Mathematical modelling reveals cellular dynamics within tumour spheroids

Tumour spheroids are widely used as an in vitro assay for characterising the dynamics and response to treatment of different cancer cell lines. Their popularity is largely due to the reproducible manner in which spheroids grow: the diffusion of nutrients and oxygen from the surrounding culture medium, and their consumption by tumour cells, causes proliferation to be localised at the spheroid boundary. As the spheroid grows, cells at the spheroid centre may become hypoxic and die, forming a necrotic core. The pressure created by the localisation of tumour cell proliferation and death generates an cellular flow of tumour cells from the spheroid rim towards its core. Experiments by Dorie et al. showed that this flow causes inert microspheres to infiltrate into tumour spheroids via advection from the spheroid surface, by adding microbeads to the surface of tumour spheroids and observing the distribution over time. We use an off-lattice hybrid agent-based model to re-assess these experiments and establish the extent to which the spatio-temporal data generated by microspheres can be used to infer kinetic parameters associated with the tumour spheroids that they infiltrate. Variation in these parameters, such as the rate of tumour cell proliferation or sensitivity to hypoxia, can produce spheroids with similar bulk growth dynamics but differing internal compositions (the proportion of the tumour which is proliferating, hypoxic/quiescent and necrotic/nutrient-deficient). We use this model to show that the types of experiment conducted by Dorie et al. could be used to infer spheroid composition and parameters associated with tumour cell lines such as their sensitivity to hypoxia or average rate of proliferation, and note that these observations cannot be conducted within previous continuum models of microbead infiltration into tumour spheroids as they rely on resolving the trajectories of individual microbeads.


