Global stability of an age-structured virus dynamics model with Beddington-DeAngelis infection function

1.   Introduction

The emergence of drug resistance in cancer cell populations during anti-cancer therapies is a major issue that is impeding the treatment of human cancers ^[1]. Most cancers are shown to develop resistance to targeted therapies or chemotherapies, even highly drug-sensitive tumors may show initial responses but develop resistance gradually in the process of treatment, presenting major challenges to personalized medicine ^[2,3]. It is confirmed that the emergence of drug resistance in cancer cell population stems from intratumor heterogeneity, that is, the coexistence of cellular populations bearing different genetic or epigenetic alterations within the same individual tumors ^[4,5]. The interplay between intratumor heterogeneity and the tumor microenvironment drives the emergence of resistance, which can be described by the tumor evolution process. Under environmental pressure, subpopulations of cancer cells expressing phenotypes adapted to the local environment have selective advantages, so that they emerge and expand at the expense of others ^[6]. The evolutionary adaption is driven by the genetic modifications (e.g., mutations) and perturbations of the tumor microenvironment induced by anti-cancer therapy agents. When cancer undergoes selective pressures imposed by anti-cancer therapy, the more stress-resistant phenotypes emerge and expand, thereby cancer adapts to the pharmacologic pressures and affects clinical outcomes. Therefore, it is important to consider the evolution of tumors in implementation of therapeutic strategies ^[7,8].

Certain combination therapies are shown to be more effective than single drugs, and some preventive combination therapies could block the emergence of resistance ^[9,10,11]. The scheduling of combination therapies is important for controlling different cancer cell subpopulations, and the therapies implemented at each clinical time point must be determined in advance carefully. Adaptive strategies targeting different branches of tumor cells with alternating administration of corresponding drugs could foster competition between drug-sensitive and drug-resistant subpopulations ^[12,13,14]. The goal of adaptive therapy design is to keep a controllable and stable tumor load in cancer treatment by allowing for sensitive cells to survive at a certain level in order to suppress proliferation of the less adapted resistant populations ^{[13,15,16,17,18,19,20]}.

Mathematical modeling with evolutionary game theory and evolutionary adaptive dynamics has been widely applied for the study of cancer evolution and design of adaptive therapy schedules ^{[14,15,16,17,18,21]}. Evolutionary game theory is a general modeling framework to study evolutionary dynamics of frequency-dependent selection among populations ^{[22,23,24,25,26,27]}. Individuals with fixed strategies in the population, as players in a game, interact randomly with other individuals. The payoff is interpreted as the fitness which relies on the relative proportions (frequencies) of the different phenotypes in the population. Strategies with higher fitness do well and reproduce faster whereas those with lower fitness do poorly and are outcompeted.

Applying frequency-dependent competition models and evolutionary game theory, West J. et al. ^[13] explored sequential and concomitant therapy strategies and presented novel multidrug adaptive therapy regimens for metastatic castrate resistant prostate cancer patients, that is, patient-specific therapeutic schedules which drive tumor evolution into cycles of a periodic and controllable loop. Zhang J. et al. ^[14] integrated evolutionary dynamics into the treatment of metastatic castrate-resistant prostate cancer by an evolutionary game theory model with Lotka-Volterra equations, and explored strategies for the control of tumor progression which showed significant improvement in outcomes in pilot clinical trial. Singh G. et al. ^[20] investigated the design of adaptive therapy dosage by a Lotka-Volterra based population dynamics model of the drug-sensitive, drug-resistant cells, and transient drug-hybrid state along with phenotypic switching during adaptive therapy. Certain key parameter regimes were identified for effective adaptive therapy, and the importance of intermediate drug-hybrid state in cancer was phenomenologically explained. Gluzman M. et al. ^[24] proposed a method for systematically optimizing the roles of adaptive policies based on an evolutionary game theory model of cancer dynamics and optimized the total drug usage and time to recovery by solving a Hamilton-Jacobi-Bellman equation.

Basanta D. et al. ^[9] proposed a general evolutionary game theory framework for tumor cell evolution under the constraints of combination therapies. Applying a more specific version of the general model, the authors explored evolutionary dynamics that how chemotherapy can improve the efficacy of the p53 vaccine for small cell lung cancer. In their model, three different tumor subpopulations are considered, that is, regular tumor cells that are susceptible to both p53 vaccine and chemotherapy (sensitive cells), cells that are resistant to chemotherapy and cells resistant to p53 vaccine. By numerical simulation, the authors showed that the schedule that p53 vaccine is administered first and then change to chemotherapy results in entire chemotherapy-resistant cells, while the schedule that chemotherapy is applied first and then change to the p53 vaccine also leads to entire chemotherapy-resistant cells, driving the susceptible and p53 vaccine-resistant populations toward extinction. However, for the schedule that chemotherapy is applied first and then change to the p53 vaccine, we found that, as the simulation time step is extended to 200, p53 vaccine-resistant cells would emerge again and establish while chemotherapy-resistant population and sensitive population are eradicated. Explanation of this phenomenon needs investigation of the underlying evolutionary dynamics.

In this paper, we investigated the evolution outcomes of tumor composed of sensitive cells, chemotherapy-resistant cells and p53 vaccine-resistant cells different schedules of chemotherapy and p53 vaccine therapy by applying a replicator dynamical model with fitness defined by a payoff matrix generalized from ^[9]. We further explored the design of adaptive therapy schedules with combination of chemotherapy and p53 vaccine. The rest of this paper is organized as follows. In Section 2, we first presented a generalized replicator dynamical model, and then explored analytically and numerically the existence and local asymptotic stability of equilibria of the model. Furthermore, designs of adaptive therapy schedules based evolutionary velocity and evolutionary absorbing regions were illustrated.

Experimental studies on targeted therapy with BRAF, ALK, or EGFR kinase inhibitors showed that these targeted therapies induces the release of sceretome by drug-sensitive cancer cells, which lead to establishment of a tumor microenvironment that supports the expansion of drug-resistant clones ^[28]. In Section 3, applying the replicator dynamics similar to the model in Section 2, we explored the influence of this supporting effect to the evolution of cancer cell populations under combination of two targeted therapies, via local stability analysis and numerical simulation. The paper is ended by Section 4, where some discussion and conclusions were presented.

2.   Tumor evolution under chemotherapy and p53 vaccine

2.1.   Replicator dynamical model

We considered the evolution dynamics of three subpopulations of tumor cells, that is, regular tumor cells susceptible to both the p53 vaccine and chemotherapy (sensitive cells, $ S $), chemotherapy-resistant cells ($ C $) and p53 vaccine-resistant cells ($ I $). The dynamics and interactions of subpopulations are given by the following replicator equations:

where $ P_C(t) $, $ P_I(t) $ and $ P_S(t) $ represent the proportion of chemotherapy-resistant cells, p53 vaccine-resistant cells and sensitive cells at time $ t $, respectively, with $ P_C+ P_I+P_S = 1 $; $ W_C $, $ W_I $ and $ W_S $ indicate the fitness of chemotherapy-resistant cells, p53 vaccine-resistant cells and sensitive cells, respectively; $ \overline W $ denotes the average fitness of the entire population. To define the fitness, we assumed the competitive interactions between these three cell types are given by the following payoff matrix:

