Basic stochastic model for tumor virotherapy

The complexity of oncolytic virotherapy arises from many factors. In this study, we incorporate environmental noise and stochastic effects to our basic deterministic model and propose a stochastic model for viral therapy in terms of Ito stochastic differential equations. We conduct a detailed analysis of the model using boundary methods. We find two combined parameters, one describes possibilities of eradicating tumors and one is an increasing function of the viral burst size, which serve as thresholds to classify asymptotical dynamics of the model solution paths. We show there are three ergodic invariant probability measures which correspond to equilibrium states of the deterministic model, and extra possibility to eradicate tumor due to strong variance of tumor growth rate and medium viral burst size. Numerical analysis demonstrates several typical solution paths with biological explanations. In addition, we provide some medical interpretations and implications.


Introduction
Cancer is a genetic disease. It is caused by changes to genes, which control how cells grow and divide. A DNA change can cause genes involved in normal cells to become oncogenes. A oncogene is difficult to be turned off and so it causes cells grow without limits. When too many cells are accumulated, they form a solid tumor, which is masses of tissue. Cancer therapy is a broad area of research, which may have three subfields: immunotherapy, gene therapy, and oncolytic virotherapy. Immunotherapy relies on the concept of stimulating the body's immune system to recognize and destroy cancer cells. Cancer cells harvested from patients are grown in vitro. Then these cells are engineered to be more recognizable to immune system by some substances or genes. These altered cells are grown in vitro and killed and their contents are incorporated into a vaccine that will be administered to patients, in order to boost the patients' immune responses. But this method had limited success [1]. Gene therapy also is called gene transfer treatment which refers to the insertion of a foreign gene into the cancer cell or surrounding tissue to express specific genes such as suicide genes. This method does not rely on the immune system. Typically, replication-incompetent viruses, such as modified strains of adenovirus, have been used to deliver these genes. However, this technique met a lot of difficulties such as gene silence, the gene not being expressed for long enough time period [2].
In this paper, we will focus on mathematical analysis of the third type of cancer treatments, known as oncolytic virotherapy. Oncolytic viral therapy is considered to be a promising therapeutic strategy to treat solid tumors [3] and has shown its efficacy in clinical trials [4,5]. This treatment involves the use of oncolytic viruses, namely genetically modified viruses to selectively infect cancer cells and induce cell death through lysis and further propagation of the virus. A number of viruses such as adenovirus, ONYX-15 and CV706, herpes simplex virus 1, and wild-type Newcastle disease virus have been used for such purposes. These viruses are shown to be unharmful to normal cells and tumor specific. In contrast to gene transfer treatments which utilize replication-incompetent viruses to alter the characteristics of cancer cells, oncolytic viruses have the ability to selectively replicate within the target cancel cell, resulting in the amplification effect in areas of tumor growth, allowing for safer doses of viral agent to be used in treatment [6].
Mathematical models formulated in terms of ordinary differential equations (ODEs) have been applied to understand spreading dynamics of oncolytic viruses through tumors for nearly twenty years. The early ODE model was proposed by Wodarz [7,8], and was generalized by Dingli et al [9] later on. These models were formulated with three physical variables: uninfected tumor cells, a free virus population, and tumor cells infected by virus particles. The uninfected tumor cells were assumed to undergo logistic growth, and infected by virus particles, which multiply rapidly with infected tumor cells. Infected tumor cells were removed from the system due to natural or virus-inflicted death, resulting in new virus particles bursting to the free virus population. Motivated by experimental evidence, Bajzer et al [10] suggested that the forming of syncytia by fusing of uninfected and infected tumor cells rather than the free virus particles was the physical mechanism which drives intratumor virus spreading. Komarova and Wodarz [11] proposed and analyzed several general mathematical formulations for oncolytic virus infection in terms of systems of two ordinary differential equations, which categorized two types of virus spread, slow and fast spread. Our work [12] proposed a simple system of three ordinary differential equations to describe the interactions among uninfected tumor cells, infected tumor cells, and oncolytic viruses. Our analytic and numerical results concluded that the oncolytic viral dynamics is mainly determined by the viral burst size. To further understand the complexity of immune responses in virotherapy, we incorporated the innate immune response into our basic model for virotherapy and investigated how the innate immunity affects the outcome of virotherapy [13].
Stochastic effects are encountered in many biological and medical systems. Stochastic models may be able to capture some stochastic effects or variations in dynamics of biological and medical problems. In recent years, several attempts have been made to characterize viral dynamics for oncolytic virotherapy using stochastic differential equations (SDEs) such as Yuan and Allen [14], Kim et al [15], and Rajalakshmi et al [16,17]. Most of these stochastic models were formulated by transforming ODE systems using the method proposed in [18]. These transformed SDE models may have some computational advantages. In this study, we propose a system of stochastic differential equations for tumor virotherapy and carry out its analysis and computation based on some suggestions from research presented in articles [19,20,21].
In [12], we proposed a common basic deterministic model for oncolytic virotherapy that includes the virus burst size b explicitly as follows, (1.1) in which x stands for the uninfected tumor cell population, y the infected tumor cell population, and v the free virus population. The tumor growth is modeled by a logistic pattern with the growth rate ρ and carrying capacity of the tumor size C. The coefficient β represents the infectivity of the virus. The infected tumor cells die with a rate δy, which means the average life time of infected tumor cells is 1 δ . The viral burst size b is the number of new viruses released from a lysis of an infected tumor cell. The term γv is the clearance rate of free virus particles by various reasons including non-specific binding and generation of defective interfering particles.
There are several ways to incorporating environmental noise or stochastic effects into mathematical models. Suppose P is a population, its growth or change is modeled dP dt = f(t, P ) in the deterministic situation. To count for environmental noise and stochastic effect, we may consider that each individual in the population make almost same contribution to the stochastic effects and receive the same environmental noise. Then, we may model the environmental noise and stochastic effects of the population is proportional to the population P. In other words, the environmental noise and stochastic effects can be represented by τPξ, where ξ is the unit noise [22] and τ can be thought as a way to measure an average variation of each individual. In general, we take the noise to be white noise assumed that W 1 , W 2 , and W 3 are mutually independent Wiener processes. Such, we obtain a system of three Ito stochastic differential equations which is a basic stochastic model for oncolytic viral therapy as follows. dx = ρx 1 − x + y C − βxv dt + τ 1 x 1 − x + y C dW 1 , dy = (βxv − δy)dt + τ 2 ydW 2 , dv = (bδy − βxv − γv)dt + τ 3 vdW 3 . (1. 2) The analysis of the deterministic system (1.1) (see [12]) shows that the virus burst size b plays a crucial role in determining its dynamics. We found two important thresholds of the burst size that give a complete picture of dynamical behavior of (1.1). Our aim in this work is to analyze the SDE system (1.2) in order to find thresholds under which we can identify the extinction or persistence of the tumor cells and, furthermore, figure out how noise intensities affect the dynamics of the SDE system (1.2).
The rest of this paper is organized as follows. In Section 2, we simplify the stochastic system (1.2), introduce notations, present our results, state medical interpretations. In Section 3, we analyze our model using boundary analysis technique, and prove our results. In Section 4, by means of published data, we demonstrate typical dynamic behaviors of our stochastic model by numerical simulations and explain possible biological meanings. We also provide a brief discussion, some open problems, and possible future work. Finally, we present some basic properties of Generalized Inverse Gaussian distribution in Appendix.

