Next Article in Journal
Timelike Constant Axis Ruled Surface Family in Minkowski 3-Space
Next Article in Special Issue
Dynamics and Optimal Harvesting for Fishery Models with Reserved Areas
Previous Article in Journal
Solitons of η-Ricci–Bourguignon Type on Submanifolds in (LCS)m Manifolds
Previous Article in Special Issue
Delay Effects on Plant Stability and Symmetry-Breaking Pattern Formation in a Klausmeier-Gray-Scott Model of Semiarid Vegetation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling the Nonmonotonic Immune Response in a Tumor–Immune System Interaction

1
School of Mathematics and Statistics, Central China Normal University, Wuhan 430079, China
2
School of Information Technology, Guangxi Police College, Nanning 530028, China
3
School of Public Health, Nanjing Medical University, Nanjing 211166, China
4
College of Science and Engineering, Aoyama Gakuin University, Kanagawa 252-5258, Japan
5
Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur 208016, India
6
Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan 430079, China
*
Authors to whom correspondence should be addressed.
Symmetry 2024, 16(6), 676; https://doi.org/10.3390/sym16060676
Submission received: 5 April 2024 / Revised: 6 May 2024 / Accepted: 14 May 2024 / Published: 31 May 2024
(This article belongs to the Special Issue Mathematical Modeling in Biology and Life Sciences)

Abstract

:
Tumor–immune system interactions are very complicated, being highly nonlinear and not well understood. A large number of tumors can potentially weaken the immune system through various mechanisms such as secreting cytokines that suppress the immune response. In this paper, we propose a tumor–immune system interaction model with a nonmonotonic immune response function and adoptive cellular immunotherapy (ACI). The model has a tumor-free equilibrium and at most three tumor-presence equilibria (low, moderate and high ones). The stability of all equilibria is studied by analyzing their characteristic equations. The consideration of nonmonotonic immune response results in a series of bifurcations such as the saddle-node bifurcation, transcritical bifurcation, Hopf bifurcation and Bogdanov–Takens bifurcation. In addition, numerical simulation results show the coexistence of periodic orbits and homoclinic orbits. Interestingly, along with various bifurcations, we also found two bistable scenarios: the coexistence of a stable tumor-free as well as a high-tumor-presence equilibrium and the coexistence of a stable-low as well as a high-tumor-presence equilibrium, which can show symmetric and antisymmetric properties in a range of model parameters and initial cell concentrations. The new findings indicate that under ACI, patients can possibly reach either a stable tumor-free state or a low-tumor-presence state in the presence of nonmonotonic immune response once the immune system is activated.

1. Introduction

The immune system can recognize and control nascent tumor cells through a process called immunosurveillance [1]. This involves main antitumor immune cells such as dendritic cells, T cells, natural killer (NK) cells and tumor-associated macrophages (TAMs) to identify and eliminate abnormal or cancerous cells [2]. However, tumors can evade attacks from the immune system by various mechanisms, such as restricting antigen recognition, orchestrating an immunosuppressive microenvironment, inducing T cell exhaustion and so on [3,4]. Understanding tumor–immune system interactions is essential for developing effective immunotherapies to treat tumors. Immunotherapy is a transformative strategy that activates the body’s immune response against a tumor [5]. It can include treatments like adoptive cellular immunotherapy (ACI), oncolytic virus therapy, cancer vaccines, immune checkpoint therapy and combination therapies [6,7].
A variety of mathematical models have been proposed and analyzed to explore the dynamics of tumor–immune system interactions and have evolved to capture much more complex aspects of the immune response [8], such as healthy cells [9], natural killer cells [10], helper T cells [11], regulatory T cells [12], macrophages [13] and so on. For mathematical simplicity, a bilinear form has mainly been used to describe the growth of the immune response to tumors [14]. However, tumor–immune system interactions are very complicated, being highly nonlinear, and evolve with time to a certain extent [15]. So it is necessary to model tumor–immune system interactions using nonlinear response functions, which are widely used in the modeling of tumor–immune interactions [16], the ecological system [17], virus dynamics [18], the metabolic system [19] and so on.
Kuznetsov et al. [20] developed a classical model describing the effector cell response to the growth of an immunogenic tumor, where the rate of stimulated accumulation of effector cells attains a maximum value as the number of T cells becomes large. Kirschner et al. [21] introduced cytokine interleukin-2 (IL-2) and ACI into the tumor–immune interaction model with the Michaelis–Menten form to indicate the saturated effects of the immune response. Their model can explain both short-term oscillations in tumor sizes and long-term tumor relapse. Dritschel et al. [22] considered adding a helper T cell compartment to Kuznetsov’s model to distinguish immune cell promotion and tumor cell killing, and they implicitly incorporated tumor-modulated immunosupression by altering the functions describing immune recruitment and proliferation. Their model can exhibit the three phases (elimination, equilibrium and escape) of cancer immunoediting, periodic growth and suppression of the immune and tumor populations. Talkington [23] explored and compared four previous models with adoptive immunotherapy that predict successful tumor elimination. Song et al. [24] developed a tumor growth model with immune system response and drug therapy to study the mechanisms of chemotherapy. Bekker et al. [25] carried out a survey of mathematical models to evaluate the effects of radiotherapy on tumor–immune system interactions. Butner et al. [26] presented a comprehensive review of mathematical models on cancer immunotherapy for clinical applications. Tang et al. [27] developed a mathematical model to investigate the interactions between effector cells, tumor cells and cytokines with pulsed radio/chemo-immunotherapy. Recently, Zhang et al. [28] investigated a tumor–immune system interaction model with saturated activation of T cells and dendritic cell therapy showing rich dynamics.
We investigate the effect of nonlinear interactions between tumor and immune cells in the tumor microenvironment. Inspired by previous studies, we propose the following ordinary differential equations model with nonlinear tumor–immune interactions:
d T d t = a T ( 1 b T ) n c E T T 2 + c 2 , d E d t = s 1 d 1 E + p E H , d H d t = s 2 d 2 H + k T H ,
where T ( t ) , E ( t ) and H ( t ) represent the populations of tumor cells (TCs), effector cells (ECs) and helper T cells (HTCs), respectively. Here, both ECs and HTCs correspond to the immune system/immune response. The growth of TCs is assumed to follow the logistic form a T ( 1 b T ) , where a is the maximal growth rate and b 1 is the maximal carrying capacity of the biological environment for TCs. In the initial stage, the immune system recognizes tumor antigens, and a large number of immune cells are recruited into the tumor microenvironment. Over time, tumors can generate various mechanisms to escape immune surveillance, such as recruiting immunosuppressive cells like regulatory T cells and secreting immunosuppressive cytokines like IL-6 [29]. Furthermore, the presence of a large number of tumor cells can suppress and weaken the immune system, allowing tumors to evade immune surveillance and continue to grow [30,31]. Therefore, we assume that the inhibition of TCs by ECs follows biphasic dependence on the number of TCs using a nonmonotonic response function F ( T ) = n c T T 2 + c 2 , where F ( T ) increases first to the peak and then decreases to zero as T is sufficiently large [32]. n represents the maximal inhibition rate of TCs by ECs, and c is the number of TCs at which ECs’ inhibition is half-maximal. The inhibition rate of ECs on TCs attains a maximum value of n / 2 when T = c and decays to zero as T + . Here, we consider a weak proliferation of ECs stimulating TCs as p E H , while ignoring the ECs’ stimulation from TCs [33]. p is the activation rate of ECs by HTCs, and k is the HTCs’ stimulation rate by the presence of identified tumor antigens. s 1 refers to the constant inflow of ACI, and s 2 is the birth rate of HTCs produced in the bone marrow. d 1 and d 2 are the natural death rates of ECs and HTCs, respectively. All model parameters are positive and described in Table 1.
As far as we know, this is the first attempt to model the immune response against a tumor by a nonmonotonic response function, which can capture long-term tumor–immune system interactions. In this study, we aim to exploit the rich dynamics of model (1), with nonmonotonic immune response function and ACI, and explain the biological implications of dynamical properties. The rest of this paper is organized as follows. In Section 2, we firstly nondimensionalize system (1) by scaling the parameters. Then we study the existence and stability of possible equilibria of the dimensionless system. In Section 3, we carry out a bifurcation analysis including saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation. In Section 4, we investigate the effects of some parameters on the stability of each equilibrium and perform the numerical simulations to show several types of bifurcations and bistability phenomena. In the last section, we give a summary and discussion.

