Targeting Cellular DNA Damage Responses: Predicting in vivo treatment responses using an in vitro-calibrated agent-based mathematical model

The ATR (ataxia telangiectasia mutated and rad3-related kinase) inhibitor AZD6738 is an anti-cancer drug that potentially hinders tumour proliferation by targeting cellular DNA damage responses. In this study, we combine a systems pharmacology approach with an agent-based modelling approach to simulate AZD6738 treatment responses in silico. The mathematical model is governed by a set of empirically observable rules. By adjusting only the rules, whilst keeping the fundamental mathematical framework and model parameters intact, the mathematical model can first be calibrated by in vitro data and thereafter be used to successfully predict treatment responses in human tumour xenografts in vivo qualitatively, and quantitatively up to approximately 10 days post tumour injection.


Introduction
Two well-studied ATR inhibitors are namely 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 [2,6]. Combination treatments with AZD6738 and radiotherapy or chemotherapy have produced synergistic results in preclinical settings [2], and AZD6738 is currently being evaluated in clinical phase I/II trials [1]. The intravenous ATR inhibitor VX-970 has demonstrated tumour controlling effects in a phase I clinical trial, both as a monotherapy and in combination with the chemotherapy drug carboplatin [1]. One cancer type that commonly features disruptions in the DDR is triple negative breast cancer (TNBC) [5]. Tu et al. [5] demonstrated that VX-970 can be used as a radiosensetiser and treatment intensifier to treat TNBC in patient-derived xenografts (PDX).
In a previous mathematical study by Checkley et al. [7], tumour responses to the ATR inhibitor AZD6738 are simulated using coupled ordinary differential equations, where both AZD6738 monotherapies and combination treatments with IR are regarded. In their study, a pharmacokinetic and pharmacodynamic (PK/PD) model of tumour growth is integrated with a mechanistic cell cycle model. Their model, which is calibrated by in vitro experiments, is predictive of in vivo xenograft studies. The model is aiding in quantitatively predicting dose and scheduling responses in a clinical Phase I trial design. The in vitro and in vivo data used in our current study is gathered from this previous work by Checkley et al. [7]. The mathematical framework used in our study is an extension of a mathematical model introduced by Powathil et al. [8] that has previously been used to study tumour growth, chemotherapy responses, radiotherapy responses, drug resistance and more [8,9,10,11,12].
Mathematical models, and their corresponding in silico experiments, can be used to simulate both in vitro and in vivo tumour scenarios. However, the microenvironment in an in vitro cell culture is significantly different from the microenvironment in a solid tumour, and many details that influence tumour dynamics differ between in vitro and in vivo settings. These details include cell proliferation, oxygen distribution and drug delivery. It follows that translating quantitative in vitro findings to in vivo predictions is not straightforward. In this study we use a multiscale, hybrid, agent-based mathematical model, in which one agent corresponds to one cancer cell (in vitro) or one group of cancer cells (in vivo). The model is governed by a few observable and well documented principles, or rules. To account for differences between the in vitro and in vivo scenarios, these rules are appropriately adjusted when moving from in vitro simulations to in vivo simulations. By only adjusting the rules, whilst keeping the fundamental mathematical framework intact, the model can first be calibrated by in vitro data and thereafter be used to successfully predict treatment responses in human tumour xenografts in vivo. Since in vitro data generally is easier to produce than in vivo data, this paper describes a modelling approach that facilitates moving towards in silico informed drug development and, ultimately, clinical trials. 3 We take a minimal parameter, agent-based modelling (ABM) approach, and use a cellular automaton (CA) to model a population of cancer cells (in vitro), or a tumour (in vivo), that evolves in time and space. In order to account for differences between in vitro and in vivo scenarios, model rules are appropriately adjusted, as described throughout Section 2. Differences between the in vitro and the in vivo model are summarised in Figure 2 in Section 2.8. In the model, one agent corresponds to one cancer cell (in vitro) or one group of cancer cells (in vivo) that is distinct from other agents in the system. The behaviour and fate of each agent is governed by both intracellular and environmental dynamic variables, that are integrated using multiscale and hybrid modelling techniques. At the start of an in silico experiment, one initial agent is placed in the centre of the lattice. This initial agent divides and duplicates to gives rise to a heterogeneous population of agents. When the population has reached an appropriate size (determined by the experiments performed by Checkley at al. [7]), AZD6738 anticancer treatments commence.