Here, we assumed as in ^[9] that chemotherapy-resistant cells and vaccine-resistant cells pay a cost of resistance, $ c_C $ and $ c_I $, respectively; the cost of vaccine to sensitive cells and chemotherapy-resistant cells are $ d_I $ and $ \alpha d_I $, respectively, when it is administered; the cost of chemotherapy to sensitive cells and vaccine-resistant cells are $ d_C $ and $ \beta d_C $, respectively, when it is applied. Here, $ \alpha $ and $ \beta $ represent the extra cost of chemotherapy-resistant cells and vaccine-resistant cells being subjected to the different drugs ^[9], that is, p53 vaccine and chemotherapy, respectively ($ \alpha, \: \beta > 1 $). We also assumed that sensitive cells and chemotherapy-resistance cells get support from vaccine-resistant cells as the vaccine is applied and as they interact with vaccine-resistant cells ^[9], respectively, which is represented by a factor $ \gamma $ ($ 0 < \gamma < 1 $). The fitness of the three subpopulations are defined by

where $ \vec P = (P_C, P_I, P_S)^T $. The average fitness is given by

The fitness of population can be understood as reproduction rate of the population.

2.2.   Stability analysis of equilibria

According the properties of general replicator dynamics ^[22,23], the system (2.1) is defined on the simplex $ S_3 = \{(P_C, P_I, P_S)|\: P_C+P_I+P_S = 1, \quad 0\le P_C, P_I, P_S\le1\} $, and the edges and the interior of the simplex are invariant, respectively. The system (2.1) has seven types of equilibria $ E_i(P_C, P_I, P_S) $, $ i = 1, ..., 7 $, given as follows.

(I) The equilibria $ E_1(0, 0, 1) $, $ E_2(0, 1, 0) $ and $ E_3(1, 0, 0) $ always exist.

(II) When $ d_C+\gamma d_I < c_I+\beta d_C < d_C+ d_I $, there exists an equilibrium $ E_4(0, \bar P_I, \bar P_S) $, where $ \bar P_I = \frac{(d_C+d_I)-(c_I+\beta d_C)}{d_I(1-\gamma)} $, $ \bar P_S = 1-\bar P_I = \frac{(c_I+\beta d_C)-(d_C+\gamma d_I)}{d_I(1-\gamma)}) $; vaccine-resistant cells and sensitive cells coexist.

(III) When $ c_C+\alpha \gamma d_I < c_I+\beta d_C < c_C+\alpha d_I $, there exists an equilibrium $ E_5(\tilde P_C, \tilde P_I, 0) $, where $ \tilde P_C = \frac{(c_I+\beta d_C)-(c_C+\alpha \gamma d_I)}{\alpha d_I(1-\gamma)} $, $ \tilde P_I = 1-\tilde P_C = \frac{(c_C+\alpha d_I)-(c_I+\beta d_C)}{\alpha d_I(1-\gamma)} $; chemotherapy-resistant cells and vaccine-resistant cells coexist.

(IV) When $ d_C+d_I = c_C+\alpha d_I $, there exists a line of equilibria $ E_6(\hat P_C, 0, \hat P_S) $, where $ \hat P_C + \hat P_S = 1 $; chemotherapy-resistant cells and sensitive cells coexist.

(V) When $ d_C+\gamma d_I < c_I+\beta d_C < d_C+d_I $ and $ (c_C+\alpha d_I) - (c_I+\beta d_C) = \alpha[(d_C+ d_I)-(c_I+\beta d_C)] $, there exists a line of equilibria $ E_7\left(P_C^*, P_I^*, P_S^*\right) $, where $ P_I^* = \frac{(d_C+d_I)-(c_I+\beta d_C)}{d_I(1-\gamma)} $, $ P_C^*+P_S^* = 1-P_I^* $; chemotherapy-resistant cells, vaccine-resistant cells and sensitive cells coexist.

Next we considered the local stability of these equilibria by linearizing the system (2.1) at the equilibria. The stability of equilibria are shown in the following theorems.

Theorem 2.1. The equilibria $ E_1(0, 0, 1) $, $ E_2(0, 1, 0) $ and $ E_3(1, 0, 0) $ of the system (2.1) always exist and their stability conditions are given as follows:

(1) If $ d_C+d_I < {\rm min}\{c_C+\alpha d_I, c_I+\beta d_C, 1\} $, then $ E_1(0, 0, 1) $ is LAS;

(2) If $ c_I+\beta d_C < {\rm min}\{ c_C+\alpha \gamma d_I, d_C+\gamma d_I, 1\} $, then $ E_2(0, 1, 0) $ is LAS;

(3) If $ c_C+\alpha d_I < {\rm min}\{ d_C+d_I, c_I+\beta d_C, 1\} $, then $ E_3(1, 0, 0) $ is LAS.

Here, LAS means locally asymptotically stable.

Proof. (1) The Jacobian matrix of the system (2.1) at $ E_1(0, 0, 1) $ is given by

The corresponding eigenvalues are $ \lambda_1 = (d_{C}+d_{I})-(c_{C}+\alpha d_{I}) $, $ \lambda_2 = (d_{C}+d_{I})-(c_{I}+\beta d_{C}) $, $ \lambda_3 = d_{C}+d_{I}-1 $. Thus, if $ d_C+d_I < {\rm min}\{c_C+\alpha d_I, c_I+\beta d_C, 1\} $, then $ \lambda_i < 0 $, $ i = 1, 2, 3 $, which implies that the equilibrium $ E_1(0, 0, 1) $ is locally asymptotically stable^[29,30].

(2) The Jacobian matrix of the system (2.1) at the equilibrium $ E_2(0, 1, 0) $ is

The corresponding eigenvalues are $ \lambda_1 = (c_{I}+\beta d_{C})-(c_{C}+\alpha \gamma d_{I}) $, $ \lambda_2 = (c_{I}+\beta d_{C})-1 $, $ \lambda_3 = (c_{I}+\beta d_{C})-(d_{C}+\gamma d_{I}) $, which are negative if $ c_I+\beta d_C < {\rm min}\{ c_C+\alpha \gamma d_I, d_C+\gamma d_I, 1\} $. Hence, the equilibrium $ E_2(0, 1, 0) $ is locally asymptotically stable if $ c_I+\beta d_C < {\rm min}\{ c_C+\alpha \gamma d_I, d_C+\gamma d_I, 1\} $.

(3) The Jacobian matrix of the system (2.1) at $ E_3(1, 0, 0) $ reads

The corresponding eigenvalues are $ \lambda_1 = (c_{C}+\alpha d_{I})-1 $, $ \lambda_2 = (c_{C}+\alpha d_{I})-(c_{I}+\beta d_{C}) $, $ \lambda_3 = (c_{C}+\alpha d_{I})-(d_{C}+d_{I}) $. We see that $ c_C+\alpha d_I < {\rm min}\{ d_C+d_I, c_I+\beta d_C, 1\} $ implies $ \lambda_i < 0 $, $ i = 1, 2, 3 $, so that the equilibrium $ E_3(1, 0, 0) $ is locally asymptotically stable.

Theorem 2.2. For the system (2.1),

(1) If $ d_C+\gamma d_I < c_I+\beta d_C < d_C+d_I $, the equilibrium $ E_4(0, \bar P_I, \bar P_S) $ exists, and furthermore if $ c_I+\beta d_C < 1 $ and $ \alpha [(d_C+d_I)-(c_I+\beta d_C)] < (c_C + \alpha d_I)-(c_I+\beta d_C) $, then $ E_4(0, \bar P_I, \bar P_S) $ is LAS, where $ \bar P_I = \frac{(d_C+d_I)-(c_I+\beta d_C)}{d_I(1-\gamma)} $, $ \bar P_S = 1-\bar P_I = \frac{(c_I+\beta d_C)-(d_C+\gamma d_I)}{d_I(1-\gamma)} $;

