Authors: Chau Hoang (1Department of Mathematical Sciences, New Mexico State University Las Cruces, New Mexico 88001, USA), Tuan Anh Phan (2Institute for Modeling Collaboration and Innovation, University of Idaho Moscow, Idaho 83844, USA), Cameron J. Turtle (3Sydney Medical School, Faculty of Medicine and Health, The University of Sydney Camperdown, NSW, 2006, Australia), Jianjun Paul Tian (1Department of Mathematical Sciences, New Mexico State University Las Cruces, New Mexico 88001, USA)
Categories: Article
Source: Mathematical biosciences
Authors: Chau Hoang, Tuan Anh Phan, Cameron J. Turtle, Jianjun Paul Tian
Based on a deterministic and stochastic process hybrid model, we use white noises to account for patient variabilities in treatment outcomes, use a hyperparameter to represent patient heterogeneity in a cohort, and construct a stochastic model in terms of Ito stochastic differential equations for tesing the efficacy of three different treatment protocols in CAR T cell therapy. The stochastic model has three ergodic invariant measures which correspond to three unstable equilibrium solutions of the deterministic system, while the ergodic invariant measures are attractors under some conditions for tumor growth. As the stable dynamics of the stochastic system reflects long-term outcomes of the therapy, the transient dynamics provide chances of cure in short-term. Two stopping times, the time to cure and time to progress, allow us to conduct numerical simulations with three different protocols of CAR T cell treatment through the transient dynamics of the stochastic model. The probability distributions of the time to cure and time to progress present outcome details of different protocols, which are significant for current clinical study of CAR T cell therapy.
Chimeric antigen receptor (CAR) T cell immunotherapy has been regarded as a major advance in the fight against cancers, especially those associated with the hematopoietic system [1]. This is a special adoptive cellular therapy in that T lymphocytes are taken from the patient’s blood, genetically modified to recognize specific antigens expressed by the tumor, expanded in vitro, and infused into the patient, usually after lymphodepletion chemotherapy. Due to the efficacy of CAR T cell immunotherapy, the Food and Drug Administration has approved 6 CAR T cell products for use in hematologic malignancies [2, 3]. However, temporary reductions in tumor burdens have been observed in patients treated with CAR T cells although overall response rate and complete response rate are high in many CAR T cell drug clinical trials [4]. This requires a complete understanding of the treatment process after CAR T cell infusion. In particular, the dynamics of CAR T cells after infusion, and the interactions among CAR T cells, cancer cells, and other relevant cells. Those are suitable for mathematical modeling once relevant key factors and interactions among components are identified. It is an opportunity for mathematical analysis and numerical simulations to make contributions to understanding CAR T cell immunotherapy.
There are several mathematical models about CAR T cell treatment at different levels. At the cellular and subcellular levels, mathematical models [5, 6, 7, 8] have been developed to better understand the downstream signaling mechanisms of CAR T-cell activation. These studies identified the competitive inhibition mechanism of CAR protein phosphorylation, the influence of co-stimulatory domain on overall CAR protein phosphorylation and how the cell-to-cell heterogeneity in terms of protein expression levels can influence the CAR T-cell activation response. These studies may provide insight into strategies to design novel CARs with optimized signaling. At the cell population level, mathematical models [9, 10, 11, 12, 13, 14, 15, 16] have been developed to analyze the complex interactions among the tumor, immune cells including CAR T cells, and microenvironmental factors. Each of these models has its specific focus and application in studies of the therapeutic outcomes at different patient conditions and doses. These models have limited input of clinical data. It is known that having adequate clinical data is a challenge for mathematical modeling, in general, in particular, for CAR T-cell modeling, since the field is still in its infancy and there is a lack of consensus on key factors driving therapeutic efficacy and safety. However, this provides an opportunity for mathematical analysis and numerical simulations to explore possible important factors at the cell population level.
There is a mathematical framework that describes dynamics and interactions among three cell populations [17], which are endogenous or normal T cells (N), CAR T cells (C), and tumor cells (B), in order to investigate the likelihood of successful treatment of CAR T cell therapy under different scenarios. The framework is a combination of a deterministic system of three ordinary differential equations and a birth-death stochastic process model in which the deterministic system is operated when the tumor cell population is above a certain threshold and the stochastic process is used to update the system when the tumor is below the threshold while the other populations are updated by the deterministic system in parallel. The deterministic system is given by (1)dNdt=rNNln(KNN+C),dCdt=rCCln(KCN+C),dBdt=rBB−γBBCkB+C.
To better fit the clinical data, Kimmel et al [17] calibrated the deterministic model (1) by considering the growth rate of the CAR T cell population to be functionally dependent on immune reconstitution reduction (2)rC:=rC(T)=ρC+b(N+C−KN)2a(N+C)2+(N+C−KN)2, where T=N+C is the total T cell count. The baseline net growth rate ρC sets the time scale at which birth and death events occur on average in the CAR T count. These birth and death events depend on signals provided by all T signalling inefficiency factor a in rC and immune reconstitution impact b in rC. The deterministic system (1) provides a general framework in which two types of T cells compete with each other. CAR T cells rapidly grow at net growth rate rC(T) but experience an emerging growth rate disadvantage compared to normal T cells growing at rate rN. This can be explained by a feature in immune reconstitution after both types of T cell proliferate, but only counts of normal T cells are supported by differentiation from stem and progenitor cells, which leads to the fact that the carrying capacity KC of CAR T cells is naturally lower than normal T cells’s carrying capacity KN. Meanwhile, the tumor cells grow autonomously at a net growth rate rB and are eradicated by CAR T cells at rate γB, proportional to the number of CAR T cells and inversely proportional to killing rate saturation parameter kB. The stochastic process model is constructed in terms of continuous time birth and death processes of three cell populations in which a system of master equations of transition probabilities among three states N, C, and B is developed (for details of this stochastic formulation, see Supplementary Material in [17]).
With the deterministic setting considered alone, there is no positive equilibrium state for the deterministic system (1) and all the equilibrium states on the boundary, which represent tumor eradication, are unstable. So the deterministic model alone cannot adequately capture the events that the tumor cell population can become arbitrarily small and the system spends a long period of time in a regime near tumor eradication, where random cell death events could lead to the elimination of the tumor. This is the reason why the authors proposed and analysed the above hybrid framework with the goal of improving the outcomes of CAR T cell therapy. By defining cure and progression as stochastic events, they pointed out that cure events occurred early and were narrowly distributed, while progression events occurred late and were more widely distributed in time.
The success of CAR T cell treatment may be subject to environmental noises and stochastic effects that come from complicated environments within the patient’s body, fluctuations in health conditions and ages among individual patients, and measurement errors in collecting data in patients. One of the patient’s internal environments that plays a crucial role in influencing the behavior and effectiveness of CAR T cells and the tumor killing rate of CAR effector is the hostile tumor microenvironment. Solid tumors often have a dense stroma that physically blocks CAR T cell infiltration and trafficking into the tumor [18]. The high interstitial pressure also prevents extravasation of CAR T cells from blood vessels into the tumor [19]. Besides the tumor microenvironment, CAR T cell proliferation can be weakened by the low arginine microenvironment, reducing their efficacy in clinical trials against haematological and solid malignancies [20]. Furthermore, modifying gut microbiome environment in patients may have a great effect on the success of CAR T cell treatment. Specifically, antibiotics, those disrupting gut bacteria like piperacillin/tazobactam, before CAR T cell treatment were associated with worse outcomes [21]. In addition to the patient’s internal environments, fluctuations in health conditions and ages among individual patients are also two key factors that limit the efficacy of CAR T cell treatment. For example, the autologous T cells are starting materials used to make CAR T cells from each individual patient. Differences in patients’ prior treatments, disease status, age, etc. affect the quality and composition of T cells in the starting materials, which can impact CAR T cell manufacturing [22]. It is clear that the potency and viability of final CAR T cell product in manufacturing heavily rely on data collection in patients, which is often subject to measurement errors.
These above stochastic factors may lead to variability in the outcomes of CAR T cell treatments in that the therapy works for some patients but not for others. To date, besides deterministic models, hybrid models, and stochastic process models mentioned above, there is no stochastic model using stochastic differential equations to quantify the events such as poor product quality, disease progression, death from toxicity in CAR T cell therapy. Stochastic differential equations (SDEs) are suitable for numerically simulating CAR T cell treatment efficacy and variability in clinical trials. This is because there are two possible interpretations for SDE modeling. First, each solution of a SDE model could represent the potential response to treatment for a hypothetical patient. When we run a SDE model for a cancer treatment 10000 times, for example, we could think of the 10000 solution paths as 10000 different hypothetical patients. The stochastic factors in the model may capture individual variations such as genetic differences, immune system responses, or other patient-specific characteristics. So a SDE model allows us to explore the variability in treatment outcomes across a cohort of hypothetical patients. Second, each solution path of a SDE model could represent a different possible trajectory for the treatment of a single hypothetical patient. By this way, if we run a SDE model 100 times, for instance, then we could treat the 100 solution paths as 100 treatment possibilities for one single patient. Stochasticity in the model may then come from sources of uncertainty in the treatment process such as variations in drug effectiveness. Thus, modeling a treatment therapy allows us to unravel the uncertainty and variability in the treatment process of a single hypothetical patient. In this paper, we will utilize the first interpretation for our stochastic modeling and discuss the second interpretation in the Discussion section.
Based on the hybrid model, we propose a system of Ito stochastic differential equations that can capture the dynamics of three given cell types in CAR T cell therapy. In our model, we use white noise to represent environmental noises and stochastic effects, and we incorporate white noise into some important parameters of the deterministic model. Suppose we would like to account for the variability of the parameter α. Then we replace α with α¯=α+τdWdt where W=W(t) is a standard Brownian Motion and τ stands for the noise intensity. Now, the new parameter α¯ is no longer a constant but a stochastic process. Within a sufficiently small interval [0, t], we can approximate dWdt by W(t)−W(0)t−0=W(t)t∼N(0,1t), which is the normal distribution with mean 0 and variance 1t. To count for patient heterogeneity in simulations, we utilize a hyperparameter σ in the following way. When we run our model, the value of α will be chosen from the normal distribution N(α,σ2α2), which is independent of N(0,1t). Since the value of W(t)t is chosen from N(0,1t), at each time t the value of α¯ will be chosen from N(α,σ2α2+τ2/t). For σ=0, we are able to generate a cohort of patients with the same health conditions (e.g. the same tumor growth rates). For σ>0, we can generate a cohort of diverse patients with different health conditions (e.g. different tumor growth rates or CAR T cell growth rates).
Based on our arguments above and our previous work [26, 27, 28], we incorporate white noise into the baseline growth rate of CAR T cells ρC, tumor growth rate rB, and the tumor-killing rate γB as follows. ρC→ρC+τ1dW1dt, rB→rB+τ2dW2dt, γB→γB+τ3dW3dt, where W1, W2 and W3 are mutually independent Wiener processes, τ1, τ2 and τ3 are intensities of corresponding noises. Then, a system of three Ito stochastic differential equations is obtained as follows, (3)dN=rNNln(KNN+C)dt,dC=rCCln(KCN+C)dt+τ1Cln(KCN+C)dW1,dB=(rBB−γBBCkB+C)dt+τ2BdW2−τ3BCkB+CdW3.
We conduct a detailed analysis of this system and a series of numerical simulations with different treatment protocols. As we know the deterministic system (1) has three equilibrium points (0,0,0), (0,KC,0), and (KN,0,0) on the boundary of the invariant domain. Since those equilibrium solutions are locally unstable, the deterministic model does not have a cure state built in. In our stochastic model (3), it has three ergodic invariant measures which correspond to the three equilibrium solutions of the deterministic model. When the tumor growth rate is smaller than a quantity of a combination of noise intensities and carrying capacities, the two invariant ergodic measures corresponding to equilibrium solutions (0,0,0) and (0,KC,0) are attractors. When the tumor growth rate is even smaller than its variance, the invariant ergodic measure corresponding to the equilibrium solution (KN,0,0) is a global attractor. These invariant ergodic measures may represent some cure states. However, when these conditions are not satisfied, we perform numerical simulations with three different protocols where simulated patient population has a size of 10000, (i) performing lymphodepletion and then one dose of CAR T cell infusion, (ii) performing lymphodepletion, and the first dose of CAR T cell infusion, then the second dose of CAR T cell infusion after 15 days, (iii) performing the first lymphodepletion, and the first dose of CAR T cell infusion; the second lymphodepletion after 10 days since the first dose, and the second dose of CAR T cell infusion at day 15. It should be noticed that second doses of CAR T cells is rarely used in clinical trials because of immune rejection. We define two stopping times (random variables), the time to cure - the first time that the tumor cells go below 1 cell, the time to progress - the first time that the tumor size go to 120% of its initial size. One set of numerical simulations shows the following results. In the first protocol (i), 47.87% of the patients get their tumor cell below 1 cell during 30 to 135 days after CAR T cell infusion while 52.13% of the patients get their tumor progression during 180 to 550 days after CAR T cell infusion. In the second protocol (ii), 47.28% of the patients get their tumor cell below 1 cell during 30 to 120 days after CAR T cell infusion while 52.72% of the patients get their tumor progression during 180 to 540 days after CAR T cell infusion. In the third protocol (iii), 48.33% of the patients get their tumor cell below 1 cell during 30 to 135 days after CAR T cell infusion while 51.67% of the patients get their tumor progression during 170 to 550 days after CAR T cell infusion. Effects of patient variability and patient heterogeneity on probability of cure and distribution of time to cure are tested using mutiple cohorts of 100 simulated patients with different values of hyperparameter σ and the noise intensity τ3. Another set of numerical simulations shows that, if the time between the first and second dose of CAR T cell infusions is longer, the probability of the tumor size goes below 1 cell will be greater with the second lymphodepletion; however, the probability of the tumor size goes below 1 cell is not increasing if the second lymphodepletion is not administered.
The rest of this article is organized as follows. In Section 2, we present main results of our study and give some medical interpretations. In Section 3, we give detailed proofs about our analytical results. In Section 4, we conduct numerical simulations to validate our analytical results and test outcomes of three different treatment protocols. We close our presentation with some discussion in Section 5.
We begin with ensuring that the noises are incorporated appropriately, which means the system possesses realistically reasonable solutions. Since all cell populations are always non-negative and lymphocyte counts cannot exceed their carrying capacities, we consider the invariant domain D={(N,C,B)⊤∣0≤N≤KN,0≤C≤KC,B≥0} for the stochastic system (3). Assume that we are working on a complete probability space (Ω, 𝓕, {𝓕t}t≥0, ℙ) with the filtration {𝓕t}t≥0 satisfying the usual condition. Let U(t)=(N(t),C(t),B(t))⊤, t≥0, be the process given by solution of the system (3). We use 𝓛 to denote the infinitesimal operator of the process U, and ℙu and Eu to denote the probability law and expectation on Ω, respectively, when the solution path starts at u=(n,c,b)⊤.
For simplicity, we let dW=(dW1,dW2,dW3)⊤ and denote the drift term and the diffusion term of the system (3) by f(U)=[rNNln(KNN+C)rCCln(KCN+C)rBB−γBBCkB+C],g(U)=[000τ1Cln(KCN+C)000τ2B−τ3BCkB+C], then, we can rewrite the system (3) into its matrix form dU=f(U)dt+g(U)dW.
The following theorem certifies the rationality of the SDE system (3) in which its solution globally stays in the domain D for any initial value in D.
For any initial value U(0)=u=(n,c,b)⊤∈ℝ+3:={(N,C,B)⊤∣N≥0,C≥0,B≥0} there exists almost surely continuous solution U(t) to the system (3) that remains in ℝ+3 for all time t≥0, and U(t) is a strong Markov process that satisfies Feller property. Furthermore, if U(0)∈D∘={(N,C,B)⊤∣0<N<KN,0<C<KC,B>0}. then U(t)∈D∘ for all t≥0 almost surely.
Since the ODE system (1) has 3 equilibrium points (0,0,0)⊤, (0,KC,0)⊤, (KN,0,0)⊤ that all stay on the boundary of the domain D, to investigate the cell dynamics in long-term treatment, we perform boundary analysis for the system (3) in Section 3. Note that both types of T cells proliferate and expand, but only normal T cells get supported in reconstitution from stem cells. As a result, CAR T cells are outcompeted by normal T cells in the long term, resulting in lower adapting ability. Therefore, we proceed our analysis under the assumption of KN>KC Besides, it can be seen that the tumor cell population completely depends on the population size of CAR T cells, while normal T cells and CAR T cells compete each other; thus, we can focus on investigating cell kinetics from the first two equations of the system (3), and then apply the results to the last equation in (3) to derive the dynamics of tumor cell population.
When the initial total lymphocyte count is zero, that is, N(0)=0 and C(0)=0, the solution U(t) of the system (3) approaches a.s. (0,0,0) provided that rB<12τ22. This means that there is a unique ergodic invariant probability measure μ0=δ0∗×δ0∗×δ0∗ for the system (3) on the boundary ∂D if rB<12τ22.
When there are only CAR T cells in the initial lymphocyte count, i.e. N(0)=0 and 0<C(0)<KC, the system (3) possesses two ergodic invariant probability measures μ0 and μ1=δ0∗×δKC∗×δ0∗ on the boundary ∂D in which the solution U(t) either converges weakly to μ0 with positive probability or converges weakly to μ1 with positive probability provided rB<γBKCkB+KC+12τ22+12τ32KC2(kB+KC)2.
When there are only initially normal T cells in the lymphocyte count, i.e. 0<N(0)<KN and C(0)=0, the solution U(t) converges weakly to the unique ergodic invariant probability measure μ2=δKN∗×δ0∗×δ0∗ provided that rB<12τ22.
When there are both normal T cells and CAR T cells in initial lymphocyte count, i.e. 0<N(0)<KN and 0<C(0)<KC, assuming that rB<12τ22, then three ergodic invariant probability measures on the boundary μ0, μ1, μ2 exist. Both μ0 and μ1 are repellers while μ2 is a local attractor. The following theorem guarantees that μ2 is actually a global attractor and, as a consequence, there is no positive invariant probability measure for the system (3) in D∘.
Assume KN>KC, for any initial value (n,c)⊤∈(0,KN)×(0,KC), the solution (N(t),C(t))⊤ to the subsystem consisting of the first two equations of (3) converges a.s. to (KN,0)⊤.
For the study of numerical simulations, we define two random variables, the time to cure - the first time that the tumor cells go below 1 cell, the time to progress - the first time that tumor size progresses to 120% of its initial size. We numerically simulate the protocols of CAR T cell treatments in the current practice and summarize the results of numerical simulations as follows.
Using parameter values estimated from clinical data [17], the model predicts long-term and transient dynamics of CAR T cells, endogenous T cells, and tumor cells. With 10000 simulated patients, three different protocols of the treatments have similar outcomes, 47.28% to 48.33% of the patients get their tumor reduced to below 1 cell first time during 30 to 135 days after CAR T cell infusion, while 51.67% to 52.72% of the patients have their tumor progressed to 120% of the initial size first time during 170 to 550 days after CAR T cell infusion; the time between the first and second dose of CAR T cell infusion positively impacts the probability of cure.
Our stochastic model has three ergodic invariant probability measures, μ0=δ0∗×δ0∗×δ0∗, μ1=δ0∗×δKC∗×δ0∗, and μ2=δKN∗×δ0∗×δ0∗, which correspond to equilibrium solutions of the deterministic system. For the deterministic system, the three equilibrium solutions are unstable while for our stochastic system, there are conditions under which those invariant probabilities are attractors. There are two critical values for the tumor growth rate, denoted by R1=12τ22 and R2=γBKCkB+KC+12τ22+12τ32KC2(kB+KC)2. When R1<rB<R2, μ0 and μ1 are local attractors. When rB<R1, μ0 and μ1 are repellers, and μ2 is a global attractor. It is interesting that the long-term behavior of the solutions is determined by the tumor growth rate, and two critical values of the tumor growth rate are given by its noise intensity, the CAR T cell carrying capacity, CAR T cell king rate and its noise intensity. If the tumor grows too slow that its growth rate is small comparing with its noise, then endogenous T cells take over and the tumor is eradicated. If the tumor growth rate is smaller than a combination of its noise intensity, CAR T cell carrying capacity, and CAR T cell killing rate, then CAR T cell may take over and the tumor is eradicated. This case seems reasonable since the tumor growth is controlled by CAR T cells. However, if the tumor growth rate does not satisfy these conditions, the system demonstrates transient behaviors. The transient dynamics is also of significance, which is used to numerically simulate different treatment protocols.
The patient variability and patient heterogeneity are encoded as noises and a hyperparameter, respectively, in our stochastic model. It should be noticed that noises, or patient variability, account for the deviation of observations in clinical dataset from the central tendency (such as mean or median) while the hyperparameter, or patient heterogeneity, refers to the diversity or differences within a clinical dataset. The results from numerical simulations when changing noises are different from those when changing the hyperparameter (see Figure 5). Those results do not follow any pattern, which are totally different from results of the hybrid model in [17] where increases in patient heterogeneity will drop probability of cure but drastically increase the range of possible times to cure.
In this section, we give a detailed proof of Theorem 2.1 in subsection 3.1. Then we conduct a detailed boundary analysis for the system (3) in subsection 3.2. Finally, proof of Theorem 2.2 is given in subsection 3.3 by using control theory and support theorem of diffusion processes.
Since the drift term f(U) and the diffusion term g(U) of the system (3) are continuously differentiable and hence locally Lipschitz continuous functions on {(N,C)⊤∣N+C>0,−kB<C<∞}×ℝ, there exists a unique locally almost surely continuous solution U(t) to system (3) up to the explosion time. The explosion time is defined as follows τe:=inf{t>0|min{N(t),C(t),N(t)ln[KNN(t)+C(t)],C(t)ln[KCN(t)+C(t)],C(t)kB+C(t),B(t)}=−∞ ormax{N(t),C(t),N(t)ln[KNN(t)+C(t)],C(t)ln[KCN(t)+C(t)],C(t)kB+C(t),B(t)}=∞}.
However, notice that C(t)ln[KCN(t)+C(t)]→−∞ happens if and only if N(t)+C(t)→0 and −∞<C(t)<0. It implies that N(t)+C(t)→0 and 0<N(t)<∞, i.e. N(t)ln[KNN(t)+C(t)]→∞. ln Following similar argument, it can be checked that C(t)ln[KCN(t)+C(t)]→∞ occurs at the same time with N(t)ln[KNN(t)+C(t)]→−∞. Therefore, we can simplify the definition of the explosion time to τe:=inf{t>0|min{N(t),C(t),C(t)ln[KCN(t)+C(t)],C(t)kB+C(t),B(t)}=−∞ ormax{N(t),C(t),C(t)ln[KCN(t)+C(t)],C(t)kB+C(t),B(t)}=∞}.
Also, the local Lipschitz continuity of f(U) and g(U) implies that the solution U(t), t∈(0,τe) is a strong Markov process that satisfies Feller properties. Now, we will first show that if the process starts with any nonnegative initial value, then it remains nonnegative almost surely for up to the explosion time. We denote u:=U(0)=(n,c,b)⊤ as the initial values of U(t). By the first equation of (3), we can solve for N(t)=nexp{∫0trNln[KNN(s)+C(s)]ds}a.s.
The second equation of (3) implies (4)C(t)=cexp{∫0t[rCln[KNN(s)+C(s)]−12τ12(ln[KCN(s)+C(s)])2]ds +τ1∫0tln[KCN(s)+C(s)]dW1(s)}a.s.
And the last equation of (3) follows B(t)=bexp{∫0t[rB−γBC(s)kB+C(s)−12τ22−12τ32(C(s)kB+C(s))2]ds +τ2W2(t)−τ3∫0tC(s)kB+C(s)dW3(s)}a.s.
It is easy to see that for all t∈(0,τe), if n, c, and b are zero, so are N(t), C(t), and B(t) correspondingly. And if these initial values n, c, and b are positive, then N(t)>0, C(t)>0, and B(t)>0 almost surely. Note that the scenario of n=0 and c=0 can happen at the same initial time, because (n,c)→(0,0) in ℝ+2 is equivalent to n+c→0, we can check that lim(n,c)→(0,0)cln(KCn+c)=0 and lim(n,c)→(0,0)nln(KNn+c)=0, which implies that these terms cannot blow up. From all arguments above, we conclude that if (n,c,b)⊤∈ℝ+3 then (N(t),C(t),B(t))⊤∈ℝ+3 for all t∈(0,τe).
Next, we show that if u∈D∘, then U(t)∈D∘ for all t∈(0,τe). Consider equation dN˜=rNN˜ln(KNN˜)dt with initial value N˜(0)=n=N(0)∈(0,KN). We can solve explicitly this equation for N˜(t)=KN(KNn)−e−rNt<KN.
Moreover, by first equation of (3) and the comparison theorem for ODE, we have 0<N(t)≤N˜(t)<KN for all t∈(0,τe) almost surely. Now consider the function with 0<c<KC V1(n,c,b)=KC−c−ln(KC−c).
Note that ∂V1∂c=−1+1KC−c and ∂2V1∂c2=1(KC−c)2. Then, by Ito’s formula, we have 𝓛V1=(−1+1KC−c)rCcln(KCn+c)+12τ12c2ln2(KCn+c)1(KC−c)2 =−rCcln(KCn+c)+rCclnKC−ln(n+c)KC−c+12τ12c2ln2(KCn+c)1(KC−c)2 ≤−rCcln(KCn+c)+rCclnKC−lncKC−c+k112τ12c2(lnKC−lncKC−c)2, where k1>0 is a constant satisfying ln2(KCn+c)≤k1(lnKC−lnc)2. Recall that rC=ρC+b(n+c−KN)2a(n+c)2+(n+c−KN)2; thus, for c∈(0,KC), we have rC≤ρC+b. Moreover, by Mean Value Theorem, there exists a value c¯∈(c,KC) that satisfies KC−lnc=1c¯(KC−c). It follows that 𝓛V1≤(ρC+b)KCln(KN+KCKC)+(ρC+b)cc¯+12k1τ12(cc¯)2 ≤(ρC+b)KCln(KN+KCKC)+ρC+b+12k1τ12.
Here, let us have some quick discussions about continuous function f(u)=u−1−lnu with u>0. It can be checked that on (0,∞), the function f(u) gets its minimum value of 0 when u=1; hence, u−lnu≥1 for all u>0. Following from this calculus result, for c∈(0,KC), we have V1(n,c,b)≥1. Furthermore, if we let K1=(ρC+b)KCln(KN+KCKC)+ρC+b+12k1τ12, then it yields 𝓛V1≤K1V1.
On the other hand, consider ξk=inf{t∈(0,τe)∣V1(U(t))≥k}, k∈ℕ. By applying Ito’s formula to V1(U(ξk∧t)), for fixed t∈(0,τe), we have EuV1(t):=EuV1(U(ξk∧t))=V1(u)+Eu∫0ξk∧t𝓛V1(U(s))ds ≤V1(u)+K1∫0ξk∧tEuV1(U(s))ds.
Since ξk∧s<ξk for all s∈(0,ξk∧t), V1(U(ξk∧s))<k<∞ which implies that ln(KC−C(ξk∧s))<∞ and so C(ξk∧s)<KC for all s∈(0,ξk∧t). By above argument, for all s∈(0,ξk∧t) 𝓛V1(U(ξk∧s))≤K1V1(U(ξk∧s)).
So EuV1(U(ξk∧t))≤V1(u)+K1∫0tEuV1(U(ξk∧s))ds.
By Gronwall’s inequality (see Theorem 8.1 p.45 in [23]), we have EuV1(t)≤V1(u)eK1t. Moreover, EuV1(t)=∫ΩV1(U(ξk∧t,ω)dℙu(ω)≥∫{ξk≤t}V1(U(ξk))dℙu≥kℙu{ξk≤t}.
Therefore, for any k∈ℕ, ℙ{ξk≤t}≤V1(u)eK1tk. It follows that ℙ{ξk>t}≥1−V1(u)eK1tk.
Furthermore, since ξk>t implies V1(U(s))<k for all s∈(0,t), it implies that ℙu{V1(U(s))<k,s∈(0,t)}≥ℙu{ξk>t}≥1−V1(u)eK1tk.
This result holds for all k≥1. Thus, we can let k→∞ and obtain ℙu{V1(U(s))<∞,s∈(0,t)}=1.
Since V1(U(s))<∞ implies C(s)<KC for all s∈(0,t), it follows that ℙu{C(s)<KC,s∈(0,t)}=1.
Note that we have worked with arbitrary t∈(0,τe); thus, it yields C(t)<KC for all t∈(0,τe) almost surely.
Finally, to complete this proof, it remains to show that τe=∞. Since we showed that C and N are bounded, only the explosion of B does matter. Thus, we can proceed with claiming that τe=τ∞:=inf{t∈(0,τe)∣B(t)=∞}=∞. We consider V2(N,C,B)=ln(1+B)on(0,KN)×(0,KC)×(0,∞).
By using Ito’s formula, we have 𝓛V2(t):=𝓛V2(U(t))=(rB−γBC(t)kB+C(t))B(t)1+B(t)−12τ22(B(t)1+B(t))2 −12τ32(B(t)1+B(t))2(c(t)kB+C(t))2 ≤(rB−γBC(t)kB+C(t))B(t)1+B(t) ≤K2,forsomepositiveconstantK2.
Let τk=inf{t∈(0,τe)∣B(t)>k}, it is clearly that τk increases to τ∞. Recall that we proved that if b=B(0)≥0 then B(t)≥0 for t∈(0,τe) almost surely. Therefore, for B(0)≥0, we have τe=τ∞. Now, we fix t>0 and apply Ito’s formula to get EuV2(τk∧t):=EuV2(U(τk∧t))=V2(u)+Eu∫0τk∧t𝓛V2(U(s))ds ≤K3+K2(τk∧t)(supposeV2(u)=K3) ≤K3+K2t.
Moreover, EuV2(τk∧t)≥∫{τk≤t}V2(U(τk))dℙu≥ln(1+k)ℙ{τk<t}. Thus, as k→∞ ℙu{τk∣<t}≤EuV2(τk∧t)ln(1+k)≤K3+K2tln(1+k)→0almostsurely.
Since it works for arbitrary t>0, we can conclude that ℙu{τ∞<∞}=0. Thus, τe=τ∞=∞ almost surely. This completes the proof.
We consider three cases for boundary analysis.
Suppose N(0)=0, C(0)=0. By using the results proved in the theorem 2.1, it follows that N(t)=0 and C(t)=0 for all t>0 almost surely. As a result, the last equation of (3) implies that dB=rBBdt+τ2BdW2. Apply Ito’s formula, we obtain B(t)=B(0)exp{(rB−12τ22)t+t2W2(t)}.
If rB<12τ22 then limt→∞(rB−12τ22+τ2W2(t)t)t=−∞. This leads to limt→∞B(t)=0 almost surely. Otherwise, B(t)→∞ as t→∞. Therefore, μ0:=δ0∗×δ0∗×δ0∗ is an ergodic invariant probability measure for the system (3) provided rB<12τ22.
Suppose N(0)=0 and 0<C(0)<KC. It follows that N(t)=0 almost surely and the second equation of (3) yields dC=(ρC+b(C−KN)2aC2+(C−KN)2)Cln(KCC)dt+τ1Cln(KCC)dW1.
For α∈(0,KC) fixed, we consider S(C)=∫αCexp{−∫αy2[ρC+b(μ−KN)2aμ2+(μ−KN)2]μlnKCμdμ/τ12μ2ln2KCμ}dy =∫αCexp{−∫αy2τ12[ρC+b(μ−KN)2aμ2+(μ−KN)2]⋅1μ⋅1lnKC−lnμdμ}dy
We can observe the behavior of S(C) when C is getting close to either 0 or KC. First, if C is in a neighborhood of KC such that α<C<KC, then the dummy variable y in the right-hand side of S(C) will be in (α, KC). So, for μ∈(α,y), since ρC+b(μ−KN)2aμ2+(μ−KN)2≥ρC, we obtain S(C)≤∫αCexp{−∫αy2ρCτ12⋅1μ⋅1lnKC−lnμdμ}dy=C1∫αC(lnKC−lnμ)2ρC/τ12dμ for some constant C1>0, which implies that limC→KCS(C)<∞. Second, if C is in a neighborhood of 0 such that 0<C<α, then the dummy variable y in the right-hand side of S(C) is in (0, α). Hence, for μ∈(y,α), since ρC+b(μ−KN)2aμ2+(μ−KN)2≤ρC+b, we get for some constant C2>0 S(C)≥∫αCexp{−∫αy2(ρC+b)τ12⋅1μ⋅1lnKC−lnμdμ}dy=C2∫αC(lnKC−lnμ)2(ρC+b)/τ12dμ
To obtain the behavior of the last integral above as C is close to 0, we need to compare the integrand (lnKC−lnμ)2(ρC+b)/τ12 with y1/2 when y is getting close to 0. So we look at the following limit limy→0+ln[y1/2(lnKC−lny)2(ρC+b)τ12]=limy→0+[12lny+2(ρC+b)τ12ln(lnKC−lny)] =limy→0+[12+2(ρC+b)τ12⋅ln(lnKC−lny)lny]lny.
By L’Hospital rule, we get limy→0+ln(lnKC−lny)lny=limy→0+−1y(lnKC−lny)1y=limy→0+−1lnKC−lny=0, it follows that lim limy→0+ln[y1/2(lnKC−lny)2(ρC+b)τ12]=−∞, i.e. limy→0+y1/2(lnKC−lny)2(ρC+b)τ12=0.
Thus, there exists 0<δ<α such that for any 0<y<δ, we have (lnKC−lny)2(ρC+b)τ12<1y1/2. Hence, for any 0<C<δ, as C→0+, S(C)=C2∫αC(lnKC−lny)2(ρC+b)τ12dy≥C2∫αC1ydy=2C2(C−α)→−2C2α.
Therefore, limC→0+S(C)>−∞. Then, by item 3 of Theorem 3.1 p.447 in [24], we can conclude that limt→∞C(t) exists a.s. and ℙ{limt→∞C(t)=0}+ℙ{limt→∞C(t)=KC}=1.
Hence, the solution to the equation (4) is not recurrent. It may either approach 0 or approach KC with positive probability. If the solution C(t) of (4) converges to 0 then the same argument as in Case 1 implies that the system (3) has the ergodic invariant probability measure μ0 provided rB<12τ22. When the solution C(t) of (4) approaches KC then B(t) behaves the same as the long-term behavior of the solution to the equation dB=(rBB−γBBKCkB+KC)dt+τ2BdW2−τ3BKCkB+KCdW3, which follows that B(t)=B(0)exp{[rB−γBBKCkB+KC−12τ22−12τ32(KCkB+KC)2]t +τ2W2(t)−τ3KCkB+KCW3(t)}.
If rB≤γBBKCkB+KC+12τ22+12τ32(KCkB+KC)2, then B(t)→0 as t→∞. Otherwise, B(t)→∞ as t→∞. Thus the system (3) possesses the ergodic invariant probability measure μ1:=δ0∗×δKC∗×δ0∗ provided rB≤γBBKCkB+KC+12τ22+12τ32(KCkB+KC)2.
Suppose C(0)=0 and 0<N(0)<KN. It follows that C(t)=0 almost surely, and the first equation of (3) becomes dN=rNNln(KNN)dt which yields N(t)=KN(N(0)KN)e−rNt→KN as t→∞. Then the long-term behavior of B(t) is the same as the behavior of the solution to equation dB=rBBdt+τ2BdW2, and so B(t)=B(0)exp{(rB−12τ22)t+τ2W2(t)}. If rB<12τ22, then limt→∞B(t)=0 almost surely. Otherwise, B(t)→∞ as t→∞. Hence, we have μ2:=δKN∗×δ0∗×δ0∗ is an ergodic invariant probability measure provided rB<12τ22.
From now on, we assume that 0<N(0)<KN, 0<C(0)<KC, and B(0)>0. Furthermore, suppose that KN>KC and rB<12τ22. Then all three invariant probability measures on the boundary μ0, μ1, and μ2 exist. Before giving a detailed proof of Theorem 2.2, we prove the following lemma.
Both measures μ0 and μ1 are repellers, while μ2 is a local attractor.
Proof. Since N(0)>0 and C(0)>0, it follows that N(t)>0 and C(t)>0 for all t>0 almost surely. We assume that U(t)=(N(t),C(t),B(t))⊤ stays in a small neighborhood of (0,0,0)⊤ for a long time. Since μ0 is ergodic, by strong law of large number, we can compute the Lyapunov exponents of μ0 along the solution components λ1(μ0)=limt→∞1tlnN(t)=limt→∞1t∫0trNln[KNN(s)+C(s)]ds=∫∂DrNln(KNN+C)μ0(dNdCdB).
As N+C is really small and close to 0, ln(KNN+C)>0 almost surely. Therefore, λ1(μ0)>0 almost surely. Similarly, λ2(μ0)=∫∂D[rCln(KCN+C)−12τ12ln2(KCN+C)]δ0∗(dN)δ0∗(dC)δ0∗(dB)<0, λ3(μ0)=∫∂D[rB−γBCkB+C−12τ22−12τ32(CkB+KC)2]δ0∗(dN)δ0∗(dC)δ0∗(dB) =rB−γB0kB+0−12τ22−12τ32(0kB+0)2=rB−12τ22<0.
Thus μ0 is a repeller. Follow the same fashion, we can check the Lyaponov exponents of μ1 along the solution components λ1(μ1)=∫∂DrNln(KNN+C)δ0∗(dN)δKC∗(dC)δ0∗(dB)=rNlnKNKC>0, λ2(μ1)=∫∂D[rCln(KCN+C)−12τ12ln2(KCN+C)]δ0∗(dN)δKC∗(dC)δ0∗(dB) =[ρC+b(KC−KN)2aKC2+(KC−KN)2]ln(KC0+KC)−12τ12ln2(KC0+KC)=0, λ3(μ1)=∫∂D[rB−γBCkB+C−12τ22−12τ32(CkB+C)2]δ0∗(dN)δKC∗(dC)δ0∗(dB) =rB−γBKCkB+KC−12τ22−12τ32(KCkB+KC)2<0.
Hence, μ1 is a repller. Finally, the Lyaponov exponents of μ2 along the solution components are λ1(μ2)=∫∂DrNln(KNN+C)δKN∗(dN)δ0∗(dC)δ0∗(dB)=rNln(KNKN+0)=0, λ2(μ2)=∫∂D[rCln(KCN+C)−12τ12ln2(KCN+C)]δKN∗(dN)δ0∗(dC)δ0∗(dB) =ρCln(KCKN+0)−12τ12ln2(KCKN+0)=ρCln(KCKN)−12τ12ln2(KCKN)<0, λ3(μ2)=∫∂D[rB−γBCkB+C−12τ22−12τ32(CkB+C)2]δKN∗(dN)δ0∗(dC)δ0∗(dB) =rB−γB0kB+0−12τ22−12τ32(0kB+0)2=rB−12τ22<0.
Hence, μ2 is a local attractor. □
Now we give a full proof of Theorem 2.2.
Proof of Theorem 2.2. We restrict ourselves to the system of the first two equations of (3) (5)dN=rNNln(KNN+C)dt,dC=[ρC+b(N+C−KN)2a(N+C)2+(N+C−KN)2]Cln(KCN+C)dt+τ1Cln(KCN+C)dW1, and show that its equilibrium point (KN,0)⊤ is globally stochastically asymptotically stable. We define V(t,v)=(N(t,v),C(t,v))⊤ is the solution to the system (5) with initial value v=(n,c)⊤. First of all, we prove that as the solution V(t,v) starts close to (KN,0)⊤ it will approach (KN,0)⊤ with large probability, i.e. (KN,0)⊤ is a locally stochastically asymptotically stable. Indeed, by computing partial derivatives of f1(n,c):=rNnln(KNn+c) at the equilibrium point (KN,0)⊤, we get ∂f1∂n|(KN,0)⊤=(rNln(KNn+c)−rNnn+c)|(KN,0)⊤=−rN, ∂f1∂c|(KN,0)⊤=(−rNnn+c)|(KN,0)⊤=−rN.
Thus, since f1(KN,0)=0, the Taylor expansion of f1(n,c) in a small neighborhood of (KN,0)⊤ is f1(n,c)=−rN(n−KN)−rNc+o((n−KN)2+c2), where lim(n,c)⊤→(KN,0)⊤o((n−KN)2+c2)(n−KN)2+c2=0. On the other hand, by proof of Lemma 3.1, λ2(μ2)=ρCln(KCKN)−12τ12ln2(KCKN)<0, so we can find a sufficiently small p>0 such that θ1:=ρCln(KCKN)−12τ12(1−p)ln2(KCKN)<0.
We consider the function V3(n,c)=(n−KN)2+cp. It can be seen that V3 is a positive-definite decrescent function on (0,KN)×(0,KC). Moreover, by Ito’s formula, 𝓛V3=2(n−KN)[−rN(n−KN)−rNc+o((n−KN)2+c2)] +pcp[ρC+b(n+c−KN)2a(n+c)2+(n+c−KN)2]ln(KCn+c)−12τ12p(1−p)cpln2(KCn+c) ≤−2rN(n−KN)2−2rNc(n−KN)+o((n−KN)2+c2) +pcp[ρCln(KCKN)−12τ12(1−p)ln2(KCKN) +ln(n+cKN)(12τ12(1−p)ln(KC2KN(n+c))−ρC) +b(n+c−KN)2a(n+c)2+(n+c−KN)2ln(KCn+c)
Let θ2(n,c)=ln(n+cKN)(12τ12(1−p)ln(KC2KN(n+c))−ρC) +b(n+c−KN)2a(n+c)2+(n+c−KN)2ln(KCn+c).
It can be checked that lim(n,c)⊤→(KN,0)⊤θ2(n,c)=0. Thus, for a ϵ1>0 small enough such that θ1+ϵ1<0, there exists a δ1>0 such that (n,c)⊤∈(KN−δ1,KN)×(0,δ1) implies θ2(n,c)<ϵ1. It follows that θ1+θ2(n,c)<θ1+ϵ1<0. On the other hand, choose ϵ2>0 such that ϵ2<min{2rN,−θ1−ϵ1}. Notice that lim(n,c)⊤→(KN,0)⊤o((n−KN)2+c2)(n−KN)2+c2=0; thus, we can find a δ2>0 so that if (n,c)⊤∈(KN−δ2,KN)×(0,δ2) then o((n−KN)2+c2)<ϵ2[(n−KN)2+c2]. Moreover, since c=0(cp) for c→0, there exists δ3>0 such that if c∈(0,δ3) then c<pcp. Now choose δ<min{δ1,δ2,δ3} so that θ1+ϵ1+2rNδ+ϵ2<0. Then, for any (n,c)⊤∈(KN−δ,KN)×(0,δ), 𝓛V3≤−2rN(n−KN)2−2rNc(n−KN)+o((n−KN)2+c2)+pcp(θ1+θ2) ≤−(2rN−ϵ2)(n−KN)2+2rNcδ+ϵ2c2+pcp(θ1+θ2) ≤−(2rN−ϵ2)(n−KN)2+pcp(θ1+ϵ1+2rNδ+ϵ2).
Choose θ3=min{2rN−ϵ2,−(θ1+ϵ1+2rNδ+ϵ2)}>0. It yields 𝓛V3≤−θ3V3, which implies that 𝓛V is a negative-definite function on (KN−δ,KN)×(0,δ). Therefore, by theorem 2.3 page 112 in [23], the trivial solution (KN,0)⊤ to the system (5) is locally stochastically asymptotically stable, i.e. the following statement For any ϵ>0, there exists a δ>0 so that for all (n,c)⊤∈(KN−δ,KN)×(0,δ), (6)ℙ{limt→∞V(t,v)=(KN,0)⊤}≥1−ϵ.
Second, we construct a compact set K˜⊂(0,KN)×(0,KC) so that the solution V(t,v) is recurrent relative to K˜, i.e. as the solution V(t,v) starts in (0,KN)×(0,KC) it will revisit K˜ infinitely many times in finite times almost surely. Indeed, consider function V4(n,c)=n+c, it can be easily seen that V4 is nonnegative and twice differentiable function on (0,KN)×(0,KC). We can compute that LV4=rNnln(KNn+c)+rCcln(KCn+c).
Choose K˜={(n,c)⊤∈(0,KN)×(0,KC)∣n+c≤KN}. Then, for any (n,c)⊤∈K˜c, we have KC<KN<n+c, which yields ln(KNn+c)<0 and ln(KCn+c)<0; hence, LV4<0 on K˜c. By theorem 3.9 page 89 in [25], it follows that V(t,v) is recurrent relative to K˜.
Third, we rewrite the system (5) in the Stratonovich form and we consider the corresponding control system (7)N˙ϕ=rNNϕln(KNNϕ+Cϕ), C˙ϕ=[ρC+b(Nϕ+Cϕ−KN)2a(Nϕ+Cϕ)2+(Nϕ+Cϕ−KN)2]Cϕln(KCNϕ+Cϕ) −12τ12Cϕln2(KCNϕ+Cϕ)+τ1Cϕln(KCNϕ+Cϕ)ϕ.
We let Vϕ(t,v)=(Nϕ(t,v),Cϕ(t,v))⊤ be the solution to the system (7) with control ϕ and initial value v=(n,c)⊤. We will show the following For given δ>0, there are a control ϕ and a time T>0 such that Vϕ(T,v)∈(KN−δ,KN)×(0,δ) for all v∈(0,KN)×(0,KC). To prove this statement, we need two following claims.
For any n0∈(0,KN) and c0, c1∈(0,KC) with c0>c1, there exist a control ϕ and a time T>0 such that |Nϕ(t,n0,c0)−n0|<ϵ and Cϕ(T,n0,c0)=c1 for all t∈(0,T).
Indeed, let f1(n,c)=rNnln(KNn+c) and f2(n,c)=[ρC+b(n+c−KN)2a(n+c)2+(n+c−KN)2]cln(KCn+c)−12τ12cln2(KCn+c).
We consider p1:=sup{|f1(n,c)|,|f2(n,c)|:|n−n0|≤ϵ,c1≤c≤c0}, and choose ϕ(t):=p2/lnKCN¯(t)+C¯(t) with p2 satisfying 0>c1−c0>ϵ(τ1p2c0p1+1) in which (N¯(t),C¯(t))⊤ solves the system N˙=f1(N,C), C˙=f2(N,C)+τ1cp2.
Note that τ1p2c0+p1<0, thus C˙ϕ(0,n0,c0)=f2(n0,c0)+τ1p2c0≤p1+τ1p2c0<0. Moreover, as Cϕ is continuous with respect to time, we can find a T0>0 such that C˙ϕ(t,n0,c0)<0 for all t∈(0,T0). It implies that Cϕ is decreasing from c0 over the period of time (0,T0). Now, assume that there exists a time t∈(0,ϵp1∧T0) that |Nϕ(t,n0,c0)−n0|>ϵ. Then, by Mean Value Theorem, for some η∈(0,t) we have ϵ<|Nϕ(t,n0,c0)−n0|=|Nϕ(t,n0,c0)−Nϕ(0,n0,c0)| =|N˙ϕ(η,n0,c0)|t=|f1(Nϕ(η,n0,c0),Cϕ(η,n0,c0))|t≤p1⋅ϵp1=ϵ, which yields a contradiction. Therefore, |Nϕ(t,n0,c0)−n0|≤ϵ for all t∈(0,ϵp1∧T0). Next assume that Cϕ(t,n0,c0)>c1 for all t∈(0,ϵp1∧T0). Then Cϕ(ϵp1∧T0,n0,c0)=limt→ϵp1∧T0Cϕ(t,n0,c0)=:c2>c1. However, by Mean Value Theorem, for some η∈(0,ϵp1∧T0) ϵ(τ1p2c0p1+1)<c1−c0≤Cϕ(ϵp1∧T0,n0,c0)−Cϕ(0,n0,c0) =C˙ϕ(η,n0,ϕ0)(ϵp1∧T0)≤(p1+τ1p2c0)ϵp1=ϵ(τ1p2c0p1+1), which is a contradiction. Thus, there exists a T∈(0,ϵp1∧T^0) so that Cϕ(T,n0,c0)=c1. This completes the proof of Claim 1.
For any n0, n1∈(0,KN) with n0<n1, there exist a c0∈(0,KC), a control ϕ, and a time T>0 such that Nϕ(T,n0,c0)=n1 and Cϕ(t,n0,c0)=c0 for all t∈(0,T).
Indeed, we consider the ODE system N˙=f1(N,C) and C˙=0 with initial value (n0, c0) where c0 will be determined later. Let (N˜(t),C˜(t)) be the solution to this system. It can be seen that C˜(t)=c0 for all t≥0 and N˜˙(t)=rNN˜(t)ln(KNN˜(t)+c0). With the feedback control ϕ satisfying [ρC+b(N˜(t)+C˜(t)−KN)2a(N˜(t)+C˜(t))2+(N˜(t)+C˜(t)−KN)2]C˜(t)−12τ12C˜(t)ln(KCN˜(t)+C˜(t))+τ1C˜(t)ϕ(t)≡0,
we have (Nϕ(t,n0,c0),Cϕ(t,n0,c0))≡(N˜(t),C˜(t)) for all t≥0. Now we choose c0>0 sufficiently small so that KN>n1+c0. It follows that infn0≤n≤n1f1(n,c)≥rNn0ln(KNn1+c0)>0.
Therefore, there exists a T>0 such that Nϕ(T,n0,c0)=n1 and Cϕ(t,n0,c0)=c0 for all t∈[0,T]. Hence Claim 2 is proved.
Finally, we will use the strong Markov property of the solution V(t,v) and the support theorem to prove that (KN,0)⊤ is a global asymptotically stable equilibrium to the system (5). Let Uδ:=(KN−δ,KN)×(0,δ). For any (n,c)⊤∈K˜, due to the statement and the support theorem of diffusion processes, there exist Tn,c>0 and pn,c>0 such that ℙ{(N(Tn,c,n,c),C(Tn,c,n,c))⊤∈Uδ}>2pn,c.
As the process (N(t),C(t))⊤ has the Feller property, we can find a neighborhood Vn,c of (n,c)⊤ such that, for all (n′,c′)⊤∈Vn,c, we get ℙ{(N(Tn,c,n′,c′),C(Tn,c,n′,c′))⊤∈Uδ}>pn,c.
On the other hand, since K˜ is a compact set, there exists a finite number of neighborhood Vni,ci(i=1,⋯,m) such that K˜=∪i=1mVni,ci Let T∗:=max{Tni,ci:i∈{1,…,m}}andp∗:=min{pni,ci:i∈{1,…,m}}.
For any (n,c)⊤∈(0,KN)×(0,KC), let τδn,c:=inf{t>0:(N(t,n,c),C(t,n,c))⊤∈Uδ}.
It is straightforward that, for any (n,c)⊤∈K˜, the event τδn,c<T∗ is followed from the event that there exists an i∈{1,…,m} such that (N(Tni,ci,n,c),C(Tni,ci,n,c))⊤∈Uδ, it follows that ℙ{τδn,c<T∗}≥ℙ{(N(Tni,ci,n,c),C(Tni,ci,n,c))⊤∈Uδ}>p∗>0.
Because the process (N(t,n,c),C(t,n,c))⊤ is recurrent relative to K˜, so we can define a sequence of finite stopping times ξ0=0, ξk=inf{t>ξk−1+T∗:(N(t,n,c),C(t,n,c))⊤∈K˜}(k≥1). Then we consider the event Ak={(N(t,n,c),C(t,n,c))⊤∉Uδforallt∈[ξk,ξk+T∗]}. It is clear that ℙ(Akc)=ℙ{(N(t,n,c),C(t,n,c))⊤∈Uδforsomet∈[ξk,ξk+T∗]}=ℙ(τδn¯,c¯<T∗)>p∗ where (n¯,c¯)⊤=(N(ξk,n,c),C(ξk,n,c))⊤∈K˜. Thus, ℙ(Ak)≤1−p∗. Next, by using the strong Markov property of (N(t,n,c),C(t,n,c))⊤, ℙ(n,c)⊤(A1∩A2)=ℙ(n,c)⊤(A1)⋅ℙ(N(ξ2,n,c),C(ξ2,n,c))⊤(A2)≤(1−p∗)2.
By induction, it yields ℙ(∩k=1nAk)≤(1−p∗)n→0 as n→∞, i.e. ℙ(∩k=1∞Ak)=0. It follows that ℙ(∪k=1∞Akc)=1. Thus, ℙ(τδn,c<∞)=1. Again, by strong Markov property of solution V(t,v), it follows from the fact τδn,c<∞ wp1 and (6) that for any (n,c)⊤∈(0,KN)×(0,KC) ℙ(limt→∞(N(t,n,c),C(t,n,c))⊤=(KN,0)⊤)≥1−ϵ.
Therefore this completes the proof of Theorem 2.2. □
In this section, we will conduct numerical simulations in details to demonstrate the usefulness of our stochastic model in predicting the likelihood of successful treatment of CAR T cell therapy in different protocols. We will firstly demonstrate the long-term behaviors of our model, and then conduct numerical simulations with three protocols. The parameter values and initial value are taken from the previous study [17] to simulate the dynamics of the stochastic system (3): KN=2.50×1011cells,KC=6.96×1010cells rN=1.70×10−1day−1, ρC=2.51×10−2day−1,a=4.23×10−1,b=5.25×10−1day−1, rB=(1∼50)×10−2day−1,γB=1.15day−1,kB=2.024×109cells, N(0)=3×109cells,C(0)=1.8×108cells,B(0)=9.486×1010cells.
Units of normal T cells (N), CAR T cells (C), and tumor cells (B) are measured as absolute lymphocyte counts and cell numbers in all figures below. For the time, it is counted as days right after CAR T Cell infusion. We fix noise intensities in all simulations as τ1=0.1, τ2=0.2, and τ3=0.1. We also fix the hyperparameter σ=0.05 for three parameters ρC, rB, and γB except in Figure 5. For each simulation, the values of ρC, rB, and γB are drawn from the normal distribution with standard deviation of σ of their mean values.
First of all, we demonstrate numerically the result in Theorem 2.2 that there is no positive equilibrium state for the stochastic system (3) in two cases. For the first case, the tumor net growth rate rB is taken to be smaller than 12τ22, and so tumor cells are completely eradicated with probability 1. Panels (a) and (b) of Figure 1 show that the endogenous T cell population gets its peak, the carrying capacity KN, around the day of 20, and the CAR T cell population achieves its carrying capacity KC in less than 10 days and starts decaying after that time. Tumor cells nearly get eradicated around the day of 100 due to the abundance of CAR T cells until day of 600. In this case, as time goes by, the normal T cell population approaches its maximum value and CAR T cells and tumor cells both decay to 0. The second case is that when the tumor net growth rate rB is taken to be much larger than 12τ22. Panels (c) and (d) in Figure 1 indicate that the dynamics of normal T cells and CAR T cells are pretty much the same as those in panels (a) and (b), while tumor cells decrease during the first 50 days and then progress right after that. This case shows that, as time goes by, the normal T cell population achieves its peak, CAR T cells eventually decay to 0, and tumor cells are likely proliferate and finally blow up almost surely. We can see from these two cases that, besides the tumor-eliminating effect of CAR T cells, noise also plays an important role in the dynamics of tumor cells and hence in the success of CAR T cell therapy. In the first case, when noise dominates the tumor net growth, it will drive tumor cells to decay in a long time period even though CAR T cells play limited role. For the second case, when the tumor net growth rate exceeds noise, it seems that noise does not help eradicate the tumor and tumor cells progress quickly as time goes by.
In reality, it is impractical to conclude that either noise is much larger than the tumor net growth rate so that noise suppresses tumor growth or noise is so small that the tumor net growth rate is greater than noise intensity and hence drives tumor cells to rapidly proliferate. It may happen that tumor net growth rate rB approximates more or less 12τ22. Especially, when rB is slightly larger that 12τ22, the ergodic invariant probability measures μ0 and μ2 do not exist. Then either the event that tumor cells are eradicated or the event that tumor cells progress to 120% of its initial size might occur with positive probability. Hence, the long-term behavior of the stochastic system (3) may not be appropriate to study the outcome of CAR T cell therapy. However, like the hybrid model in [17], to test the efficacy of different treatment strategies, we can use transient dynamics of the stochastic system (3) to investigate the behavior of two random times - the first time that tumor cells go below 1 cell (which is considered as tumor eradication), called the time to cure, and the first time that tumor cells progress to 120% of initial size, called the time to progression - when the baseline CAR T cell growth rate ρC, tumor net growth rate rB, and tumor killing rate γB are chosen from normal distribution with standard deviation of σ=5% of their mean values. We estimate the distribution of the time to cure and time to progress by simulating 10000 solution paths of the stochastic system (3) with the same initial value and then calculating probabilities from these solution paths in three different
Note that we take the median tumor net growth rate rB=10×10−2 for these three protocols. In the first protocol, each patient is treated by lymphodepletion before day 0 and one dose of CAR T cells at day 0. We run 10000 solution paths with the value of σ=0.05, which represent a cohort of 10000 simulated patients with different health conditions, of the stochastic system (3) with the same initial values. A patient is regarded to be cured if the tumor sample path of that patient firstly goes below 1 cell (Figure 2(c)) and a patient is considered to be progressed if the tumor sample path of that patient does not go below 1 cell but firstly exceeds 1.2 times of its initial tumor burden (Figure 2(d)). Figure 2(a) shows the percentage of the number of cured patients and the time it takes for patients to be cured ranges from day 30 to day 120. Figure 2(b) shows the percentage of the number of progressed patients and the time it takes for patients to be progressed ranges from day 180 to day 550.
For the second protocol, each patient undergoes lymphodepletion before day 0 and is treated by first dose of CAR T cells at day 0, and at day 15, receives the second dose of CAR T cells without lymphodepletion. We run 10000 solution paths of the stochastic system (3) with the same initial values and the value of σ=0.05. Note that each solution path, which is representative of a patient, includes a sample path of normal T cells, a sample path of CAR T cells, and a sample path of tumor cells. Figure 3(c) shows a sample path of normal T cells and a sample path of CAR T cells. For the sample path of the CAR T cell population, the amount of the second dose of CAR T cells, which is the same as that of the first dose, is added into the CAR T cells at day 15. We count the number of patients who are cured (meaning the sample path of tumor cells firstly goes below 1 cell) and the number of patients who progressed (meaning the sample path of tumor cells does not go below 1 cell but firstly exceeds 120% of its initial tumor burden) among 10000 simulated patients. Figures 3(a) and 3(b) show the percentage of the number of cured patients and progressed patients, respectively, in the second protocol.
In the third protocol, each patient is treated by lymphodepletion before day 0 and the first dose of CAR T cells at day 0, then treated by second lymphodepletion at day 10 and the second dose of CAR T cells with the same amount as the first dose at day 15. We run 10000 solution paths of the stochastic system (3) with the same initial values and the value of σ=0.05. Note that each solution path, which is representative of a patient, includes a sample path of normal T cells, a sample path of CAR T cells, and a sample path of tumor cells. Figure 4(c) shows a sample path of normal T cells and a sample path of CAR T cells in which CAR T cells are eradicated to 0 by the second lymphodepletion at day 10 and then the second dose of CAR T cells is infused at day 15, while normal T cells are dropped dramatically by second lymphodepletion at day 10 and recovered up to its initial value at day 15. Figures 4(a) and 4(b) show the percentage of the number of cured patients and progressed patients, respectively, in the third protocol.
The probability of cure and probability of progress in Figures 2, 3, and 4 are computed by adding all the percentages of the number of cured patients and progressed patients, respectively.
In the first protocol, most patients are cured between days 30 and 80 (Figure 2a). It is rare for us to find late cure events up to or after day 135. Meanwhile, the time to progression, defined as 120% of the initial tumor burden, occurs anywhere between days 180 and 525 (Figure 2b). The reason for this big difference in time between two random times is due to the time period when CAR T cells decline is longer than the time period when CAR T cells expand. For the second protocol, the distributions of both time to cure and time to progress look like the same as those in the first protocol. However, the efficacy of the second dose of CAR T cell infusion can be raised up for later time points from day 80 to day 90 (Figure 3a). This indicates that, without lymphodepletion, the second dose of CAR T cell infusion could only have a fairly measurable effect in patients which exhibit some response but not cure. In the third protocol, we can see that it has a huge impact on the patient survival. In Figure 4a, the distribution of the time to cure concentrates densely from day 35 to day 70 and then spreads out up to day 135. It can be explained that the second lymphodepletion at day 10, after the first dose of CAR T cell infusion, kills most of normal T cells and CAR T cells; when the second dose of CAR T cells is treated at day 15, it takes a while for normal T cells to out-compete CAR T cells, so the CAR T cell population has more time to get its peak and maintain its large amount for a while to interact and kill tumor cells before normal T cells achieve its carrying capacity.
It is obvious that, from our simulations, both parameter σ and noise intensities τi(i=1,2,3) account for the patient heterogeneity and patient variabilities in treatment outcomes, respectively. By our theoretical analysis in Section 3 and numerical simulations in Figure 1, noise plays a crucial role in determining the success of CAR T cell treatment. To test whether variability in timing to cure or progression is also shaped by the hyperparameter σ and noise, we perform 100 simulations for the stochastic system (3) using the first protocol with different values of σ and τ3 both ranging from 0.001 to 0.15 and then estimate the probability of cure and the distribution of time to cure for each value of σ and τ3 separately. Figure 5(a) shows the probability of cure ranges between 25% and 60% but does not show any pattern. Similarly, Figure 5(c) indicates that the probability of cure ranges from 40% to 62% but still does not show any pattern. In Figure 5(b) and Figure 5(d), as increasing either σ or τ3, the possible range of times to cure seems unpredictable, which contrasts with results in [17]. This demonstrates that both hyperparameter σ and noises play a key role in predicting the success of CAR T cell treatment.
Finally, to confirm necessity of lymphodepletion for the treatment with the second dose of CAR T cell infusion, we analyze numerically the change in probability of cure as a function of the time of the second dose of CAR T cell infusion with lymphodepletion 5 days before and without lymphodepletion. Here we chose all parameter values as we introduced at the beginning of this section except for a low tumor killing rate γB=0.9 since it would not almost surely lead to cure if only one dose of CAR T cells was infused. At each time of the second dose of CAR T cell infusion, 100 stochastic simulations (corresponding to a cohort of 100 simulated patients) of the stochastic system (3) are run, and then its probability of cure is computed by adding all the percentages of the number of cured patients whose tumor sample paths firstly go below 1 cell among 100 patients simulated under two treatment protocols. In Figure 6a, it can be seen that the probability of cure nearly approaches 1 with increasing time of the second dose of CAR T cell infusion since the first infusion. This is because of the fact that CAR T cells have an early selective advantage against a lowered amount of normal T cells by lymphodepletion, which leads to a second peak of CAR T cells (Figure 4c). On the other hand, a second infusion of CAR T cells without lymphodepletion is able to be effective in elevating its probability of cure from 0 to about 10–30% (Figure 6b). Thus, to improve success with a second infusion of CAR T cells, it is necessary to reset normal T cells by lymphodepletion. This indicates that a second, lower dose of lymphodepletion alone may be good to get maximal benefits for patients if it does not wipe out all CAR T cells but lowers overall T cell density. This observation suggests that lymphodepletion is a key factor to increase the chance of successful outcomes of CAR T cell therapy.
Our stochastic model is based on the hybrid model [17]. The hybrid model is a mixture of deterministic ODE model with birth-death stochastic process. In general, a deterministic model is based on understanding of mechanisms of relevant components in a system. The hybrid model builds a competition for the total immune cell density between endogenous T cells and CAR T cells and a saturation killing functional relation by CAR T cells into it. It considers tumor growth follows the linear growth pattern. For constructing stochastic models in terms of stochastic differential equations, we consider both external and internal noises in the system. In this study, we consider the noises arising with CAR T cell expansion rate after infusion, tumor growth rate, and the rate of CAR T cells killing tumor cells. Those noises may represent variability of patients. Our stochastic model has three invariant ergodic measures which corresponds to the three equilibrium solutions of the deterministic model. When the tumor growth rate is smaller than the quantity R2, the two invariant ergodic measures corresponding to equilibrium solutions (0,0,0) and (0,KC,0) are attractors. When the tumor growth rate is even small, say smaller than the quantity R1, the invariant ergodic measure corresponding to the equilibrium solution (KN,0,0) is a global attractor. These invariant ergodic measures may represent some cure states from the long-term viewpoint. The stochastic model overcomes a drawback of the hybrid model which lacks stable cure states. However, the tumor growth rate alone plays an important role in determining long-term outcomes of the treatment. The killing rate by CAR T cells and capacity of CAR T cells only indirectly impact tumor growth by the quantity R2. This provides some rooms for further improvement of model construction.
Our numerical simulations in Section 4 have both similarities and differences compared to those of the hybrid model in [17]. Both models yielded similar distributions of time to cure and time to progress in the first protocol (see Figure 2 and Figure 4 in [17]). Both models demonstrated that lymphodepletion is an important component of the success of CAR T cell therapy when the second dose is administered (see Figure 6 and Supplementary Figure 5 in [17]). The difference in numerical results between the two models lies in the patient heterogeneity and the patient variability in treatment outcomes, which are shaped by the hyperparameter σ and how stochasticity is incorporated into the two models (see Figure 5 and Figure 3 in [17]). The main reason for this difference relies on how noise is simulated in each model. The hybrid model in [17] updates its solution at the next time point by using a Poisson process with a random time step, which is an exponential distribution with a rate parameter related to the solution at the previous time point. This time step represents the minimum amount of time for the next interaction among cell populations to occur. So there is no controlled noise intensity for each parameter in the hybrid model in [17]. Meanwhile, we used white noise and three different noise intensities to capture the variability in three important parameters, which are picked up from a normal distribution with variance changing over time. This explains why the probability of cure and the distribution of time to cure did not follow any pattern when we changed the hyperparameter σ and the noise intensity τ3 in Figure 5.
As mentioned in Introduction part, two ways of interpreting the solutions of a SDE model can be used to evaluate a cancer therapy efficacy and variability in treatment outcomes. One can consider SDE solutions as either treatment outcomes of different simulated patients or different treatment possibilities of a single simulated patient. Utilizing which interpretation is relevant depends on the specific details of the model and how it is set up. In this study, we aimed to use the first interpretation for several reasons. First, our stochastic model was based on the hybrid model in [17], which also utilized the first interpretation, and so we can compare our numerical results with theirs. Second, running simulations for our model with the first interpretation is less expensive than that with the second interpretation. Furthermore, in our model, we incorporated white noises with large values of noise intensities and a hyperparameter into three important parameters that can account for tumor heterogeneity across a cohort of patients.
However, the second interpretation, where different solution paths can be explained as different possible treatment outcomes for one single patient, can be used for our model as follows. When a patient is injected with a CAR dose, treatment outcome of this patient may vary due to the baseline growth rate of CAR T cells, ρC, and the rate at which tumor is killed by CAR T cells, γB. These two parameters can represent CAR dose effectiveness and may be subject to uncertainty. So we can incorporate white noise with small values of noise intensities into these two parameters to capture the variability in treatment outcomes of one patient and use hyperparameter to catch up different tumor growth rates among a cohort of patients. This new SDE model with the second interpretation is worthy of future investigation.
In this study, we aim to use the transient dynamics of a stochastic differential equation system to numerically simulate different treatment protocols in practice because it offers some promising advantages. First, stochastic differential equations capture the randomness and variability inherent in biological systems, making them suitable for modeling complex and dynamic processes. Second, transient dynamics enable the incorporation of patient-specific variability, which is crucial in numerically simulating diverse patient populations. Third, transient dynamics help in identifying critical time points during the course of treatment, where certain events or changes in the system are more likely to occur. This information can guide the timing of interventions and the monitoring of treatment efficacy. In addition to these advantages, we also address some drawbacks of this framework. The first issue is the identifiability of stochastic differential equations, relying on robust biological and clinical data. Obtaining precise estimates for the parameters of a stochastic model can be challenging, and uncertainties in parameter values may affect the reliability of predictions. The second challenge is validating the predictions of a transient stochastic model against real-world clinical data, which can be a daunting task. The inherent stochasticity and variability may make it difficult to discern whether observed deviations from clinical data are due to model limitations or genuine biological variability.
As demonstrated in this study, stochastic differential equations offer a powerful framework for capturing the inherent variability and uncertainties in clinical data, providing insights and understanding into the mechanistic function of different new treatment protocols. However, as clinical data continue to expand in the near future, this framework would exhibit some of its limitations, such as accurate parameterization and validating predictions against real-world clinical data, as mentioned above. To overcome these limitations, using machine learning techniques, which can learn patterns from data and make accurate predictions, can be a remedy. In a clinical context, machine learning alone is capable of precisely predicting individual patient outcomes from past observations using a patient database. Nevertheless, it can potentially identify which of the existing treatments is most adequate but is intrinsically unable to recommend new treatment protocols or provide accurate predictions for new treatments. Stochastic differential equations can help overcome the limitations of the machine learning approach. Thus, the integration of stochastic differential equations and machine learning presents a promising approach for studying the efficacy of CAR T cell therapy as clinical data continue to increase in the future. We will continue to study this approach in the future.