Bridging in vitro and in vivo research via an agent-based modelling approach: predicting tumour responses to an ATR-inhibiting drug

Translating quantitative information between in vitro and in vivo research remains a scien-tiﬁcally and ﬁnancially challenging step in preclinical drug development processes. However, well-developed in silico tools can be used to facilitate this in vitro to in vivo translation, and we here propose using a data-driven, agent-based model to bridge the gap between in vitro and in vivo research. The agent-based model used in this study is governed by a set of empirically observable rules, and by adjusting only the rules when moving between in vitro and in vivo simulations, whilst keeping the fundamental mathematical model and parameters intact, the agent-based model can ﬁrst be parameterised by in vitro data and thereafter be used to predict in vivo treatment responses. As a proof-of-concept, this modelling approach is here validated against data pertaining to LoVo cells subjected to the ATR (ataxia telangiectasia mutated and rad3-related kinase) inhibiting drug AZD6738, but the modelling approach has the potential to be expanded to other applications. In this article we also highlight how agent-based models, that are currently underutilised in pharmaceutical contexts, can be used in preclinical drug development.


Bridging in vitro and in vivo research
Mathematical models, and corresponding in silico tools, can be used to simulate both in vitro and in vivo scenarios that involve cancer cell populations, or tumours, and their responses to anti-cancer treatments [1,2,3,4,5]. However, cancer cells in an in vitro cell culture experience a microenvironment that is significantly different from the microenvironment experienced by cancer cells in a solid tumour in vivo. As these microenvironments influence cell proliferation and the delivery of oxygen, drug and nutrient molecules to cells, it follows that the dynamics of a cancer cell population in vitro differs from the dynamics of a solid tumour in vivo. Consequently, translating data obtained by in vitro experiments into quantitative information that can guide or predict in vivo experiments remains a challenging, but important, step in drug development processes.
Agent-based models (ABMs) are used in many applications in mathematical biology but are underutilised in the context of pharmaceutical drug development [6]. An ABM consists of multiple, distinct agents that may interact with each other and with their microenvironment. In this study, we introduce a novel modelling approach that uses an agent-based mathematical model to bridge the gap between in vitro and in vivo research, as is conceptually illustrated in Figure 1. In the ABM at the core of this modelling approach, an agent consists of one cancer cell or a group of cancer cells, where the behaviour and fate of each agent is governed by a set of empirically observable and well-established modelling rules that incorporate both intracellular and microenvironmental dynamic variables, as is thoroughly described throughout Section 2. To account for differences between in vitro and in vivo scenarios, the modelling rules are appropriately adjusted when moving between in vitro and in vivo simulations. By only adjusting the rules, whilst keeping the fundamental mathematical model and parameters intact, when moving between in vitro and in vivo simulations, the mathematical framework can first be parameterised by in vitro data and thereafter be used to predict in vivo treatment responses. Figure 1: A schematic of the mathematical modelling approach used in this study. An agentbased mathematical model, that distinguishes between in vitro and in vivo modelling rules, is used to bridge the gap between in vitro and in vivo research. The mathematical model is first parameterised by in vitro data and is thereafter used to predict in vivo outcomes.
As a proof-of-concept of this modelling approach, we here simulate LoVo (human colon carcinoma) cells subjected to the anti-cancer drug AZD6738. The in vitro and in vivo data used in this study are gathered from previous work by Checkley et al. [7]. The ABM used in this study 2 is an extension of a cellular automaton (CA) introduced by Powathil et al. [8].

DNA damage response inhibiting drugs
The deoxyribonucleic acid (DNA) in human cells is perpetually exposed to, potentially harmful, influences that can be derived from both exogenous and endogenous sources and events [9,10]. Exogenous sources include ultraviolet radiation, ionising radiation and chemotherapeutic drugs, whilst erroneous DNA replication is an example an endogenous event yielding DNA damage [9]. Regardless of the source, a multitude of intracellular events are triggered when the DNA in a cell becomes damaged. Cells may, for example, respond to DNA damage by activating DNA repair mechanisms, cell cycle arrest or, in cases of severe DNA damage, apoptosis [11]. Such cellular responses to DNA damage are mainly governed by the DNA damage response (DDR), which comprises a complex network of signalling pathways [11]. The DDR has many functionalities and, amongst other things, it monitors DNA integrity and repairs DNA damage in order to maintain genomic stability in cells. The DDR also governs DNA replication, cell cycle progression, and apoptosis [9,12]. When DNA repair in a cell is needed, the DDR activates relevant effector proteins [9]. Included in the group of DDR-associated effector proteins are approximately 450 proteins [12], out of which the two main regulators for cell cycle checkpoints are ataxia telangiectasia mutated kinase (ATM) and ataxia telangiectasia mutated and rad3related kinase (ATR) [10]. ATM and ATR belong to the enzyme family phosphatidyilinositol-3-OH-kinases (PI3K), and they both play central roles when cells respond to DNA damage [11]. In this work, we study the effects of an anti-cancer drug, namely AZD6738, that works by inhibiting ATR activity.
DNA lesions in form of single-strand breaks are a common result of replication stress, and the repair of single-strand DNA breaks is mainly attributed to ATR activity. A drug that inhibits ATR activity consequently inhibits the repair of single-strand DNA breaks post replication stress. Cancer cells are associated with high replication stress and consequently ATR inhibitors have, during the last decade, been explored as anti-cancer agents [9,11,13]. With the premise that inhibiting DNA damage responses should increase the effect of some other main therapy, DDR inhibitors have been explored as both radiotherapy and chemotherapy treatment intensifiers [11,13]. Two well-studied ATR inhibitors are AZD6738 and VX-970. AZD6738 is an oral ATR inhibitor, and its anti-tumour potential has been demonstrated in preclinical vitro and in vivo xenograft studies for various ATM deficient cell lines, including ATM deficient lung cancer, chronic lymphocytic leukemia and metastatic adenocarcinoma of the colon [7,10,14]. Combination treatments that combine AZD6738 with either radiotherapy or chemotherapy have produced synergistic results in preclinical settings [10], and AZD6738 is currently being evaluated in clinical phase I/II trials [9,13]. VX-970 is an intravenous ATR inhibitor [15] that has demonstrated tumour controlling effects in a phase I clinical trial, both as a monotherapy and in combination with the chemotherapy drug carboplatin [9]. A summarising table of clinical trials involving ATR-inhibitors can be found in an article by Mei et al. [13].