(2) If $ c_C+\alpha \gamma d_I < c_I+\beta d_C < c_C+\alpha d_I $, the equilibrium $ E_5(\tilde P_C, \tilde P_I, 0) $ exits, if in addition, $ c_I+\beta d_C < 1 $ and $ \alpha [(d_C+d_I)- (c_I+\beta d_C)] > (c_C+\alpha d_I) - (c_I+\beta d_C) $, then $ E_5(\tilde P_C, \tilde P_I, 0) $ is LAS, where $ \tilde P_C = \frac{(c_I+\beta d_C)-(c_C+\alpha \gamma d_I)}{\alpha d_I(1-\gamma)} $, $ \tilde P_I = 1-\tilde P_C = \frac{(c_C+\alpha d_I)-(c_I+\beta d_C)}{\alpha d_I(1-\gamma)} $.

Proof. (1) The fitness of the three subpopulations at the equilibrium $ E_4(0, \bar P_I, \bar P_S) $ are given by

The Jacobian matrix of the system (2.1) at $ E_4(0, \bar P_I, \bar P_S) $ reads

One of the corresponding eigenvalues is $ \lambda_1 = \overline{W}_{C}-\overline{W}_{I} = (c_I+\beta d_C)-(c_{C}+\alpha d_{I}) + \alpha [(d_C+d_I))-(c_I+\beta d_C)] $, which is negative if $ \alpha [(d_C+d_I)-(c_I+\beta d_C)] < (c_C + \alpha d_I)-(c_I+\beta d_C) $. The other two eigenvalues ($ \lambda_2, \lambda_3 $) have negative real parts if the following conditions are satisfied:

which are equivalent to

respectively, noticing that $ \bar P_I+\bar P_S = 1 $, $ \overline{W}_I = \overline{W}_S $, $ \bar P_I, \bar P_S, d_I > 0 $ and $ 0 < \gamma < 1 $. Hence, $ \overline W_{S} = 1-(c_I+\beta d_C) > 0 $ ensures the two eigenvalues $ \lambda_2 $ and $ \lambda_3 $ have negative real parts. Therefore, if $ c_I+\beta d_C < 1 $, $ \alpha [(d_C+d_I)-(c_I+\beta d_C)] < (c_C + \alpha d_I) - (c_I+\beta d_C) $, and the equilibrium $ E_4(0, \bar P_I, \bar P_S) $ exists, then it is locally asymptotically stable.

(2) The Jacobian matrix of the system (2.1) at $ E_5(\tilde P_C, \tilde P_I, 0) $ is given by

where

We see that one corresponding eigenvalue is $ \lambda_1 = \tilde W_{S}-\tilde W_{I} = (c_I+\beta d_C)-(d_{C}+d_{I}) + \frac{1}{\alpha}[(c_{C}+\alpha d_{I})-(c_{I}+\beta d_{C})] $, so that $ (c_C + \alpha d_I)-(c_I+\beta d_C) < \alpha [(d_C+d_I)-(c_I+\beta d_C)] $ implies $ \lambda_1 < 0 $. The other two eigenvalues $ \lambda_2, \lambda_3 $ have negative real parts, if

which are equivalent to

respectively, noticing that $ \tilde P_C+\tilde P_I = 1 $, $ \tilde {W}_C = \tilde {W}_I $, $ \tilde P_I, \tilde P_C, d_I > 0 $ and $ 0 < \gamma < 1 $. Therefore, $ \tilde W_{C} = 1-(c_I+\beta d_C) > 0 $ guarantees $ {\rm Re}(\lambda_i) < 0 $, $ i = 2, 3 $. In summary, if the equilibrium $ E_5(\tilde P_C, \tilde P_I, 0) $ exits, and in addition, $ c_I+\beta d_C < 1 $ and $ (c_C + \alpha d_I)-(c_I+\beta d_C) < \alpha [(d_C+d_I)-(c_I+\beta d_C)] $, then $ E_5(\tilde P_C, \tilde P_I, 0) $ is locally asymptotically stable.

Theorem 2.3. If and only if $ d_C+d_I = c_C+\alpha d_I $, there exists a line of equilibria, $ E_6(\hat P_C, 0, \hat P_S) $, where $ \hat P_C+\hat P_S = 1 $, and they are locally stable if $ c_{C}+\alpha d_{I} < \min\{c_I+\beta d_C, \: 1\} $.

Proof. The fitness of the subpopulations at the equilibrium $ E_6(\hat P_C, 0, \hat P_S) $ are

Note that $ \hat W_C = \hat W_S $, if $ E_6(\hat P_C, 0, \hat P_S) $ exists. The Jacobian matrix of the system (2.1) at $ E_6(\hat P_C, 0, \hat P_S) $ is given by

The corresponding eigenvalues are $ \lambda_1 = 0 $, $ \lambda_2 = \hat W_I-\hat W_C = (c_{C}+\alpha d_{I})-(c_{I}+\beta d_{C}) $, $ \lambda_3 = -\hat W_C = (c_{C}+\alpha d_{I})-1 $. Hence, if $ c_{C}+\alpha d_{I} < \min\{c_{I}+\beta d_{C}, 1\} $, then $ \lambda_2, \lambda_3 < 0 $, so that $ E_6(\hat P_C, 0, \hat P_S) $ is locally stable when it exists (shown numerically in the next subsection).

Theorem 2.4. If $ d_C+ \gamma d_I < c_I+\beta d_C < d_C+d_I $ and $ (c_C+\alpha d_I)-(c_I+\beta d_C) = \alpha[(d_C+ d_I)-(c_I+\beta d_C)] $, there exists a line of equilibria $ E_7\left(P_C^*, P_I^*, P_S^*\right) $, where $ P_I^* = \frac{(d_C+d_I)-(c_I+\beta d_C)}{d_I(1-\gamma)} $, $ P_C^*+P_S^* = 1-P_I^* $; they are locally stable if $ c_I+\beta d_C < 1 $.

Proof. The Jacobian matrix of the system (2.1) at $ E_7\left(P_C^*, P_I^*, P_S^*\right) $ reads

where $ W^* = W_C^* = W_I^* = W_S^* = =1-c_{I}-\beta d_{C} $. Using the elementary transformation $ T^{-1}J_7T $, with

the matrix $ J_7 $ is transformed into

where $ a = \alpha (1-\gamma)d_{I}P_{C}^*P_{I}^*+(1-\gamma)d_{I}P_{I}^*P_{S}^*-P_{C}^*W^*-P_{S}^*W^* $. The two similar matrices $ J_7 $ and $ J^* $ have the same eigenvalues, one of which is $ \lambda_1 = 0 $. The real parts of the other two eigenvalues $ \lambda_2, \lambda_3 $ are negative as the following conditions hold:

which are equivalent to

respectively, where $ P_C^*, P_I^*, P_S^*, \alpha, d_I > 0 $ and $ 0 < \gamma < 1 $. Therefore, $ W^* = 1 -c_{I}-\beta d_{C} > 0 $ permits the other two eigenvalues have negative real parts, and the equilibrium $ E_7\left(P_C^*, P_I^*, P_S^*\right) $ is locally stable (shown numerically in the next subsection), where $ P_I^* = \frac{(d_C+d_I)-(c_I+\beta d_C)}{d_I(1-\gamma)} $, $ P_C^*+P_S^* = 1-P_I^* $.