2. Stability Analysis

In this section, we firstly nondimensionalize system (1), and then we study the existence and stability of possible equilibria of the dimensionless system.

2.1. Dimensionless Model

The first step is to nondimensionalize system (1). Let
t ^ = t d 1 , x = b T , y = b E , z = p H d 1 .
Then we obtain
d x d t ^ = a ^ x ( 1 x ) n ^ x y x 2 + c ^ 2 , d y d t ^ = s 1 ^ y + y z , d z d t ^ = s 2 ^ d 2 ^ z + k ^ x z ,
where
a ^ = a d 1 , n ^ = n c b d 1 , c ^ = b c , s 1 ^ = s 1 b d 1 , s 2 ^ = s 2 p d 1 2 , d 2 ^ = d 2 d 1 , k ^ = k b d 1 .
Without any loss of generality, we can ignore the hat and write the dimensionless model as follows:
d x d t = a x ( 1 x ) n x y x 2 + c 2 , d y d t = s 1 y + y z , d z d t = s 2 d 2 z + k x z ,
with initial conditions x ( 0 ) 0 , y ( 0 ) > 0 , z ( 0 ) > 0 , where x ( t ) , y ( t ) and z ( t ) denote the dimensionless densities of TCs, ECs and HTCs.

2.2. Existence of Equilibria

In this section, we study the existence of the equilibrium points of system (3). We set d x d t = 0 , d y d t = 0 and d z d t = 0 in system (3), and we have
0 = a x ( 1 x ) n x y x 2 + c 2 , 0 = s 1 y + y z , 0 = s 2 d 2 z + k x z .
Assuming x = 0 yields the tumor-free equilibrium, namely,
E 0 = ( x 0 , y 0 , z 0 ) = 0 , s 1 d 2 d 2 s 2 , s 2 d 2 ,
which exists when s 2 < d 2 .
Putting x 0 and eliminating z in (4), we find following two equations:
y = a ( 1 x ) ( x 2 + c 2 ) n f ( x ) , y = s 1 1 s 2 d 2 k x = s 1 ( d 2 k x ) d 2 s 2 k x g ( x ) .
If the tumor-presence equilibrium E = ( x , y , z ) = ( x , s 1 ( d 2 k x ) d 2 s 2 k x , s 2 d 2 k x ) exists, then x is the positive root of the equation f ( x ) = g ( x ) and x should satisfy 0 < x < min { 1 , d 2 s 2 k } x ˜ and s 2 < d 2 , since y and z are both positive. It is straightforward to show that these steady-state solutions lie at the intersection of the two curves y = f ( x ) and y = g ( x ) for x > 0 only.
We differentiate f ( x ) and g ( x ) with respect to x in (5) as follows:
f ( x ) = a n ( 3 x 2 + 2 x c 2 ) , g ( x ) = s 1 s 2 k ( d 2 s 2 k x ) 2 .
From the expression of f ( x ) in (5), we have f ( 0 ) = a c 2 n > 0 , f ( x ) 0 when x 1 , and f ( x ) < 0 when x > 1 . From the expression of f ( x ) in (6), we can obtain the following cases.
Case (a): When c 3 3 , f ( x ) decreases in the range x ( 0 , 1 ) (see Figure 1a).
Case (b): When 0 < c < 3 3 , f ( x ) has two zero points x 1 = 1 1 3 c 2 3 ( 0 < x 1 < 1 3 ) and x 2 = 1 + 1 3 c 2 3 ( 1 3 < x 2 < 2 3 ) , where 1 3 is the inflection point. So f ( x ) decreases when x ( 0 , x 1 ) or ( x 2 , 1 ) and increases when x ( x 1 , x 2 ) (see Figure 1b). By direct calculations, we have
f ( x 1 ) = 2 a 27 n ( 2 + 1 3 c 2 ) ( 1 + 3 c 2 1 3 c 2 ) 0 , 8 a 27 n , f ( x 2 ) = 2 a 27 n ( 2 1 3 c 2 ) ( 1 + 3 c 2 + 1 3 c 2 ) 4 a 27 n , 8 a 27 n .
Figure 1. Two possible shapes of f ( x ) . We fix a = 10 , n = 3 , s 1 = 0.31 , s 2 = 0.5 , d 2 = 1 , k = 0.19 .
Figure 1. Two possible shapes of f ( x ) . We fix a = 10 , n = 3 , s 1 = 0.31 , s 2 = 0.5 , d 2 = 1 , k = 0.19 .
Symmetry 16 00676 g001
From the expression of g ( x ) in (5), we have g ( 0 ) = s 1 d 2 d 2 s 2 > 0 when d 2 > s 2 . The asymptotic line of g ( x ) is x = d 2 s 2 k x 0 . As x , g ( x ) s 1 . From the expression of g ( x ) in (6), we have g ( x ) > 0 and g ( x ) = 2 s 1 s 2 k 2 ( d 2 s 2 k x ) 3 > 0 for 0 < x < x ˜ . So f ( x ) is a convex function in the interval ( 0 , x ˜ ) .
The above analysis can be summarized into the following results:
Case (a): c 3 3 .
A1. If f ( 0 ) g ( 0 ) , then there exists one intersection point of f ( x ) and g ( x ) .
A2. If f ( 0 ) < g ( 0 ) , then there is no intersection point of f ( x ) and g ( x ) .
Case (b): 0 < c < 3 3 .
In order to find the intersection point of f ( x ) and g ( x ) , we study the critical point in which f ( x ) and g ( x ) are tangent. We set
f ( x ) = g ( x ) , i . e . , a n ( 3 x 2 + 2 x c 2 ) = s 1 s 2 k ( d 2 s 2 k x ) 2 ,
and
f ( x ) = g ( x ) , i . e . , a ( 1 x ) ( x 2 + c 2 ) n = s 1 ( d 2 k x ) d 2 s 2 k x .
We suppose that x + is a root of the equation in (7), and we show the existence of x + as follows. Firstly, we study the interaction point of f ( x ) and g ( x ) .
Case (1): If f ( 1 3 ) < g ( 1 3 ) , i.e., a n ( 1 3 c 2 ) < s 1 s 2 k ( d 2 s 2 1 3 k ) 2 , then there is no intersection point of f ( x ) and g ( x ) (see Figure 2a).
Case (2): If f ( 1 3 ) = g ( 1 3 ) , i.e., a n ( 1 3 c 2 ) = s 1 s 2 k ( d 2 s 2 1 3 k ) 2 , then there exists one intersection point of f ( x ) and g ( x ) (see Figure 2b).
Case (3): If f ( 1 3 ) > g ( 1 3 ) , i.e., a n ( 1 3 c 2 ) > s 1 s 2 k ( d 2 s 2 1 3 k ) 2 , then there exist two intersection points of f ( x ) and g ( x ) (see Figure 2c). We suppose that these two intersection points are x 1 + and x 2 + ( x 1 , 2 + ). Then we have x 1 < x 1 + < 1 3 < x 2 + < x 2 . When x is a root of both Equations (7) and (8), f ( x ) and g ( x ) will be tangent.
Figure 2. The relative positions of f ( x ) and g ( x ) depending on k. We fix a = 10 , n = 3 ,   s 1 = 0.31 , s 2 = 0.5 , d 2 = 1 , c = 0.5 .
Figure 2. The relative positions of f ( x ) and g ( x ) depending on k. We fix a = 10 , n = 3 ,   s 1 = 0.31 , s 2 = 0.5 , d 2 = 1 , c = 0.5 .
Symmetry 16 00676 g002
Next, we discuss the intersection point of f ( x ) and g ( x ) as follows.
B1. When f ( 0 ) g ( 0 ) , i.e., a c 2 n s 1 d 2 d 2 s 2 , there exists at least one intersection point (since when x d 2 s 2 k , g ( x ) + ).
Case B1.1. When f ( 1 3 ) < g ( 1 3 ) , then f ( x ) and g ( x ) only have one intersection point. In fact, if f ( x ) and g ( x ) have two intersection points x 3 and x 4 (that is, F ( x 3 ) = F ( x 4 ) where F ( x ) = f ( x ) g ( x ) ), then by intermediate value theorem, there exists η ( x 3 , x 4 ) , such that F ( η ) = 0 , i.e., f ( η ) = g ( η ) . It is a contradiction.
Case B1.2. When f ( 1 3 ) = g ( 1 3 ) , there exists at least one intersection point.
Case B1.3. When f ( 1 3 ) > g ( 1 3 ) , note that x 1 , 2 + satisfies (7), and we have
f ( x 1 , 2 + ) = a ( 1 x 1 , 2 + ) ( ( x 1 , 2 + ) 2 + c 2 ) n and g ( x 1 , 2 + ) = s 1 ( d 2 k x 1 , 2 + ) d 2 s 2 k x 1 , 2 + , then we have the following cases:
(a)
If f ( x 1 + ) > g ( x 1 + ) , there exists one intersection point (see Figure 3a).
(b)
If f ( x 1 + ) = g ( x 1 + ) , there exist two intersection points (see Figure 3b).
(c)
If f ( x 1 + ) < g ( x 1 + ) , f ( x 2 + ) > g ( x 2 + ) , there exist three intersection points (see Figure 3c).
(d)
If f ( x 2 + ) = g ( x 2 + ) , there exist two intersection points (see Figure 3d).
(e)
If f ( x 2 + ) < g ( x 2 + ) , there exists one intersection point (see Figure 3e).
B2. When f ( 0 ) < g ( 0 ) , i.e., a c 2 n < s 1 d 2 d 2 s 2 , we have the following cases:
Case B2.1. When f ( 1 3 ) g ( 1 3 ) , i.e., a n ( 1 3 c 2 ) < s 1 s 2 k ( d 2 s 2 1 3 k ) 2 , there is no intersection point (see Figure 3f).
Case B2.2. When f ( 1 3 ) > g ( 1 3 ) , i.e., a n ( 1 3 c 2 ) > s 1 s 2 k ( d 2 s 2 1 3 k ) 2 (we know f ( x 1 ) < g ( x 1 ) , f ( x 1 + ) < g ( x 1 + ) ).
(a)
If f ( x 2 + ) < g ( x 2 + ) , there is no intersection point (see Figure 3g).
(b)
If f ( x 2 + ) = g ( x 2 + ) , there exists one intersection point (see Figure 3h).
(c)
If f ( x 2 + ) > g ( x 2 + ) , there exist two intersection points (see Figure 3i).
Figure 3. The interactions of f ( x ) and g ( x ) , illustrating the possible number of equilibrium points. We choose a = 10 , n = 3 , s 1 = 0.23 , s 2 = 0.5 , d 2 = 1 .
Figure 3. The interactions of f ( x ) and g ( x ) , illustrating the possible number of equilibrium points. We choose a = 10 , n = 3 , s 1 = 0.23 , s 2 = 0.5 , d 2 = 1 .
Symmetry 16 00676 g003

