Enhancing oncolytic virotherapy: Observations from a Voronoi Cell-Based model

Journal of Theoretical Biology 485 (2020) 110052

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/jtb

Enhancing oncolytic virotherapy: Observations from a Voronoi Cell-Based model Adrianne L Jenner a,∗, Federico Frascoli b, Adelle C.F. Coster c, Peter S. Kim d a

School of Mathematics and Statistics, University of Sydney, Sydney, NSW, Australia Department of Mathematics, Faculty of Science, Engineering and Technology, Swinburne University of Technology, Melbourne, VIC, Australia c School of Mathematics and Statistics, University of New South Wales, Sydney, NSW, Australia d School of Mathematics and Statistics, University of Sydney, Sydney, NSW, Australia b

a r t i c l e

i n f o

Article history: Received 7 November 2018 Revised 13 September 2019 Accepted 14 October 2019 Available online 15 October 2019 Keywords: Oncolytic virus Tumour Voronoi tessellation Agent-based model Cell-based model Cellular automaton Mathematical modelling

a b s t r a c t Oncolytic virotherapy is a promising cancer treatment using genetically modified viruses. Unfortunately, virus particles rapidly decay inside the body, significantly hindering their efficacy. In this article, treatment perturbations that could overcome obstacles to oncolytic virotherapy are investigated through the development of a Voronoi Cell-Based model (VCBM). The VCBM derived captures the interaction between an oncolytic virus and cancer cells in a 2-dimensional setting by using an agent-based model, where cell edges are designated by a Voronoi tessellation. Here, we investigate the sensitivity of treatment efficacy to the configuration of the treatment injections for different tumour shapes: circular, rectangular and irregular. The model predicts that multiple off-centre injections improve treatment efficacy irrespective of tumour shape. Additionally, we investigate delaying the infection of cancer cells by modifying viral particles with a substance such as alginate (a hydrogel polymer used in a range of cancer treatments). Simulations of the VCBM show that delaying the infection of cancer cells, and thus allowing more time for virus dissemination, can improve the efficacy of oncolytic virotherapy. The simulated treatment noticeably decreases the tumour size with no increase in toxicity. Improving oncolytic virotherapy in this way allows for a more effective treatment without changing its fundamental essence. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Oncolytic viruses are promising treatment agents designed to eliminate cancer cells. These viruses are genetically engineered to preferentially infect, replicate within and kill tumour cells. This genetic modification localises the infection to the tumour site and leaves nearby healthy cells unaffected. Recent clinical and experimental trials from a range of genetically modified cancer-killing viruses have shown increasing promise (Aghi and Martuza, 2005; Jebar et al., 2015; Lawler et al., 2017; Russell et al., 2012). However, oncolytic virotherapy still faces many challenges, especially those related to the rapidity of treatment decay. The rapid decay in concentration of viral particles due to clearance and dispersion at the tumour site shortens the window of effectiveness for oncolytic virotherapy. Additionally, the inability to efficiently distribute the viruses within solid tumours represents a

∗

Corresponding author. E-mail addresses: [email protected] [email protected] (P.S. Kim). https://doi.org/10.1016/j.jtbi.2019.110052 0022-5193/© 2019 Elsevier Ltd. All rights reserved.

(A.L.

Jenner),

significant barrier to the success of clinical trials (Liu et al., 2007; Parato et al., 2005). The relatively static viral distribution within a tumour is caused primarily by two factors: the non-uniformity of the tumour structure and the increase in viral clearance as a function of the number of infected tumour cells. Some studies have tried to improve the efficacy of oncolytic virotherapy by combining it with treatments to disrupt the tumour structure and reduce viral clearance, including degradation of the extracellular matrix (ECM) with Relaxin (Ganesh et al., 2007; Kim et al., 2011a) and Anti-VEGF therapies (Kottke et al., 2010). Gel-based mediums and nanoparticles are currently being investigated as a way of enhancing and controlling delivery (Yokoda et al., 2017). These modifications have the potential to delay the release of the viral particles able to cause infection of the tumour. In this study, we investigate in silico how modifying virus particles to delay viral infection, for example by coating the virus particles in alginate (a hydrogel polymer used in a range of cancer treatments), could overcome the effects of viral clearance and inhomogeneous infection and diffusion. Mean-field mathematical models of an oncolytic virus interacting with cancer cells have been shown to provide in-

2

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 1. Initial Voronoi tessellation. Healthy cells are coloured pale pink and tumour cells are bright green. The boundaries for each cell are represented by a solid line formed from a Voronoi tessellation of the lattice points, which are small dots in the centres of the cells in (a). In (b), the tessellation is overlaid with the network of connected lattice points obtained using a Delaunay triangulation. The neighbourhood of interaction is indicated in blue for one point in the lattice. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

sight into a range of treatment perturbations, see for example (Dingli et al., 2006; Jenner et al., 2018b; 2018c; Komarova and Wodarz, 2010; Wodarz, 2001). Nonetheless, for aggressive tumours, stochasticity in tumour characteristics and behaviours can be the dominant driver of cancer progression, and mean-field models are unable to fully capture this process. In this work we use an agentbased approach like that of Wodarz et al. (2012), using an offlattice framework for the formation of a tumour in a 2-D setting where cells are designated edges from a Voronoi tessellation and viruses move freely across the tessellation of cells. Researchers have used Voronoi tessellations successfully in tumour histopathological image analysis (Haroske et al., 1996), and the use of a Voronoi Cell-Based model (VCBM) allows for a versatile representation the interaction between cancer cells and virus particles. In this study, we derive a VCBM for cancer cells treated with oncolytic virotherapy and investigate treatment perturbations that could overcome obstacles of this therapy. Two key areas are investigated: the dependence of outcome on (1) the multiplicity and location of treatment injections and (2) a virus modification that allows for further intratumoural dissemination by delaying infection. Both of these key areas are investigated on three different tumour shapes: circular, rectangular and irregular. We also investigate the effect of modelling viral movement using subdiffusive as opposed to classical diffusion. Overall, our findings suggest a new method to improve treatment dissemination and antitumour effectiveness: delayed viral infection. 2. Biological model development Agent-based models can be used to simulate mechanical and physiological phenomena in cells and tissues (Van Liedekerke et al., 2015). In off-lattice agent-based models, interactions between cells are usually described by forces or potentials, and position changes in cells can be obtained by solving an equation of motion (Metzcar et al., 2019; Van Liedekerke et al., 2015). The Voronoi Cell-Based model (VCBM) we have designed is an off-lattice agent-based model that mimics tumour formation and treatment with an oncolytic virus. In the model, cells are generated from a set of points with boundaries obtained from a Voronoi tessellation. Viruses are modelled as a separate agent-based population that diffuses across the Voronoi tessellation of cells. The model evolution is driven by the virus and cell characteristics and interactions. There are four key areas to the biological assumptions for the VCBM we have designed: the underlying tumour growth biology, the oncolytic virus

life cycle, the movement of viruses, and the interaction of oncolytic viruses and tumour cells. 2.1. Underlying tumour growth model There are two main cell types used to model tumour growth in the VCBM: healthy tissue cells and tumour cells. Each cell’s position in 2-D space is defined by a point and when all points are connected they form a lattice. The Voronoi tessellation of the lattice is used to define the edges of a particular cell in the VCBM and determine the neighbourhood of interaction for a particular point in the lattice, see Fig. 1(a). The Voronoi tessellation is generated by determining the region of space where the Euclidean distance to a point is less than the distance to any other point in the lattice. The boundary of a particular cell is the line equidistant from that cell’s point and another point in the lattice, such that the set of cells generated by all the points of the lattice forms the tessellation. Voronoi cells on the boundary of the tessellation will have infinite area by definition. To avoid any interference from these boundary cells, the initially generated grid of points is always generated large enough so that the boundary cells do not influence the dynamics of the model. The advantage of the chosen lattice topology is that cells are not fixed in space, and are not inherently confined to a particular arrangement. We consider a finite domain of tumour cells and the surrounding environment and employ Dirichlet boundary conditions (i.e. cell states are fixed beyond the boundary). The initial grid of points is larger than the necessary domain for the dynamics seen in all simulations in this manuscript. Initially, the points in the lattice are arranged so that the corresponding Voronoi cells form a hexagonal tessellation (see Fig. 1(a)) analogous to other work in the literature (Buijs et al., 2004; Lobo, 2014). Many mechanisms promote cell movement in a tumour, such as pressure-driven motility, chemotactic migration or ECM-directed motion (Paul et al., 2017; Yamaguchi et al., 2005). In this model, we assume cell movement is governed by pressure-driven motility that arises from the introduction of new tumour cells through proliferation into an already crowded domain. To model this in the VCBM, the spatial relationship between points in the lattice is defined by a network of springs and modelled using Hooke’s Law. All mature cells are initially assigned an equal spring rest length with their neighbours. When a mature cell proliferates, the daughter cells are assigned a shorter spring rest length which increases over time to the mature cell spring rest length. This formalism en-

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

3

Fig. 2. Oncolytic virus life cycle. Virus particles infect either uninfected or infected tumour cells. Once inside a cell, virus particles undergo replication for a period of time. Eventually they lyse the cell, causing it to burst, and release new viral progeny that will infect other tumour cells. These virus particles either decay or are cleared by the immune system. During this process, only uninfected cells continue to proliferate.

courages movement in the lattice as the newly proliferated cells apply a time-dependent force to their neighbouring cells as they grow into mature cells. The neighbourhood of interaction for a particular Voronoi cell is defined as the neighbouring cells that share a connecting edge with it, see Fig. 1(b). These points are determined by taking a Delaunay triangulation of the lattice and finding the set of cells that are connected in the triangulation. Cells in the neighbourhood of interaction for a particular Voronoi cell are the ones that can influence the movement of that cell at each time step. A known hallmark of cancer is rapid cell proliferation; however, for tumour cells to divide there must be sufficient surrounding space and nutrients. Pressure from closely packed neighbouring cells restricts the access of oxygen and nutrients, hence cells towards the centre of an enlarging tumour receive a smaller level of nutrients than those at the periphery and tend to form a quiescent (non-proliferative) tumour cell population (Sherar et al., 1987). To account for this, we assumed that the probability of a tumour cell proliferating would be a function of its distance from the edge of the tumour. This was previously shown to be a realistic assumption for modelling the proliferation of brain tumour cells by Kansal et al. (20 0 0). Additionally, to facilitate the rapid formation of a tumour in a static tissue environment, we assume healthy cells divide much slower than tumour cells, so for the timescale and extent of the model we ignore the effects of healthy cell proliferation. 2.2. Oncolytic virus life cycle In general, the success of most oncolytic virotherapies relies on the inherent ability of viruses to infect, replicate and lyse tumour cells (Kim et al., 2006; Russell et al., 2012). More than one virus can infect a single cell (Phan and Wodarz, 2015; Syverton and Berry, 1947), but we assume that the multiplicity of infection does not affect the replication rate (Jenner et al., 2018a). A summary of the infection process of an oncolytic virus and corresponding death of a tumour cell is summarised in Fig. 2. The immune system is stimulated by the presence of virusinfected cells, initiating the clearance of extracellular virus particles (Janeway et al., 2005). In addition, viruses decay from the tu-

mour site by entering the tissue or blood vessels. In this model, the clearance of extracellular virus particles by immune cells or the decay of viruses is not modelled explicitly. Instead, the combined effects of immune clearance and viral decay are assumed proportional to the number of infected cells in a neighbourhood of a given virus particle. The neighbourhood of a given virus particle is defined as the quadrant of the tumour that the virus is located within, where horizontal and vertical axis defining the quadrants intersect at the tumour centre. 2.3. Movement of viruses The movement of virus particles within a tumour is governed by its structure (Wang and Yuan, 2006). Extracellular transport of viruses is significantly hindered by the tumour cells and the extracellular matrix, often because the size of virus particles are close to or larger than the space between fibres in the extracellular matrix (Jain, 1997; Juweid et al., 1992; Wang and Yuan, 2006). As viruses cannot move in a self-directed manner, chemotaxis or homogeneous directed motion does not generally occur (McKerrow and Salter, 2002; Wang and Yuan, 2006). Continuous spatial models for oncolytic virotherapy either do not explicitly model viral movement (Wein et al., 2003; Wodarz et al., 2012) or model viral movement by classical diffusion (Friedman et al., 2006; Mok et al., 2009). Agent-based models that consider individual virus particles have similarly used lattice random walks to model individual particle movement through a tumour (Paiva et al., 2011). However, since virus particles can struggle to diffuse due to the dense extracellular matrix and disorganised structure of the tumour cells, classical diffusion is not neccessarily the best way to model viral movement. Virus particles commonly get stuck for long periods of time, especially at their initial entry site (Kim et al., 2006), and this dynamic is analogous to a population of particles diffusing anomalously. Diffusive motion of viruses can be crucial to the outcome of oncolytic virotherapy as particles cannot disseminate and infect cells within the tumour if they get stuck at the periphery or in dense areas of the ECM. To model the possible crowding or trapping of viruses in our VCBM, we have used subdif-

