MODELING THE SPATIOTEMPORAL DYNAMICS OF VIROTHERAPY AND IMMUNE RESPONSE AS A TREATMENT FOR CANCER

Oncolytic virotherapy is a promising treatment for cancer, which uses viruses that selectively target and kill cancer cells. In this article, we developed a system of partial differential equations, which describes the spatiotemporal dynamics of cancer cells under virotherapy and immune response. We carried out numerical simulations of the system and analyzed the important biological parameters. The numerical results suggest that high viral infection plays an important role in improving the treatment outcomes. Moreover, strong immune response can improve the result of virotherapy along with a high immune killing rate of cancer cells. This model could be validated by using patients’ data; and thus it may help oncologists choose the optimal therapy to eradicate cancer.


INTRODUCTION
Recently, there has been a remarkable development in treatments for cancer. The efficacy of a treatment depends on its ability to eliminate the tumor in a short time, without damaging the surrounding healthy tissue. In the field of genetic engineering, researchers have been able to find a new treatment for cancer with genetically altered viruses [1]. Oncolytic virus is a promising treatment for cancer because it specifically infects the cancer cells and multiplies within them, eventually leading to its lysis. New virus molecules are then released to infect the neighboring cancer cells [2]. Virotherapy is administered into the body via different delivery methods, for example, intratumoral and intravenous [3][4][5][6][7][8]. Intratumoral delivery is restricted to easily reachable solid tumors, whereas intravenous delivery is a good option for advanced, metastatic, or unreachable cancer [9]. Clinical trials have shown success in the treatment of cancer using these two delivery methods [10][11][12][13].
By mathematical modeling we can describe the interaction between oncolytic virotherapy and tumor growth. Moreover, numerical analysis can predict the optimal strategies to ensure the effectiveness of virotherapy. Wu et al. [14] built a mathematical model consisting of a system of partial differential equations (PDEs), describing tumor growth during virotherapy.
They discussed three strategies for injecting the virus: (1) uniform injection (injecting the entire tumor), (2) core injection, and (3) rim injection. They concluded through the analysis of a simplified ordinary differential equation (ODE) model that the virus injected into the tumor has the ability to control its growth. Friedman and Tao [15] developed a PDE model similar to the previous one, where the virus is injected into the tumor, but included a diffusion term for viral density. They analyzed the model to find the effect of injecting the virus on reducing tumor size and found that if viral density is large enough, it minimizes tumor size in a very short time.
One of the most important obstacles that limits the efficiency of virotherapy is the presence of immune response because it destroys both viruses and infected cancer cells (blocking the spread of viral infection) [16]. Some studies have shown how the evolution of nanotechnology has contributed to overcoming this problem. Encapsulating the oncolytic virus by nanomaterials leads to stabilization of the drug in the bloodstream [17]. There are mathematical modeling approaches to understand the immune response against virotherapy. Wein et al. [18] added into their earlier model [14] an immune response that kills virus-infected cancer cells, and they used recent preclinical and clinical data to validate the model and estimate key parameter values. They concluded that an effective immune response contributes to the removal of the virus before complete elimination of the tumor. Therefore, complete control of the tumor is achieved using the oncolytic virus together with immune suppressors. Later, Tao and Guo [19] developed the previous model by adding a diffusion term for both viral density and immune response density [18]. Through their analysis of the model, they found a threshold for the strength of immune response to control tumor growth during virotherapy. Another model, which describes the effect of immune response on infected cancer cell population and virus population, is proposed by Phan and Tian [20]. The analysis of the model shows that if the burst size of the virus is not too big, then the immune response may complicate virotherapy by creating more equilibria.
There is still a growing need to understand how to improve the efficiency of virotherapy by combining it with other treatments. Therefore, some mathematical models have examined the effect of virotherapy in combination with immune suppressors and inhibitors [21,22]. Friedman et al. [21] presented a mathematical model for spherical glioma, a type of cancer, undergoing virotherapy by injecting into its center. The model also includes the effect of cyclophosphamide (a chemotherapy drug that suppresses the immune response). They concluded that tumor size would decrease to being very small even without cyclophosphamide if the burst size is big. Furthermore, repeated cyclophosphamide treatment decreases the density of uninfected tumor cells and thus reduces the risk of a secondary tumor. Tao and Guo [22] developed a mathematical model describing the effect of combined treatment using an oncolytic virus and mitogen-activated protein kinase (MEK) inhibitors on tumor size. MEK inhibitors enhance the coxsackie-adenovirus receptor to ensure the entry of the virus into the tumor cell, but it can also prevent the proliferation of the virus. As a result, they studied the optimal dose and timing of MEK inhibitors. They found that the injection of the virus and MEK inhibitors at the same time achieved more effectiveness for the combined treatment.
Our goal is to understand the dynamics of virotherapy in the presence of immune response.
Cancer cells have the ability to escape from the immune system [16]; however, the presence of viruses could stimulate immune cells and help them recognize and, thus, attack cancer cells [23]. We developed the model by Phan and Tian [20] by considering the effect of the immune system on uninfected cancer cells (to be published). We also added a death term for these uninfected cancer cells. The model consists of four nonlinear ODEs, describing the interaction between tumor growth and virotherapy, along with the effect of immune response. In this article, we develop this ODE model by incorporating the spatial evolution of cells because the virus spreads into the tumor. Therefore, we added a diffusion term for viral density. The main purpose of this article is to determine the optimal strategy for viral treatment. The model would be solved numerically for different values of parameters, where the virus is given at a constant rate for a certain period of time. The latter can be accomplished by nanotechnology [24]. In Section 2, we present the model, and then in Section 3, we carry out numerical solutions of the model for different parameter values. Finally, discussion and future work are presented in Section 4.