Model and Method
An ABM, specifically a CA, is used to model a population of cancer cells (in vitro), or a solid tumour (in vivo), that evolves in time and two spatial dimensions. The model describes the behaviour of cancer cells using a set of modelling rules. In order to account for differences between in vitro and in vivo scenarios, these rules are appropriately adjusted when moving between in vitro and in vivo simulations, as is described throughout Section 2. Taking a minimal parameter approach, we aim to use as few rules and parameters as possible to capture the nature of the regarded system. We here chose to include model rules and parameters that pertain to the cells' doubling time and cell cycle (see Section 2.2), cell proliferation on the lattice (see Section 2.3), the distribution of oxygen and drugs across the lattice (see Sections 2.4 and 2.5 respectively) and cellular responses to local oxygen and drug concentrations (see Sections 2.4 and 2.6 respectively). In this work, details concerning nutrient distribution and its effect on tumour growth are not included. Instead, under a simplifying modelling assumption, the diffusion of oxygen forms a surrogate for the distribution of nutrients. Differences between the in vitro and the in vivo modelling rules are pictorially summarised in Section 2.8, and in vitro-calibrated model parameters are listed in Table 1.
The in vitro and in vivo data used in this study are gathered from previous work by Checkley et al. [7]. In the regarded in vitro experiments, populations of LoVo cells were plated and subjected to AZD6738, where population sizes of up to roughly 4000 cells were reported [7]. In the in vivo experiments, LoVo cells were subcutaneously injected in flanks of female Swiss nude mice in order to produce human tumour xenografts, and AZD6738 treatments started when the tumours had reached a volume of 0.2-0.3 cm 3 [7]. Here, we regard treatment responses in terms two dynamic variables: population or tumour size and percentage of DNA-damaged (i.e. γH2AX-positive) cells. The in vitro and in vivo data used in our current study are available in the Supplementary Material (SM).

Cellular automaton lattice
In the model, one agent corresponds to one cancer cell (in vitro) or, due to computational costs, one group of cancer cells (in vivo). The behaviour and fate of each agent is governed by a set of rules that incorporate both intracellular and environmental dynamic variables using multiscale modelling techniques [16]. At the start of an in silico experiment, one initial agent is placed in the centre of the lattice. This initial agent produces daughter agents and ultimately gives rise to a heterogeneous population of agents. When the population has reached an appropriate size (chosen to match the in vitro and in vivo data), AZD6738 anticancer treatments commence. The CA lattice is a square lattice, and every lattice point is either empty or occupied by one agent. If a lattice point is empty, it consists of extracellular solution (in vitro) or only extracellular matrix (ECM) (in vivo). The ECM comprises multiple components such as collagen, elastin and fibronectin but, as under a simplifying modelling assumption, we do not distinguish between these components in the model [17]. In the in vitro simulations, the dispersion of any molecules across the lattice is modelled as instantaneous, and thus the extracellular solution is considered to render the entire lattice homogeneous in terms drug and oxygen concentrations at all times. In the in vivo simulations, however, drug and oxygen molecules are modelled as gradually diffusing over the tumour and its local environment, and consequently the in vivo lattice will be heterogeneous in terms of drug and oxygen concentrations. Oxygen and drug distribution across the lattice are further discussed in Sections 2.4 and 2.5 respectively. Differences between the simulated in vitro an in vivo lattices are described below.
In vitro lattice: Cell populations evolve on a two-dimensional square lattice with 100 × 100 lattice points, where the spacing in both spatial directions, x 1 and x 2 , corresponds to one cell diameter.
In vivo lattice: Approximating the tumour as spherical, we simulate (only) a central cross section of the tumour as an, approximately circular, disk of cells living on a two-dimensional square lattice. This lattice is specifically anL ×L = 500 × 500 square lattice, with a spacing in both spatial directionsx 1 andx 2 equal to 40µm. The dimensions are chosen in order to allow our agent-based model to simulate the required physical dimensions, whilst keeping computational costs low. Post simulation time, the two-dimensional cross section of cells is extrapolated to represent a three-dimensional tumour-spheroid. This disk-to-spheroid extrapolation process is outlined in the SM.