The conditions for existence and local stability of the equilibria $ E_i $, $ i = 1, ..., 7 $, were summarized in Table 1, with introduction of the notations $ f_C: = c_C+\alpha d_I $, $ f_I : = c_I+\beta d_C $, $ f_S: = d_C+d_I $, $ f_C^*: = c_C+\gamma \alpha d_I $ and $ f_S^*: = d_C+\gamma d_I $.

Notice that the existence and local stability of the evolutionary stable points $ E_i $, $ i = 1, ..., 7 $, are mutually exclusive. Here, we explain some of them, except for the obvious ones. For example, $ E_1(0, 0, 1) $ and $ E_5(\tilde P_C, \tilde P_I, 0) $ cannot be bistable: The existence condition of $ E_5 $, $ f_I < f_C $, and its LAS condition, $ f_C-f_I < \alpha (f_S-f_I) $, together imply that $ f_I < f_S $, under which, $ E_1 $ is unstable. Similarly, $ E_3(0, 0, 1) $ and $ E_4(0, \bar{P_I}, \bar{P_S}) $ (or $ E_7(P_C^*, P_I^*, P_S^*) $) cannot be bistable. The stability condition of $ E_2(0, 1, 0) $, $ f_I < f_S^* $, and the stability condition of $ E_6(\hat{P_C}, 0, \hat{P_S}) $, $ f_S < f_I $, cannot hold simultaneously since $ f_S^* < f_S $ for $ 0 < \gamma < 1 $; hence $ E_2 $ and $ E_6 $ cannot be bistable.

The replicator system (2.1) has the property that if one cell subpopulation does not exist initially, then it cannot emerge forever. In this case, the dynamics of the system is determined by the interaction between other two subpopulations. For example, there are three different possibilities for convergence of the solutions with initial values $ P_S = 0 $, $ 0 < P_C, P_I < 1 $, that is, evolutionary stable points $ E_2(0, 1, 0) $, $ E_3(1, 0, 0) $ and $ E_5(\tilde P_C, \tilde P_I, 0) $; the conditions for convergence to the three evolutionary stable points are given in Table 2. The convergence possibilities of solutions of (2.1) in the other two cases and the corresponding conditions are also shown in Table 2.

2.3.   Numerical simulation results

2.3.1.   Evolution of cancer cell populations to evolutionary stable points

The evolution trends of the three cancer cell subpopulations (i.e., chemotherapy-resistant, vaccine-resistant, sensitive cells) under different stability conditions are illustrated in a simplex $ \{(P_C, P_I, P_S)|P_C+P_I+P_S = 1, \quad 0\le P_C, P_I, P_S\le 1\} $ which is presented in a triangle under trilinear coordinate ^[22,23](see Figure 1). Each point $ P $ in the simplex represents a particular structure of the population, $ (P_C, P_I, P_S) $. At three vertices $ S $, $ I $ and $ C $ of the simplex, $ (P_C, P_I, P_S) = (0, 0, 1) $, $ (0, 1, 0) $, $ (1, 0, 0) $, respectively (see Figure 1(A)), that is, only one subpopulation is present (with $ 100\% $) while the other two subpopulations become extinct in each case. The edges in the simplex are the sets of points where two subpopulations coexist and another one becomes extinct. For example, the coordinate of a point $ P $ in the edge $ \overline{SC} $ is $ (P_C, P_I, P_S) $ with $ P_I = 0 $, $ P_C > 0 $, $ P_S > 0 $, where $ P_C $ and $ P_S $ show the distances from $ P $ to the side line $ \overline{IS} $, $ \overline{CI} $, respectively. The trilinear coordinates of points in the other two edges are given similarly. The interior of the simplex is the set of points $ (P_C, P_I, P_S) $ with $ P_C, P_I, P_S > 0 $; the trilinear coordinates $ (P_C, P_I, P_S) $ of an inner point $ P $ in the simplex are given by the distances from $ P $ to the three edges; for instance, the distance from inner point $ P $ to the edge $ \overline{SC} $ represents the $ P_I $ value of the point $ P $ (see Figure 1(A)).

Figure 1(B)–(H) demonstrated the evolution dynamics of the three cancer cell subpopulations toward evolutionary stable points $ E_i $, $ i = 1, 2, ..7 $, respectively. Figure 1(B) showed the evolution toward evolutionary stable point $ E_1(0, 0, 1) $ under the condition $ f_S < \min\{f_C, f_I, 1\} $, where the parameter values were chosen as in Table 3. The arrows indicated the evolution trend and the red solid circles denoted the stable equilibria. Note that the population with initial sensitive cell proportion $ P_S > 0 $ would evolve toward the evolutionary stable point $ E_1(0, 0, 1) $, whereas if $ P_S = 0 $ initially, the population evolved toward the evolutionary stable points $ E_2(0, 1, 0) $, $ E_3(1, 0, 0) $ or $ E_5(\tilde{P_C}, \tilde{P_I}, 0) $ depending on the parameter values and the conditions as in Table 2; in Figure 1(B), the chosen parameters satisfied $ f_C < \min\{f_I, 1\} $, so that the population with initial $ S $ cell proportion $ P_S = 0 $ evolved toward $ E_3(1, 0, 0) $ which was presented by the green solid square. In a similar way, the evolution toward evolutionary stable points $ E_i $, $ i = 2, ..., 7 $, were depicted in Figure 1(C)–(H), under different parameter values shown in Table 3 satisfying the conditions in Table 1, respectively. The red lines in Figure 1(G)–(H) presented the equilibria lines $ E_6(\hat P_C, 0, \hat P_S) $ with $ \hat P_C+\hat P_S = 1 $, and $ E_7(P_C^*, P_I^*, P_S^*) $ with $ P_I^* = \frac{(d_C+d_I)-(c_I+\beta d_C)}{d_I(1-\gamma)} $ and $ P_C^*+P_S^* = 1-P_I^* $, respectively, which are stable but not asymptotically stable.

Note that $ f_S: = d_C+d_I $ was the collective cost of chemotherapy ($ d_C $) and vaccine ($ d_I $) to sensitive cells as they encounter another sensitive cell or chemotherapy-resistance cells; $ f_S^*: = d_C+\gamma d_I $ indicated the collective cost of chemotherapy ($ d_C $) and vaccine ($ d_I $) to sensitive cells as they play with vaccine-resistant cells; $ f_C: = c_C+\alpha d_I $ denoted the collective cost with the resistant cost of chemotherapy-resistant cells ($ c_C $) and the cost of vaccine to chemotherapy-resistant cells ($ \alpha d_I $); $ f_C^*: = c_C+\gamma\alpha d_I $ was the collective cost with the resistant cost of chemotherapy-resistant cells ($ c_C $) and the cost of vaccine to chemotherapy-resistant cells ($ \alpha d_I $) as they encounter vaccine-resistant cells; $ f_I: = c_I+\beta d_C $ represented the collective cost with the resistant cost of vaccine resistant cells ($ c_I $) and the cost of chemotherapy to vaccine-resistant cells ($ \beta d_C $).

We found that if the collective cost for sensitive cells $ (f_S) $ is less than the collective cost for chemotherapy-resistant cells ($ f_C $) and the collective cost for vaccine-resistant cells ($ f_I $), and also less than one, that is, $ f_S < \min\{f_C, f_I, 1\} $, then the populations containing sensitive cells initially would evolve to the population with $ 100\% $ of sensitive cells (i.e., $ E_1(0, 0, 1) $, see Figure 1(B)); the populations absent of sensitive cells initially would evolve to population with $ 100\% $ of chemotherapy resistant cells (i.e., $ E_3(1, 0, 0) $), as the collective cost for chemotherapy-resistant cells ($ f_C $) is lower than one and also less than the collective cost for vaccine resistant cells ($ f_I $).

When only the chemotherapy is applied (i.e., $ d_C > 0 $, $ d_I = 0 $), the evolutionary stable points $ E_1(0, 0, 1) $, $ E_3(1, 0, 0) $ and $ E_6(\hat P_C, 0, \hat P_S) $ are possible to be stable; if the monotherapy of vaccine is implemented (i.e., $ d_C = 0 $, $ d_I > 0 $), the evolutionary stable points $ E_1(0, 0, 1) $, $ E_2(0, 1, 0) $ and $ E_4(0, \bar P_I, \bar P_S) $ are possible to be stable. We explored the effects of timing for the effectiveness of treatment when the drugs are administered sequentially as in ^[9]. Figure 3 showed the evolution of cancer cell population when vaccine is applied first followed by chemotherapy (Figure 3(A)–(B)) or when chemotherapy is applied first followed by vaccine (Figure 3(C)–(E)). We saw that if the drug is changed to another one, then the corresponding resistant strain would grow up again or the sensitive cells would grow up. As the vaccine is applied first and changed to chemotherapy later on, the cell population with initial distribution of large sensitive cell proportion and very small resistant cell proportions evolved toward the tumor entirely consisting of either sensitive cells or chemotherapy resistant cells, depending on the chemotherapy dosage; the sensitive cell population would grow fast under low dose of chemotherapy (see Figure 3(A)), while high dose of chemotherapy could keep the sensitive cells' growth under control (see Figure 3(B)). Under a very strict dosing schedule that satisfies $ d_C = c_C < 1 $, the sensitive cells and chemotherapy-resistant cells could coexist. Similarly, as chemotherapy is applied first and then changed to vaccine, the cell population evolved toward the tumor consists entirely of either sensitive cells (see Figure 3(C)) or vaccine resistant cells (see Figure 3(D)), or coexistence of them (see Figure 3(E)), relying on the dosage of vaccine. In a long run, it is impossible to maintain both the sensitive cells and vaccine resistant cells under control with the implementation of chemotherapy as shown in ^[9], which could only be obtained in a short time scale.

2.3.2.   Adaptive therapy schedules

We further investigated the design of adaptive therapy schedules with the chemotherapy and p53 vaccine, that is, a periodic treatment cycle that traps the tumor into periodic and repeatable temporal dynamics of tumor composition. Evolutionary velocity is introduced for designing adaptive therapy schedules ^[13]. A fast dynamic region is needed for control of high proportion of resistant cells, whereas a slow dynamic region may be beneficial to the treatment holidays ^[13].

The evolutionary velocity of the system (2.1) evolving toward different evolutionary stable points (same as in Figure 1) was demonstrated in Figure 2, where the total velocity was calculated by L2-Norm of vector $ \left(\frac{\text{d}P_C}{\text{d}t}, \frac{\text{d}P_I}{\text{d}t}, \frac{\text{d}P_S}{\text{d}t}\right) $ for the system (2.1). It can be observed that the evolution velocity near the evolutionary stable point is very low; the evolutionary velocity at an evolutionary stable point is the lowest (equal to 0). For example, near the vertex $ S $ in Figure 2(A), that is, near the evolutionary stable point $ E_1(0, 0, 1) $, the evolution velocity is very small; the evolution velocity is also very small near the vertices $ C $ and $ I $, where the initial sensitive cell proportion is zero. A low velocity indicates a slow change to the tumor composition. In a similar way, the evolution velocity for the case with stable $ E_i $, $ i = 2, 3, ..., 7 $, was given in Figure 2 (B)–(G), respectively. In order to control the tumor progression, keeping the velocity of evolution towards evolutionary stable points at a relatively low level would be beneficial to the design of adaptive therapy schedule, which means drug dosages should be changed to avoid high evolution velocity.

The relative subpopulation velocity for chemotherapy-resistant cells ($ C $), vaccine-resistant cells ($ I $) and sensitive cells ($ S $) were shown in Figure 4 under the condition for stability of $ E_i $, $ i = 1, 3, 4, 5 $, where positive (negative) velocity indicated growth (decline) of the corresponding subpopulation. In order to examine the slow velocity (close to zero) regions, that is, slow dynamics regions, we derived the relative evolution velocities of each subpopulations by rescaling the positive and negative velocities separately.

Under each combination treatment $ (d_C, d_I) $ for a sufficiently long time, the tumor composition $ (P_C, P_I, P_S) $ would evolve to one of the evolutionary stable points $ E_i $, $ i = 1, ..., 7 $, relying on the condition satisfied by $ d_C $ and $ d_I $ (see Table 4). Under the treatment strategy associated with $ E_3 $, the velocity of the chemotherapy-resistant cell population is positive for most of the state space (see Figure 4(A2)). To control the chemotherapy-resistant cell population, it is necessary to switch to a new treatment after some period of the treatment, for example, a treatment schedule associated with $ E_4 $, where the velocity of the chemotherapy-resistant cell population is negative for almost all the state space (see Figure 4(A3)). Further, under the treatment schedule associated with $ E_1 $, the velocity of the sensitive cell population is positive for most of the state space (see Figure 4(C1)), therefore, to control sensitive cell population, it is necessary to change to another new treatment schedule after some period of the treatment, for example, a treatment schedule associated with $ E_5 $, where the velocity of sensitive cell population is negative in some parts of the state space (see Figure 4(C4)). Under the treatment schedule associated with $ E_4 $, the vaccine-resistant cells have positive evolution velocities in most of the state space (see Figure 4(B3)), so it is also necessary to change to another treatment schedule after some period of the treatment to control the vaccine-resistant cells. It is possible to design a cycle of the treatment schedules which could maintain growth of the three subpopulations under control.

A design of a treatment schedule with sequential application of the three treatment strategies associated with ($ E_3 $, $ E_5 $, $ E_4 $) (see Table 4) according to the relative velocities was displayed in Figure 5(A)–(B). There existed at least one periodic cycle in the interior of the absorbing state region (region in orange) whose boundary is given by the trajectories (dodger blue lines with arrows) connecting the evolutionary stable points (red solid circles). Under the same treatment strategies ($ E_3 $, $ E_5 $, $ E_4 $), it is impossible for the tumor composition to go out of the region, once it is driven into the absorbing regions. There may be other periodic cycles resulting from sequential administration of treatment strategies ($ E_3 $, $ E_5 $, $ E_4 $). The outer rim of the absorbing region is the largest cycle under the treatment regimen. Similarly, the sequential implementation of the treatment strategies associated with ($ E_1 $, $ E_5 $, $ E_3 $) could also be scheduled so that it results in periodic cycle evolution of the tumor compositions (see Figure 5(C)–(D)). Under these treatment schedules, the proportion of the drug-resistant cells and the sensitive cells are all under control with periodic oscillations, and the cancer could be treated as a chronic disease.

3.   The auxiliary effect of sensitive cells on drug-resistant cells

