The Effect of the Infection Rate on Oncolytic Virotherapy

Oncolytic viruses have become a novel therapeutic tool for various cancer treatments. Several naturally occurring oncolytic viruses and engineered oncolytic viruses are developed for oncolytic virotherapies. Although we have a good understanding on molecular mechanisms of viral replication and virus-induced cell lysis at the cellular level, it is unclear how oncolytic viruses and cancer cells interact as a population. Several mathematical models of oncolytic virotherapy have been developed to advance the understanding of dynamic interaction between oncolytic viruses and cancer cells. Many authors investigated the effect of the virus replication on dynamics of cancer cell population and proposed that the bursting rate of viruses is an important factor for successful oncolytic virotherapy. In this study, we investigate the effect of infection rate of oncolytic viruses on an oncolytic virotherapy model. Particularly, we focused on studying the relationship between two control parameters, bursting rate and infection rate of the virus, to generate the patterns from equilibrium steady state to periodic solutions. Based on the model, the interaction between cancer cells and oncolytic viruses shows an intriguing two-dimensional bifurcation, showing three parameter regions (equilibrium steady state, damped oscillations and oscillations). Our result suggests that both infection rate and bursting rate are crucial properties of oncolytic viruses to design a successful oncolytic virotherapy.


Introduction
Oncolytic virotherapy (OV) is a novel treatment by which oncolytic replication competent viruses are either inoculated to a cancer patient or directly injected into the tumor to kill cancer cells without causing harm to healthy normal cells. Upon introduction, viruses selectively infect cancer cells by binding to a specific receptor protein present in cancer cells. The antitumor activity of OVs can be achieved either by direct lysis of cancer cells or by triggering a host immune response [1][2][3]. Upon the lysis of infected cancer cell, new virus particles burst out and infect adjacent cancer cells. Recent studies have shown that measles, adenovirus, Newcastle-disease virus (NDV), reovirus, herpes simplex virus 1 (HSV-1), and coxsackievirus are relatively non-toxic and tumor specific [4,5]. In addition to naturally occurring viruses, viruses have been engineered to increase the specificity. For instance, an adenovirus H101 was the first engineered virus that have been approved by a government regulatory agency [6].
Many mathematical models of OV have been developed to understand complex and dynamic interactions between tumor cells and viruses. In addition, mathematical models can give an insight to design treatment regimens for a cancer intervention with minimal adverse effects. Many studies have focused on examining how variation in viral and host parameters influence the outcomes of treatments [7][8][9][10][11][12][13][14]. Wodarz et al. proposed a mathematical model to investigate the effect of virus-specific lytic cytotoxic T lymphocytes on the dynamics of OV [8]. The ordinary differential equations (ODEs) were used to describe the interactions between two types of cancer cells populations (uninfected cancer cells and infected cancer cells) according to predator-prey model including virus population. There are also mathematical models using partial differential equations (PDEs) [12][13][14].
It has been suggested that the ability of the viral replication, represented as a bursting rate, is a key factor for successful OV. In addition, the oscillations of cancer cell and virus populations are an intrinsic property of OV [11]. Compared to the ability of virus replication, little attentions have been given to the effect of infectivity, represented as an infection rate, of the virus. Furthermore, it will be crucial to understand the relationship between bursting rate and infection rate of viruses for the development of optimal OVs. In this paper, we investigated how the variableinfection rates impact-on the cancer cell population using dynamical system (phase space analysis) and bifurcation approaches. We also analyzed the population dynamics of oncolytic viruses and cancer cells while considering (or applying) both bursting rate and infection rate in a simple OV model. Our results suggest that 1) the infection rate coefficient is a bifurcation parameter; existence of two threshold values, 2) these threshold values of infection rate are dependent on the bursting rate values and 3) there are a set of values or parameter region in choosing the bursting rate and infection rate to the therapeutic way to succeed or fail.

Model
The OV model is a three dimensional [11] Where (t) , (t) and (t) represent populations of uninfected cancer cells, infected cancer cells and free viruses, respectively. The is the growth rate of cancer cell and being the carrying capacity of a tumor. The term 1 − describes the logistic growth rate of an uninfected cancer cell population (t) . The constant value represents the infection rate of the virus in the infected cancer cells so that the term describes the rate of infected cancer cells by free viruses (t).
represents the death rate of infected cancer cells. The is the bursting rate of free virus particles. The term is the rate of elimination of free virus particles by various causes including nonspecific binding, generation of defective interfering particles, and host immune system.
For non-dimensionalization, we set τ = , = , y = , z = . Then Equation (1) become We have the following model by setting the parameters; = + − * − All parameters are described in Table 1.

Analysis and Stability of Equilibrium
There are three fixed points (equilibrium points); two equilibrium solutions -(0, 0, 0) and -(1, 0, 0) in the positive invariant domain D and the other one -( * , * , * ) in either the negative or the positive domain depending on the parameter values.
Theorem Therefore, E • E > E , 2E((I + E) + K > 0 , which implies E > 0. Thus, the complex roots have negative real parts, which means E 3 is stable.

Numerical Simulation
We used the Runge-Kutta 2nd order method [15] to compute the numerical solutions in MATLAB (The Mathworks, Natick, MA). The small-time step ∆t = 0.05 was used to check the accuracy of the numerical method. For numerical simulations, we set the parameters a = 0.31 (growth rate of cancer cell from the experimental data [16]) and d = 0.44 with initial conditions x(0) = 0.5, y(0) = 0, and z(0) = 1.5. Both the infection rate (b) and the bursting rate of virus (c) are considered as variables.

Infection Rate Coefficient Is a Control Parameter in Changing the Structure of Cancer Cell Population
We investigated the effect of the infection rate of the virus on the cancer cell population. We set the busting rate (c = 5; different values of c will be discussed later) and increased the infection rate constant from 0 to 1 with step size 0.01. Our numerical result shows that the increased infection rate (b) results in the altered the pattern of cancer cell population over time. The cancer cell population exhibits three patterns over the relative time depending on the infection rate coefficient ( Figure 1A). For the infection rate b = 0.1, the cancer cell population approaches to the free virus equilibrium -(1, 0, 0) in phase portrait ( Figure 1B) or maximal capacity (green curve in Figure 1a). For b = 0.4, it shows fluctuations over short time period and then converges to the relative population less than the maximal capacity (damped oscillation, Figure 1C). For b=0.64, the cancer cell population converges to the minimum population 0.172 (data not shown). For b = 0.7, the population starts to oscillate over time (red curve in Figure 1a) and the trajectory of populations circulates around a stable limit cycle in a dynamical system ( Figure 1D). As the infection rate increases, the amplitude of oscillation of the cancer cell population also increases. This result indicates that the infection rate coefficient should be considered as a control parameter that can change the dynamic structure of the relative cancer cell population over time as it changes.

Transcritical and Hopf Bifurcations
Qualitative changes in the stability of fixed points depending on a parameter are called bifurcation in the dynamics 17 . The control parameter values where bifurcations occurred are called bifurcation points. A bifurcation diagram of cancer cell population can be obtained by changing the parameter of virus infection rate. We computed the relative cancer cell population density with different values of b over the relative time. The bifurcation diagram for the cancer cell population over the infection rate is shown in Figure 2. There are two bifurcations at b = b 1 and b 2 where the shape of the graph changes. In our numerical result, we calculated two values b 1 = 0.12 and b 2 = 0.64 (Figure 2A). We also calculated the Lyapunov exponent to determine whether cancer cell population for b > b 2 shows periodic solutions or non-periodic solutions ( Figure 2B). The cancer cell population density in this parameter interval shows periodic patterns since the Lyapunov exponent is less than 0 for b > b 2 . For 0 < b ≤ 0.12, the cancer cell population reaches its maximum population (capacity). This indicates that the viral population diminished due to low strength of infectivity of the virus so that the cancer cell population approaches to the maximum. In this case, it can be concluded that the virotherapy fails in that parameter regions. For 0.12 < b ≤ 0.64, the cancer cell population shows oscillations for a short time period and then approaches to the values less than 1. Particularly, the population decreases to the relative population of 0.172 at b=0.64, which means that the virotherapy is partially successful in the interval and the cancer cell population can be reduced to minimum population of 0.172. For b > 0.64, the population start to show oscillations over the relative time. Populations of the infected cancer cells and viruses also oscillate (data not shown). Oscillations of cancer cells and viruses suggests that the interaction between cancer cells and viruses are the same dynamics as the prey-predator system (coexistence in dynamical system).
A brief bifurcation diagram in terms of fixed points, discussed in section 2.2, is shown in Figure 2C. A horizontal line is used to describe the fixed points of -(0, 0, 0) and -1, 0, 0 and a curve is used to describe the parameter dependent fixed point -0 , 0 , 0 . Here, solid and dotted lines represent stable and unstable, respectively. For 0 J * = * , -is a stable fixed point and bothandare saddle (unstable), so that their populations approach to the free virus equilibrium valuein dynamical system. Particularly,is in the negative invariant domain that is outside of our interest region. However, -is in the positive invariant domain D and it is in the same position asat the bifurcation point b = b 1 . In other words, -moves closely to from the negative invariant domain to the positive invariant domain D as the b value increases from 0 and then is in the position ofat b=b 1 . For b > b 1 , the stability of andchanged from stable to unstable and from unstable to stable, respectively. This indicates a transcriticial bifurcation occurs at b=b 1 . For * J * = * , bothandbecome unstable fixed (saddle) points butbecomes a stable spiral fixed point, so that the trajectories circulate about and then gradually converge to -. For b > b 2 , all fixed points are unstable fixed points and the eigenvalues of Jacobian matrix athas a family of complex numbers with positive real part. The trajectories rotate around a stable limit cycle near the fixed pointand the populations show oscillations, which implies that there is a Hopf bifurcation occurs at b=b 2 .

Bursting Rate-dependent Bifurcation Values
We investigated the effect of the correlation between the bursting rate and the infection rate of the virus on the cancer cell population. We simulated the bifurcation diagrams for the infection rate coefficient with various fixed bursting rate coefficients to examine how the bifurcations points (b 1 and b 2 ) change ( Figure 3). For the fixed value of c=4, the cancer cell population undergoes the first (transcritical) bifurcation at b 1 =0.146 and the second (hopf) bifurcation at b 2 =0.82 ( Figure 3A). The other values are given by (c, b 1 , b 2 ) = (7, 0.073, 0.45), (15, 0.031, 0.205), and (30, 0.015, 0.102) ( Figure 3B-3D). These results provide increasing the bursting rate of the virus results in decreasing the bifurcation points (b 1 and b 2 ). Those bifurcation values are dependent on the bursting rate values, which implies that different values of bursting rate induce different bifurcation points in the infection rate. These results leave questions about how the correlation between the bursting rate of the virus and the infection rate of the virus affects the success or failure of the oncolytic virotherapy.

Parameter Regions from Equilibrium Steady State to Periodic Solutions
The oncolytic model undergoes two bifurcations with respect to two parameters; the bursting rate and the infection rate of the virus. The transcritical bifurcation occurs at (b, c) = (b 1 , c 1 ) where the equation form of the bifurcation point clearly can be formulated from the stability of the fixed point -1, 0, 0 and is given by *1 :5 . In other words, the set of bifurcation points (b 1 , c 1 ) are the set of points on the graph of + 6 6 in the b-c plane. Therefore, the set of parameters (b, c), which satisfies the inequality + = 6 6 , represents the region where the cancer cell population converges to the maximum carrying capacity and the virotherapy is a failure. Otherwise, the therapy is partially successful, or all populations coexist showing oscillations over the relative time. Figure 4A shows the two-dimensional bifurcation diagram. The diagram shows three different colored regions which represent the following; 1) The blue region is the set of parameters (b, c) at which the relative cancer cell population converges to the maximal capacity, 2) the yellow one is the region where the population exhibits damped oscillations and converges to the population less than the maximal capacity and 3) the red one is the region where the population shows oscillations over the time. Our numerical result provides the critical set of values of bifurcation points (b 1 , c 1 ) and (b 2 , c 2 ) which are located on the border between two different colored regions. For example, the transcritical bifurcation occurs at the border between blue and yellow color and the Hopf bifurcation occurs on the border between yellow and red colored regions.