2.3. Stability of Equilibria

In order to investigate the local stability of the equilibria of system (3), we linearize the system and obtain the characteristic equation evaluated at a typical equilibrium E ¯ = ( x ¯ , y ¯ , z ¯ ) .
det a 2 a x ¯ n y ¯ ( c 2 x ¯ 2 ) ( x ¯ 2 + c 2 ) 2 λ n x ¯ x ¯ 2 + c 2 0 0 z ¯ 1 λ y ¯ k z ¯ 0 k x ¯ d 2 λ = 0 .
At E 0 , characteristic Equation (9) becomes
λ a n y 0 c 2 λ s 2 d 2 1 ( λ + d 2 ) = 0 .
Three eigenvalues are λ 1 = a n y 0 c 2 , λ 2 = s 2 d 2 1 < 0 , λ 3 = d 2 < 0 . Note that the feasibility of E 0 is s 2 < d 2 . It is easy to show that λ 1 < 0 if a < n s 1 d 2 ( d 2 s 2 ) c 2 , i . e . , c < n s 1 d 2 a ( d 2 s 2 ) . Therefore, we have the following result.
Theorem 1.
System (3) has one tumor-free equilibrium E 0 , which is locally asymptotically stable (LAS) if s 2 < d 2 and a < n s 1 d 2 ( d 2 s 2 ) c 2 .
Note that n y = a ( 1 x ) ( x 2 + c 2 ) , y = s 1 1 z , s 2 z = d 2 k x for E = ( x , y , z ) , so the feasibility of equilibrium components requires that x < 1 , z < 1 .
At E , characteristic Equation (9) becomes
λ 3 + b 1 λ 2 + b 2 λ + b 3 = 0 ,
where
b 1 = a x ( 3 x 2 + 2 x c 2 ) x 2 + c 2 + 1 z + s 2 z , b 2 = a x ( 3 x 2 + 2 x c 2 ) x 2 + c 2 ( 1 z + s 2 z ) + ( 1 z ) s 2 z , b 3 = a x ( 3 x 2 + 2 x c 2 ) x 2 + c 2 ( 1 z ) s 2 z + k z n x y x 2 + c 2 .
By some algebraic calculations, we find
b 1 > 0 1 z + s 2 z > a x ( 3 x 2 + 2 x c 2 ) x 2 + c 2 ,
b 3 > 0 k z n y > a ( 3 x 2 + 2 x c 2 ) ( 1 z ) s 2 z ,
and
b 1 b 2 b 3 > 0 a 2 x 2 ( 3 x 2 + 2 x c 2 ) 2 ( x 2 + c 2 ) 2 ( 1 z + s 2 z ) a x ( 3 x 2 + 2 x c 2 ) x 2 + c 2 ( 1 z + s 2 z ) 2 + ( 1 z + s 2 z ) ( 1 z ) s 2 z k z n x y x 2 + c 2 > 0 .
By Routh–Hurwitz criteria, E is stable if the conditions (12)–(14) hold. Therefore, we have the following results.
Theorem 2.
System (3) has at most three tumor-presence equilibria E , which are LAS if the conditions (12)–(14) hold.