Cell cycle model
In order to capture the influence of ATR and the ATR inhibitor AZD6738 on the cell cycle, we use a probabilistic, rule-based cell cycle model adapted from previous mathematical (nonagent-based) work by Checkley et al. [7]. As illustrated in Figure 2, this cell cycle model can be represented as a graph with nodes (cell cycle phases or states) that are connected via various paths (phase/state transitions). A cell can be in an undamaged state (G1, S or G2/M), a replication stress-induced DNA damaged state (D-S) or a dead state, where the cause of cell death is unrepaired replication stress. As ATR is active in the checkpoint in the intra-S phase of the cell cycle, both under undamaged circumstances and in response to DNA damage [11], ATR inhibition will inhibit the cell from progressing to the G2/M state in the mathematical cell cycle model. A cell can take different possible paths through the cell cycle graph, and every time that paths fork, random number generation (from a uniform distribution) determines which path will be taken. Every cell commences its life in the G1 state, but thereafter a cell can enter either the S state or the damaged S (D-S) state. The probability that a cell enters the D-S state is denoted Π D−S and is calibrated by in vitro data [7]. If a cell enters the D-S state, it has a chance to repair itself and enter the S state. If there is no drug in the system, this repair is always achieved, however the repair path is inhibited by the presence of the drug AZD6738. The higher the drug-concentration, the more unlikely it is that a cell in the D-S state will successfully repair itself to the S state. If a cell in the D-S state fails to repair, it is sentenced to die. Whether a cell in state D-S repairs or dies is decided by comparing a random number, generated from a uniform distribution, to the cell's survival probability, which is influenced by the local drug concentration C(x, t), as is described in detail in Section 2.6. A cell that has successfully reached the S state continues to the G2/M state, after which it duplicates and starts over in the G1 state again, ready to perform another cell cycle.
In order to allow for asynchronous populations, each agent i on the lattice is assigned an shown as nodes in the graph. Viable (undamaged or damaged) states are shown in circles, whilst the dead state is shown as a cross. Paths illustrate transitions between states, and symbols next to the paths denote the probabilities that the corresponding paths will be taken. The dashed path can be inhibited by an ATR inhibiting drug, such as AZD6738.
individual doubling time τ i , where τ i is a random number generated from a normal distribution with mean value µ and standard deviation σ. Each agent is attributed an individual cell cycle clock, that determines when the agent should progress to a subsequent state in the cell cycle model. Progression to a subsequent state occurs once an agent has spent a certain fraction of its doubling time in its current state. The fraction of the doubling time spent in the G1, S (including D-S) and G2/M states are respectively denoted Θ G1 , Θ S and Θ G2/M , where these values are approximate and chosen from literature to match values for typical human cells with a rapid doubling time of 24 hours so that Θ G1 = 11/24, Θ S = 8/24 and Θ G2/M = 5/24 [18]. The fraction of an agent's doubling-time spent in the D-S state, Θ D−S , is on the other hand fitted by in vitro data produced by Checkley et al. [7], as outlined in the SM. In vitro and in vivo cell cycle modelling rules are described below.
In vitro cell cycle model rules: One agent corresponds to one cancer cell that is assigned an individual doubling time τ i . The cell cycle path taken by cell i is governed by random number generations specific to that cell.
In vivo cell cycle model rules: One agent comprises a group of identical cancer cells. Each agent is assigned an individual doubling time, τ i , and thus all cells belonging to agent i progress simultaneously and uniformly through the cell cycle model. Random number generations specific to agent i determine which path the agent takes through the cell cycle.

Cell proliferation
When an agent has completed the mitosis state in the cell cycle model, a secondary agent, namely a daughter agent, is produced. Each daughter agent is placed on a random lattice point in the (approximately circular) neighbourhood of its parental agent. To accomplish circular-like growth, the model stochastically alternates between placing daughter agents on Moore and von Neumann neighbourhoods of parental agents, as is pictorially described in the SM. A daughter agent is allowed to be placed on, up to, a νth order neighbourhood of its parental agent, but lower order neighbourhoods (i.e. neighbourhoods closer to the parent) are prioritised and populated first. Modelling rules concerning in vitro and in vivo cell proliferation are outlined below.
In vitro proliferation rules: In the experimental in vitro setup, there is no spatial constraint or nutrient deficiency that is inhibiting cell division within the time-course of the experiment. Consequently cells are allowed to divide freely in the in vitro model and we set ν to be equal to infinity in the in vitro case (with the restriction that agents can not be placed outside the lattice in the in silico implementation).
In vivo proliferation rules: In vivo tumours typically consist of a core with quiescent cells and a shell of proliferating cells. To accommodate for this, a daughter agent (representing a group of daughter cells) is allowed to be placed on up to a third order (approximately circular) neighbourhood of its parental agent, so thatν = 3, in accordance with previous mathematical models [8]. For the in vivo experiment regarded in our current study,ν = 3 matches the experimental data. However, for other experiments, the value ofν may be adjusted to fit the specific cell-line and modelling scenario at hand. When an agent is in the G1 phase of the cell cycle, it scans its environment to see if it has enough resources, in terms of space and nutrients, to commence the process of producing a daughter cell. If not, the cell enters the quiescent phase [19]. Thus in the model, when an agent is in the G1 phase, it continues to progress through the cell cycle model, provided that some free space is available on the lattice within its νth neighbourhood. If this is not the case, the agent exits the cell cycle to enter a quiescent state G0. Should neighbourhood space be made available again, here as a result of anticancer targeting, quiescent agents may re-enter the cell cycle.

