A mathematical model for tumor growth and treatment using virotherapy

1.   Introduction

Oncolytic virotherapy is a cancer treatment that involves injecting cancerous tumor cells with a virus to both infect and break down those cells without destroying healthy cells ^[1]. This treatment also works to jump-start the body’s defenses by stimulating the immune system ^[2]. Oncolytic virotherapy attacks cancer as a virus would normally attack the body and works without chemotherapy agents or any kind of radiation. Because of this, it is not vulnerable to the same drug and radiation resistance that current commonly used treatments experience. Due to the way in which virotherapy works, this type of treatment can be applied in tandem with other treatments; it can be administered before or after surgery or between radiation or chemotherapy appointments ^[2]. The average length of virotherapy treatment is three years with scheduled monitoring, and oncolytic virotherapy avoids the detrimental side effects other common cancer treatments such as chemotherapy and radiation tend to exhibit ^[3]. In order for a virus to be acceptable for oncolytic virotherapy, it must be capable of replication and selective infection ^[4].

In recent years, scientists have been looking into the possibility of a single-shot cure, a threshold limit for vascular delivery, and an alternate way for viruses to target the cancer cells ^[5]. In 2017, information was released on the first study of herpes simplex virus-1, HSV-1, being used in children and young adults for its oncolytic properties ^[6]. Researchers from both Nationwide Children’s Hospital and Cincinnati Children’s Hospital Medical Center have together completed the first phase 1 trial of a mutated HSV-1 virus, HSV1716 ^[6]. They have determined that after the completion of phase 1, the HSV1716 is both endurable and nontoxic ^[7]. As of 2018, viruses from nine families have progressed to clinical trials of oncolysis. While these viruses have shown encouraging results, their efficacy as a single agent therapy is limited ^[8]. Current research is exploring how oncolytic viruses can support immunotherapy, particularly in cancers susceptible to checkpoint inhibitors ^[9].

Since the onset of oncolytic virotherapy, mathematicians have utilized experimental results and analytical techniques to build mathematical models which could be studied to determine key model parameters, as well as short and long-term dynamics of such a treatment approach. Previous mathematical studies have incorporated a virotherapy treatment within their models using a constant source rate; these same papers have focused on the dynamics of the infected and uninfected cell populations in their main equations, without a free virus equation ^[10,11]. Some studies have included an immune response in their system of differential equations ^[11,12], while others have neglected to include an immune response to the cancerous cells ^{[13,14,15,16]}. The aim of this paper is to study the long-term dynamics of a system of nonlinear differential equations that describes the role of virotherapy on tumors and the impact of immune response specific to fighting cancer.

2.   Model development

In 2015, Kim, Crivelli and Choi et al. ^[17] studied the short-term dynamics of a model describing the interaction between an oncolytic virus and tumor cells. Their model explicitly utilizes a free virus population that includes a virotherapy treatment term as well as growth due to infected tumor cell lysis:

where $U, I, V, T,$ and $A$ represent uninfected tumor cells, infected tumor cells, virions, T-cells, and APCs ^[17].

In ^[17], Kim et al. utilize experimental data to fit parameter values to their proposed model, and then vary treatment strategies to determine the effects of various dosage, treatment schedules, and targeted viruses have on the short-term behavior of the tumor cell populations. The authors conclude that the most important factors in controlling short-term tumor growth are the immune response and the virus burst size.

The goal of this paper is to use build upon the work of Kim et al. ^[17] to study the long-term dynamics of oncolytic virotherapy on a tumor cell population. The use of an exponential growth rate in ^[17] allows for the model to simulate up to at most sixty days. Since we are interested in the long-term dynamics of virotherapy, we propose the following modification to the model presented in ^[17]. First, we include logistic growth in place of exponential growth to allow for long-term population dynamics. Secondly, we simplify the immune response to study the response of effector T-cells on the infected tumor cell population. Here, we assume that the effector T-cells are recruited at a rate proportional to the infected tumor cell population and the effector T-cells decrease the infected tumor cell population according to the law of mass action. The effector T-cell recruitment is consistent with the assumptions made in Kim et al. ^[17], while relaxing the frequency-dependent impact on tumor cells facilitates the global study of long-term dynamics. Furthermore, these assumptions are well-documented in mathematical models, and we refer the reader to ^[18] for further reading on the modified terms.