4

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

fusion, a type of anomalous diffusion. Anomalous diffusion is a diffusion process whose variance scales non-linearly with time. The analogy between anomalous diffusion and diffusive motion of macromolecules has previously been discussed by Höfling and Franosch (2013). The movement of viral particles in the VCBM is modelled using a continuous-time random walk (CTRW), where a stable distribution is used to determine the waiting times between an individual particle’s consecutive movements. We choose to approximate anomalous diffusion of a population of virus particles on discrete-time intervals. While we have not proved these processes converge, the desired qualitative behaviour is exhibited by our model. There are many other factors that can influence an individual virus’s movement, for example fluid flow through the tumour and drag forces from the movement of cells. We assume that subdiffusion is able to capture the overall dynamics exhibited by viral diffusion through the tumour, including the influence of cells and fluid. We also have assumed that viruses move freely in space and are not confined to the cell lattice. 2.4. Interaction of oncolytic viruses and tumour cells To model the interaction between an oncolytic virus and a growing tumour, we consider four different sub-types of tumour cells: uninfected tumour cells, virus-infected tumour cells, dead tumour cells and empty space. Once an uninfected tumour cell is infected by a virus, it becomes a virus-infected cell. Virus-infected tumour cells do not proliferate in the model as viruses typically replicate their genomes and generate new progeny by deregulating cell-cycle checkpoint controls and modulating cell proliferation pathways (Bagga and Bouchard, 2014). Hence, the likelihood of a virus-infected cell proliferating is low and the effects are taken to be negligible in our model. Naturally, viruses are unable to distinguish between tumour cells and healthy cells and can infect both; however, the effects of viral infection of healthy cells are negligible as a result of their genetic modification. As such, we assume that healthy cells do not become infected by viruses. Additionally, as we are primarily interested in the influence of oncolytic viruses, we assume the only mode of cell death is through viral lysis. Therefore, healthy cells do not die as viral particles do not replicate within them. After a period of time, infected cells undergo lysis and become dead tumour cells. When a cancer cell dies from viral-induced cell lysis, its remnants disintegrate over a period of time. To simulate this in the model, cells can slowly take up space previously held by the dead cell. Once a dead cell has disintegrated, it turns into empty space. In our model, there are cells designated as empty space. These cells are removed from the lattice and do not contribute to the force calculation for any individual cell. These cells are only part of the model to keep the size of each living cell bounded in the visualisation. 3. Model implementation At any given time, each cell is endowed with one of five possible states: uninfected tumour cell, virus-infected tumour cell, dead tumour cell, empty space or normal healthy cell. Uninfected tumour cells can proliferate, move or become infected cells. Virusinfected tumour cells can move or die. Dead cells slowly disintegrate into empty space and do not move. Healthy cells only move over the time-scale of the investigation. Details for the implementation of viral movement, viral clearance, cell movement, cell proliferation, cell infection and cell death, are given below. Since we have assumed healthy cells can only move and do not proliferate, they are unable to regenerate and proliferate back into the empty space left by dead tumour cells in the time frame of

the simulations. This is biologically plausible if we consider for example breast cancer: once a breast tumour has been resected, patients are often left with soft tissue defects and disfigurations from the nearby tissue’s inability to regenerate itself (Stosich and Mao, 2005). 3.1. Viral movement In our VCBM, movement of virus particles through a tumour is captured by a Pearson-gamma random walk, independent of the cellular Voronoi lattice, with random waiting times between consecutive movements drawn from a heavy-tailed distribution of the form P(W > w ) ∼ w−1−α where α ∈ (0, 1). Trajectories may be simulated on a discrete-time grid by drawing a waiting time W for each particle after a single step from the stable distribution with stability parameter α using

sin(α (V + π /2 )) W = cos(V )1/α



cos(V − α (V + π /2 )) E

 1−αα

,

(1)

where V is uniformly distributed on the interval (−π /2, π /2 ) and E is exponentially distributed with unit rate parameter. For the original derivation of Eq. (1) see the work by Carnaffan and Kawai (2017) or Janicki and Weron (1993). The uncorrelated and unbiased Pearson-gamma random walk is used to simulate the movement of viruses after they have waited the appropriate amount of time. This random walk has numerous applications in the fields of physics, biology and ecology (Bartumeus et al. (2008); Codling et al. (2008); Kiefer and Weiss (1984)) and models independent steps with gamma distributed lengths and uniform orientations (Le Caër, 2010). To simulate the viral motion, we use the following algorithm. Initially, each virus particle is assigned a waiting time W, drawn from the distribution in Eq. (1). Once the virus has waited the appropriate number of time intervals, the step length of the displacement of the virus particle is drawn from a gamma distribution with mean rμ and variance rσ . The angle the virus rotates relative to its previous position is a random variable drawn from the uniform distribution [0, 2π ). Whilst other distributions could have been used, the choice of the gamma distribution was motivated, in common with other studies, by its strictly positive bell-shape, definite average and ‘tunable’ characteristics, (Frank, 2009; Jenner et al., 2018a). After each step in the virus particle’s motion, a new waiting time W is then drawn from Eq. (1). In Fig. 3, we have compared the density of viral particles after 200 hours, with and without waiting times W (Fig. 3(a), (b), and (c) respectively). It is known that the variance of the distribution of the anomalously diffusing population scales as a power-law, proportionally to tα . Smaller values of α (corresponding to heavier tails in the waiting time distribution) result in the slower spreading of viruses, while as α → 1, linear scaling of variance with time is recovered as the regularity of long trapping events decreases (Carnaffan and Kawai, 2017; Janicki and Weron, 1993). As a result, as the value of α is increased in Fig. 3(a) and (b), the spread in the histograms is noticeably increased. To provide more insight into how subdiffusive viral motion differs from a continuous random walk, in Fig. 3(d) we have plotted the mean-squared displacement (MSD), corresponding to the density of viral particles after 200 hours in Fig. 3(a), (b), and (c). It is evident from Fig. 3(d) that viral movement without waiting times results in a linear MSD as a function of time. Whereas, viral particles with waiting times between consecutive movement generated from Eq. (1) with stability parameter of α = 0.6 and α = 0.8 result in a nonlinear MSD over time. Initially the displacement of viral particles from the initial seeding location increases quickly for subdiffusive movement. Then, as time goes on, we see a decrease

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

5

Fig. 3. Spatial histograms of the distribution of 30 0 0 virus particles initially located at the origin after 200 hours, where particles are diffusing with waiting times given by Eq. 1 with (a) α = 0.6, (b) α = 0.8 and (c) no waiting times. The corresponding mean-squared displacement (MSD) of the virus particles in (a), (b) and (c) are plotted in (d), and the approximate slope, , of the MSD at 200 hours has been labelled.

in how rapidly the mean displacement of the virus particles from their initial position is increasing. As estimation of α is difficult Burnecki et al. (2015)) and data for the displacement of viral particles within a tumour is not available. We, therefore, fixed α = 0.6 based on Ernst et al.’s calculation for particle movement in mucus (Ernst et al., 2017). For more thorough insight into the time averaging of CTRW with a broad distribution of waiting times, see the work by Neusius et al. (2009). 3.2. Viral clearance and decay To simulate clearance of the virus due to immune stimulation, the number of infected cells in a region close to each virus is considered. To do this, the tumour is split into four sections which are formed from vertical and horizontal axes intersecting at the centre of the tumour. The number of infected cells in each section, Ii , is calculated at each time step. The probability that an individual particle is cleared or decays into the tissue or blood vessels pc is given by

Ii pc = pδ , IT

(2)

 where IT is the total number of infected cells, i.e. IT = 4i=1 Ii , and pδ is the clearance and decay constant. This constant pδ is the dimensionless probability of a cell decaying in a given time step t + t. Modelling viral clearance in this way simulates the increased likelihood of virus clearance or decay with an increase in the proportion of infected cells in a region close of the virus (Filley and Dey, 2017; Janeway et al., 2005). Additionally, it accounts for the division of resources when multiple infected sites exist in different quadrants. For all simulations, the clearance and decay con-

Fig. 4. Schematic illustrating the connection between points k, j and i in the lattice at a fixed time t. Springs connect points in the lattice and the movement of each point depends on the force derived from Hooke’s Law, assuming that motion is overdamped due to strong friction. The spring rest length between points k and j is sk,j , and points k and i is sk,j , and sk,j < sk,i .

stant was pδ = 1/3, fixing the maximum probability of clearance or decay in a single time step as 1/3. 3.3. Cell movement Cell movement in the tumour microenvironment is assumed to be a result of neighbouring cell pressure caused by the proliferation of tumour cells. To model this, the position of each cell (except dead cells and empty cells) is updated by calculating the effective displacement of the cell’s lattice point using Hooke’s Law. The force exerted on the kth cell is modelled as the force exerted on that cell’s point in the lattice. Force is modelled as a

6

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 5. Schematic to illustrate how cell-to-cell adhesion is assumed to be negligible after the cells have reached a distance greater than s + al .

network of damped springs connecting the kth point to its neighbouring points. The spring connecting the kth and jth point has a rest length sk,j (t), which can vary over time t. All points in the lattice are connected in this way and the spring rest lengths between points can be unique. In Fig. 4, we present an example of a possible set-up between three points k, j and i. In this example the spring rest length sk,j is shorter than the spring rest length sk,i . This variable spring rest length is used to model the growth and decay of newly proliferated cells and cells in the process of disintegrating. Following the implementation in (Meineke et al., 2001; Murray et al., 2009; Osborne et al., 2017), the displacement of the kth point on the lattice is given by

mk

 d2rk = Fk,I j + FkV , 2 dt

(3)

where mk is the mass of the kth point, rk is its spatial position, Fk,I j is the interaction force between a pair of neighbouring points, FV k

is the viscous force acting on the kth point, and the sum is taken over neighbouring points to k in the lattice. The total interaction force FkI (t ) acting on the kth point at time t is equal to the sum of all forces from the springs of all points i connected to k:



Fk,I j = μ

j=k

 rk,i (t )   s (t ) − ||rk,i (t )|| , ||rk,i (t )|| k,i

(4)

∀i

where μ is the spring constant, rk,i (t ) is the vector from the kth to the ith point at time t, sk,i is the spring rest length from the kth to the ith point at time t and ||rk,i (t )|| is the L2-norm of the vector rk,i (t ), see Fig. 4. Eq. 3 can then be simplified using two key assumptions. The first is that the viscous force Fk,V j , i.e. point-point, point-medium and point-matrix interactions, can be modelled by assuming that the drag on the kth point is independent of the springs and is proportional to its velocity, with constant of proportionality η. Secondly, the points are assumed to be in a relatively dissipative environment, so point motion can be approximated as being overdamped due to strong friction. Hence

d2r mk 2k dt

∼ 0.

(5)

Thus

FkI = −FKV = ηvk ,

(6)

where vk is the velocity of the kth point. Approximating this velocity over a small time interval t, we obtain

r (t + t ) − r (t ) FkI ≈ η . t

rk (t + t ) = rk (t ) + +λ

1 Fk (t )t = rk (t )

η

 rk,i (t )   s (t ) − ||rk,i (t )|| , ||rk,i (t )|| k,i

(8)

∀i

j=k

FkI (t ) =

Fig. 6. Schematic for the probability of a particular cell proliferating given a particular distance d from the edge of the tumour, see Eq. 9. The maximum radial distance for which proliferation occurs, dmax , separates the tumour into proliferating and non-proliferating sections, with the cells inside the shaded circle having a distance greater than dmax from the edge, and hence being unable to proliferate.

(7)

Thus the effective displacement of the kth point within a small time interval t in the overdamped limit is

where rk (t ) is the position of the kth point in the lattice at time t and η is the damping constant. Cell mobility is described by the ratio λ = μ/η, which is known to influence the velocity of the relaxation process (Meineke et al., 2001). Many diverse molecules take part in the cell-to-cell and cellto-ECM adhesion in cancer (Okegawa et al., 2004). Here, adhesion effects between neighbouring cells are modelled simply using a linear force and cut-off distance al . When the Euclidean distance between points on the lattice of neighbouring cells is longer than sk,l + al , no interaction takes place, see Fig. 5. 3.4. Cell proliferation To model cell proliferation, a cell’s distance to the nutrient source and local spatial limitations are considered. We assume the distance d from a cell to the nutrient source, is the distance between the cell and its closest peripheral tumour cell. The maximum radial distance that still allows for a cell to obtain nutrients from its surroundings is accounted for by dmax . If d > dmax , then the cell does not proliferate, essentially becoming a quiescent cell. The probability of a cell dividing based on the nutrients it receives is