Introduction
By the time tumours are clinically detectable in vivo they are typically highly heterogeneous in terms of their spatial composition [1]. Tumours contain multiple cell types, including stromal cells (e.g., fibroblasts) and immune cells (e.g., macrophages, T cells) and their growth is sustained by an irregular network of tortuous and immature blood vessels which deliver vital nutrients such as oxygen to the tumour cells. When characterising tumour cell lines or testing new cancer treatments it is important to have a reproducible experimental assay. In such situations, tumour spheroids are widely used due to the predictable manner in which they grow [2].
Tumour spheroids are clusters of tumour cells whose growth in vitro is limited by the diffusion of oxygen and other nutrients, such as glucose, from the surrounding medium into the spheroid centre. Other factors which may limit the growth of tumour spheroids include intercellular communication, contact sensing, pH levels and/or the circadian clock. In small spheroids, all cells receive sufficient nutrients to proliferate and exponential growth ensues. As a spheroid increases in size, nutrient levels at its centre decrease and may eventually become too low to support cell proliferation, driving cells to halt division and become quiescent. Slower growth of the spheroid will occur until nutrient levels at its centre fall below those needed to maintain cell viability, leading to the formation of a central necrotic core containing dead cells. Growth will continue until the spheroid reaches an equilibrium size at which the proliferation rate of nutrient-rich cells in the outer shell of the spheroid balances the degradation rate of necrotic material at the spheroid centre [2][3][4]. During necrosis, the cell membrane collapses causing rapid ejection of cell constituents into extracellular space [5], leading to a reduction in cell size as liquid matter disperses into the spheroid.
A wide range of models have been developed to describe the growth and mechanical properties of tumour spheroids [6][7][8] and organoids [9,10] and their response to treatment [11,12]. The simplest models, which include logistic growth and Gompertzian growth, recapitulate the characteristic sigmoid curve describing how the total spheroid volume changes over time [13][14][15]. These phenomenological models are, however, unable to describe the internal spatial structure of tumour spheroids. More detailed mechanistic models relate the internal spatial structure of the spheroids to the supply of vital nutrients such as oxygen and glucose [16][17][18][19][20], and may be adapted to include the effect of anti-cancer treatments. While some models of spheroid growth account explicitly for factors such as glucose, ATP, pH, and contact inhibition of cell proliferation (e.g., [21]), it is common in mathematical models of tumour spheroids to simplify these complex metabolic processes while retaining the qualitative behaviour of the experimental observations. Most models therefore represent oxygen, glucose and other nutrients via a single diffusible species described variously as "oxygen" or "nutrient", which is assumed to be vital for the survival and proliferation of tumour cells (e.g., [22][23][24]).
Agent-based models (ABMs), which resolve individual cells, can also be used to model tumour spheroids. ABMs are often multiscale, linking processes that act at the tissue, cell and subcellular scales. For example, the cell cycle dynamics of individual cells may be modelled via ordinary differential equations (ODEs) at the subcellular scale, may depend on local levels of tissue scale quantities such as oxygen concentration, and may influence cell scale processes such as cell proliferation. ABMs are termed 'hybrid' if they combine different modelling approaches. For example, a reaction-diffusion equation describing the spatial distribution of oxygen within a tumour spheroid may be coupled to a stochastic, rule-based cellular automata (CA) model governing the dynamics of individual tumour cells [25]. ABMs can be formulated using on-and off-lattice approaches. On-lattice approaches include rule-based CA models (e.g., [26,27]) in which each lattice site is typically occupied by at most one cell, and the cellular Potts model [28][29][30] where individual cells may occupy multiple lattice sites. CA models are generally unable to capture realistic cell shapes or intercellular forces. However, while cellular Potts models permit more realistic cell geometries, they are also constrained by a lattice and do not allow full consideration of mechanical effects.
A weakness of on-lattice models is that cell locations are restricted to discrete lattice sites. By contrast, off-lattice ABMs allow cells to move in a continuous manner through space. Examples of off-lattice ABMs include those which track cell centres (cell-centre approaches) and those which track the cell boundaries (vertex-based approaches). We refer the interested reader to [31] for a comparison of five ABMs (CA, cellular Potts models, overlapping spheres [32], Voronoi tesselation [33] and vertex-based methods).
While many groups develop their own ABMs (e.g., [34]), increasing numbers are using open source software specifically designed for simulating ABMs. The Chaste framework [35,36] is designed to implement a wide range of ABMs. PhysiCell [37] utilises BioFVM [38] to obtain efficient simulations involving large numbers of diffusing substrates such as oxygen. Morpheus [39] focusses on user-friendliness, with a GUI designed to bypass many of the coding challenges associated with developing agent-based models. Like Chaste, it can implement models using a range of on-lattice or off-lattice frameworks. CompuCell3D [40] provides an intuitive way to implement models using the cellular Potts framework [28]. Other software tools that implement agent-based models include CellSys [41], Biocellion [42], HAL [43] and Timothy [44]. Advantages of these software tools are that code can be more effectively reused and benchmarked, and errors in model implementation are more likely to be identified. Taken together, these frameworks also provide multiple ways of implementing ABMs, so that researchers can choose the framework (or frameworks) best suited to the questions they seek to address.
When developing theoretical models of tumour spheroid growth, a key consideration is the experimental data available to validate and/or parameterise the model. Typically, spheroid experiments generate dynamic data showing how the total tumour volume changes over time. These data may be supplemented by spatially-resolved images of spheroid composition at discrete timepoints [3,45]. In a series of papers, Dorie et al. adopted an alternative approach [46,47]. They added inert microbeads to the outer edge of well-developed tumour spheroids and collected time-series data showing how the spatial distributions of the microbeads changed over time, as they moved radially inwards, towards the centres of the spheroids. This data is consistent with "passive infiltration", in which inert particles are advected into the spheroid by tumour cells which are moving down pressure gradients caused by spatial variation in cell proliferation and death. This cellular flow, induced by mechanical stresses within a tumour spheroid, has been proposed as a mechanism by which drugs may exploit advection to enter a tumour more efficiently [48]. Passive infiltration of this type differs from advection by fluid flow, which has previously been implicated as a mechanism by which chemotherapeutic drugs may enter, or be inhibited from entering, tumours [49,50].
Several authors have developed continuum models to describe Dorie et al.'s experiments [51][52][53]. These continuum models focus on cell populations and, as such, do not resolve individual cells. McElwain and Pettet [51] showed that a possible cause of bead internalisation is pressure gradients caused by differential rates of cell proliferation and death between the proliferative rim and the necrotic core. A modified version of this model [53] also distinguishes between proliferating and quiescent cell populations. Both models assume that dead cells are instantaneously removed from the spheroid and do not occupy space. They also require chemotactic movement of tumour cells in response to the oxygen gradient to reproduce the results from [46]. Thompson and Byrne [52] developed an alternative continuum model to explain the observed infiltration patterns, arguing that differences in infiltration can be explained by non-uniform death and proliferation in the spheroids. Their model assumes for simplicity that the tumour spheroids do not contain a necrotic core. While it reproduces the observed data, their model is unable to reproduce predicted infiltration patterns unless the microbeads are assumed initially to lie strictly inside the tumour spheroid.
In order to reproduce additional data from Dorie et al. [46] describing the infiltration of tumour cells labelled with tritiated thymidine ( 3 H-labelled cells), McElwain and Pettet [51] assumed that tumour cells move, via chemotaxis, up spatial gradients in the oxygen concentration, and neglected proliferation of the labelled cells. Thompson and Byrne [52] assumed that labelled cells proliferate, but as with their microbead model the labelled cells had to be initially placed within the tumour spheroid to reproduce infiltration results. By varying chemotaxis coefficients between the subpopulations, the authors reproduced the observed infiltration patterns. As they use a continuum framework, these models cannot resolve individual microbeads or labelled cells.
Inspired by the original experiments of Dorie et al. [46], and the increasing use of ABMs, we develop an off-lattice hybrid ABM to describe the growth of tumour spheroids and their infiltration by microbeads and 3 H-labelled tumour cells. This enables us to resolve the trajectories of individual cells and microbeads within the model, something which continuum models are unable to do. We show further how observations of microbead trajectories can be used to identify the composition of a simulated spheroid and estimate parameters associated with cell cycle duration and cell responses to hypoxia. After verifying that the model reproduces the spatio-temporal dynamics of tumour spheroids growing in free suspension in vitro, we use it to simulate their infiltration by inert microbeads and labelled tumour cells (see Supporting Information, S5 Appendix: Replication of 3 H-labelled cells experiments), obtaining good agreement with Dorie et al.'s experimental results. A parameter sensitivity analysis reveals how the growth rate, size and spatial structure of the spheroids change as we vary key model parameters. We show how spheroids with the same equilibrium size may differ in their spatial organisation. We conclude by showing how dynamic data describing the trajectories of individual microbeads, which cannot be resolved using continuum models, can be used to infer the composition of simulated tumour spheroids, and also to estimate model parameters pertaining to tumour cell proliferation rates and cell sensitivity to hypoxia.