Experimental studies on targeted therapy with BRAF, ALK, or EGFR kinase inhibitors showed that these targeted therapies induces the release of sceretomes by drug-sensitive cancer cells, which lead to establishment of a tumor microenvironment that supports the expansion of drug-resistant clones ^[28]. Applying the replicator dynamics similar to the system (2.1), we explored the influence of this supporting effect of sensitive cells to drug-resistant cells for the evolution of cancer cell populations under combination of two targeted therapies. Here, we considered the targeted therapy with kinase inhibitors (KI) such as vemurafenib and dabrafenib (BRAF inhibitors) ^[28], which we denoted as KI-1 and KI-2, respectively. Assuming the support of sensitive cells to drug-resistant cells under targeted therapy by reduction to the drug cost of KI-1-resistant cells and KI-2-resistant cells with factor $ \rho_1 $ and $ \rho_2 $, respectively, we have the following payoff matrix

where $ 0 < \rho_1, \rho_2 < 1 $. Similar to the system (2.1), $ c_1 $ and $ c_{2} $ are the cost of resistant to two targeted therapies, KI-1 and KI-2, respectively; $ d_1 $ and $ d_2 $ are the cost of KI-1 and KI-2 to sensitive cells; $ \alpha $ and $ \beta $ ($ \alpha, \beta > 1 $) represent the extra cost of KI-1-resistant and KI-2-resistant cells being subjected to KI-2 and KI-1 drugs, respectively. In this case, the fitness of the subpopulations reads

where $ \vec P = (P_1, P_2, P_3) $, and $ P_i $, $ i = 1, 2, 3 $, denote the the proportions of KI-1-resistant cells, KI-2-resistant cells and sensitive cancer cells, respectively, which satisfy $ P_1+ P_2+P_3 = 1 $. The corresponding average fitness is

The evolution of the three subpopulations is given by the replicator dynamics:

In the following, we investigated the evolution dynamics of this replicator model.

3.1.   Local stability analysis

Similar to the system (2.1), the system (3.3) is defined on the simplex $ S_3 = \{(P_1, P_2, P_3)|\: P_1+P_2+P_3 = 1, \quad 0\le P_1, P_2, P_3\le1\} $, and the edges and the interior of the simplex are invariant, respectively. There are seven types of equilibria for the system under certain conditions, which are $ E_1(0, 0, 1) $, $ E_2(0, 1, 0) $, $ E_3(1, 0, 0) $, $ E_4(0, \bar{P_2}, \bar{P_3}) $, $ E_5(\tilde P_1, 0, \tilde P_3) $, $ E_6(\hat P_1, \hat P_2, 0) $ and $ E_7(P_1^*, P_2^*, P_3^*) $, where $ \bar{P_2} = \frac{(d_1+d_2)-(c_2+\rho_2\beta d_1)}{\beta(1-\rho_2)d_1} $, $ \bar{P_3} = 1- \bar{P_2} $; $ \tilde P_1 = \frac{(d_1+d_2) - (c_1+\rho_1\alpha d_2)}{\alpha(1-\rho_1)d_2} $, $ \tilde P_3 = 1-\tilde P_1 $; $ \hat P_1+\hat P_2 = 1 $; $ P_3^* = \frac{(c_1+\alpha d_2) -(d_1+d_2)}{\alpha(1-\rho_1)d_2} $, $ P_1^*+P_2^* = 1-P_3^* $.

The conditions for existence and local stability of these equilibria are given in Table 5, with the notations $ g_1: = c_1+\alpha d_2 $, $ g_2 : = c_2+\beta d_1 $, $ g_3: = d_1+d_2 $, $ g_1^*: = c_1+\alpha \rho_1 d_2 $ and $ g_2^*: = c_2+\beta \rho_2 d_1 $. Note that the existence and local stability of the equilibria $ E_i $, $ i = 1, ..., 7 $, are also mutually exclusive.

The detailed proof of the stability of these equilibria are shown in the following theorem.

Theorem 3.1. For the system (3.3),

(1) $ E_1(0, 0, 1) $ always exists and is LAS if $ d_1+d_2 < {\rm min}\{c_1+\rho_1\alpha d_2, c_2+ \rho_2\beta d_1, 1\} $;

(2) $ E_2(0, 1, 0) $ always exists; it is LAS if $ c_2+\beta d_1 < {\rm min}\{ c_1+\alpha d_2, d_1+ d_2, 1\} $;

(3) $ E_3(1, 0, 0) $ always exists; it is LAS if $ c_1+\alpha d_2 < {\rm min}\{ d_1+d_2, c_2+\beta d_1, 1\} $.

(4) If $ c_2+\beta d_1 > d_1+d_2 > c_2+\rho_2\beta d_1 $, the equilibrium $ E_4(0, \bar P_2, \bar P_3) $ exists, and furthermore if $ d_C+d_I < 1 $ and $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] < \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $, then $ E_4(0, \bar P_2, \bar P_3) $ is LAS, where $ \bar{P_2} = \frac{(d_1+d_2)-(c_2+\rho_2\beta d_1)}{\beta(1-\rho_2)d_1} $, $ \bar{P_3} = \frac{(c_2+\beta d_1) -(d_1+ d_2)}{\beta(1-\rho_2)d_1} $;

(5) If $ c_1+\alpha d_2 > d_1+d_2 > c_1+\rho_1\alpha d_2 $, the equilibrium $ E_5(\tilde P_1, 0, \tilde P_3) $ exits, if in addition, $ d_1 + d_2 < 1 $ and $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] > \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $, then $ E_5(\tilde P_1, 0, \tilde P_3) $ is LAS, where $ \tilde P_1 = \frac{(d_1+d_2) - (c_1+\rho_1\alpha d_2)}{\alpha(1-\rho_1)d_2} $, $ \tilde P_3 = \frac{(c_1+\alpha d_2)-(d_1+ d_2)}{\alpha(1-\rho_1)d_2} $;

Proof. (1) The Jacobian matrix of the system (3.3) at $ E_1(0, 0, 1) $ is given by

The corresponding eigenvalues are $ \lambda_1 = (d_1+d_2)-(c_1+\rho_1\alpha d_2), \lambda_2 = (d_{1}+d_{2})-(c_{2}+\rho_2\beta d_{1}), \lambda_3 = (d_{1}+d_{2})-1 $. Thus, if $ d_1+d_2 < {\rm min}\{c_1+\rho_1\alpha d_2, c_{2}+\rho_2\beta d_{1}, 1\} $, then $ \lambda_i < 0 $, $ i = 1, 2, 3 $, which implies that the equilibrium $ E_1(0, 0, 1) $ is locally asymptotically stable.

(2) The Jacobian matrix of the system (3.3) at the equilibrium $ E_2(0, 1, 0) $ is

The corresponding eigenvalues are $ \lambda_1 = (c_{2}+\beta d_{1})-(c_{1}+\alpha d_{2}), \lambda_2 = (c_{2}+\beta d_{1})-1, \lambda_3 = (c_{2}+\beta d_{1})-(d_{1}+d_{2}) $, which are negative if $ c_{2}+\beta d_{1} < {\rm min}\{ c_{1}+\alpha d_{2}, d_{1}+d_{2}, 1\} $. Hence, the equilibrium $ E_2(0, 1, 0) $ is locally asymptotically stable if $ c_{2}+\beta d_{1} < {\rm min}\{ c_{1}+\alpha d_{2}, d_{1}+d_{2}, 1\} $.

(3) The Jacobian matrix of the system (3.3) at $ E_3 (1, 0, 0) $ reads