pd = p0 1 −

d dmax



,

(9)

where p0 is a proliferation constant. Note that p0 is dimensionless as pd is the dimensionless probability of a cell proliferating in a given time step t + t. Fig. 6 is a schematic illustrating how dmax segregates the tumour into a rim of proliferating and nonproliferating cells. Additionally, to account for the spatial limitation on proliferation, tumour cells only divide if there is least rmin space between a tumour cell and any cell in its neighbourhood of interaction. This formulation was used by Kansal et al. (20 0 0) to model the proliferation of brain tumour cells. A similar probability calculation for cell proliferation can also be seen in the work on invasive tumour growth by Jiao and Torquato (2011). Only uninfected cells proliferate in our framework and the division is symmetric.

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

7

Table 1 Parameters in the model and their meanings. Cell parameters gage Time steps taken for a daughter cell to grow to adult size Probability constant for proliferation p0 d Distance from a tumour cell to the nearest cell on the tumour edge Radial distance that nutrient reaches by diffusing from the tumour edge dmax Minimum distance between neighbouring cells for proliferation to occur rmin Age a cell needs to reach before it can proliferate again page Distance spring length of a dead cell decreases at each time step dfrac Time steps taken for dead cell to disintegrate dage Cell motility parameters s Spring rest length μ Spring constant η Damping constant Adhesion distance between two Voronoi cell points al Virus parameters Time steps taken from infection to when the infected cell bursts iage Probability of infection occurring pi Mean distance of viral movement per time step rμ Standard deviation for virus displacement rσ ν Number of new virus agents created through lysis

cells are placed at a random angle of orientation from the original position of the mother cell at a random distance less than or equal to s/gage apart. Once a cell has proliferated, it takes page time steps before the daughter cells will proliferate again, including the gage time steps for the daughter cells to mature and reposition. 3.5. Cell infection

Fig. 7. Schematic for cell motility after proliferation. The proliferation of cell k into two new cells k and l results in a spring rest length sk,l = s/gage between daughter cells, which increases over time to s (the mature cell spring rest length).

In each time-step, all uninfected and infected tumour cells are checked to see whether they will be infected. If there is at least one virus agent within the perimeter of a cell, then the cell becomes infected with probability pi . If there are ς T virus particles within the perimeter of a cell, then a random number of virus particles ς , drawn from a uniform distribution of the number of viruses able to infect that cell, i.e. ς ∼ U(0, ς T ), infects the cell. These virus agents are then removed from the extracellular virus population and the tumour cell becomes an infected cell. 3.6. Cell death Once inside a tumour cell, viral particles undergo replication for iage time steps, after which the cell will burst and release ν new virus particles. These new particles are placed randomly within the perimeter of the cell that burst. The cell that bursts then becomes a dead cell. We assume that dead cells take dage time steps to disintegrate. To simulate disintegration, at each time increment, the spring rest lengths of a dead cell to each of its neighbours, sk,i , reduces by dfrac = sk,i /dage .

Fig. 8. Illustration of the spring rest length for two daughter cells after cell division. After a cell divides, the two daughter cells k and l will have a spring rest length sk,l = s/gage . This spring rest length, sk,l , increases as a function of time since division until gage time steps have passed, at which point the cell has the ’adult’ cell rest length s. The rest length between the daughter cells is then fixed until one of the cells proliferates after page time steps have passed. See Fig. 7 for a schematic detailing the spring rest lengths after cell proliferation.

If a cell proliferates, the framework also allows the addition and rearrangement of lattice points (and the associated VCBM) as a cell divides into two daughter cells. When a uninfected cell at position k divides, a new lattice point l is created, so that k and l are now the points associated with the two daughter cells, see Fig. 7. To simulate the enlargement and repositioning of the daughter cells (Ghaffarizadeh et al., 2018; Mumenthaler et al., 2013), the resting spring length of the connection between k and l is taken to be linearly increasing from s/gage to the mature resting spring length, s, over a time gage as indicated in Fig. 8. Note the two daughter

4. Parameter optimisation and sensitivity All major parameters in the model are collated in Table 1. The parameters relating to cell state characteristics were optimised using time-series measurements for the growth of cervical cancer SKOV3 cells in vivo (Kim et al., 2011a). The model was assumed to be updated on a time step of 4 h for the data optimisation and all subsequent numerical simulations. A summary of the viral characteristic parameters obtained from the literature is also presented in Section 4.2. 4.1. Optimising cell state characteristics To verify the robustness and predictive capability of our model, we investigated the growth of cervical cancer in mice. The volume of cervical cancer SK-OV3 cells has been measured in three mice

8

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052 Table 2 Cell parameter values. The calibrated parameters are obtained by tuning the model to the experimental measurements of SK-OV3 cell growth seen in Fig. 9 and the fixed parameters are taken from the initial hexagaonal tessellation and the literature.

Fig. 9. Model calibration for in vivo cervical cancer SK-OV3 cell growth. Individual mouse tumour cell numbers recorded by Kim et al. (2011a) are plotted as grey circles. Overlayed in light blue dotted lines are 15 model simulations of tumour growth with the mean of these simulations in black. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Calibrated parameters

SK-OV3 data

gage (hour) p0 dMAX ( μm) page (hour) Fixed parameters rmin (μm) s (μm) μ(μg/hour2 ) η (μg/hour−1 ) al (μm) dage (hour)

12 0.7 4 28 17.28 18.50 0.01 0.133 0.15 8

Table 3 Virus characteristic parameter values. These parameters were fixed to values found in the literature. Parameter

Value

Source

pi iage (hour) ν (no. of virus particles) rμ ( μm) rσ ( μm)

0.9 24 5 0.05 0.1

Jenner et al. (2018a) Jenner et al. (2018a) Jenner et al. (2018c) Mok et al. (2009) Mok et al. (2009)

fore given in terms of hours. We confirmed that the dynamics do no change for different time steps (1 hour, 2 hours and 4 hours). 4.2. Viral characteristic parameters

Fig. 10. Schematic for invasive cell proliferation. The designated invasive cell (blue) that is on the edge of the tumour (green cells) proliferates in the direction of its healthy cell neighbour closest to the vector joining the centre of the tumour and the invasive cell. One daughter cell is designated an invasive cell and kills and takes the position of the healthy cell, and the other daughter (now a non-invasive tumour cell) cell remains in the original position of the invasive cell. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

over time by Kim et al. (2011a). SK-OV3 cells have an average diameter of 14.1 μm (Chen et al., 2011) and, interpolating the calculation of Monte (Del Monte, 2009), this equates to approximately 4.8 × 108 cells per mm3 . Due to computational and geometric constraints in our model, rather than representing a single cell in our 2D setting, each Voronoi cell was taken to represent the average characteristics of 1010 SK-OV3 cells, and the data scaled accordingly by a factor of 1010 . The VCBM was then used to determine the growth of SK-OV3 cells, with the trend showing good correspondence with the data, see Fig. 9 and Table 2. We approximated the time taken for a cell to disintegrate after lysis as dage = 8 hours based on experiments of Monastra et al. (1996) which measured the nucleus disintegration in lysed cells over a period of 8–20 h. The parameters rmin and s are taken from the original hexagonal lattice set up (Fig. 1). Hooke’s-law related parameters μ, η and λ are taken from Meineke et al. (2001), but scaled to our time step of 4 hours, and al is arbitrary. Note that the volume of cervical cancer SK-OV3 cells was measured over 25 days, and the parameters returned were therefore a function of days. In the Results section, we have plotted the number of Voronoi cells and considered the tumour growth as a function of hours. The parameter values presented in Table 2 are there-

The remaining parameter values in the model were approximated from literature. The probability pi of infection was approximated by the rate of infection reported by Jenner et al. (2018a), assuming that the infection rate can be modelled by an exponential probability distribution. The decay rate was approximated using the average duration of cell lysis of 24 hours (Jenner et al., 2018a). The number of new viruses created, ν , was approximated from Jenner et al. (2018c), increased slightly, because we are considering individual virus agents rather than a mean-field approximation to the viral dynamics. The mean and variance for viral step length were approximated from the diffusion coefficient of HSV particles reported at 5 × 10−10 cm2 s−1 by Mok et al. (2009). See Table 3 for a summary. 4.3. Model simulation Cancer cell proliferation, spatial limitations and obstructions influence the shape of a tumour. To understand the sensitivity of the treatment to tumour shape, we investigate the outcome of virotherapy on three generic tumour shapes: circular (Fig. 11(a)), rectangular (Fig. 11(c)), and irregular (Fig. 11(e)). To generate each shape, the same basic underlying VCBM rules have been used with some additions in the case of rectangular and irregular tumour growth. Circular tumours can be generated directly using the model described in the previous section If there are no spatial limitations for proliferating cancer cells in vivo, a roughly circular tumour will form, similar to the cross-section of a hanging drop tumour spheroid (Weiswald et al., 2015), see Fig. 11(a). Note the apparent difference between the exponential growth seen in Fig. 9 and the more linear growth seen in Fig. 11(a). This difference is due to the slower growth of the tumour on the time scale of hours as opposed to days, and also due to the smaller number of cells.

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

9

Fig. 11. Representative tumour shapes considered for treatment with an oncolytic virus (a)-(b) circular, (c)-(d) rectangular, and (e)-(f) irregular. The corresponding number of tumour cells as a function of time from 12 simulations has been plotted for each shape. Typical evolution plots for each of these shapes can be found in Figs. A.28, Fig. A.29 and A.30 (Appendix). Note the different scale for (f) due to the extremely fast growth of the irregular tumours.

When there is an obstruction above and below an initial seeding of tumour cells in vivo, a rectangular tumour shape will form. This obstruction can be considered as stiffer stromal tissue, similar to that seen in breast ductal carcinomas in situ (DCIS) (Ghaffarizadeh et al., 2018; Kim et al., 2011b), which have approximately rectangular cross sections. These tumours form in a rectangular shape due to spatial limitations above and below the tumour. To encourage rectangular tumour formation in our VCBM, a horizontal impenetrable boundary is positioned above and below the

initial grouping of tumour cells. This was simulated by placing a dense horizontal line of points into the lattice (not plotted) at the position of the impenetrable boundary, and then requiring that all healthy cells above and below the two horizontal boundaries are unable to move. In certain cancers, cells on the periphery of a tumour can become invasive cells, allowing them to degrade the ECM and remove nearby cells, forming irregular branches (Jiao and Torquato, 2011). To investigate the effect of oncolytic virotherapy on a tumour of

10

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 12. To compare the differences in outcome measures between the circular, rectangular and irregular control tumour growths in Fig. 11 at (a) 160hours and (b) 280hours, a one-way ANOVA was performed followed by a Tukey–Kramer posthoc test for significance. The results of these tests are in Table B.4 (Appendix). Horizontal lines linking groups indicate statistically significant differences between the groups, p < 0.05. The horizontal positions of the box lines for each sample data are the 25th, 50th and 75th percentiles, the top and bottom of the angled lines represent the upper and lower 95% confidence limit for each sample median and the whiskers represent the maximum and minimal value of the sample data. Note that only the circular and rectangular tumour growths are compared in (b) since the irregular tumour was simulated only until 160 hours.

this geometry, we generate what we call an irregular tumour, which is a circular tumour that has grown invasive branches. To generate an irregular tumour, we designate a collection of tumour cells on the tumour periphery as invasive cells. These cells are inherently tumour cells and will undergo infection and movement in the same way as tumour cells. Their additional invasive properties come into effect when the invasive cell undergoes proliferation. When an invasive cell satisfies the probability to proliferate designated by Eq. (9) and there is at least one neighbouring healthy cell, the invasive cell divides with one invasive daughter cell killing and occupying the position of the healthy cell and the other daughter cell occupying the original position (and reverting to a non-invasive tumour cell). Furthermore, in the case of more than one healthy neighbour, the invasive tumour cell does not split in a random direction, but rather kills and occupies the position of the healthy cell neighbour which is the furthest from the tumour centre. In this way, invasive cells move in a direction that maximises the nutrient and oxygen supply for the tumour (Jiao and Torquato, 2011). In our model we assume nutrients are flowing in from the boundary of the domain, so these invasive cells move in a direction away from the tumour. See Fig. 10 for a summary of the generation of an irregular tumour. Representative time evolutions for each of the tumour shapes (circular, rectangular and irregular) can be found in Figs. A.28, A.29 and A.30 (Appendix). While the shapes generated in this article have not been directly matched to experimental tumour images, they represent shapes occurring in tumour formation as discussed above. Over time, healthy cells are surrounded by the tumour cell population and this is reminiscent of the true biological scenario, since healthy cells are regularly found within tumours (Park et al., 20 0 0). The one-way analysis of variance (ANOVA) with posthoc Tukey–Kramer analysis, see Fig. 12, shows that the circular and rectangular tumours are approximately the same size at 160 hours, with the irregular tumour being significantly larger (p < 0.05). By 280 hours, a difference can be seen between the circular and rectangular tumours, with the circular tumour becoming significantly larger than the rectangular. Using the aforementioned rules for oncolytic viruses and corresponding parameter values in Tables 1 and 3, a representative