Oxygen distribution and influence on cells
Tumour growth and treatment responses are highly influenced by intratumoural oxygen levels [20,21,22] and severely hypoxic (cancer) cells proliferate slower than do well-oxygenated cells [19].
In vitro oxygen distribution and responses: In the mathematical in vitro model, all cells are assumed to be well-oxygenated in accordance with the experimental in vitro setup performed by Checkley et al. [7]. Consequently, neither oxygen dynamics nor cellular responses to oxygen levels are incorporated in the in vitro model.
In vivo oxygen distribution and responses: Within solid tumours, oxygen concentrations typically vary and hypoxic regions are common tumour features [22,23,24]. Avoiding complicated details of vasculature in the model, we here approximate oxygen as diffusing in from 'outside the boundary of the tumour'. Oxygen dynamics across the CA lattice is here described using a mechanistic diffusion equation, where the oxygen concentration in locationx at time t is denoted by K(x, t) where The first term in Equation 1 describes oxygen diffusion across the CA lattice, the second term is an oxygen supply term and the third term describes oxygen uptake by cells. Accordingly, D K (x, t) is the oxygen diffusion coefficient, and r K and φ K are supply and consumption coefficients respectively. The diffusion coefficient for oxygen is known from literature to be 2.5 × 10 −5 cm 2 s −1 [8]. Assuming that oxygen diffuses slower inside the tumour than outside the tumour, the oxygen diffusion coefficient is divided by a factor 1.5 if there is a cell in locationx at time t. The binary factor m(x, t) is 1 if the regarded locationx is outside the tumour boundary at time t and 0 otherwise, i.e. m(x, t) is 1 if the regarded lattice point is not occupied by an agent nor completely surrounded by agents. Similarly, the binary factor cell(x, t) is 1 if there is a viable (non-dead) cell in locationx at time t, and 0 otherwise [8]. Equation 1 is coupled with no-flux boundary conditions, thus the total amount of oxygen in the system will fluctuate over time [25]. A scaled oxygen variableK(x, t) is introduced in order to express oxygenation levels in percentages (%) between 0% and 100%. This scaled oxygen value is computed at every unique time step t u by where maxx ,tu K(x, t u ) denotes the maximum occurring K(x, t u )-value at the time point t u and h is a scaling factor that is included in order to classify cells with oxygen levels of 10%, or less, as hypoxic [8,23]. Low cellular oxygen levels have been shown to delay cell cycle progression by inducing arrest in, particularly, the G1 phase of the cell cycle [19]. Consequently, in our model, hypoxic cells display arrest (i.e. delay) in the G1 phase of the cell cycle. In mechanistic Tyson-Novak type cell cycle models [26,27,28], the cell cycle is governed by a system of ordinary differential equations (ODEs) in which the G1 phase can be inherently elongated under hypoxic conditions by incorporating hypoxia-induced factors into the ODEs [8]. In the mathematical model discussed in this paper, however, we use agent-attributed clocks to model cell cycle progression and thus, in order to achieve a longer G1-phase under hypoxic conditions, we introduce a G1 delay factor (G1DF) where otherwise. ( The G1DF is an approximation for how much the G1 phase is expanded in time as a function of cellular oxygen concentration. It is matched to fit data points extracted from a previous mathematical study by Alarcon et al. [19], in which a Tyson-Novak cell cycle model is extended to incorporate the action of p27, a protein that is up-regulated under hypoxia and delays cell cycle progression. Data-fitting yields the parameter values a 1 = 0.9209, a 2 = 0.8200, a 3 = −0.2389 [23]. Thus the fraction of an agent's doubling time spent in the G1 state is G1DF (K(x, t)) · Θ G1 , where G1DF (K(x, t)) = 1 for normoxic cells.

Drug distribution across the lattice
Drug distribution significantly varies between in vitro and in vivo settings. In the regarded in vitro setup, the drug concentration can be regarded as homogeneous, whilst heterogeneous drug concentrations must be accounted for in vivo.
In vitro drug distribution: In the in vitro experiments performed by Checkley et al. [7], plated cell populations of roughly 1000 cells were treated with AZD6738 in the solvent dimethylsulfoxide (DMSO). In the mathematical model, we approximate the drug distribution across the CA lattice to be instantaneous (occurring at treatment time T 0 ) and homogeneous. We furthermore assume that the drug has a half-life time that exceeds the time course of the experiment, and note that there is no other drug elimination from the in vitro system. In our mathematical model, this is equivalent to there being no drug decay or elimination, hence the drug concentration C(x, t), in locationx at time t is simply given by where C denotes the applied drug concentration (in units of molarity).
In vivo drug distribution: In the in vivo experiments performed by Checkley et al. [7], AZD6738, or vehicle in the control case, were administered via oral gavage once per day to female Swiss nude mice. In the mathematical in vivo model, we consider the drug to diffuse through the tumour from its surrounding, creating a drug gradient within the tumour. In the mathematical model, this drug dynamics is modelled using a partial differential equation (PDE), where the concentration of AZD6738 at locationx at time t is denoted by C(x, t) such that where D AZD is the diffusion coefficient of the drug AZD6738, and the supply coefficient p(x, t) is greater than zero at drug administration times only for lattice points outside the tumour. Assuming first order kinetics for drug elimination, the drug decay constant η AZD is matched to the reported half-life time of 6 hours for AZD6738 in vivo [29,30]. Note that the drug decay term here represents all drug elimination from the system, both metabolic and that caused by excretion.
The diffusion rate of a drug is predominantly affected by the molecular size of the drug. More specifically, the diffusion coefficient of a drug is inversely proportional to the square root of the molecular weight of the drug, so that large molecules diffuse more slowly than do small molecules [30]. Using this assumption, the drug diffusion coefficient D AZD is set in relation to the oxygen diffusion coefficient D 0 2 , as is done in previous mathematical studies [8]. Thus the relationship between the diffusion coefficients corresponds to the square of the inverse relationship between the molecular weights, such that where the molecular weights are collected from the PubChem database [31]. Details regarding pharmacokinetics are outside the scope of this study, bioavailability is instead calibrated using extreme case drug scenarios, as described in the SM.

Drug responses
AZD6738 inhibits DNA repair from the D-S state to the S state in the cell cycle model, as illustrated in Figure 2, and, in our model, maximal drug effect corresponds to complete repair inhibition. The drug effect is modelled using an agent-based adaptation of the sigmoid Emax model [32], in which the drug effect on a cell in positionx at time t is given by E max denotes the maximal drug effect, here corresponding to complete repair inhibition (E max = 1), EC 50 denotes the drug concentration required to achieve half of the maximal drug  effect (0.5 · E max ) and γ is the Hill-exponent [32]. EC 50 and γ are fitted from the in vitro data, as outlined in the SM. When an agent is scheduled to progress from the D-S state in the cell cycle, it has a probability Π rep ∈ [0,1] to repair, where Π rep is determined by the local drug concentration so that Note that in the absence of drugs, the repair probability is 1. When a cell dies, it is transformed into a membrane-enclosed 'cell-corpse' [30]. In the in vivo setting, this cellular debris is digested by macrophages but in the in vitro setting such 'cell-corpses' linger on the lattice during the course of the experiment. Post the lethal event (i.e. the D-S to S repair failure) a cell is declared 'dead' in the model after a time T L→D has passed (where L stands for 'lethal event' and D stands for 'death'). The parameter T L→D is calibrated by in vitro experiments. The differences between modelling rules for in vitro and in vivo drug responses are described below.
In vitro drug responses: After failure to repair from the D-S state, a cell (i.e. and agent) is considered to be dead after a time T L→D has passed. However, a dead cell is never physically removed from the lattice.
In vivo drug responses: An agent (i.e. a group of cells) is declared to be dead and removed from the lattice after an amount of time T L→D post the lethal event (failure to repair).

Parameters
The parameters used in the mathematical model are calibrated by in vitro data, this calibration process is described in the SM. In the context of quantitative pharmacology, knowledge about a model's robustness is crucial [33], therefore we have provided results from the uncertainty and sensitivity analysis in the SM. We performed three different uncertainty and sensitivity analyses techniques, suitable for agent-based models with stochastic elements, namely (i) consistency analysis, (ii) robustness analysis and (iii) Latin hypercube analysis [34,35]. Detailed descriptions on how to perform and interpret these techniques are available in an introductory uncertainty and sensitivity analyses review [34]. In accordance with the performed consistency analysis, we run 100 simulations per in silico experiment in order to formulate results (in terms of mean values and standard deviations) that mitigate uncertainty originating from intrinsic model stochasticity.

Differences between in vitro and in vivo modelling rules
Differences between the in vitro and in vivo modelling rules in the ABM are pictorially summarised in Figure 3. A note on the simplifying modelling assumptions that we have used in this study is provided in the SM.

Implementation
The mathematical model is implemented in an in-house C++ [36] framework, in which PDEs are solved using explicit finite difference methods. Simulation cell-maps are visualised using ParaView [37] and uncertainty and sensitivity analyses are performed in MATLAB [38].

Results
The mathematical framework is first calibrated by in vitro data, and is thereafter used to predict treatment responses in human tumour xenografts in vivo.

Simulating in vitro experiments
The mathematical framework is parameterised using in vitro data produced by Checkley et al. [7]. In the in vitro experiments, populations of LoVo (human colon carcinoma) cells were exposed to the ATR inhibiting drug AZD6738. The in silico results in Figure 4 show the evolution of the in vitro cell population over time in terms of percentage of DNA damaged, i.e. γH2AX-positive, cells (Figure 4 Left) and in terms of cell count (Figure 4 Right). AZD6738 drugs are given at 0 hours, when the cell population has reached a size of approximately 1000 cells. Simulated response curves for six different drug concentrations, including the zero-drug concentration control case, are shown. Also shown in Figure 4 are simulation standard deviations and in vitro data [7].
Using a minimal-parameter modelling approach, the mathematical framework is calibrated to fit in vitro data points without introducing any variable model parameters. This calibration process is described in the SM. Our results demonstrate that, post in vitro parameterisation, our mathematical framework is able to capture the qualitative nature of in vitro LoVo cell population growth and drug (AZD6738) responses. Figure 4 (Left) demonstrates that the model is able to qualitatively reproduce the asymptotic fraction of DNA damaged cells in the system but fails to match early in vitro data points. The sensitivity analysis demonstrates that the treatment timing (in relation to the overall cell cycle phase composition of the cancer cell population) notably influences treatment responses in terms of percentage of γH2AX-positive cells. This indicates that in the in vitro experiments, a large fraction of the cells would have been in the DNA damaged D-S state when, or shortly after, the drug was applied. The experimental error bars in Figure 4 (Right) and the numerical cell count data, available in Table 1 in the SM, demonstrate that the doubling time of the cell population drastically decreased towards the end of the in vitro experiment and, consequently, our agent-based model was not able to replicate cell count data at 72 hours as the modelling rules and parameters were not updated over time.

Simulating in vivo experiments
Post in vitro calibration, the mathematical framework is used to simulate the in vivo experiments performed by Checkley et al. [7] in which LoVo xenografts, that are injected in mice flanks, are treated with AZD6738 once daily for 14 days. The results in Figure 5 show AZD6738 drug responses in terms of the percentage of DNA damaged (γH2AX-positive) cells ( Figure 5 Left) and tumour volume (Figure 5 Right). Simulated response curves to three different drug doses (0, 25 and 50 mg/kg) and in vivo data are provided in Figure 5. Figure 5 (Right) demonstrates that our simulated results qualitatively agrees with the in vivo results reported by Checkley et al. [7] for approximately 12 days post tumour injection for control case tumours, and for approximately 8 days post tumour injection for tumours subjected to drugs. This can be explained by the fact that the behaviour of the agents in our current model does not change over time, when in fact tumours are highly adaptable and responsive to external pressures. It follows that details pertaining to tumour growth and drug sensitivity may vary over time, and in future work the agent-based model used in this study can be updated to incorporate variable modelling rules and parameters.

Discussion
In silico results obtained in this study were compared to in vitro and in vivo data and can, furthermore, be compared to previous mathematical modelling results produced by Checkley et al. [7]. In their study, Checkley et al. [7] modelled tumour responses to AZD6738 using coupled ordinary differential equations, where a pharmacokinetic/pharmacodynamic (PK/PD) model of tumour growth was integrated with a mechanistic cell cycle model. Their model is predictive of in vivo xenograft studies and is being used to quantitatively predict dose and scheduling responses in a clinical Phase I trial design [7]. Our modelling results qualitatively agree with those produced by Checkley et al. [7], although two different modelling approaches have been taken: Checkley et al. [7] regard the tumour as one entity with different compartments whilst we here use a bottom-up modelling approach and regard the tumour as consisting of multiple, distinct agents.
Moving drug-response investigations from in vitro to in vivo settings is a key step involved in the process of moving a drug from bench-to-bedside. However, in vivo data are often sparse, as gathering in vivo data is associated with practical, financial and ethical constraints. Plentiful and adaptable in silico data are, on the other hand, easy to produce, and can thus be used as epistemic complements to sparse in vivo data. Well-formulated in silico tools can be extended to investigate various dose-schedule scenarios in order to guide in vitro and in vivo experiments. Such in silico experiments may provide a testbed for simulating various mono and combination therapies. We here propose creating ABM in silico tools in which modelling rules are based on "fundamental" principles that describe how cancer cells in a system behave (where it is up to the modeller to decide which principles should be considered "fundamental" in the specific modelling scenario at hand). The ABM considered in this study is an extension of a mathematical model that has previously been used to study tumour growth and treatment responses to chemotherapy, radiotherapy, hyperthermia and hypoxia-activated prodrugs [8,39,40,41,42,23]. In recent years, several ABMs have been developed for the purpose of describing various aspects of cancer dynamics [43], and it should be noted that the modelling approach proposed in Figure 1 is not conceptually limited to usage with the ABM described in this study. The choice of ABM should be influenced by the research question at hand, the desired level of model details and the available data. Examples of data-driven ABMs are available in a recent review article by Chamseddine and Rejniak [44] discussing hybrid models, and hybrid modelling techniques, used in the field of mathematical oncology today.
Data-driven modeling, exploitation of existing data and proof-of-concept studies are important steps involved in current and future procedure for enabling mathematical modeling in systems medicine, as argued in a report by Wolkenhauer et al. [45]. Despite the fact that mathematical modelling is becoming increasingly popular in the pharmaceutical industry, there are not that many ABMs present in the pharmaceutical scene [6]. We argue that this is a missed opportunity in the context of oncology, as ABMs naturally capture the heterogeneous nature of tumours, which is known to complicate treatments. As multiscale ABMs organically enable the integration of data across various scales in time and space, it follows that they are useful to the interdisciplinary team that wishes to combine data and knowledge from its team members. Following interdisciplinary collaborations between clinicians, biologists and mathematicians, mathematical modelling may be used to enable in silico informed drug development.

Competing Interests
No competing interests to declare.

Experimental Data
The experimental in vitro and in vivo data used in our current study are gratefully gathered from a previous study performed by Checkley et al. [1]. In vitro data are listed in Table 1 and in vivo data are listed in Tables 2 and 3.

Model Parameterisation
Using a minimal-parameter approach, seven model parameters are calibrated using the in vitro data previously produced by Checkley et al. [1], as listed in Table 1 in Section 2.7 in the main paper. Parameter sensitivity is explored in the sensitivity analysis (Section 5 of the supplementary material). The calibration process is outlined in Sections 2.1 through to 2.4. The in vivo calibration is described in Section 2.5.

Cell doubling
In the model, the doubling time of a cell i is denoted τ i , where τ i is stochastically picked from a normal distribution with mean value µ and standard deviation σ. Thus µ corresponds to the average cell doubling time and σ corresponds to how synchronised the cells are. If σ is zero, then all cells have perfectly synchronised cell cycles and duplicate at the same time. Higher σ values achieve less synchronised cell cycles amongst cells and smoother cell count growth curves over time. The control case (i.e. no drug) cell count data is used to estimate µ and σ.
By observing the control case cell count data in Table ??, we note that the cell population roughly doubles between 2 and 24 hours (indicating an average cell doubling time, or µ, of approximately 22 hours, which is less than 24 hours). Furthermore, the cell population also roughly doubles between 24 and 48 hours (indicating that µ is approximately 24 hours). However, in the last 24-hour interval, between 48 and 72 hours, the control population increases by less than 5% (indicating that µ is more than 24 hours). From these three observations, we choose to make the modelling assumption that the average doubling time for cells should be around 24 hours, and σ-values in the parameter range [22,26] hours are investigated in silico. Due to the synchronised nature of the cell count data, σ-values between 0 and 2.5 hours were investigated in silico, where σ = 0 h corresponds to completely synchronised cells and σ = 2.5 h achieves a smooth cell count growth curve. After an iterative process of tuning parameters and running in silico experiments, the calibrated values are set to be µ = 24 hours and σ = 0.5 hours. As discussed in the main article, the ABM can be improved to better fit wet lab data by including variable parameter values or rules, that are updated over time. However, in the current stage of our work, we decided to fix µ and σ.

Cell cycle progression
The in vitro data provides information on how many cells are in the damaged S state via the biomarker γH2AX. For the control case, the number of γH2AX positive cells in our mathematical model depends on two variables: (1) the probability (Π D−S ) that a cell enters the D-S state and (2) the amount of time (Θ D−S · τ i ) spent in the D-S state prior to repairing. Recall that Θ D−S is the fraction of a cell's doubling time (τ i ) spent in the D-S state. As a first step, in silico experiments are performed in which we find various parameter pairs (Π D−S , Θ D−S ) that agree with the control data. We thereafter note that the in vitro drug effect saturates for concentrations 3, 10 and 30 µM and assume that the maximal dose (30 µM) yields 100% D-S to S repair inhibition. Thus a second step we test the variable pairs (Π D−S , Θ D−S ) for this 'maximal drug and no repair' scenario in silico, and we match these in silico results to the 30 M in vitro data. Here, we only use data from early time points (time < 12 hours) in order to avoid the influence that dying cells have on the data and model outputs. After iterative in silico testing, the variable pair (Π D−S , Θ D−S ) that best fits these both extreme cases is Π D−S = 0.75 and Θ D−S = 0.03. The first extreme case refers to the 'no drug' in silico experiment matched to the in vitro control data, where we assume that all D-S cells repair to state S. The second extreme case refers to the 'maximum drug' in silico experiment matched to the 30µM control data, where we assume that no D-S cells repair to state S.

Drug response
Drug effects are modelled using the sigmoid E-max model [2], where the drug effect E is a function of the drug concentration C, so that where E max denotes the maximal drug effect. Here we set E max = 1 to corresponds to total D-S to S repair inhibition. EC 50 denotes the drug concentration that achieves half of the maximal drug effect and γ is the Hill-coefficient. If drug effect is plotted over time, the EC 50 -value determines the asymptotic behaviour of the effect whilst the γ-value determines how quickly the asymptotic value is reached.
From the in vitro data, we note that the drug concentration 1 µM achieves roughly half of the total drug effect in terms of γH2AX-positive cells. (Note from Table ?? that when the drug concentration is 10 µM or 30 µM, the percentage of γH2AX-positive cells is roughly 67% at 72 hours, and when the drug concentration is 1 µM, the percentage of γH2AX-positive cells is roughly 33% at 72 hours. Furthermore the lower drug concentration of 0.3 µM yields a percentage of roughly 7% γH2AX-positive at 72 hours, and the higher drug concentration of 3 µM yields a percentage of roughly 64% γH2AX-positive at 72 hours. Consequently, we use 0.3 µM as a lower bound and 3 µM as an upper bound for the parameter range in which we seek EC 50 , and EC 50 values in the range [0.3 µM, 3 µM] are investigated with various Hill coefficients to fit in vitro data for all (non-control) drug concentrations. In order to avoid the impact that dying cells have on the data used parameterise EC 50 and γ, only early in vitro data (time < 12 hours) is used to guide the calibration. After iterative in silico testing, the best variable pair (EC 50 , γ) is determined to be EC 50 = 1 µM and γ = 2.

Cell death
In the in vitro experiments, cells that are damaged (but not yet dead) are γH2AX-positive. In the model, the time it takes between the 'lethal event' (i.e. a cell's failure to repair) and a cell being 'dead' is denoted T L→D and is matched from the in vitro experiment. After noting the asymptotic behaviour of the in vitro data, both in terms of cell damage and cell count, we estimate that the rate of cell elimination should roughly correspond to the rate of cell production, and thus T L→D should be in the same order of magnitude as the doubling time. Consequently, values of T L→D between 0 and 2 τ i are explored in silico after which T L→D = τ i is chosen as it best matches the in vitro data for all tested (non-control) drug concentrations.

In vivo calibration
For the control case, the in vivo model is directly calibrated by the in vitro data, and no further calibration is needed. For drug concentrations larger than 0 µM, we use the in vivo data for the highest administered drug dose to calibrate the model in order to disregard details concerning pharmacokinetics and bioavailability. In future work, our model can be integrated with pharmacokinetic modelling techniques. 5

Proliferation of Agents on the Lattice
When an agent divides, it places a daughter agent on the lattice. To achieve an approximately circular population growth on the square lattice, the placement of daughter agents alternates between following the von Neumann (Figure 1(a)) and Moore (Figure 1(a)) convention. The lattice points that are directly adjacent to the lattice point occupied by the parental agent (i.e. the agent that is dividing and thus producing a daughter cell) constitute the 1st order (von Neumann or Moore) neighbourhood of the parental cell, as is illustrated in Figure 1. If a parent is scheduled to place a daughter agent on the lattice according to the von Neumann/Moore convention, then a daughter agent will be placed on a random lattice point in the 1st order von

Cross-Section to Tumour Spheroid Extrapolation
When implementing our mathematical in vivo model, only a central cross-section of the tumour is actually simulated in silico and post simulation time this cross-section area (that is approximately circular) is extrapolated to a tumour volume (that is approximately spherical). From the extrapolated tumour spheroid, the two outputsX 1 (percentage of γH2AX-positive cells) and X 2 (tumour volume) are gathered. This is done by using simulated areas to compute the total tumour volume, and the quiescent tumour volume, Quiescent Tumour Volume = 4π 3 From the above, the volume of cycling, or proliferating cells, is obtained by Cycling Tumour Volume = Total Tumour Volume − Quiescent Tumour Volume.
Now the outputX 1 can be computed where,

A Note on Simplifying Modelling Assumptions
Note that an agent in the in vivo setting can, in general, be chosen to comprise either one cancer cell or a group of cancer cells. Note also that modelling the in vivo tumour as a (spatially) twodimensional disk means that the distribution of nutrients and drugs is modelled across a twodimensional (rather than a three-dimensional) space. Likewise, the parameterν, that governs the allowed distance between a parental agent and its daughter agents, only concerns proliferation on a two-dimensional plane. By increasingν, the tumour will grow quicker and comprise a higher fraction of cycling (to quiescent) agents. Thus both the parametersν and µ influence the rate of tumour growth in the in vivo simulation, and should ideally be fitted to match detailed in vivo data. In the current study, we choose to keep µ at the in vitro-calibrated value and thereafter fitν to match the in vivo data. A modeller can chose to use a three dimensional spatial domain, and thus explicitly model a tumour spheroid instead of a tumour cross-section. Computational costs, fineness of available data, and the desired level of simulation details should be used to guide the choice of agents and spatial domain.

Uncertainty and Sensitivity Analyses
To evaluate the in silico findings obtained in our in vitro study, three uncertainty and sensitivity analyses techniques are performed. The three techniques are namely: (1) Consistency Analysis, which is used to determine how many in silico runs should be performed before defining results in terms of statistical metrics in order to mitigate uncertainty originating from intrinsic model stochasticity, (2) Robustness Analysis, which investigates model sensitivity to local parameter perturbations and (3) Latin Hypercube Analysis, which investigates model sensitivity to global parameter perturbations. These techniques are thoroughly described in a review/methodology paper [3] which provides background information concerning the origin of these techniques as well as detailed information on how to implement them. To perform uncertainty and sensitivity analyses we need to specify a set of inputs and outputs. Here, the output variables are X 1 : the percentage of γH2AX-positive (i.e. damaged) cells at the end time of the experiment (72 hours), and X 2 : the cell count (i.e. the number of non-dead cells) at the end of the experiment. The input variables are the seven model parameters listed in Table 1, in the main article, that we calibrate using in vitro data. These inputs are namely µ, σ, Π D−S , Θ D−S , EC 50 , γ and T L→D .

Consistency Analysis
Results from the Consistency Analysis are provided in Figures 2, 3, 4, 5, 6 which show thê A-measures, in both computed and scaled forms, for the distribution sizes n = 1, 5, 50, 100, 300 respectively. By observing Figures 2 through to 6, it is clear that the statistical significance decreases with increasing distribution size n, as is shown in Figure 7 and Table 6.1 which show the maximal scaledÂ-values for all tested distribution sizes. These results demonstrate that the distribution size n = 100 is the smallest tested distribution size that yields a small statistical significance (i.e. a maximum scaledÂ-value smaller than 0.56) for both regarded output variables X 1 and X 2 . From this we decide to base every in silico result (here in terms of mean values and standard deviations) on 100 simulation runs.

Robustness Analysis
We use Robustness Analysis to investigate how sensitive the output is to local parameter perturbations, that is to say when input parameters are varied one at a time. Figures 8,9,10,11,12,13,14 provide boxplots andÂ-measures that demonstrate the effect that local perturbations of the input variables µ, σ, Π D−s , Θ D−S , EC 50 , γ and T L→D respectively have on the output variables X 1 and X 2 . Key findings are listed below, discussing the impact of one input parameter at a time.
Remarks regarding input parameter µ: Figure 8 shows that, for small parameter perturbations, increasing the average doubling times of cells, µ, overall decreases the percentage of γH2AX positive cells and increases the cell count, however this decrease/increase is not linear. This indicates that the results of the in vitro simulation (and of the in vitro experiment nonetheless) are sensitive to the timing of the drug administration. In other words, Robustness Analyses demonstrates that treatment responses depend on how many cells are in the susceptible cell-cycle state at time of drug administration.
Remarks regarding input parameter σ: Figure 9 demonstrates that the level of cell cycle synchronisation amongst cells, quantified by the input σ, affects in silico outputs for small parameter perturbations. The results indicate that for highly asynchronised cells (i.e. high σ-values) the smoother growth curves yield higher cell counts at certain time-points (such as the end time 72 hours) and a lower percentage of γH2AX-positive cells. As discussed in the remark above, the timing between cell cycles and drug administration affect treatment responses.
Remarks regarding input parameter Π D−S : Figure 10 illustrates that increasing the probability that a cell enters the damaged S state, i.e. the variable Π D−S , increases the percentage of γH2AX cells and decreases the cell count, as expected.
Remarks regarding input parameter Θ D−S : Figure 11 shows how the amount of time that damaged cells spend in the D-S state before attempting to repair, and thus the Θ D−S -value, affects the output. Results show that the percentage of γH2AX positive cells increases with increasing values of Θ D−S , as more damaged cells will accumulate in the D-S state. However, this does not affect the probability of cells repairing, so the cell count is not as sensitive to small perturbations of Θ D−S . The value of ΘD − S implicitly affects the measured cell count at the end time of the experiment as a decreased/increased ΘD − Svalue yields a slightly decreased/increased time lag between a cell entering the D-S state and dying.
Remarks regarding input parameter EC 50 : Figure 12 demonstrates that output variables are highly sensitive to perturbations of EC 50 . Increasing EC 50 results in a higher percentage of γH2AX positive cells and a lower cell count. Thus the input parameter EC 50 should be regarded as a highly influential on quantitative results.
Remarks regarding input parameter γ: Figure 13 illustrates that output variables measured at the end time of the experiment are not very sensitive to small perturbations of γ. This can be understood as the γ parameter inherently corresponds to 'how quickly' a drug achieves asymptotic behaviour in the E max model, the model used in our mathematical framework to express cellular drug response.
Remarks regarding input parameter T L→ : Figure 14 shows how output variables change as a result of perturbations to the input variable T L→ , that describes how long it takes for a cell that has failed to repair to die (i.e. how long a cell with a 'death-sentence' is picked up as γH2AX positive in the in vitro experiment). Results show that both the percentage of γH2AX positive cells and the cell count increases with increasing values of T L→ , as dying cells will categorised as γH2AX positive longer before being categorised as dead. When calibrating the model, we avoid the effect of this input parameter by only regarding in vitro data at time points that are early enough to correspond to systems with no (or a negligible amount of) dead cells.

Latin Hypercube Analysis
Latin Hypercube Analysis is here used to investigate how sensitive output responses are to global parameter perturbations. We here investigate parameter values within ranges that we consider to be 'plausible' from the calibration process and the Robustness Analysis. Figures 15,16,17,18,19,20,21 provide scatter-plots that demonstrate correlations between the output variables X 1 and X 2 and the input variables µ, σ, Π D−s , Θ D−S , EC 50 , γ and T L→D respectively. The Pearson Product Moment Correlation Coefficients between the various input-output pairs are listed in Remarks regarding input parameter µ: Figure 15 and the first column in Table 6.3 show that, for the allowed parameter range, µ and X 1 are moderately, negatively correlated as the correlation coefficient is -0.48 and the scatterplot displays an overall trend of the output (X 1 ) decreasing with increasing values of the input µ. The relationship between µ and the other output variable X 2 is, on the other hand, negligible. We explain this by the fact that treatment responses are sensitive to the timing of the drug administration, but there is a time-lag T L→D between a cell's lethal event (failure to repair) and its death. As damaged (but not dead) cells are included in the cell count, the (µ,X 2 )-relationship is more weakly linearly correlated than the (µ,X 1 )-relationship.
Remarks regarding input parameter σ: Figure 16 and the second column in Table 6.3 demonstrate that the linear relationships between input variable σ and the output variables X 1 and X 2 are both negligible, within the regarded input parameter value range.
Remarks regarding input parameter Π D−S : Figure 17 and the third column in Table 6.3 indicate that the relationships between the input variable Π D−S and the output variables X 1 and X 2 are, respectively, positively and negatively weakly linearly correlated. This agrees with the intuitive notion that if the probability that a cell enters the D-S state increases, cell damage (X 1 ) increases whilst the cell count (X 2 ) decreases.
Remarks regarding input parameter Θ D−S : Figure 18 and the fourth column in Table 6.3 show that the input variable Θ D−S is has a negligible linear correlation with the output variables X 1 and X 2 .
Remarks regarding input parameter EC 50 : Figure 19 and the fifth column in Table 6.3 demonstrate that the input variable EC 50 impacts the output responses more than do other input variables, within the regarded ranges for input variables. EC 50 is negatively, moderately linearly correlated with X 1 and EC 50 is strongly, positively linearly correlated with X 2 . These relationships are visually apparent in the regarded scatterplots.
Remarks regarding input parameter γ: Figure 20 and the sixth column in Table 6.3 indicate negligible linear correlations between the input parameter γ and both output variables X 1 and X 2 .
Remarks regarding input parameter T L→D : Figure 21 and the last column in Table 6.3 demonstrate that the input variable T L→D is positively, weakly, linearly correlated with the output X 1 , whilst the linear correlation between T L→D and X 2 is negligible.