MATHEMATICAL MODEL
Our mathematical model describes the dynamics of tumor growth under oncolytic virotherapy in the presence of immune response. In particular, we present a PDE model that takes into account the spatiotemporal evolution of cancer cells. We assume that the domain is spherically symmetric, and thus, the variables depend on time t and radial distance r. The virus is delivered through a blood vessel and diffuses into the tumor, with no flux at the boundary (right). Also, the virus is released at a constant rate for a certain time interval. We hypothesize that cancer cells depend on the closest blood vessel, which has the radius r b . Therefore, we estimate the radius of the spherical region supported by the blood vessel by r b √ BV F , where BVF is the blood volume fraction [25,26].
Supposing that the oncolytic virus diffuses into the tumor, infects the cells, and then becomes inactive, the cancer cells can be divided into uninfected and infected cells. We assume that the uninfected cells grow logistically and die with a constant rate. The infected cells would eventually die, thus releasing new viruses (the virus replicates inside the infected cell). Moreover, we assume that the virus would either die naturally or by the immune system. The presence of cancer cells (uninfected and infected) stimulates the immune system, thus causing their death by natural killer cells. Finally, the immune cells die naturally. Therefore, we have the following system of PDEs: where the homogenous initial and boundary conditions are as follows: The model variables are x = x(r,t) is the density of uninfected tumor cells, y = y(r,t) is the density of infected tumor cells, v = v(r,t) is the density of oncolytic viruses, and z = z(r,t) is the density of immune cells.
The parameters and their description, values, and units are summarized in Table 1. Note that our model is similar to the one in [20], where the new terms involve the parameters D, α, d, and s 2 . Also, we consider the spatial dependence of variables because of the diffusion of the drug.
Note that v 0 can depend on the time at which the virus is injected according to a specific dosing and timing regimen (e.g., bolus treatment). Here, we only consider a drug that is continuously administered for a period of time at a constant dose v 0 , as suggested by Borges et al. [27] for chemotherapy. This can be done by nanotechnology, which reinforces the stability of the drug in the bloodstream [17].
We calculate f x that represents the ratio of the uninfected (viable) cancer mass to its initial mass. This is done by integrating the density x at each time step over the spherically symmetric domain around the blood vessel, as follows: where Similarly, f y , f v , and f z are calculated.
Before varying the parameters and considering the effect of the treatment on cancer growth, we solve the model (1) with the conditions (4), (5) and parameter values from Table 1, as shown in Figure 1. The time t varies from 0 to 504, which represents virotherapy given for 3 weeks.
The plots in Figure 1 Table 1. In (a)-(d), the variables are plotted after the 2nd, 4th, and 6th day of the treatment as indicated in the legend (t is given in hours). In (e), the ratio of the uninfected cancer mass to its initial mass is plotted for the whole treatment which is 3 weeks.  Table 1). In all simulations the initial and boundary conditions are as given in (4), (5).

Parameter analysis.
In this section, we vary the parameters in the model. We focus on the infection rate of the virus (β ), burst size (b), immune killing rate of uninfected cancer cells (α), clearance rate of viruses (γ), and stimulation rate of immune cells by uninfected cells (s 2 ).
The parameters β , b, and γ depend on the type of the virus, and thus, by varying them we consider different types of treatments. Moreover, α and s 2 are related to the immune response to cancer (uninfected cells). Therefore, by varying these parameters we consider cases where the immune system is stimulated by the presence of the virus or by, for example, monoclonal antibodies [32] and vaccine [33].
First, in Figure 2 Then, after a long time, f x would reach smaller values than those in the other two cases (red and blue curves). Perhaps as viruses die fast (which means less viruses), the immune system kills uninfected cells. Finally, we vary the rate at which uninfected cells stimulate the immune system (s 2 ), as demonstrated in Figure 2(e). It is expected that as we increase this rate, f x would decrease and reach zero in a short time, as shown in the figure (yellow curve). This is because the strong response of the immune system would lead to the killing of uninfected cells.

DISCUSSION AND FUTURE WORK
Virotherapy is a promising therapeutic approach for the treatment of cancer that uses genet-  [32] and vaccine [33]. Another combined therapy could include virotherapy and chemotherapy [36] or radiotherapy [37].

Conflict of Interests
The author(s) declare that there is no conflict of interests.