Cellular Automaton Lattice
Every lattice point in the CA is either empty or occupied by one agent. If a lattice point is empty, it consists of extracellular solution (in vitro) or extracellular matrix (ECM) in vivo). The ECM comprises multiple components such as collagen, elastin and fibronectin but we do not distinguish between these components in the mathematical model [13]. Differences between the simulated in vitro an in vivo lattices are described bellow.
In vitro lattice: In the in vitro experiments performed by Checkley et al. [7], populations of LoVo (human colon carcinoma) cells were plated, and population sizes of up to roughly 4000 cells were investigated. In our corresponding mathematical in vitro model, we regard cells on a two-dimensional square lattice with L × L = 100 × 100 lattice points, where the spacing in both spatial directions, x 1 and x 2 , corresponds to one cell diameter so that ∆x 1 = ∆x 2 = 20µm [8]. Thus the lattice simulates a physical space of 2 × 2 mm 2 [8].
In vivo lattice: In the in vivo experiments performed by Checkley et al. [7], LoVo cells were subcutaneously injected on flanks on Female Swiss nude mice in order to produce human tumour xenografts. Treatments started when the tumours had reached a volume of 0.2-0.3 cm 3 [7]. Approximating the tumour as spherical, we simulate (only) a central cross section of the tumour as an, approximately circular, plane 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, where the tilde-notations over the variables denote the in vivo setting. The dimensions are chisen in order to allow our agent-based model to simulate these physical dimensions, whilst keeping computational costs low, we let one agent correspond to one group of cancer cells. Post simulation time, the two-dimensional cross section of cells is extrapolated to represent a three-dimensional tumour-spheroid. This extrapolation process is outlined in the Supplementary Material. Paths illustrate transitions between states, and symbols next to the paths denote the probability that the respective path will be taken. The teal, dashed path indicates a replication path that can be inhibited by ATR-inhibitors, such as AZD6738, as is further described in Section 2.6.

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 Checkley et al. [7]. In this model, a cell progresses through various states in the cell cycle, where the states correspond to different cell cycle phases. As is illustrated in Figure 1, 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. The cause of cell death is here unrepaired replication stress. A cell can take different possible paths through the cell cycle, and every time the cell cycle path forks, stochastic 'dice-rolls' determine which path will be taken. Every cell commences its life in the G1 state of the cell cycle, 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 fitted from in vitro data reported by Checkley et al. [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 AZD6738. The higher the drug-concentration is, 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 can not repair, it is sentenced to die. Whether a cell in state D-S repairs or dies is decided by a stochastic 'dice-roll' influenced by the local drug concentration C(x, t), as is further described 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 cycle.
Each agent i on the lattice is assigned an individual doubling time τ i , where τ i is stochastically selected from a normal distribution with mean value µ and standard deviation σ. A normal distribution is chosen in order to model the in vitro cells as being almost in sync, this is motivated by observations of the in vitro data, listed in the Supplementary Material, which implies fairly synchronised cell cycles amongst cells. Each agent is also attributed an individual cell cycle clock, that determines when the agent should progress to a subsequent state in the cell cycle. Every agent that is placed on the lattice commences its life in the G1 state, and progression to a subsequent state (i.e. S or D-S) is scheduled to occur once an agent has spent a certain fraction of its doubling time in the G1 state. The fraction of the doubling time spent in the G1, S (including D-S), 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 [14]. 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]. In order to account for differences between in vitro and in vivo systems, cell cycle model rules are adjusted as described below.
In vitro cell cycle model rules: Each cell has an individual doubling time, and individually (stochastically) determines which paths to take between the various states in the cell cycle model.
In vivo cell cycle model rules: Each agent comprises a group of identical cells that simultaneously and uniformly progress through the states in the cell cycle model.