model evolution is shown for a circular tumour with viral treatment in Fig. 13 (See also Fig. A.31)". All simulations presented in the following results section use the same VCBM virus rules, unless specified otherwise. Since the size of the adenovirus is approximately 90–100 nm (Appert et al., 2012), this means that based on the size of an SK-OV3 cell, an adenovirus is 0.6%-0.7% of the cell’s size. Therefore we have assumed virus transport mechanisms are not influenced by whether they are in a cell-occupied area or an area with no cells. 5. Results: Simulating alternative treatment protocols With the advancement of oncolytic viruses to clinical development, delivery is a major barrier of success. Traditionally, viral therapy is administered by either intravenous or intratumoral injection (Wang and Yuan, 2006). Irrespective of the administration protocol, host immunity, tumour microenvironment and abnormal vascularity contribute to inefficient virus delivery (Yokoda et al., 2017). Two major therapy perturbations are examined in the following subsections: the configuration of the viral injections entry sites and the effects of delaying the infection of cancer cells by viral particles. These modifications to the current therapy are examined in detail on the three tumour shapes generated. The impact of modelling viral dissemination with subdiffusion is then investigated. In this study, we have chosen to investigate oncolytic virotherapy on small tumours, aiming to improve the diffusive viral properties of this therapy for early-stage cancers. 5.1. Dependence of treatment outcome on entry site configuration Intratumoural injections are regularly used to administer oncolytic virotherapy; however, the possible dependence of therapy outcome on the position of initial treatment injection has not been investigated systematically. Using our model, we determined whether there is an optimal injection configuration for the three tumour shapes: Circular (Fig. 14), rectangular (Fig. 16) and irregular (Fig. 18). The varying injection configurations considered are represented pictorially by enumerated and coloured virus shapes with the resulting number of tumour cells as a function of time also shown. Across all injection configurations, the total dosage was the

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

11

Fig. 13. Representative evolution of the VCBM model in a circular configuration. The four snapshots above represent equal intervals of model dynamics. Pale pink cells represent healthy cells, dark green cells represent tumour cells, light green cells represent quiescent cells, bright pink cells represent infected tumour cells and grey cells represent dead cells with empty space shaded in light grey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

same. In the case of more than one injection, the dosage was split evenly amongst the injections. Increasing the number of injection sites and distance from the tumour centre improved the overall treatment efficacy for circular tumours, see Fig. 14. A comparison of the tumour size at 280 h was undertaken using a one-way ANOVA and Tukey–Kramer posthoc test, see Fig. 15 and Table B.5. The least effective modalities were the single injections (A, B and C), and those that were centrally applied (D and G). These were not statistically significantly different to each other, and the tumour size under these treatment protocols was significantly (p < 0.05) larger than that for the other multiple off-centre injections, whether applied uniformly (E and F) or not (H and I). The application of multiple off-centre injections was significantly more effective when they were applied uniformly (E and F) compared to non-uniformly (H and I), although the radius at which they were applied did not appear to have a significant effect at 280 hours, with both E and F not significantly different to each other, and similarly H and I not being significantly different. It is possible that the multiple off-centre injections performed better

than the central applications as the virus did not need to spread as far to reach the growing boundary the tumour. Similarly, off-centre single injections allowed large sections of the growing boundary to increase without treatment. In rectangular tumours, single injections had a more varied effect on tumour growth than circular tumours, see Fig. 16. The rectangular tumours were constrained along symmetric impenetrable boundaries which forced the tumour to grow only in the direction of the major axis, i.e. the axis through the centre of the tumour (horizontal in Fig. 16). Injections were given either along the constrained-edges or along the major axis in three positions: at the tumour periphery, at the tumour centre or halfway between the periphery and the centre. Different treatment protocols were compared using an ANOVA and Tukey–Kramer posthoc test of tumour size at 280 hours (significance for p < 0.05), Fig. 17 and Table B.6 (Appendix). The untreated tumour was significantly larger than those resulting from all treatment protocols. Single injections along the major axis were found to be significantly better (smaller tumour size) if placed more centrally (A and B, not significantly different to each other) compared to a single injec-

12

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 14. The effects of the multiplicity and configuration of the injection site on a circular tumour. The initial injection configurations considered are represented in Figs. 14(a), (c) and (e) by the coloured regions, depicting viral particles. The corresponding total number of tumour cells over time for each injection type is plotted in Figs. 14(b), (d) and (f), respectively for 12 simulations. For reference, untreated tumour growth is plotted for 12 simulations in grey. Note quiescent cells are not plotted.

tion closer to the tumour periphery (C, worse than A and B). This is most likely a result of intratumoural single injections accessing more of the tumour than those at the edge. There were examples of individual tumour growths after a single injection that stabilised (see Fig. 16(a)); however, these occurred rarely. Splitting the treatment into two major axis injections showed a significant improvement with two mid-axis injections being more effective than a single injection similarly placed on one side alone (D significantly less than B). Similarly, double injections at both growing faces of the tumour (E) were significantly better than a single injection at only one edge (C). If the treatment was applied at the constrained edge the outcome depended on axial position. At the centre of the tumour, a single injection on the major axis

was more effective than splitting it into two symmetric injections at the constrained edge (A significantly less than F). Mid-axis, a single injection at the constrained edge was worse than a single major axis injection (B was significantly less than H); however, splitting the treatment into two injections at the constrained edge was more effective than a single axial application (G was significantly less than B). Interestingly, the most effective treatments were D and G, both mid-axial applications. The application of two mid major-axis injections (D) was significantly more effective than the dual injections on one side (and at the constrained edge) (G), perhaps suggesting that targeting both growing edges of the tumour is advantageous. Additionally keeping the treatment a little back from

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

13

Fig. 15. To compare the differences in outcome measures for the circular tumour growths at 280hours for injection configurations A to I in Fig. 14, a one-way ANOVA was performed followed by a Tukey–Kramer posthoc test for significance. The results of these tests are in Table B.5 (Appendix). There were no significant differences between the tumour size in the case of single injections and centrally applied injections (A, B, C, D and G). Similarly the tumour size for off-centre, uniform application of injections (E and F) was not significantly different, as was the case for off-centre, non-uniform injections (H and I). All other differences were statistically significant with p < 0.05.

the growing edge may allow the virus to be more effective than applications at the periphery, due to the interplay of diffusion to uninfected cells and the time for viral replication. Profile G however only controls one of the proliferating directions of the tumour whereas D targets both growing domains. Different treatment protocols in irregular tumours (Fig. 18) were again assessed by ANOVA and Tukey–Kramer posthoc testing of the tumour size at 160 h (Fig. 19 and Table B.7 in the Appendix). All outcomes were statistically different (p < 0.05) except B and C. The best result by far in the treatment of an irregular tumour occurs if every branch off-shoot of the tumour is identified and targeted, splitting the virus into multiple injections at each growing branch tip (D). If the tips cannot be targeted, but the branches identified, then the next best (although a far worse outcome) is to target the root of the off-shoot at the periphery of the tumour bulk (E). Interestingly this performs better than injections part-way along the branches (F). In reality, it is unlikely that every offshoot can be identified and located. Targeting all the virus application as a single injection in the centre resulted in the worst outcome, indicating that distribution of the initial application is warranted. Considering splitting the treatment into three injections, if some growing branches can be identified, the treatment could be applied at two growing offshoots and also the tumour bulk (B). Unsurprisingly this performs significantly worse than the profiles targeting each off-shoot (E, F and D); however, it is equivalent to three injections in the tumour bulk (C). Both B and C performed better than three injections at the periphery of the tumour bulk (A). 5.2. Treatment with delayed initial viral-infection While optimising the injection configuration can improve the treatment outcome, the efficacy of virotherapy can also be improved by altering the virus’s delivery mechanism. Traditionally, oncolytic viruses are administered by either intravenous or intratumoural injection; however, there are a number of novel approaches for enhancing oncolytic virotherapy delivery, such as gelbased mediums and nanoparticles (Yokoda et al., 2017). Gel-based mediums have been used to improve virotherapy, chemotherapy and immunotherapy, see for examples (Choi et al., 2013; Jung et al., 2017; 2014; Oh et al., 2017; Smith et al., 1999; Wade et al., 2017). In each of these examples, the gels are

loaded with the therapeutic agents and are injected adjacent to the tumour providing sustained therapy release and therapeutic efficacy. Nanoparticles have also been investigated for their use as a viral DNA and RNA delivery system (Riley and Vermerris, 2017). The advantage of these vectors is they have a decreased immune response and may be fuctionalised for different purposes through their exterior design (Riley and Vermerris, 2017). Nanoparticle physical properties, such as encapsulation and coating with biodegradable polymers, can be used to provide controlled viral release and diminish systemic infectivity, maintaining an elevated local concentration (Yokoda et al., 2017). For example, Tseng et al. (2018) designed a carrier loaded with virus that oxidised from contact with lactate (present in the tumour microenvironement), releasing the viral contents. Another example of controlled viral delivery was by Choi et al. (2015), who developed nanoparticles deployed with a synthesised pH-sensitive polymer that released the viral particle in an acidic environment. Polymer-based nanomaterials have also drawn interest as potential delivery vectors, such as Polyethyleneglycol (PEG), which have the ability to shield nanoparticles from the extracellular environment preventing subsequent clearance (Riley and Vermerris, 2017). In all of these therapies, the treatment activity was extended and controlled resulting in an improvement of therapy effectiveness. The current investigations point to the possibility that modifying viral vectors to delay the initial infection time would allow for further diffusion through the tumour bulk prior to the initial infection and activation of clearance, resulting in a more effective therapy. While this idea is yet to be tested experimentally, there are a few examples of how it could be implemented either using gel-based coating or nanoparticle technology. Alginate, a naturally occurring biopolymer, has several unique properties that have enabled it to be used for delivery of a variety of biological agents, including viruses (Choi et al., 2013; Qurrat-ul Ain et al., 2003). The ability to encapsulate viral particles in alginate microbeads has been tested as a vaccine delivery system (Kwok et al., 1989; Wee et al., 1995). Choi et al. (2013) showed that the biological activity of viral particles loaded in alginate gel is prolonged compared with naked virus (Choi et al., 2013). The microenvironment around the gel-encapsulated virus also provides protection from clearance by the immune system over an extended time period (Choi et al., 2013; Muruve et al., 1999; Ruzek et al., 2002).

14

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 16. The effects of the multiplicity and configuration of the injection site on a rectangular tumour. The initial injection configurations considered are represented in Fig. 16(a), (c) and (e) by the coloured regions, depicting viral particles. The corresponding total number of tumour cells over time for each injection type is plotted in Fig. 16(b), (d) and (f), respectively for 12 simulations. For reference, untreated tumour growth is plotted for 12 simulations in grey. Note quiescent cells are not plotted.

Our hypothesis is treatment efficacy can be improved by developing a mechanism that can delay infection of the virus. To investigate, we used the VCBM framework and considered the effects of a fixed delay before the initial virus infection. We additionally assumed that the modification allows the virus to diffuse, but does not promote immune clearance. As the size of an adenovirus is 0.6%-0.7% the size of an SK-OV3 cell, we have assumed that the modification will not influence the size significantly enough to warrant a change in the transport properties and rates. The effects of delayed viral treatment on circular, Fig. 21, rectangular, Fig. 23, and invasive tumours, Fig. 25, were investigated. In each of these figures, the predicted number of tumour cells under treatment with the original oncolytic virus is overlaid with that of

the delayed virus for the various delay times. Four different delays of 40, 52, 60 and 80 h were simulated. These delays could be achieved by controlling pore size, degradation rate and release kinetics of alginate (Gombotz and Wee, 1998). Delaying infection allows the treatment to disseminate further into the tumour, see Figs. 20, 22 and 24. In Figs. 22 and 24, it is clear that the variance of the position of the viral particles initially is smaller than that of the delayed viral particles before their initial infection, after their corresponding wait periods. The treatment was administered in three injections on the periphery for circular tumours (profile F, Fig. 14(c)), two injections on the short ends for rectangular tumours (profile E, Fig. 16(c)) and three injections on the tumour bulk periphery for irreg-

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

