Examining the efficacy of localised gemcitabine therapy for the treatment of pancreatic cancer using a hybrid agent-based model

The prognosis for pancreatic ductal adenocarcinoma (PDAC) patients has not significantly improved in the past 3 decades, highlighting the need for more effective treatment approaches. Poor patient outcomes and lack of response to therapy can be attributed, in part, to a lack of uptake of perfusion of systemically administered chemotherapeutic drugs into the tumour. Wet-spun alginate fibres loaded with the chemotherapeutic agent gemcitabine have been developed as a potential tool for overcoming the barriers in delivery of systemically administrated drugs to the PDAC tumour microenvironment by delivering high concentrations of drug to the tumour directly over an extended period. While exciting, the practicality, safety, and effectiveness of these devices in a clinical setting requires further investigation. Furthermore, an in-depth assessment of the drug-release rate from these devices needs to be undertaken to determine whether an optimal release profile exists. Using a hybrid computational model (agent-based model and partial differential equation system), we developed a simulation of pancreatic tumour growth and response to treatment with gemcitabine loaded alginate fibres. The model was calibrated using in vitro and in vivo data and simulated using a finite volume method discretisation. We then used the model to compare different intratumoural implantation protocols and gemcitabine-release rates. In our model, the primary driver of pancreatic tumour growth was the rate of tumour cell division. We were able to demonstrate that intratumoural placement of gemcitabine loaded fibres was more effective than peritumoural placement. Additionally, we quantified the efficacy of different release profiles from the implanted fibres that have not yet been tested experimentally. Altogether, the model developed here is a tool that can be used to investigate other drug delivery devices to improve the arsenal of treatments available for PDAC and other difficult-to-treat cancers in the future.


