A mathematical model for tumor growth and treatment using virotherapy

: We present a system of four nonlinear di ﬀ erential equations to model the use of virotherapy as a treatment for cancer. This model describes interactions among infected tumor cells, uninfected tumor cells, e ﬀ ector T-cells, and virions. We establish a necessary and su ﬃ cient treatment condition to ensure a globally stable cure state, and we additionally show the existence of a cancer persistence state when this condition is violated. We provide numerical evidence of a Hopf bifurcation under estimated parameter values from the literature, and we conclude with a discussion on the biological implications of our results.


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.

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 longterm 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 welldocumented 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, β to represent the rate of uninfected tumor cells becoming infected, γ U to represent the rate of decay of uninfected cells via T-cells, γ I to represent the rate of decay of infected cells via T-cells, δ I to represent the rate of decay for infected cells, c to represent the rate of T-cell growth via infected tumor cells, δ E to represent the rate of decay for effector T-cells, δ V to represent the rate of decay for virions, α 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 , and δV = δ V r .We drop all the tildes from our notation and arrive at the following non-dimensionalized model: (2.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.

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, lim sup t→∞ U(t) ≤ 1.Using this upper bound, we can derive an upper bound for the infected tumor cell population.From (2.1), we have: It

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 * = N δ V .Thus we find a cure state equilibrium point of the form P 0 = (0, 0, 0, N δ V ).
If U * 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 δ 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 δ V > N. We summarize these results in the following theorem: Theorem 1.

Stability of the cure state
In this section, we explore the stability of the cure state equilibrium P 0 = 0, 0, 0, N δ V .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.

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 Thus, we find the local stability condition for the cure state to be N > δ V .If this condition is not met, the cure state is unstable.We demonstrate this result with numerical simulations in Section 4.

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 lim inf We next seek conditions to ensure U(t) → 0 as t → ∞.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 We note that this is the same condition for local stability of the cure state.Under this condition, it then follows that U(t) → 0 as t → ∞.
We next explore the asymptotic behavior of I: Thus U(t) → 0 implies I(t) → 0 as well.Similarly, for E and V we have: Hence I(t) → 0 implies E(t) → 0 and V(t) → N δ V .These results are summarized as follows: Theorem 2. If N > δ V , then the cure state P 0 = 0, 0, 0, N δ V is globally asymptotically stable.

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 < δ 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.

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 × 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 β = 8.9 × 10 −4 by the total cell population; thus we similarly adjust the value of β to align with the chosen carrying capacity.The corresponding non-dimensionalized parameter values are found in Table 2.  2, Theorem 2 guarantees a globally stable cure state for Ñ > 7.41935 (corresponding to N > 8.01123 × 10 11 virions), while the Routh-Hurwitz criteria above predicts a locally stable cancer persistence state for 0.000819195 < Ñ < 7.41935 (corresponding to 8.84546×10 7 < N < 8.01123×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 Ñ < 0.000819195 (corresponding to N < 8.84546 × 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 Ñ = 0.000819195, (correspondingly N = 8.84546 × 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: Then the dynamics undergo a Hopf bifurcation at (x 0 , µ 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 , µ 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 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 ±ω 0 i.Let q be the eigenvector of A corresponding to the eigenvalue ω 0 i.Let p be the adjoint eigenvector such that A T p = −ω 0 ip and p, q = 1.If I denotes the identity matrix, then the first Lyapunov coefficient l 1 (P * ) of the system dx dt = F(x) at P * is given by l 1 (P * ) = 1 2ω 0 Re[ p, C(q, q, q) − 2 p, B(q, A −1 B(q, q)) + p, B(q, (2iω 0 I n − A) −1 B(q, q)) ] 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 Ñ = 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 (λ 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 ω 0 = 1.188783 and q . To calculate p, we derive the eigenvalues and corresponding eigenvectors of A T : λ i = −11.0387,−1.73456, ±1.18878i; We then find that p = When we compute p•q, we see that it is not equal to 1, but is equal to −0.000313395−0.000303225iso we must scale q as follows:
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 * ): 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 × 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.

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 2 1 + a 2 3 a 0 .If N 0 < N < δ 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.
varied With the non-dimensionalized parameters from Table