Discussion
OVs are promising new therapeutics for various types of cancers. There have been multiple clinical trials using oncolytic viruses [18]. For the OV, replication competent viruses are used such as reoviruses [19]. The infection oncolytic viruses can kill cancer cells by triggering a cytotoxicity directly in the infected cancer cells or antitumoral immune responses or both [18,20]. For an effective and safe treatment with oncolytic viruses, it is important to understand the behaviors of cancer cells and viruses. Mathematical modelling is a useful method to study the dynamic interactions between cancer cells and viruses. Although it is a simplified model excluding many other variations that can happen in clinical settings, our model supports the idea that not only the bursting rate of viruses but also infection rate of viruses must be considered for OV. One of barriers in OV is finding the optimal number of viruses to exterminate cancer cells without adverse effects such as immune hypersensitivity or spontaneous infection of normal cells [21]. Our model can be useful to determine the optimal number of viruses to be used for OV based on the bursting rate and infection rate that can be expressed as a plaque-forming unit (pfu). Both the bursting rate and pfu of a specific virus can be determined by in vitro experiments.
In clinical setting, there can be many other factors to affect the interaction between cancer cells and viruses. One of factors have not been considered in our model is host immune response. The host immune system either can facilitate OVs by triggering cytotoxic response to kill infected cancer cells or can interfere OVs by eliminating viruses. Intrinsic properties of virus can also determine the lifespan or death rate of free virion. Based on figure 5, it will be important to develop advanced models to consider the death rate of viruses when assessing cancer cell-virus infections.

Conclusions
In this paper, we studied the effect of the infection rate of the virus on the cancer cell population using the bifurcation diagrams of the cancer cell population density and the dynamical system (phase-space analysis). Our results show that there are two bifurcation threshold values of infection rate, which are the bursting rate dependent. The first bifurcation threshold value (b 1 ) can be obtained from the stability of the fixed point at -1, 0, 0 and be expressed as *1 :5 , which has an inverse relationship with the bursting rate of the virus and is proportional to the death rate of free virus particles. However, it is difficult to compute the second threshold value explicitly and express it using parameter values. For this issue, we provided the two-dimensional bifurcation diagram to show the set of parameters where three different patterns are generated. Our result gives an insight into understanding the relationship between two control parameters at which the cancer cell population exhibits the patterns from equilibrium steady-state solution to periodic solutions. It provides the second threshold value (b 2) has the same functional inverse relationship with c 2 . The two-dimensional bifurcation diagram provides a promising result in determining optimal parameters for the successful virotherapy.