Thus, we propose the following model describing the interactions between uninfected tumor cells, $U$, infected tumor cells, $I$, effector T-cells, $E$, and virions, $V$.

Here, we use $r$ to represent the growth rate of uninfected tumor cells, $k$ to represent the total carrying capacity of tumor cells, $\beta$ to represent the rate of uninfected tumor cells becoming infected, $\gamma_U$ to represent the rate of decay of uninfected cells via T-cells, $\gamma_I$ to represent the rate of decay of infected cells via T-cells, $\delta_I$ to represent the rate of decay for infected cells, $c$ to represent the rate of T-cell growth via infected tumor cells, $\delta_E$ to represent the rate of decay for effector T-cells, $\delta_V$ to represent the rate of decay for virions, $\alpha$ to represent the number of virions released via infected cell lysis, and $N$ to represent the virotherapy dosage.

To simplify later calculations, we non-dimensionalize our model by setting $\tilde{t} = rt, \tilde{U} = \frac{U}{k}, \tilde{I} = \frac{I}{k}, \tilde{E} = \frac{\gamma_U}{r} E, \tilde{V} = \frac{\beta}{r} V, \tilde{\gamma} = \frac{\gamma_I}{\gamma_U}, \tilde{\delta_I} = \frac{\delta_I}{r}, \tilde{c} = \frac{\gamma_u k c}{r^2}, \tilde{\delta_E} = \frac{\delta_E}{r}, \tilde{N} = \frac{\beta N}{r^2}, \tilde{\alpha} = \frac{\beta k \alpha}{r},$ and $\tilde{\delta_V} = \frac{\delta_V}{r}.$ We drop all the tildes from our notation and arrive at the following non-dimensionalized model:

3.   Analysis

3.1.   Preliminary results

Assuming (2.1) is subject to non-negative initial conditions, we note that the system is invariant in the non-negative orthant. Additionally, because the associated vector field is continuously differentiable, there exists a unique solution to (2.1) under non-negative initial conditions by the Picard–Lindelöf Theorem.

3.1.1.   Boundedness

In order to confirm that our model does not predict uninhibited cell growth, we ensure that our cell populations are bounded above. For uninfected tumor cells:

Thus, $\limsup\limits_{t \rightarrow \infty}U(t) \leq 1$. Using this upper bound, we can derive an upper bound for the infected tumor cell population. From (2.1), we have:

It follows that $\limsup\limits_{t \rightarrow \infty}I(t) \leq 1$. Utilizing the asymptotic upper bound of the infected tumor cells, for effector T-cells we have $\frac{dE}{dt} \leq c - \delta_E E.$ By standard comparison theory, it follows that $\limsup\limits_{t \rightarrow \infty}E(t) \leq \frac{c}{\delta_E}$. Similarly, for the virion population, we have $\limsup\limits_{t \rightarrow \infty}V(t) \leq \frac{N + \alpha \delta_I}{\delta_V}.$

3.1.2.   Existence of equilibria

To establish equilibria of (2.1), we must solve the following system of equations:

If $U^{*} = 0$, it readily follows that $I^{*} = E^{*} = 0$ and $V^{*} = \frac{N}{\delta_V}$. Thus we find a cure state equilibrium point of the form $P_0 = (0, 0, 0, \frac{N}{\delta_V})$.

If $U^{*}\neq 0$, we can use (3.1), (3.3), and (3.4) to solve for $U^{*}, E^{*},$ and $V^{*}$ in terms of $I^{*}$:

Substituting these expressions into (3.2) leaves us with a polynomial in $I^{*}$, denoted by $f(I^{*})$:

The number of internal equilibria is thus determined by the number of solutions to $f(I^{*}) = 0$. We first note that $f(I^{*})$ is a quadratic function and that the coefficient on $I^{*^2}$ is negative. We also note that the constant term is positive if and only if $\delta_V > N$. By Descartes' Rule of Signs, it follows that there exists one unique positive real root for $f(I^{*})$. Since $I^{*}$ is positive and real, $U^{*}, E^{*},$ and $V^{*}$ must also be positive and real. We conclude that there exists a unique cancer persistence state of the form $P^* = (U^{*}, I^{*}, E^{*}, V^{*})$ when $\delta_V > N.$ We summarize these results in the following theorem:

Theorem 1.

1. There exists a unique cure state of the form $P_0 = (0, 0, 0, \frac{N}{\delta_V})$.

2. When $N < \delta_V$, there exists a unique cancer persistence state of the form $P^* = (U^{*}, I^{*}, E^{*}, V^{*})$.

3.2.   Stability of the cure state

In this section, we explore the stability of the cure state equilibrium $P_0 = \left(0, 0, 0, \frac{N}{\delta_V}\right).$ We note that the nonzero virion population in the cure state results from assuming a continuous constant dosage treatment. Furthermore, the lack of effector T-cells in the cure state represents there no longer being a need for an immune response due to cancer clearance.

3.2.1.   Local stability of the cure state

We first consider the local stability of the cure state equilibrium $P_0.$ Recall that our non-dimensionalized model is

Evaluating the Jacobian matrix at $P_0$ yields

with eigenvalues $\dfrac{ - N + \delta_V }{\delta_V}, - \delta_I, - \delta_V, -\delta_E$. Thus, we find the local stability condition for the cure state to be $N > \delta_V$. If this condition is not met, the cure state is unstable. We demonstrate this result with numerical simulations in Section 4.

3.2.2.   Global stability of the cure state

We next explore the global stability of $P_0$. We begin by establishing a lower bound on the virion population $V$:

Using standard comparison theory, this implies that $\mathop {\lim \inf }\limits_{t \to \infty } V(t) \ge \frac{N}{{{\delta _V}}}$.

We next seek conditions to ensure $U(t) \rightarrow 0$ as $t \rightarrow \infty$. Defining the Lyapunov-type function $L = U$, we compute its derivative along trajectories of our system:

Using the lower bound on $V$, we then have

Thus $\frac{dL}{dt} < 0$ if $1 - \frac{N}{\delta_V} < 0$, or equivalently if $N > \delta_V$. We note that this is the same condition for local stability of the cure state. Under this condition, it then follows that $U(t) \rightarrow 0$ as $t \rightarrow \infty$.

We next explore the asymptotic behavior of $I$:

Thus $U(t) \rightarrow 0$ implies $I(t) \rightarrow 0$ as well. Similarly, for $E$ and $V$ we have:

Hence $I(t) \rightarrow 0$ implies $E(t) \rightarrow 0$ and $V(t) \rightarrow \frac{N}{\delta_V}$. These results are summarized as follows:

Theorem 2. If $N > \delta_V$, then the cure state $P_0 = \left(0, 0, 0, \frac{N}{\delta_V}\right)$ is globally asymptotically stable.

3.3.   Stability of the cancer persistence state

We now explore the long-term dynamics of the model when the stability condition for the cure state is violated; that is, $N < \delta_V$. Recall that in this case, the cancer persistence equilibrium $(U^*, I^*, E^*, V^*)$ exists in the nonnegative orthant. Although full stability analysis of this equilibrium proves intractable, we can glean some useful bounds on our treatment term $N$ from a local stability analysis. Substituting the cancer persistence equilibrium into the Jacobian matrix yields:

After algebraic simplification, the characteristic polynomial can be written in the form:

where $a_i > 0$ for $i = 0, 1, 2, 3$. Our characteristic polynomial has the following coefficients:

The Routh-Hurwitz criterion provide necessary and sufficient conditions for the cancer state to be locally stable, normally with the conditions

While condition (3.5) is satisfied, establishing (3.6) and (3.7) in general is intractable due to the nonlinear nature of the model. However, we note that (3.7) implies (3.6), and we will show in the numerical example below that (3.7) establishes a lower bound for stability of cancer persistence on the dosage, $N$.

4.   Numerical results

We continue our exploration of the stability of the cancer persistence state numerically, with parameter values chosen from literature as seen below in Table 1. We utilize the parameter values derived from experimental results in ^[17] for our numerical simulations, in which the authors perform various simulations that elicit different immune responses. Since we are only considering the effector T cell immune response in our model, we choose parameter values from ^[17] that are present with no APC response. To address the change from exponential growth to logistic growth, we use a carrying capacity of $k = 3 \times 10^9.$ This is equivalent to a tumor having a volume of 3000 $mm^3,$ which Kim et al. ^[17] found experimentally to be the lethal size of tumors in mice. Finally, in ^[17], the authors divide the infection rate $\beta = 8.9 \times 10^{-4}$ by the total cell population; thus we similarly adjust the value of $\beta$ to align with the chosen carrying capacity. The corresponding non-dimensionalized parameter values are found in Table 2.

With the non-dimensionalized parameters from Table 2, Theorem 2 guarantees a globally stable cure state for $\tilde{N} > 7.41935$ (corresponding to $N > 8.01123 \times 10^{11}$ virions), while the Routh-Hurwitz criteria above predicts a locally stable cancer persistence state for $0.000819195 < \tilde{N} < 7.41935$ (corresponding to $8.84546 \times 10^7 < N < 8.01123 \times 10^{11}$ virions). The expected behavior of the model with representative dosages chosen from these ranges can be seen in Figures 1 and 2.

For representative dosages satisfying $\tilde{N} < 0.000819195$ (corresponding to $N < 8.84546 \times 10^7$ virions), we observe sustained oscillations in each population as seen in Figure 3.

This suggests the existence of a Hopf bifurcation at the critical value $\tilde{N} = 0.000819195,$ (correspondingly $N = 8.84546 \times 10^7$ virions.) The following theorem from ^[19] provides a method of establishing the existence of a Hopf bifurcation and analyzing its stability by calculating the first Lyapunov coefficient:

Theorem 3. Suppose that the system $\dot{x} = f(x, \mu), x \in \mathbb{R}^N, \mu \in \mathbb{R}$ has an equilibrium $(x_0, \mu_0)$ and the following properties are satisfied:

(H1) The Jacobian $D_x f \Big|_{(x_0, \mu_0)}$ has a simple pair of pure imaginary eigenvalues $\lambda(\mu_0)$ and $\overline{\lambda(\mu_0)}$ and no other eigenvalues with zero real parts,

(H2) $\frac{d}{d\mu} \mathit{\mbox{Re}}(\lambda(\mu)) \Big|_{\mu = \mu_0} \neq 0.$

Then the dynamics undergo a Hopf bifurcation at $(x_0, \mu_0)$ resulting in periodic solutions. The stability of the periodic solutions is given by the sign of the first Lyapunov coefficient of the dynamics $l_1(x_0, \mu_0).$ If $l_1 < 0,$ then these solutions are stable limit cycles and the Hopf bifurcation is supercritical, while if $l_1 > 0,$ the periodic solutions are repelling.

To numerically calculate the first Lyapunov coefficent $l_1(P^*),$ we utilize the following theorem from ^[20]:

Theorem 4. Let $\frac{dx}{dt} = F(x)$ be a differential system having $P^*$ as an equilibrium point. Consider the third order Taylor approximation of $F$ around $P^*$ given by

Assume that $A$ has a pair of purely imaginary eigenvalues $\pm \omega_0 i$. Let $q$ be the eigenvector of $A$ corresponding to the eigenvalue $\omega_0 i$. Let $p$ be the adjoint eigenvector such that $A^T p = - \omega_0 i p$ and $\langle \overline{p}, q \rangle = 1$. If $I$ denotes the identity matrix, then the first Lyapunov coefficient $l_1 (P^*)$ of the system $\frac{dx}{dt} = F(x)$ at $P^*$ is given by

where

We begin numerically establishing the existence of a Hopf bifurcation by recalling the Jacobian evaluated at $P^* = (U^*, I^*, E^*, V^*)$:

Using the non-dimensionalized parameter values found in Table 2 together with the critical dosage value $\tilde{N} = 0.000819195$, we first calculate the values of $U^{*}, I^{*}, E^{*},$ and $V^{*}$. We find the approximate values of ($0.265432, 0.0000609925, 0.733598, 0.000909817$) and substitute them into the Jacobian matrix to find

This gives the following eigenvalues and corresponding eigenvectors $(\lambda_i, v_i)$ of the system:

With these eigenvalues and eigenvectors, we note that (H1) of Theorem 3 is satisfied. To establish the transversality condition (H2), we use the following result from ^[20]:

with $A$, $p$ and $q$ as defined in Theorem 4. Using the eigenvalues and eigenvectors calculated above, we compute $\omega_0 = 1.188783$ and $q = \begin{bmatrix}-0.0466716 + 0.208015 i \\ 0.0000812302 + 0.0000855293 i\\0.977010 \\ 0.00121312 + 0.000926622 i \end{bmatrix}$. To calculate $p,$ we derive the eigenvalues and corresponding eigenvectors of $A^T:$

We then find that $p = \begin{bmatrix}-0.000162668 - 0.000728537 i\\-0.999377\\-0.0000417853 + 0.00012728 i\\-0.0348568 - 0.00555895 i\end{bmatrix}$.

When we compute $\overline{p} \cdot q$, we see that it is not equal to $1$, but is equal to $-0.000313395 - 0.000303225i$ so we must scale $q$ as follows:

Finally, we compute

and so

Thus, we have satisfied the transversality condition of Theorem 3 and therefore have established the existence of a Hopf bifurcation.

To determine the stability of the Hopf bifurcation, we proceed with calculations needed for $l_1 (P^*)$:

$(2i \omega_0 I_4 - A)^{-1} =$

$\left[\begin{array}{cccc} -0.0042 - 0.4720 i & 133.179 - 148.759 i & 0.04497 + 0.0243 i & 2.7793 - 6.1957 i \\ -0.00012 - 0.00004 i & 0.0602 - 0.3429 i & 0.00002 - 5.2976 \times 10^{-6} i & -0.0016 - 0.0118 i\\ -0.4547 + 0.5099 i & -1465.14 - 1039.52 i & 0.1725 - 0.4270 i & -58.333 - 18.5147 i \\ -0.0017 + 0.00004 i & -0.5909 - 4.3056 i & 0.0002 - 0.0001 i & 0.0584 - 0.1727 i \end{array} \right],$

Substituting the above matrices and vectors into $l_1 (P^*)$ yields

Since $l_1 (P^*) < 0$, we conclude that the Hopf bifurcation that occurs at the critical dosage $N = 8.84546 \times 10^7$ virions is supercritical, implying that the cancer persistence equilibrium bifurcates into a stable limit cycle. This agrees with the behavior observed in Figure 3.

5.   Conclusions

The model of Kim et al. ^[17] uses experimental data to simulate tumor progression for up to sixty days due with the use of an exponential growth term for the uninfected tumor cells. In this work, we modify the model presented in ^[17] to allow for long-term dynamics through logistic growth. In doing so, we find thresholds for dosages of an effective virotherapy treatment that are sufficient to reach long term tumor eradication. Conversely, when the virotherapy protocol is not strong enough to ensure tumor eradication, our model gives two possibilities. The first possibility is a stable cancer persistence state where the tumor may shrink, but is never eradicated. In such a case, this model predicts that virotherapy could be useful as a neoadjuvant therapy in preparation for surgery or radiotherapy treatment. The second possibility is periodic cancer recurrence that may indicate further progression of the tumor or metastasis.

In this latter case, to numerically demonstrate the existence of a Hopf bifurcation, we compute the first Lyapunov coefficient using parameter values gathered from the literature. $l_1 (P^*)$ was found to be less than zero, which allows us to conclude that a supercritical Hopf bifurcation exists at the critical dosage and thus the resulting limit cycle is stable. This calculation, combined with the Routh–Hurwitz criteria, motivates the following conjecture for bounds on the specific range of dosage to guarantee a stable cancer persistence state:

Let $N_0$ denote the smallest positive root of $a_3 a_2 a_1 = a_1^2 + a_3^2 a_0$. If $N_0 < N < \delta_V$, then the cancer persistence state $(U^{*}, I^{*}, E^{*}, V^{*})$ is globally asymptotically stable. If $N < N_0$, there exists a globally stable limit cycle.

While Kim et al. ^[17] found that the most important factors in controlling short term tumor growth were the immune response and the virus burst size, our model suggests that the virotherapy dosage and the infection rate of the virus are key parameters to ensure long term tumor eradication. However, as treatment was modeled using a constant dosage, future work should include how more realistic periodic treatment schedules influence the resulting stability analysis. Additionally, since sustained oscillations of tumor load are not typically clinically observed, finding criteria for the nonexistence of limit cycles in higher dimensional nonlinear models such as the one presented in this paper will be useful for building larger feasible parameter spaces.

Acknowledgments

This manuscript benefited from the comments of two anonymous reviewers. This material is based upon work supported by Winthrop University's Ronald E. McNair Program and by an Institutional Development Award (IDeA) from the National Institute of General Medical Sciences (2 P20 GM103499 15) from the National Institutes of Health.

Conflict of interest

All authors declare no conflicts of interest in this paper.