15

Fig. 17. To compare the differences in outcome measures for the rectangular tumour growths at 280hours for injection configurations A to H in Fig. 16, a one-way ANOVA was performed followed by a Tukey–Kramer posthoc test for significance. The results of these tests are in Table B.6 (Appendix). All comparisons were statistically significant with p < 0.05, apart from groups (A, B & E), (C, F & H) and (D & G). The horizontal positions of the box lines for each sample data are the 25th, 50th and 75th percentiles, the top and bottom of the angled lines represent the upper and lower 95% confidence limit for each sample median and the whiskers represent the maximum and minimal value of the sample data.

Fig. 18. The effects of the multiplicity and configuration of the injection site on an irregular tumour. The initial injection configurations considered are represented in Fig. 18(a) and (c) by the coloured regions, depicting viral particles. The corresponding total number of tumour cells over time for each injection type is plotted in Fig. 18(b) and (d), respectively for 12 simulations. Note the shorter timescale, and greater growth compared to Figs. 14 and 16. For reference, untreated tumour growth is plotted for 12 simulations in grey. Note quiescent cells are not plotted.

ular tumours (profile A, Fig. 18(a)). For circular and rectangular tumours, Figs. 21 and 23, the modified delayed virus resulted in lower tumour cell numbers at 280 h than that of the non-delayed virus. In Fig. 26 one-way ANOVA analyses for the number of tumour cells for systems with different delays demonstrate there were statistically significant differences between these simulations. Each injection had the same amount of

virus particles, irrespective of the tumour shape or delay length. In Fig. 26, a comparison of the treatment effectiveness at discrete time points is presented. The number of tumour cells as a function of the virus delay at 10 0, 20 0 and 280 h have been summarised for circular and rectangular tumours, with differences assessed using one-way ANOVA and posthoc Tukey–Kramer analyses.

16

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 19. To compare the differences in outcome measures for the irregular tumour growths at 160hours for injection configurations A to G in Fig. 18, a one-way ANOVA was performed followed by a Tukey–Kramer posthoc test for significance. The results of these tests are in Table B.7. All injection configurations were statistically significant (p < 0.05) except for B & C. The horizontal positions of the box lines for each sample data are the 25th, 50th and 75th percentiles, the top and bottom of the angled lines represent the upper and lower 95% confidence limit for each sample median and the whiskers represent the maximum and minimal value of the sample data.

Fig. 20. Virus diffusion for (a) non-modified and (b) modified virus particles at 52 hours. A comparison of the two cases shows that the region within which the virus particles have diffused is similar in both cases; however, there are fewer virus particles in the non-modified case due to immune clearance.

In the circular tumours, the delay of the onset of viral infection initially allowed the tumour cell numbers to increase rapidly compared to the tumours undergoing non-delayed viral treatment, Fig. 21. At the onset of the infection a dramatic drop in tumour cell numbers was observed. For shorter delay times, Fig. 21(a) and (b), the tumour numbers were similar under the delayed and nondelayed treatments following the onset time. For treatments with longer delays, Fig. 21(c) and (d), the tumour numbers under the delayed treatment dropped below that of the non-delayed treatment at the onset time and the trend of the subsequent increase was on average below that of the non-delayed treatment with significant differences at 10 0, 20 0 and 280 h (p < 0.05, ANOVA with Tukey–Kramer posthoc Fig. 26 and Table B.8). This illustrates that while the addition of a delay in the viral infection allows the tumour to grow to a larger size, the treatment’s effectiveness is increased due to its ability to disseminate further into the tumour, see Fig. 22. It is clear from the distribution of viral particles at the end of each delay period, that the ability to disseminate further aids in the reduction of the overall tumour growth.

For rectangular tumours, delaying the onset of viral infection results in tumour growths with equivalent slopes to the non-delay treatment and different tumour sizes at 280 h (p < 0.05, ANOVA with Tukey–Kramer in Fig. 26 and Table B.9 in the Appendix). In the case of rectangular tumours, Fig. 23, delaying the viral infection means that once the virus particles infect, the size of the tumour is lower at that point than the size of the tumour under treatment with non-delayed virus. So while the tumour growths have the same rate, the size of the tumour is smaller for the delayed viral treatment. Again, from the distribution of viral particles at the end of each delay period, see Fig. 24, the increased dissemination of the virus aids in the reduction of the overall tumour growth. In the case of irregularly shaped tumours, Fig. 25, whilst the dramatic drop in tumour cell numbers was observed upon the onset time in the delayed treatments, this did not decrease the numbers below that of the non-delayed treatments, and the long term trends were the same. This may be caused by the rapid growth of the irregular tumours, which made the delayed approach ineffective. Therefore, while modifying viral particles to delay their in-

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

17

Fig. 21. Size of circular tumour under treatment with a delayed oncolytic virus. The number of tumour cells predicted by our model over time is plotted under treatment with the original oncolytic virus (light blue) and with the delayed virus (dark blue) applied using profile F, Fig. 14. Twelve model simulations were considered for each case. The wait times for the delayed virus cases were (a) 40, (b) 52, (c) 60 and (d) 80 hours. Note the dramatic drop in tumour cell numbers upon initial viral infection in each case. In (e), the ratio of the number of extracellular virus particles to uninfected tumour cells is plotted for each treatment as a function of time, and in (f) the number of infected cells for each treatment is plotted as a function of time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

fection could significantly improve therapy, the model simulations improve therapy, the model simulations showed that it depends heavily on tumour shape. It is clear from Fig. 25 that the modified treatments in irregular tumours have no significant effect without further analysis. On average, the delayed virus was more effective for circular and rectangular tumours than the non-delayed treatment, or the control case when no treatment was administered, irrespective of the length of the delay, see Fig. 26. It was clear that at 100 hours, the choice of delay had an effect on the number of tumour cells,

with a delay of 52 hours resulting in the smallest average tumour cell number for both circular and rectangular tumours; however, at 200 hours, there was less sensitivity to the initial delay. This is confirmed by the ANOVA and posthoc Tukey–Kramer test, where p < 0.05 between the different delays at 100 hours but not at 200 or 280 hours. This illustrates that in the short term, the length of the delay before the initial viral infection can play a significant role in the size of the tumour. To determine whether this result solely relied on the fact that the tumour size at the start of a delayed virus’s infection is larger, we simulated the two types of treatment

18

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 22. The corresponding distribution of viral particles for Fig. 21 has been plotted at the moment before they start infecting for the case of (a) no delay, (b) 40 hour delay, (c) 60 hour delay and (d) 80 hour delay.

on the same size tumour at the start of their infection time, see Fig. C.33 (Appendix). We found that since the delayed virus had disseminated further, the treatment performed much better than the non-delayed virus, given its inability to disseminate as far.

5.3. Impact of subdiffusive viral movement on therapy effectiveness In the VCBM, viral movement follows a subdiffusive process to account for the hindrance of extracellular transport in the tumour microenvironment. As current models for oncolytic virotherapy in the literature have only considered classical diffusion, we investigated the dependence of the previous section’s results on the underlying model for viral movement. All previous simulations were obtained using waiting times between successive individual viral movements drawn from the distribution W, Eq. (1), with stability parameter α = 0.6. This value of α was chosen based on the work by Ernst et al. (2017) and is classified as strong subdiffusion (Burnecki et al., 2015). To examine how the treatment is affected by the type of subdiffusion, the model was simulated with weak subdiffusion, i.e. α = 0.8. Additionally, to assess how treatment performs under more classical diffusion, three Pearson-gamma random walks were considered by removing the waiting times between successive viral movements, see Fig. 27. For the strong and weak subdiffusion models, the step length for each virus particle was drawn from a gamma distribution with mean rμ = 0.05 and variance rσ = 0.1. As demonstrated in Fig. 3(d), the long-term slopes of the MSD at 200 hours, , for the two subdiffusion models were quite different with  = 0.06 for strong subdiffusion and  = 0.11 for weak subdiffusion.

To compare the effect of subdiffusion or a Pearson random walk on the efficacy of treatment, the values of rμ and rσ were chosen so that the long-term slope of the MSD for the Pearson random walk would match the long-term slope of the MSD for either the strong or weak subdiffusion. The equivalent Pearson random walk for strong subdiffusion was obtained by drawing step lengths from the gamma distribution with rμ = 0.12 and rσ = 0.05. Comparatively, the Pearson random walk for weak subdiffusion was obtained by drawing step lengths from the gamma distribution with rμ = 0.25 and rσ = 0.15. We also chose to simulate a Pearson random walk with the same rμ and rσ as for the weak and strong subdiffusion models, this had a long-term MSD slope of  = 0.76, see Fig. 3(d). The effect of the weak and strong subdiffusion and the comparative Pearson random walk on the efficacy of treatment for circular and rectangular tumours is plotted in Fig. 27. The model was simulated for the injection protocol F on circular tumours, Fig. 14, and injection protocol E on rectangular tumours, Fig. 16. The treatment was simulated with and without the delay modification of 40 hours, see Figs. 21 and 23. The mean and standard deviation of the number of tumour cells for 12 simulations for each of the different viral movement models is plotted. From Fig. 27, it is clear that the different viral movement models affect the number of tumour cells and, therefore, the treatment efficacy. The Pearson-gamma random walks with low MSD slope, , result in a practically ineffective treatment. Increasing the stability parameter from α = 0.6 to α = 0.8, improves the treatment, demonstrating that weak subdiffusion results in a more effective therapy than strong subdiffusion. The most effective treatment is achieved in the absence of waiting times and with a large MSD slope of  = 0.76. Overall, qualitatively similar results are seen for

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

19

Fig. 23. Size of rectangular tumour under treatment with a delayed oncolytic virus. The number of tumour cells predicted by our model over time is plotted under treatment with the original oncolytic virus (light blue) and with the delayed virus (dark blue) applied using profile E, Fig. 16. Twelve model simulations were considered for each case. The wait times for the delayed virus cases were (a) 40, (b) 52, (c) 60, and (d) 80 hours. Note the dramatic drop in tumour cell numbers upon initial viral infection. In (e), the ratio of the number of extracellular virus particles to uninfected tumour cells is plotted for each treatment as a function of time, and in (f) the number of infected cells for each treatment is plotted as a function of time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

both the circular and rectangular tumours, suggesting that there is a relationship between the model for viral movement and the effect of virotherapy. Future work will be dedicated to experimentally quantifying this relationship. 6. Discussion The rapid clearance of viral particles is a major obstacle in the effectiveness of oncolytic virotherapy. Viral particles are cleared by the immune system, reducing both the number of particles acting and the window of time within which the treatment persists. In this article, we developed a Voronoi Cell-Based

model (VCBM) for the interaction between a growing tumour and an oncolytic virus treatment and investigated ways to optimise the treatment protocol. We have found that by optimising the injection site configuration and modifying the viruses to delay their infection, it may be possible to improve the efficacy of this therapy with a particular focus on small early stage tumours. The VCBM in this article is malleable and translatable to other cancer types both by extension to 3D and tuning to other growth rates and formation characteristics. Additionally, the VCBM could serve as an excellent visualisation tool for biologists and clinicians to use to pinpoint the depth of the treatment spread and

20

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 24. The corresponding distributions of viral particles for Fig. 23 before they can initially infect for the case of (a) no delay, (b) 40 hour delay, (c) 60 hour delay and (d) 80hour delay.

its effectiveness, as exemplified by Fig. 14. One difference of this model to previous spatial models for oncolytic virotherapy is that viruses are assumed to follow subdiffusive motion generated from a Pearson-gamma random walk with waiting times between successive movements drawn from a stable distribution. This particular modelling has been chosen specifically to capture the trapping phenomena and the inability of viruses to diffuse through the tumour microenvironment, Kim et al. (2006), something that is not captured through the use of classical diffusion (Friedman et al., 2006; Mok et al., 2009) or simple lattice random walks (Paiva et al., 2011). There are two primary protocols for administering an oncolytic virus, either intratumorally or intravenously. When treatment is administered intravenously, it is challenging to predict where the treatment will enter a tumour and it is usually at multiple sites on the tumour periphery. Alternatively, when treatment is administered intratumourally, it is possible to designate the entry site location to some extent. In Figs. 14, 16 and 18, we show how the application location is a crucial determinant of treatment efficacy for circular, rectangular and irregular tumours. The optimal injection profile for each shape varied significantly. In general, multiple injections performed better than a single injection; however, this was dependent on injection locations. To determine which injection profiles were significantly different, one-way ANOVA and posthoc Tukey–Kramer analyses were conducted on the end-point of each simulations, see Figs. 15, 17 and 19 and Tables B.5, B.6 and B.7). For a circular tumour, Fig. 14, single injections of oncolytic virus particles gave rise to the highest average tumour size over time.