3. Local Bifurcation Analysis

In this section, we carry out a bifurcation analysis including saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation.

3.1. Saddle-Node Bifurcation

Consider the dimensionless model (3). y- and z-nullclines are given by
z = s 2 g 1 ( x ) , y = s 1 g 1 ( x ) g 1 ( x ) s 2 ,
where g 1 ( x ) = d 2 k x . Then, the nontrivial x-nullcline is given by
y = a n ( 1 x ) f 1 ( x ) ,
where f 1 ( x ) = x 2 + c 2 . If E ( x , y , z ) denotes a typical coexistence equilibrium point, then x is a positive root of the equation
Φ ( x ) a ( 1 x ) n s 1 g 1 ( x ) f 1 ( x ) ( g 1 ( x ) s 2 ) = 0 .
Let us assume that the system undergoes a saddle-node bifurcation at E , then x should be a double root of the equation Φ ( x ) = 0 , which requires
Φ ( x ) = 0 , Φ ( x ) = 0 .
We can prove that the Jacobian matrix of system (3) evaluated at E should have a simple zero eigenvalue when the conditions in (16) are satisfied. The Jacobian matrix for system (3) evaluated at E is given by
J = x a + n y f 1 ( x ) f 1 2 ( x ) n x f 1 ( x ) 0 0 z 1 y z g 1 ( x ) 0 g 1 ( x ) ,
where
y = s 1 g 1 ( x ) g 1 ( x ) s 2 , z = s 2 g 1 ( x ) .
Now we calculate the determinant of J as follows:
det ( J ) = x a + n y f 1 ( x ) f 1 2 ( x ) ( 1 z ) g 1 ( x ) + n x f 1 ( x ) y z g 1 ( x ) = x a + n s 1 g 1 ( x ) f 1 ( x ) f 1 2 ( x ) ( g 1 ( x ) s 2 ) ( g 1 ( x ) s 2 ) + n x s 1 s 2 g 1 ( x ) ( g 1 ( x ) s 2 ) .
And we obtain from Φ ( x ) = 0 that
a + n s 1 g 1 ( x ) f 1 ( x ) f 1 2 ( x ) ( g 1 ( x ) s 2 ) = n s 1 s 2 g 1 ( x ) ( g 1 ( x ) s 2 ) 2 .
Hence, we can find
det ( J ) = x n s 1 s 2 g 1 ( x ) ( g 1 ( x ) s 2 ) 2 ( g 1 ( x ) s 2 ) + n x s 1 s 2 g 1 ( x ) ( g 1 ( x ) s 2 ) = 0 .
Next, we have to verify the transversality conditions. The eigenvectors associated with the zero eigenvalue of the matrices J and ( J ) T are given by
V = 1 s 1 s 2 g 1 ( x ) ( g 1 ( x ) s 2 ) 2 s 2 g 1 ( x ) g 1 2 ( x ) , W = 1 n x g 1 ( x ) f 1 ( x ) ( g 1 ( x ) s 2 ) n s 1 x g 1 ( x ) f 1 ( x ) ( g 1 ( x ) s 2 ) 2 .
Here, we assume k as the bifurcation parameter. Then k S N is a threshold of k such that the condition (16) is satisfied and E is evaluated at k S N . Following the same notations as in [34], we have to check the transversality conditions.
First, we calculate
W T F k ( E ; k S N ) = W T 0 0 z x = n s 1 s 2 ( x ) 2 f 1 ( x ) ( g 1 ( x ) s 2 ) 2 0 ,
where
F x a ( 1 x ) n y x 2 + c 2 s 1 y + y z s 2 z ( d 2 k x ) ,
and F k is the partial derivative of F with respect to k.
After some lengthy algebraic calculation, we can verify that
W T D 2 F ( E ; k S N ) ( V , V ) 0 .
Hence, the saddle-node bifurcation is nondegenerate. Therefore, we have the following result.
Theorem 3.
System (3) undergoes a saddle-node bifurcation at E for k = k S N , where k S N satisfies the condition (16).

3.2. Transcritical Bifurcation

Consider the dimensionless model (3). The axial equilibrium point is given by E 0 = 0 , s 1 d 2 d 2 s 2 , s 2 d 2 , and the Jacobian matrix of system (3) evaluated at E 0 is given by
J 0 = a n s 1 d 2 ( d 2 s 2 ) c 2 0 0 0 s 2 d 2 1 s 1 d 2 d 2 s 2 k s 2 d 2 0 d 2 .
E 0 is feasible for d 2 > s 2 , and hence, two eigenvalues of J 0 are λ 2 = s 2 d 2 1 < 0 and λ 3 = d 2 < 0 , and the third eigenvalue is λ 1 = a n s 1 d 2 ( d 2 s 2 ) c 2 . Boundary equilibrium point E 0 is stable for a < n s 1 d 2 ( d 2 s 2 ) c 2 and is unstable if the inequality is reversed. E 0 undergoes a transcritical bifurcation at a T C = n s 1 d 2 ( d 2 s 2 ) c 2 , and E bifurcates from E 0 for a > a T C . λ 1 is a simple eigenvalue of J 0 when a > a T C .
The eigenvectors associated with the zero eigenvalue of J 0 and J 0 T are given by
V = ( d 2 s 2 ) 2 k s 1 s 2 1 ( d 2 s 2 ) 2 s 1 d 2 2 , W = 1 0 0 .
Following the same notations as in [34], we have to check the transversality conditions. First, we calculate
W T F a ( E 0 ; a T C ) = W T x x 2 0 0 E 0 = 0 ,
where F is defined in (19) and F a is the partial derivative of F with respect to a.
Next, we can calculate
W T D F a ( E 0 ; a T C ) V = ( d 2 s 2 ) 2 k s 1 s 2 0 ,
and
W T D 2 F ( E 0 ; a T C ) ( V , V ) = 2 ( d 2 s 2 ) 2 k s 1 s 2 a k s 1 s 2 + n c 2 0 .
Hence, all the transversality conditions are satisfied and the transcritical bifurcation is nondegenerate. Therefore, we have the following result.
Theorem 4.
System (3) exhibits a transcritical bifurcation for a > a T C , where a T C = n s 1 d 2 ( d 2 s 2 ) c 2 .

3.3. Hopf Bifurcation