Results and interpretations
First of all, for simplicity, we non-dimensionalize the system (   All parameters are positive. Assume that we are working on a complete probability space Ω, ℱ, ℱ t t ≥ 0 , ℙ with a filtration ℱ t t ≥ 0 satisfying the usual condition. The process given by the solution to the system (2.2) will be denoted by u or u(t) = (x(t), y(t), v(t)), t ≥ 0.
We denote the drift term and the diffusion term of the system (2. where F u is the gradient of F and F uu is the Hessian matrix of F. We use ℙ u or ℙ x, y, v to denote the probability law on Ω when the solution path starts at u = (x, y, v) and E u or E x, y, v is the expectation corresponding to ℙ u .
Based on the results (see Theorem 1.1 and Theorem 1.3 in [23]) about asymptotic behaviors of stochastic Kolmogorov systems in non-compact domains, we derive a sufficient and almost necessary condition to determine the extinction and persistence of populations of uninfected tumor cells, infected tumor cells, and free viruses. However, these results cannot be applied directly to our model because the drift term of the system (2.2) is not in Kolmogorov form as that of the system in [23]. To apply these results, we need to change variables. In view of (2.2), set the transformation of variables x = x, y = y, v = yz, or z = v y , and use Ito's formula, we get  The following theorem, that will be proved in Section 3, guarantees the global nonnegativity of the solution of the system (2.3) for any positive initial value. For any initial value (x(0), y(0), z(0)) ∈ ℝ + 3 : = (x, y, z): x ≥ 0, y ≥ 0, z ≥ 0 , there exists a unique a.s. continuous global solution (x(t), y(t), z(t)), t ≥ 0, that remains in ℝ + 3 a.s. In particular, if x(0) = 0 then x(t) = 0 for all t > 0 a.s. and if x(0) > 0 then x(t) > 0 for all t > 0 a.s. Similarly, if y(0) = 0 then y(t) = 0 for all t > 0 a.s. and if y(0) > 0 then y(t) > 0 for all t > 0 a.s. Finally, if z(0) ≥ 0 then z(t) > 0 for all t > 0 a.s. Furthermore, the solution (x(t), y(t), z(t)) is a strong Markov process that possesses the Feller property.
K θ (w) , and K θ (·) is the modified Bessel function of the third kind with index θ which is given by With these parameters and their thresholds, we give the complete picture of the stochastic dynamics of the system (2.3). Our main result is stated in the following theorem that will be proved in Section 3.
• If λ > 0 then the solution u(t) is strongly stochastically persistent in the sense that the solution converges to its unique invariant probability measure μ 3 supported by ℝ + 3,°.
From the transformation of variables, the information about the system (2.2) can be obtained. We write them as the interpretation of our main theorem. We will give some medical interpretation of these results and compare with our study in [12].

Interpretation 2.1.
Consider the non-dimensionalized uninfected tumor cell population x(t), infected tumor cell population y(t), and free virus population v(t) start in ℝ + 3,°: = (x, y, z): x > 0, y > 0, v > 0 such that x + y ≤ 1, which corresponds to the system (2.2). Then, according to the thresholds ζ and λ, we can describe how each population will evolve as follows. Case 1.-When ζ < 0, the tumor cannot be eradicated completely a.s.
Using the transformation of variables or directly deduce, we have a similar proposition as 2.1. probability measures μ 1 = δ 0 * × δ 0 * × δ 0 * , μ 2 = δ 1 * × δ 0 * × δ 0 * , and μ 3 , which correspond to Q 1 , In our study [12], we obtained asymptotic properties of the system (1.1). There are three equilibrium solutions Q 1 , Q 2 , and Q 3 . 2. However, we would like to understand these results from the original system or how environmental noise and stochastic effects change the dynamical behaviors of the original system (1.1). Then, we need to understand how these two combined parameters connect to original parameters and their biological meanings. We have a proposition about the parameter λ.

Proposition 2.3.
The parameter λ = ab R θ (w) − 1 − τ 2 2 2 is an increasing function of the virus burst size b. Also consider λ is a function of noise intensities τ 2 and τ 3 , and set λ: = lim τ 2 , τ 3 (0, 0) λ. Then λ = 0 if and only if b = b s 1 = 1 + c a , λ < 0 if and only if b < b s 1 , and λ > 0 if and only if b > b s 1 ; or simply, λ also is a increasing function of b.
The parameter ζ combines infected tumor cell lysis rate δ, virus degradation rate γ, and their stochastic variation τ 2 and τ 3 , which describes possibilities if the tumor can be eradicated. More specifically, ζ = 2c − 2 − τ 2 will disappear, and only tumor cell population is left a.s., or the system approaches the invariant probability measure μ 2 = δ 1 * × δ 0 * × δ 0 * . If b is greater than the threshold b s 1 which corresponds to λ > 0, three populations will coexist, or the system approaches the invariant probability measure μ 3 , in which we may say viral therapy achieve some partial success.
When ζ > 0 which means which corresponds to λ < 0, the noise intensity τ 1 or tumor cell variance τ 1 2 is smaller than the double of the scaled tumor cell growth rate, then the viral treatment will complete fail. When b is greater than the threshold b s 1 which corresponds to λ > 0, and the noise intensity tumor cell growth rate, the viral therapy will eradicate the tumor. A medical implication could be that viral therapy can success without too big virus burst size.

Analysis of the model
This section is devoted to proving results in Section 2.
This completes the proof. □

Phan and Tian Page 11
Math Biosci Eng. Author manuscript; available in PMC 2022 February 25.

Proof of Theorem 2.2
Before giving the detailed proof of the main theorem 2.2, we analyze solutions of the system    where C 3 is some positive constant. The integrand can be written as Phan  Clearly, there are k 1 and k 2 in Z + such that Hence s(0+) = −∞ and s(∞) = ∞. So z(t) oscillates between 0 and ∞, and thus (3.3) has a unique invariant measure π 2 ~ GIG(θ, χ, ψ), which is the generalized inverse Gaussian distribution with parameters θ ∈ ℝ, χ > 0, and ψ > 0 (see the Appendix), whose density Then, we have the following claim.
In order to prove Claim 3.1, we will utilize the non-negative semi-martingale convergence theorem (see Theorem 3.1 in [24]) and the following two lemmas, whose proofs will be given in the end of this subsection.  Then But this contradicts the fact that lim Note that, since 0 ≤ x(t) ≤ 1 for all t ≥ 0 a.s. by Lemma 3.1, we can show that E u (x(t) + y(t)) ≤ 1. This implies that E u y(t) ≤ 1 and hence, by Markov's inequality, we can prove that y(t) is bounded a.s. Thus 1 + c T u 2 for some constant K 4 > 0. When z is large enough, we can choose γ b > 0 sufficiently small This shows that where ‖u‖ := |x| + |y| + |z|. Moreover, it is easy to compute diag g 1 , g 2 , g 3 ΓΓ T diag g 1 , g 2 , g 3 = which is positive definite for all (x, y, z) ∈ ℝ + 3,° satisfying x + y ≤ 1. Thus, Assumption 1.1 in [23] is fulfilled for the system (2.3).

Author Manuscript
Author Manuscript
We complete the proof. □ Proof of Lemma 3.2.: By Ito's formula, since 0 ≤ x(t) ≤ 1 for all t ≥ 0 a.s. (by Lemma 3.1), Taking expectation both sides yields here we have used the equality E u x 2 (t)z 2 (t) ≥ E u x(t)z(t) 2 . This implies that Therefore for all t ≥ 0 we obtain So, Lemma 3.2 is proved.  virus particles over relative infected tumor cells. In our simulation, time is scaled or relative time since T = δt. In [27] and [12], the parameter b, the burst size of free virus particles, plays a pivotal role in determining the success of glioma virotherapy. So we will simulate the trajectories of the system (2.2) with the initial value (0.5, 0.5, 1.5) and the modified system (2.3) with the initial value (0.5, 0.5, 3) as b and noise intensities are varied while all the other parameters are fixed.
Example 1.-We illustrate the situation when λ < 0. In the figures 1 and 2, we take b = 5, τ 1 = 0.2, τ 2 = 0.3, and τ 3 = 0.2. By computation, θ = 7.3077, w = 22.8191, R θ (w) = 1.39, and hence λ = −0.0141 < 0. Figure 1 indicates that the relative uninfected tumor cells increase to the relative maximal cell density of the tumor, which is 1; the relative infected tumor cells decay to zero; and the ratio of relative free virus particle over relative infected tumor cells reaches an equilibrium state, which explains why the relative free viruses are wiped out as in Figure 2. These two pictures verify the conclusion in Case 1 of Theorem 2.2. In terms of biological meaning, since the burst size b is not big enough, the number of new viruses released from a lysis of an infected cell is inconsiderable when compared with the number of free viruses dying out. Because of that, in the early stage, the population of free virus particles increases by contribution from lysis of some infected tumor cells but later the number of free viruses decrease and decay to zero. The decrease in free viruses leads to decrease in infected tumor cells and hence the infected also decay to zero. Then the uninfected becomes less and less infected by free viruses, and finally increases to its carrying capacity. Therefore, virotherapy fails. phenomenon can be explained as follows. When the burst size of free viruses is increased, say, to 10, the number of new viruses from a lysis of an infected tumor cell becomes significant. Then, the dying-out infected cells contribute much to the number of free viruses within the tumor. The number of free virus particles is big enough to prevent the growth of uninfected tumor cells. Some of the uninfected getting infected by free viruses becomes the infected, while some of them keep growing. Three populations interact in the mutually coexistent way. This shows that injecting free viruses with stronger burst size into the tumor yields better treatment.
If the burst size b is doubled to 20 while the noise intensities are the same as in Example 1 and 2, Figure 5 indicates that all solution paths still persist and oscillate most of the time around an equilibrium state. The difference is that the tumor load, which is the total number of the uninfected and the infected cells, is much smaller than the tumor load with the burst size b = 10. When we increase the burst size b to 40 and noise intensities are kept the same, the solution behaves differently. Figure 6 shows that all solution paths represent a pulsating Phan  oscillatory. The minimum of the uninfected tumor population can reach a very small value comparing with the maximum tumor size. In this case, the tumor may be regarded to be undetectable and then we consider that the tumor is eradicated. This phenomenon becomes more visible when the burst size is taken very large, say b = 80, as illustrated in Figure 7.
Thus, the viral treatment can be seen to be some success.

Discussion
In this paper, our basic virotherapy model of stochastic type is able to predict the dynamics of viral therapy based on the viral burst size b and noise intensities. We found thresholds ζ and λ that provide conditions for various outcomes of our stochastic model. The parameter ζ combines infected tumor cell lysis rate δ, virus degradation rate γ, and their stochastic variation τ 2 and τ 3 , which describes possibilities if the tumor can be eradicated by viral therapy. The parameter λ is a differential function of the viral burst size b and the noise intensities τ 2 and τ 3 , and is increasing as b increases. We elaborate some medical implications of these parameters and medical outcomes theoretically in Section 2.
We also numerically demonstrate these dynamical outcomes and present more biological explanations in Subsection 4.1. We compare our stochastic model and its deterministic counterpart. Equilibrium states of deterministic models correspond to ergodic invariant probability measures of stochastic models. However, our stochastic model demonstrate some new features. For the deterministic model, there is no possibility to eradicate the tumor, but for the stochastic model, there is a case where the tumor can be eradicated. This is due to introducing a big variance τ 1 2 of tumor cell growth rate.
There are several interesting questions arisen in our study. For two ergodic invariant probability measures on the boundary, we obtain their explicit probability distributions, so that we can compute their expectations which correspond equilibrium solutions of the deterministic system. For the ergodic invariant probability measure supported by the interior of the domain, we are unable to find its probability distribution explicitly in this study although we know it correspond the coexistence equilibrium solution of the deterministic model. One question is to find the explicit expression of this probability distribution. A second question is about Hopf bifurcations. In the deterministic model, when the viral burst size passes through the second threshold value b s 2 , there is a Hopf bifurcation with appearance of three families of periodic solutions. We know that λ is an increasing function of the viral burst size b, and its threshold b s 1 also serve well for classification of solutions to stochastic model. We ask if there is a Hopf bifurcation for the stochastic model when the viral burst size b passes through b s 2 . This would be very interesting from both theoretical and practical viewpoint.
One of the major challenges in current medical practice of oncolytic viral therapy is to get insight into the complexity of the immune responses. Understanding the dynamics of oncolytic virotherapy in the presence of immune responses is a considerable need. The innate immune response has a tendency to reduce the efficacy of oncolytic viral treatment by lowering new virus multiplication and blocking the infection spreading, while the stimulated adaptive immune response tends to reduce tumor cells. So the extension of our stochastic model to incorporate the innate and adaptive immune systems is expected. We plan to conduct these studies in the future.