Modelling and Simulation of Complex Biological SystemsView this Special Issue
Mathematical Modelling of the Inhibitory Role of Regulatory T Cells in Tumor Immune Response
The immune system against tumors acts through a complex dynamical process showing a dual role. On the one hand, the immune system can activate some immune cells to kill tumor cells (TCs), such as cytotoxic T lymphocytes (CTLs) and natural killer cells (NKs), but on the other hand, more evidence shows that some immune cells can help tumor escape, such as regulatory T cells (Tregs). In this paper, we propose a tumor immune interaction model based on Tregs-mediated tumor immune escape mechanism. When helper T cells’ (HTCs) stimulation rate by the presence of identified tumor antigens is below critical value, the coexistence (tumor and immune) equilibrium is always stable in its existence region. When HTCs stimulation rate is higher than the critical value, the inhibition rate of effector cells (ECs) by Tregs can destabilize the coexistence equilibrium and cause Hopf bifurcations and produce a limit cycle. This model shows that Tregs might play a crucial role in triggering the tumor immune escape. Furthermore, we introduce the adoptive cellular immunotherapy (ACI) and monoclonal antibody immunotherapy (MAI) as the treatment to boost the immune system to fight against tumors. The numerical results show that ACI can control TCs more, while MAI can delay the inhibitory effect of Tregs on ECs. The result also shows that the combination of both immunotherapies can control TCs and reduce the inhibitory effect of Tregs better than a single immunotherapy can control.
Tumors can be benignant (not cancerous), premalignant (precancerous), and malignant (cancerous). Every year millions of people suffer with cancer and die from this disease throughout the world . It is important to understand tumor’s mechanisms of establishment and destruction, cell-mediate immunity with cytotoxic T lymphocytes (CTLs), and natural killer cells (NKs), generally called effector cells (ECs) that are cytotoxic to tumor cells (TCs), and play a basic role in immune response against tumors [2, 3]. Moreover, efficient antitumor immunity requires the action of helper T cells (HTCs), which can directly activate naive CD T cells to differentiate into CTLs [4–6]. Recently, it has been reported that regulatory T cells (Tregs) can inhibit CTLs and promote the escape of TCs . Tregs suppress immune cells, and when the war between T cells and infection is over, the Tregs signal to stop . Cancer immunotherapy fights against cancer by strengthening the body’s immune system, but the involvement of Tregs inhibits the immune response and turns off the anticancer effect. Tregs inhibition is important in the dynamics of the tumor immune system, which is one motivation of this work.
Adoptive T cell immunotherapy (ACI) as a common immunotherapy involves injecting adoptive T cells directly into tumor patients [9–11]. Its advantages are good destruction of tumor and persistence, while its disadvantages are serious toxic and side effects. The monoclonal antibody immunotherapy (MAI) is the immune checkpoint inhibitor [12, 13], which has the advantage of removing the suppression state of the immune system and restoring the immune function of the body to TCs and the disadvantage of having serious immune-related adverse reactions.
Tregs have become an important target in tumor immunotherapy because of their contribution to tumor immune escape. Cytotoxic T lymphocyte antigen 4 (CTLA-4) is a marker that is expressed on the surface of activated T cells and transmits inhibitory signals in the immune response [14–16]. Blocking CTLA-4 can reduce the inhibitory activity of Tregs, and the anti-CTLA-4 humanized monoclonal antibody Ipilimumab and Tremelimumab are used to treat advanced melanoma and malignant mesothelioma, respectively . Similar to CTLA-4, programmed death receptor 1 (PD-1) can also promote the activation and development of Tregs [18–21]. Blocking PD-1 can prevent the development of Tregs and prevent the conversion of HTCs into Tregs . Currently, OPDIVO (Nivolumab), an anti-PD-1 monoclonal antibody, has been approved by the US FDA for the treatment of melanoma, renal cell carcinoma, and non-small cell lung cancer [23–25]. Establishment of a mathematical model to study the immunotherapy on the reduction of Tregs inhibition has both theoretical and practical significance.
In order to describe the mechanisms of host’s own immune response to against TCs, various types of mathematical models have been proposed [26–43]. The modelling of the tumor immune system described by ordinary differential equations (ODEs) has a long history, which can be traced back to the classic research of Stepanova in 1980 . In 1994, Kuznetsov et al. established the famous two-dimensional ODEs model, postulating that tumor growth follows the Logistic growth pattern. They evaluated the parameters of the model by fitting experimental data from mice . In 2003, Stolongo-Costa et al. assumed that TCs follows the exponential growth pattern and constructed a two-dimensional ODEs model. They analyzed the basic properties of the model and provided conditions for stability of the tumor-free equilibrium, explaining its epidemiological significance . In 2004, Galach simplified Kuznetsov’s system to account for the effect of immune delay on the tumor immune system . In 2014, Dong et al. constructed a three-dimensional ODEs model focusing on the effects of HTCs on the tumor immune system .
In 1998, Kirschner and Panetta generalized Kuznetsov–Taylor model and illustrated the dynamics between TCs, ECs, and IL-2. They firstly introduced ACI into their model which can explain both short-term tumor oscillations in tumor sizes as well as long-term tumor relapse . In 2003, in order to study the role of cytokine therapy in the activation of the immune system, Stolongo-Costa et al. introduced cycle therapy term and established a cycle immunotherapy model. They obtained some thresholds of the frequency and intensity of immunotherapy . In 2006, de Pillis et al. constructed the six-dimensional ODEs model to investigate the effects of combined chemotherapy and immunotherapy on tumor control. They briefly analyzed the nature of the model and discussed the optimal treatment using optimal control theory . In 2008, Bunimovich-Mendrazitsky et al. established a pulsed differential equation model with Bacillus Calmette-Guerin tumor immunotherapy. They obtained the critical threshold and pulse frequency of BCG injection dose that could successfully treat superficial bladder cancer . In 2012, Wilson and Levy established a mathematical model containing Tregs. They studied the absence of treatment, vaccine treatment, anti-TGF treatment, and combination vaccine and anti-TGF treatment, as well as sensitivity analysis of some important parameters . In 2018, Radunskaya et al. established a mathematical model with blood, spleen, and tumor compartments to study PD-L1 inhibitors in the role of tumor immunotherapy. The model was used to fit parameters with the experimental data. The results showed that increasing the resistance of PD-L1 doses can greatly improve the clearance rate of tumor .
This paper investigates the role of Tregs in the tumor immune system. Therefore, we incorporate the fourth population of Tregs into the previous system in . For the mathematical simplicity, a bilinear term also has been used to describe the interactions between immune response and tumor. To our knowledge, HTCs can recognize TCs and promote the growth of ECs. And ECs can provide direct protective immunity by attacking TCs. When there are more HTCs and ECs, in order to maintain immune homeostasis, the body will produce corresponding Tregs to suppress ECs, and Tregs originating from both HTCs and ECs. Then, we establish a four-dimensional ODEs model described as below:where , , , and represent the populations of TCs, ECs, HTCs, and Tregs, respectively. The first equation describes the rate change of TCs population. The tumor follows logistic growth dynamics with growth rate , and the maximum capacity is . represents the loss rate of TCs by ECs interaction. The second equation describes the rate change of the ECs population. is the mortality rate of ECs. is the activation rate of ECs by HTCs, and is the inhibition rate of Tregs on ECs. The third equation describes the rate change of the HTCs population, is birth rate of HTCs produced in the bone marrow, and HTCs have a natural lifespan of an average days. is HTCs stimulation rate by the presence of identified tumor antigens. The fourth equation gives the rate change of the Tregs population, and are the activation rates of Tregs by ECs and HTCs, respectively. represents per capita decay rate of Tregs. A diagram of the various interactions between these cell populations is shown in Figure 1.
We nondimensionalize model (1) by taking the following scaling:and we choose the scaling . By replacing by , we obtain the following scaled model:with initial conditions
Here, , and denote the dimensionless densities of TCs, ECs, HTCs, and Tregs populations, respectively.
2. Model Analysis
2.1. Well Posedness of Model (3)
Proof. Since the right-hand side of model (3) is completely continuous and locally Lipschitz on the interval , there exists a constant such that the solutions of model (3) with initial conditions (4) are existent and unique on the interval , where .
From model (3) with initial conditions (4), we obtainIt is obvious to see that the solutions of model (3) with initial conditions (4) are nonnegative for all . The proof is complete.
2.2. Existence of Equilibria
Putting yields the tumor-free and ECs-free equilibrium, namely,which always exists.
Putting yields the tumor-free equilibrium, namely,which exists when and .
Putting yields the tumor-dominant equilibrium, namely,which exists when .
Putting and eliminating in (6), we havewhere .
Then, we have the coexistence equilibriumwhich exists when , where is the positive root of (10).
Below, we consider the existence condition of : Case (1): , i.e., . Since and , we have , which contradicts with . Hence, does not exist in this case. Case (2): , i.e., .
If , then (10) does not have a positive root. Hence, does not exist in this case.
If , then (10) has one positive root , where , which contradicts with . Hence, does not exist in this case.
If , then (10) has two positive roots, where and . Since , which contradicts with . If , i.e., and , there exists a unique coexistence equilibrium , where .
The existence conditions for each equilibrium are given in Table 1.
2.3. Stability of Equilibria
In order to investigate the local stability of the above equilibria of system (3), we linearize the system and obtain Jacobian matrix at each equilibrium :
The corresponding characteristic equation is
Theorem 1. System (3) always has one tumor-free and ECs-free equilibrium , which is unstable.
Proof. At , characteristic (13) becomesIt can be seen that one eigenvalue is positive. Hence, is unstable.
Theorem 2. System (3) has one tumor-free equilibrium , when and , which is locally asymptotically stable (LAS) if the inequality holds.
Proof. At , characteristic (13) becomesThen, one root of characteristic equation is . It is easily noted that as , so has solutions with negative real parts. Therefore, is LAS when , that is, .
Theorem 3. System (3) has one tumor-dominant equilibrium , when , which is LAS if the inequality holds.
Proof. At , characteristic (13) becomesThe three roots of characteristic equation are , and , if . So, is LAS when , which is equivalent to .
At , characteristic (13) becomeswhereNote that and are necessary for the existence of , then we have . By the Routh–Hurwitz criterion, the roots of (17) have only negative real parts if and only ifHence, we obtain the sufficient conditions for stability of :where
2.4. Numerical Simulations
In this section, we choose some suitable parameters in (3) to simulate numerically the theoretical conclusions obtained in the previous sections by using the Matlab software package MATCONT . We select the following parameters’ set :
Note that . By calculations, we have
And we find the stability region of (see Figure 2). is LAS in regions I and II. is unstable in region III.
Let us denote the point in plane as . Case (a): we choose a point in the region I; then, system (3) has one interior equilibrium: The eigenvalues of Jacobian matrix of (12) are so is stable, as shown in Figure 3(a). Case (b): we choose a point in the region II; then, system (3) has one interior equilibrium: The eigenvalues of Jacobian matrix of (12) are so is stable, as shown in Figure 3(b). Case (c): we choose a point in the region II, then system (3) has one interior equilibrium: The eigenvalues of Jacobian matrix of (12) are so is stable, as shown in Figure 3(c). Case (d): we choose a point in the region III; then, system (3) has one interior equilibrium:
Below, we perform numerically bifurcation analysis of against for different values of . Case (1): we choose (see Figure 4(a)). We obtain the following result.
Proposition 2. When , exists and is unstable, exists and is LAS, exists and is unstable, and does not exist. When , exists and is unstable, exists and is unstable, exists and is unstable, and exists and is LAS. When , exists and is unstable, does not exist, exists and is unstable, and exists and is LAS. When , exists and is unstable, does not exist, exists and is LAS, and does not exist. Case (2): we choose (see Figure 4(b)).We obtain the following result.
Proposition 3. When , exists and is unstable, exists and is LAS, and does not exist. When , where is a Hopf bifurcation, exists and is unstable, exists and is LAS, and exists and is LAS. When , exists and is unstable, exists and is unstable, and exists and is unstable. When , exists and is unstable, does not exist, and exists and is unstable. When , exists and is unstable, does not exist, and does not exist.
Consider the case where HTCs stimulation rate () is low by the presence of identified tumor antigens (). When the inhibition rate of Tregs to ECs is lower than , the solution of system (3) approaches implying that ECs can still effectively remove TCs. When , the solution of system (3) approaches showing the coexistence of TCs and immune cells, which means that the patient can survive with tumors. When , TCs escape the control of the immune system and develop into malignant tumors.
Next consider the case where HTCs stimulation rate () is high by the presence of identified tumor antigens (). When , the solution of system (3) approaches implying that TCs can be effectively removed by ECs. When , the solution of system (3) approaches implying that TCs can still be controlled by the immune system. When , the system has a Hopf bifurcation point and induces a limit cycle (see Figure 3(e)). In the biological sense, it can be understood that the number of TCs presents a periodic change.
3. Treatment Model
In order to investigate well the effect of Tregs in tumor immune response under the treatment, we follow the way in  to introduce the constant treatment term into the second equation of (1). Since CTLA-4 and PD-1 can inhibit the development of Tregs and prevent the transformation of HTCs into Tregs , the effect of CTLA4 or PD-1 on Tregs can be considered. Then, we shall establish a five-dimensional ODEs model described as below:where represents the concentration of monoclonal antibody in human body at time , represents the treatment term of introducing LAK and TIL into the region of tumor localization, represents the inhibition rate of monoclonal antibody on Tregs, represents the amount of monoclonal antibody entering human body at time , and represents the attenuation coefficient of monoclonal antibody.
We scale those new parameters in model (32) as follows:and we choose the scaling . By replacing by , we obtain the following scaled model with treatment:with initial conditions
Here, denotes the dimensionless concentration of monoclonal antibody.
3.1. Model Analysis
Next, we discuss the effect of three types of immunotherapy on tumors such as Case I: ACI model () Case II: MAI model () Case III: combined immunotherapy model ()
We set in (34), and we have
3.1.1. ACI Model
When , we put and (36) becomes
Putting yields the tumor-free equilibrium, namely,where , . always exists.
Putting and eliminating in (37), we have Case (a): if (where , due to and , then (39) has one positive root in the interval . Therefore, system (34) has a unique coexistence equilibrium: We show the existence conditions of in the following: Case (1): . Since , we have . Since , we have the following two cases: Case where , i.e., . Since , we have . That is a contradiction. Case where , i.e., . Since , we have implying . If , i.e., , we have , there exists . Case (2): . Since , we have . Since , we have the following two cases: Case where , i.e., . Since , we have , there exists . Case where , i.e., . Since , we have implying . If , i.e., , we have , there exists . Case (b): If , then (39) can have two positive roots and in the interval , where and are the roots of (39). Therefore, system (34) can have two coexistence equilibria:
In the following, we show the existence condition of and .
Differentiating (39), we havewhere and . If and , then (43) has two positive roots and . And if , then (43) has one positive root . Case (1): . Since , we have . Since , we have the following two cases: Case where , i.e., . Since , we have . If , then there exist and . Case where , i.e., . Since , we have implying . Then, we have . If , then there exist and . Case (2): .
Since , we have . Since , we have the following two cases:
Case where , i.e., . Since , we have . That is a contradiction.
Case where , i.e., . Since , we have implying . Then, we have . If , then there exist and .
Theorem 5. System (34) with always has one tumor-free equilibrium , which is LAS if the inequality holds.
Proof. At , the characteristic equation becomesand one root of characteristic equation is . It is easily noted that as , so the solutions of have always negative real parts. Therefore, is LAS if and only if , that is, if and only if .