Basic stochastic model for tumor virotherapy

1.   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,

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 $\rho$ and carrying capacity of the tumor size $C$. The coefficient $\beta$ represents the infectivity of the virus. The infected tumor cells die with a rate $\delta y$, which means the average life time of infected tumor cells is $\frac{1}{\delta}$. The viral burst size $b$ is the number of new viruses released from a lysis of an infected tumor cell. The term $\gamma 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 $\frac{dP}{dt} = f(t, P)$ in the deterministic situation. To count for environmental noise and stochastic effects, 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 $\tau P\xi$, where $\xi$ is the unit noise ^[22] and $\tau$ can be thought as a way to measure an average variation of each individual. In general, we take the noise to be white noise $\xi = \frac{dW}{dt}$, where $W = W(t)$ is the standard Wiener Process. So, we obtain a Ito stochastic differential equation $\frac{dP}{dt} = f(t, P)+\tau P\frac{dW}{dt}$, or $dP = f(t, P)dt+\tau PdW$ as a stochastic model for the population $P$. We may call the noise added this way the linear noise. For our model (1.1), we will incorporate linear noise to the infected tumor cell population and free virus population. that is, we will add $\tau_2y\frac{dW_2}{dt}$ and $\tau_3v\frac{dW_3}{dt}$ to the second and the third equations of (1.1), respectively. However, for the uninfected tumor cell population, we incorporate environmental effects into per capital growth rate $\rho$. That is, we replace $\rho$ by $\rho +\tau_{1}\frac{dW_{1}}{dt}$, where $\tau_{1}$ represent the strength of the noise contributed by each tumor cells. It should be 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.

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.

2.   Results and interpretations

First of all, for simplicity, we non-dimensionalize the system (1.2) by setting $T = \delta t$, $x = C{\bar x}$, $y = C{\bar y}$, $v = C{\bar v}$, $r = \frac{\rho}{\delta}$, $a = \frac{\beta C}{\delta}$, $c = \frac{\gamma}{\delta}$, $\tau_1 = {\bar \tau}_1$, $\tau_2 = {\bar \tau}_2$, and $\tau_3 = {\bar \tau}_3$. Then (1.2) becomes

Dropping all bars over the parameters and variables and writing $T$ as $t$, we obtain

All parameters are positive. Assume that we are working on a complete probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\ge 0}, \mathbb{P})$ with a filtration $\{\mathcal{F}_t\}_{t\ge 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\ge 0$. We denote the drift term and the diffusion term of the system (2.2) by

Let $L$ be the infinitesimal generator of the process $u$ and, for any smooth enough functions $F: \mathbb{R}^3_{+}: = [0, \infty)^3\to \mathbb{R}$, the generator $L$ acts as

where $F_u$ is the gradient of $F$ and $F_{uu}$ is the Hessian matrix of $F$. We use $\mathbb{P}_u$ or $\mathbb{P}_{x, y, v}$ to denote the probability law on $\Omega$ when the solution path starts at $u = (x, y, v)$ and $\mathbb{E}_u$ or $\mathbb{E}_{x, y, v}$ is the expectation corresponding to $\mathbb{P}_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 = \frac{v}{y}$, and use Ito's formula, we get

Then (2.2) is changed to

We still denote by $u(t) = (x(t), y(t), z(t))$ the solution process of the system (2.3). The drift term and the diffusion term of (2.3) are also denoted by

The following theorem, that will be proved in Section 3, guarantees the global non-negativity of the solution of the system (2.3) for any positive initial value.

Theorem 2.1. For any initial value $(x(0), y(0), z(0))\in \mathbb{R}_{+}^3: = \{(x, y, z):x\ge 0, y\ge 0, z\ge 0\}$, there exists a unique a.s. continuous global solution $(x(t), y(t), z(t))$, $t\ge 0$, that remains in $\mathbb{R}_{+}^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)\ge 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.

Our analysis in Section 3 shows that there are only two ergodic invariant measures

of (2.3) on the boundary $\partial \mathbb{R}_{+}^3$. Here $\delta_0^{*}$ and $\delta_1^{*}$ are Dirac measures with mass at $0$ and $1$, respectively. The invariant measure $\pi_1$ has the inverse gamma distribution: $\pi_1\sim\text{IG}\left(\frac{2(c-1-\tau_2^2)}{\tau_2^2+\tau_3^2}+1, \frac{2b}{\tau_2^2+\tau_3^2}\right)$. The invariant measure $\pi_2$ has the generalized inverse Gaussian distribution: $\pi_2\sim\text{GIG}(\theta, \chi, \psi)$, where $\theta = \dfrac{2(1+\tau_2^2-a-c)}{\tau_2^2+\tau_3^2}-1$, $\psi = \dfrac{4a}{\tau_2^2+\tau_3^2}$, and $\chi = \dfrac{4b}{\tau_2^2+\tau_3^2}$.

To classify solutions of the system (2.3), we define two combined parameters as follows.

where $w = \dfrac{4\sqrt{ab}}{\tau_2^2+\tau_3^2}$, $R_{\theta}(w) = \frac{K_{\theta +1}(w)}{K_{\theta}(w)}$, and $K_{\theta}(\cdot)$ is the modified Bessel function of the third kind with index $\theta$ 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.

Theorem 2.2. Assume that the initial values $u = (x, y, z)$ are in $\mathbb{R}_{+}^{3, \circ}: = \{(x, y, z):x > 0, y > 0, z > 0\}$ such that $x+y\le 1$. The complete classification of solutions of the system (2.3) is as follows.