The corresponding eigenvalues are $ \lambda_1 = (c_{1}+\alpha d_{2})-1, \lambda_2 = (c_{1}+\alpha d_{2})-(c_{2}+\beta d_{1}), \lambda_3 = (c_{1}+\alpha d_{2})-(d_{1}+d_{2}) $. We observe that $ c_{1}+\alpha d_{2} < {\rm min}\{ d_{1}+d_{2}, c_{2}+\beta d_{1}, 1\} $ implies $ \lambda_i < 0 $, $ i = 1, 2, 3 $, so that the equilibrium $ E_3(1, 0, 0) $ is locally asymptotically stable.

(4) The Jacobian matrix of the system (3.3) at $ E_4(0, \bar P_2, \bar P_3) $ is given by

where $ \overline{W}_1 = 1-c_{1}-\alpha d_{2} + \alpha (1-\rho_1)d_{2}\bar P_3 $ and $ \overline{W}_2 = \overline{W}_3 = 1-d_1-d_2 $. One of the corresponding eigenvalues is $ \lambda_1 = \overline{W}_1 - \overline{W}_2 = (d_1+d_2)- (c_{1}+\alpha d_{2}) + \alpha (1-\rho_1)d_{2} \frac{(c_2+\beta d_1)-(d_1+ d_2)}{\beta(1-\rho_2)d_1} $, which is negative if $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] < \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $. The other two eigenvalues ($ \lambda_2, \lambda_3 $) have negative real parts if the following conditions are satisfied:

which are equivalent to

respectively, noticing that $ \bar P_2+\bar P_3 = 1 $, $ \overline{W}_2 = \overline{W}_3 $, $ \bar P_2, \bar P_3, d_1 > 0 $ and $ 0 < \rho_2 < 1 $. Hence, $ \overline W_{3} = 1-(d_1+d_2) > 0 $ ensures the two eigenvalues $ \lambda_2 $ and $ \lambda_3 $ have negative real parts. Therefore, if $ d_1+d_2 < 1 $, $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] < \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $, and the equilibrium $ E_4(0, \bar P_I, \bar P_S) $ exists, then it is locally asymptotically stable.

(5) The Jacobian matrix of the system (3.3) at $ E_5(\tilde P_1, 0, \tilde P_3) $ reads

where $ \tilde W_1 = \tilde W_3 = 1-d_1-d_2 $ and $ \tilde W_2 = 1-c_{2}-\beta d_{1} + (1-\rho_2)\beta d_{1} \frac{(c_1+\alpha d_2)-(d_1+ d_2)}{\alpha(1-\rho_1)d_2} $.

We obtain that one corresponding eigenvalue is $ \lambda_1 = \tilde W_{2}-\tilde W_1 = (d_1+d_2)- (c_{2} + \beta d_1)+\beta (1-\rho_2)d_{1} \frac{(c_1+\alpha d_2)-(d_1+ d_2)}{\alpha(1-\rho_1)d_2} $, so that $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] > \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $ implies $ \lambda_1 < 0 $. The other two eigenvalues $ \lambda_2, \lambda_3 $ have negative real parts, if

which are equivalent to

respectively, where $ \tilde P_1, \tilde P_3, d_2 > 0 $ and $ 0 < \rho_1 < 1 $. Therefore, $ \tilde W_{1} = 1-d_1-d_2 > 0 $ guarantees $ {\rm Re}(\lambda_i) < 0 $, $ i = 2, 3 $. In summary, if the equilibrium $ E_5(\tilde P_1, 0, \tilde P_3) $ exits, and in addition, $ d_1+d_2 < 1 $ and $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] > \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $, then $ E_5(\tilde P_1, 0, \tilde P_3) $ is locally asymptotically stable.

Theorem 3.2. For the system (3.3),

(1) If and only if $ c_1+\alpha d_2 = c_2+\beta d_1 $, there exists a line of equilibria, $ E_6(\hat P_1, \hat P_2, 0) $, where $ \hat P_1+\hat P_2 = 1 $, and they are locally stable if $ c_1+\alpha d_2 < d_1+d_2 < 1 $.

(2) If $ c_1+\alpha \rho_1 d_2 < d_1+d_2 < c_1+\alpha d_2 $, $ c_2+\beta \rho_2 d_1 < d_1+d_2 < c_2+\beta d_1 $ and $ \alpha(1-\rho_1)d_2[(c_2+\beta d_1)-(d_1+d_2)] = \beta(1-\rho_2)d_1[(c_1+\alpha d_2)-(d_1+d_2)] $, there exists a line of equilibria $ E_7\left(P_1^*, P_2^*, P_3^*\right) $, where $ P_3^* = \frac{(c_1+\alpha d_2) -(d_1+d_2)}{\alpha(1-\rho_1)d_2} = \frac{(c_2+\beta d_1) -(d_1+d_2)}{\beta(1-\rho_2)d_1} $, $ P_1^*+P_2^* = 1-P_3^* $; they are locally stable if $ c_2+\beta d_1 < 1 $.

Proof. (1) The fitness of the subpopulations at the equilibria $ E_6(\hat P_1, \hat P_2, 0) $ are

Note that $ \hat W_1 = \hat W_2 $, if the equilibria $ E_6(\hat P_1, \hat P_2, 0) $ exist. The Jacobian matrix of the system (3.3) at $ E_6(\hat P_1, \hat P_2, 0) $ is given by

The corresponding eigenvalues are $ \lambda_1 = 0 $, $ \lambda_2 = \hat W_3- \hat W_1 = (c_{1}+\alpha d_{2})-(d_{1}+ d_{2}) $, $ \lambda_3 = -\hat W_1 = (c_{1}+\alpha d_{2})-1 $. Hence, if $ c_{1}+\alpha d_{2} < \min\{d_{1}+ d_{2} , 1\}$, then $ \lambda_2, \lambda_3 < 0 $, so that the equilibria $ E_6(\hat P_1, \hat P_2, 0) $ are locally stable when they exist.

(2) The Jacobian matrix of the system (3.3) at the equilibria $ E_7\left(P_1^*, P_2^*, P_3^*\right) $ is given by

where $ W^* = W_1^* = W_2^* = W_3^* = 1-c_{2}-\beta d_{1} $. Using the elementary transformation $ T^{-1}J_7T $, with

the matrix $ J_7 $ is transformed into

where $ a = \alpha (1-\rho_1)d_2P_1^*P_3^*+\beta (1-\rho_2)d_1P_2^*P_3^*-(P_1^*+P_2^*)W^* $. The two similar matrices $ J_7 $ and $ J^* $ have the same eigenvalues, one of which is $ \lambda_1 = 0 $. The real parts of the other two eigenvalues $ \lambda_2, \lambda_3 $ are negative as the following conditions hold:

which are equivalent to

respectively, where $ P_1^*, P_2^*, P_3^*, \alpha, \beta, d_1, d_2 > 0 $ and $ 0 < \rho_1, \rho_2 < 1 $. Therefore, $ W^* = 1 -c_{2}-\beta d_{1} > 0 $ permits the other two eigenvalues have negative real parts. As a result, the equilibria $ E_7\left(P_1^*, P_2^*, P_3^*\right) $ are locally stable, where $ P_3^* = \frac{(c_1+\alpha d_2) -(d_1+d_2)}{\alpha(1-\rho_1)d_2} = \frac{(c_2+\beta d_1) -(d_1+d_2)}{\beta(1-\rho_2)d_1} $, $ P_1^*+P_2^* = 1-P_3^* $.

