Tumour-immune dynamics with an immune checkpoint inhibitor

The use of immune checkpoint inhibitors is becoming more commonplace in clinical trials across the nation. Two important factors in the tumour-immune response are the checkpoint protein programmed death-1 (PD-1) and its ligand PD-L1. We propose a mathematical tumour-immune model using a system of ordinary differential equations to study dynamics with and without the use of anti-PD-1. A sensitivity analysis is conducted, and series of simulations are performed to investigate the effects of intermittent and continuous treatments on the tumour-immune dynamics. We consider the system without the anti-PD-1 drug to conduct a mathematical analysis to determine the stability of the tumour-free and tumorous equilibria. Through simulations, we found that a normally functioning immune system may control tumour. We observe treatment with anti-PD-1 alone may not be sufficient to eradicate tumour cells. Therefore, it may be beneficial to combine single agent treatments with additional therapies to obtain a better antitumour response.


Introduction
The immune system works to differentiate between healthy cells and abnormal cells. Immune cells, known as T cells, fight the abnormal cells while leaving the healthy cells untouched. In accomplishing this, the use of 'checkpoints', molecules on immune cells that enhance or suppress an immune response, is employed (Mahoney, Freeman, & McDermott, 2015).
Two important factors in the tumour-immune response are programmed death-1 (PD-1) and its ligand PD-L1. PD-1 is a protein expressed on activated T cells. PD-L1, found mainly on various types of tumour cells and T cells (Maute et al., 2015;Talay, Shen, Chen, & Chen,2009), is one of PD-1's ligands. When PD-1 binds to PD-L1, the PD-1-PD-L1 complex is formed, allowing the cell expressing PD-L1 to be unrecognized as a danger when detected. In response to an immune attack, tumour cells may overexpress PD-L1 and PD-L2 which can bind to PD-1 receptors on T cells, reducing or eliminating the effectiveness of T cells' attacks (He, Hu, Hu, & Li, 2015).
Recently, instead of directly killing tumour cells, scientists have worked to exploit the immune system through targeting the checkpoints and thus allowing the T cells to recognize tumour cells (Mahoney et al., 2015). Drugs targeting PD-1 or PD-L1 can prevent the formation of the PD-1-PD-L1 complex. In particular, many PD-1 inhibitors (anti-PD-1s), including pembrolizumab and nivolumab, have been developed and shown positive results in over 270 international clinical trials (Liu & Wu, 2017). A 2012 clinical trial demonstrated the value of PD-1 inhibitors, as treated patients saw increases to their overall survival (Topolian et al., 2012). Therapies involving immune checkpoint inhibitors are being developed to treat a wide array of tumour types.
A mathematical model of the combination of anti-PD-1 and a tumour vaccine (GVAX) was first developed in Lai andFriedman in 2017 (Lai &. According to clinical trial data, the use of PD-1 inhibitors alone is not ideal due to high costs and low antitumour effects. Instead, single agent treatments in combination with therapies such as chemotherapy, radiotherapy and other forms of immunotherapy have shown greater promise through ongoing clinical trials (Hamanishi et al., 2016). Combination therapies are proposed with the objective to increase the antitumour response effects. Lai and Friedman examine a combination treatment with a system of 13 partial differential equations incorporating dendritic cells, activated CD4+ T cells, activated CD8+ T cells, tumour cells, necrotic tumour cells, cytokines interleukin-2 (IL-2) and interleukin-12 (IL-12), PD-1, PD-L1, PD-1-PD-L1, anti-PD-1 and the tumour vaccine GVAX (Lai & Friedman, 2017). Their model examines the synergy of the two drugs, suggesting that the immune checkpoint inhibitor and tumour vaccine work better to reduce the tumour volume in combination than individually. Their results are verified through a comparison with data collected from mice experiments. An optimal combination level of each drug is found, which lies within the bounds of the maximum tolerated dose allowed in clinical trials.
Although we appreciate the need for a combination of treatments as addressed by Lai and Friedman, we wish to solely examine the effects of a single treatment. We hope an in-depth understanding of the immune checkpoint inhibitor therapy will enhance our appreciation and comprehension of the combined treatment. The model in Lai and Friedman (2017) is a complex system of partial differential equations which forbids a comprehensive mathematical exploration and an insightful biological interpretation of the model. Our main objective is to simplify their model into a tractable system of ordinary differential equations to help gain a better mathematical and biological understanding of its rich dynamics. We consider the case when only one treatment, anti-PD-1, is applied. Our model depicts the tumour-immune interactions following the activation process of T cells by interleukins. We compare the treatment in the case when continuous drug and periodic injections are applied. In the Lai and Friedman paper (Lai & Friedman, 2017), most of the main dependent variables in the model go to steady state. We take advantage of these steady state values in our simplification, reducing the system to three nonlinear ordinary differential equations. The effectiveness of intermittent vs. continuous therapy will be examined in detail. Specifically, a study will be conducted into how altering the Table 1. Description of variables in the model system (1)-(3).

Variable
Meaning Unit C Density of tumour cells g/cm 3 T Density of activated T cells g/cm 3 P Free PD-1 concentration g/cm 3 A Anti-PD-1 concentration g/cm 3 Q PD-1-PD-L1 concentration g/cm 3 Figure 1. The tumour-immune interactions following activation of T cells. Sharp arrows indicate proliferation/activation. Blocked arrows indicate killing/blocking. Dashed lines indicate proteins on the tumour or T cells. Tumour cells express PD-L1. Activated T cells express PD-1 and PD-L1. Activated T cells kill tumour cells. When PD-1 and PD-L1 bind together, the complex PD-1-PD-L1 is formed, which inhibits the function of activated T cells. Anti-PD-1 binds to PD-1, inhibiting the formation of the PD-1-PD-L1 complex.
dosage frequency, while maintaining the same total dosage, impacts the tumour volume. We examine the drug-free system and perform an analysis on the existence and local stability of the model equilibria. Moreover, we illustrate our findings through selective simulation and bifurcation results.

Model formulation
The variables of the model, their meanings and their respective units are listed in Table 1. The essence of our simplified mathematical model is captured by the following dramatically simplified tumour-immune interactions network following activation of T cells shown in Figure 1.
Our model of tumour-immune interaction with immune checkpoint inhibitor anti-PD-1 treatment is a vastly simplified and slightly modified version of the model presented in Lai and Friedman (2017). Our model consists of three nonlinear ordinary differential equations of the following form: dT dt = λ TI 12 T N I 12 K I 12 + I 12 activation by IL-12 where P is defined as: The function for suppression of T cell activation and proliferation by the PD-1-PD-L1 complex is given by: As in Lai and Friedman (2017), we assume the dissociation (d Q Q) and association (α PL PL) of the complex PD-1-PD-L1 are fast so that we can apply the usual quasi-steady-state argument to obtain α PL PL = d Q Q. Thus, As in Lai and Friedman (2017), we assume the level of free PD-L1 can be represented by since PD-L1 is expressed on both T cells and tumour cells, and upregulated on tumour cells, indicated by C > 1. Thus The parameters, their meanings and estimates of their values are given in Table 2. A more detailed discussion of the parameter estimates is given in Appendix 1. The cytokines IL-2 and IL-12 were estimated from their steady-state values. In our model, we consider the process after the T cell population has been activated. Our hope is to focus primarily on the interactions between the T cells, cancer cells and PD-1-PD-L1 complexes. As a result, we base the values of IL-2 and IL-12 densities off the control case simulated in Figure 2 of Lai and Friedman (2017). The simulation shows the IL-2 and IL-12 densities are not very dynamic, as they both move within a tight range. Since the densities tend to a steady-state in Lai and Friedman's simulations, we assume a steady-state in our model and take IL-2 = 2.37 · 10 −11 g/cm 3 and IL-12 = 1.5 · 10 −11 g/cm 3 . Equation (1) models the tumour cell population, where logistic growth is used to denote the net growth of the tumour cells in the presence of an anti-PD-1 immune checkpoint inhibitor (Kirschner & Panetta, 1998). Through interactions between tumour cells and activated T cells, the loss of tumour cells occurs at rate η. We assume that the PD-1/PD-L1 interaction reduces the number of activated effective T-cells with negligible inhibition on the immune cytotoxicity. Equation (2) captures the rate of change of the activated (effective) T cell population. IL-12 and IL-2 play a major role in activation and differentiation of the population of T-cells, shaping the immune response. The first source term accounts for the activation of naive T cells (T N ) by IL-12. The other source term represents the proliferation of activated T cells, as stimulated by IL-2. Michaelis-Menten forms are used to capture the saturation in the immune response (Kirschner & Panetta, 1998). Larger values of PD-1-PD-L1 concentration inhibit the activation and proliferation of T cells, reducing the source terms by a factor of 1 1+Q/K TQ . Natural death of T cells occurs at rate d T . The anti-PD-1 treatment is modeled by Equation (3), where A(t) = γ A corresponds to an intravenous, continuous injection. Intradermal injections are also considered with varying dosage frequencies, as discussed in Numerical Simulations. PD-1 binds to anti-PD-1 at a rate of μ PA . We assume there is no dissociation of PD-1 with anti-PD-1 once bound due to the very quick nature of the association and dissociation. The rates of association and dissociation are relatively small values and lead to an assumption that the variables can be approximated by their levels at quasi-steady state. Anti-PD-1 naturally decays at rate d A . Finally, Equation (4) describes the density of the free PD-1. Newly activated T cells have a fixed amount of PD-1 per T cell, ρ P . As anti-PD-1 is injected and the drug binds to PD-1, PD-1 is depleted at a rate of γ A.

Treatment free model dynamics
It is easy to see that the functions contained in the full system (1)-(3) are differentiable, which ensures its solutions with positive initial values exist and are unique by a direct application of standard differential equation theory. Observe that the full system (1)-(3) consists of three highly nonlinear differential equations that prevent us from conducting a comprehensive mathematical analysis. In order to appreciate the full model complexity and to gain an in-depth understanding of the natural tumour-immune interaction via the lens of a mathematical model, we examine the case when no anti-PD-1 treatment is applied (A = 0) in the full system (1)-(3). In this case, Equation (4) becomes P = ρ P T and Q = βT(T + C C), β = κρ P .
The reduced system takes the form: where we use the assumptions on the complex-forming process to define F(C, T): M and N are two positive parameters defined as: M = λ TI 12 T N I 12 K I 12 + I 12 ,  (2017) T N Density of naive T cells 6 · 10 −4 g/cm 3 Lai and Friedman (2017) I 12 IL-12 concentration 1.5 · 10 −10 g/cm 3 Lai and Friedman (2017) K I 12 Half-saturation of IL-12 1.5 · 10 −10 g/cm 3 Lai and Friedman (2017) λ TI 2 Activation rate of T cells by IL-2 0.5 day −1 Lai and Friedman (2017) I 2 IL-2 concentration 2.37 · 10 −11 g/cm 3 Lai and Friedman (2017) Half-saturation of IL-2 2.37 · 10 −11 g/cm 3 Lai and Friedman (2017 We seek to examine the local and global stability for our system. In doing so, we first need to establish positivity and boundedness of the solutions with positive initial values. Proposition 3.1: Solutions of (5) that start positive remain positive and bounded.

Proof:
We first seek to show the proof of positivity, examining T first. We assume T(t 0 ) ≥ 0. Then for T(t) < 0 for some t > 0, we would require dT dt ≤ 0 when T = 0. However, we have dT dt | T=0 = M > 0 since all parameters are positive. Thus T(t) > 0, ∀t > 0.
In examining the positivity of C, we assume C(t), T(t) exist on [0, t 1 ]. By a standard separation of variables, We can then conclude Since solutions remain nonnegative, we now seek to prove boundedness. We start with boundedness of C using a comparison argument: In addition, we see that lim Therefore C is eventually bounded above by C K + a for any positive constant a.
We then look to prove that T is bounded. First, let G(C, T) = (M + NT)F(C, T). Then By a standard comparison argument, Then T ≤ max{T 0 , α d T } and T is bounded. The reduced system (5) may contain two types of equilibria: the tumour-free equilibria E * 0 = (0, T * 0 ) and tumorous equilibria E * 1 = (C * , T * ), C * > 0, T * > 0. Before examining stability, we must establish the existence and uniqueness of the tumour-free equilibrium Proposition 3.2: The system (5) has a unique tumour-free equilibrium is a strictly increasing function of parameter M. Proof: The T * values corresponding to C * = 0 are the roots of the equation: Since p(0) = −M < 0 and lim T→∞ p(T) = +∞, we see that p(T) = 0 has a unique positive root which is a strictly increasing function of M.
If d T < N, with M = 0, then p(T) = 0 has three real roots: However, when M > 0, the entire curve is shifted down, as shown in Figure 2. There is only one real positive root T * 3 , and T * 3 increases strictly as M increases. We now establish the existence and uniqueness of the tumorous equilibrium E * 1 = (C * , T * ), C * > 0, T * > 0.

Proposition 3.3: The system (5) has a unique tumorous equilibrium E
Proof: It is easy to see that These together yield The number of positive roots of Equation (8) can be determined by the Descartes' rule of signs. For example, if C C K η λ C < 1 and N d T < 1, then the polynomial f (T * ) in Equation (8) has only one change of signs in coefficients and hence has a single positive root; while if (8) has two changes of signs in coefficients and hence may have two positive roots or none.
We also need to ensure the positivity of C * . We have We are now able to perform a stability analysis with our reduced system. The Jacobian for the system (5) is given by: We first examine the stability of the tumour-free equilibrium E * 0 = (0, T * 0 ). Proposition 3.4: Assume that d T > N. Then for system (5) the following statements are true: For part (a) and the first part of (b), the Jacobian matrix for the tumour-free equilibrium Since we know F(C, T) ≤ 1 and F(C, T) is a decreasing function of C and T, ∂ ∂T F(C, T), ∂ ∂C F(C, T) < 0 and the eigenvalues are: Thus, we have local asymptotic stability for E * 0 . Similarly, if λ C > ηT * 0 , then E * 0 is a saddle point.
We can biologically interpret several of the conditions for the local stability. We have that λ C < ηT * 0 implies that during the onset of tumour initiation, resident T cells can kill the cancer cells faster than the cancer cells can multiply. Thus, it is unsurprising to find that this condition is necessary for the local stability of the tumour-free equilibrium. The condition d T > N can be interpreted as the death rate of the T cells being greater than the stimulation rate of T cells by IL-2.
We next examine the stability of the tumorous equilibrium E * 1 = (C * , T * ). Proposition 3.5: Assume that C C K η λ C < 1, N d T < 1 and f λ C η < 0, then the tumorous steady state E * 1 of system (5) is unique and is locally asymptotically stable. Proof: The Jacobian matrix for the tumorous equilibrium E * 1 is: Assume that d T > N. Then the trace for the tumorous equilibrium is Therefore, τ < 0. The determinant is given by: We can rearrange the terms in : Observe that ∂Q ∂T > 0 and ∂Q ∂C > 0. We see that > 0 if and only if which is equivalent to Clearly this is true if η C < 2 λ C C K or 1 λ C η C C K < 2. The first condition of Proposition 3.3 ensures this is true. Therefore, > 0. Together with Proposition 3.3, we have the following local stability result for the tumorous equilibrium E * 1 . We now establish global asymptotic stability of the tumorous steady state E * 1 of system (5). To this end, we need to identify conditions that can rule out the existence of nontrivial positive periodic solutions from system (5). Proposition 3.6: The system (5) has no nontrivial positive periodic solutions provided that d T > N.

Proof:
We employ the Dulac criterion to show there are no nontrivial positive periodic orbits. Let h(C, T) = 1/C. We have: Thus, the Dulac criterion ensures there will be no nontrivial positive periodic solutions for system (5).
We are ready to state the following global stability results for system (5). Theorem 3.7: Assume that d T > N. Then the following are true: (i) If λ C < ηT * 0 and no tumorous equilibrium exists, then the tumour-free equilibrium E * 0 of system (5) is globally asymptotically stable. (ii) If C C K η λ C < 1 and f λ C η < 0, then the tumorous steady state E * 1 of system (5) is globally asymptotically stable.
Proof: Recall solutions of system (5) with respect to positive initial values are positive and bounded. By Poincare-Bendixson Theorem, bounded positive solution must tend to a nonnegative equilibrium, a nontrivial positive periodic solution, a homoclinic cycle or a heteroclinic cycle. In case of (i), the system (5) admits no positive steady states, nontrivial positive periodic solutions, homoclinic cycles or heteroclinic cycles. This implies that all positive solutions of system (5) tend to the tumour-free equilibrium E * 0 . In case of (ii), the system (5) admits a unique positive steady state E * 1 , which is locally asymptotically stable. This eliminates the existence of heteroclinic cycles since system (5) admits only two equilibria. It is easy to see that if E * 0 is locally asymptotically stable, then S148 E. NIKOLOPOULOU ET AL.  there can not be any homoclinic cycles originating from E * 0 . Observe that the T axis is invariant. If there is any homoclinic cycle originating from E * 0 , then the cycle will contain at least one branch of its stable manifold, which is unbounded, a contradiction to the fact that positive solutions of the system (5) are bounded. Together, these arguments exclude the possibility of the existence of any homoclinic cycle originating from E * 0 . By Poincare-Bendixson Theorem, the tumorous steady state E * 1 of system (5) attracts all its positive solutions.

Sensitivity analysis
By the method outlined in Marino, Hogue, Ray, and Kirschner (2008), we performed a sensitivity analysis using partial rank correlation coefficients (PRCC) via Latin hypercube sampling (LHS) to identify the main drivers of the tumour burden. We used LHS to produce 1000 sets of randomly determined parameters, which were then used to compute the PRCC and corresponding p-values. Parameters were tested within the ranges given in Table 2. For parameters without ranges, sampling was completed within a range from 1/10 to twice the values in Table 2. The sensitivity analyses with respect to the tumour cell density C and the activated T cell density T at day 200 are shown in Figures 3 and 4, respectively.
The density of tumour cells is most sensitive to changes in the tumour cell growth rate (λ C ), kill rate of tumour cells by T cells (η) and death rate of T cells (d T ) (p-value < .01 for each). We notice that C, the tumour cell density, is highly positively correlated with the tumour cell growth rate (λ C ) and the death rate of T cells (d T ), and negatively correlated with the kill rate of tumour cells by T cells (η).
The sensitivity analysis conducted with respect to T, the activated T cell density, reveals the tumour cell growth rate (λ C ), death rate of T cells (d T ) and kill rate of tumour cells by T cells (η) are the main parameters influencing T (p-value < .01 for each). As could be expected, the activated T cell density is negatively correlated with the growth rate of the tumour cells (λ C ) and death rate of T cells (d T ). The kill rate of the tumour cells by the T cells (η) captures the most highly positively correlated parameter.
In making biological sense of the results, we note that the reaction of immune cells (T cells) and tumour cells can either result in the death of tumour cells or the inactivation of the immune cells. Our model solely represents the former case. Thus, the density of tumour cells and activated T cells are sensitive to d T , the death rate of T cells, which causes the T cells to decline. As the activated T cells decline, there are fewer reactions, and thereby less deaths of tumour cells. The kill rate of tumour cells by T cells (η) increasing or decreasing reflects a stronger or weaker immune response. As η increases, the tumourimmune interactions result in more deaths of tumour cells. When the tumour cells are decreased, less PD-L1 is in the system to bind to PD-1. Thus, less PD-1-PD-L1 complexes are formed to inhibit the activation of T cells. Therefore, η is positively correlated with the activated T cell density and negatively correlated with the density of tumour cells. Finally, the production parameter for the tumour cells, λ C , is positively correlated with the density of tumour cells. As the number of tumour cells, and thereby the amount of PD-L1, increases, more T cells are being inactivated as a result of the increased PD-1/PD-L1 bindings. Therefore, λ C is negatively correlated with the activated T cell density.
The bifurcation plot of d T , as shown in Figure 5, suggests the existence of a critical value for the maximum death rate of T cells d T , namely d T ≈ N, since N = 0.25. When d T < N, this causes an increase in the number of activated T cells, resulting in the killing of more tumour cells. However, when d T > N, the tumour cells dramatically increase towards the carrying capacity. The plot suggests our mathematical results can be improved upon to find sharper conditions, as d T > N is not required for the tumour-free equilibrium E * 0 to be stable. We can now discern that Proposition 3.4 and 3.5 provide sufficient, but not necessary conditions for the stability of the tumour-free and tumorous equilibria. Further exploration of the bifurcation plots for different parameters illustrate the approach of the same steady state.

Numerical simulations
All simulations of the mathematical model were performed by MATLAB using the parameter estimates given in Table 2 and the initial conditions (in units of g/cm 3 ) C(0) = 0.3968, T(0) = 6 · 10 −3 , and A(0) = 0 (Lai & Friedman, 2017). We examine the cases when no drug, a continuous drug and injections of varying periodicity are applied.
We first explore the control case when no drug is applied ( Figure 6). We keep the death of the T cells small enough so that the results closely resemble the behaviour seen for those with a normally functioning immune system. The simulation is run for a duration of 150 days to determine the behaviour of certain cell densities/concentrations. We observe that the number of activated T cells increases up to a certain level, and later stabilizes, approaching a steady state. The density of tumour cells decreases significantly until reaching a disease free equilibrium. Since the death rate of T cells is small, tumour is controlled to a great extent. Consequently, people with a strong immune system have a chance to eradicate tumour. However, for those with a weakened immune system, tumours are significantly less responsive to treatment. We next consider intermittent anti-PD-1 therapy compared to continuous treatment. We seek to determine how the densities of tumour and activated T cells are affected by varying dosage frequencies and maintaining the dosage level. The time period between injections is varied from 3 days to 150 days. Simulations are run for 365 days to reflect a clinically realistic treatment period. Figure 7 suggests a periodicity in the immune response, resulting from the treatment acting as a forced oscillator. Each time the anti-PD-1 is administered, the level of activated T cells increases due to less PD-1-PD-L1 complexes being formed and thereby less inhibition of T cell activation. As a result, the tumour cells sharply decrease. Following the decline, the tumour cells logistically grow, returning to their previous level. A continuous therapy offers the best control for the tumour cells. Figure 8 suggests that when treatments are less frequent (every 60, 90, 150 days), the tumour cells periodically approach the same density as when no drug is applied before decreasing with the next treatment. The more frequent treatments (every 3, 10, 30 days) periodically approach a lower density, closer to that of the continuous drug. Thus, the simulations suggest that more frequent treatments produce more effective results. As the tumour-immune reactions and PD-1-PD-L1 complex formations are fast, the drug is entirely eliminated from the system before the next treatment begins. When the time between treatments is smaller, the drug remains continuously in the system and is thus more effective. An interesting area we wish to explore with our model is the behaviour of the tumourimmune system when the treatment is abruptly ended (switch on/off treatment). We have chosen a period of 180 days (approximately half of a year) to apply the treatment, after which the treatment will be terminated, as shown in Figure 9. Once the treatment has ended, the density of the tumour cells and activated T cells eventually approach the same steady state as if no drug had ever been applied. Similar results were found when the time interval before treatment termination was increased.

Discussion
Tumour immunotherapy is one of the current focuses in oncology. One common form of immunotherapy employs the use of immune checkpoint inhibitors, such as anti-PD-1, in a variety of tumours, including melanoma and lung cancer. A growing number of clinical trials using drugs with this mode of action seek to improve survival rates and quality of life. Current trials suggest that the use of a single agent such as anti-PD-1 produces a low tumour-immune response. Combinations of anti-PD-1 with other types of treatments are the recommended protocol to produce the most effective results (Hamanishi et al., 2016). The primary purpose of this work is to gain insight into the tumour-immune response while highlighting immune checkpoint therapy. Due to lack of data, many of the parameter values are difficult to estimate. Therefore, we developed a simplified model in an effort to explore interesting mathematical results. In particular, we assume that association and dissociation of PD-1/PD-L1 are fast enough to be at a steady state at the scale of tumour growth. Despite this simplification, our simulations reflect expected dynamics. Although we cannot currently test the predictive capability of our model because of sparse data, we can use the results to predict general behaviours of the tumour-immune response to treatment.
Our approach contrasts and complements that of Lai and Friedman (2017). Their model is a system of 13 partial differential equations with 48 parameters. Their goal was to explore the synergy of two drugs (anti-PD-1 and GVAX). One of Lai and Friedman's major contributions is a detailed parameterization. Unfortunately, many of these parameters are unexplored in the empirical literature, usually due to overwhelming technical challenges. Because of these uncertainties, we arrived at quite wide ranges for many of the parameters. In some cases, these ranges were not in agreement with Lai and Friedman. For example, since many tumour cells express PD-L1, they set C , the expression of PD-L1 in tumour cells vs. T cells, to be 0−0.01. However, we concluded that C = 1−100 based on literature indicating PD-L1 is upregulated in some tumour types (Spranger et al., 2013). Due to a lack of data, our model cannot be validated, so there is limited predictive capability. However, it provides insight into how the tumour and T cell densities depend on parameter values. Our sensitivity analysis using PRCC via LHS showed λ C , the tumour cell growth rate, η, the kill rate of tumour cells by T cells, and d T , the death rate of T cells, were the most sensitive parameters and are thus the main drivers of the tumour burden. Furthermore, the bifurcation analysis for d T confirmed the findings of the sensitivity analysis. The bifurcation diagram suggests that the stability of the tumour-free and tumorous equilibria changes as d T varies. These results are biologically significant. If the lifetime of T cells lasts less than four days, the tumorous equilibrium becomes stable, resulting in a significant disease. Thus, our model suggests that cancers are likely to be more common and significant for patients with weakened immune systems, whether from age, emotional or physical stress. Our model further suggests that, if caught early enough before the immune system becomes too weak, it is possible to eradicate the tumour with minimal treatment.
The simulations suggest anti-PD-1 is most effective when it is applied continuously. Without a continuous or frequent injection of the drug, the end results are as if no drug was applied. However, even when applied continuously, it is not a sufficient treatment to eradicate the tumour cells. It would be most beneficial to pair anti-PD-1 with another treatment. Previous anti-PD-1 trials reached similar conclusions (Antonios et al., 2016). When combined with a dendritic cell vaccination, the treatment produced a substantial increase in survival. However, when anti-PD-1 alone was used, the benefit was not notable.
Recent clinical trials suggest that combining the immune checkpoint inhibitor with an alternative treatment is the most effective way to treat tumours (Antonios et al., 2016). Anti-PD-1 is not fully effective due to a lack of T cells (Kleponis, Skelton, & Zheng, 2015). Although a single-agent checkpoint inhibitor increases T cell density, our simulations suggest anti-PD-1 alone cannot raise the activated T cells to the level needed to eradicate the tumour cells (Figure 7). When combined with a vaccine or therapy that increases the amount of T cells, anti-PD-1 has a more substantial anti-tumour immune response, as it can help more T cells remain active in killing tumour cells (Lai & Friedman, 2017).
Our stability analysis of the drug-free system established requirements for global stability of the tumorous and tumour-free equilibria. Two conditions are required to eliminate the tumour: λ C < ηT * 0 and d T > N. Biologically, we see that λ C < ηT * 0 means that when the tumour is initially growing, the T cells can kill the cancer cells faster than the cancer cells multiply. The condition d T > N can be interpreted as the death rate of the T cells being greater than the stimulation rate of T cells by IL-2. Thus, the mathematical analysis suggests that even patients with weakened immune systems, and thereby shorter lifetimes of T cells, can see tumour elimination for slow-growing tumours. Though the condition d T > N complements the mathematical analysis, the bifurcation diagram suggests our mathematical results still have room for improvement. This condition proves to be sufficient but not necessary for determining the stability of the tumour-free and tumorous equilibria. Further exploration is required to determine necessary conditions. Future work will include two competing tumour species: one sensitive and the other insensitive to treatment. Such a model promises to help determine the most effective frequency and dose of treatment, corresponding to patient-specific initial conditions, to delay the onset of treatment resistance. To investigate the most effective reduction of the tumour burden, it would be beneficial to examine the effects of adding additional treatments, like chemotherapy or radiation, to the model. Hence λ C = 0.18 − 0.67 day −1 . The carrying capacity of tumour cells can be estimated using the assumption that a tumour reaching 1 gram net weight contains 10 9 cells (Del Monte, 2009). By Szomolay, Eubank, Roberts, Marsh, and Friedman (2012), the total capacity of living and dead tumour cells can be estimated to be 9.45 · 10 8 cells/cm 3 . Then 9.45 · 10 8 cells 1 cm 3 · 1 gram 10 9 cells = 0.945 g/cm 3 By Friedman and Hao (2017), the carrying capacity of tumour cells is 0.8 g/cm 3 . Hence C K = 0.8 − 0.945 g/cm 3 . By Lai and Friedman (2017), the killing rates of tumour cells by CD4+ and CD8+ T cells are estimated to be 11.5 day −1 · cm 3 /g and 46 day −1 · cm 3 /g, respectively. Thus, the killing rate of tumour cells by T cells is η = 11.5 + 46 = 57.5 day −1 · cm 3 /g. Equation (2): By Lai and Friedman (2017) the activation rate of CD4+ and CD8+ T cells by IL-2 are both estimated to be 0.25 day −1 . Hence, the activation rate of T cells by IL-2 is λ TI 2 = 0.25 + 0.25 = 0.5 day −1 . Similarly, following the same process for the activation rate of T cells by IL-12, we have λ TI 12 = 4.66 + 4.15 = 8.81 day −1 . We consider the density of naive T cells as T N = 4 · 10 −4 + 2 · 10 −4 = 6 · 10 −4 g/cm 3 , as the density of naive CD4+ and CD8+ T cells are estimated to be 4 · 10 −4 g/cm 3 and 2 · 10 −4 g/cm 3 , respectively (Lai & Friedman, 2017).
By Grossman and Berke (1980), the lifetime of T cells is assumed to be 50 days (d T = 1/50 = 0.02 day −1 ). By McDonagh and Bell (1995), the lifespan for CD8+ T cells is given to be 68 days (d T = 1/68 = 0.014 day −1 ). When fit to experimental data of a BCL 1 tumour in the spleens of chimeric mice, the death rate of T cells was estimated as d T = 0.0412 day −1 (Kuznetsov et al., 1994). In Eikenberry et al. (2009), the death rate of T cells is estimated to be d T = 0.0 − 0.05 day −1 . Collecting the findings, we take d T = 0.0 − 0.05 day −1 .
In Cheng et al. (2013), experiments found an average of 9282 PD-L1 on a single activated T cell, with a general range of about 8000 − 11500 PD-L1 on a single activated T cell. Following the same procedure as in Lai and Friedman (2017), we estimate ρ L as follows: where m L and m T are the mass of one PD-L1 protein and the mass of one T cell (T 1 or T 8 ). We assume that one T cell is approximately 10 microns in diameter, with the same density as water at 39 • C, 0.9926 g/cm 3 . Then we can calculate the mass of one T cell: m T = 4π(0.005) 3 cm 3 3 · 0.9926 g cm 3 ≈ 5.24 · 10 −10 g.
Thus we assume L ≈ 2 · 10 −6 g/cm 3 at steady state. Equation (4): In Cheng et al. (2013), experiments showed an average of 3096 PD-1 on a single activated cell, with a general range of 2000 − 5000 PD-1 on a single activated cell. Following the same procedure as in Lai and Friedman (2017), we estimate ρ P as follows: where m P and m T are the mass of one PD-1 protein and the mass of one T cell. By Agata et al. (1996), a flow cytometry shows the mass of one PD-1 protein is m P = 50 − 55 kDa (8.3 · 10 −20 − 9.13 · 10 −20 g). As outlined above, m T = 5.24 · 10 −10 − 10 −9 g. We then calculate: