Mathematical Modelling of the Inhibitory Role of Regulatory T Cells in Tumor Immune Response

)e 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 Tcells’ (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. )is 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. )e numerical results show that ACI can control TCs more, while MAI can delay the inhibitory effect of Tregs on ECs. )e 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.


Introduction
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 [1]. 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 CD8 + T cells to differentiate into CTLs [4][5][6]. Recently, it has been reported that regulatory T cells (Tregs) can inhibit CTLs and promote the escape of TCs [7]. Tregs suppress immune cells, and when the war between T cells and infection is over, the Tregs signal to stop [8]. 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][10][11]. Its advantages are good destruction of tumor and persistence, while its disadvantages are serious toxic and side effects. e 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][15][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 [17]. Similar to CTLA-4, programmed death receptor 1 (PD-1) can also promote the activation and development of Tregs [18][19][20][21]. Blocking PD-1 can prevent the development of Tregs and prevent the conversion of HTCs into Tregs [22]. 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][24][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][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. e 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 [26]. In 1994, Kuznetsov et al. established the famous two-dimensional ODEs model, postulating that tumor growth follows the Logistic growth pattern. ey evaluated the parameters of the model by fitting experimental data from mice [27]. In 2003, Stolongo-Costa et al. assumed that TCs follows the exponential growth pattern and constructed a two-dimensional ODEs model. ey analyzed the basic properties of the model and provided conditions for stability of the tumor-free equilibrium, explaining its epidemiological significance [28]. In 2004, Galach simplified Kuznetsov's system to account for the effect of immune delay on the tumor immune system [29]. In 2014, Dong et al. constructed a three-dimensional ODEs model focusing on the effects of HTCs on the tumor immune system [4].
In 1998, Kirschner and Panetta generalized Kuznetsov-Taylor model and illustrated the dynamics between TCs, ECs, and IL-2. ey firstly introduced ACI into their model which can explain both short-term tumor oscillations in tumor sizes as well as long-term tumor relapse [11]. 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 Fcos 2 ωt and established a cycle immunotherapy model. ey obtained some thresholds of the frequency and intensity of immunotherapy [28]. In 2006, de Pillis et al. constructed the six-dimensional ODEs model to investigate the effects of combined chemotherapy and immunotherapy on tumor control. ey briefly analyzed the nature of the model and discussed the optimal treatment using optimal control theory [30]. In 2008, Bunimovich-Mendrazitsky et al. established a pulsed differential equation model with Bacillus Calmette-Guerin tumor immunotherapy. ey obtained the critical threshold and pulse frequency of BCG injection dose that could successfully treat superficial bladder cancer [31]. In 2012, Wilson and Levy established a mathematical model containing Tregs. ey 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 [8]. 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. e model was used to fit parameters with the experimental data. e results showed that increasing the resistance of PD-L1 doses can greatly improve the clearance rate of tumor [32].
is paper investigates the role of Tregs in the tumor immune system. erefore, we incorporate the fourth population of Tregs into the previous system in [8]. 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. en, we establish a four-dimensional ODEs model described as below: where T(t), E(t), H(t), and R(t) represent the populations of TCs, ECs, HTCs, and Tregs, respectively. e first equation describes the rate change of TCs population. e tumor follows logistic growth dynamics with growth rate a, and the maximum capacity is 1/b. n represents the loss rate of TCs by ECs interaction. e second equation describes the rate change of the ECs population. d 1 is the mortality rate of ECs. p is the activation rate of ECs by HTCs, and q is the inhibition rate of Tregs on ECs. e third equation describes the rate change of the HTCs population, s 2 is birth rate of HTCs produced in the bone marrow, and HTCs have a natural lifespan of an average 1/d 2 days. k 2 is HTCs stimulation rate by the presence of identified tumor antigens. e fourth equation gives the rate change of the Tregs population, r 1 and r 2 are the activation rates of Tregs by ECs and HTCs, respectively. d 3 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 T 0 � E 0 � H 0 � R 0 � 10 6 . By replacing τ by t, we obtain the following scaled model: with initial conditions Here, x, y, z, and u denote the dimensionless densities of TCs, ECs, HTCs, and Tregs populations, respectively.

Well Posedness of Model (3).
e following proposition establishes the well posedness of model (3) with initial conditions (4).
If Δ < 0, then (10) does not have a positive root. Hence, P * does not exist in this case.
If Δ > 0, then (10) has two positive roots, where x * − . e existence conditions for each equilibrium are given in Table 1.

Stability of Equilibria.
In order to investigate the local stability of the above equilibria P 0 , P 1 , P 2 , and P * of system (3), we linearize the system and obtain Jacobian matrix at each equilibrium P(x, y, z, u): e corresponding characteristic equation is Theorem 1. System (3) always has one tumor-free and ECsfree equilibrium P 0 , which is unstable.
We summarize the above results in Table 1

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 [45]. We select the following parameters' set [4]: Note that ω * 2 � βδ 2 � 0.00011. By calculations, we have Table 1: e existence and stability conditions of each equilibrium.

Existence conditions Stability conditions
Complexity And we find the stability region of P * � (x * , y * , z * , u * ) (see Figure 2). P * is LAS in regions I and II. P * is unstable in region III.
Case (a): we choose a point Q 1 � (0.48, 0.0001) in the region I; then, system (3) has one interior equilibrium: e eigenvalues of Jacobian matrix of (12) are so P * 1 is stable, as shown in Figure 3 e eigenvalues of Jacobian matrix of (12) are so P * 2 is stable, as shown in Figure 3 so P * 3 is stable, as shown in Figure 3 so P * 4 is unstable, as shown in Figure 3(d).
Below, we perform numerically bifurcation analysis of x against θ for different values of ω 2 .

Proposition 2.
When 0 < θ < θ 0 , P 0 exists and is unstable, P 1 exists and is LAS, P 2 exists and is unstable, and P * does not exist. When θ 0 < θ < θ 1 , P 0 exists and is unstable, P 1 exists and is unstable, P 2 exists and is unstable, and P * exists and is LAS. When θ 1 < θ < θ 2 , P 0 exists and is unstable, P 1 does not exist, P 2 exists and is unstable, and P * exists and is LAS. When θ > θ 2 , P 0 exists and is unstable, P 1 does not exist, P 2 exists and is LAS, and P * does not exist. Figure 4(b)).
We obtain the following result. Proposition 3. When 0 < θ < θ 0 , P 0 exists and is unstable, P 1 exists and is LAS, and P * does not exist. When θ 0 < θ < θ 4 , where θ 4 � 0.4657 is a Hopf bifurcation, P 0 exists and is unstable, P 1 exists and is LAS, and P * exists and is LAS. When θ 4 < θ < θ 1 , P 0 exists and is unstable, P 1 exists and is unstable, and P * exists and is unstable. When θ 1 < θ < θ 3 , P 0 exists and is unstable, P 1 does not exist, and P * exists and is unstable. When θ > θ 3 , P 0 exists and is unstable, P 1 does not exist, and P * does not exist.
Consider the case where HTCs stimulation rate (ω 2 ) is low by the presence of identified tumor antigens (ω 2 < ω * 2 ). When the inhibition rate of Tregs to ECs θ is lower than θ 0 , the solution of system (3) approaches P 1 implying that ECs can still effectively remove TCs. When θ 0 < θ < θ 2 , the solution of system (3) approaches P * showing the coexistence of TCs and immune cells, which means that the patient can survive with tumors. When θ > θ 2 , TCs escape the control of the immune system and develop into malignant tumors.
Next consider the case where HTCs stimulation rate (ω 2 ) is high by the presence of identified tumor antigens (ω 2 > ω * 2 ). When θ < θ 0 , the solution of system (3) approaches P 1 implying that TCs can be effectively removed by ECs. When θ 0 < θ < θ 4 , the solution of system (3) approaches P * implying that TCs can still be controlled by the immune system. When θ > θ 4 , 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.

Treatment Model
In order to investigate well the effect of Tregs in tumor immune response under the treatment, we follow the way in [11] to introduce the constant treatment term s 1 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 [12], the effect of CTLA4 or PD-1 on Tregs can be considered. en, we shall establish a five-dimensional ODEs model described as below: where W(t) represents the concentration of monoclonal antibody in human body at time t, s 1 represents the treatment term of introducing LAK and TIL into the region of tumor localization, m represents the inhibition rate of monoclonal antibody on Tregs, s 3 represents the amount of monoclonal antibody entering human body at time t, and d 4 represents the attenuation coefficient of monoclonal antibody.
We scale those new parameters in model (32) as follows: and we choose the scaling W 0 � 10 6 . By replacing τ by t, we obtain the following scaled model with treatment: with initial conditions Here, v denotes the dimensionless concentration of monoclonal antibody.

Model Analysis.
By a similar proof of Proposition 1, we can obtain the well posedness of model (34) with initial conditions (35) as follows.

MAI Model.
When σ 1 � 0, σ 3 > 0, (36) becomes Equation (49) can be simplified as follows: where δ 31 � δ 3 δ 4 + ξσ 3 /δ 4 . e analysis of (50) is similar to that of (12), and the equilibria of (34) can be obtained as follows: Here, x * ∈ (0, μ) and satisfies the equation e Jacobian matrix of system (34) at any equilibrium F(x, y, z, u, v) is as follows: e corresponding characteristic equation is Substituting δ 31 for δ 3 in (12), we have us, the stability analysis of F is similar to that of P. Substituting δ 31 for δ 3 in (20) to get new conditions for the stability of F * : where We set Complexity and we obtain the following results.
For system (34) at N * � (x * , y * , z * , u * , v * ), Jacobian matrix is given as follows: e corresponding characteristic equation is us, the stability analysis of N * is similar to that of E * . Set θ 9 � ρδ 31 /c 2 and substitute δ 31 for δ 3 in (47) to get new conditions: where erefore, we can obtain the following result.

Numerical Simulations
. We conduct numerical simulations of the ACI model. We choose ω 2 � 0.0004 and the other parameter values are given in (22). We find the stability region of E * � (x * , y * , z * , u * ), as shown in Figure 5(a). E * is stable in region I and unstable in region II. Case so E * is stable, as shown in Figure 6 so E * is unstable, as shown in Figures 6(c) and 6(d).
By comparing the curves in Figure 5(d), we see that, as σ 1 increases gradually, the stable region of tumor-free equilibrium E 0 of system (34) gradually increases (the intersection point of the curve and the x-coordinate gradually moves to the right) and the stability region of equilibrium E * in system (34) gradually increases (the blue curve gradually moves upward).
is indicates that increasing the injection volume of adoptive T cells can not only delay the inhibitory effect of Tregs on tumor immune response but also help the immune system to remove more TCs. It also helps the immune system to control more TCs, keeping them at a stable state.
By comparing the curves in Figure 7(d), we can find that, with the gradual increase of σ 3 , the stability region of the tumour-free equilibrium F 1 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right). is shows that the increase of antibody injection quantity can help to slow down Tregs inhibition of tumor immune responses. e stability region of equilibrium F * of system (34) gradually increases (the blue curve gradually moves upward), which means that increasing the amount of antibody injected can help the immune system to control more TCs. By comparing Figures 5 and 7, we find that the effect of MAI is better than that of AIC in delaying the inhibitory effect of Tregs on tumor immune response (at the same injection dose, the intersection point of curves and x-coordinate in Figure 7 moves to the right more widely than that in Figure 5). AIC is more effective than MAI in controlling TCs (the blue  16 Complexity curve in Figure 7 moves up more than the blue curves in Figure 5 at the same injection dose).