Case 1. When $\zeta < 0$, there is only one ergodic invariant measure $\mu_2$ for solutions of (2.3) on the boundary $\partial \mathbb{R}_{+}^3$.

● If $\lambda < 0$ then $x(t)$ converges to $1$ a.s., $y(t)$ converges to $0$ a.s., and $z(t)$ converges a.s. to $\pi_2$ weakly.

● If $\lambda > 0$ then the solution $u(t)$ is strongly stochastically persistent in the sense that the solution converges to its unique invariant probability measure $\mu_{3}$ supported by $\mathbb{R}_{+}^{3, \circ}$.

Case 2. When $\zeta\ge 0$, there are two ergodic invariant measures $\mu_1$ and $\mu_2$ for solutions of (2.3) on the boundary $\partial \mathbb{R}_{+}^3$.

● If $\lambda < 0$ and $\tau_1 < \sqrt{2r}$, then $x(t)$ converges to $1$ a.s., $y(t)$ converges to $0$ a.s., and $z(t)$ converges a.s. to $\pi_2$ weakly.

● If $\lambda < 0$ and $\tau_1 > \sqrt{2r}$, then solutions starting near the interior of $\mathit{\text{supp}}(\mu_2)$ will tend to stay close and concentrate on $\mathit{\text{supp}}(\mu_1)$.

● If $\lambda > 0$ and $\tau_1 < \sqrt{2r}$, then the solution $u(t)$ is strongly stochastically persistent.

● If $\lambda > 0$ and $\tau_1 > \sqrt{2r}$, then $x(t)$ and $y(t)$ both converge a.s. to $0$.

Proposition 2.1. The deterministic part of the system (2.3) has three possible nonnegative equilibrium solutions, $E_{1} = (0, 0, \frac{b}{c-1})$, $E_{2} = (1, 0, \frac{1}{2a}((1-a-c)+\sqrt{(1-a-c)^{2}+4ab}))$, and $E_{3} = (\frac{1}{a{\bar z}}, \frac{r(1-/a{\bar z})}{r+a{\bar z}}, {\bar z})$ where ${\bar z} = \frac{b-1}{c}$. The ergodic invariant measures $\mu_{1}$ and $\mu_{2}$ correspond to $E_{1}$ and $E_{2}$, respectively, in the sense that the means of the distributions of $\mu_{1}$ and $\mu_{2}$ approaches $E_{1}$ and $E_{2}$, respectively, when $(\tau_{1}, \tau_{2}, \tau_{3})$ approaches $(0, 0, 0)$.

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 $\mathbb{R}_{+}^{3, \circ}: = \{(x, y, z):x > 0, y > 0, v > 0\}$ such that $x+y\le 1$, which corresponds to the system (2.2). Then, according to the thresholds $\zeta$ and $\lambda$, we can describe how each population will evolve as follows.

Case 1. When $\zeta < 0$, the tumor cannot be eradicated completely a.s.

● If $\lambda < 0$ then $(x(t), y(t), v(t))$ converges to ${\bar \mu}_{2} = \delta_1^{*}\times\delta_0^{*}\times\delta_0^{*}$ a.s.

● If $\lambda > 0$ then $(x(t), y(t), v(t))$ is strongly stochastically coexistence in the sense that $(x(t), y(t), v(t))$ converges to a unique invariant probability measure ${\bar \mu}_{3}$ supported by $\mathbb{R}_{+}^{3, \circ}$.

Case 2. When $\zeta\ge 0$, there is some possibilities to eradicate the tumor by oncolytic viruses.

● If $\lambda < 0$ and $\tau_1 < \sqrt{2r}$, then $(x(t), y(t), v(t))$ converges to ${\bar \mu}_{2} = \delta_1^{*}\times\delta_0^{*}\times\delta_0^{*}$ a.s.

● If $\lambda > 0$ and $\tau_1 < \sqrt{2r}$, then the solution $(x(t), y(t), v(t))$ is strongly stochastically coexistence.

● If $\lambda > 0$ and $\tau_1 > \sqrt{2r}$, then $(x(t), y(t), v(t))$ converge to ${\bar \mu}_{1} = \delta_0^{*}\times\delta_0^{*}\times\delta_0^{*}$ a.s.

Using the transformation of variables or directly deduce, we have a similar proposition as 2.1.

Proposition 2.2. The deterministic part of the system (2.2) has three equilibrium solutions, $Q_{1} = (0, 0, 0)$, $Q_{2} = (1, 0, 0)$, and $Q_{3} = (\frac{1}{a{\bar z}}, \frac{r(1-/a{\bar z})}{r+a{\bar z}}, \frac{r(1-/a{\bar z})}{a+r/{\bar z}})$. The system (2.2) has three ergodic invariant probability measures ${\bar \mu}_{1} = \delta_0^{*}\times\delta_0^{*}\times\delta_0^{*}$, ${\bar \mu}_2 = \delta_1^{*}\times\delta_0^{*}\times\delta_0^{*}$, and ${\bar \mu}_{3}$, which correspond to $Q_{1}$, $Q_{2}$, $Q_{3}$, respectively.

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}$. $Q_{1}$ is always unstable for any positive values of parameters. $Q_{2}$ is globally asymptotically stable when the virus burst size $b$ is smaller a threshold value $b_{s_{1}}$, while it is unstable if $b$ is greater than $b_{s_{1}}$. There is a second threshold value of the viral burst size $b_{s_{2}}$, and under the second threshold value and other conditions, $Q_{3}$ is locally asymptotically stable. The system (1.1) undergoes Hopf bifurcations with three families of periodic solutions when the virus burst size passes the second threshold value $b_{s_{2}}$. It is interesting that $Q_{3}$ can be approximated by $(O(\frac{1}{b}), O(\frac{1}{b}), \frac{r}{a})$ when the viral burst size $b$ is very big. After incorporating environmental noise and stochastic effects into the system (1.1), there are three invariant probability measures in which solutions will approach them under various conditions. We have two combined parameters $\zeta$ and $\lambda$ to describe asymptotical properties of solutions to the systems (2.2) or (2.3) as in Theorem 2.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 $\lambda$.

Proposition 2.3. The parameter $\lambda = \sqrt{ab}\, R_{\theta}(w)-1-\frac{\tau_2^2}{2}$ is an increasing function of the virus burst size $b$. Also consider $\lambda$ is a function of noise intensities $\tau_2$ and $\tau_3$, and set ${\bar \lambda}: = \lim_{(\tau_2, \tau_3)\to(0, 0)}\lambda$. Then ${\bar \lambda} = 0$ if and only if $b = b_{s_1} = 1+\frac{c}{a}$, ${\bar \lambda} < 0$ if and only if $b < b_{s_1}$, and ${\bar \lambda} > 0$ if and only if $b > b_{s_1}$; or simply, ${\bar \lambda}$ also is a increasing function of $b$.

The parameter $\zeta$ combines infected tumor cell lysis rate $\delta$, virus degradation rate $\gamma$, and their stochastic variation $\tau_{2}$ and $\tau_{3}$, which describes possibilities if the tumor can be eradicated. More specifically, $\zeta = 2c-2-\tau_{2}^{2}+\tau_{3}^{2} = 2\frac{T_{\delta}}{T_{\gamma}}-2-\tau_{2}^{2}+\tau_{3}^{2}$, where $T_{\delta}$ is the average life time of infected tumor cells, and $T_{\gamma}$ is the average life time of free viruses in tumor tissue. $\zeta < 0$ means $\frac{T_{\delta}}{T_{\gamma}}+\frac{\tau_{3}^{2}}{2} < \frac{\tau_{2}^{2}}{2}+1$. We may interpret that, if the ratio between the life time of infected tumor cells to the life time of free viruses is small and stochastic effects of viruses also is small comparing with stochastic effects of infected tumor cells, it is impossible to eradicate the tumor for viral therapy. However, in this situation, the viral therapy may partly success which depends on $\lambda$, or implicitly the outcome of the virotherapy depends on the virus burst size $b$. As in the deterministic model (1.1), if $b$ is smaller than the threshold $b_{s_{1}}$ which corresponds to $\lambda < 0$ (it is deduced from continuity of $\lambda$ as a function of $b$, $\tau_{2}$, and $\tau_{3}$), then the infected tumor cell population and virus population will disappear, and only tumor cell population is left a.s., or the system approaches the invariant probability measure ${\bar \mu}_{2} = \delta_1^{*}\times\delta_0^{*}\times\delta_0^{*}$. If $b$ is greater than the threshold $b_{s_{1}}$ which corresponds to $\lambda > 0$, three populations will coexist, or the system approaches the invariant probability measure ${\bar \mu}_{3}$, in which we may say viral therapy achieve some partial success.

When $\zeta > 0$ which means $\frac{T_{\delta}}{T_{\gamma}}+\frac{\tau_{3}^{2}}{2} > \frac{\tau_{2}^{2}}{2}+1$. We may interpret that, if the ratio between the life time of infected tumor cells to the life time of free viruses is big and stochastic effects of viruses also is big comparing with stochastic effects of infected tumor cells, there is some possibilities to eradicate the tumor by viral therapy. In this case, there is a third threshold value for stochastic variations of tumor cell growth $\tau_{1}$ that comes to play some roles. This value is $2r = 2\frac{\rho}{\delta}$, scaled tumor cell growth rate. When $b$ is smaller than the threshold $b_{s_{1}}$ which corresponds to $\lambda < 0$, the noise intensity $\tau_{1}$ or tumor cell variance $\tau^{2}_{1}$ is smaller than the double of the scaled tumor cell growth rate, then the viral treatment will completely fail. When $b$ is greater than the threshold $b_{s_{1}}$ which corresponds to $\lambda > 0$, and the noise intensity $\frac{1}{2}\tau^{2}_{1}$ is not strong or smaller than scaled tumor cell growth rate, the system eventually will have three populations coexist, where the viral therapy reaches partial success. However, unlike in the corresponding deterministic model (1.1), when $b$ is greater than the threshold $b_{s_{1}}$ which corresponds to $\lambda > 0$, and the noise intensity $\frac{1}{2}\tau^{2}_{1}$ is strong or greater than scaled 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.

3.   Analysis of the model

This section is devoted to proving results in Section 2.

3.1.   Proof of Theorem 2.1

Consider the system (2.3). Since the drift term $f(u)$ and the diffusion term $g(u)$ are locally Lipschitz continuous, there exists a unique local a.s. continuous solution $u(t)$ up to the explosion time

Also, the solution $u(t) = (x(t), y(t), z(t))$, $t\in(0, \tau_e)$, is a strong Markov process (see ^[24]). Denote by $(x, y, z)$ the initial value of $u(t)$. First, we will show that if $(x, y, z)$ is in $\mathbb{R}_{+}^3$ then $u(t)$ is also in $\mathbb{R}_{+}^3$ for all $t\in(0, \tau_e)$ a.s.

From the equation of $x(t)$, we get

So, if $x = 0$, then $\mathbb{P}_{0, y, z}\{x(t) = 0\, \, \forall\, t\in(0, \tau_e)\} = 1$ for all $y\ge 0$ and $z\ge 0$; if $x > 0$, then $\mathbb{P}_{x, y, z}\{x(t) > 0\, \, \forall\, t\in(0, \tau_e)\} = 1$ for all $y\ge 0$ and $z\ge 0$.

The second equation of (2.3) implies

If $y = 0$, then $\mathbb{P}_{x, 0, z}\{y(t) = 0\, \, \forall\, t\in(0, \tau_e)\} = 1$ for all $x\ge 0$ and $z\ge 0$; if $y > 0$, then $\mathbb{P}_{x, y, z}\{y(t) > 0\, \, \forall\, t\in(0, \tau_e)\} = 1$ for all $x\ge 0$ and $z\ge 0$.

The last equation of (2.3) follows

where

This implies that if $z\ge 0$, then $\mathbb{P}_{x, y, z}\{z(t) > 0\, \, \forall\, t\in(0, \tau_e)\} = 1$ for all $x\ge 0$ and $y\ge 0$.

Hence, we have shown that if $x\ge 0$, $y\ge 0$, and $z\ge 0$ then $x(t)\ge 0$, $y(t)\ge 0$, and $z(t)\ge 0$ for all $t\in(0, \tau_e)$ a.s.

Next, we show that $\tau_e = \infty$ a.s. Consider $V(x, y, z) = x+y+\ln(1+z)$. By Ito's formula, for all $t\in(0, \tau_e)$ we get

Let $\tau_n: = \inf\{t\in[0, \tau_e):x(t) > n\, \, \text{or}\, \, y(t) > n\, \, \text{or}\, \, z(t) > n\}$. Clearly, $\tau_n$ increases to $\tau_{\infty}$ as $n\to\infty$ where

Since $\tau_{\infty}\le\tau_e$ a.s., it suffices to prove that $\mathbb{P}_{x, y, z}\{\tau_{\infty} = \infty\} = 1$. Fix $t > 0$, Ito's formula for $V$ implies

where $K: = V(x, y, z)$. On the other hand,

Thus

Since $t > 0$ is arbitrary, $\mathbb{P}_{x, y, z}\{\tau_{\infty} < \infty\} = 0$ and hence $\tau_{\infty} = \infty$ a.s.

This completes the proof.

3.2.   Proof of Theorem 2.2

Before giving the detailed proof of the main theorem 2.2, we analyze solutions of the system (2.3) on the boundary $\partial \mathbb{R}_{+}^3$ firstly. When $x(0) = 0$, $x(t) = 0$ for all $t\ge 0$ a.s. If $x = 0$ then the system (2.3) becomes

The second equation for $y(t)$ implies

So $Y(t)\to 0$ a.s. for all $Y(0) = y(0)\ge 0$. Consider the last equation for $Z$

Fix $\alpha_1 > 0$, consider

where $C_1$ is some positive constant. Rewrite the integrand as

Since there exists a $k\in \mathcal{Z}_{+}$ such that $-\frac{2(1+\tau_2^2-c)}{\tau_2^2+\tau_3^2}-k < -1$, $s(0+) = -\infty$. If $\zeta = 2c-2-\tau_2^2+\tau_3^2 < 0$, then $-\frac{2(1+\tau_2^2-c)}{\tau_2^2+\tau_3^2}+1 < 0$, and, so $s(\infty) < \infty$. By the item 2 of Theorem 3.1 on page 447 in ^[25], $\lim\limits_{t\to\infty}Z(t) = \infty$ a.s. In this case, (3.1) does not have any invariant measure. If $\zeta\ge 0$, then $-\frac{2(1+\tau_2^2-c)}{\tau_2^2}+1\ge 0$, and this implies that $s(\infty) = \infty$. Then $Z(t)$ oscillates between $0$ and $\infty$. Then (3.1) has a unique invariant measure $\pi_1\sim\text{IG}\left(\frac{2(c-1-\tau_2^2)}{\tau_2^2+\tau_3^2}+1, \frac{2b}{\tau_2^2+\tau_3^2}\right)$ (the inverse gamma distribution with parameters $\frac{2(c-1-\tau_2^2)}{\tau_2^2+\tau_3^2}+1$ and $\frac{2b}{\tau_2^2+\tau_3^2}$).

When $x(0) > 0$, $x(t) > 0$ for all $t > 0$ a.s. If $y(0) = 0$, then the second equation of (2.3) implies $y(t) = 0$ for all $t > 0$ a.s. So, when $y = 0$, the equation for $x$ becomes

Fix $\alpha_2 > 0$, we compute

where $C_2$ is some positive constant. Clearly, $s(1-) < \infty$. Since $\lim_{y\to 0+}(\frac{1}{y}-1)^{2r/\tau_1^2} = \infty$, for any $M > 0$, there exists a $0 < \delta < \alpha_2$ so that, if $0 < y < \delta$, then $(\frac{1}{y}-1)^{2r/\tau_1^2}\ge\frac{M}{\delta C_2}$. But we have

Letting $M\to\infty$, it gives $s(0+) = -\infty$. This means that $\lim_{t\to\infty}\tilde{x}(t) = 1$ a.s. for any $\tilde{x}(0) = x(0) > 0$. When $x = 1$ and $y = 0$, the last equation for $z$ becomes

Fix $\alpha > 0$, consider

where $C_3$ is some positive constant. The integrand can be written as

Clearly, there are $k_1$ and $k_2$ in $\mathcal{Z}_{+}$ such that

Hence $s(0+) = -\infty$ and $s(\infty) = \infty$. So $\tilde{z}(t)$ oscillates between $0$ and $\infty$, and thus (3.3) has a unique invariant measure $\pi_2\sim\text{GIG}(\theta, \chi, \psi)$, which is the generalized inverse Gaussian distribution with parameters $\theta\in \mathbb{R}$, $\chi > 0$, and $\psi > 0$ (see the Appendix), whose density takes the form

where $\theta: = \dfrac{2(1+\tau_2^2-a-c)}{\tau_2^2+\tau_3^2}-1$, $\chi: = \dfrac{4b}{\tau_2^2+\tau_3^2}$, and $\psi: = \dfrac{4a}{\tau_2^2+\tau_3^2}$; $K_{\theta}(\cdot)$ is the modified Bessel function of the third kind with index $\theta$. By law of large numbers,

in which $R_{\theta}(w): = \dfrac{K_{\theta+1}(w)}{K_{\theta}(w)}$ and $w: = \dfrac{4\sqrt{ab}}{\tau_2^2+\tau_3^2}$.

In summary, on the boundary $\partial \mathbb{R}_{+}^3$

● If $\zeta = 2c-2-\tau_2^2+\tau_3^2 < 0$, then the system (2.3) has only one invariant measure $\mu_2: = \delta_1^{*}\times\delta_0^{*}\times\pi_2$.

● If $\zeta\ge 0$ then the system (2.3) has two invariant measures $\mu_1: = \delta_0^{*}\times\delta_0^{*}\times\pi_1$ and $\mu_2$.

Note that

We define a combined parameter as the threshold

and define the family of the random normalized occupation measures

Then, we have the following claim.

Claim 3.1. Assume that $\lambda < 0$. For any initial value $u = (x, y, z)$ in $\mathbb{R}_{+}^{3, \circ}$ satisfying $x+y\le 1$, if $\Pi_t(\cdot)$ converges weakly to $\mu_2$ a.s. and $y(t)$ converges a.s. to $0$ exponentially fast with the rate $\lambda$, then $x(t)$ converges a.s. to $1$ and $z(t)$ converges weakly to $\pi_2$.

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.

Lemma 3.1. If $u(t) = (x(t), y(t), z(t))$ is the solution of (2.3) with the initial value $u = (x, y, z)$ satisfying $x > 0$, $y > 0$, $z > 0$, and $x+y\le 1$ then $0 < x(t) < 1$ for all $t > 0$ a.s.

Lemma 3.2. Suppose the assumption of Lemma 3.1 is satisfied. Then $\limsup\limits_{t\to\infty} \mathbb{E}_u x(t)z(t) < \infty$.

Proof of Claim 3.1. From the equation for $x(t)$ of (2.3), we get

Denote

Clearly, $A_t$ and $U_t$ are continuous adapted ($\mathcal{F}_t$-measurable) increasing processes with $A_0 = U_0 = 0$. $M_t$ is a local martingale with $M_0 = 0$ and $1-x(t)\ge 0$ a.s. by Lemma 3.1.

Then, we show that $\lim_{t\to\infty}A_t < \infty$ a.s. Since $\lim\limits_{t\to\infty}\dfrac{\ln y(t)}{t} = \lambda < 0$, there is a $\Theta > 0$ such that $t\ge\Theta$ implies $y(t)\le\exp\{\frac{\lambda t}{2}\}$. But, then for $t\ge\Theta$

which follows that $\lim\limits_{t\to\infty}\int_0^t rx(s)y(s)ds < \infty$ a.s. On the other hand, by Lemma 3.2, we can use the Markov inequality to show $N: = \sup\limits_{t\ge 0}x(t)z(t) < \infty$ a.s. So the same argument as above implies $\lim\limits_{t\to\infty}\int_0^t ax(s)y(s)z(s)ds < \infty$ a.s. Therefore $\lim\limits_{t\to\infty}A_t < \infty$ a.s.

By Theorem 3.9 on page 14 in ^[24], $\lim\limits_{t\to\infty}(1-x(t)) < \infty$ a.s. and $\lim\limits_{t\to\infty}\int_0^t x(s)[1-x(s)]ds < \infty$ a.s.

If $x(t)[1-x(t)]$ did not converge a.s. to $0$, then there would be an $\Omega_1\subseteq\Omega$ with $\mathbb{P}(\Omega_1) > 0$ so that $\liminf\limits_{t\to\infty}x(t, \omega)[1-x(t, \omega)] = p(\omega) > 0$ for all $\omega\in\Omega_1$. Fix $\omega\in\Omega_1$, there exists a $T: = T(\omega) > 0$ so that $t\ge T$ implies $x(t, \omega)[1-x(t, \omega)] > \frac{1}{2}p(\omega)$. Hence

Then $\Omega_1\subseteq\Omega_2$, where $\Omega_2 = \left\{\omega; \, \int_0^{\infty}x(s, \omega)[1-x(s, \omega)]ds = \infty\right\}$. This implies that $\mathbb{P}(\Omega_2) > 0$. But this contradicts the fact that $\lim\limits_{t\to\infty}\int_0^t x(s)[1-x(s)]ds < \infty$ a.s. Therefore

Since $\Pi_t(\cdot)$ converges weakly to $\mu_2$, there exists a sequence $\{t_k\}_{k\ge 1}$ such that $t_k\uparrow\infty$ and

where $P(t, u, \cdot)$ is the transition probability of the solution $u(t)$ of the system (2.3). In other words, $\mathbb{E}_u x(t_k)\to 1$. Combining this fact with (3.4), we can conclude that $x(t)$ converges a.s. to $1$. Moreover, since

$\mathbb{E}_u|z(t_k)-\tilde{z}(t_k)|\to 0$. As $\tilde{z}(t)$ converges weakly to $\pi_2$, so is $z(t)$.

Now, we give a proof of our main theorem 2.2. Notice that Assumptions 1.1-1.5 and Theorems 1.1 and 1.3 mentioned in the proof are referred in ^[23], since our proof is based on Theorem 1.1 and Theorem 1.3 there.

Proof of Theorem 2.2. First, we denote $x_1 = x$, $x_2 = y$, $x_3 = z$, $f_1 = r(1-x-y)-ayz$, $f_2 = axz-1$, $f_3 = \frac{b}{z}+1+\tau_2^2-ax-c-axz$, $g_1 = 1-x-y$, $g_2 = g_3 = 1$, and

It is clear that $f_i$ and $g_i$ ($i = 1, 2, 3$) are locally Lipschitz. For $c = (0, 0, 1)^T$, $u = (x, y, z)^T$, and $\gamma_b > 0$, we have

Note that, since $0\le x(t)\le 1$ for all $t\ge 0$ a.s. by Lemma 3.1, we can show that $\mathbb{E}_u(x(t)+y(t))\le 1$. This implies that $\mathbb{E}_uy(t)\le 1$ and hence, by Markov's inequality, we can prove that $y(t)$ is bounded a.s. Thus

for some constant $K_4 > 0$. When $z$ is large enough, we can choose $\gamma_b > 0$ sufficiently small so that

This shows that

where $\|u\|: = |x|+|y|+|z|$. Moreover, it is easy to compute

which is positive definite for all $(x, y, z)\in \mathbb{R}_{+}^{3, \circ}$ satisfying $x+y\le 1$. Thus, Assumption 1.1 in ^[23] is fulfilled for the system (2.3).

Next, let $\mathcal{M}$ be the set of ergodic invariant measures of the system (2.3) supported by the boundary $\partial \mathbb{R}_{+}^3$. Consider two cases.

Case 1. $\zeta = 2c-2-\tau_2^2+\tau_3^2 < 0$. There is only one ergodic invariant measure on $\partial \mathbb{R}_{+}^3$, which is $\mu_2 = \delta_1^{*}\times\delta_0^{*}\times\pi_2$. Observe that

Then $\mathcal{M} = \{\mu_2\}$, and so $\text{Conv}(\mathcal{M}) = \{\mu_2\}$ (the convex hull of $\mathcal{M}$, that is the set of probability measure $\pi$ of the form $\pi(\cdot) = \sum_{\mu\in\mathcal{M}}p_{\mu}\mu(\cdot)$ with $\sum_{\mu\in\mathcal{M}}p_{\mu} = 1$, $p_{\mu}\ge 0$). If $\lambda > 0$, then Assumption 1.2 holds and thus, by Theorem 1.1, the solution of (2.3) is strongly stochastically persistent. If $\lambda < 0$ then Assumption 1.3 holds. Note that

So

and hence $\mathcal{M}^2: = \mathcal{M}\backslash\mathcal{M}^1 = \emptyset$. This means that Assumption 1.5 is satisfied. Furthermore, since $\sum_{i = 1}^3 g_i^2 = 2+(1-x-y)^2$ is bounded, for any $0 < \delta_1 < 1$ we have

Thus Assumption 1.4 is fulfilled. By Theorem 1.3, for any initial value $(x, y, z)$ in $\mathbb{R}_{+}^{3, \circ}$ satisfying $x+y\le 1$, $\Pi_t(\cdot)$ converges weakly to the unique invariant measure $\mu_2$ and $\lim_{t\to\infty}\frac{\ln y(t)}{t} = \lambda$ w.p.1. By Claim 3.1, $x(t)$ converges a.s. to $1$ and $z(t)$ converges weakly to $\pi_2$ a.s.

Case 2. $\zeta\ge 0$. There are two ergodic invariant measures on $\partial \mathbb{R}_{+}^{3}$, which are $\mu_1$ and $\mu_2$. It is straightforward to see that

Then $\mathcal{M} = \{\mu_1, \mu_2\}$ and so

If $\tau_1 < \sqrt{2r}$ and $\lambda > 0$ then $\lambda_1(\mu_1) > 0$ and $\lambda_2(\mu_2) > 0$. It is clear that, for any $\mu\in\text{Conv}(\mathcal{M})$, $\max\{\lambda_1(\mu), \lambda_2(\mu), \lambda_3(\mu)\} > 0$. Then Assumption 1.2 holds and hence, by Theorem 1.1, the solution of (2.3) is strongly stochastically persistent.

If $\tau_1 < \sqrt{2r}$ and $\lambda < 0$ then $\mathcal{M}^1 = \{\mu_2\}$. Since $\mathcal{M}_{\mu_2} = \{\mu_1\}$ and $\text{Conv}(\mathcal{M}_{\mu_2}) = \{\mu_1\}$, $\max\{\lambda_1(\mu_1), \lambda_3(\mu_1)\} > 0$. Hence Assumption 1.3 holds. Since $\mathcal{M}^2 = \{\mu_1\}$, $\text{Conv}(\mathcal{M}^2) = \{\mu_1\}$ and so $\max\limits_{i = 1, 2, 3}\lambda_i(\mu_1) > 0$. Then Assumptions 1.4 and 1.5 are satisfied. By Theorem 1.3 and Claim 3.1, $x(t)$ converges a.s. to $1$, $y(t)$ converges a.s. to $0$, and $z(t)$ converges weakly to $\pi_2$ a.s.

If $\tau_1 > \sqrt{2r}$ and $\lambda > 0$ then, since $\max\limits_{i\in I_{\mu_1}^c}\lambda_i(\mu_1) < 0$ and $\mathcal{M}_{\mu_1} = \emptyset$, $\mathcal{M}^1 = \{\mu_1\}$ and so $\mathcal{M}^2 = \{\mu_2\}$. As $\max\limits_{i = 1, 2, 3}\lambda_i(\mu_2) = \lambda > 0$, so Assumption 1.5 are satisfied. With the same argument as Case 1, we also have Assumption 1.4 is fulfilled. Thus, by Theorem 1.3, we can conclude that $x(t)$ and $y(t)$ both converge a.s. to $0$ with the rates $\lambda_1(\mu_1)$ and $\lambda_2(\mu_1)$, respectively.

Lastly, if $\tau_1 > \sqrt{2r}$ and $\lambda < 0$, then $\max\limits_{i\in I_{\mu_1}^c}\lambda_i(\mu_1) < 0$ and $\max\limits_{i\in I_{\mu_2}^c}\lambda_i(\mu_2) < 0$. However, the condition $(1.7)$ in Assumption 1.3 does not hold because $\max\limits_{i\in I_{\mu_1}}\lambda_i(\mu_1) < 0$ and $\max\limits_{i\in I_{\mu_2}}\lambda_i(\mu_2) = 0$. This means that $\mu_1$ and $\mu_2$ are not repellers. In this case, solutions starting near the interior of $D_{\mu_2}$ will stay close and concentrate on $D_{\mu_1}$.

This completes the proof.

Proof of Lemma 3.1.

Take $V(x, y, z) = 2-x-y-\ln(1-x)$ for $x > 0$, $y > 0$, $z > 0$ and $x+y\le 1$. By Ito's formula, for all $x > 0$, $y > 0$, $z > 0$, and $x+y\le 1$ we have

for some suitable positive constant $K_1$. Let $\zeta_k = \inf\{t\ge 0:V(x(t), y(t), z(t))\ge k\}$ and fix $t > 0$. Then Ito's formula implies

By Gronwall's inequality, $\mathbb{E}_{x, y, z}V(t)\le V(x, y, z)\exp\{K_1t\}$. But, since

$\mathbb{P}_{x, y, z}\{\zeta_k\le t\}\le\dfrac{V(x, y, z)\exp\{K_1t\}}{k}$ for all $k\ge 1$ and hence

for all $k\ge 1$. On the other hand, since $\zeta_k > t$ implies $V(x(s), y(s), z(s)) < k$ for all $s\in[0, t]$,

for all $k\ge 1$. Letting $k\to\infty$ yields

As $V(x(s), y(s), z(s)) < \infty$ implies $1-x(s) > 0$, so $\mathbb{P}_{x, y, z}\{0 < x(s) < 1\, \, \forall\, s\in[0, t]\} = 1$. Since $t > 0$ is arbitrary, $\mathbb{P}_{x, y, z}\{0 < x(s) < 1\, \, \forall\, s\ge 0\} = 1$.

We complete the proof.

Proof of Lemma 3.2.

By Ito's formula, since $0\le x(t)\le 1$ for all $t\ge 0$ a.s. (by Lemma 3.1),

Taking expectation both sides yields

here we have used the equality $\mathbb{E}_ux^2(t)z^2(t)\ge(\mathbb{E}_ux(t)z(t))^2$. This implies that

Therefore for all $t\ge 0$ we obtain

So, Lemma 3.2 is proved.

3.3.   Proofs of Propositions

Proof of Proposition 2.1.

$\mu_1 = \delta_0^{*}\times\delta_0^{*}\times\pi_1$, where the invariant measure $\pi_1$ has the inverse gamma distribution with parameters $\alpha = \frac{2(c-1-\tau_2^2)}{\tau_2^2+\tau_3^2}+1$ and $\beta = \frac{2b}{\tau_2^2+\tau_3^2}$. The mean is given by $\frac{\beta}{\alpha -1} = \frac{b}{c-1-\tau_{2}^{2}}$. It is clearly that when $\tau_{2}$ approaches zero, the mean approaches $\frac{b}{c-1}$ which is the third coordinate of the equilibrium $E_{1}$.

$\mu_2 = \delta_1^{*}\times\delta_0^{*}\times\pi_2$, where the invariant measure $\pi_2$ has the generalized inverse Gaussian distribution with parameters $\theta = \dfrac{2(1+\tau_2^2-a-c)}{\tau_2^2+\tau_3^2}-1$, $\psi = \dfrac{4a} {\tau_2^2+\tau_3^2}$, and $\chi = \dfrac{4b}{\tau_2^2+\tau_3^2}$. From A.17 on page 172 in ^[26], we know the mean of this distribution $\pi_2$ is $R_{\theta}(w)\sqrt{\frac{b}{a}}$, where $w: = \sqrt{\psi\chi} = \frac{4\sqrt{ab}}{\tau_2^2+\tau_3^2}$, $R_{\theta}(w): = \frac{K_{\theta+1}(w)}{K_{\theta}(w)}$, and $K_{\theta}(w)$ is the modified Bessel function of the third kind with index $\theta$.

Since $\frac{\theta}{w} = \frac{2(1+\tau_3^2-a-c)}{4\sqrt{ab}}-\frac{\tau_2^2+\tau_3^2}{4\sqrt{ab}}$, $\lim_{(\tau_2, \tau_3)\to(0, 0)}\frac{\theta}{w} = \frac{1-a-c}{2\sqrt{ab}}$. From the reference ^[26], we get

Due to the asymptotic expansion of $D_{\theta}(w)$ as $w\to\infty$ (see A.22 on page 173 in ^[26])

it is clear that $D_{\theta}(w)$ approaches $1$ as $(\tau_2, \tau_3)$ approaches $(0, 0)$. Hence

Proof of Proposition 2.3.

We first show the threshold $\lambda = \sqrt{ab}\, R_{\theta}(w)-1-\frac{\tau_2^2}{2}$ is an increasing function of the virus burst size $b$. Let $h(b) = \sqrt{ab}\, R_{\theta}(w)$ where $w = \frac{4\sqrt{ab}}{\tau_2^2+\tau_3^2}$. Since $h'(b) = \frac{2a}{\tau_2^2+\tau_3^2}\left[R^2_{\theta}(w)-2\frac{\theta}{w}R_{\theta}(w)-1\right]$, $h'(b)\ge 0$ is equivalent to $R_{\theta}(w)\ge\frac{\theta}{w}+\sqrt{(\frac{\theta}{w})^2+1}$. By the integral representation of $D_{\theta}(w)$ (see A.29 on page 176 in ^[26]), we have

It follows that

Therefore $h(b)$ is increasing w.r.t. $b$.

From the proof of Proposition 2.1, $\lim_{(\tau_2, \tau_3)\to(0, 0)}\lambda = \lim_{(\tau_2, \tau_3)\to(0, 0)}\sqrt{ab}R_{\theta}(w)-1-\frac{1}{2}\tau_{2}^{2} = \frac{1-a-c+\sqrt{(1-a-c)^2}+4ab}{2}-1 = {\bar \lambda}$. It is easy to see that ${\bar \lambda} = 0$ gives $b = 1+\frac{c}{a}$ which is the threshold $b_{s_{1}}$ for the deterministic system (1.1).

4.   Numerical simulation and discussion

4.1.   Numerical simulation with biological explanation

This section is devoted to demonstrate our main analytical results in Section 2. Data values from our previous research (see ^[27] and ^[12]) are used to estimate parameter values and simulate our stochastic model. The maximal radius of the tumor in mice brain, which is considered to be dead from the tumor, is $5$ millimeters. Because our SDE model neglects spatial variations, tumor size is converted into cell numbers using the constant of cell density $K = 10^6$ per cubical millimeter. After non-dimensionalization, the parameter values are $r = 0.36$, $a = 0.11$, and $c = 0.44$. For simplicity, we carry out numerical simulations based on the non-dimensionalized SDE systems (2.2) and (2.3). The quantities $x$, $y$, and $v$ are, respectively, the portion of uninfected tumor cells, infected tumor cells, and free virus particles over the maximal cell density of the tumor $C$. These quantities are not absolute numbers but relative numbers. We just call them relative uninfected tumor cells and so on in the figures below. Notice that the quantity $z$ in the system (2.3) is the ratio of relative free virus particles over relative infected tumor cells. In our simulation, time is scaled or relative time since $T = \delta 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 $\lambda < 0$. In the figures 1 and 2, we take $b = 5$, $\tau_1 = 0.2$, $\tau_2 = 0.3$, and $\tau_3 = 0.2$. By computation, $\theta = 7.3077$, $w = 22.8191$, $R_{\theta}(w) = 1.39$, and hence $\lambda = -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.

Example 2. We consider the situation when $\lambda > 0$. Take $b = 10$ and noise intensities have the same values as in Example 1. By computation, we can obtain $\lambda = 0.2832 > 0$. The second conclusion in Case 1 of Theorem 2.2 indicates that relative populations of uninfected, infected tumor cells and free viruses persist strongly and finally settle down into a coexistence equilibrium state. Both Figures 3 and 4 show that, as time goes by, each solution path oscillates most of the time around a positive equilibrium point (which is the positive equilibrium point of the corresponding ODE system, $E_{3}$ or $Q_{3}$). Biologically, this 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 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.

4.2.   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 $\zeta$ and $\lambda$ that provide conditions for various outcomes of our stochastic model. The parameter $\zeta$ combines infected tumor cell lysis rate $\delta$, virus degradation rate $\gamma$, and their stochastic variation $\tau_{2}$ and $\tau_{3}$, which describes possibilities if the tumor can be eradicated by viral therapy. The parameter $\lambda$ is a differential function of the viral burst size $b$ and the noise intensities $\tau_2$ and $\tau_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 $\tau_{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 corresponds to 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 $\lambda$ 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.

Acknowledgments

JPT would like to acknowledge grant supports from National Science Foundation (DMS-1446139) and National Institutes of Health (U54CA132383) which also supported TAP in some semesters and summers during the grant periods.

Conflict of interest

All authors declare that there are no conflicts of interest regarding the publication of this manuscript.

Appendix

We state some properties of the Generalized Inverse Gaussian Distribution GIG($\theta, \chi, \psi$) without proof (see ^[26] for more details), that is needed in Section 3, where $\theta\in \mathbb{R}$, $\chi > 0$, and $\psi > 0$. Its probability density function is given by

where $K_{\theta}(\cdot)$ is the modified Bessel function of the third kind with index $\theta$, given by

Let $w = \sqrt{\chi\psi}$ and $\eta = \sqrt{\dfrac{\chi}{\psi}}$, then the above probability density function takes the form

Moments of random variable $X\sim\text{GIG}(\theta, w, \eta)$ are given by