Cell Proliferation
When an agent has completed the mitoses 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 the mother agent, where up to ν order of neighbourhoods are regarded. After all lattice points in neighbourhood n are occupied, daughter cells may be placed in neighbourhood n + 1. To accomplish circular-like growth, the model stochastically alternates between placing daughter agents on Moore and von Neumann neighbourhoods of the mother agent. Modelling rules for the in vitro and in vivo scenarios are outlined below.
In vitro proliferation rules: In the experimental in vitro setup, there is no spatial constraint or lack of nutrient deficiency 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 ν in vitro = L/2, where L corresponds to the size of the lattice, in order to restrict daughter cells to the lattice.
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 neighbourhood of the 'mother agent', so thatν = 3, in accordance with previous mathematical models [8]. For the in vivo experiment regarded 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 [15]. Thus in the model, when an agent is in the G1 phase, it continues to progress through the cell cycle model, provided that free space is available on the lattice within in 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, 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. [16,17,18]. Hypoxic (cancer) cells proliferate slower than do well-oxygenated cells [15] and, hypoxic tumour regions express reduced sensitivity to radiotherapy and a plethora of chemotherapeutic drugs [19]. In our model, the hypoxic cells (i.e. cells with a partial pressure of oxygen (pO 2 ) value of 10 mmHg or less [20]) display arrest (i.e. delay) in the G1 phase of the cell cycle [21]. As is described bellow, details regarding oxygen dynamics are included in the in vivo model, but not in the in vitro model.
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 low oxygen levels are incorporated in the mathematical in vitro model.
In vivo oxygen distribution and responses: Within solid tumours, oxygen concentrations typically vary and hypoxic regions are common tumour features [21,18,22]. Avoiding complicated details of vasculature in the model, we here approximate the oxygen as diffusing in from 'outside 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 D K (x, t) denotes the oxygen diffusion coefficient, and r K and φ K are source 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 over cells than in the ECM, 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 at time t and 0 otherwise. Similarly, the binary factor cell(x, t) is 1 if there is a viable cell in locationx at time t, and 0 otherwise [8]. Equation 1 is coupled with no-flux boundary conditions so that the total amount of oxygen in the system fluctuates over time [23]. A scaled oxygen variableK(x, t) is introduced in order to express oxygenation levels in units of mmHg.K(x, t) is obtained bŷ where maxx ,t K(x, t) denotes the maxima occurring K(x, t)-value at time t [23] and h is a scaling factor. The scaling factor is included in order to achieve a lattice on which lattice points have oxygenation levels ranging between 0 and 100 mmHg.
Low cellular oxygen levels have been shown to delay cell cycle progression by inducing arrest in, particularly, the G1 phase of the cell cycle [15]. In mechanistic Tyson-Novak type cell cycle models [24,25,26], the cell cycle is governed by a system of ordinary differential equations (ODEs) in which the G1 phase is inherently elongated under hypoxic conditions by incorporating hypoxiainduced factors into the ODEs [8]. In this model, we use a clock to model cell cycle progression and thus we introduce a G1 delay factor (G1DF) in order to achieve a longer G1-phase under hypoxic conditions where The G1DF is an approximation for how much the G1 phase is expanded in time as a function of oxygen pressure. It is matched to fit data points extracted from a previous mathematical study by Alarcon et al. [15], 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 [21]. 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, as described below. 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 thus approximate the drug distribution across the lattice to be instantaneous (occurring at treatment time T 0 ) and homegeneous. 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 (unless drug wash-out is applied). In our mathematical model, this is equivalent to there being no drug decay. Hence in the mathematical in vitro model, 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 smooth drug gradient within the tumour. In the mathematical model, this drug dynamics is modelled using a mechanistic partial differential equation (PDE), where the concentration of the drug AZD6738 in lattice point n at locationx at time t is denoted by C(x, t) and where D AZD is the diffusion coefficient of the drug AZD6738, and the production coefficient p(x, t) is greater than zero at drug administration times only for lattice points outside the tumour, i.e. in the extracellular matrix. The value of p(x, t) corresponds to the amount of administered drug. 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 [27]. 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 predominately 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 [28]. Using this assumption, the drug diffusion coefficient is set in relation to the oxygen diffusion coefficient, as is done in previous mathematical studies [8]. Thus the relationship between the diffusion coefficients for the drug (AZD) and oxygen (O 2 ) corresponds to the square of the inverse relationship between the corresponding molecular weights, such that where the molecular weights are collected from the PubChem database [29]. Details regarding pharmacokinetics are outside the scope of this study, bioavailability is instead calibrated using the extreme scenario in which the maximum drug dose is administered in vivo.

Drug Responses
AZD6738 inhibits the repair from the D-S state to the S state in the cell cycle, as illustrated in Figure 1, and maximal drug effect corresponds to complete repair inhibition. The drug effect is modelled using an agent-based adaptation of the sigmoid Emax model [30], 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 concentration required to achieve half of the maximal drug effect,  0.5 · E max and γ is the Hill-exponent [30]. EC 50 and γ are fitted from the in vitro data. When an agent is scheduled to progress from the D-S state in the cell cycle, it has a probability Π rep to repair which is determined by the local drug concentration such 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' [28]. In the in vivo setting, the cellular debris is digested by macrophages but in the in vitro setting such 'cell-corpses' linger on the lattice. 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: The survival probability of a cell is modelled using the Emax model. After failure to repair from the D-S state, a cell 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's survival probability is modelled using the Emax model. An agent 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
sec:ddrparams The parameters used in the mathematical framework are calibrated by in vitro data and are listed in Table 1. The in vitro and in vivo data produced by Checkley et al. [7] is listed in the Supplementary Material. The Supplementary Material also includes information regarding the model calibration process. In the context of quantitative pharmacology, knowledge about a model's robustness is crucial [31]. Accordingly, results from the uncertainty and sensitivity analysis are also included in the Supplementary Material. Three different uncertainty and sensitivity analyses techniques, suitable for agent-based models with stochastic elements, are used as described in the Supplementary Material. These three techniques are namely (i) Consistency Analysis, (ii) Robustness Analysis and (iii) Latin Hypercube Analysis [32,33].

Differences Between in vitro and in vivo Modelling Rules
The differences between the in vitro and in vivo rules used in the mathematical framework are pictorially summarised in Figure 2.

Implementation
The mathematical model is implemented in an in-house C/C++ framework in order to produce in silico results. Partial differential equations are solved using explicit finite difference methods.
Simulations are visualised using ParaView [34]. Results reported in Section 3 are based on 100 simulation runs as, according to results from the uncertainty analysis, 100 runs are sufficient for mitigating uncertainty originating from intrinsic model stochasticity. Uncertainty and sensitivity analyses are performed using MATLAB [35].

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. By only adjusting the modelling rules in the mathematical framework, whilst keeping the in vitro calibrated model parameters intact, the mathematical framework is able successfully predict treatment responses in human tumour xenografts in vivo.

Simulating in vitro Experiments
The mathematical framework used in this study is calibrated by results from an in vitro experiment performed by Checkley et al. [7] in which populations of LoVo (human coloncarcinoma) cells are exposed to the ATR-inhibiting drug AZD6738. The in silico results in Figure 3 simulate the evolution of the in vitro cell population over time in terms of cell damage (Figure 3 Left) and in terms of cell count (Figure 3 Right). AZD6738 drugs are given at 0 hours, when the cell population has reached a size of approximately 1000 cells. Response curves for six different drug concentrations, including the zero-drug concentration control case, are shown. Each response curve is based on a mean values from 100 simulation runs, according to the uncertainty analysis described in the supplementary material. Also shown in Figure 3 are simulation standard deviations and in vitro results produced by Checkley et al. [7]. Using a minimal-parameter modelling approach, the mathematical framework is calibrated to best fit all in vitro data points without introducing any auxiliary scaling variables or similar. The calibration process is described in the Supplementary Material. Our results demonstrate that, after calibration, our mathematical framework is able to capture in vitro LoVo cell population growth and drug (AZD6738) responses.