Model overview
We develop a hybrid agent-based mathematical model to describe the in vitro growth of tumour spheroids in response to an externally-supplied nutrient, here taken to be oxygen. We use the model to simulate spheroid infiltration by inert microbeads and 3 H-labelled cells. Our model is implemented within Chaste (Cancer, Heart and Soft Tissue Environment, available at https://www.cs.ox.ac.uk/chaste/), an open source simulation package designed to solve computationally demanding, multiscale problems that arise in biology and medicine [35,36]. We choose this framework because it provides an efficient means of implementing off-lattice ABMs, and has previously been used to simulate multicellular spheroids [35,54]. The extensions to the framework described in this paper will be made available in a subsequent release of the Chaste software.
A schematic highlighting the key features of our hybrid ABM is presented in Fig 1. Individual cells are represented using an off-lattice, cell-centre model, and cell movement is determined via consideration of the force balances on each cell and assuming that inertial effects can be neglected (Panel C of Fig 1 indicates those forces which act on cells). Interactions between cells are modelled by connecting cell centres with springs (Fig 1, Panel C), which simulate both intercellular adhesion and volume exclusion. Through appropriate choice of the spring lengths associated with each cell, the model includes a notion of cell size which can be adjusted to account for size differences between developed cells and those which have just proliferated, or which are decaying due to necrosis.
We distinguish two types of agent: tumour cells and microbeads. Tumour cell behaviours (e.g., cell cycle progression, quiescence and cell death) depend on the local oxygen concentration, which is determined by a reaction-diffusion equation accounting for oxygen consumption by viable tumour cells and oxygen diffusion from the spheroid boundary towards its centre. The flowchart at the bottom of Fig 1 shows how simulation results are generated. We note that while this model can simulate spheroid growth in three dimensions, here we restrict attention to 2D simulations to reduce computational time.
In the rest of this Section, we introduce the agent-based model we have developed to simulate tumour spheroid growth and microbead infiltration. We first describe the PDE used to calculate the oxygen distribution throughout the spheroid, and then discuss the impact oxygen concentration has on tumour cell behaviour. We outline the rules used to implement proliferation and death, and the forces which act on different types of agent to determine their movement. Finally, we summarise the different simulations conducted using this model and the rules used to initialise them.

Oxygen distribution
While we model cells and microbeads as discrete entities, we assume that the concentration of oxygen ω(x, t) is continuous and can be described via a reaction-diffusion equation, with oxygen consumption by live tumour cells modelled by placing point sinks at the centres of viable cells. Written in dimensional form, the equation governing the spatio-temporal evolution of the oxygen concentration is thus for x 2 O, where x i is the location of viable cell i, the parameter D ω is the assumed constant diffusion coefficient of oxygen, κ is the oxygen consumption rate and O is the simulation domain, which we take to be a square domain large enough to fully enclose the spheroid. δ(x) is the delta function (δ(x) = 1 when x = 0; δ(x) = 0 otherwise). Eq (1) is solved subject to Dirichlet boundary conditions, which are prescribed on the domain boundary, and suitable initial conditions. We assume that oxygen is maintained at a constant level in the culture medium surrounding the tumour spheroid and, hence, by Overview of our agent-based model for the growth of tumour spheroids. A: As oxygen ω diffuses from the outer spheroid boundary, it is consumed by live cells. Consequently, the oxygen concentration at the centre decreases as the spheroid increases in size. We use the oxygen distribution to distinguish up to four different regions, or compartments, within a spheroid: a welloxygenated rim, where ω q � ω � 1, contains proliferating cells; a quiescent compartment, where ω h � ω < ω q , contains nonproliferating viable cells; a hypoxic compartment, where ω � ω h , contains non-proliferating viable cells which will become necrotic if they remain hypoxic for longer than a prescribed time period; and a necrotic compartment, containing dead cells which degrade over time. B: Schematic showing how the way in which cells switch between different compartments depends on the local oxygen concentration. C: Schematic showing the forces which act on individual cells to determine cell movement. All nodes experience continuity that the oxygen concentration on the spheroid boundary is also maintained at this constant value, ω 1 .
Since the timescale for oxygen diffusion (seconds) is much shorter than the timescale for cell proliferation (hours), when we nondimensionalise Eq (1), we make the standard, quasisteady state approximation (see S3 Appendix: Non-dimensionalisation of oxygen equation for details) and solve instead We solve Eq (2) on a regular tetrahedral finite element mesh which spans O, a square domain large enough to contain all the cell centres. We fix ω = ω 1 at any nodes of the mesh which lie outside the spheroid.

Tumour cell phenotypes
While microbeads do not consume oxygen and are unaffected by the oxygen concentration, the phenotype of tumour cells depends on their local oxygen concentration. We distinguish four types of tumour cell behaviour by introducing the following phenotypes (see Fig 1): • A tumour cell is proliferative if its local oxygen concentration exceeds a threshold value, ω q .
• If ω(x, t) � ω q , then the cell becomes quiescent and immediately pauses its cell cycle. If ω increases above ω q then the cell immediately becomes proliferative and resumes its cell cycle.
• If the oxygen concentration falls below a second, hypoxic threshold, 0 � ω h � ω q , then the cell immediately becomes hypoxic. A hypoxic cell can re-enter the quiescent compartment if its oxygen concentration rises back above ω h . However, if a cell remains hypoxic for sufficiently long (t i hours) then it will irreversibly become necrotic [55].
• A necrotic cell is dead, and no longer consumes oxygen. Necrotic cells continue to occupy space for approximately � t hours before being removed from the simulation.
We also simulate cells labelled with tritiated thymidine ( 3 H-labelled cells) to replicate experimental results from Dorie et al. [46] (see Supporting Information, S5 Appendix: Replication of 3 H-labelled cells experiments). In this model 3 H-labelled cells are assumed to behave in the same way as other tumour cells. They adopt the same phenotypes as unlabelled cells in response to local environmental cues. The only difference between 3 H-labelled cells and other tumour cells is the presence of a label which which is passed on to their daughter cells, enabling their lineage to be tracked.

Tumour cell proliferation and death
For each individual cell i, proliferation and death are determined by two subcellular dependent variables: the cell cycle time denoted by T i , which tracks a cell's progress through the cell cycle, and the hypoxia time denoted byT i which determines whether a cell is sufficiently hypoxic to become necrotic. These internal timers move at a rate which depends on the local oxygen spring forces due to interactions with their neighbours, a random force which represents local fluctuations in the cell environment, and a drag force which resists cell movement. Boundary nodes also experience a surface tension force which is directed inwards, towards the spheroid centroid, and which resists spheroid expansion. D: Flowchart summarising how the ABM is updated on each timestep-see main text and Supporting Information, S2 Appendix: Algorithm for updating the cell cycle, for details. https://doi.org/10.1371/journal.pcbi.1007961.g001 concentration ω(x, t). Pseudocode describing how the cell cycle is updated is presented in the Supporting Information, S2 Appendix: Algorithm for updating the cell cycle.
Cell proliferation. At birth, the cell cycle time of cell i is initialised such that T i = 0, and the cell is assigned a cell cycle duration τ i drawn from a uniform distribution U(0.75τ, 1.25τ) where the parameter τ defines the average cell cycle length. The distribution was chosen to be sufficiently wide to ensure that cell cycles do not become artificially synchronised over time, while ensuring that τ remains a good descriptor of the average cell cycle duration. If cell i is at position x at time t then its cell cycle evolves as follows: where ω(x, t) is the local oxygen concentration at time t and location x, H is the Heaviside step When ω � ω q , the cell cycle pauses and the cell remains dormant until either ω increases above the threshold (and progress through the cell cycle resumes) or the cell becomes necrotic (details of this process are described in the next section). When T i = τ i , cell division occurs and a daughter cell is placed half a cell diameter away from the parent cell centre in a randomly chosen direction. Both cells are assigned new cell cycle durations. Their cell cycle times are reset to 0 and evolve according to Eq (3). The resting spring length of the new cells is adjusted to account for the reduced size of the new cells (for details, see the description of the spring force laws below). Cell death. Hypoxic cells at locations where ω � ω h undergo necrotic cell death if they remain hypoxic for longer than a threshold timet i . In a manner similar to that used to model cell cycle progress,t i is drawn from a uniform distribution Uð0:75t; 1:25tÞ wheret is the average duration for which a cell is hypoxic before it becomes necrotic. As for the parameter τ, the distribution aroundt enables us to account for stochastic fluctuations in cell properties and also to suppress the emergence of artificial oscillations in the number of necrotic cells caused by multiple cells simultaneously becoming necrotic.
Each cell is assigned an internal hypoxia time,T i , which progresses when a cell is hypoxic and evolves as follows: withT i ¼ 0 at the onset of hypoxia. If the cell moves or the oxygen distribution changes so that ω(x i , t) > ω h then we setT i ¼ 0, indicating that the cell has received sufficient nutrient to prevent cell death. A cell becomes necrotic whenT i ¼t i . It is then irreversibly marked for cell death. Once a tumour cell has become necrotic, it is no longer viable and no longer progresses through the cell cycle. It continues to occupy space, but reduces in size until it is removed from the simulation over a period of � t i hours where � t i is drawn from a uniform distribution Uð0:75 � t; 1:25 � tÞ and � t is the average duration of necrosis. As with the cell cycle duration and the time threshold which triggers necrosis, the range of this distribution is estimated and reflects variation in the process of cell degradation. Details on how size reduction is implemented are included below.

Force balance
We use Newton's second law to derive the equations of motion for cells and microbeads. In the over-damped limit, we neglect inertial effects and obtain the following force balances for cell i and microbead j respectively: In Eqs (5) and (6) we assume that the drag forces on cell i and microbead j are proportional to their velocities, the constant ν denoting the drag coefficient. We denote the mechanical force by F m i , random forces by F r i , and surface tension forces by F s i . Mechanical and random forces act on cells and microbeads, whereas surface tension forces only act on cells. Functional forms for these forces are introduced below.
The timestep dt used to generate numerical solutions is taken to be 1/120 of a dimensionless time unit, which is equivalent to 30 seconds (see Supporting Information, S7 Appendix: Table of parameters).

Mechanical forces, F m i (cells and microbeads.)
The mechanical spring force acting on a cell or microbead is the net force exerted on it by its neighbours. We assume that cells/microbeads i and j only interact if the distance between their centres is less than a fixed value, R int . In more detail, and following the overlapping sphere approach outlined in [32,33,56,57], if |x i − x j | < R int then the interaction force between cells/ microbeads i and j is parallel to the vector x i − x j connecting their centres. The magnitude of the force depends on the sizes of the cells/microbeads and the distance between them. While agents in this model are represented as points, each point has an associated size which is implemented by adjusting the resting spring lengths for each cell/microbead. The resting spring length between two nodes, s i,j , is the sum of the equilibrium springs for each cell (s i,j = s i + s j ). For most cells i, s i = R Cell is a constant which is approximately equal to the radius of a cell. For newborn and necrotic cells, cell growth or shrinkage may mean that s i < R Cell . These processes are described below.
If the distance between the cell centres of cells i and j is greater than s i,j then the cells experience an attractive force representing intercellular adhesion, but if the distance is less than s i,j then the force is instead repulsive and models volume exclusion. The net force acting on a cell or microbead i at location x i due to mechanical forces is the sum of the contributions over all cells and microbeads j within radius R int : where F m i;j is the mechanical force between cells i and j. This force always points in the direction of the vector between the cells. The magnitude of F m i;j is defined as follows: where x = |x i − x j | − s i,j is the overlap between cells i and j, the parameter μ represents the spring stiffness and the parameter λ determines the strength of intercellular adhesion between neighbouring cells. A sketch showing how jF m i;j j changes as x varies is shown in the Supporting Information, S1 Appendix: Spring force magnitude. The adhesive force defined in Eq (8) grows stronger as the cell centres draw closer, since more intercellular bonds form as the surface area of contact between the cells increases [32,58].
Inert microbeads also interact with neighbouring cells and microbeads. Since microbeads do not adhere to each other we modify the mechanical force described in Eq (8) when i and j are both microbeads as follows: where μ bead is the spring stiffness associated with microbeads. Beads therefore resist being compressed, but do not adhere to other beads.
With the exception of newly divided cells and necrotic cells, we assume that s i = R Cell for cells and microbeads (as the radius of a microbead is comparable to that of a cell). For convenience, we rescale lengths according to this lengthscale, assuming that 1 cell diameter = 2R Cell = 20μm. When a cell divides, the radius of the new cells is set to s i ¼ R Cell 2 and increases linearly over the course of one hour until s i = R Cell . For necrotic cells, s i decreases linearly over the course of � t i hours until it reaches 0 and the cell is removed from the simulation. The associated spring constant of springs attached to a necrotic cell is reduced linearly at the same rate, representing a weakening intercellular force between a necrotic cell and neighbouring cells as the necrotic cell degrades. These reductions in cell size are incorporated in the pseudocode in the Supporting Information, S2 Appendix: Algorithm for updating the cell cycle.

Random forces, F r i (cells and microbeads)
We assume that cells and microbeads experience random forces due to heterogeneity in the surrounding environment and that the random force, F r i ¼ ðF r x ; F r y Þ, acting on cell/microbead i during the timestep dt is given by: In Eq (10), D is a diffusion coefficient and ξ = (ξ x , ξ y ) where ξ x and ξ y are random variables drawn from a standard normal distribution.

Surface tension forces, F s i (boundary cells only)
Cells on the spheroid boundary experience a force, F s i , of the form wherex i is a unit vector pointing in the direction of the line connecting cell i on the boundary with the spheroid centroid, and the parameter β determines the strength of the surface tension force. We define boundary cells as those which belong to the α-shape of the set of cell centres, where α = R Cell [59]. α-shapes generalise the concept of a convex hull to permit concave boundaries, and can be understood as the set of points which make contact with a ball of radius α rolled around the edge of the spheroid.

Simulations
We use our ABM to perform three types of simulations: (i) growth of a tumour spheroid, (ii) microbead infiltration into a well-developed tumour spheroid, and (iii) infiltration of 3 Hlabelled cells into a well-developed tumour spheroid (see Supporting Information, S5 Appendix:

Replication of 3 H-labelled cells experiments).
Simulations of spheroid growth are initialised by uniformly distributing 300 cells in a circle of radius 5 cell diameters. The cells are then are allowed to grow for 300 hours, which is sufficient for them to attain a steady state for most parameter regimes examined here. At the start of each simulation, all cells are randomly assigned a cell cycle time T i from a uniform distribution U(0, 0.75τ). This ensures that the cell cycles are not synchronised. When performing parameter sensitivity analyses, three model parameters were varied: the average cell cycle length, τ, the oxygen threshold for quiescence, ω q , and the oxygen threshold for hypoxia, ω h . We focus on these parameters as they have a significant effect on tumour spheroid composition, and are known to vary between different tumour cell lines.
For microbead infiltration simulations, 100 microbeads were distributed randomly around the spheroid edge after 300 hours, and the simulation allowed to evolve for a further 100 hours, reproducing the conditions under which Dorie et al. [46] tracked microbeads experimentally. For 3 H-labelled cell experiments, 50 cells on the spheroid boundary were randomly selected after 300 hours and marked with a label which was transmitted on cell division but did not affect cell behaviour.
For details regarding the parameter values used here, we refer the reader to the Supporting Information, S7 Appendix: Table of parameters. Throughout this paper, when estimates of simulated tumour cell phenotype compartment widths are stated for a particular simulation realisation, this width is calculated as the difference in radial distance from the spheroid centroid between the innermost and outermost cells with that phenotype.

Units and non-dimensionalisation
We refer to non-dimensionalised times, distances and concentrations. Lengths are nondimensionalised with the diameter of a cell, so a dimensionless distance of 1.0 corresponds to approximately 20 μm. Times are rescaled so one dimensionless time unit corresponds to 1 hour. Concentrations are rescaled with the externally-supplied concentration of oxygen, so a dimensionless concentration of 1.0 equates to approximately 3.298 × 10 −6 m 3 kg −1 in dimensional units (using Henry's law to convert from partial pressure to concentration and following [12,60]).

Results
We use our hybrid agent-based model to simulate tumour spheroids whose growth dynamics and spatial structure replicate those of spheroids cultured in vitro. We first demonstrate that this model reproduces the sigmoidal growth curves observed in in vitro experiments. we then show that when microbeads are added to the in silico spheroids their infiltration patterns are qualitatively similar to those described by Dorie et al. [46]. By performing a systematic parameter sensitivity analysis, we show how varying the oxygen thresholds ω q and ω h and the average cell cycle length τ can generate in silico tumour spheroids with similar growth dynamics and different spatial compositions. Finally, we demonstrate that resolving the trajectories of individual microbeads within our agent-based model provides sufficient information to predict spheroid composition and infer certain model parameter values.

Model Validation
We first show that our ABM qualitatively reproduces the dynamics and changing spatial structure that characterises spheroid growth in vitro. Four parameters relating to the proliferation and death rates of tumour cells were varied: the oxygen thresholds at which cells become quiescent or hypoxic, ω q and ω h , the average cell cycle duration τ and the average duration of hypoxia before cell death,t (ranges given in Supporting Information, S7 Appendix: Table of parameters). Typical simulation results generated using representative parameter values (ω q = 0.5, ω h = 0.3, τ = 16,t ¼ 8) are presented in Fig 2. We observe a short burn-in period of approximately 2-4 hours, during which the cells move from their initial positions to form a densely packed spheroid. Thereafter the tumour radius grows rapidly, slowing when nutrient levels at its centre become too low to support cell proliferation, and saturating when the rate at which nutrient-rich cells proliferate balances the rate at which necrotic cells are degraded [3]. For the parameter values chosen here, the tumour reached its steady state approximately 150

PLOS COMPUTATIONAL BIOLOGY
hours after beginning the simulation, with a steady state radius comparable to those observed in the experiments by Dorie et al. [46]. Fig 3 demonstrates that the model reproduces the qualitative behaviour of the radial distribution of infiltrating microbeads which was observed in [46]. These distributions resemble a wave which becomes increasingly dispersed as it moves radially inwards from the spheroid edge. The parameter values used in Fig 3 (ω q = 0.7, ω h = 0.1, τ = 8,t ¼ 16) were chosen to produce a distribution which matches that described in [46]. We remark that spheroids generated using different parameter values may produce microbead distributions which are not visually distinguishable (see Supporting Information, S4 Appendix: Microbead infiltration histograms for other parameter combinations). In addition to reproducing the dynamics of infiltrating microbeads in the experimental data described in [46], this model can also reproduce the observed distributions of 3 H-labelled cells (see Supporting Information, S5 Appendix: Replication of 3 H-labelled cells experiments). In particular, we note that, as for microbead infiltration patterns, 3 H-labelled cell infiltration patterns are strongly affected by variation of the parameters which control tumour cell proliferation and death. We can identify regimes in which the distribution of 3 H-labelled tumour cells closely resembles that described by Dorie et al. In contrast to previous models of 3 H-labelled cell infiltration, no additional mechanisms are required to describe differences in the distributions of 3 H-labelled cells and microbeads. Instead, differences in their distributions can be attributed to proliferation of the 3 H-labelled cells which causes the peak of the distribution to remain localised near the spheroid boundary where oxygen levels and, hence, cell proliferation rates are maximal. (For further details, see Supporting Information, S5 Appendix: Replication of 3

Spheroids with the same equilibrium size may have different spatial compositions
We systematically vary three key model parameters: the average cell cycle length, τ, and the oxygen thresholds for quiescence and hypoxia, ω q and ω h . We focus on these parameters as they are directly linked with the spheroid growth rate, the proportion of cells which are part of the proliferative rim, and the rate at which non-proliferating cells undergo cell death. As these parameters vary, we can generate spheroids which exhibit a range of behaviours (see Fig 4).
Panel A of Fig 4 shows spheroid growth curves for simulations generated from three different parameter sets. Each simulated spheroid reaches a steady state radius of approximately 14 cell diameters. Simulations i and ii are both generated with ω q = 0.6, and different values of τ and ω h . While the cells in Spheroid ii proliferate more often than those in Spheroid i, this effect is offset by the change in ω h which causes cells in Spheroid ii to become necrotic at higher oxygen levels than those in Spheroid i. This difference can be seen in Fig 4, Panel B, which shows that Spheroid ii has a more pronounced necrotic core and thinner quiescent region than Spheroid i.
Similarly, Spheroid i and Spheroid iii have the same average cell cycle length τ = 31.25, but differ in both ω q and ω h . While ω q has been lowered from ω q = 0.6 to ω q = 0.47, they realise similar equilibrium sizes because ω h increases from ω h = 0.34 to ω h = 0.42. These parameter changes cause Spheroid iii to have a thicker proliferating rim and narrower quiescent region than Spheroid i.
By appropriately tuning the proliferation rate and oxygen thresholds at which cells become quiescent or die, we can generate synthetic spheroids of the same equilibrium size which possess different internal compositions. This suggests that spheroid composition cannot be predicted by observing the overall growth dynamics alone. Panels C and D of Fig 4 show the effect of varying τ and ω q on the steady state size and composition of tumour spheroids, averaged For each fixed value of ω h , lowering ω q causes the proliferating rim to become wider and the spheroid to become larger. Decreasing the average cell cycle length τ causes the spheroid radius to increase, but this is due mainly to increases in the hypoxic and necrotic volumes of the spheroid. Results comparing the model with the 3 H-labelled data from [46] can also be found in the Supporting Information, S5 Appendix: Replication of 3

Hlabelled cells experiments.
https://doi.org/10.1371/journal.pcbi.1007961.g003 Cells are coloured according to their oxygen concentration and phenotype. By varying the thresholds ω q and ω h , and the growth rate τ, we can generate spheroids with similar growth dynamics and steady state radii which have different internal compositions. Spheroids i and ii have the same ω q , but differ in ω h and τ. Spheroids i and iii have the same τ, but differ in ω q and ω h . Spheroid i (blue line) contains a large number of quiescent cells and has a slow growth rate, but has a low oxygen threshold before cells become hypoxic. Cells in Spheroid ii (yellow line) proliferate more quickly than in Spheroid i, but since ω h is higher the spheroid has a larger necrotic core at steady state. Spheroid iii (green line) also has only a narrow quiescent region, but possesses a thicker proliferative rim than spheroid ii. C and D: Average spheroid sizes and compositions change as ω q and τ vary, averaged over 50 repetitions for each parameter set (witht ¼ 8 and ω h = 0.1 (C) or ω h = 0.3 (D)). Scale bars are 10 cell diameters in length. This result indicates that data on tumour size alone are not predictive of tumour composition. Additional data describing spatial structure and heterogeneity in tumour composition could improve the predictive power of models, particularly their ability to simulate responses to treatments such as radiotherapy which affect hypoxic and well-oxygenated cells differently.

Bead trajectories change as spheroid composition and growth dynamics vary
One advantage of using an agent-based framework rather than a continuum model to simulate passive infiltration [51][52][53] is that the trajectories of individual microbeads can be tracked in addition to the bead distribution (given by the infiltration histograms in Fig 3). We can use these trajectories to better understand why microbead distributions disperse while travelling radially inwards, but also to move beyond the original experiments in [46] to use beads to predict both spheroid composition and parameters associated with individual spheroids.
Panel A of Fig 5 shows the radial trajectory of 15 randomly selected beads from spheroids simulated with the same parameter set (ω q = 0.5, ω h = 0.3, τ = 16,t ¼ 8). Microbeads are coloured according to the oxygen concentration at their position. Fig 5 shows that microbead trajectories typically have two distinct phases. Initially, microbeads remain in the proliferating rim close to the spheroid boundary where their movement is characterised by random movement generated by tumour cell proliferation. During cell division, the new daughter cell is placed randomly in a neighbourhood of the parent cell. This position may be closer to the spheroid centre than the parent cell, or closer to the spheroid edge. Since the new location is selected randomly, the forces which act on a microbead that is close to a dividing cell due to the presence of the new cell are as likely to push the bead radially inwards or outwards. Since we expect the forces directed outwards to be comparable with those directed inwards over time, we expect radial movement in this regime to be dominated by Brownian motion and neglect movement in the tangential direction in our analysis. On leaving the proliferating region, the microbeads typically move on a linear trajectory towards the spheroid centre. Although the time microbeads remain in the proliferative rim varies, once they enter the central region they appear to move radially inwards at the same velocity. For each bead i, we define its radial infiltration velocity, V i r , to be its average radial velocity from when the bead enters the second regime until it reaches the necrotic core. Similarly, its waiting time, T i wait , is the time taken for a microbead to leave the outer proliferating rim and to enter the inner region. We denote the average of these quantities as V r ¼ 1 where n is the number of beads.
We use observations of the microbead velocities to identify the distance from the spheroid boundary at which cell movement transitions from Brownian dominated to advection dominated. The velocities from the first 5 hours of simulation are used to estimate D est , the estimated diffusion coefficient for Brownian motion. The expected displacement of a particle undergoing 1D Brownian motion for a time t is ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2D est t p , and we use this distance to estimate the location of the threshold separating the two regimes. As this value is close to 0 when t is small, we offset this boundary to ensure that microbeads do not cross this threshold during the first timestep. If a spheroid has radius R at time t after adding microbeads, we estimate that the threshold between regimes is Threshold ¼ 0:9R À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2D est t p . This estimated boundary is shown as a dashed black line in Panel A of Fig 5. Panel Panel D of Fig 5 shows how varying ω q and τ affects T wait , the mean time that a microbead spends in the random motion dominated regime. Reducing τ lowers T wait , indicating that bead movement across this part of the spheroid is driven by cell proliferation. Reducing ω q also increased the time spent near the spheroid edge for small τ, although this effect is less prominent as τ increases. . We use a 2D metric, as our simulations are in 2D; a similar 1D metric based on the relative widths of the tumour spheroid and quiescent compartment produces similar results. Panels B and C illustrate that the quiescent proportion can be estimated accurately by observing T wait and V r for a given simulation.

Inferring spheroid composition and parameters from microbead trajectories
In panel B of Fig 6 we present results from simulations involving a range of spheroids generated using randomly selected values of ω q and τ. The markers are coloured according to the quiescent proportion of each simulated spheroid, and are scattered according to the values of T wait and V r observed in each simulation. Spheroids with similarly sized quiescent regions are found close together, suggesting that measurements of individual microbead velocities and waiting times could be used to predict the composition of a tumour spheroid. The same effect can be observed for spheroids with similarly sized proliferating, hypoxic and necrotic regions (see S6 Appendix: Predicting spheroid composition from microbead observations).
We applied a k-nearest neighbours classifier [61] (with k = 10) to predict the composition of the spheroid based on the observed values of T wait and V r . The simulations were split into two groups to use as training and testing data for the classifier (100 data points for training, 48 data points for prediction, shown as circles and crosses respectively in Fig 6, panel B). Panel C of Fig 6 shows that the predictions of this classifier are accurate, with a high R 2 value (R 2 = 0.86) indicating that the predictions are well correlated with the true quiescent proportions of the simulated spheroids in the validation set. Similarly, it is possible to predict the size of the proliferating, hypoxic and necrotic regions (see S6 Appendix: Predicting spheroid composition from microbead observations).
As well as predicting the spheroid composition, T wait and V r can also be used to infer parameters associated with spheroid growth. Panel D of Fig 6 shows scatterplots in which the true values of τ and ω q for simulations in the validation dataset are compared with the predictions from the nearest neighbours classifier. Observing the trajectories of individual beads allows accurate prediction of the average cell cycle length (R 2 = 0.88) and oxygen threshold for quiescence (R 2 = 0.85) of simulated tumour cells.

Discussion
We have developed a hybrid, off-lattice ABM for oxygen-limited spheroid growth. Our model reproduces the sigmoidal growth dynamics seen in in vitro spheroids as well as their spatial structure, consisting of an outer rim of proliferating cells and a central necrotic core (Fig 2). Our simulations reveal that it is possible to generate synthetic spheroids with similar growth dynamics but different spatial compositions of proliferating, quiescent, hypoxic and necrotic cells (Fig 4). These changes in spheroid composition are driven by variations in cell-scale behaviours such as the average length of the cell cycle or the oxygen thresholds at which cells are affected by hypoxia.
The model describes the passive infiltration of microbeads into tumour spheroids and reproduces the distributions of microbeads described by Dorie et al. [46] (Fig 3). It also reproduces the dynamics of infiltrating 3 H-labelled cells described in [46], and can explain differences in the spatial distributions of the microbeads and 3 H-labelled cells without invoking additional mechanisms such as chemotaxis (see Supporting Information, S5 Appendix: Replication of 3 H-labelled cells experiments). In contrast to previous continuum models describing this data, this model can distinguish the trajectories of individual microbeads. This allows us to derive experimentally measurable properties of microbead infiltration, T wait and V r , which Observations of T wait and V r can be used to predict spheroid composition and simulation parameters. A: Schematic showing how the composition of the spheroid was defined in order to make predictions. The proportion of the spheroid which was quiescent is defined as the area of the quiescent region, A q , divided by the total spheroid area A total . (Comparable results are obtained when defining the quiescent proportion of the spheroid as the width of the quiescent annulus divided by the total spheroid radius.) B: Proportion of the spheroid radius accounted for by the quiescent compartment for simulations with randomly selected ω q and τ (with fixed ω h = 0.3 andt ¼ 8). Approximately two-thirds of the simulations (100 out of 148, marked as circles) were used to train a k-nearest neighbours classifier to predict the composition of spheroids in the remaining simulations (48 out of 148, marked as crosses). C: Results of using 10-nearest neighbours classification to predict the proportion of each tumour spheroid which was quiescent based on observed values of T wait and V r for microbeads in each simulation. D: permit accurate prediction of both model parameters and the spheroid composition. T wait describes the average length of time which a microbead spends close to the spheroid edge, and correlates with the average length of the cell cycle, τ (Fig 5). V r is the speed at which microbeads move radially inwards once they have moved sufficiently far from the spheroid boundary to enter the advection dominated regime inside the spheroid, and correlates with both τ and the oxygen threshold at which cells become quiescent, ω q (Fig 5). Combinations of V r and T wait can be used to accurately predict the values of τ and ω q used to simulate a spheroid (Fig  6), as well as the proportion of the spheroid consisting of proliferating, quiescent, hypoxic or necrotic cells (Fig 6 and Supporting Information, S6 Appendix: Predicting spheroid composition from microbead observations).
The model enables detailed simulation of the forces acting on microbeads during passive migration. In common with other agent-based models, it includes a large number of parameters. The results presented here were obtained by fixing many of the parameters at values specified in S7 Appendix: Table of parameters. While some of these parameters have been informed by the literature, this model has identified those that should be informed by further experiments to identify the relationship between spheroid simulations and in vitro spheroids. We restrict our analysis to 2D tumour spheroids, for several reasons. A single realisation of a 2D simulation may take several hours to run on one core of a desktop PC; by contrast a 3D simulation of a spheroid containing tens of thousands of cells would take several days to complete. This computational effort would make it prohibitive to conduct in 3D the parameter sensitivity analyses performed in this paper, particularly as multiple simulations are required for each parameter set. Further, the 2D simulations described here are sufficiently detailed to identify qualitative differences in passive migration patterns as parameters relating to tumour cell proliferation and death are varied. Identifying the relationships between models of tumour spheroid growth implemented in 2D or 3D, or implemented within different software frameworks, remains an open problem.
We note that in the necrotic core of some spheroids observed in in vitro experiments, loss of volume in dead cells due to fluid leakage can cause "cracking", or fluid filled voids, when coupled with intercellular adhesion [62,63]. While this behaviour has been reported in existing agent-based models of tumour spheroids [32,37], we note that for our simulations this effect is only observed when the magnitude of the surface tension force β is close to zero. As predicted by Landman and Please [64], spheroids featuring this "void" effect frequently produce travelling wave solutions which do not reach a steady state when this force resisting spheroid growth is sufficiently small. Extensions of our model could include the response of tumour spheroids to treatments such as radiotherapy or chemotherapy which are most effective at targetting proliferating cells. Our agent-based model could be used to examine the role that tumour composition plays during these therapies, as previous continuum models predict that treatment response is highly dependent on the internal composition of tumours [12]. Advection has also been identified as being of importance in models of drug treatment efficiency (e.g., [49,50]). Advection-diffusion equations have also been used to model the growth of vascular tumours [65], and the development of tumour metastasis [66].
This model could also be extended to investigate tumour-immune interactions in which immune cells, such as macrophages, infiltrate into clusters of tumour cells. Previous continuum models of macrophage infiltration into tumour spheroids and avascular tumours Comparison of parameter values for (left) τ and (right) ω q based on 10-nearest neighbours classification with the values used to generate the synthetic data.
https://doi.org/10.1371/journal.pcbi.1007961.g006 emphasise the role of chemotaxis [67][68][69][70][71] or paracrine signalling [72,73]. However, to date the influence of passive migration on the infiltration of immune cells into spheroids and avascular tumours has not been considered. Much focus has been placed on the impact that immune cells have on tumour cells, but the impact of solid tumours on the spatial distributions of immune cells has attracted less attention.
In future work we could also extend the model to account for other metabolites and waste products which impact tumour spheroid growth, such as glucose, lactate and pH levels [21]. These model extensions could be used to simulate tumour responses to metabolic inhibitors or treatments that are affected by pH.
We have used an agent-based framework to simulate passive migration into tumour spheroids. This migration is caused by the collective movement of tumour cells, which is itself a consequence of the spatial composition of a tumour spheroid. This model shows that passive migration is strongly influenced by spheroid composition, and shows how microbead infiltration experiments similar to those conducted by Dorie et al. [46] could be used to deduce this composition.