The lack of response to the single injection protocol could be due to single injections promoting an increased multiplicity of infection with a subsequent enhancement of viral trapping that limits the breadth of spread irrespectively of the location on the tumour. In contrast to the single injections, the mean and standard deviation of a circular tumour size at 280 hours was noticeably lower for three injections along the radius of the tumour or at the periphery (Fig. 15). Increasing the injections allows for an improved diffusivity and virus to cell contact rate. The best performing injection protocol on circular tumours for a window of 50 hours was three radial and rotationally symmetric injections. After 50 hours, the efficacy of central tumour injections was overtaken by injections on the tumour periphery. We hypothesise that this effect may be due to injections on the periphery of circular tumours restricting cell growth on the boundaries, where the growth rate is usually the fastest. The optimal treatment injection configuration for rectangular tumours, Fig. 16, was not dissimilar to that of circular tumours: on average, increasing the injection multiplicity improved the treatment efficacy. However, there were some cases where a single injection close to the edge of a tumour caused a reduction of tumour growth similar to multiple injections, see B and E in Fig. 16(b) and (d) and Fig. 17. This stochasticity of tumour response could be explained by the fact that the tumour’s primary growth occurs in a horizontal direction. Since ductal carcinomas in situ (DCIS) have approximately rectangular cross sections, this could suggest that treating these cancers with virotherapy would result in a wide variety of unpredictable responses, depending on the location and multiplicity of injection.

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

21

Fig. 25. Invasive tumour under treatment with a delayed oncolytic virus. The number of tumour cells predicted by our model over time is plotted under treatment with the original oncolytic virus (light blue) and with the delayed virus (dark blue) for profile A, Fig. 16. Twelve model simulations were considered for each case. The wait time for the delayed virus cases was (a) 40, (b) 52, (c) 60, and (d) 80 hours. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Two intratumoural injections positioned mid-way along the major axis of the rectangular tumour (profile D, Fig. 16) significantly reduced the tumour size. However, administering treatment in this way could be difficult. An alternative is to administer two injections mid-way along the contrained edges above and below the tumour (profile H, Fig. 16). In real tumours, this could be a simpler way of administering treatment as the tumour itself does not need to be penetrated. Statistical comparisons (Fig. 17 and Table B.6) showed that the effects of these injection profiles are similar and thus the simpler application may be preferable. For an irregular tumour shape, it is shown in Fig. 18 that the most successful treatment outcome in the first 60 hours was obtained when the treatment was administered at the periphery of the tumour bulk as opposed to the invasive branches. The multitude of injections significantly affects this outcome: when ten injections are given at the periphery of the tumour, the treatment did considerably better than when three injections were administered at the periphery (see Fig. 19). Additionally, comparing this option to one single intratumoural injection, it fared the worst out of all the possible profiles. After 60 hours, the optimal injection configuration was the one that eliminated all invasive cells on the tumour spines. We hypothesise that by eradicating the invasive cells the tumour growth is reduced significantly so that, while treatment takes longer to be effective, the overall tumour growth rate is reduced. It is also interesting to note that there is a pronounced oscillation in tumour cell numbers soon after injection for both the circular and rectangular tumours, see Figs. 14, 16, 21 and 23. There are two reasons why we see this oscillation in the model. As this behaviour is observed predominantly after intratumour injections

(profile E and H in Fig. 14, and profile D in Fig. 16), the oscillations may be caused by the increase in free-space within the tumour allowing cells to proliferate that were previously too confined. The second reason is, in the case of the delayed virus, a significant amount of cells are infected once the virus delay has passed, and these cells lyse simultaneously, allowing nutrients to access cells, that previously may have had a low probability of proliferating. While optimising the treatment injection configuration can help to improve the standard of the current treatment, there is more we can do to increase the efficacy. To tackle the diffusivity obstacle presented by the tumour microenvironment, oncolytic viral particles can be modified to delay their infection for a specific period of time. Delaying the onset of viral infection within a tumour allows for further infiltration of the tumour prior to the onset of immune clearance. This delay could be achieved through alginate coating of the viral vectors, or developing gel-based delivery systems, or nanocarriers. For circular tumours, the delayed virus caused a significant reduction in the number of tumour cells, see Figs. 21 and 26. At 100 hours, the length of the delay had a noticeable effect on the tumour cell number with a delay of 80 hours resulting in a significantly less effective treatment than all other delays. While early-on a delay of 52 hours resulted in the lowest average tumour number, by 280 hours the statistical analysis revealed there was no significant difference between the lengths of the delays. We suggest that this reflects a relationship between the length of the delay and the time it takes for the effects of the delayed treatment to stabilise. The success of delayed virus on circular tumours can be attributed to the diffusive properties of the virus. Since the infection of the tumour cells is initially delayed, the virus is able to

22

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. 26. Comparison of treatment effectiveness at discrete time points for different infection delays. Individual, mean and standard deviation measurements corresponding to the number of tumour cells from the simulations Figs. 21 and 23 after (a) & (d) 100hours (b) & (e) 200hours and (c) & (f) 280 hours is plotted for circular tumours (a)-(c) and rectangular tumours (d)-(e). To compare the differences in outcome measures for the treatments, one-way ANOVA with Tukey–Kramer posthoc analyses were undertaken. The results of these tests can be found in Tables B.8 and B.9 (Appendix). Horizontal lines indicate statistically significant differences between the groups, p < 0.05. The control result without treatment (from Fig. 11(b) and (d)) is included in each figure for reference (but not in the statistical analysis). Note the vertical axis break in the top row and the different vertical scales.

disseminate further into the tumour before it infects the cells. Additionally, due to the coating of the virus, the immune system is unable to clear the virus during the period of time before its initial infection. Therefore, since the mean-squared displacement of the virus will be larger for a virus that has had a delayed infection, more virus undergoes the first round of infection. As a result, more secondary virus is created from the first round of infection than the non-delayed virus’s first round of infection, see Figs. 21(e) and 23(e). We also hypothesise that the initially large number of infected cells in the delayed virus treatment compared to the non-delayed virus (Figs. 21(f) and 23(f)) is one of the main drivers of the success of this treatment. These large numbers of infected cells result

in a large number of virus particles (Figs. 21(e) and 23(e)) which allows the treatment to reduce the tumour growth more significantly. Interestingly, for all simulated treatments the infected cell number and ratio of viruses to tumour cells appears to stabilise to a similar population size or growth pattern, suggesting that it is the initial response of treatment due to the delay that drives the overall difference seen in tumour growth. For rectangular tumours, the modified virus, irrespective of the delay, had a considerable effect on reducing the tumour burden, see Figs. 23 and 26. By 200 hours, each delayed virus reduced the tumour numbers to approximately the same level and at 280 hours, there was no statistically significant difference between the tumour sizes. The only major difference in the delay was how

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

23

Fig. 27. Viral movement modelled using weak and strong subdiffusion (SD) or a Pearson-gamma random walk (RW). Subdiffusive motion was generated with waiting times drawn from W, Eq. (1), with stability parameter α = 0.6 for strong subdiffusion and α = 0.8 for weak subdiffusion. These had a long-term MSD slope of  = 0.06 and  = 0.11 respectively (see Fig. 3(d)). Two Pearson random walks were generated to have equivalent  to weak and strong subdiffusion by manipulating the mean rμ and variance rσ of the step-lengths (i.e. rμ = 0.12, rσ = 0.05 for strong subdiffusion and rμ = 0.25, rσ = 0.15 for weak subdiffusion). A Pearson random walk with the equivalent rμ and rσ as the subdiffusion models was also simulated and this had corresponding  = 0.76. The number of tumour cells predicted by our model over time is plotted for a circular tumour (a) - (b) and rectangular tumour (c) - (d) under treatment with the different viral movement models. Each solid line and shaded region corresponds to the mean and standard deviation of the number of tumour cells from 12 simulations for a virus without any delayed infection (a) & (c) and with a delayed infection of 40 hours (b) & (d).

many cells were initially removed by the virus, but even this was not significant. There appeared to be no optimal delay in this case as all worked just as well as each other. This may be due to rectangular tumours having a maximum tumour growth rate, since they are bounded above and below by an impenetrable boundary. The overall growth rate of rectangular tumours was less affected by the delayed viral infection treatment compared to that of the circular tumours. Since the virus diffuses radially from the initial injection and rectangular tumours grow primarily horizontally, the effectiveness of the delayed virus, due to its ability to disseminate further into the tumour, could be reduced accordingly. In the case of irregular tumour, the modified delayed virus was no more effective than the original non-delayed form. In a few cases, the non-delayed form resulted in lower tumour sizes than the delayed form. We hypothesise this is due to the aggressive nature of the tumour shape. Invasive tumours are able to grow away from the treatment sites in a multitude of ways, evading the viruses swiftly and efficiently. In general, the chosen delay does not contribute significantly to the outcome for circular and rectangular tumours, with all delays resulting in a smaller tumour cell numbers on average than the treatments with no delay or for a tumour growing without treatment. For circular and rectangular tumours, modifying viral particles to delay their infection can produce a notable advantage to

therapy; however, for irregular tumours the benefit remains unclear. This model is presented in 2-dimensions and with the addition of a third dimension, the drop in tumour cell numbers could be significantly larger. Also, administering additional injections of the delayed-treatment could be worth investigating, especially in respect to the timing of the subsequent injections. Currently, there is limited data for individual viral movement through tumour tissue. As such, we assumed that the stability parameter for the subdiffusion model could be approximated by α = 0.6, similar to particles diffusing through mucus (Ernst et al., 2017). From numerical simulations of different values for the stability parameter in Fig. 3, we can see that increasing α to α = 0.8 reduces the trapping of the virus at the initial injection site and increases the mean-squared displacement (MSD) of the viral particles. This then improved the overall efficacy of the oncolytic virus for circular and rectangular tumours (Fig. 27), suggesting that a more effective treatment can be obtained by increasing virus dissemination. To quantify the impact of subdiffusion on the efficacy of treatment, the VCBM was simulated with viral movement following Pearson-gamma random walks with equivalent MSD slopes to strong and weak subdiffusion (i.e.  = 0.06 and  = 0.11 respectively), see Fig. 27. Treatment with virus movement following these random walks was less effective than their subdiffusion counter-

24

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

parts. Most likely this is due to the linear time-dependence of the Pearson random walks MSD, see Fig. 3(d). While the slopes of the MSD at 200 hours for the Pearson random walks and weak/strong subdiffusion were equivalent, the MSD curves were not. As subdiffusion has a non-linear time-dependent MSD, the initial MSD would be larger for the virus population moving with subdiffusion as opposed to the Pearson random walk with the same MSD slope at 200 hours. Overall, this suggests that the initial displacement of the virus may be what determines the long-term outcome. Simulating viral movement following a Pearson-gamma random walk with a long-term MSD slope of  = 0.76 noticeably improves treatment (Fig. 27). This  value is considerably larger than the previous subdiffusion and Pearson-gamma random walks considered. It is surprising, therefore, that treatment does not improve more significantly as the virus is able to diffuse further. This suggests that increasing viral dissemination through the tumour, for example by using a substance that can break down the ECM like Relaxin (Ganesh et al., 2007; Kim et al., 2011a), may have a limiting impact on improving the treatment efficacy. As such, further treatment modifications such as additional injections may be needed to improve treatment. It is worth noting that the results of this article are dependent on the underlying dynamics of the model. To investigate the possible validity of these results on a more heterogeneous sampling of tumour types and viral treatment modalities, it would be worth considering more extensively the different viral dissemination profiles, tumour growth dynamics and viral clearance rates. It has also previously been reported that the injection timing and volume can play a role in the effectiveness of oncolytic virotherapy, (Dingli et al., 2006; Jenner et al., 2018c), and this would be worth investigating in future work. In the next stage of this work, it also would be worth quantifying how the immune system dynamics may influence the viral clearance and whether the delayed virus is truly more effective (Fig. A5).