Combined Immunotherapy Model
We conduct numerical simulations of the combined immunotherapy model. We choose ω 2 � 0.0004, ξ � 0.05, and δ 4 � 0.25, and the other parameter values are given in (22). We find the stability region of N * � (x * , y * , z * , u * , v * ), as shown in Figure 9(a). When σ 3 increases gradually, the stability region of N * also increases gradually.
(91) e eigenvalues of Jacobian matrix of (64) are so N 0 is stable, as shown in Figure 10 so N 0 is unstable. We choose combined immunotherapy parameter σ 1 � 0.4 and σ 3 � 0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 9(c)). By some calculations, we have θ 8 � αδ 31  . e bifurcation diagram of x with respect to θ (the stable steady state is represented by the blue curve, and the unstable one corresponds to the red curve). σ 1 δ 2 δ 31 /α(αc 1 δ 2 + σ 2 c 2 ) � 0.646, and θ 9 � ρδ 31 / c 2 � 0.791. And Hopf bifurcation point appears at θ * CI � 0.727. We can obtain the following result. Proposition 7. When 0 < θ < θ 8 , N 0 exists and is LAS, and N * is nonexistent. When θ 8 < θ < θ * CI , N 0 exists and is unstable, and N * exists and is LAS. When θ * CI < θ < θ 9 , N 0 exists and is unstable, and N * exists and is unstable. When θ > θ 9 , N 0 exists and is unstable, and N * is nonexistent.
For combined immunotherapy, by comparing the curves in Figure 9(d), we know that, with the increase of σ 1 and σ 3 , the stability region of the tumour-free equilibrium of N 0 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right), and the stability region of equilibrium N * of system (34) gradually increases (the blue curve gradually moves upward). By comparing Figure 9 with Figures 5 and 7, it can be seen that combined immunotherapy has better effects on delaying the inhibitory effect of Tregs on tumor immune response and helps the immune system to control more TCs than ACI or MAI.

Discussion and Conclusion
Tregs-mediated tumor immune escape is one of the core mechanisms of tumor immune regulation. And Tregs have been found to mediate tumor evasion and immune escape in many different solid tumors [46]. e study on Tregs has a very high research value and application prospect in the immunotherapy of tumors. If the activity of Tregs is controlled or blocked during the tumor immune response, or a barrier is set to prevent Tregs from migrating into the tumor microenvironment, then the effect of tumor immunotherapy can be improved [47]. First, we developed a mathematical model to study the inhibitory role of Tregs in the tumor immune system. For the lower recognition of tumor antigens by the immune system, the stronger the inhibition effect of Tregs on ECs, TCs can easily escape the control of the immune system (see Figure 4(a)). When the immune system is highly sensitive to tumor antigens, the immune system activates ECs; the stronger the inhibition effect of Tregs on ECs, the more complicated interactions between TCs and immune cells (see Figure 4(b)).
Second, we incorporated the previous mathematical model with three types of immunotherapy to obtain ACI model, MAI model, and combined both ACI and MAI model. rough the theoretical analysis and numerical simulations, we found that ACI can control more TCs, but have no obvious effects on reducing the inhibitory effect of Tregs on ECs (see Figure 5). MAI can effectively reduce the inhibitory effect of Tregs on ECs, but cannot control more TCs (see Figure 7). However, combination immunotherapy with ACI and MAI is more effective than single immunotherapy. It can not only significantly reduce the inhibitory effect of Tregs on ECs but also help the immune system to kill TCs to the maximum extent (see Figure 9). erefore, we recommend the use of combined immunotherapy in the treatment of tumors. Besides, clinical trials are needed to further evaluate the safety and efficacy of combined immunotherapy.
is paper focused on the general process of tumor immune response with negative feedback. Using the mathematical model, it is possible to simulate the state of tumors in the immune system at different inhibition states. e results of the study can contribute to the understanding of tumor immunity; at the same time, it also provides new ideas for the treatment of tumors. However, due to the complexity and heterogeneity of tumor microenvironment, there is still a certain gap between the mathematical model and the description of the real interactions between the tumor and immune system. erefore, specific tumor microenvironment and heterogeneous tumoral populations should be considered in practical application to make the model more realistic [34]. In addition to the bilinear incidence model considered in this paper, nonbilinear model with saturation incidence should be employed in the further study [35]. Besides, other important factors such as immune activation delays [36][37][38]48], stochastic effects [39,40,49], and impulsive perturbations [41,42,50] can be considered in the modelling of the tumor immune system. e tumor immune response dynamics in vivo is very complex and not well understood primarily because the measurements of the necessary parameters are difficult in vivo [43].

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.