In order to discuss how system (3) undergoes Hopf bifurcation near the tumor-presence equilibrium E , we choose k as a bifurcation parameter and define
ϕ ( k ) = b 1 ( k ) b 2 ( k ) b 3 ( k ) ,
where b 1 , b 2 and b 3 are given in (11). Without loss of generality, we suppose that k = k H such that
ϕ ( k H ) = ( b 1 ( k ) b 2 ( k ) b 3 ( k ) ) | k = k H = 0 .
Under the assumption b 1 b 2 b 3 = 0 , b 1 > 0 and b 2 > 0 , we have that characteristic Equation (9) has one negative eigenvalue b 1 and two purely imaginary eigenvalues ± i ω , where ω > 0 is a real number and satisfies ω 2 = b 2 . By the assumption that k = k H and b 2 ( k H ) > 0 , characteristic Equation (9) becomes
( λ 2 + b 2 ) ( λ + b 1 ) = 0 ,
which has three roots
λ 1 = i b 2 ( k H ) , λ 2 = i b 2 ( k H ) , λ 3 = b 1 ( k H ) .
When 0 < | k k H | 1 , the roots can be expressed as the general form
λ 1 ( k ) = μ ( k ) + i ν ( k ) , λ 2 ( k ) = μ ( k ) i ν ( k ) , λ 3 ( k ) = b 1 ( k ) < 0 ,
where μ ( k H ) = 0 , ν ( k H ) = b 2 ( k H ) .
Next, we verify the transversality condition. Substituting λ 1 ( k ) = μ ( k ) + i ν ( k ) into the characteristic Equation (9), calculating the derivative with respect to k at k H and separating the real and imaginary parts, we achieve that
( b 2 ( k H ) 3 ν 2 ( k H ) ) μ ( k H ) 2 b 1 ( k H ) ν ( k H ) ν ( k H ) + b 3 ( k H ) ν 2 ( k H ) b 1 ( k H ) = 0 2 b 1 ( k H ) ν ( k H ) μ ( k H ) + ( b 2 ( k H ) 3 ν 2 ( k H ) ) ν ( k H ) + ν ( k H ) b 2 ( k H ) = 0 .
Since ν ( k H ) = b 2 ( k H ) , by eliminating ν ( k H ) in (26), we obtain that
d μ ( k ) d k | k = k H = ( b 2 ( k H ) 3 ν 2 ( k H ) ) ( b 3 ( k H ) ν 2 ( k H ) b 1 ( k H ) ) + 2 b 1 ( k H ) ν 2 ( k H ) b 2 ( k H ) ( b 2 ( k H ) 3 ν 2 ( k H ) ) 2 + 4 ν 2 ( k H ) b 1 2 ( k H ) = b 1 ( k H ) b 2 ( k H ) + b 2 ( k H ) b 1 ( k H ) b 3 ( k H ) 2 ( b 1 2 ( k H ) + b 2 ( k H ) ) = 1 2 ( b 1 2 ( k H ) + b 2 ( k H ) ) d ϕ ( k ) d k | k = k H .
Therefore, if there exists k H > 0 such that d ϕ ( k ) d k | k = k H 0 , then the transversality condition d μ ( k ) d k | k = k H 0 can hold. Moreover, since λ 3 ( k ) = b 1 ( k ) < 0 , Hopf bifurcation can occur at k = k H by the Hopf bifurcation theorem [35,36]. Therefore, we have the following result.
Theorem 5.
Suppose the tumor-presence equilibrium E exists. If there exists k H > 0 such that b 2 ( k H ) > 0 , ϕ ( k H ) = 0 and d ϕ ( k ) d k | k = k H 0 , system (3) undergoes a Hopf bifurcation near E at k = k H .

4. Numerical Simulations

In this section, we perform numerical simulations to show how parameters k and c affect the dynamics of system (3). Here, the tumor-free equilibrium is named E 0 . Furthermore, three tumor-presence equilibria are named as follows: low tumor equilibrium E 1 , intermediate tumor equilibrium E 2 and high tumor equilibrium E 3 . We choose the following parameter set [11,20,22]:
a = 10 , n = 3 , s 1 = 0.31 , s 2 = 0.5 , d 2 = 1 .
First, we provide the biparametric bifurcation diagram showing the Bogdanov–Takens (BT) bifurcation point where the Hopf bifurcation curve meets with the saddle-node bifurcation curve (see Figure 4). Furthermore, two saddle-node bifurcation curves intersect at the cusp point (CP). On one branch of the Hopf bifurcation, curve we find a generalized Hopf bifurcation point, which is marked as GH.
Next, we provide the bifurcation diagrams of variable x with respect to k by choosing the six different fixed values of c = 0.35 , 0.43 , 0.44 , 0.45 , 0.46 , 0.464 . We find that system (3) can undergo saddle-node bifurcation and Hopf bifurcation (see Figure 5).
Case 1: c = 0.35 (satisfying f ( 0 ) < g ( 0 ) , see Figure 5a).
  • There is one saddle-node bifurcation point named k S N = 0.0589 . For any k > 0 , the tumor-free equilibrium E 0 always exists, which is stable. When 0 < k < k S N , there exist two tumor-presence equilibria: unstable E 2 and stable E 3 . Furthermore, we find the bistability of E 0 and E 3 .
Case 2: c = 0.43 (satisfying f ( 0 ) < g ( 0 ) , see Figure 5b).
  • There is one saddle-node bifurcation point named k S N = 0.2538 and one Hopf bifurcation point named k H = 0.2398 . For any k > 0 , E 0 always exists, which is stable. When 0 < k < k S N , there exist E 2 and E 3 , where E 2 is unstable. E 3 is stable when 0 < k < k H and unstable when k H < k < k S N . Furthermore, we find that the Hopf bifurcating limit cycle disappears through a homoclinic bifurcation k H C L .
Case 3: c = 0.44 (satisfying f ( 0 ) > g ( 0 ) , see Figure 5c).
  • There is one saddle-node bifurcation point named k S N = 0.2835 and a Hopf bifurcation point named k H = 0.261 . For any k > 0 , E 0 always exists, which is unstable. When 0 < k < k S N , there exist E 1 , E 2 and E 3 , where E 1 is stable and E 2 is unstable. E 3 is stable when 0 < k < k H and unstable when k H < k < k S N . When k > k S N , there exists E 1 , which is stable.
Case 4: c = 0.45 (satisfying f ( 0 ) > g ( 0 ) , see Figure 5d).
  • There are two saddle-node bifurcation points named k S N 1 = 0.1774 and k S N 2 = 0.3171 and one Hopf bifurcation point named k H = 0.284 . For any k > 0 , E 0 always exists, which is unstable. When 0 < k < k S N 1 , there exists E 3 , which is stable. When k S N 1 < k < k S N 2 , there exist E 1 , E 2 and E 3 , where E 1 is stable and E 2 is unstable. E 3 is stable when k S N 1 < k < k H and unstable when k H < k < k S N 2 . When k > k S N 2 , there exists E 1 , which is stable. By fixing k = 0.19 ( k S N 1 < 0.19 < k H ) , we can obtain E 1 ( 0.1314 , 0.6363 , 0.5128 ) , E 2 ( 0.189 , 0.644 , 0.5186 ) , E 3 ( 0.6593 , 0.7236 , 0.5716 ) and find the bistability of E 1 and E 3 .
Case 5: c = 0.46 (satisfying f ( 0 ) > g ( 0 ) , see Figure 5e).
  • There are saddle-node bifurcation points named k S N 1 = 0.347 and k S N 2 = 0.3586 and two Hopf bifurcation points named k H 1 = 0.306 and k H 2 = 0.3927 . For any k > 0 , E 0 always exists, which is unstable. When 0 < k < k S N 1 , there exists E 3 . E 3 is stable when 0 < k < k H 1 , and unstable when k H 1 < k < k S N 1 . When k S N 1 < k < k S N 2 , there exist E 1 , E 2 and E 3 , where they are all unstable. When k > k S N 2 , there exists E 1 . E 1 is unstable when k S N 2 < k < k H 2 and stable when k > k H 2 . By fixing k = 0.39 ( k H 1 < 0.39 < k H 2 ) , we can obtain E 1 = ( 0.1598 , 0.6641 , 0.5332 ) and find stable periodic solutions around E 1 (see Figure 6). Furthermore, we find two homoclinic bifurcation thresholds k H C L 1 and k H C L 2 in this case.