Simulating in vivo Experiments
Post the in vitro calibration, our mathematical framework is used to simulate 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 4 simulate AZD6738 drug responses in terms of the percentage of damaged cells (Figure 4 Left) and tumour volume ( Figure  4 Right). Response curves to three different drug doses (0, 25 and 50 mg/kg) are shown, where the curves represent mean values for 100 simulation runs. Simulation standard deviations and in vivo data are also provided in Figure 4. These result graphs demonstrate that our mathematical framework quantitatively agrees with the in vivo results reported by Checkley et al. [7] for up to approximately 10 days post tumour injection. This empirically demonstrates the predictive ability of our mathematical framework and modelling approach.  Simulation mean values for 100 simulation runs (coloured lines), simulation standard deviations (black error bars) and in vivo data with standard errors (coloured error bars) [7] are shown. 13 We demonstrated that the mathematical framework used in this study can successfully simulate growth and drug (AZD6738) responses both in cancer cell populations (in vitro) and in tumour xenografts (in vivo). This versatile mathematical framework is an extension of a multiscale, hybrid, agent-based, on-lattice model that has previously been used to study tumour growth and treatment responses to chemotherapy, radiotherapy, hypoxia-activated prodrugs and more [8,9,10,11,12]. By only adjusting the rules in the mathematical framework, whilst keeping model parameters intact, the mathematical model can first be calibrated by in vitro data and thereafter be used to successfully predict treatment responses in human tumour xenografts in vivo. After comparing in vivo-scenario in silico results to in vivo data, one could chose to further fine-tune the modelling parameters in order to obtain an even better fit.
Data-driven modeling, exploitation of existing data and proof-of-concept studies are important steps involved in current and future procedure for enabling multiscale modeling in systems medicine, as argued in a report by Wolkenhauer et al. [36]. In the data-driven, minimalparameter, 'proof-of-concept' study discussed in this article, we demonstrate that our mathematical framework is able to produce in silico results that match both in vitro and in vivo data produced in a previous study by Checkley et al. [7]. In their study, Checkley at al. [7] used a pharmacokinetic/pharmacodynamic (PK/PD) model to predict in vitro and in vivo AZD6738 treatment responses. The mathematical results produced by Checkley et al.'s [7] PK/PD model agree with the mathematical result produced by our multiscale, hybrid, agent-based, on-lattice model agrees. Thus our mathematical framework is validated by, not only in vitro and in vivo data, but also by another mathematical model that is based on modelling techniques different from ours. Despite the fact that mathematical modelling is, in general, is becoming increasing popular in the pharmaceutical industry, there are not that many agent-based models present on the pharmaceutical scene [37]. We argue this is a missed opportunity in the context of oncology, as agent-based models naturally capture the heterogeneous nature of tumours, that is known to complicate treatment matters. Moreover, multiscale models enable the integration of knowledge and data on different scales in time and space. By formulating agent-based modelling rules using 'low' principles (our ideological equivalent to 'first principles' in this context) regarding how the agents in the system behave, the produced in silico data can be handled as classical in vitro and in vivo data post simulation [37].
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. Data from in vivo experiments is often sparse, as gathering in vivo data is associated with practical, financial and ethical constraints. Plentiful and adaptable in silico data are however easy to produce, and thus sparse in vivo data can be complemented and validated by in silico data. Thus mathematical frameworks, such as the one used in this study, can be used as an epistemic contribution to sparse experimental data. Since our in vitro-calibrated mathematical framework is able to predict in vivo responses, it can be used as a tool in the process of moving a drug-response investigation from an in vitro to an in vivo setting. Furthermore, the corresponding in silico experiments can be extended to investigate various dose-schedule scenarios in order to, not only validate, but also, guide in vitro and in vivo experiments. Ultimately, following interdisciplinary collaborations between clinicians, biologists and mathematicians, mathematical oncology may be used to further personalised cancer care in clinical settings according to bench-to-bedside and blackboard-to-bedside approaches [7,11,36,31]. As multiscale, agent-based models 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 the knowledge and data of its team members. 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.

Calibrating the model using in vitro data
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 described later on in 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 choice of using a normal distribution from which to pick τ i is motivated by the observation that the cells in the in vitro experiment are 'fairly synchronised', as can be seen by observing the cell count data listed in Table 1. The control case (i.e. no drug) cell count data is used to estimate µ and σ.
By observing the control case cell count data, we conclude that the average doubling time for cells should be between 22 and 26 hours, hence σ-values in the parameter range [22,26] hours are investigated in silico. Note that in the last 24-hour interval, between 48 and 72 hours, the control population increases by less than 5%. However, according to cell count data from earlier time poins, we assume that cells are on the brink of cell division at the end time of the experiment and thus we disregard the influence of the 72-hour data point in the overall model calibration. 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.

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 m ax 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, and further that the drug concentrations 0.3 and 3 µM respectively achieve less than, and more than, the half of the total drug effect. With this in mind, EC 50 values over 0.3 and under3 µ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 5 0, γ) 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 siico 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.

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,X 2 = Total Tumour Volume = 4π 3 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,

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 1, 2, 3, 4, 5 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 1 through to 5, it is clear that the statistical significance decreases with increasing distribution size n, as is shown in Figure 6 and

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 7,8,9,10,11,12,13 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 7 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 8 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 9 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 10 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 11 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 12 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 13 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 14, 15 Remarks regarding input parameter µ: Figure 14 and the first column in Table 4.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 15 and the second column in Table 4.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 16 and the third column in Table 4.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 17 and the fourth column in Table 4.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 18 and the fifth column in Table 4.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 19 and the sixth column in Table 4.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 20 and the last column in Table 4.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.