Introduction
Inoperable pancreatic ductal adenocarcinoma (PDAC) has a dismal prognosis, with a median survival of 3−5 months for untreated disease [1]. Treatment of PDAC with the chemotherapeutic agent gemcitabine can achieve clinical benefit and symptom improvement in 20−30% of patients [1,2], although PDAC is still regarded as a chemotherapy-resistant tumour [3,4]. Gemcitabine is designed to target and kill cancer cells by incorporating into the DNA strand of a PDAC cell allowing only one deoxynucleotide to be incorporated, which prevents strand elongation [5,6], resulting in cell cycle arrest and subsequent cell death [7,8]. Despite gemcitabine being established as a standard treatment for advanced PDAC over 20 years, most subsequent large phase III studies have not shown significantly improved survival benefit [9]. Overall prognosis for PDAC has seen little improvement in the last 3 decades, largely due to drug resistance and poor intratumoural drug accumulation.
The majority of chemotherapeutics, gemcitabine included, are administered systemically via bolus or infusion intravenous administration. This often results in significant systemic toxicity, with only a fraction of the injected dose reaching the tumour. As such, there has been a growing interest in the development of localised targeted delivery systems which can modify the bio-distribution of drugs and achieve local drug accumulation in the tumour tissue [10][11][12] (Fig 1). For example, drug-eluting polymeric implants are designed to deliver high concentrations of chemotherapeutic drugs directly at the tumour site over some period of time, overcoming transport and tissue barriers as well as limiting off-target toxicities [13]. Biodegradable implants, can be designed to provide sustained drug release over weeks or months, avoiding repeated external drug dosing, clinic visits and other surgical interventions. The characteristics of these devices make local delivery especially attractive for chemotherapeutics with a narrow therapeutic window or short in vivo half-life [14], such as gemcitabine. However, there is still some discussion around the optimal way in which the implanted devices should release drug, for example whether the rate at which drug leaves the device is variable or fixed.
Drug-loaded polymeric fibres can be prepared by various cross-linking methods and allow for drug molecules to be released in a controlled manner depending on the cross-linking type and methods [15]. Previously, Wade et al. [14] showed that wet-spun gemcitabine-loaded alginate fibres inhibited ex vivo PDAC spheroid growth and reduced PDAC cell viability compared to gemcitabine delivered as a free drug. In a subsequent study, Wade et al. [13,16] showed that a coaxial fibre formulation, in which the alginate was encased by a polycaprolactone (PCL) shell, demonstrated significant in vivo antitumour efficacy; however, it is not possible to conclude experimentally whether an alternative release-profile of gemcitabine may be more effective. Fortunately, computational and mathematical modelling is well situated as a predictive tool for quantifying the efficacy of alternative drug release profiles and drug administration patterns.
Mathematical models have been used to help understand formation and treatment of a range of different cancers for some time now [17][18][19][20][21][22]. In particular, agent-based models (ABM) have been used extensively in cancer modelling as they allow for the consideration of spatial and phenotypic heterogeneity [23][24][25][26][27][28][29][30][31] which are known to be major drivers of variations in treatment outcomes. In ABMs, the likelihood of events, such as cell proliferation, movement, death or mutation are modelled as probabilities, allowing the simulation to evolve stochastically in time. ABMs have been used to contribute understanding to tumour growth dynamics [24,[26][27][28], such as angiogenesis [31] and cell cycling [32]; the impacts of certain treatments [30,33,34]; and the likelihood of tumour recurrence [23]. ABMs have also been used to model chemotherapy treatment and movement through a tumour [35][36][37][38][39][40][41]. For example, Tang et al. [38] developed a 3D computational model of tumour growth under treatment with chemotherapy, where cells are modelled as discrete agents whose proliferation is Sustained-delivery implants are a promising treatment methodology over conventional single free-drug intravenous or intrantumoural injections. A hypothetical comparison of drug concentrations at the tumour site under these two protocols is pictured. Systemic injections of anti-cancer drugs often result in only a fraction of the drug arriving at the tumour site followed by a rapid decrease of drug concentration at the tumour site. In comparison, sustained-release mechanisms deliver drug over a prolonged period resulting in a durable drug presence at the tumour site. Created using biorender.com. dependent on interstitial pressure and chemotherapy is captured by a partial differential equation (PDE). With their model they investigated the effects of the tumour microenvironment on the intratumoural distribution of chemotherapy and concluded that the main driver was the pressure difference within capillary blood and extracellular space. While important for intravenously administered chemotherapy, their model is too complex for considering an implanted degradable chemotherapy device. Furthermore, there is yet to be developed an ABM for the release of chemotherapy treatment through sustained-release or degradable devices. Insights on therapeutic failure in immunotherapy and virotherapy have also been provided through an ABM software known as PhysiCell [42][43][44][45]. There have been ABMs developed that specifically focus on pancreatic cancer growth [46,47]; however, an ABM describing pancreatic cancer growth and treatment with a degradable polymer implant has not yet been developed.
For some time, mathematical models of degradable drug delivery mechanisms have been used to assist in the understanding of polymer degradation, hydrolysis kinetics and the subsequent effect of drug release on the applied system [10,[48][49][50][51][52][53][54][55]. Using mass-balance kinetic equations, McGinty et al. [52] investigated the extent to which variable porosity drug-eluting coatings can provide better control over drug release using transport diffusion equations.
Their results indicate that the contrast in properties of two layers can be used as a means of better controlling the release, and that the quantity of drug delivered in early stages can be modulated by varying the initial drug distribution. More recently, Spiridonova et al. [56] fitted drug release from polymer microparticles and investigated the effect of size distribution on diffusional drug release from sustained-delivery systems using a system of PDEs. Whilst useful for capturing the drug delivery mechanism, most models of drug-loaded polymers such as these have not examined the influence of changes to drug release profiles on antitumour efficacy or how intratumoural stochasticity impacts drug delivery.
In this work, we have developed a hybrid mathematical and computational model of PDAC tumour growth and death from treatment with gemcitabine released from a polymeric fibre. We extended a previously published ABM known as a Voronoi cell-based model (VCBM) [57] to model tumour cell growth and death and coupled this with a PDE model for gemcitabine release from polymeric implants. In vitro drug release curves were used to optimise the PDE formulation describing how gemcitabine is released from fibres. A numerical simulation was then used to initialise the parameters in the ABM using in vivo control PDAC tumour growth measurements. The potential impact of these fibres on tumour growth and cell death was then investigated with the VCBM-PDE model and improvements on drug release kinetics and fibre placement were suggested. We quantified the impact of varying sustained-release profiles including a constant release, exponential release, and sigmoidal release from these devices. In this way, we could conclude whether the rate at which drug was released from the implant had any impact on treatment effectiveness. The model was developed as a tool that can be applied to interrogate other cancer therapies using polymeric implants with the goal to improve treatment response for PDAC patients.

Experimental methods
Ethics statement. All animal experiments were conducted in accordance with the NHMRC Australian Code for the Care and Use of Animals for Scientific Purposes, which requires 3R compliance (replacement, reduction and refinement) at all stages of animal care and use, and the approval of the Animal Ethics Committee of the University of Wollongong (Australia) under protocol AE18/13.
Fibre fabrication and characterisation. Full details for the fabrication and characterisation of alginate fibres loaded with or without gemcitabine are described in Wade et al. [13,14].
Briefly, gemcitabine-loaded alginate fibres had a uniform surface area from 50−120 μm in diameter. Experimentally, fibres displayed different drug release profiles and total drug released depending on the concentration of polymer used. Primarily, the effect of changing the alginate percentage and PCL presence was seen in the speed of release of the drug and time to full drug release; see [14] for more details. To investigate this further in this work, using a computational model we consider alternative release profiles not yet fabricated to examine their effectiveness, these include constant, exponential and sigmoidal release gradients. These are considered purely in silico. Fibre diameter also varied depending on the materials used [14].
Fibre gemcitabine release kinetics. Full details for the experiments measuring gemcitabine release can be found in Wade et al. [14] with brief details here. Gemcitabine-loaded fibres were added to 2mL of simulated body fluid (SBF), Ph 7.4 and incubated at 37˚C. At various time points (10, 30, 60, 90 min hourly for 10h and then daily for 3 weeks), buffer solution (200μL) was removed for analysis of gemcitabine release and replaced with fresh SBF. The amount of drug released from alginate fibres was determined using high performance liquid chromatography (HPLC). The amount of gemcitabine released (μg) was calculated by interpolating AUC values from the standard curve using Empower Pro V2 (Waters) software.
Implant toxicity in vitro. Gemcitabine loaded fibres were tested for their cytotoxicity against human pancreatic cancer cells (Mia-PaCa-2) cells over 72h. Cells were incubated with 0.5 cm lengths of gemcitabine loaded or non-drug loaded fibre formulation before an endpoint MTS cell viability assay was performed. Results are displayed as a percentage of an untreated control. Experiments were performed in triplicate. Full details for the toxicity experiments can be found in Wade et al. [13].
In vivo Mia-PaCa-2 cell growth. Animals were subcutaneously inoculated with 100μL suspension of 1×10 6 Mia-PaCa-2 cells in PBS. Tumour volume measurements began when tumours reached a volume of 200 mm 3 using where width is the longest tumour diameter measurement and length is the tumour measurement along the axis perpendicular to this (Fig H in S1 Text). Tumour volume was measured daily for a duration of approximately 33 days. Full details for this experiment can be found in Wade et al. [13]. All animal experiments were conducted in accordance with the National Health and Medical Research Council (NHMRC) Australian Code for the Care and Use of Animals for Scientific Purposes, which requires 3R compliance (replacement, reduction, and refinement) at all stages of animal care and use, and the approval of the Animal Ethics Committee of the University of Wollongong (Australia) under protocol AE18/13.

Mathematical methods
The model developed for the release of gemcitabine from alginate fibres and the impact on a growing PDAC tumour is formulated in two parts below. The first describes the PDE describing the concentration of gemcitabine in the tumour microenvironment (TME) and surrounding tissue over time. The second describes the VCBM [57] that captures the way tumour cells proliferate, move and undergo apoptosis from gemcitabine. We chose to model the tumour growth and treatment in 2 dimensions as the in vivo tumour growth measurements were taken as only width and height measurements and we assumed any third-dimension effects were analogous to what happens in 2 dimensions. All parameters introduced for the model are summarised in Tables A-E in S1 Text and a schematic for the model is in Fig 2. Model of gemcitabine. To capture the concentration of gemcitabine in the tumour microenvironment, we first considered a 2D rectangular domain with boundary B (Fig A in  S1 Technical Supplementary Information). Inside this domain, is implanted a gemcitabine drug-loaded fibre which is represented by a vertical line source (panel A of Fig A in S1 Technical Supplementary Information and Fig 2A). Gemcitabine diffuses from the line source at some time-dependent rate that decreases as the polymeric fibre degradation slows. The The concentration of drug in the TME was modelled in a 2 dimensional domain bounded by B, where C(x,y,t) was the concentration in the TME at position (x,y). The fibre implant was then placed at a position x = x F and modelled as a line source. To capture the diffusion of drug from the fibre, we modelled the concentration of gemcitabine inside the fibre F(r,y,t) at radial position r and domain position y where the continuity condition in Eq (3) required equal concentrations at the fibre boundary and at the immedicate local microenvironment, i.e. F(r total ,y,t) = C(x F ,y,t). (B) The concentration of gemcitabine inside the polymeric fibres was modelled by radially symmetric diffusion Eq (2) using a finite volume method (FVM) discretisation and considering the 2D cylindrical cross section of the fibres which have length L and radius r total . The fibre was discretised into concentric annuli F m,j at annulus m and cross section j, (i = 0,1,. . .,M) and the concentration of drug in each annulus F m,j was modelled by considering drug diffusion across the bounadaries (e.g. F m−1,j and F m+1,j flow into F m,j and vice vera). The full discretisation is presented in S1 Technical Supplementary Information. (C) Modelling assumptions for the VCBM were that cancer cells (orange) proliferate and some are able to cause epithelial to mesenchymal transtion and become invasive. We model this transition by assuming cells differentiate into a mesenchymal cancer cell (MCC) with one daughter cell placed on a neighbouring healthy cell. These MCCs cause the break down of surrounding tissue (i.e. replace healthy neighbouring cells with their progeny). Cancer cells can then die through gemcitabine uptake from their local environment. For more details see Fig B in S1 Text. (D) Individual cells were modelled as cell centres connected by springs [57]. The proliferation of a cell introduced a new cell into the lattice network which caused the rearrangement of the cells in the lattice with movement governed by Hooke's law. (E) To simulate the gemcitabine concentration in the TME, Eq (1), we introduced a FVM discretisation, where the gemcitabine concentration was defined at discrete volumes centered around points in the discretisation. Cells could take up drug from the nearest grid point to their centre, and this concentration was used to determine their likelihood of drug-induced cell-death.
https://doi.org/10.1371/journal.pcbi.1010104.g002 gemcitabine concentration in the domain outside the fibre is diffusing and decaying. PDAC cells in the domain are also taking up gemcitabine, removing it from the concentration in the domain. Inside the fibre, we model the diffusion of drug as radially symmetric (Fig 2B).
We denote the concentration of drug in the TME at position (x,y) by C(x,y,t) and model this concentration by zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl { cell uptake þ dðx À x F ÞJðy; tÞ; zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl fflffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl fflffl { fibre release where D is the diffusion coefficient in the TME, and λ is the decay rate of the drug. To model cancer cells taking up gemcitabine, we used δ(x) which is the Dirac delta function in onedimension, where (x i , y i ) is the ith cancer cell's Voronoi centre position in the domain (Fig A  in S1 Text), and W i is the cell's volume. Pancreatic cancer cells take up drug in the domain at a rate ν c . Cell uptake was modelled by point sinks analogous to that in PhysiCell and BioFVM [43,58], where cells are considered discrete "point masses" in the domain that take up drug from a single rectangular discretised voxel weighted by the local concentration of drug. We then used a line source at x = x F , y 0 �y�y 0 +L to model the release of gemcitabine from the polymeric fibre, where y 0 is the location of the bottom of the fibre and L is the fibre length ( To derive the flux of drug from the line source, we first assumed that the release of drug from the fibre would be time dependent. As such, we chose to explicitly model a concentration of drug diffusing inside the fibre. We denote the concentration of gemcitabine at radial position r and location (x F , y) by F(r,y,t) (Fig B in S1 Technical Supplementary Information and Fig 2A). We model the diffusion and movement of drug inside the fibre assuming radial symmetry. We assumed that diffusion in the radial direction is significantly faster than along the fibre since the radius of the fibre r total is significantly less than the length of the fibre L (Figs A and B in S1 Technical Supplementary Information). This gives where D F (t) is the time-dependent diffusion of drug inside the fibre. We imposed the continuity condition so that the diffusion of drug out of the fibre at the line source will depend on the location (x F , y) and local exterior concentration. The flux out of the line source J(y,t) in Eq (1) can then be approximated from the release of drug across the boundary of the fibre: : This term is derived by converting the flux out of the radial fibre into the flux represented by the line source in Eq (1) and converting to a concentration per surface area where h is the depth of the rectangular region (presumed thin, see Fig A in S1 Technical Supplementary Information). Both Eq (3) and Eq (4) are necessary boundary conditions for Eq (1) and Eq (2). In this way, we assume the concentration is continuous and the flux of the fibre is equal to the flux into the TME, equivalent to a conservation of mass.
The diffusivity of the drug, D F (t), is modeled by the function where k controls the gradient of the drug release rate(i.e. how quickly the fibre swells), D const is the constant decay rate from the fibre and � is a tuning constant to provide a finite initial diffusion coefficient, i.e. D F (0) = k/�+D const . We expect D F (0) to be initially large (>1) since the polymeric fibre is hydrophilic and drug would immediately diffuse out of the fibre. In addition, some drug is never properly loaded into the fibre and can be released instantaneously. The formalism in Eq (5) was broadly chosen to capture the rapid decline in release as the polymeric fibre degrades. It is possible to model the breakdown of the drug release mechanisms to include device swelling and degradation and for examples of this see [56,[59][60][61].
No-flux boundary conditions on B, the exterior of the TME, are imposed: where C 0 is the amount of drug in μg, the denominator is the volume of the fibre and there is no drug initially in the domain B. We assume the location of the fibre is fixed in space over the course of the simulation and is not affected by cells around it. For more details on the derivation of the model see S1 Technical Supplementary Information. We solved Eqs (1)-(4) numerically using a Finite Volume approximation. In particular, the diffusion of drug within the fibre, Eq (2), was solved through discretising the cross section of a fibre into annuli (see Fig 2B and S1 Technical Supplementary Information). The model is solved using a finite volume method (FVM) discretisation, for examples of this form of discretisation in cancer growth and treatment see [62][63][64][65][66][67][68][69][70]. As part of the investigations of these fibres in this work, we chose to quantify the impact of varying the drug release profile to a constant release, exponential release, or sigmoidal release profiles. To do this, we replaced the drug diffusion in Eq (2) by the corresponding profile being examined and imposed drug conservation on the total drug released. More details on this can be found in S1 Technical Supplementary Information. Voronoi Cell-Based Model (VCBM) of pancreatic tumour growth. Agent-based models (ABMs) are primarily used to simulate heterogeneity that arises through stochasticity in cellular interactions. We present an ABM to capture the 2D formation of a pancreatic tumour in the pancreas. Our model extends a Voronoi cell-based model (VCBM) for tumour growth already published in [57]. The model describes how individual cells behave over time by considering their behaviour to be a stochastic process. We consider a set of (x,y) points in the domain B as representatives of cell centres and then overlay this with a Voronoi tessellation to define individual cell boundaries. A Voronoi tessellation defines the region of space where the Euclidean distance to a point is less than the distance to any other point e in the lattice. Voronoi tessellations have been used to model tissue and cancer cell dynamics for some time [71][72][73][74][75]. Using a Voronoi tessellation for the ABM allows cell morphology to be heterogeneous and not fixed, and the morphology can change with cell movement. The model is solved on a time increment of 1hr to account for the fact that cellular interactions are slow in comparison to drug diffusion (Fig C in S1 Technical Supplementary Information). To model pancreatic tumour formation, we assumed the primary functions of pancreatic tumour cells were movement and proliferation. Below are details of the cell types, the model for cell movement and proliferation, a description of the dynamics of tumour mesenchymal cells, the model for cell death and details of how the domain changes as the tumour grows.
PDAC cells can acquire mesenchymal-like phenotype properties through a process known as epithelial-mesenchymal transition (EMT) [76][77][78][79]. In the EMT process, epithelial elements undergo cytoskeleton remodelling and migratory capacity acquisition due to the loss of intracellular contacts and polarity [77]. This enables the formation of mesenchymal-like cancer cells (MCCs) which have enhanced migratory capacities and invasiveness, as well as elevated resistance to apoptosis [78]. Since there is evidence that EMT plays an important role in PDAC progression [76][77][78][79], we have introduced this cell type into the model.
We Each cell has an initial location in domain B such that for cell i its cell centre would be at The initial tissue comprises of healthy cells, arranged so that the corresponding Voronoi cells form a hexagonal tessellation, analogous to other work in the literature [80,81]. To initialise the tumour formation, we removed a healthy cell from the centre of the domain and replaced it with a pancreatic tumour cell (Fig A in S1 Text). These pancreatic tumour cells could proliferate, die from gemcitabine, or form MCCs. Once formed, these MCCs then move and proliferate until they die. Healthy cells are assumed to be able to move or be removed from the domain by being replaced by MCCs (Fig B in S1 Text).
Cell movement is governed by pressure-driven motility, modelled using Hooke's law [57]. Each cell's (x,y) position is updated by calculating the effective displacement of the cell's lattice point by the sum of the forces exerted on that cell, where force is modelled as a network of damped springs connecting a cell to its nearest neighbours (defined by a Delaunay triangulation). To determine the movement of an individual cell, we use the vector representation for that cell's (x,y) position from the origin r For cell i, the displacement of this point in time Δt cells is given by where r ! i ðtÞ is the vector representation for the ith cell's centre position in the lattice at time t, λ m is a damping and mobility constant, r ! i;j is the vector between cell i and j, s i,j is the spring rest length (equilibrium distance) between cell i and j. Cells j are the cells connected to cell i in a Delaunay triangulation, i.e. in a neighbourhood of cell i. The introduction of new cells in the lattice through proliferation introduces new spring connections and shortens or extends others, promoting the movement of cells in the environment (Fig 2D).
Tumour cell proliferation was assumed to be a function of the cell's distance, d neut , to the nutrient source (tumour periphery, i.e. nearest healthy cell centre, see Fig C in S1 Text). The maximum radial distance for nutrient-dependent cell proliferation is d max . Cells that are a further distance from the nutrients than d max enter a quiescent (non-proliferative state), forming what is commonly known as a necrotic core. The probability of a cell dividing p d in time step Δt cells is given by where p 0 is a proliferation constant derived based on the maximum rate of cell proliferation ϕ (i.e. p 0 ¼ 1 À expðÀ �Dt cells Þ � �Dt cells ). The formalism in Eq (8) is similar to what was used by Kansal et al. [82], and Jiao and Tarquato [83]. Mechanical feedback between cells has been shown to be a regulatory mechanism for growth control of tissue cells in vitro and in vivo [84]. A cancer cell's ability to proliferate, while generally considered to be dysregulated, has been shown to slow down and eventually arrest with compression from neighbouring cells [85,86]. To model this, we introduced a constraint on cell proliferation: if all cells within a cell's neighbourhood (i.e. connected to that cell by a Delaunay triangulation) are within s/p age of the cell, then the cell will not proliferate. If a cell i proliferates, a new lattice point j is created and the two cells are placed at a distance s/p age from the original proliferating cells position at a rotation θ~2 U(0,2π] (Fig 2D). To simulate the enlargement and repositioning of the daughter cells, the resting spring length of the connection between i and j linearly increases over time from s/p age to the mature resting spring length s as was formulated in our previous work [57]; see Fig C in S1 Text. Once a cell has proliferated, it takes g age time steps before the daughter cell will try to proliferate again, accounting for G1 phase of the cell cycle where the cell transitions from mitosis M to DNA synthesis S [57]. It is well known that tumours contain highly heterogeneous populations of cells that have distinct reproductive abilities. To account for heterogeneity in the cell cycling, cells sampled the age at birth from a Poisson distribution with mean 50.
MCCs are created at the boundary of the tumour with probability p MCC . These cells are created from tumour cell differentiation. In this way, a single tumour cell can either proliferate into two new tumour cells or undergo differentiation into a tumour cell and an MCC. We model the invasive property of MCCs by placing the newly differentiated MCC daughter cell at the position of a neighbouring healthy cell, removing that healthy cell from the domain; see Fig B in S1 Text. As such, through their creation, these MCCs contribute to the degradation of the healthy tissue surrounding the tumour.
As in [87][88][89][90], we assumed that cancer cells die from gemcitabine contact at a rate described by the Michaelis-Menten term where δ m is the maximum death rate due to the drug, C i,j is the concentration of drug at the grid position (i,j) in the FVM discretisation closest to the cell's centre (Fig 2E and S1 Technical Supplementary Information), and IC 50 is the concentration at which half the effect of the drug is attained. From this, the probability of an individual cell dying can be determined by assuming Prob(cell death) = 1 À expðÀ bDt cells Þ � bDt cells . While we chose not to model explicitly the resistance to gemcitabine that cancer cells can develop [3,4], we believe that by modelling cell death probabilistically we can capture some of the heterogeneity that may exist intratumourally. If a cell dies, then its phenotype changes to be a dead cell and it takes d age hours to disintegrate. To simulate disintegration, at each time increment the spring rest lengths of a dead cell to each of its neighbours, s i,j , decreases by s i,j /d age . Once the dead cell has disintegrated, it is removed from the simulation and the space left behind becomes "empty space" which is plotted in white. Cells are free to move into this space over time following Hooke's law.
As the tumour grows, the model domain expands. To reduce computational cost, new healthy cells are added to the domain only when a tumour cell's radial distance from empty space is <10 model units (Fig D in S1 Text). A summary of the overall model rules is provided in a decision tree diagram in Fig E in S1 Text.

Numerical simulations and parameter estimation
The VCBM-PDE model was written in C++ and simulations called through Matlab 2021b by creating a definition file for the C++ library using clibgen and build in Matlab 2021b. Code for the model at the various stages (e.g. fibre, single injections) can be found on github (https:// github.com/AdrianneJennerQUT/hybrid-VCBM-of-gemcitabine-and-pancreatic-cancer). Full details on all aspects of the code can be found in S1 Code Documentation.
An approximation for tumour volume was then determined from the 2D simulations using the same formula as the calibre measurements, multiplied by a scalar σ: where width is the longest distance of a cell on the periphery from the centre and length is the distance of the farthest cell from the centre on the radial axis perpendicular to the radial axis of the longest distance (Fig H in S1 Text) where σ unit length of the model is equivalent to 1 mm. This calculation choice was made to closely resemble the tumour volume calculation with calipers done in vivo. As the size of the computational domain was smaller than the size of the real tumour, the length unit was scaled by σ, which scaled the unit length in the VCBM domain to a comparable mm unit measurement that reduced the computational cost. We chose σ = 0.1728. We performed a convergence test for the VCBM to determine a threshold for the minimum number of simulations of the model to have average convergence at day 33 (Fig F in S1 Text). We found that from 500 simulations onward, the average tumour volume at day 33 had converged to some mean value. As such, we picked n = 500 to be the number of simulations of the model we would run at each investigation.
To fit the parameters in the model, we first consider only the drug release and diffusion compartment of the model and using in vitro fibre release data obtain the relevant parameters for this compartment. Following this, we fit the cell proliferation and death rates based on in vitro cell count and viability measurements. Lastly, we estimate the remaining VCBM cell proliferation parameters by using a Latin Hypercube Sampling (LHS).
All fitting was undertaken using lsqnonlin in Matlab 2021b using pdepe and ode45 to simulate the model. Parameters in the model were fit using experimental data or estimated from the literature. To fit the parameters relating to drug release from the fibre we used the in vitro drug release experiments. We simplified the model to consider only one cross section, i.e. F m,j = F m , since the outside concentration of drug was independent of location in the absence of cells in the in vitro experiment.
To estimate parameters for the pancreatic cell growth kinetics, we did a large Latin Hypercube sample of the parameter space and determined parameters that resulted in a minimal least squares distance to the in vivo control tumour growth measurements. To do this, we set an initial seed and obtained 1000 samples of the parameter space and simulated the model 30 times for each sample each time calculating the residual between the data and the simulation. The parameter sample with the lowest average was then fixed for the remainder of the model simulations. It should be noted that this method does not provide us with a fitted parameter value but was a numerical way of obtaining a reliable estimate for these parameters. Other parameters were either fixed to previous values in the literature or estimated based on previous work. See Tables A-E in S1 Text for a full summary of all parameter values and relevant references.

Calibration of drug release kinetics and drug-induced cell death to in vitro measurements
Gemcitabine-loaded fibres were placed in a solution bath and the resulting cumulative concentration of gemcitabine measured (Fig 3A). To obtain a model describing the release rate of the drug from the fibre, we fitted parameters from Eqs (1)-(4) to these in vitro measurements for the release of gemcitabine from 3% alginate 15% PCL fibres [14]. Fitting the release curve parameters k, d const , C 0 and A out gave the fit in Fig 3B and Table A in S1 Text. Overall, the model was able to obtain the fit to the data and followed the trend which showed a rapid initial release of gemcitabine followed by a steady-state threshold. We validated the model's predictive capability by also fitting gemcitabine release from 1% and 2% alginate fibres (Fig G in S1 Text).

parameter values in
To assess the efficacy of the drug on inducing death in PDAC cells, cell viability studies were performed using Mia-PaCa-2 cell lines. To model these experiments, we considered a simplified deterministic and spatially independent version of our model with only live cancer cells P L (t), dead cancer cells P D (t) and a concentration of drug C(t):  Table B in S1 Text. https://doi.org/10.1371/journal.pcbi.1010104.g003 where ϕ is the exponential proliferation rate of cancer cells in vitro, δ m is the death rate of cancer cells by gemcitabine, IC 50 is the drug's half effect concentration, and λ is the decay rate of the drug (Fig 3C). To first determine the proliferation rate of pancreatic cancer cells in vitro, an exponential growth curve was fit to cell count measurements for Mia-PaCa-2 cells [91] (Fig  3D, parameter values Table B in S1 Text) using simple exponential growth (i.e. setting C(0) = 0 in Eq (9)). Fixing this growth rate and the estimate for the decay rate of drug, we then determined the antitumour efficacy of gemcitabine-loaded fibres in the cell viability experiments. Cells were treated with aliquots of simulated body fluid from gemcitabine-loaded fibres that had been incubating for 24, 48 or 72 h (Fig 3E). To simulate these experiments, the model is solved piecewise such that μ(t) = δ(t−t aliquot )C(t aliquot ), where t aliquot are the times of the drug administrations. An approximation for the concentration of drug at each time point, C(t aliquot ), can be determined using the calibrated PDE model for drug release from the fibres. Fitting the drug-induced death rate and IC 50 gave a good approximation to the data (Fig 3F). The resulting parameter values from the fit of the model can be found in Table B in S1 Text.

Calibration and sensitivity of pancreatic tumour growth
The VCBM simulation of pancreatic tumour growth in the absence of treatment depicts invasive and disorganised movement of cancer cells into surrounding healthy tissue (Fig 4A). To calibrate tumour growth parameters in the model, we used an exhaustive numerical search of the parameter space using a Latin Hypercube Sampling for g age , d max , p 0 and p MCC , where we were minimising the least squares of the simulation with the in vivo tumour volume of Mia-PaCa-2 cells over 33 days (Fig 4B and Table C in S1 Text). To obtain an understanding of the stochasticity in our model, we fixed the parameter values obtained and we simulated the model 500 times and plotted the tumour volume over 33 days. From Fig 4B, while the growth is varied at points, there are no distinct outliers or unusual tumour growth rates, and the standard deviation throughout the entire period of observation remains small. In addition, the simulations sit within the in vivo tumour growth measurements from Wade et al. [13] for untreated pancreatic cancer growth. The histogram for the number of MCCs across the simulations (Fig 4C) shows only a small number of MCCs are created over the 33 days of growth, which is realistic when considering the ratio between a single cell agent in the model and a real cell in a biological tumour and matches findings that MCCs will compose only a small subset of the tumour [92][93][94].
To analyse the drivers of pancreatic tumour growth dynamics in our model, we conducted a detailed sensitivity analysis. A systematic multi-parameter sensitivity analysis was performed for p = [p 0 , p MCC , d max , g age , p age ] using weighting identified by Wells et al. [95] (Fig 4D). This sensitivity analysis can identify combinatorial influences of multiple parameters and elucidate systemic features of the model. The average tumour volume predicted by the model at day 33 for 500 simulations was recorded for each parameter set. Pairs of parameters were varied, with each cell of Fig 4D depicting  The time taken for a cell to prepare for mitosis, g age , has the greatest impact on final tumour volume (Fig 4D). Increasing g age decreases tumour volume and conversely a decrease in g age increases the final volume. As a result, the model predicts that if cells take longer to move through the cell cycle and undergo mitosis this will result in a smaller tumour volume. Reducing the maximum distance, d max , from the periphery at which a cell can still proliferate decreases the final tumour volume. This is to be expected, as reducing the proliferating cell rim (through decreasing the distance from the periphery for which cells can proliferate) will reduce the number of cells available to proliferate and subsequently reduce the tumour volume. Decreasing the value of d max only appears to have a significant impact on the final  Fig I in S1 Text). (B) Mia-PaCa-2 tumour volume over 33 days measured in vivo in mice (black, n = 4). Overlaid is the tumour volume from the VCBM simulation (pink, n = 500) with parameters from Table C in S1 Text. (C) MCC counts in the VCBM simulations (n = 500). (D) Sensitivity analysis of control tumour growth. Maximum tumour volume over 33 days for perturbations of parameters with weights of 0.25, 0.75, 1.25, 1.75 and 2.25, and spatial plots of large and small tumours simulated using the depicted weightings. In the heatmap, each pixel represents 500 averaged simulations with two parameters. In the boxes, the parameters vertically and horizontally in the grid are the weightings in ascending order, with each pixel being a "coordinate" representing the weighting for each parameter and the result from 500 averaged tests. Diagonal pixels only use individual parameters with different weightings. Legend for cell colouring: cancer cell (orange) healthy cell (grey), MCC (blue). https://doi.org/10.1371/journal.pcbi.1010104.g004

PLOS COMPUTATIONAL BIOLOGY
tumour volume when the weighting applied is �50%. In comparison with d max and g age , the tumour volume is insensitive to changes in both the probability of a cell proliferating if it has reached mitosis, p 0 , and the probability of a new pancreatic cancer stem cell being created, p MCC . The time taken for a cell to reach adult size (when it can proliferate), g age , similarly has a negligible impact on the tumour volume.

Intratumoural implantation provide an alternate effective protocol
As previously shown by Wade et al [13][14][15][16], pancreatic cancer growth was inhibited under administration of gemcitabine as a free drug. This growth inhibition was furthered when gemcitabine was administered in loaded alginate fibres. Before quantifying the efficacy of varying gemcitabine-loaded fibre characteristics, we first looked to evaluate the impact of varying protocols for single point free-drug injections (Fig 1) of gemcitabine on the tumour volume. Simulating single point free-drug injections with the VCBM-PDE is a simplification of the full model presented in Eqs (1)-(4) where F(r,y,t) = 0. More details on this can be found in S1 Technical Supplementary Information. We considered free-drug injections of gemcitabine as administered along a radial axis of the tumour in either a single dose or four free-drug injections which are rotationally symmetric (Fig 5A). In the case of the four injections, the total dosage is spread across the injections so that the total amount of drug administered is conserved.
Simulations of the model under the different injection protocols can be found in Figs 5B and 5C, and J in S1 Text. The sensitivity of parameter values governing tumour volume were again probed, now under a single administration of gemcitabine at the centre of the tumour (Figs 5D and 5E, and K in S1 Text). The same trends with g age and d max were observed; however, an additional parameter, which represents the concentration the drug required to have an impact on the tumour volume, IC 50 , was found to influence the volume under further perturbations of the parameter value (Fig 5E). As expected, a lower value of IC 50 , which indicates that a smaller concentration of the drug is required for it to influence cancerous cells, leads to a lower tumour volume, while an increased value of IC 50 leads to a higher tumour volume when compared to original estimate for IC 50 .
To determine the effect of injection placement on tumour volume over time, we considered two protocols. The first was a single injection into the tumour of concentration C 0 , and the second was four injections into the tumour, each of concentration C 0 /4 placed at equal rotations around the tumour. For each of these protocols we considered five potential placements of each injection at a distance d m from the centre: a central injection (d m = 0), and injections d m = 0.9 mm from the centre, d m = 1.7 mm from the centre, d m = 2.5 mm from the centre and d m = 3.5 mm from the centre (Fig 5A). For example, for the second protocol where four injections are being placed, we could choose a distance d m = 2.5mm from the centre of the tumour for the placement of the injections and the four injections were then placed at equal rotations around the tumour centre. For each of these placements, 500 simulations were run over 33 days and both the number of tumour cells and the tumour volume over time were measured (Fig 5F). For a single injection, distance did not impact the effectiveness of the injection and the tumour volume is qualitatively similar. The tumour volume was more significantly affected by distance in the case of four injections (Fig 5F), with free-drug injections further away from the centre of the tumour performing worse than those intratumoural injections. Primarily, single free-drug injections implanted peritumourally may encourage branching of external tumour structures in the model, and hence increase the calculated volume as it is based on the maximum distance from the centre of the tumour to the edge. Tumours with invasive edges or branching formations have been shown in some cases to be more challenging to treat in the long term [83,96]. From our model, this suggests that avoiding injection protocols that induce branching is crucial. While we present an approximation for tumour volume and placement of injections in units relevant to in vivo models (i.e. mm 3 and mm respectively), more work needs to be done to validate that the efficacy of treatment predicted by the model would map to the human scale.

Fibre location and release kinetics are a major driver of tumour arrest or tumour growth
Using the VCBM-PDE, we analysed the impact of varying the position of the fibre and the initial drug concentration on the tumour growth dynamics (Fig 6A). We introduced three classifications for the tumour growth dynamics: tumour eradication (i.e. a tumour volume <1mm 3 ) tumour stabilisation, i.e. a tumour volume at day 33 less than the initial tumour size (� 100mm 3 ), and tumour growth, i.e. a tumour volume on day 33 greater than the initial tumour volume. Large concentrations of gemcitabine loaded into the fibre positioned at d m = 3.5 mm or d m = 4.3 mm from the tumour centre were unable to stabilise or eradicate the tumour, also known as tumour arrest (Figs 6B and 6C, and M in S1 Text). Once the fibre was positioned closer to the tumour centre (�1.7 mm) lower concentrations of drug were sufficient to result in stabilisation of the tumour growth (Fig 6B). It was only with high drug concentration and centered fibres that we saw complete tumour eradication (Fig L in S1 Text).
There are large variations in the response of tumour growth to the different protocols, suggesting that tumour stabilisation or arrest might be achievable for some tumours whereas others might experience tumour growth even in the presence of drug-loaded fibre.
To then analyse the effects of changes to the drug release profile on the tumour growth, we investigated four different release profiles: constant release, exponential release, sigmoidal Emax/Imax release profiles [97][98][99] (See Section TS3 in S1 Technical Supplementary Information,). Each of these release profiles were parameterised by a release rate γ and for the Emax and Imax curves a half-effect term η. The different release profiles were tested with the fibre placed either centrally (intratumourally) (Fig 6D) or on the periphery of the tumour (peritumourally) (Fig 6E). The four different release profiles (constant, exponential, sigmoid emax, sigmoid imax) were tested with 8 different release rates. For each parameter value, 500 simulations were run over 33 days, with an initial amount of 500 μg of gemcitabine.
For fibres positioned in the centre (Fig 6D), it is possible to eradicate the tumour with all release profiles considered given a small enough value of γ. In comparison, none of the drug release profiles resulted in tumour eradication when positioned peripherally (Fig 6E). Although, of note, the exponential release profile performed best out of the four release cases in both the centered and peripheral fibre configurations. Overall, the results suggest that it is possible to reduce tumour size with peripheral fibre injections irrespective of the shape of the release curve.

Discussion
PDAC is a difficult-to-treat cancer with a poor prognosis. Novel therapeutic interventions are desperately needed to improve patient survival. While chemotherapy drugs, such as gemcitabine, have shown durable efficacy for pancreatic cancer, there has been little to no improvement in patient survival in the last 30 years [100]. PDACs are notorious for a dense fibrotic stroma that is interlaced with ECM [101] and is a major cause of therapeutic resistance [102]. One way of improving drug retention at the tumour site, and by consequence increase tumour eradication and patient survival, is through sustained-delivery devices (Fig 1). Polymeric fibres loaded with gemcitabine have shown increased therapeutic efficacy over conventional treatment delivery. To further analyse the potential of these novel therapeutic implants, we have designed a hybrid Voronoi cell-based model (VCBM)-partial differential equation (PDE) model to describe pancreatic tumour formation in healthy pancreatic tissue and the resulting effect of gemcitabine on the tumour tissue when delivered locally. With this model, we considered both the impact of a single fibre implanted with varying drug release profiles and hypothesised alternative and more effective treatment protocols.
The model was calibrated to in vitro and in vivo data. A limitation of this calibration was the lack of spatial data available to calibrate the parameters in the model. To try and quantify the impact of this limitation we performed a range of sensitivity analyses at different stages of the model investigation. The parameter sensitivity analysis then revealed that the fundamental driver of tumour growth in our model was the rate of cell mitosis. The idea that the cell cycling time is a fundamental part of tumour progression has been found in other mathematical models [103], suggesting that the model's sensitivity in terms of tumour volume is in line with other models in the literature. It is also known that molecules can modulate the cell cycle of cancer cells, changing the cancer aggressivity. For example, melatonin is a hormone known for its antitumour efficacy as it significantly increases the duration of the cell cycle of human breast cancer cells [104]. Given a heterogeneous cohort of individuals with varying degrees of tumour growth rates, our model suggests that the driver of these differences is most likely the cell cycling rate. Drugs targeting this should, therefore, be considered.
Depending on the cancer type, administering an intratumoural injection of a drug can be extremely difficult and administering treatments on the periphery can be an easier course of action. Simulating the model, we found that intratumoural administration of gemcitabineloaded fibres significantly outperforms peritumoural administration. However, there is a threshold distance from the tumour to achieve an effective treatment, beyond which placing fibres further into the tumour bulk sees no added benefit. There is a clear benefit to increasing the dosage multiplicity and spreading the administered drug out amongst the tumour compared to a single high dose. Tumour volume was most significantly decreased when four freedrug point injections were administered compared to a single free-drug point injection. This proposes the existence of a potential threshold above which increasing the multiplicity of dosages or dosage size has a negligible effect over spreading out the dosages.
The location of the fibre and the total drug concentration in the fibre was a major driver of tumour eradication. For fibres located within the centre of the tumour with a significantly high drug concentration, it was possible to completely eradicate the tumour. Moving the fibre farther away from the centre, we found that there was no concentration of drug that would inhibit growth. This suggests that a large amount of drug from the implants is lost to the surrounding tissue, and this has detrimental effects on the efficacy of these devices. Fortunately, simulations show there is a minimal concentration of drug necessary for stabilisation, allowing these predictions to be used a way to guide dosage so that toxicity is minimised, and efficacy is maximised.
Interestingly, we found that changing the release profile of the drug from the device had a similar effect for fast enough release curves, i.e. γ between 10 −3 and γ −1 . This suggests that there is some margin of error for the creation of these devices and provides some promise for their optimisation.
More recently, research has been focused on combining gemcitabine with other drugs to improve its efficacy in the sustained-release devices. For example, nanoparticle albuminbound paclitaxel (nab-paclitaxel) administered in combination with gemcitabine [9] is one of the standard of care treatment regimens that has shown an increase in overall survival in patients with advanced PDAC, as shown in a Phase I/II clinical trial [9]. A phase III clinical trial showed that gemcitabine and erlotibin also significantly increased overall survival in advanced PDAC patients compared to gemcitabine alone [105,106]. This drug combination has been investigated in an analogous sustained-release system by Wade et al. [13] and using this model, future work hopes to examine how the release profile of combination drugs such as this might be optimised.
Due to wanting to reduce the computational complexity of the VCBM, we made some simplifying assumptions that have introduced limitations into our model. To avoid simulating excessively large numbers of cells, we have chosen to scale the spatial unit appropriately so that we simulate on the order of~10 6 cells. An improvement for this model, could be to parallelise the agent update step to increase the speed of the simulation. In addition, we consider only a 2-dimensional cross section of the tumour, which is a simplification given tumour's grow in 3-dimensional environments. While this is a simplification, since the measurements obtained by Wade et al. [13] only measure tumour dimensions in a 2-dimensional cross-section, i.e. width and length, we feel it is appropriate to model growth in 2 dimensions. In addition, since the large effects felt in 3 dimensions will be based on surrounding organ tissue, we feel that since we model neighbouring tissue as having a homogenous effect on tumour growth, there would be no significant impact of extending our model to 3 dimensions. Lastly, we model cell uptake by point sink terms; however, a cell would uptake drug across its surface area through drug molecule binding and internalisation. It would be possible to model this by extending the framework from a single point uptake to a uniform uptake across a cell's defined Voronoi cell region.
There are considerable avenues for future extensions of this work, and we feel the platform we have built is easily extendable by other computational oncologists. In particular, future modelling could extend the model to account for the dense fibrotic nature of PDAC [101,102] and investigate the impact the release and delivery of drug. In addition, the model could be used to simulate the efficacy of dual drug-loaded polymer and verify whether improvements on the current treatment protocol exist. There are many applications of degradable polymeric drug delivery systems in cancer therapy [10], for example, Rezk et al. [10] developed a pH-sensitive polymeric carrier to study the local delivery of anticancer drug bortezomib. They fitted the release profile of the drug from their carrier system to a mathematical formalism. Using our pancreatic cancer growth VCBM, it would be possible to feed in their drug release mechanism and simulate the efficacy under alternative protocols and predict the remaining tumour volume. Lastly, while we did not consider gemcitabine resistance in our model, it does occur in PDAC [3,4]. A simple extension of the model could consider the impact of resistance on the performance of therapy like other works on resistance of chemotherapeutics using mathematical models [107,108].

Conclusion
Treatment for cancers with a poor prognosis, such as PDAC, are in vital need of novel therapeutic approaches that provide sustained, heightened, localised drug concentrations. The computational platform developed in this work can recapitulate spatially heterogeneous tumour growth and treatment with the chemotherapy drug gemcitabine. Investigating the efficacy of gemcitabine released from a degradable polymeric fibre implant, we can suggest that a minimum dosage for maximum efficacy exists based on the location of the device within the tumour. Furthermore, certain release profiles are more effective than others, suggesting that the way in which drug is released from these devices is crucial to improving patient treatment. Moving forward, a study of this form could be used to help inform experimental design and be integrated into future device development.