Figure 6. (a) The time evolution curves of x , y , z of system (3) around E 1 = ( 0.1598 , 0.6641 , 0.5332 ) , which show periodic oscillations. (b) The 3D phase portrait depicts a stable limit cycle corresponding to (a).
Figure 6. (a) The time evolution curves of x , y , z of system (3) around E 1 = ( 0.1598 , 0.6641 , 0.5332 ) , which show periodic oscillations. (b) The 3D phase portrait depicts a stable limit cycle corresponding to (a).
Symmetry 16 00676 g006
Case 6: c = 0.464 (satisfying f ( 0 ) > g ( 0 ) , see Figure 5f).
  • There are two Hopf bifurcation points named k H 1 and k H 2 . For any k > 0 , E 3 always exists. E 3 is stable when 0 < k < k H 1 or k > k H 2 and unstable when k H 1 < k < k H 2 .
Next, we provide the bifurcation diagrams of variable x with respect to c by choosing the four different typical cases of k = 0.1 , 0.19 , 0.3 , 0.4 . We find that system (3) can appear with saddle-node bifurcation, transcritical bifurcation and Hopf bifurcation (see Figure 7). As k increases, upper saddle-node bifurcation point c S N 1 will disappear, which means that E 2 and E 3 will merge, and lower saddle-node bifurcation point c S N 2 will disappear, which means that E 1 and E 2 will merge. Furthermore, upper Hopf bifurcation point c H 2 and lower Hopf bifurcation point c H 1 will appear one after another, causing the stability switch of E 3 . As k increases, the curve of x of the tumor-presence equilibrium goes up on the right side, and the upper branch goes down. In the case of k = 0.19 from Figure 7b, we find two bistable scenarios: E 0 and E 3 ( c H < c < c T C , see Figure 8a,b), and E 1 and E 3 ( c T C < c < c S N 2 , see Figure 8c,d), where c H = 0.4076 , c T C = 0.4316 , c S N 2 = 0.4506 . For larger values of k (e.g., k = 0.4 ), there is only one tumor-presence equilibrium E as c varies (see Figure 7d). Furthermore, we find two Hopf bifurcation points at c H 1 = 0.464 and c H 2 = 0.5381 .
When k = 0.19 , we identify five distinct regions separated by four critical points as c varies: two saddle-node bifurcation points c S N 1 = 0.4059 and c S N 2 = 0.4506 , one Hopf bifurcation point c H = 0.4076 and one transcritical bifurcation point c T C = 0.4316 (see Figure 7b). E 0 always exists. Furthermore, E 0 is stable when 0 < c < c T C and unstable when c > c T C . The qualitative behavior in each region as c varies is summarized as follows:
(1)
0 < c < c S N 1 : there is no tumor-presence equilibrium.
(2)
c = c S N 1 : one saddle-node bifurcation point at which E 2 and E 3 will appear. In this case, the functions f ( x ) and g ( x ) are tangent.
(3)
c S N 1 < c < c H : there exist E 2 and E 3 . Both E 2 and E 3 are unstable.
(4)
c H < c < c T C : there exist E 2 and E 3 . E 2 is unstable and E 3 is stable.
(5)
c = c T C : one transcritical bifurcation point at which E 0 loses the stability.
(6)
c T C < c < c S N 2 : there exist E 1 , E 2 and E 3 . Both E 1 and E 3 are stable, while E 2 is unstable.
(7)
c = c S N 2 : one saddle-node bifurcation at which E 1 and E 2 will merge. In this case, the functions f ( x ) and g ( x ) are tangent.
(8)
c > c S N 2 : there exists E 3 , which is stable.

5. Conclusions and Discussion

In this paper, we proposed a tumor–immune system interaction model with nonmonotonic immune response function in the presence of a large number of tumors. We studied the existence and stability of tumor-free equilibrium E 0 and tumor-presence equilibrium E of system (3). Our results showed that E 0 undergoes a transcritical bifurcation at a T C = n s 1 d 2 ( d 2 s 2 ) c 2 and E bifurcates from E 0 for a > a T C . We found that system (3) displays a saddle-node bifurcation in which the sudden appearance of an equilibrium point which splits into two equilibria can be found. We also present a biparametric bifurcation diagram to show the appearance of the BT bifurcation point and cusp bifurcation point (see Figure 4). In addition, Hopf bifurcations at E were detected by choosing k (see Figure 5) and c (see Figure 7) as the bifurcation parameter, respectively. The Hopf bifurcation is supercritical and can induce stable periodic solutions (see Figure 6). Interestingly, we also found the bistability of two tumor-presence equilibria named the low-tumor-presence equilibrium and the high-tumor-presence equilibrium (see Figure 8). Two types of bistability phenomena show symmetric and antisymmetric properties in a range of model parameters and initial cell concentrations, which reveals that patients with ACI can possibly reach either a stable tumor-free state or a low-tumor-presence state when the nonmonotonic immune response is properly activated.
The two-dimensional dynamic model proposed by Kuznetsov et al. [20] can also have at most three tumor-presence equilibria with the bistability of small and large tumor-presence equilibria, but the model has no closed orbits, shown by the Dulac–Bendixson criterion [37]. The three-dimensional dynamic model proposed by Dong et al. [11] can also exhibit a Hopf bifurcation inducing stable periodic orbits, but the model has only a unique tumor-presence equilibrium. However, the consideration of nonmonotonic immune response function in our model (3) can capture both the bistability phenomenon and periodic orbits, showing rich dynamics. Comparing the monotonic immune response function to, for example, the saturated type used in the modeling of the immune system against the tumor [21], the nonmonotonic immune response function adopted by this study made the theoretical analysis more challenging while also bringing new, interesting findings. In this study, we only considered the constant input of ACI treatment as one immunotherapy. In further studies, we can introduce other tumor therapies (radiotherapy, chemotherapy, CAR T cell therapy and oncolytic virotherapy) into the model (1) to evaluate the effects of combination treatments for tumor management [38]. In addition, due to proliferation, differentiation and transport of immune cells, the immune system needs some time to be activated, which also plays an important role in tumor–immune system interactions [39].

Author Contributions

Conceptualization, Y.T., M.B. and Y.D.; methodology, Y.L., Y.M., C.Y., Z.P., Y.T., M.B. and Y.D.; software, Y.L., Y.M., M.B. and Y.D.; investigation, Y.L., Y.M., C.Y., Z.P. and Y.T.; writing—original draft preparation, Y.L., M.B. and Y.D.; writing—review and editing, Y.L., Y.M., C.Y., Z.P., Y.T., M.B. and Y.D.; supervision, M.B. and Y.D.; project administration, Y.D.; funding acquisition, Y.T. and Y.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (No. 12371488), and by the Japan Society for the Promotion of Science “Grand-in-Aid 20K03755”.

Data Availability Statement

Data are contained within the article.

Acknowledgments