further advance vector therapy by maximizing safety, efficacy and duration of transgene expression. We have shown that the treatment injection site configuration plays a significant role in the overall treatment outcome and found optimal injection locations as a function of tumour shape. We believe that the suggested optimisation of treatment injection configuration is worth investigating further and verifying experimentally as it could prove to be a better treatment protocol. In this article, we discuss how a modified delayed-virus is able to reduce the tumour size in circular and rectangular tumours, but has limited to no effect on irregular tumours. This novel addition to therapy could prove to be a simple yet effective way to improve the efficacy of oncolytic virotherapy for some tumour types. As we continue to apply our Voronoi cellular automaton and agent-based modelling to oncolytic virotherapy, we plan to improve the current model to account for other more complicated cancer therapies. Extending the current model to three-dimensions could allow for other tumour characteristics to be included in the model, such as vascularisation, angiogenesis and necrosis. The extension to such models does, however, increase the computational cost, as both the number of cells and time period over which they require simulation would necessarily be much larger than the current investigation. In the last decade, the immune system and its interaction with cancer therapies and cancer cells have become a topic of much discussion. The current study provides the framework for models of realistic tumour geometries and for the explicit inclusion of immune cells enabling a deeper understanding of immune involvement in oncolytic virotherapy. Acknowledgments The authors received support through an Australian Postgraduate Award (ALJ) and an Australian Research Council Discovery Project DP180101512 (ACFC, FF and PSK). We would also like to thank CO Yun for providing the published data in Fig. 9.

7. Conclusion Appendix A. Typical time evolution of the model The theoretical perspective presented in this paper provides valuable insight into the biological process behind cancer formation and treatment using oncolytic viruses. The development of an optimised oncolytic virus and an effective delivery system would

Typical evolutions of the model for three tumour shapes: circular, rectangular and irregular are presented in Figs. A.28, A.29 and A.32 respectively.

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. A1. Typical evolution of the cellular automaton for a circular tumour.

25

26

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. A2. Typical evolution of the cellular automaton for a rectangular tumour.

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. A3. Typical evolution of the cellular automaton for an irregular tumour.

27

28

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. A4. Typical evolution of the cellular automaton for a circular tumour with a non-delayed viral treatment.

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

29

Fig. A5. Typical evolution of the cellular automaton for a circular tumour with a delayed viral treatment.

Appendix B. Results of the ANOVA and posthoc Tukey–Kramer tests

Appendix C. Impact of initial tumour size on delayed viruses effectiveness

In the following tables are the results of the one-way ANOVA and posthoc Tukey–Kramer tests for the simulations in Section 5.1 and 5.2. The F-statistic and p-value for the one-way ANOVA are placed at the top of the table and the individual pvalues from the posthoc Tukey–Kramer of the grouped comparison are in the rows and columns of the table. The threshold for statistical significance is p < 0.05.

To determine whether the results in Figs. 21 and 23 might depend on the fact that the tumour size at the start of a delayed virus’s infection is larger, we simulated the two types of treatment on the same size tumour at the start of their infection time, see Fig. C.33. This was done by allowing the delayed treatment to disseminate as the tumour grew for the first 40 hours, while the nondelayed treatment was stagnated as the tumour grew. In this way, at 40 hours when both virus treatments started taking effect, the tumour was at the same size. It is clear from the simulations in Fig. C.33 that tumour size is most influenced by the delay allowing for treatment dissemination and not the initial size of the tumour at the time of infection.

Table B1 One-way ANOVA and posthoc Tukey–Kramer test results for the control tumour growths (circular, rectangular and irregular) at 160 hours, see Fig. 13. At 280 hours, the circular and rectangular tumour growths gave F-statistic = 801.02 and p-value = 10−19 . F statistic =28371.64, p-value = 10−54

Circular Rectangular

Rectangular

Irregular

0.6

10−10 10−10

30

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052 Table B2 One-way ANOVA and posthoc Tukey Kramer test results for the injection configurations A to I applied to circular tumours in Fig. 15 at 280 hours. F statistic = 143.7, p-value = 10−56 A Control A B C D E F G H

10

B −6

C −6

10 0.7251

D −6

10 0.3040 0.9997

E −6

10 0.9999 0.3363 0.0828

F −6

10 10−7 10−7 10−7 10−7

G −6

10 10−7 10−7 10−7 10−7 0.9982

H −6

10 0.9942 0.9973 0.8881 0.8734 10−7 10−7

I −6

10−6 10−7 10−7 10−7 10−7 10−5 10−6 10−7 0.9999

10 10−7 10−7 10−7 10−7 10−6 10−7 10−7

Table B3 One-way ANOVA and posthoc Tukey Kramer test results for the injection configurations A to H applied to circular tumours in Fig. 17 at 280 hours. F statistic = 54.68, p-value = 10−33

Control A B C D E F G

A

B

C

D

E

F

G

H

10−8

10−8 0.9955

10−5 10−4 0.0065

10−8 10−8 10−8 10−8

10−8 0.9890 0.6894 10−6 10−7

10−5 10−4 0.0133 1 10−8 10−5

10−8 10−4 10−5 0.9996 0.2542 0.0140 10−8

10−4 10−5 10−4 10−8 10−8 10−7 0.9961 10−8

Table B4 One-way ANOVA and posthoc Tukey Kramer test results for the injection configurations A to G applied to circular tumours in Fig. 19 at 280 hours. F statistic = 849.59, p-value = 10−78

Control A B C D E F

A

B

C

D

E

F

G

10−8

10−8 10−6

10−8 0.0019 0.8276

10−8 10−8 10−8 10−8

10−8 10−8 10−8 10−8 10−8

10−8 10−8 0.0037 10−5 10−8 10−8

10−8 10−7 10−8 10−8 10−8 10−8 10−8

Table B5 One-way ANOVA and posthoc Tukey–Kramer test results for the delayed virus treatment of circular tumours, see Figs. 21 and 26. The analysis of numerical simulations was conducted at 100 hours, 200 hours and 280 hours.

Table B6 One-way ANOVA and posthoc Tukey–Kramer test results for the delayed virus treatment of rectangular tumours, see Fig. 23 and 26. The analysis of numerical simulations was conducted at 100 hours, 200 hours and 280 hours. At 100 hours F statistic = 32.75, p-value = 10−14

No delay 40 hours 52 hours 60 hours

40 hours

52 hours

60 hours

80 hours

10−6

10−8 10−4

10−8 10−4 1

10−6 1 10−4 10−4

At 200 hours F statistic = 14.24, p-value = 10−8

At 100 hours 40 hours

F statistic = 21.11, p-value = 10−10

No delay 40 hours 52 hours 60 hours

40 hours

52 hours

60 hours

80 hours

10−6

10−8 0.6398

10−8 0.4209 0.9966

0.2313 0.0068 10−5 10−5

No delay 40 hours 52 hours 60 hours

10

−7

40 hours

No delay 40 hours 52 hours 60 hours

60 hours

80 hours

0.001

10−6 0.6378

10−6 0.4202 0.9967

10−5 0.6878 1 0.9927

At 280 hours F statistic = 8.47, p-value = 10−5 40 hours No delay 40 hours 52 hours 60 hours

0.0172

52 hours 0.0036 0.9836

60 hours −5

10 0.2410 0.5334

80 hours 10−4 0.6148 0.8987 0.9638

10 0.9927

60 hours −6

10 0.9830 1

80 hours 10−7 1 0.9927 0.9830

F statistic = 8.22, p-value = 10−5

F statistic = 12.07, p-value = 10−7 52 hours

−6

At 280 hours

At 200 hours

40 hours

52 hours

No delay 40 hours 52 hours 60 hours

10

−4

52 hours −4

10 0.9996

60 hours

80 hours

0.002 0.9402 0.9810

10−4 1 0.9996 0.9402

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Fig. C1. Circular tumour size under treatment with a delayed infecting oncolytic virus. The number of tumour cells predicted by our model over time is plotted under treatment with the original oncolytic virus (light blue) and with the delayed virus (dark blue) with a wait time of 40 hours. To investigate the dependence of the results in Figs. 21, and 23 on the initial tumour size, both viruses were only able to infect from 40 hours, but the non-delayed virus was unable to disseminate. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

References Aghi, M., Martuza, R.L., 2005. Oncolytic viral therapies-the clinical experience. Oncogene 24 (52), 7802–7816. Qurrat-ul Ain, S., GK Khuller, S., Garg, S.K., 2003. Alginate-based oral drug delivery system for tuberculosis: pharmacokinetics and therapeutic effects. J. Antimicrobial Chemother. 51 (4), 931–938. Appert, J., Raynor, P.C., Abin, M., Chander, Y., Guarino, H., Goyal, S.M., Zuo, Z., Ge, S., Kuehn, T.H., 2012. Influence of suspending liquid, impactor type, and substrate on size-selective sampling of ms2 and adenovirus aerosols. Aerosol Sci. Technol. 46 (3), 249–257. Bagga, S., Bouchard, M.J., 2014. Cell cycle regulation during viral infection. Cell Cycle Control: Mech. Protocols 165–227. Bartumeus, F., Catalan, J., Viswanathan, G., Raposo, E., Da Luz, M., 2008. The influence of turning angles on the success of non-oriented animal searches. J. Theoret. Biol. 252 (1), 43–55. Buijs, J.O.d., Musters, M., Verrips, T., Post, J.A., Braam, B., Van Riel, N., 2004. Mathematical modeling of vascular endothelial layer maintenance: the role of endothelial cell division, progenitor cell homing, and telomere shortening. Am. J. Physiol.-Heart Circul. Physiol. 287 (6), H2651–H2658. Burnecki, K., Kepten, E., Garini, Y., Sikora, G., Weron, A., 2015. Estimating the anomalous diffusion exponent for single particle tracking data with measurement errors-an alternative approach. Scientif. Rep. 5, 11306. Carnaffan, S., Kawai, R., 2017. Solving multidimensional fractional Fokker–Planck equations via unbiased density formulas for anomalous diffusion processes. SIAM J. Scientif. Comput. 39 (5), B886–B915. Chen, Y.-C., Lou, X., Ingram, P., Yoon, E., 2011. Single cell migration chip using hydrodynamic cell positioning. In: International Conference on MicroTAS, pp. 1409–1411. Choi, J.-W., Jung, S.-J., Kasala, D., Hwang, J.K., Hu, J., Bae, Y.H., Yun, C.-O., 2015. Ph-sensitive oncolytic adenovirus hybrid targeting acidic tumor microenvironment and angiogenesis. J. Control. Release 205, 134–143. Choi, J.-W., Kang, E., Kwon, O.-J., Yun, T., Park, H.-K., Kim, P.-H., Kim, S., Kim, J., Yun, C.-O., 2013. Local sustained delivery of oncolytic adenovirus with injectable alginate gel for cancer virotherapy. Gene Therapy 20 (9), 880. Codling, E.A., Plank, M.J., Benhamou, S., 2008. Random walk models in biology. J. R. Soc. Interf. 5 (25), 813–834. Del Monte, U., 2009. Does the cell number 109 still really fit one gram of tumor tissue? Cell Cycle 8 (3), 505–506. ´ K., Russell, S.J., Bajzer, Ž., 2006. Mathematical modelDingli, D., Cascino, M.D., Josic, ing of cancer radiovirotherapy. Math. Biosci. 199 (1), 55–78. Ernst, M., John, T., Guenther, M., Wagner, C., Schaefer, U.F., Lehr, C.-M., 2017. A model for the transient subdiffusive behavior of particles in mucus. Biophys. J. 112 (1), 172–179. Filley, A.C., Dey, M., 2017. Immune system, friend or foe of oncolytic virotherapy? Front. Oncol. 7, 106. Frank, S.A., 2009. The common patterns of nature. J. Evol. Biol. 22 (8), 1563–1585. Friedman, A., Tian, J.P., Fulci, G., Chiocca, E.A., Wang, J., 2006. Glioma virotherapy: effects of innate immune suppression and increased viral replication capacity. Cancer Res. 66 (4), 2314–2319. Ganesh, S., Edick, M.G., Idamakanti, N., Abramova, M., VanRoey, M., Robinson, M., Yun, C.-O., Jooss, K., 2007. Relaxin-expressing, fiber chimeric oncolytic adenovirus prolongs survival of tumor-bearing mice. Cancer Res. 67 (9), 4399–4407. Ghaffarizadeh, A., Heiland, R., Friedman, S.H., Mumenthaler, S.M., Macklin, P., 2018.

31

Physicell: an open source physics-based cell simulator for 3-d multicellular systems. PLoS Comput. Biol. 14 (2), e1005991. Gombotz, W.R., Wee, S., 1998. Protein release from alginate matrices. Adv. Drug Delivery Rev. 31 (3), 267–285. Haroske, G., Dimmer, V., Steindorf, D., Schilling, U., Theissig, F., Kunze, K.D., 1996. Cellular sociology of proliferating tumor cells in invasive ductal breast cancer.. Anal. Quant. Cytol. Histol. 18 (3), 191–198. Höfling, F., Franosch, T., 2013. Anomalous transport in the crowded world of biological cells. Rep. Progr. Phys. 76 (4), 046602. Jain, R.K., 1997. Delivery of molecular and cellular medicine to solid tumors. Adv. Drug Delivery Rev. 26 (2–3), 71–90. Janeway, C.A., Travers, P., Walport, M., Shlomchik, M., 2005. Immunobiology: the immune system in health and disease. Garland Science New York. Janicki, A., Weron, A., 1993. Simulation and Chaotic Behavior of Alpha-stable Stochastic Processes, 178. CRC Press. Jebar, A.H., Errington-Mais, F., Vile, R.G., Selby, P.J., Melcher, A.A., Griffin, S., 2015. Progress in clinical oncolytic virus-based therapy for hepatocellular carcinoma. J. General Virol. 96 (7), 1533–1550. Jenner, A., Yun, C.-O., Yoon, A., Kim, P.S., Coster, A.C., 2018a. Modelling heterogeneity in viral-tumour dynamics: the effects of gene-attenuation on viral characteristics. J. Theoret. Biol.. Jenner, A.L., Coster, A.C., Kim, P.S., Frascoli, F., 2018b. Treating cancerous cells with viruses: insights from a minimal model for oncolytic virotherapy. Lett. Biomath. 1–20. Jenner, A.L., Yun, C.-O., Kim, P.S., Coster, A.C., 2018c. Mathematical modelling of the interaction between cancer cells and an oncolytic virus: insights into the effects of treatment protocols. Bull. Math. Biol. 1–15. Jiao, Y., Torquato, S., 2011. Emergent behaviors from a cellular automaton model for invasive tumor growth in heterogeneous microenvironments. PLoS Comput. Biol. 7 (12), e1002314. Jung, B.-K., Oh, E., Hong, J., Lee, Y., Park, K.D., Yun, C.-O., 2017. A hydrogel matrix prolongs persistence and promotes specific localization of an oncolytic adenovirus in a tumor by restricting nonspecific shedding and an antiviral immune response. Biomaterials 147, 26–38. Jung, S.-H., Choi, J.-W., Yun, C.-O., Yhee, J.Y., Price, R., Kim, S.H., Kwon, I.C., Ghandehari, H., 2014. Sustained local delivery of oncolytic short hairpin rna adenoviruses for treatment of head and neck cancer. J. Gene Med. 16 (5–6), 143–152. Juweid, M., Neumann, R., Paik, C., Perez-Bacete, M.J., Sato, J., van Osdol, W., Weinstein, J.N., 1992. Micropharmacology of monoclonal antibodies in solid tumors: direct experimental evidence for a binding site barrier. Cancer Res. 52 (19), 5144–5153. Kansal, A.R., Torquato, S., Harsh, G., Chiocca, E., Deisboeck, T., 20 0 0. Simulated brain tumor growth dynamics using a three-dimensional cellular automaton. J. Theoret. Biol. 203 (4), 367–382. Kiefer, J.E., Weiss, G.H., 1984. The pearson random walk. In: AIP Conference Proceedings, 109. AIP, pp. 11–32. Kim, J.-H., Lee, Y.-S., Kim, H., Huang, J.-H., Yoon, A.-R., Yun, C.-O., 2006. Relaxin expression from tumor-targeting adenoviruses and its intratumoral spread, apoptosis induction, and efficacy. J. Natl. Cancer Inst. 98 (20), 1482–1493. Kim, P.-H., Sohn, J.-H., Choi, J.-W., Jung, Y., Kim, S.W., Haam, S., Yun, C.-O., 2011a. Active targeting and safety profile of peg-modified adenovirus conjugated with herceptin. Biomaterials 32 (9), 2314–2326. Kim, Y., Stolarska, M.A., Othmer, H.G., 2011b. The role of the microenvironment in tumor growth and invasion. Progr. Biophys. Mol. Biol. 106 (2), 353–379. Komarova, N.L., Wodarz, D., 2010. ODE models for oncolytic virus dynamics. J. Theoret. Biol. 263 (4), 530–543. Kottke, T., Hall, G., Pulido, J., Diaz, R.M., Thompson, J., Chong, H., Selby, P., Coffey, M., Pandha, H., Chester, J., et al., 2010. Antiangiogenic cancer therapy combined with oncolytic virotherapy leads to regression of established tumors in mice. J. Clin. Investig. 120 (5), 1551–1560. Kwok, K., Groves, M., Burgess, D., 1989. Sterile microencapsulation of BCG in alginate-poly-l-lysine by an air spraying technique. In: Proceedings International Symposium Control Release Bioactive Material, 16, pp. 170–171. Lawler, S.E., Speranza, M.-C., Cho, C.-F., Chiocca, E.A., 2017. Oncolytic viruses in cancer treatment: a review. JAMA Oncol. 3 (6), 841–849. Le Caër, G., 2010. A pearson random walk with steps of uniform orientation and dirichlet distributed lengths. J. Stat. Phys. 140 (4), 728–751. Liu, T.-C., Galanis, E., Kirn, D., 2007. Clinical trial results with oncolytic virotherapy: a century of promise, a decade of progress. Nat. Rev. Clin. Oncol. 4 (2), 101. Lobo, E.P., 2014. Modelling the Role of Interclonal Cooperativity During Early Carcinogenesis. University of Sydney Ph.D. thesis. McKerrow, J.H., Salter, J., 2002. Invasion of skin by schistosoma cercariae. Trends Parasitol. 18 (5), 193–195. Meineke, F.A., Potten, C.S., Loeffler, M., 2001. Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Proliferation 34 (4), 253–266. Metzcar, J., Wang, Y., Heiland, R., Macklin, P., 2019. A review of cell-based computational modeling in cancer biology. JCO Clin. Cancer Inform. 2, 1–13. Mok, W., Stylianopoulos, T., Boucher, Y., Jain, R.K., 2009. Mathematical modeling of herpes simplex virus distribution in solid tumors: implications for cancer gene therapy. Clin. Cancer Res. 15 (7), 2352–2360. Monastra, G., Cabrelle, A., Zambon, A., Rosato, A., Macino, B., Collavo, D., Zanovello, P., 1996. Membrane form of tnfα induces both cell lysis and apoptosis in susceptible target cells. Cellular Immunol. 171 (1), 102–110. Mumenthaler, S., D’Antonio, G., Preziosi, L., Macklin, P., 2013. The need for integrative computational oncology: an illustrated example through MMP-mediated tissue degradation. Front. Oncol. 3, 194.

32

A.L. Jenner, F. Frascoli and A.C.F. Coster et al. / Journal of Theoretical Biology 485 (2020) 110052

Murray, P.J., Edwards, C.M., Tindall, M.J., Maini, P.K., 2009. From a discrete to a continuum model of cell dynamics in one dimension. Phys. Rev. E 80 (3), 031912. Muruve, D.A., Barnes, M.J., Stillman, I.E., Libermann, T.A., 1999. Adenoviral gene therapy leads to rapid induction of multiple chemokines and acute neutrophil-dependent hepatic injury in vivo. Human Gene Therapy 10 (6), 965–976. Neusius, T., Sokolov, I.M., Smith, J.C., 2009. Subdiffusion in time-averaged, confined random walks. Phys. Rev. E 80 (1), 011109. Oh, E., Oh, J.-E., Hong, J., Chung, Y., Lee, Y., Park, K.D., Kim, S., Yun, C.-O., 2017. Optimized biodegradable polymeric reservoir-mediated local and sustained co-delivery of dendritic cells and oncolytic adenovirus co-expressing il-12 and gm-csf for cancer immunotherapy. J. Control. Release 259, 115–127. Okegawa, T., Pong, R.-C., Li, Y., Hsieh, J.-T., et al., 2004. The role of cell adhesion molecule in cancer progression and its application in cancer therapy. Acta Biochimica Polonica-English Edition- 51, 445–458. Osborne, J.M., Fletcher, A.G., Pitt-Francis, J.M., Maini, P.K., Gavaghan, D.J., 2017. Comparing individual-based approaches to modelling the self-organization of multicellular tissues. PLoS Comput. Biol. 13 (2), e1005387. Paiva, L., Martins, M., Ferreira, S., 2011. Questing for an optimal, universal viral agent for oncolytic virotherapy. Phys. Rev. E 84 (4), 041918. Parato, K.A., Senger, D., Forsyth, P.A., Bell, J.C., 2005. Recent progress in the battle between oncolytic viruses and tumours. Nat. Rev. Cancer 5 (12), 965. Park, C.C., Bissell, M.J., Barcellos-Hoff, M.H., 20 0 0. The influence of the microenvironment on the malignant phenotype. Mol. Med. Today 6 (8), 324–329. Paul, C.D., Mistriotis, P., Konstantopoulos, K., 2017. Cancer cell motility: lessons from migration in confined spaces. Nat. Rev. Cancer 17 (2), 131. Phan, D., Wodarz, D., 2015. Modeling multiple infection of cells by viruses: challenges and insights. Math. Biosci. 264, 21–28. Riley, M., Vermerris, W., 2017. Recent advances in nanomaterials for gene deliveryâa review. Nanomaterials 7 (5), 94. Russell, S.J., Peng, K.-W., Bell, J.C., 2012. Oncolytic virotherapy. Nature Biotechnol. 30 (7), 658–670. Ruzek, M.C., Kavanagh, B.F., Scaria, A., Richards, S.M., Garman, R.D., 2002. Adenoviral vectors stimulate murine natural killer cell responses and demonstrate antitumor activities in the absence of transgene expression. Mol. Therapy 5 (2), 115–124. Sherar, M., Noss, M., Foster, F., 1987. Ultrasound backscatter microscopy images the internal structure of living tumour spheroids. Nature 330 (6147), 493–495. Smith, J.P., Kanekal, S., Patawaran, M.B., Chen, J.Y., Jones, R.E., Orenberg, E.K., Ning, Y.Y., 1999. Drug retention and distribution after intratumoral chemother-

apy with fluorouracil/epinephrine injectable gel in human pancreatic cancer xenografts. Cancer Chemother. Pharmacol. 44 (4), 267–274. Stosich, M.S., Mao, J.J., 2005. Stem cell-based soft tissue grafts for plastic and reconstructive surgeries. In: Seminars in plastic surgery, 19. Copyright© 2005 by Thieme Medical Publishers, Inc., 333 Seventh Avenue, New, pp. 251–260. Syverton, J.T., Berry, G.P., 1947. Multiple virus infection of single host cells. J. Exper. Med. 86 (2), 145–152. Tseng, S.-J., Kempson, I.M., Huang, K.-Y., Li, H.-J., Fa, Y.-C., Ho, Y.-C., Liao, Z.-X., Yang, P.-C., 2018. Targeting tumor microenvironment by bioreduction-activated nanoparticles for light-triggered virotherapy. ACS Nano 12 (10), 9894–9902. Van Liedekerke, P., Palm, M., Jagiella, N., Drasdo, D., 2015. Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results. Comput. Particle Mech. 2 (4), 401–444. Wade, S.J., Zuzic, A., Foroughi, J., Talebian, S., Aghmesheh, M., Moulton, S.E., Vine, K.L., 2017. Preparation and in vitro assessment of wet-spun gemcitabine-loaded polymeric fibers: towards localized drug delivery for the treatment of pancreatic cancer. Pancreatology 17 (5), 795–804. Wang, Y., Yuan, F., 2006. Delivery of viral vectors to tumor cells: extracellular transport, systemic distribution, and strategies for improvement. Ann. Biomed. Eng. 34 (1), 114–127. Wee, S., Gombotz, W., Fanslow, W., 1995. Evaluation of alginate microbeads for intranasal delivery of ovalbumin. In: Proceedings International Symposium Control Release Bioact Mater, 22, pp. 566–567. Wein, L.M., Wu, J.T., Kirn, D.H., 2003. Validation and analysis of a mathematical model of a replication-competent oncolytic virus for cancer treatment: implications for virus design and delivery. Cancer Res. 63 (6), 1317–1324. Weiswald, L.-B., Bellet, D., Dangles-Marie, V., 2015. Spherical cancer models in tumor biology. Neoplasia 17 (1), 1–15. Wodarz, D., 2001. Viruses as antitumor weapons. Cancer Res. 61 (8), 3501–3507. Wodarz, D., Hofacre, A., Lau, J.W., Sun, Z., Fan, H., Komarova, N.L., 2012. Complex spatial dynamics of oncolytic viruses in vitro: mathematical and experimental approaches. PLoS Comput. Biol. 8 (6), e1002547. Yamaguchi, H., Wyckoff, J., Condeelis, J., 2005. Cell migration in tumors. Curr. Opin. Cell Biol. 17 (5), 559–564. Yokoda, R., Nagalo, B.M., Vernon, B., Oklu, R., Albadawi, H., DeLeon, T.T., Zhou, Y., Egan, J.B., Duda, D.G., Borad, M.J., 2017. Oncolytic virus delivery: from nano-pharmacodynamics to enhanced oncolytic effect. Oncolytic Virother. 6, 39.