3.2.   Influence of the supportive effect on tumor evolution

We investigated the effects of $ \rho_1 $ and $ \rho_2 $ on the evolution of cancer cell population. Figure 6 illustrated the evolution of cancer cell populations as the supports of sensitive cells to KI-1-resistant cells and KI-2-resistant cells increased, that is, as $ \rho_1 $ and $ \rho_2 $ decreased, respectively. Different shaded colors indicated evolution of cancer cell populations to different evolutionary stable points.

We chose a set of parameters so that the cell population evolves to the $ 100\% $ sensitive cell population (see Figure 6) as the supports of sensitive cells to KI-1-resistant cells and KI-2-resistant cells are very small (i.e., $ \rho_1 $ and $ \rho_2 $ are close to 1). As the supportive effect to KI-1-resistant cells increases, the KI-1-resistant cells grow up and the cell population evolves to the equilibrium point where sensitive cells and KI-1-resistant cells coexist (see the subplots in the first column of Figure 6). Similarly, as the supportive effect to KI-2-resistant cells increases, the KI-2-resistant cells establish and the cell population evolves to the coexistence of sensitive cells and KI-2-resistant cells (see the subplots in the first row of Figure 6). The three subpopulations could coexist (see subplots shaded with green in Figure 6) at very specific situations when $ \rho_1 $ and $ \rho_2 $ satisfy $ \rho_2 = 1-(1-\rho_1)[\alpha d_2(c_2+\beta d_1-d_1-d_2))]/[\beta d_1(c_1+\alpha d_2 -d_1-d_2)] $, $ \rho_1 < (d_1+d_2-c_1)/(\alpha d_2) $ and $ \rho_2 < (d_1+d_2-c_2)/(\beta d_1) $. The corresponding bifurcation diagrams are shown in Figure 7.

Figure 7(A)–(C) demonstrated the bifurcations with respect to the proportions of sensitive cells ($ P_3 $), KI-1-resistant cells ($ P_1 $) and KI-2-resistant cells ($ P_2 $), respectively. The color bars indicated the corresponding proportions of the three subpopulations. With $ \rho_1 $ and $ \rho_2 $ in the region $ S_1 $, that is, $ \rho_1 > (d_1+d_2-c_1)/(\alpha d_2) $ and $ \rho_2 > (d_1+d_2-c_2)/(\beta d_1) $, the cell population evolves to the evolutionary stable point $ E_1(0, 0, 1) $, so that sensitive cells establish. In the region $ S_2 $ with $ \rho_1 < (d_1+d_2-c_1)/(\alpha d_2) $ and $ \rho_2 > 1-(1-\rho_1)[\alpha d_2(c_2+\beta d_1-d_1-d_2))]/[\beta d_1(c_1+\alpha d_2 -d_1-d_2)] $, the KI-1-resistant cells establish and coexist with sensitive cells at $ E_5(\tilde P_1, 0, \tilde P_3) $. Similarly, the KI-1-resistant cell population evolves to coexists with sensitive cells as $ \rho_1 $ and $ \rho_2 $ are in the region $ S_3 $ where $ \rho_2 < (d_1+d_2-c_2)/(\beta d_1) $ and $ \rho_2 < 1-(1-\rho_1)[\alpha d_2(c_2+\beta d_1-d_1-d_2))]/[\beta d_1(c_1+\alpha d_2 -d_1-d_2)] $. On the slanting line segment $ \rho_2 = 1-(1-\rho_1)[\alpha d_2(c_2+\beta d_1-d_1-d_2))]/[\beta d_1(c_1+\alpha d_2 -d_1-d_2)] $, the three subpopulations evolve to coexist at $ E_7(P_1^*, P_2^*, P_3^*) $.

4.   Discussion and conclusions

The emergence and growth of drug-resistant tumor cell subpopulations is a major problem associated with cancer treatments. Drug resistance is often mitigated by application of combination therapies. In the present paper, we investigated the evolution of cancer cell populations under different schedules of combination treatment with chemotherapy and p53 vaccine, by constructing a nonlinear replicator dynamical system model of sensitive cancer cells, chemotherapy-resistant cells and p53 vaccine-resistant cells based on the payoff matrix proposed by Basanta D. et al. ^[9]. By analyzing the local stability of evolutionary stable points of the model and determining corresponding stability conditions, we demonstrated that cancer cell populations could evolve to a pure cancer cell population of sensitive cells (equilibrium $ E_1(0, 0, 1) $), chemotherapy-resistant cells (equilibrium $ E_3(1, 0, 0) $), or vaccine-resistant cells (equilibrium $ E_2(0, 1, 0) $), respectively, or coexistence of sensitive cells and p53 vaccine-resistant cells (equilibrium $ E_4(0, \bar P_I, \bar P_S) $), or coexistence of chemotherapy-resistant cells and p53 vaccine-resistant cells (equilibrium $ E_5(\tilde P_C, \tilde P_I, 0) $) under certain combination strategies of chemotherapy and p53 vaccine (see Table 4 and Figure 1). The coexistence of sensitive cells and chemotherapy-resistant cells (equilibria $ E_6(\hat P_C, 0, \hat P_S) $), or coexistence of the three subpopulations (equilibria $ E_7(P_C^*, P_I^*, P_S^*) $) could establish only under very strict conditions of combination therapy. Under monotherapy with chemotherapy, the population evolved to 100% of sensitive cells or 100% of chemotherapy-resistant cells (see Figure 3 (A)–(B)), while under monotherapy of p53 vaccine, the population evolved to 100% of sensitive cells or vaccine-resistant cells, or coexistence of them (see Figure 3(C)–(E)).

We further explored the design of adaptive therapy schedules with combination of chemotherapy and p53 vaccine, based on the evolutionary velocity and evolutionary absorbing regions under different treatment scenarios. We demonstrated the design of treatment schedules with sequential application of three treatment strategies associated with $ (E_3, E_5, E_4) $ or $ (E_1, E_5, E_3) $ which showed periodic cycle evolution the proportions of three tumor cell subpopulations being under control (see Figure 5).

The experimental research on targeted therapy revealed that sensitive cells can support the survival of drug-resistant cells by secreting cytokines. Through establishment a new replicator dynamical model based on the previous model, we explored the supportive effect of sensitive cells on two kind of BRAF inhibitor-resistant cells, that is, KI-1-resistant (vemurafenib-resistant) cells and KI-2-resistant (dabrafenib-resistant) cells. The influence of the supporting effects of sensitive cells on the evolutionary outcomes of the cancer cell populations are demonstrated by analysis of the existence and local asymptotic stability of the equilibrium points of the model. With an increase of the supporting effect to KI-1-resistant cells, the cell population could evolves from 100% of sensitive cells toward coexistence of sensitive cells and KI-1-resistant cells, and similarly with the improvement of the supporting effect to KI-2-resistant cells, cell populations evolved toward coexistence of sensitive cells and KI-2-resistant cells (see Figure 6). The evolutionary bifurcation diagram were shown with the the strength of supporting effects to KI-1-resistant cells and KI-2-resistant cells. At the bifurcation boundary, the three subpopulations could coexist (see Figure 7).

The evolution outcome of the three subpopulations of sensitive cells and two types of resistant cells, under corresponding combination therapy were explored in this paper by replicator dynamical modeling and theoretical analysis of the model. The conclusions about the evolution outcome and the design of the adaptive therapy schedules still need experimental and clinical verification.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 12171478).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.