A Cellular Automata Model of Oncolytic Virotherapy in Pancreatic Cancer

Oncolytic virotherapy is known as a new treatment to employ less virulent viruses to specifically target and damage cancer cells. This work presents a cellular automata model of oncolytic virotherapy with an application to pancreatic cancer. The fundamental biomedical processes (like cell proliferation, mutation, apoptosis) are modeled by the use of probabilistic principles. The migration of injected viruses (as therapy) is modeled by diffusion through the tissue. The resulting diffusion–reaction equation with smoothed point viral sources is discretized by the finite difference method and integrated by the IMEX approach. Furthermore, Monte Carlo simulations are done to quantitatively evaluate the correlations between various input parameters and numerical results. As we expected, our model is able to simulate the pancreatic cancer growth at early stages, which is calibrated with experimental results. In addition, the model can be used to predict and evaluate the therapeutic effect of oncolytic virotherapy.


Introduction
Oncolytic virotherapy is a novel cancer treatment where natural or genetically modified viruses infect cancer cells and then self-replicate until the host cancer cells lysis (see Fig. 1). Ruptured cancer cells release chemicals like tumor antigens, which make cancer cells easily recognizable by the immune system. Moreover, the released viruses can infect more cancer cells to trigger a chain reaction and effectively act as a followup treatment. As early as in 1912, De Pace (1912) observed a tumor regression after inoculation of an attenuated rabies vaccine in a patient with uterine cervical carcinoma. Later on, an animal-based test (Levaditi and Nicolau 1922) and a human trial (Pack 1950) were conducted in 1920 and 1940, respectively, where both experiments yielded an obvious partial tumor regression (Kasuya et al. 2005). In the subsequent decades, more works (Kirn 2001;Gil et al. 2013Gil et al. , 2014 demonstrated that oncolytic virotherapy lead to tumor attenuation. Some milestones in the development of oncolytic virotherapy are shown in Fig. 2 (De Pace 1912;Martuza et al. 1991;Xia et al. 2004;Fukuhara et al. 2016).
Pancreatic ductal adenocarcinoma (PDA), recognized as the most common pancreatic cancer, is a lethal disease due to late detection, a low resectability rate, medication resistance and poor prognosis (Korc 2007;Gore and Korc 2014). Currently, pancreatic cancer is the seventh leading cause of cancer death worldwide and its 5-year survival rate is less than 5% (Feig et al. 2012;Bray et al. 2018). Compared with other types of cancer, PDA has more cancer-associated fibroblasts (CAFs) resulting in abundantly desmoplasitic stroma that constitutes up to 90% of a solid tumor volume (Moir et al. 2015;Öhlund et al. 2017). The profuse desmoplasia in stroma produced by CAFs acts as a physical barrier for drug delivery and leads to medication resistance (Provenzano et al. 2012;Jacobetz et al. 2013). However, CAFs make cancer cells more susceptible to be infected by oncolytic viruses. Ilkow et al. (2015) experimentally demonstrated that the cross-talk between cancer cells and CAFs facilitates the oncolytic virus-based therapies. Therefore, oncolytic virotherapy offers an avenue for the treatment of pancreatic cancer.
The ideal oncolytic virus should be able to selectively replicate itself in cancer cells without damaging normal somatic cells (see Fig. 1). A couple of studies (Sunamura et al. 2004;Kasuya et al. 2005;Wennier et al. 2011) summarized advantages and disadvantages of various replication-competent oncolytic viruses proposed for pancreatic cancer therapies, e.g., adenoviruses, herpesviruses, poxviruses, parvoviruses, reoviruses and paramyxoviruses. A few types of viruses have been tested in animalbased xenograft models; however, only a few kinds of viruses have reached clinical trials. In particular for pancreatic cancer, relevant studies are rare, among which Fu et al. (2006) observed that a type of oncolytic virus produced antitumor effects in human pancreatic cancer xenografts. In addition, Sunamura et al. (2004) carried out adenovirus therapy in immunodeficient mice with human pancreatic cancer xenografts that resulted in a remarkable inhibition of tumor growth under consecutive injections of the virus. Typically, if animal testing is successful, the new drug will reach clinical trials that are classified into four phases: (1) phase I is a safety test for healthy volunteers; (2) phase II demonstrates whether a drug can have any efficacy against the disease; (3) phase III checks in a randomized multicenter tests if a drug has the right therapeutic effect; and (4) phase IV post-marketing surveillances its efficacy and side effects after extensive use. Regarding the clinical applicability, Kasuya et al. (2005) stated that a clinical trial of viruses adenovirus ONYX-015 (phase I and II) has been conducted in pancreatic cancer patients, where half of the patients (phase II) exhibited either tumor reduction or stabilization. In contrast, a phase I trial of the efficacy of several oncolytic herpes viruses (such as G207, 1716 and OncoVEX GM-CSF) have tested against various tumors and the herpesvirus exhibited a good tolerance at all dosages. Although oncolytic virotherapy has been proposed for decades, a thorough understanding of the interactions of virus, tumor and microenvironment in vivo is still needed to be further researched, like a proper viral dose for a specific virus, how to control of virulence and etc. Therefore we develop a three-dimensional (3D) spatial Markov Chain cellular automata model to mimic pancreatic tumor (pancreatic cancer at early stages) progression and subsequent regression under the interference of oncolytic virotherapy. In addition, the model can be used to quantify the impact of virotherapy with different viral doses, viral infectivity and levels of immunity in patients with pancreatic cancer.
Cellular automata models are lattice-based models that facilitate analysis of the spatiotemporal dynamics based on the interplay between cells and their microenvironment. The cellular automata model has been introduced as a computer model of self-reproduction by John von Neumann and Stanislaw Ulam (Langton 1984). In the past decades, cellular automata models, in addition to self-reproduction, have been extended to other model applications successfully within a wide range of spectrum from biology, physics, chemistry to other sciences (Deutsch et al. 2005). Regarding the cancer modeling in cellular automata, Reis et al. (2009) proposed a model that could capture the Gompertzian behavior of tumor growth. Hatzikirou et al. (2008) developed a model of tumor invasion dynamics. In addition, a couple of studies demonstrated applications of cellular automata in cancer therapy, e.g., radiotherapy (Enderling et al. 2010), and chemotherapy (Pourhasanzade and Sabzpoushan 2019). In terms of computational models on oncolytic virotherapy, Wodarz et al. (2012) studied the distinct patterns of oncolytic viral spreading through tumor cell population both experimentally and computationally. Furthermore, Paiva et al. (2009) developed a multiscale cellular automata model for oncolytic virotherapy, where authors found that a solid tumor can be either eradicated completely or keeps on growing. The resulting behavior depends on the input parameters, which represent the biological circumstances around the tumor. Malinzi et al. (2017) developed a partial differential equation-based model for oncolytic virotherapy. For the one-dimensional case, they assessed the stability of traveling wave solutions and the rate of tumor progression by means of the minimal traveling wave speed. One of their most important conclusions is that the combination of chemotherapy and virotherapy most successfully removes tumors.
The 2D model in Paiva et al. (2009) revealed an oscillatory behavior of cancer cells and virus population, which hints a strong host immune response is necessary. As an extension of this model, we develop a 3D model to phenomenologically show pancreatic cancer initiation at early stages and its subsequent regression under oncolytic virotherapy, where the various levels of host immune responses and toxicities of residual viruses are taken into account. To the best of our knowledge, it is the first description of a hybrid cellular automata model with an application to pancreatic cancer. In addition, this model also shows the pattern of virus spreading by using partial differential equations in a spatial domain. One of the advantages is its efficacy in evaluation of therapeutic outcomes based on injected dose of viruses, viral infection potent and personal immune strength. Next to the development of the model, we calibrate the model to experimental outcomes and we carry out an uncertainty quantification by using Monte Carlo simulations in terms of assessment of the likelihood of success of the treatments. Therefore, therapeutic effects as well as toxicities of residual viruses are predicted, which is expected to be helpful for viral administration and the patientspecific treatment.

Mathematical Formalism
Cellular automata models consist of a class of lattice-based models, where lattice approaches are classified as: (1) a single lattice site is occupied by one cell only; (2) a single lattice site is occupied by a cluster of cells; (3) one cell takes many lattice sites. They are all capable of investigating biological processes with single cell or multiple cells, where division, death or other biomedical phenomena are modeled by stochastic processes (Van Liedekerke et al. 2015). In contrast, the first two categories are typically used to describe volume effects, whereas the last category is able to capture the morphological evolution of cells.
In our simulations, each lattice represents a volume element filled with multiple cells in a bounded 3D computational domain ⊂ R 3 , which is divided into a set of lattice points N = {1, . . . , n}. The boundary of is denoted by . The lattice point i has a finite number of discrete states S i that indicates the state of cells in the corresponding volume, which reads as (1) Assigning an initial state for each lattice point, and subsequently adjusting the state of the specified lattice at position x i = [x i , y i , z i ] at subsequent times is correlated with the states of its neighborhood marked in Fig. 3. Subsequently, we consider fundamental biological processes like cell division, mutation, infection and death modeled as stochastic processes. The likelihood of changing a state of a lattice i to another state satisfies a memoryless exponential distribution which is given by f (λ i , t) t during a time period (t 0 , t 0 + t). Here λ i is the probability rate per unit of time, which is determined by biological mechanisms like cell division, mutation, infection and death. The f (λ i , t) reads as, and hence the transition probability P within a time interval of length t is given by In the work by Vermolen and Pölönen (2020), it is proved that the likelihood of changing state depends on a simple binary states of the neighbors, which has been applied to modeling the progression of skin cancer. All lattice points in the domain are initialized to the epithelial cell state {S i = 1} or unoccupied state {S i = 0}. During the growth of epithelial cells, cancer mutation happens as a result of exposure to carcinogenic factors like genetic inheritance, chemical carcinogens, electromagnetic radiation, or viral infection. The very complicated biological process is simplified by the application of a transition probability from epithelial cells to uninfected cancer cells, that is from {S i = 1} to {S i = 2} over a time interval t with a likelihood following Here λ mu represents the mutation probability rate per unit of time and the second part of Eq. (4) reflects that this transition is irreversible. Since both epithelial cells and cancer cells are able to proliferate and to migrate, lattice points are allowed to change their states from {S i = 0} to {S i = 1, 2}. The likelihood for these transitions over a time interval t are given by where λ pro denotes the probability rate of transition from 'not occupied by any cells' to being 'occupied by either epithelial cells or uninfected cancer cells'. Note that the probability rate λ pro is determined by the number of neighbors that are in cancer/ epithelial cell state. That is λ pro = λ max n 12 h , where n 12 denotes the number of lattice points that are either in state 1 or in state 2. The distance between two lattice points is represented by h and λ max is a constant to regulate the overall growth rate of cells (including epithelial and uninfected cancer cells). Whether a free lattice point will be occupied by multiple epithelial cells or uninfected cancer cells depends on the states of the surrounding lattice points. Consider an unoccupied node i at time t, that is S i (t) = 0. We denote the number of neighboring lattice points that are in state 1 by n 1 . We further denote the likelihood that the lattice point i, given that it changes state, by α = α 0 · n 1 n 12 , where α 0 ∈ [0, 1] is a constant. This fraction α is used to determine the transition probability of node i, which is given by Apoptosis is programmed death of cells, however, uninfected cancer cells are able to proliferate uncontrollably and to resist cell apoptosis. In the current model, oncolytic viruses are incorporated to infect and damage cancer cells. We let cancer cells jump from state {S i = 2} (uninfected cancer state) to state {S i = 3} (infected state) as soon as the virus concentration exceedsĉ. Hence, we have If the viral concentration does not exceed the thresholdĉ, then we disregard the release of viruses. Therefore, the likelihood of a lattice point i to be infected by viruses, which is the state transition from {S i = 2} to {S i = 3}, depends on the released viruses from the neighborhood and the local concentration of viruses. Subsequently, infected cancer cells are, like epithelial cells, subject to possible cell death. Hence a node i is allowed to change from a cell state {S i = {1, 3}} to an unoccupied state {S i = 0}, which is given by the following likelihood Here λ de denotes the probability rate that an infected cancer or an epithelial cell dies.
Oncolytic virotherapy is initiated when the fraction of tumor constitutes up to 100% of the computational domain (that is the tissue). In animal-based experiments, the viruses are given by injections (Sunamura et al. 2004). Therefore, we consider one or multiple injections, denoted by V(t), as source points at position x p by using the Dirac delta function δ(x) at time t. After the injection, the delivery of viruses is simulated by the reaction-diffusion equation written as where c, D and γ denote the concentration, diffusivity and injection rate of virus in this domain. Note that u(x) is utilized to model an increase in production of viruses released by necrotic cancer cells. Moreover, T represents the mass transfer rate coefficient between the computational domain and its environment. Since viruses infect cancer cells and copy themselves until host cells lysis, infected cancer cells act as sources where newborn viruses originate. We define ic (t) to denote the portion of the computational domain that is occupied by virally infected cancer cells. The function u(x) in Eq. (9) increases as where β and N v denote the proliferation rate of virus and a burst size of viruses from a necrotic cancer cell. The function u(x) is a hypothetic function that accounts for the regeneration of virus particles as long as the cancer cell has been infected. We assume that the carrying capacity of the viral particles is determined by the availability of limited amounts of resources in the cancer cell. This assumption translates into the above logistic differential equation. This function can be revised easily if experimental observations require this.
Since most clinical data are not available and some parameters have even never been measured, we estimate some of the input parameters, which are listed in Table 1.

Discretization
We first consider a 2D or 3D that is occupied by uninfected cancer cells, where a dose of oncolytic viruses is injected into the computational domain as part of the therapy. Subsequently, viruses diffuse and thereby spread throughout the domain. The change in concentration of viruses is modeled by the reaction diffusion equation Note that the above equation is a combination of Eqs. (9) and (10), in which we only take one injection source into account. The numerical method is described for the case of only one injection source since our simulations are only carried out for one source only. We note that incorporating multiple injection sources would not change the numerical method conceptually.
To solve the problem in 2D, the Laplace operator is discretized by the finite difference method (FDM) as, where h is the distance between adjacent lattice points. If the computational domain is extended to 3D, the above equation needs to be revised to If a lattice point on boundary x=0 , then the point (−h, y, z) is assumed as a virtual point out of the computational domain. The Robin boundary condition in Eq. (11) is dealt with the FDM as, and thereby the viral concentration at the virtual point is calculated by, Analogously, the viral concentration on the other virtural points can be obtained by Similarly, if a lattice point is located on boundary x=x n , y=y n or z=z n , the viral concentration at the corresponding virtual point is estimated by (17) Furthermore, the injection of viruses is simulated by a point source that is mathematically inspired by the Dirac Delta function δ(x), which is mollified by using the normal distribution, where ε and d, respectively, denote the source width and the dimensionality.

Time Integration
To update the concentration of virus on each lattice at the next time step, an IMplicit-EXplicit(IMEX) time integration is utilized, where the linear parts and nonlinear parts are treated by a Euler backward method and a Euler forward method, respectively. Thereby the concentration of virus c is updated by Note that this IMEX approach avoids the need of inner iterations to solve a nonlinear problem at each time step. We have used the finite difference scheme on several mesh sizes and time steps and we find that enlarging the resolution dose not yield any significant changes in the approximation for the current numerical setting (see Table  2).

Monte Carlo Simulation
Monte Carlo simulations are widely used in many quantitative probabilistic and statistical investigations that permeates much of finance, engineering and contemporary sciences (Kroese et al. 2014). To obtain quantities of interest, such as (cumulative) probability distributions of output variables and correlations, Monte Carlo simulations enable random sampling of input parameters from predefined probability distributions and extensive repetitive experiments. Due to the variety of viruses and variations from patient to patient, many variables can hardly be determined or measured. For instance, the dose of a virus γ may depend on its effectiveness and toxicity, which varies among viruses. In addition, some variables, such as the concentration threshold at the time of viral infectionĉ, the reproductive rate of the virus in cancer cells α and human immune strength β, may all depend on patient lifestyle, gender and genetic pattern, and hence vary from patient to patient. However, the above-mentioned variables may be quantitatively correlated to viral treatment outcomes, and thereby Monte Carlo simulations are performed on X ∈ {γ,ĉ, α, β}. We assume that X follows a normal distribution X ∼ N(μ, σ 2 ) with the mean value μ and the standard deviation σ . Therefore, the stochastic variable X with a number of samples N s follows Taking the computational time into account, Monte Carlo simulations are performed in two-dimensional simulations. Furthermore, Monte Carlo algorithms tend to be scalable and rely less on computational dimensionality. Referring to our previous work in (Chen et al. 2018a), the accuracy of Monte Carlo simulation is proportional to the reciprocal of the square of the number of samples N s , therefore 5000 samples are chosen to guarantee a small error.

Numerical Results
To perform the simulations, the most numerical parameters have been listed in Table 2.

Cancer Progression
We consider a cubic domain to represent the tissue in the pancreas. The 3D domain has been divided into N l lattice points. Each lattice point is occupied by multiple where Vol( ) is the domain volume and N c (t) denotes the number of lattice points in cancer state {S i = 2} at time t. To model cancer mutation occurring at the edge of a tissue or organ and its competitive growth with epithelial cells, a small number of lattice points are generated randomly in one octant of the domain only (see Fig. 4a). Those lattice points are initiated with epithelial cell state indicated by blue color. Due to mutation, several lattice points change their states from {S i = 1} to {S i = 2} that are visualized by the red dots (see Fig. 4a). Typically, normal cells stop dividing once they contact with each other during division as a result of contact inhibition, which can prevent excessive proliferation. Contrarily, mutated cancer cells often show uninterrupted growth that is called 'autonomous growth'. Moreover, cancer cells disperse more easily and invade the neighborhood tissue. Therefore, they have a larger growth and division rate despite limited space and nutrient supply. In the current simulation, the probability of mutation, proliferation and death is based on Eq.
(3) and several consecutive snapshots are shown in Fig. 4. After 1400 hours (approximate 58 days), cancer cells occupy the entire computational domain and its corresponding growth curve with the respect of time indicated in red color is shown in Fig. 5a. According to Eq. (8), the growth of cancer cells is influenced by α, which can be decided by α = α 0 · n 1 n 12 . To investigate the impact of α 0 on the tumor growth curve, multiple values (i.e., 0.75, 0.95, 0.98, 1) are used and the results show that growth of tumor volume slows down with the increase of α 0 value. Therefore, small variations of α 0 value facilitate our numerical model fitting experimental results.
As an example, Fig. 5b shows experimental results of tumor growth curve from the work by Durrant et al. (2015), where pancreatic cancer cells are inoculated into immunodeficient mice, where the inoculation site is subcutaneous. Implanted cancer cells are allowed to grow during two weeks before the initiation of gemcitabine drug treatment (100 µg/kg) and its growth curve is indicated by the blue curve in Fig. 5b. As a control experiment, the black curve in Fig. 5b exhibits the growth of inoculated  , where control and gem in the legend denote tumor growth without drug and with gemcitabine drug, respectively. In the simulation with curve 1, λ max and α 0 are equal 1 × 10 −3 and 0.94, respectively. However, to calibrate the model to curve 2, λ max decreases to 5.5 × 10 −4 and α falls to 0.85

Oncolytic Virotherapy
Oncolytic virotherapy has been recognized as a promising cancer treatment approach. We first develop a phenomenological model of oncolytic virotherapy in 3D, where the intratumoral injection of the virus is taken into consideration. The spread of viruses is simulated by the reaction-diffusion equation that is solved by using the FDM method. As a result, the diffusion of viruses in two different situations at time t = 50 h are compared in Fig. 6. Due to the very different doses of virus administration (Aghi and Martuza 2005;Wollmann et al. 2012), we assume that the injection is carried out during a time span of 0.5 h with a total dose of approximate 3.6 × 10 5 pfu viral particles. Figure 6a shows viruses spread with the absence of cancer cells and new breeding viruses, whereas Fig. 6b gives the distribution of viruses at time t = 50 h with the viral infection and newly generated viruses. A few isosurfaces are plotted with a color bar indicating the concentration of viruses. In contrast, viruses remain in the core of the computational domain and the highest concentration of viruses is up to 9.4 × 10 2 pfu mm 3 in Fig. 6a. This is mainly due to a slow viral diffusivity (Bajaj et al. 2001) and insufficient viruses supply. The isosurface in Fig. 6b indicates that a small amount of viruses has spread near the boundary. Note that irregular isosurface in grey color has a concentration value of slightly less than 100 pfu mm 3 , since the new burst size of viruses from a lysed cancer cell is 100 pfu mm 3 in the current model. Due to viral infection, the highest concentration of viruses in the core is 1.69 × 10 2 pfu mm 3 .
In addition, viral diffusion with the time evolution (when t = 1, 10, 50, 100 h) of each situation is shown by slices in Figs. 7 and 8, respectively. The slices are taken from the angle of a z-axis top view, which is located in the middle of the computational domain. Figure 7 presents a slow and relatively smooth diffusion phenomenon, with no viruses on the computational boundary at t = 100 h. However, viral spread in Fig. 8 is faster as a result of the supply new breeding viruses from necrotic cancer In the presence of cancer cells, viral infection ensues, viruses replicate leading to rupture of cancer cells, which then releases the viruses. The viruses are thus found also at distant locations. The isosurface in grey color has a concentration value of slightly less than 100 pfu mm 3 cells and seems to be more random, which we think to be more in line with the infection and spread of viruses in reality. Eventually, some viruses spread to the edge and are dissipated from the border of the domain. Dissipated or remaining viruses after treatment might be removed by immune cells or, in worse cases, be virulent to healthy tissue. Therefore, it is vital to assess the toxicity of the remaining viruses after treatment.
To visualize the modeling progression of oncolytic virotherapy, some consecutive snapshots are shown in Fig. 9. An extension of the model in Fig. 4, where the computational domain = 15 × 15 × 15 mm 3 has been divided into N l = 30 × 30 × 30 lattices. The domain is filled with cancer cells at time t = 0 h, see Fig. 9a. Typically, viruses are injected intratumorally if cancer occurs under the epidermis (Rehman et al. 2016), otherwise, intravenous injection is the main approach for virotherapy (Aghi and Martuza 2005). However, the intravenous injections could cause many viruses to infect other tissue outside of the tumor, or be removed by the immune system or be dissipated before reaching the cancer area. To make the problem tractable, we consider one dose of intratumoral injection, which is given at the center of the domain. Subsequently, internal cancer cells will start to get infected, indicated in black color, and subsequently die, which is indicated by the white color, see Figs. 9b, c. Local cells at lattice point i may get infected once the local concentration of viruses exceeds the threshold, which isĉ = 10 pfu mm 3 . Afterward, infected cancer cells (black color in Fig. 9) are able to die (unoccupied grid nodes in white color) and release new breeding viruses with a burst size N s = 100 pfu mm 3 . The dead lattice sites at state S = 0 are reminiscent to a wound. Then a chain reaction is triggered such that the virotherapy speeds up. Since internal lattice points are released after the death of infected cancer cells, we suppose that the normal constitutive cells around the cancer region will migrate to this area and fill the wound by proliferation, see Fig. 9d-f.
The model describes an ideal virus type with a small dose to kill cancer cells, however, the role of the viral dose remains unclear. Since some viruses, like NDV, lead to a significant therapeutic benefit at high doses, whereas other viruses do not (Aghi and Martuza 2005). However, different types of viruses have different side effects (Yu and Fang 2007). The risk could be tiny symptoms, such as flu or fever (Cripe et al. 2015), and also could be severe like fatal muscle toxicity or neurotoxicity (Russell et al. 2012). Therefore, the evaluation of residual viruses after treatment is crucially important. In our model, according to the boundary condition in Eq. (9), viruses will dissipate from the boundary to other tissues or organs. Thereby we estimate the remaining viruses in the modeled area and ignore the dissipated viruses when t = 100  Figure 10a, b show changes in total viruses and cancer volume over time in the domain, respectively. At the beginning, a total dose 0.18 × 10 5 pfu (injection rate γ = 0.5 × 10 4 pfu/h) is given (see the enlarged view in Fig. 10a), where the domain is fully occupied by uninfected cancer cells with a volume as large as 3375 mm 3 (see Fig.  10b). Once cancer cells get infected by viruses, successful viruses begin to replicate themselves until the host cancer cells rupture, which results in a significant increase in viral quantity and decrease in cancer volume. When time approaches 80 h, the number of viruses in the domain have accumulated to a peak (see Fig. 10a), whereas most cancer cells are damaged (see Fig. 10b). Furthermore, Fig. 10a shows that the Fig. 9 Consecutive snapshots of oncolytic virotherapy. The blue, red and black colors are visualized as epithelial, cancer and infected cancer cells, respectively. In addition, white color means the dead cells or unoccupied lattice points. A small scale of cancerous tissue that returns to normal tissue by cell reproduction or migration under the oncolytic virotherapy after t = 110 h  Changes in viral quantity as the evolution of time with an injection rate γ = 0.5 × 10 5 pfu/h; c Changes in viral quantity as the evolution of time with an injection rate γ = 1 × 10 5 pfu/h number of viruses gradually decreases after 80 h and this is mainly because a fraction of viruses escapes from the domain boundary. Note that there is a minor decline in the number of viruses (see the enlarged view in Fig. 10a), which may be due to the fact that the actual number of viruses present exceeds the carrying capacity of viruses. In order to investigate whether there is a maximum capacity of viruses in a limited domain, various injection rates (i.e., 1 × 10 4 , 0.5 × 10 5 , 1 × 10 5 pfu/h) are compared. The results given in Fig. 11 show that a larger viral dose leads to a greater decline in total particles after injection and do not affect the eventual result because of the maximum capacity of viruses in the computational domain. The results suggest that if a certain threshold is exceeded for the amount of injected viruses, then its temporal evolution is more or less the same. Kelly and Russell (2007) showed that immunosuppressed patients normally have a better therapeutic benefit than those who have an intact immune system in oncolytic virotherapy. However, a defective immune system would lead to a large number of viruses, which is associated with unacceptable toxicity in most cases (Russell et al.  2012). To make our model applicable to a wide range of virus species, the antiviral immune response is incorporated in the Monte Carlo simulations as one of the input variables, which we sample from a probability distribution. Therefore, Eq. (9) is revised slightly to

Monte Carlo Simulations
The injection is modeled by Dirac delta function δ(x) with source points at position x p . Moreover, β and N v define the proliferation rate and a burst size of viruses from a necrotic cancer cell. On the domain boundary ∂ , T again represents the mass transfer rate. Note that ηc represents the neutralization process by immune cells where η denotes the neutralization rate. Therefore, the antiviral immune strength is investigated by variation of the η parameter. Since the appropriate dose of a specific virus is still unclear, the total dose of viral injection is considered by varying the injection rate γ . Moreover, the infection thresholdĉ is used to evaluate the ability of viral infectivity regarding its impact on the final total particles of the remaining viruses and cancer area. To perform the Monte Carlo simulations, 5000 samples are chosen for the injection rate γ , infection thresholdĉ and immune strength η, where sampling parameters follow the normal distribution. The mean and variance of the sampling parameters have been listed in Table 3. Taking CPU time into consideration, we limit each simulation up to 50 h and then compare total particles of the remaining viruses and cancer area in the computational domain. Based on 5000 samples, Fig. 12a, b show the histograms of the total particles of the remaining viruses and cancer area, respectively, in 2D simulations with a total area of the domain of 225 mm 2 (15 mm × 15 mm). Of 5000 samples, 700 simulations end with few residual viruses, see Fig. 12a and thereby there are around 700 cases with a cancer area above 200 mm 2 in Fig. 12b. This may be caused by a combination of low injection rate, high infection rate and a strong antiviral immune response.
Since the simulation is limited to 50 h, most cases end with a large cancer area compared with the original area, which is from 100 to 200 mm 2 . Correspondingly, there is a large portion of simulations that have the remaining viral quantity ranging from 1E 4 to 3E 4 pfu at t = 50 h. To see the correlations between variables and the numerical results, several scatter plots are shown in Figs. 13, 14 and 15. The role of viral dose is tested by using the injection rate γ in Fig.13, which shows that there is no obvious correlation between the injected virus dose and the remaining viral quantity and cancer area. This is probably because of an insufficient simulation time period or the maximum capacity of viruses in a limited domain (see Fig. 11 as an illustration for this claim). Since the viral particles are not necessarily contained in the domain of computation, but are leaving the considered region as a result of diffusion (possibly also as a result of being transported by the blood vessel network). The 'loss' of particles is modeled by the Robin boundary condition. This is the reason why that once all cancer cells have been neutralized, like happens in some of the Monte Carlo simulations, the viral particles concentration will tend to zero in the computational domain, and hence a flattening behavior takes place along the 'no residual viruses axis', see Fig. 13a. The flattening behavior that is shown in Fig. 13b reflects the Monte Carlo runs for cases that the viral particles are already transported away from the domain of computation and hence the cancer cells are not reached by which the viral doses are not able to reproduce sufficiently and hence the entire region (15 × 15 × 15mm 2 = 225mm 2 ) remains full with cancer cells.
In contrast, the infection threshold, which is used to represent the ability of viral infectivity, shows a significant correlation with the remaining viral quantity and cancer  Fig. 14. The higher the threshold value, the higher the concentration of the virus is needed to infect the cancer cells, which hints at a lower ability of the viral infectivity. From Fig. 14a, the number of residual viruses decreases with increasing infection threshold since viruses with low infection ability are not able to damage cancer cells, but they can be eliminated by the antiviral immune response. Therefore, antiviral immune and insufficient newborn viruses facilitate cases with few residual viruses after 50 h. Certainly, the cancer area will not reduce significantly if the viral infection ability is weak. Viruses with a good infectivity (ĉ ≤ 5) are able to neutralize cancer cells within a period of t = 50 h. Based on Eq. (22), the term ηc(r) reflects the immune strength, therefore, the immune strength is investigated through variation of immune reduction rate η. Note that the value of c(r) is quite large and thereby the η is chosen very small (from 0 to 0.06) to guarantee a simulation with a likelihood of success. A large η denotes a strong antiviral immune response that would result in the death of most viruses. According to Fig. 15a, in the case of immunodeficiency, the residual viruses could accumulate to a large amount, while the amount falls as the immunity increases. When the antiviral response is strong like η > 0.04, viruses will be eliminated by immune cells completely in the domain. On the contrary, the cancer area has declined with the intervention of residual viruses if the antiviral immune is defective (see when η approaches to 0). However, the cancer area is more likely to be large in size when the immune response is strong, like η > 0.04. This indicates that patients with a weaker immune response may benefit from a larger reduction of the tumor size. However, at the same time, patients with a weak immune system are sensitive to have large amounts of residual viruses in their bodies.

Conclusions
Many animal-based experiments and clinical trials yielded a noticeable tumor attenuation by using oncolytic viruses (De Pace 1912;Kasuya et al. 2005). However, currently, viruses are not deemed as a useful means to stop or inhibit cancer since there is no effective way to control the virulence and retaining their replication capability in cancer cells (Kelly and Russell 2007). We have developed a cell-based model in pancreatic cancer at early stages (Chen et al. 2018b), which is subsequently extended to therapy model in (Chen et al. 2019). However, compared to classical treatments for pancreatic cancer like surgery, radiotherapy, chemotherapy, virotherapy has its own limitations, which needs further scientific assessment. In particular, the limitations include antiviral immune responses, inefficient delivery of virus as well as the poor virus spread in tumor area (Wennier et al. 2011). Therefore, more research in terms of oncolytic virotherapy is needed.
In the present study, we develop a 3D cellular automata model for oncolytic virotherapy. As we expected, the model is able to simulate cancer progression at early stages, which include the biological processes such as mutation, proliferation and death. Within 1400 h (appropriate 58 days), cancerous cells mutate from healthy somatic cells and then colonize the computational domain as big as 15 × 15 × 15 mm 3 . Certainly, the model is scalable and the speed of cancer progression can be adjusted by variation of input parameters. Therefore, different growth trends have been compared and one numerical result of our model could fit experimental results very well. Subsequently, oncolytic virotherapy is phenomenologically simulated in the same domain that is completely occupied by cancer cells. The migration and proliferation of the virus is modeled by using a reaction-diffusion equation, which is solved by a FDM method. Since viruses specifically infect and damage cancer cells, the model predicts cancer attenuation as time evolves. Eventually, normal somatic cells fill in the gap through migration and proliferation.
In addition, Monte Carlo simulations are performed in a 2D model to quantitatively investigate the correlations between several input variables and numerical results. Among 5000 samples, there are 700 simulations ending with few residual viruses and large cancer area, which dues to failed virotherapy probably as results of the extreme parameter values. The results indicate an insignificant correlation between the injection dose of viruses and simulated results (total residual viruses and cancer area), and that is probably because of an improper value range of injection rate, an insufficient simulation time period or a limited computational domain. However, we believe that this result is acceptable, since some virus species, such as NDV, show a high correlation between given doses and therapeutic benefits, whereas others do not (Aghi and Martuza 2005). Further, the viral infection threshold has a significant correlation with total amount of remaining viruses and with the final cancer area, which means that viruses with low viral infectivity likely allow a large cancer area with just few viruses left. Moreover, the anti-viral immune response presents an obvious correlation with the numerical results. Specifically, most simulations end up with relatively fewer viruses if the anti-viral reaction is strong and thereby the corresponding residual cancer area is also larger.
Due to gene mutation (i.e., RAS, TP53), the anti-viral infection ability of cancer cells is weakened, which gives the oncolytic viruses a chance (Chiocca 2002) to infect the cancer cells. Research in oncolytic viruses is not limited to cancer therapy research, also in studies that combine other treatments, such as immunotherapy (Bommareddy et al. 2018) and chemotherapy (Molnar-Kimber et al. 2002). In order to optimize the viral therapy in terms of fighting cancer, and leaving as few viral particles post-therapy as possible, further experimental studies are necessary. The required quantification in order to optimize viral therapy implies that mathematical modeling is a necessary and very helpful step.