We would like to thank the four anonymous referees for their valuable comments which helped us to improve the manuscript significantly.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Vesely, M.D.; Schreiber, R.D. Cancer Immunoediting: Antigens, mechanisms and implications to cancer immunotherapy. Ann. N. Y. Acad. Sci. 2013, 1284, 1–5. [Google Scholar] [CrossRef] [PubMed]
  2. Yang, Y.L.; Yang, F.; Huang, Z.Q.; Li, Y.Y.; Shi, H.Y.; Sun, Q.; Xu, F.H. T cells, NK cells, and tumor-associated macrophages in cancer immunotherapy and the current state of the art of drug delivery systems. Front. Immunol. 2023, 14, 1199173. [Google Scholar] [CrossRef] [PubMed]
  3. Beatty, G.L.; Gladney, W.L. Immune escape mechanisms as a guide for cancer immunotherapy. Clin. Cancer Res. 2015, 21, 687–692. [Google Scholar] [CrossRef] [PubMed]
  4. Kim, S.K.; Cho, S.W. The evasion mechanisms of cancer immunity and drug intervention in the tumor microenvironment. Front. Pharmacol. 2022, 13, 868695. [Google Scholar] [CrossRef] [PubMed]
  5. Tan, S.; Li, D.; Zhu, Z. T cells, Cancer immunotherapy: Pros, cons and beyond. Biomed. Pharmacother. 2020, 124, 109821. [Google Scholar] [CrossRef] [PubMed]
  6. Farkona, S.; Diamandis, E.P.; Blasutig, I.M. Cancer immunotherapy: The beginning of the end of cancer? BMC Med. 2016, 14, 73. [Google Scholar] [CrossRef] [PubMed]
  7. Waldman, A.D.; Fritz, J.M.; Lenardo, M.J. A guide to cancer immunotherapy: From T cell basic science to clinical practice. Nat. Rev. Immunol. 2020, 20, 651–668. [Google Scholar] [CrossRef] [PubMed]
  8. Eftimie, R.; Bramson, J.L.; Earn, D.J.D. Interactions between the immune system and cancer: A brief review of non-spatial mathematical models. Bull. Math. Biol. 2011, 73, 2–32. [Google Scholar] [CrossRef]
  9. de Pillis, L.G.; Radunskaya, A.E. A mathematical tumor model with immune resistance and drug therapy: An optimal control approach. Computat. Math. Methods Med. 2001, 3, 79–100. [Google Scholar]
  10. de Pillis, L.G.; Radunskaya, A.E.; Wiseman, C.L. A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. 2005, 65, 7950–7958. [Google Scholar] [CrossRef]
  11. Dong, Y.; Miyazaki, R.; Takeuchi, Y. Mathematical modeling on helper T cells in a tumor immune system. Discret. Contin. Dyn. Syst. Ser. B 2014, 19, 55–72. [Google Scholar] [CrossRef]
  12. Wilson, S.; Levy, D. A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy. Bull. Math. Biol. 2012, 74, 1485–1500. [Google Scholar] [CrossRef]
  13. den Breems, N.Y.; Eftimie, R. The re-polarisation of M2 and M1 macrophages and its role on cancer outcomes. J. Theor. Biol. 2016, 390, 23–39. [Google Scholar] [CrossRef] [PubMed]
  14. Shu, Y.; Huang, J.; Dong, Y.; Takeuchi, Y. Mathematical modeling and bifurcation analysis of pro- and anti-tumor macrophages. Appl. Math. Model. 2020, 88, 758–773. [Google Scholar] [CrossRef]
  15. d’Onofrio, A. Metamodeling tumor-immune systeminteraction, tumor evasion and immunotherapy. Math. Comput. Model. 2008, 47, 614–637. [Google Scholar] [CrossRef]
  16. Mahlbacher, G.E.; Reihmer, K.C.; Frieboes, H.B. Mathematical modeling of tumor-immune cell interactions. J. Theor. Biol. 2019, 469, 47–60. [Google Scholar] [CrossRef] [PubMed]
  17. Tao, Y.; Campbell, S.A.; Poulin, F.J. Dynamics of a diffusive nutrient-phytoplankton-zooplankton model with spatio-temporal delay. SIAM J. Appl. Math. 2021, 81, 2405–2432. [Google Scholar] [CrossRef]
  18. Xu, Q.; Huang, J.; Dong, Y.; Takeuchi, Y. A delayed HIV infection model with the homeostatic proliferation of CD4 +T cells. Acta Math. Appl. Sin. Engl. Ser. 2022, 38, 441–462. [Google Scholar] [CrossRef]
  19. Tao, Y.; Sun, Y.; Zhu, H.; Lyu, J.; Ren, J. Nilpotent singularities and periodic perturbation of a GIβ model: A pathway to glucose disorder. J. Nonlinear Sci. 2023, 33, 33–49. [Google Scholar] [CrossRef]
  20. Kuznetsov, V.A.; Makalkin, I.A.; Taylor, M.A.; Perelson, A.S. Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bull. Math. Biol. 1994, 56, 295–321. [Google Scholar] [CrossRef]
  21. Kirschner, D.; Panetta, J. Modeling immunotherapy of the tumor-immune interaction. J. Math. Biol. 1998, 37, 235–252. [Google Scholar] [CrossRef] [PubMed]
  22. Dritschel, H.; Waters, S.L.; Roller, A.; Byrne, H.M. A mathematical model of cytotoxic and helper T cell interactions in a tumour microenvironment. Lett. Biomath. 2018, 5, S36–S68. [Google Scholar] [CrossRef]
  23. Talkington, A.; Dantoin, C.; Durrett, R. Ordinary differential equation models for adoptive immunotherapy. Bull. Math. Biol. 2018, 80, 1059–1083. [Google Scholar] [CrossRef] [PubMed]
  24. Song, G.; Liang, G.; Tian, T.; Zhang, X. Mathematical modeling and analysis of tumor chemotherapy. Symmetry 2022, 14, 704. [Google Scholar] [CrossRef]
  25. Bekker, R.A.; Kim, S.; Pilon-Thomas, S.; Enderling, H. Mathematical modeling of radiotherapy and its impact on tumor interactions with the immune system. Neoplasia 2022, 28, 100796. [Google Scholar] [CrossRef] [PubMed]
  26. Butner, J.D.; Dogra, P.; Chung, C. Mathematical modeling of cancer immunotherapy for personalized clinical translation. Nat. Comput. Sci. 2022, 2, 785–796. [Google Scholar] [CrossRef] [PubMed]
  27. Tang, S.; Li, S.; Wang, X.; Xiao, Y.; Checke, R.A. Hormetic and synergistic effects of cancer treatments revealed by modelling combinations of radio-or chemotherapy with immunotherapy. BMC Cancer 2023, 23, 1040. [Google Scholar] [CrossRef] [PubMed]
  28. Zhang, Y.; Xie, L.; Dong, Y.; Ruan, S.; Takeuchi, Y. Bifurcation analysis in a tumor-immune system interaction model with dendritic cell therapy and immune response delay. SIAM J. Appl. Math. 2023, 83, 1892–1914. [Google Scholar] [CrossRef]
  29. Elmusrati, A.; Wang, J.; Wang, C.Y. Tumor microenvironment and immune evasion in head and neck squamous cell carcinoma. Int. J. Oral Sci. 2021, 13, 24. [Google Scholar] [CrossRef]
  30. Vinay, D.S.; Ryan, E.P.; Pawelec, G.; Talib, W.H.; Stagg, J.; Elkord, E.; Kwon, B.S. Immune evasion in cancer: Mechanistic basis and therapeutic strategies. Semin. Cancer Biol. 2015, 35, S185–S198. [Google Scholar] [CrossRef]
  31. Gonzalez, H.; Hagerling, C.; Werb, Z. Roles of the immune system in cancer: From tumor initiation to metastatic progression. Genes. Dev. 2018, 32, 1267–1284. [Google Scholar] [CrossRef] [PubMed]
  32. Serre, R.; Benzekry, S.; Padovani, L.; Meille, C.; André, N.; Ciccolini, J.; Barbolosi, D. Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy. Cancer Res. 2016, 76, 4931–4940. [Google Scholar] [CrossRef] [PubMed]
  33. Yu, M.; Dong, Y.; Takeuchi, Y. Dual role of delay effects in a tumour-immune system. J. Biol. Dyn. 2017, 11, 334–347. [Google Scholar] [CrossRef] [PubMed]
  34. Perko, L. Differential Equations and Dynamical Systems, 3rd ed.; Springer: New York, NY, USA, 2013. [Google Scholar]
  35. Ruan, S.; Wolkowicz, G.S.K. Bifurcation analysis of a chemostat model with a distributed delay. J. Math. Anal. Appl. 1996, 204, 786–812. [Google Scholar] [CrossRef]
  36. Marsden, J.E.; McKracken, M. The Hopf Bifurcation and Its Applications, 1st ed.; Springer: New York, NY, USA, 1976. [Google Scholar]
  37. Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos, 2nd ed.; Springer: New York, NY, USA, 2003. [Google Scholar]
  38. Malinzi, J.; Basita, K.B.; Padidar, S.; Adeola, H.A. Prospect for application of mathematical models in combination cancer treatments. Int. J. Med. Inform. 2021, 23, 100534. [Google Scholar] [CrossRef]
  39. Ruan, S. Nonlinear dynamics in tumor-immune system interaction models with delays. Discrete Contin. Dyn. Syst. Ser. B 2021, 26, 541–602. [Google Scholar] [CrossRef]
Figure 4. The biparametric bifurcation diagram of system (3) in ( k , c ) plane. Magenta and blue curves are saddle-node bifurcation and Hopf bifurcation, respectively. BT, CP and GH denote Bogdanov–Takens bifurcation point, cusp bifurcation point and generalized Hopf bifurcation point, respectively.
Figure 4. The biparametric bifurcation diagram of system (3) in ( k , c ) plane. Magenta and blue curves are saddle-node bifurcation and Hopf bifurcation, respectively. BT, CP and GH denote Bogdanov–Takens bifurcation point, cusp bifurcation point and generalized Hopf bifurcation point, respectively.
Symmetry 16 00676 g004
Figure 5. The bifurcation diagrams of x with respect to the parameter k of system (3) with different values of c = 0.35 , 0.43 , 0.44 , 0.45 , 0.46 , 0.464 . Solid and dashed curves represent tumor-presence equilibrium and tumor-free equilibrium, respectively. Blue and red (magenta) curves correspond to stable equilibrium and unstable equilibrium, respectively. Saddle-node bifurcation points are labeled as S N , S N 1 and S N 2 . Hopf bifurcation points are labeled as H 1 and H 2 . Homoclinic bifurcation thresholds are labeled as H C L , H C L 1 and H C L 2 . Saddle-node bifurcation of periodic orbits are denoted as S N L C 1 and S N L C 2 .
Figure 5. The bifurcation diagrams of x with respect to the parameter k of system (3) with different values of c = 0.35 , 0.43 , 0.44 , 0.45 , 0.46 , 0.464 . Solid and dashed curves represent tumor-presence equilibrium and tumor-free equilibrium, respectively. Blue and red (magenta) curves correspond to stable equilibrium and unstable equilibrium, respectively. Saddle-node bifurcation points are labeled as S N , S N 1 and S N 2 . Hopf bifurcation points are labeled as H 1 and H 2 . Homoclinic bifurcation thresholds are labeled as H C L , H C L 1 and H C L 2 . Saddle-node bifurcation of periodic orbits are denoted as S N L C 1 and S N L C 2 .
Symmetry 16 00676 g005
Figure 7. The bifurcation diagrams of x with respect to the parameter c of system (3) with different values of k = 0.1 , 0.19 , 0.3 , 0.4 . Solid and dashed curves represent tumor-presence equilibrium and tumor-free equilibrium, respectively. Blue and red(magenta) curves correspond to stable equilibrium and unstable equilibrium, respectively. Saddle-node bifurcation points are labeled as S N 1 and S N 2 . Hopf bifurcation points are labeled as H, H 1 and H 2 . Transcritical bifurcation is labeled as T C . Saddle-node bifurcations of periodic orbits are denoted as S N L C 1 and S N L C 2 .
Figure 7. The bifurcation diagrams of x with respect to the parameter c of system (3) with different values of k = 0.1 , 0.19 , 0.3 , 0.4 . Solid and dashed curves represent tumor-presence equilibrium and tumor-free equilibrium, respectively. Blue and red(magenta) curves correspond to stable equilibrium and unstable equilibrium, respectively. Saddle-node bifurcation points are labeled as S N 1 and S N 2 . Hopf bifurcation points are labeled as H, H 1 and H 2 . Transcritical bifurcation is labeled as T C . Saddle-node bifurcations of periodic orbits are denoted as S N L C 1 and S N L C 2 .
Symmetry 16 00676 g007
Figure 8. The time evolution curves of x , y , z of system (3). (a,b) show the bistability of E 0 ( 0 , 0.62 , 0.5 ) and E 3 = ( 0.63 , 0.7176 , 0.568 ) . (c,d) show the bistability of E 1 ( 0.1314 , 0.6363 , 0.5128 ) and E 3 ( 0.6593 , 0.7236 , 0.5716 ) .
Figure 8. The time evolution curves of x , y , z of system (3). (a,b) show the bistability of E 0 ( 0 , 0.62 , 0.5 ) and E 3 = ( 0.63 , 0.7176 , 0.568 ) . (c,d) show the bistability of E 1 ( 0.1314 , 0.6363 , 0.5128 ) and E 3 ( 0.6593 , 0.7236 , 0.5716 ) .
Symmetry 16 00676 g008
Table 1. Model parameters.
Table 1. Model parameters.
ParameterDescriptionReference
aProliferation rate of TCs [20]
b 1 Carrying capacity of TCs [20]
nMaximal inhibition rate of TCs by ECs [20]
cHalf-velocity constant [22]
s 1 Constant inflow of ACI [11,22]
d 1 Death rate of ECs [11,22]
pActivation rate of ECs by HTCs [11]
s 2 Birth rate of HTCs produced in the bone marrow [11,22]
d 2 Death rate of HTCs [11,22]
kHTCs’ stimulation rate by identified tumor antigens [11]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Liu, Y.; Ma, Y.; Yang, C.; Peng, Z.; Takeuchi, Y.; Banerjee, M.; Dong, Y. Modeling the Nonmonotonic Immune Response in a Tumor–Immune System Interaction. Symmetry 2024, 16, 676. https://doi.org/10.3390/sym16060676

AMA Style

Liu Y, Ma Y, Yang C, Peng Z, Takeuchi Y, Banerjee M, Dong Y. Modeling the Nonmonotonic Immune Response in a Tumor–Immune System Interaction. Symmetry. 2024; 16(6):676. https://doi.org/10.3390/sym16060676

Chicago/Turabian Style

Liu, Yu, Yuhang Ma, Cuihong Yang, Zhihang Peng, Yasuhiro Takeuchi, Malay Banerjee, and Yueping Dong. 2024. "Modeling the Nonmonotonic Immune Response in a Tumor–Immune System Interaction" Symmetry 16, no. 6: 676. https://doi.org/10.3390/sym16060676

APA Style

Liu, Y., Ma, Y., Yang, C., Peng, Z., Takeuchi, Y., Banerjee, M., & Dong, Y. (2024). Modeling the Nonmonotonic Immune Response in a Tumor–Immune System Interaction. Symmetry, 16(6), 676. https://doi.org/10.3390/sym16060676

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop