A mathematical model of oncolytic virotherapy with time delay

1.   Introduction

Viruses that selectively replicate in tumor cells have recently demonstrated their potential use in cancer treatment ^[1,2,3,4]. Replicating oncolytic viruses are able to infect and lyse cancer cells and spread through the tumor while leaving normal cells largely unharmed. A variety of viruses have shown promising results in clinical trials ^[5]. Among the oncolytic viruses with potential use for virotherapy are the adenovirus Onyx-015 ^[6], the herpes simplex virus HSV-1 ^[7], the Newcastle disease virus NDV ^[8] and M1 virus ^[9].

In order to have a complete understanding of how virus and host characteristics influence the outcome of therapy, many mathematical models have been established. For example, in 2001, Wodarz ^[10] established a mathematical model as follows.

In this model, there are three variables, uninfected tumor cells $x$, infected tumor cells $y$ by the viruses and virus-specific CTL $z_v$. The tumor cells grow in a logistic fashion at a rate $r$ and die at a rate $d$. The maximum size of space where the tumor is allowed to occupy is given by its carrying capacity $k$. The viruses spread to tumor cells at a rate $\beta$. Infected tumor cells are killed by the viruses at a rate $a$ and grow in a logistic fashion at a rate $s$. The virus-specific CTL expands in response to antigen at a rate $c_v$ and decays at a rate $b$. Hereafter, without specific indication, we denote by $x, y, z_v$ uninfected tumor cells, infected tumor cells and virus-specific CTL response, respectively.

In ^[10], the authors mainly focused on the total tumor load by analyzing each equilibrium of Eq (1.1). Certain equilibria particularly drew their attention. In the absence of the virus-specific CTL response, they found

● an infection free equilibrium $E_0 = (\frac{k(r-d)}{r}, 0, 0)$,

● a 100% virus prevalence equilibrium $E_1 = (0, \frac{k(s-a)}{s}, 0)$,

● and a coexistence equilibrium with both infected and uninfected tumor cells

In the presence of the virus-specific CTL response, they found

● an equilibrium $E_3 = (0, \frac{b}{c_v}, \frac{kc_v(s-a)-sb}{p_v kc_v})$ at which there is 100% virus prevalence in the tumor cell population,

● a coexistence equilibrium with both infected and uninfected cells $E_4(x_4, y_4, z_{v4})$, where $x_4 = \frac{r(kc_v-b)-k(c_v d+b\beta)}{rc_v}, y_4 = \frac{b}{c_v}, z_{v4} = \frac{\beta k(rc_v-b\beta-c_v d)-c_v(ar-sd)-b\beta(r-s)}{p_v c_v r}.$

At equilibrium $E_0$, the tumor is at its maximum size $\frac{k(r-d)}{r}$ without virotherapy; At $E_1$, all of the tumor cells are infected and tumor size is given by $x+y = \frac{k(s-a)}{s}$. If $[\beta k(r-d)+sd]/r > a$, then the total tumor load equals to $\frac{k(r-s+a-d)}{\beta k+r-s}$ at equilibrium $E_2$. Wodarz claimed that there is an optimal choice of virus cytotoxicity $a = \frac{s(d+\beta k)}{r+\beta k}$ which led to a minimum tumor load given by $x+y = \frac{k(r-d)}{r+\beta k}$. Using the same method, He examined equilibria $E_{3}, E_{4}$ to find the optimal choice for immune response rate $c$.

We note that in ^[10], the authors did not provide rigorous mathematical proofs of the stability of equilibria mentioned above. Some experimental evidence shows that several different viruses inhibit cellular replication after infection, see for example ^[11]. Henceforth, the growth term in infected tumors is not considered in most of mathematical models, i.e. $s = 0$ ^[5,12,13].

Recently, mathematical models with intracellular viral life-cycle have been formulated ^[14]. At the molecular level, a great deal of phenomena about intracellular viral life cycles have been found experimentally. Indeed, there are several stages in a typical viral life-cycle: viral entry, viral replication, viral shedding and viral latency. For the details of the viral life-cycle, we refer the reader to ^{[14,15,16,17,18,19]}. In ^[14], the model is formulated as

where $v$ is the free virus particle, the time period of the lytic cycle is described by $\tau$, $k$ is the maximal tumor size, $r$ is the per capita tumor growth rate. The coefficient $b$ represents the infectivity of the virus. The coefficient $a$ is the death rate of infected tumor cells, and $\gamma$ is the virus clearance rate. The parameter $\delta$ is the burst size of the virus.

We point out that the variable $y$ in Eq. (1.2) is the number of infected tumor cells in the last stage of lytic cycle, which is infective and slightly different from the one in Eq. (1.1). Also, note that Eq. (1.1) are only accurate if the death rate of newly infected cells in the first stage of the lytic cycle (eclipse phase)^[20], here referred to as $I(t)$, is zero. More generally, if $n$ denotes the death rate of these infected tumor cells, then the term $b x(t-\tau)v(t-\tau)$ in the second equation of Eq. (1.2) should be replaced by $b\exp\{-n\tau\}x(t-\tau)v(t-\tau)$, where $\exp\{-n\tau\}$ represents the survival rate during this period. Observe that then $I(t) = \int_{t-\tau}^t e^{-n(t-s)}bx(s)v(s)ds$. Technically, density dependence in tumor growth should include all tumor cells, e.g. $x(t)+I(t)+y(t)$ in the equation for $x$. However, including this term would then require including an equation for $I(t)$ in the model, a serious mathematical complication. If $\tau$ is small, $I(t)$ should also be small so we may neglect $I(t)$ in the logistic term. For these reasons, we leave the logistic term as in ^[14].

We assume that the turnover of free virus is fast compared to that of infected cells. Similar to ^[21], using a quasi-steady state assumption we assume that infected tumor cell density is proportional to virus density, thereby allowing us to drop the equation for virus. Thus, Eq. (1.2) can be reduced as below:

In oncolytic virotherapy, the effect of immune response is indispensable and can not be neglected. Experiments with injecting mutant herpes simplex virus 1 (hrR3) into glioma implanted in brains of rats show the lack of efficacy in eradicating the cancer, due to interference by the immune system ^[22]. In other words, the presence of free virus will lead to the virus-specific CTL response that will kill infected tumor cells. Similar to the idea of Wodarz ^[10], virus-specific CTL response will also be involved in our model.

According to discussions above, the schematic representation of the assumptions underlying the mathematical models is drawn as follows.

Very recently, research on increasing the oncolytic effect of M1 virus indicated that a classical protein kinase A (PKA) inhibitor, H89, inhibits the innate antiviral response ^[23]. Some other immune suppressive drugs are reported in ^[1,22]. From this point, we think that the strength, defined by $\frac{c}{d}$, of virus-specific immune response can be controlled by immune suppressive drugs. In this paper, we will support that the use of immune suppressive drugs will benefit the oncolytic effect of virus as reported in clinic results ^[23].

In the next section, we propose a new mathematical model based on the schematic diagram 1. Main mathematical results are listed in Table 2 and Table 3. Numerical simulations and biological interpretations will be given in Section 3. Section 4 is devoted to the rigorous mathematical analysis of our main results. In the last section, we will give optimal strategies in tumor therapy in different situations.

2.   Materials and method

As has been stated in the previous section, our model will have the following form

where $n$ is the death rate of infected tumor cells during the period of the lytic cycle. Explanations of other parameters in Eq. (2.1) are listed in the Table 1. As $I(t)$ is uncoupled with $x(t), \ y(t), z(t)$. We just consider the mathematical model about $x(t), \ y(t), \ z(t)$.

First, we consider equilibria with no virus-specific CTL response. There are two such equilibria: the trivial equilibrium $E_0 = (0, 0, 0)$, and the original tumor equilibrium $E_1 = (k, 0, 0)$. We claim that infections by viruses in tumor cell population occur if $R_0 > 1$, where $R_0$ is defined as

Indeed, $R_0$ can be calculated in the same way as the basic reproductive number in epidemiology model ^[26]. See also Chapt. 8 in ^[27]. By the second equation of Eq. (2.1), the life time of infected component is $\frac{1}{a}$, and every infected tumor cell produces $\frac{bx_{*}\exp\{-n\tau\}}{a}$ newly infected individuals during its life time. Note that $x_{*}$ is just the maximum carrying capacity $k$ of tumor load. So, every infected individual can produce $R_0$ infective tumor cells. Thus, if $R_0 > 1$, infected cells proliferate. Thus, virus therapy equilibrium $E_{vt}(x_{vt}, y_{vt}, z_{vt})$ exists, where

Especially, if the cytotoxicity of virus $a = 0$, the virus can attain $100\%$ prevalence in the tumor cell population. In this case, the virus therapy equilibrium $E_{vt}$ is reduced to $E_* = (0, \frac{rk}{r+bk}, 0)$. If $R_0 = 1$, equilibrium $E_{vt}$ coincides with equilibrium $E_1$.

Next, we consider virus replication in the presence of the virus-specific CTL response, i.e. $z\ne 0$. In this case, there is another equilibrium $E_{CTL} = (x_{CTL}, y_{CTL},$ $z_{CTL})$, where

Similar to the definition of $R_0$, we can define

By simple calculations, we can verify that $R_1 > 1$ is a sufficient and necessary condition for the existence of positive equilibrium $E_{CTL}$, see Theorem 4.4. It implies that the virus-specific CTL component can invade the infection component. If $R_1 = 1$, equilibria $E_{vt}$ and $E_{CTL}$ coincide. Note that if $R_0\leq 1$, then $R_1\leq0$, which implies that if the oncolytic viruses fail to spread in the tumor cell, then the virus specific CTL response will not occur.

Our main mathematical results are listed in Table 2 and Table 3. Detailed proofs are given in Section 4.

3.   Biological results

Based on the mathematical results described above, in this section we will investigate the treatment outcome and its dependence on system parameters. Recall our key parameters: virus infection rate $b$, cytotoxicity of the virus $a$, intercellular viral life-cycle length $\tau$, immune cell stimulation rate by virus $c$, and the clearance rate of immune cells $d$. We focus on the stability properties of the the equilibria $E_1$, $E_{vt}$ and $E_{CTL}$ and how much the tumor load $x+y$ is reduced for the "treatment equilibria" $E_{vt}$ and $E_{CTL}$ relative to that of $E_1$ where the load is $k$. Mathematical analysis shows that success or failure of the virus therapy is determined by $R_0 = bk\exp(-n\tau)/a$. Virus therapy occurs only if $R_0 > 1$ and virus elicits an immune response only if $R_1 > 1$ where $R_1 = \frac{crk}{d(r+bk)}(1-R_0^{-1})$.

If $R_0\leq 1$, infected tumor cells will go extinct before the virus has a chance to significantly spread, see Figure 5. If $1 < R_0\leq 3, R_1 < 1$, then infected tumor cells can invade the uninfected tumor but virus CTL response is not elicited, see Figure 6. If $R_0 > 3, R_1 < 1$, sustained oscillations occur in the tumor load, see Figure 7. However, suitable virus-specific CTL response make the periodic oscillation fade away, see Figure 8. In particular, if $R_0 > 1, R_1\leq1, a = 0$, the virus will attain the $100\%$ prevalence in the tumor population. In this case, the infected tumor population $y = \frac{rk}{r+bk}$, and the uninfected tumor population $x = 0$, see Figure 9. Parameters in Figure 2-Figure 9 are chosen as listed in Table 4.

3.1.   Virus therapy

In what follows, we are interested in the total tumor load $x+y$. In the simplest case, $R_0\leq 1$, the original tumor equilibrium $E_1 = (k, 0, 0)$ is globally stable. In this case, total tumor load is equal to $x_1+y_1 = k$. This is the original total tumor load. Next we focus on virus therapy equilibrium $E_{vt}$. It is easy to see that $E_{vt}$ exists if and only if $R_0 > 1$, i.e.

By inequality (3.1), we see that in order to make the viral therapy work, one needs to engineer the oncolytic virus with suitable properties. Among these are weaker cytotoxicity $a$, faster viral replication rate $b$, and/or shorter period of intracellular cycle $\tau$. By the definition of $R_0$, one finds

It is obvious that $x_{vt}+y_{vt} < k$ when $R_0 > 1$, where $k$ is the original total tumor load. Besides, uninfected tumor cells $x_{vt}$ will be decreasing under the conditions of weaker cytotoxicity $a$, faster viral replicating rate $b$, or shorter period of intracellular cycle $\tau$. Intuitively, these condition will benefit the infected tumor cells $y_{vt}$, conforming to biological intuition. The short period of intracellular cycle means the high effect of transformation from the uninfected tumor cells to infected tumor cells. And the weak cytotoxicity reduces the death rate of infected tumor cells. All of these conditions will benefit the growth of infected tumor cells, which leads to increased contact probability between infected tumor cells and uninfected tumor cells, and thus decreasing the uninfected tumor population. If the cytotoxicity of virus $a = 0$, then the infected tumor cells may invade the uninfected tumor cells completely, i.e. $x = 0$.

Now we consider the total tumor load $x_{vt}+y_{vt}$ in different situations.

Case 1. When the period of intracellular cycle $\tau = 0$ and $R_0 > 1$, then the viral therapy equilibrium is stable. In this case, the total tumor load $x_{vt}+y_{vt} = \frac{kr+ak}{r+bk}$. Both weak cytotoxicity and fast viral replication lead to reduced total tumor load. When the viral replication rate $b$ is fixed, the optimal cytotoxicity rate is $a = 0$, and the total tumor load is $\frac{kr}{r+bk}$. If $a$ is fixed, then the total tumor load tends to $0$ as $b$ tends to infinity, see Figure 10.

Case 2. When viral cytotoxicity is $a = 0$ and $R_0 > 1$, the $100\%$-virus prevalence equilibrium $E_*$ is stable. In this case, the uninfected tumor cell density is zero and the total tumor load is $\frac{kr}{r+bk}$, see Figure 9. Hence, increasing viral replication rate $b$ will reduce the total tumor load.

Case 3. If $\tau\neq0$, $a\neq 0$ and $1 < R_0\le 3$, then the viral therapy equilibrium is stable, see Figure 6. Recall that the total tumor load is given by

When $x_{vt}+y_{vt}$ is considered as a function of $R_0$ and $b$, it can be rewritten as

This is decreasing function of $R_0$. Thus, $\min_{R_0\in(1, 3]} x_{vt}+y_{vt} = k\frac{3r+bk}{3r+3bk}$. And $k\frac{3r+bk}{3r+3bk}$ is also a decreasing function of $b$. Thus, the minimal total tumor load is $\frac{k}{3}$ as $R_0 = 3$ and $b$ tends to $\infty$.

Case 4. When the period of intracellular cycle $\tau\neq0$, viral cytotoxicity $a\neq0$, and $R_0 > 3$, infected tumor load and uninfected tumor load oscillate about a positive mean value, see Figure 7. In fact, a Hopf Bifurcation may occur from the viral therapy equilibrium $(x_{vt}, y_{vt}, 0)$. By Eq. (3.2), rapid replication $b$, weak cytotoxicity of virus $a$ and short cycle time of the intracellular $\tau$ facilitate tumor remission, if amplitude of this periodic solution is small. Large amplitude oscillations may be dangerous since it may mislead patients to stop further therapy when the total load is of very low. However, suitable virus-specific CTL response will make this phenomenon fade away, see Figure 8. Furthermore, we conjecture that the combination of viral therapy with immunotherapy should be an effective strategy for tumor remission.

3.2.   Virus-specific CTL response

When $R_0 > 1$ and $R_1 = \frac{crk}{d(r+bk)}(1-\frac{1}{R_0}) > 1$, virus specific CTL response occurs. In this case, uninfected tumor population density is given by

and infected tumor population density is $y_{CTL} = \frac{d}{c}$ so total tumor load equals

Here, $\frac{c}{d}$ can be viewed as an index representing the strength of virus-specific immune response. Larger $\frac{c}{d}$ corresponds to stronger virus-specific immune response and note that $R_1 > 1$ implies a lower bound for it. Consequently, if virus CTL response is strong, the uninfected tumor population $x_{CTL}\approx k$, infected tumor cells $y_{CTL}\approx 0$ and total tumor population $x_{CTL}+y_{CTL}\approx k$. This implies that the oncolytic virotherapy fails due to the interference by the immune system as reported in ^[22]. On the other hand, by the definition of $R_1$, $\frac{c}{d} > \frac{r+bk}{rk}(1-\frac{1}{R_0})^{-1} > \frac{b}{r}$ holds. Therefore, when $\frac{c}{d}$ tends to $\frac{r+bk}{rk}(1-\frac{1}{R_0})^{-1}$, total tumor load will tend to its minimum value $\frac{k(r+a\exp\{n\tau\})}{r+aR_0 \exp\{n\tau\}}$, see Figure 11.

As references ^[1,22,23] showed that immune suppressive drug can decrease the virus-specific CTL response. We assume $\frac{c}{d}$ can be decreased by using immune suppressive drugs. In ^[23], the authors used the percentage of viable (uninfected) cells, given by $\frac{x_{CTL}}{x_{CTL}+y_{CTL}}$ to quantify the oncolytic effect of virotherapy. By the expressions for $x_{CTL}$ and $y_{CTL}$, we see that $y_{CTL}$ will be increasing and $x_{CTL}$ be decreasing on using immune suppressive drugs (decreasing $\frac{c}{d}$), which leads to a decrease of the percentage of uninfected cells, see Figure 12. Thus, our model corroborates the experimental results obtained in ^[23].

4.   Mathematical proofs

In this section, we give detailed proofs of our mathematical results listed in Table 2 and Table 3. Due to the biological background, we study Eq. (2.1) with nonnegative initial conditions

where $\varphi_i(\theta)\ge 0$, $\varphi_i(0) > 0$ for $\theta\in [-\tau, 0], i = 1, 2, 3$.

Let $C = C([-\tau, 0], {\bf R}^{3})$ be the Banach space of continuous mappings from $[-\tau, 0]$ into ${\bf R}^{3}$ with norm

We denote

As usual, for any continuous function $x\in C([-\tau, +\infty), {\bf R})$ and any given $t\ge 0$, $x_t\in C([-\tau, 0], {\bf R})$ is defined as $x_t(\theta) = x(t+\theta)$, for any $\theta\in [-\tau, 0]$.

Theorem 4.1. For any $(\varphi_1(\theta), \varphi_2(\theta), \varphi_3(\theta))\in C^+$, there exists unique solution to Eq. (2.1) which is nonnegative. Moreover, there exist positive constants $M_1, M_{2}, M_3$ independent of initial data, such that $x(t)\leq M_1, y(t)\leq M_{2}, z(t)\leq M_3$ for sufficiently large $t$.

Proof. By Theorem 3.4 in ^[27], the first part of Theorem 4.1 holds. By the first equation of Eq. (2.1),

So

For any fixed $\varepsilon > 0$, let $M_1 = k+\varepsilon$, then $x(t)\leq M_1$ for sufficiently large $t$.

Let $u(t) = b\exp\{-n\tau\}x(t)+(b+\frac{r}{k})y(t+\tau)$, then

Thus, $b\exp\{-n\tau\}x(t)+(b+\frac{r}{k})y(t+\tau)\leq br\exp\{-n\tau\}\frac{brk\exp\{-n\tau\}}{4 a}$. By $x(t), \ y(t)\geq 0$, there exists a positive constant $M_2 > 0$ such that $y(t)\leq M_2$ for sufficiently large $t$. Similarly, consider $v(t) = y(t)+\frac{p}{c}z(t)$, we can find that $v$ is bounded. Therefore, there is a constant $M_3$ such that $z(t)\leq M_3$ for $t$ large enough.

Now we consider the linearization of the Eq. (2.1) at equilibrium $(x_*, y_*, z_*)$ as follows.

where $X(t) = (x(t), y(t), z(t))^{T}$,

The characteristic equation of Eq. (4.1) is formulated as

Theorem 4.2. Trivial equilibrium $E_0 = (0, 0, 0)$ is always unstable and equilibrium $E_1 = (k, 0, 0)$ is asymptotically stable when $R_0 < 1$.

Proof. Let $x_* = y_* = z_* = 0$ in Eq. (4.2). We have $\lambda_1 = r, \lambda_2 = -a, \lambda_3 = -d$. So trivial equilibrium $E_0 = (0, 0, 0)$ is always unstable.

Let $x_* = k, y_* = z_* = 0$ in Eq. (4.2). Then

It is clear that $\lambda_1 = -r, \lambda_2 = -d$. The roots of equation

dominate the stability of $E_1$. As $R_0 < 1$, it is easy to check that $bk\exp\{-n\tau\}-a < 0$ and $b\exp\{-n\tau\}k > -a$. By Theorem 4.7(b) in ^[27], the original tumor load equilibrium $E_1 = (k, 0, 0)$ is asymptotically stable.

Theorem 4.3. If $R_0 > 1$ and $R_1 < 1$, then equilibrium $E_{vt}$ is locally stable when $\tau = 0$. If $\tau > 0$ and if either $1 < R_0\le 3$ or $a = 0$, then virus therapy equilibrium $E_{vt}$ is also locally stable. When $R_0 > 3$, equilibrium $E_{vt}$ undergoes a Hopf bifurcation at a critical value of the delay $\tau$.

Proof. By substituting $E_{vt}$ in characteristic Eq. (4.2), we have

It is easy to check that $\lambda = cy_{vt}-d < 0$ due to $R_1 < 1$. Thus the roots of the following equation

determine the stability of equilibrium $E_{vt}$, where

Consequently, $a_1, a_2 > 0$ as $\tau = 0$ which implies that the characteristic Eq. (4.2) has no solution with positive real part.

For $\tau > 0$, we rewrite Eq. (4.3) as below.

where

Suppose that there exists a $\tau_0 > 0$ such that Eq. (4.4) has a pair of pure imaginary roots $\pm i \omega, \omega > 0$. Then $\omega$ satisfies

Separating the real and imaginary parts, we have

Then we have

By direct calculation, we have

Substituting $x_{vt}$ and $y_{vt}$ into formula (4.9), then

Therefore, if $1 < R_0\leq 3$ or $a = 0$, there is no positive positive solutions to the Eq. (4.7). Thus, the solutions of Eq. (4.4) always have negative real parts, which implies that $E_{vt}$ is asymptotically stable.

Otherwise, if $R_0 > 3$, then there exists some $\omega^2 > 0$ satisfying Eq. (4.7). In this case, by Eq. (4.5) and Eq. (4.6), we have

where $k = 1, 2, 3, \cdots$.

Suppose that $\lambda(\tau) = \sigma(\tau)+i\omega(\tau)$ is a root of Eq. (4.7) with $\tau > \tau_0$. Differentiating Eq. (4.4) with respect to $\tau$, we obtain

So we have

By formula (4.8), we have $sign \Re \frac{d\tau}{d\lambda}|_{\tau = \tau_0} > 0$, where $\Re \frac{d\tau}{d\lambda}|_{\tau = \tau_0}$ means the real part of $\frac{d\tau}{d\lambda}|_{\tau = \tau_0}$. By Hopf bifurcation theorem (Theorem 6.1 in Chapter 6 of ^[27]), conclusions in Theorem 4.3 hold.

Theorem 4.4. If $R_0 > 1$ and $R_1 > 1$, then the CTL equilibrium $E_{CTL}$ exists and is stable for small values of the delay $\tau$.

Proof. In order to find the CTL equilibrium, we consider algebraic system

Assuming $z_{CTL}\neq0$, we get $y_{CTL} = \frac{d}{c}$ from the third equation of equation (4.11). Substituting $y_{CTL}$ into the first equation of system (4.11), we get $x_{CTL} = \frac{ckr-(r+bk)d}{cr}$. By the second equation of system (4.11), we get $z_{CTL} = \frac{bx_{CTL}\exp\{-n\tau\}-a}{p}$. In order to find positive solutions, we assume $x_{CTL}, y_{CTL}, z_{CTL} > 0$. Obviously, $y_{CTL} = \frac{d}{c} > 0$. It is easy to see that $z_{CTL} = \frac{bx_{CTL}\exp\{-n\tau\}-a}{p} > 0$ can lead to $x_{CTL} > 0$. Assume $z_{CTL} > 0$, we have

Substituting $E_{CTL}$ into determinant (4.2) and by direct calculation, we get

It is easy to see that if $\tau = 0$, then $b_1, b_2, b_3 > 0$ and $b_1b_2 > b_3$. Thus, by Hurwitz criterion, Eq. (4.12) has no solution with positive real part. Because solutions of Eq. (4.12) depend continuously on its parameters, there is no root with positive real part of Eq. (4.12) for sufficiently small $\tau$.

Theorem 4.5. If $R_0 < 1$, then the original tumor load equilibrium $E_1$ is globally stable in the first quadrant for every $\tau\geq0$.

Proof. Let

where $V_{11}(\varphi_1, \varphi_2, \varphi_3) = \varphi_1(0)-k-k\ln\frac{\varphi_1(0)}{k}$, $V_{12}(\varphi_1, \varphi_2, \varphi_3) = L\varphi_2(0)$, $V_{13}(\varphi_1, \varphi_2, \varphi_3) = \frac{pL}{c}\varphi_3(0)$, $V_{14}(\varphi_1, \varphi_2, \varphi_3) = \frac{r+bk}{k}\int_{-\tau}^{0}\varphi_1(\theta)\varphi_2(\theta)d\theta$ and $L = \frac{r+bk}{bk\exp\{-n\tau\}}$. Then,

Clearly, $V_1$ is always positive except for $x = k, y = 0, z = 0$.

Differentiating $V_{1}$ with time $t$ along Eq. (2.1), we get

Then

As $R_0 < 1$, we have $V_1'(x_t, y_t, z_t)\leq 0$. If $V_1'(x_t, y_t, z_t) = 0$, then $x = k, y = 0, z = 0$. Hence, by LaSalle's invariance principle, $E_1$ is globally stable.

Theorem 4.6. If $R_0 > 1, R_1 < 1, \tau = 0$, then equilibrium $E_{vt}$ is globally stable in the first quadrant.

Proof. Let

where,

Differentiating $V_{2i}, i = 1, 2, 3$ along Eq. (2.1),

Thus, we have that

Since $R_1 = \frac{c}{d}y_{vt} < 1$, $V_2'(x, y, z)\le 0$ and $V_2'(x, y, z) = 0$ if and only if $x = x_{vt}, z = 0$. Set

By Lemma A.16 in ^[28], we obtain $\frac{dx}{dt} = 0$ which leads to $y = y_{vt}$. Thus, the largest invariant set of (2.1) contained in $D_{00}$ is $E_{vt}$. By LaSalle's invariance principle, $E_{vt}$ is globally stable.

Theorem 4.7. If $R_0 > 1, R_1 > 1, \tau = 0$, then the CTL equilibrium $E_{CTL}$ is globally stable in the first quadrant.

Proof. Let

where,

By similar calculations, we have

Consequently,

and $V_3'(x, y, z) = 0$ if and only if $x = x_{CTL}$. Set

By Lemma A.16 in ^[28], we can see $\frac{dx}{dt} = 0$ which leads to $y = y_{CTL}$. Using the same method in the second equation of Eq. (2.1), we further obtain $z = z_{CTL}$. Thus, the largest invariant set of (2.1) contained in $D_{01}$ is $E_{CTL}$. Conclusions of Theorem 4.7 follows from LaSalle's invariance principle.

5.   Discussion and conclusion

Since oncolytic virotherapy has significant potential benefits in the fight against cancer, there has been much interest in constructing and analyzing mathematical models of the effects of virotherapy. To the best of our knowledge, Wodarz was the first to model oncolytic virotherapy using a simple ODE system and the first PDE model was considered by Wu ^[29]. Neither of these authors included an intracellular viral life-cycle delay in their models. The first to do so was the paper of Wang ^[14], although their model did not include an immune response. Thus, we formulate, to the best of our knowledge, the first time, an oncolytic virotherapy model with both virus-specific CTL response and viral life-cycle delay.

Our results provide new insight into how to engineer an effective oncolytic virus. Furthermore, we show that immune suppressive drugs can enhance the oncolytic effect (i.e. the percentage of viable tumor cells $\frac{x}{x+y}$ will decrease) which coincides with the experimental results in ^[23].

We choose several strategies to reduce the total tumor load in different situations, see Table 5. In addition, we also verify that immune suppressive drugs can enhance the oncolytic effect of virus.

The most favorable result would be that the oncolytic virus should be designed to have low toxicity, short viral life-cycle, fast replication, and virus-specific immune avoidance.

Generally, $R_0$ will determine the total tumor load as $\tau\neq0, a\neq0, R_1 < 1$. When $R_0\in(1, 3]$, the optimal choice of parameters satisfies $bk = 3a\exp\{n\tau\}$ and $b$ is big enough as $\lim\limits_{b\rightarrow\infty} x+y = \frac{k}{3}$. Thus, virus should be designed: rapid replication rate $b$ and cytotoxicity $a$ should be close to number $\frac{bk}{3\exp\{n\tau\}}$. On the other hand, if $R_0\in(3, +\infty)$, the tumor load changes periodically. When the amplitude of this periodic solution is small, virus needs to be designed as rapid replication, weak cytotoxicity of virus, and short intracellular viral life-cycle. Otherwise, suitable virus-specific immune responses can reduce the periodic oscillations ultimately causing the tumor load to tend to a constant, see Figure 13.

It is worth mentioning that dosage of medication is not considered in the present paper, one can find the successful therapy treatment in our recent papers ^[30,31] where the optimal dosage of medication is given and the total tumor load tends to $0$ under suitable conditions.

In this paper, we focus on how to design a more effective virus to be used in oncolytic virotherapy. Original tumor load equilibrium $E_1$ is globally stable as $R_0\leq 1$. Virus therapy equilibrium $E_{vt}$ is stable as $R_0 > 1, \ R_1\leq 1, \ \tau = 0$ or $R_0\in(1, 3], \ R_1\leq 1 \ \tau\neq 0$. Besides, virus therapy equilibrium can undergo a Hopf bifurcation such that the cell and virus populations oscillate in time as $R_0 > 3, \ R_1\leq 1, \ \tau\neq0$. Lastly, virus-specific CTL equilibrium $E_{CTL}$ is stable for small $\tau$, $R_0 > 1$ and $R_1 > 1$.

Acknowledgments

This work was supported by National Natural Science Foundation of China (No. 11371107, 11771104), Program for Changjiang Scholars and Innovative Research Team in University (IRT-16R16) and Simons Foundation Grant 355819.

We are very grateful to anonymous referees for their valuable comments and suggestions. The first author would like to thank the staff in the School of Mathematical and Statistical Sciences at Arizona State University for their friendly and helpful assistance during his visit. He would especially like to thank Professor Yang Kuang for his help.

Conflict of interest

The authors declared that they have no conflicts of interest to this work.