A stochastic mathematical model of 4D tumour spheroids with real-time fluorescent cell cycle labelling

In vitro tumour spheroids have been used to study avascular tumour growth and drug design for over 50 years. Tumour spheroids exhibit heterogeneity within the growing population that is thought to be related to spatial and temporal differences in nutrient availability. The recent development of real-time fluorescent cell cycle imaging allows us to identify the position and cell cycle status of individual cells within the growing spheroid, giving rise to the notion of a four-dimensional (4D) tumour spheroid. We develop the first stochastic individual-based model (IBM) of a 4D tumour spheroid and show that IBM simulation data compares well with experimental data using a primary human melanoma cell line. The IBM provides quantitative information about nutrient availability within the spheroid, which is important because it is difficult to measure these data experimentally.


Introduction
In vitro tumour spheroids are widely adopted to study avascular tumour growth and drug design [1][2][3]. Unlike simpler two-dimensional assays, tumour spheroid experiments exhibit heterogeneity within the growing population which is thought to be partly driven by spatial and temporal differences in the availability of diffusible nutrients, such as oxygen [3,4]. Historically, tumour spheroid experiments are analysed using bright-field imaging to measure the spheroid size [5,6], however, this does not reveal information about the internal structure of the spheroid. Since 2008, fluorescent ubiquitination-based cell cycle indicator (FUCCI) has enabled real-time identification of the cell cycle for individual cells within growing populations [4,7,8]. Using FUCCI, nuclei of cells in G1 phase fluoresce red, nuclei of cells in S/G2/M phase fluoresce green, and nuclei of cells in early S (eS) phase appear yellow as a result of both red and green fluorescence being active [7] (figure 1a). FUCCI provides information about spheroid size and heterogeneity of the cell cycle status (figure 1c-e). At early times the entire spheroid is composed of freely cycling cells, with a relatively even distribution of FUCCI colours, whereas at intermediate times cells in the central region becomes predominantly red, indicating G1-arrest [4]. Late time growth is characterized by the formation of a central necrotic region, indicated by the absence of fluorescence. FUCCI allows us to identify both the position of individual cells in three spatial dimensions, as well as identifying cell cycle status, leading to the notion of a four-dimensional (4D) tumour spheroid [9]. Assuming spherical symmetry, the geometry of 4D spheroids can be characterized by three radii: (i) r o ðtÞ . 0 is the outer radius; (ii) r a ðtÞ ! 0 is the arrested radius; and, (iii) r n ðtÞ ! 0 is the necrotic radius, with r o ðtÞ . r a ðtÞ ! r n ðtÞ. In figure 1e, we see that r n ðtÞ ¼ 0 for t ≤ 3, with the necrotic core forming sometime between t = 3 and t ¼ 6 days.
Continuum mathematical models of tumour spheroids have been developed, analysed, and deployed for over 50 years [10][11][12][13][14][15], and these developments have included very recent adaptations of classical models to study tumour spheroids with FUCCI [9,16,17]. Continuum modelling approaches lack the ability to track individual cells within the growing population, and typically neglect heterogeneity and stochasticity. In comparison, individual-based models (IBMs) allow us to study population dynamics in more detail, by keeping track of all individuals and explicitly capturing heterogeneity and stochasticity [18,19]. Over the last 20 years, as computing power has increased at the same time that experimental imaging resolution has improved, there has been an increasing interest in interpreting tumour spheroid experiments using IBMs [20][21][22][23][24][25], with some studies using these models to focus explicitly on how the balance of cell migration and cell proliferation impacts phenotype selection [26]. Flegg & Nataraj [27] succinctly review mathematical modelling methodologies used to interpret tumour spheroid experiments. While several previous IBMs have the ability to track the cell cycle within individuals [20][21][22][23][24][25], our aim is to track the cell cycle in a relatively minimal IBM and to quantitatively interpret this information in terms of a new set of 4D tumour spheroid experiments.
In this work, we develop an IBM of 4D tumour spheroid growth with FUCCI. The IBM describes how individual cells migrate, die and progress through the cell cycle to mimic FUCCI. Certain mechanisms in the IBM are coupled to the local availability of a diffusible nutrient. The biological fidelity of the IBM is demonstrated by qualitatively comparing simulation results with detailed experimental images at several cross sections (figure 1b). Quantitative data from the model are used to assess the spheroid population distribution, nutrient concentration, and the role of variability in spheroid growth. We extract and quantitatively compare simulation radius estimates with measurements from a series of 4D tumour spheroid experiments using a human primary melanoma cell line (figure 1). Using a careful choice of parameter values, we show that the IBM quantitatively replicates key features of 4D tumour spheroids. Experimental images of a tumour spheroid using the human melanoma cell line WM793B at days 0, 3, 6 and 10 (after formation) showing. (c) Full spheroids, viewed from above, (d ) spheroid hemispheres and (e) spheroid slices, where the cross section is taken at the equator. White dashed lines in (e) denote the boundaries of different regions, where the outermost region is the proliferative zone, the next region inward is the G1-arrested region, and the innermost region at days 6 and 10 is the necrotic core. In (a) and (d ), we use cyan colouring for dead cells, which assist in identifying the necrotic core in (d ). Spheroid outer radii are labelled alongside their corresponding time points, and scale bars represent 200 μm.

Methods
Experimental methods are described in electronic supplementary material, §S1.

Individual-based mathematical model
We simulate 4D spheroid growth inside a cubic domain, V, of side length L, where L is chosen to be large enough so that agents do not reach the boundary of the domain during the simulation, but not so large as to incur significant computational overhead (electronic supplementary material, S4.3). Biological cells are represented as discrete agents located at x n (t) = (x n (t), y n (t), z n (t)) for n = 1, 2, 3, …, N(t), where N(t) is the total number of agents at time t.

Gillespie algorithm
The IBM describes key cellular-level behaviours; namely cell cycle progression and mitosis, cell motility, and cell death, as discrete events simulated using the Gillespie algorithm [28]. Each agent has an allocated rate of cell cycle progression, dependent on its cell cycle status and the local nutrient concentration (figure 2a). Agents in each phase of the cell cycle are coloured according to FUCCI, with G1 agents coloured red, eS agents coloured yellow, and S/G2/M agents coloured green.
We make the natural assumption that biological cells require access to sufficient nutrients to commit to entering the cell cycle. Therefore, the red-to-yellow transition rate, R r (c), depends on the local nutrient concentration, c(x, t) (figure 2a). Once an agent has committed to entering the cell cycle, we assume the yellow-togreen transition takes place at a constant rate R y , and the greento-red transition, which involves mitosis, occurs at a constant rate R g (figure 2a).
The rates of agent migration and death, m(c) and d(c), respectively, are assumed to depend on the local nutrient concentration. When an agent dies, it is removed from the simulation and we record the location at which the death event occurs (figure 2c). When an agent moves or undergoes mitosis (figure 2d,e), a random direction in which the agent will migrate, or its daughter agents will disperse, is chosen (figure 2b). For an agent undergoing mitosis, the first daughter agent is placed a distance σ/2 along the randomly chosen direction, and the second daughter agent is placed at a distance σ/2 in the opposite direction, leaving the two daughter agents dispersed a distance of σ apart, where we set σ to be equal to a typical cell diameter [29] (figure 2d and table 1). When migrating, agents are displaced a distance μ along the randomly chosen direction (figure 2e). Similar to the dispersal, we simulate migration by taking the step length μ to be a typical cell diameter.
We specify the agent cycle progression rates, ð2:1Þ R y ðcÞ ¼ R y , ð2:2Þ where c(x n , t) ∈ [0, 1] is the non-dimensional nutrient concentration at the location of the nth agent; R r > 0 is the maximum red-to-yellow transition rate; m max > m min ≥ 0 are the maximum and minimum migration rates, respectively; d max > d min ≥ 0 are the maximum and minimum death rates, respectively; η 1 > 0, η 2 > 0 and η 3 > 0 are Hill function indices; and c a > 0, c m > 0 and c d > 0 are the inflection points of R r (c), m(c) and d(c), respectively (figure 2a).

Nutrient dynamics
We make the simplifying assumption that cell migration, death and progression through the cell cycle are regulated by a single diffusible nutrient, such as oxygen [4,10,12]. The spatial and temporal distribution of nutrient concentration, C(x, t), is assumed to be governed by a reaction-diffusion equation To solve this reaction-diffusion equation, we set v(x i,j,k ,t) = N i,j,k /h 3 , where N i,j,k is the number of agents within the control volume surrounding the node located at (x i , y j , z k ) and h 3 is the volume of the control volume. On the boundary, @V, we impose C = C b , where C b is some maximum far-field concentration.
Our experiments lead to spheroids of diameter 500-600 μm over a period of 10 days after spheroid formation (figure 1) (14 days after seeding). Since these length and timescales are clear, we leave the independent variables x and t in equation (2.6) as dimensional quantities. By contrast, spatial and temporal variations of C(x, t) are very difficult to measure during spheroid growth, so we non-dimensionalize the independent variable c(x, t) = C(x, t)/C b , giving Typically, the timescale of nutrient diffusion is much faster than the timescale of spheroid growth [10]. Consequently, we approximate equation (2.7) by . Therefore, we describe the spatial and temporal distribution of nutrients by solving equation (2.8) repeatedly during the simulation. This quasi-steady approximation is computationally convenient, as we describe later. We solve equation (2.8) with a finite volume method on a uniform structured mesh (electronic supplementary material, S4) with node spacing h. The IBM treats cells as point particles, which we call agents. Interactions between agents are modelled implicitly by specifying rates of migration and cell cycle progression that depend upon the local nutrient concentration. As we will show later, this very simple framework replicates our experimental observations reasonably well, however, we acknowledge there are other ways to treat agent-agent interactions, including latticebased hard-core exclusion [29] and lattice-free methods for minimizing [30] or avoiding agent overlap [23,24]. Other options for modelling agent-agent interactions include direct simulation of agent adhesion or repulsion [31]. Here, our modelling philosophy is to work with the simplest possible biologically plausible model that can describe the key features of interest. Therefore, we deliberately avoid these additional mechanisms in our IBM framework.

Simulation algorithm
We simulate spheroid growth by supposing the spheroid initially contains N(0) agents distributed uniformly within a sphere of radius r o ð0Þ . 0 ½mm. While it is experimentally relevant to assume the population is spherically symmetric at t = 0, this assumption is not necessary, and we will discuss this point later. The proportion of agents chosen to be red, yellow, or green at t = 0 can be selected arbitrarily, but we choose these proportions so that the internal structure and composition of the in silico spheroids are consistent with our in vitro measurements. We achieve this by choosing the initial numbers of red, yellow, and green agents, N r (0), N y (0) and N g (0), respectively, noting that ). The most appropriate timescale for individual cell-level behaviour is hours, however, spheroid development takes place over 10 days, so we will use a mixture of timescales to describe different features of the experiments and simulations as appropriate. We simulate spheroid growth from t = 0 to t = T h, updating the nutrient concentration at M equally spaced points in time. This means that the nutrient concentration is updated at intervals of duration t* = T/M [h]. The accuracy of our algorithm increases by choosing larger M (smaller t*), but larger M decreases the computational efficiency. We explore this tradeoff and find that setting t* = 1 h is appropriate (electronic supplementary material, S4.4). When equation (2.8) is solved for c(x, t), the value of c(x n , t) at each agent is calculated using linear interpolation. These local nutrient concentrations are held constant for each agent while resolving all the various agent-level events (cycling and proliferation, migration and death) from time t to time t + t*. After resolving the appropriate agent-level events, we update the agent density before updating the nutrient profile again. Pseudo-algorithms for the IBM are provided (electronic supplementary material, S9), and code to reproduce key results is available on GitHub https://github. com/ProfMJSimpson/4DFUCCI.

IBM image processing
To estimate r o ðtÞ, r a ðtÞ and r n ðtÞ, we apply methods described in [16,32,33] to the IBM output. Briefly, we import the agent locations from a particular cross section, and map these locations to an (L + 1) × (L + 1) pixel image, increase the size of the agents to 12 pixels in diameter, and use edge detection to identify and estimate r o ðtÞ, r a ðtÞ and r n ðtÞ (electronic supplementary material, S2). This procedure adapts the image processing approach for the experimental images so that it is applicable to the synthetic results from the IBM.

Results and discussion
We now compare and analyse images and measurements from a range of in vitro experiments and in silico simulations. All experiments use the WM793B melanoma cell line, which takes approximately 4 days to form spheroids after the initial seeding in the experiments [17]. This means that t ¼ 0 days corresponds to 4 days after seeding to give the experimental spheroids sufficient time to form. Snapshots from the IBM correspond to a single realization, however, time-series data from the IBM are reported by simulating 10 realizations of the IBM and then averaging appropriate measurements across the 10 simulations. While we take great care to analyse the simulation images and the experimental images using the same imageprocessing algorithm, the format of the experimental images is different to the format of the images produced by the IBM, and we do not explicitly account for these differences.

Qualitative comparison of experiments and simulations
We  [17,22], we aim to provide more comprehensive information about the internal structure of the spheroid by making comparisons at multiple locations. At the beginning of the experiment, in all cross sections (in vitro and in silico), we see the population is relatively uniform, with an even distribution of colours, suggesting the entire spheroid is composed of freely cycling cells. At t = 2 and t ¼ 4 days, however, we begin to see the development of heterogeneity within the growing in vitro and in silico populations, with those cells and agents at the centre of the growing spheroid predominantly red, indicating G1-arrest. By t ¼ 4 days, we see the value of comparing different cross sections, since the G1-arrest is clear in the centre of the equatorial cross section, but there is no obvious heterogeneity present across either the upper or lower cross sections. Similarly, by t ¼ 6 days, we see the formation of a necrotic core in the equatorial cross section, but this is not present at either the upper or lower cross sections. By t = 8 and t = 10 days, the spheroid has developed into a relatively complicated heterogeneous structure; the outer spherical shell contains freely cycling cells, the intermediate spherical shell contains living G1-arrested cells, and the internal region does not contain any fluorescent cells.
Overall, the qualitative match between the IBM and the experiment confirms that the IBM captures both the macroscopic growth of the entire spheroid, as well as the emergent spatial and temporal heterogeneity. We now build on this preliminary qualitative information by comparing quantitative measurements of growth of the spheroid.

Spheroid structure and nutrient profiles
Given the ability of the IBM to capture key spatial and temporal patterns of spheroid growth, cell cycle arrest, and cell royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20210903 death throughout the spheroid, we now demonstrate how to take these preliminary simulations and extract detailed quantitative data that would be difficult to obtain experimentally. Figure 4a shows a typical IBM simulation during the interval where we observe the development of internal structure. For clarity, we plot the locations of all living agents as in figure 3 together with the locations at which agents die, which is difficult to estimate experimentally, but is straightforward with the IBM. Each spheroid in figure 4a is shown with an octant removed to highlight the development of the internal structure, and for further clarity we show equatorial cross sections in figure 4b.
To quantify the internal spheroid structure we simulate 10 identically prepared realizations of the IBM and extract averaged quantitative data that are summarized in figure 4c (electronic supplementary material, S10). These data include plotting the non-dimensional nutrient concentration, c(x, t), and various normalized agent densities, ϱ( p(t), t), as a function of distance from the spheroid periphery, pðtÞ ¼ r o ðtÞ À r, where r is the distance from the spheroid centre. Hence, p(t) = 0 corresponds to the spheroid periphery, and pðtÞ ¼ r o ðtÞ corresponds to the spheroid centre. This representation of internal spheroid structure is made by assuming that the growing population remains spherically symmetric, which is a reasonable assumption since our initial condition and spheroid growth is spherically symmetric (figure 4a). Each density profile is normalized relative to the maximum value of all agent densities across all time points, so that we can compare how the density of the various subpopulations of agents and nutrient are distributed (electronic supplementary material, S10). Using the IBM we are able to describe the spatial and temporal densities of living agents in various phases of the cell cycle (G1, eS and S/G2/M) as well as G1-arrested agents. We plot each density profile as a function of the distance from the periphery as this allows us to compare various profiles as the size of the spheroid increases [9,34].
Averaged relative agent density profiles from the IBM provide quantitative information that cannot be easily obtained from experimental observations. Initially, we see the relatively evenly distributed G1, eS and S/G2/M populations become rapidly dominated by agents in G1 phase, which then form an obvious inner-most arrested region by about t = 2 days. During the interval 3 < t < 6 days, we  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20210903 observe rapid growth in the arrested population. During the later interval, 6 < t < 10 days, we see the formation of a clear necrotic core. These results indicate the spatial and temporal role of stochasticity, with the variability most evident in the G1 and arrested G1 populations at early and intermediate times. Plotting the relative agent densities in this way provides a simple approach to interpret the spatial and temporal organization of cell cycle status within the growing spheroid and visualizing the agent densities together with the non-dimensional nutrient concentration is particularly useful when this kind of information cannot be easily obtained experimentally. In particular, it is technically challenging to measure absolute concentrations of nutrient profiles during these experiments [15,35,36] and so we now focus on visualizing the nutrient concentration profile that drives this heterogeneity.
Results in figure 5 show spatial and temporal patterns in the nutrient profile, c(x, t), for a typical IBM simulation from figure 4. Figure 5a shows the three-dimensional evolution of c(x, t), with the colour bar highlighting the death and arrest thresholds, c d and c a , respectively. These three-dimensional plots show the depletion of nutrient over time in the central region of the spheroid, leading to strong spatial gradients of nutrient concentration near the edge of the growing spheroid. Profiles in figure 5b show the nutrient profile at the equatorial plane with the c(x, y, 0) = c a contour (red) and the approximate size of the necrotic core (cyan) superimposed. Simplified one-dimensional profiles of c(x, t), along x = (x, 0, 0), are shown in figure 5c, where the diameter of the growing spheroid (Àr o ðtÞ , x , r o ðtÞ) is shaded in yellow. Again, these simplified cross sections illustrate how nutrient consumption leads to the formation of spatial nutrient gradients near the outer radius of the growing spheroid. Overall, a key strength of the IBM is the ability to extract agent-level information (figure 4) as well as information about the nutrient distribution (figure 5), whereas experimental studies typically report cell-level data without explicitly showing nutrient-level information [4,6].
While it is very difficult to measure the spatial and temporal distribution of diffusible nutrient experimentally in the growing spheroid, it is possible to indirectly examine our assumption that spatial and temporal differences in cell cycle status are partly driven by the availability of oxygen. Figure 6 shows a series of spheroids stained with pimonidazole and pimonidazole-detecting antibodies, which indicate hypoxia [37]. In this series of images, we see evidence of hypoxia staining in the central region of the spheroid at t = 0, with persistent hypoxia staining adjacent to the necrotic core at later times. These results support our hypothesis that spatial and temporal differences in nutrient availability correspond with spatial and temporal differences in cell cycle status, and in this case the pimonidazole staining suggests that oxygen availability plays a role in the development of heterogeneity within the growing population. While this observation is consistent with our IBM, it does not rule out the possibility of multiple diffusible signals acting in unison, and we will discuss this possibility later.

Role of variability
Experimental images (figures 1, 3 and 6) suggest that spheroid development is variable, as we see spheroids of different  (c) Relative concentrations ϱ( p, t) of nutrient (black) and cycling red, yellow and green agents (coloured appropriately), based on distance from the periphery pðtÞ ¼ r o ðtÞÀ r, averaged over 10 identically prepared simulations. The dashed red line shows the relative density of arrested red agents, also averaged over 10 simulations with identical initial conditions. For nutrient, ϱ( p, t) = c. For agents, ϱ( p, t) is the relative agent density (electronic supplementary material, S10). Shaded areas represent plus or minus one standard deviation about the mean, and are non-zero as a consequence of stochasticity in the model, even though the 10 simulations start with identical populations and radii.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20210903 diameters at the same time points. One of the limitations of relying on experimentation alone is that it can be difficult to quantify the importance of different sources of variability, whereas this can be assessed very simply with the IBM. For example, we can simulate multiple spheroids that start from precisely the same initial condition to quantify the variability that arises due to the stochastic growth process, or we can deliberately introduce variability into the initial composition of the spheroid to explore how this variability evolves during spheroid growth for a suite of simulated spheroids.
Simulation data in figure 7a show the temporal evolution of various agent subpopulations, including the total number of living agents, dead agents, G1, eS, S/G2/M and G1-arrested agents. Each profile shows the mean number of agents obtained by simulating 10 identically initialized spheroids with r o ð0Þ ¼ 245 mm, which matches the average spheroid diameter at t = 0 days in the suite of in vitro experiments. The variability in these profiles is quantified by calculating the sample mean and sample standard deviation and shading the region corresponding to the sample mean  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20210903 plus or minus one sample standard deviation, and we see that, at this scale, the variability is barely noticeable. By contrast, results in figure 7b show equivalent data from a suite of simulations where the initial density of agents in the spheroid is held constant, but the initial radius of the 10 simulated spheroids is deliberately varied to mimic the observed variability in our experiments. The initial radius in each simulation corresponds to one of 10 particular experimental measurements (figure 7). The mean of these 10 initial radii measurements is r o ð0Þ ¼ 245 mm, which is precisely the same as the initial radii for the simulations in figure 7a. Comparing results in figure 7a,b shows that the average population profiles are very similar, but the variability is strikingly different. This exercise shows that quantifying the variability in spheroid size at the beginning of the experiment is the key to understanding and predicting the variability in spheroid composition and size at the end of the experiment. The clear differences we see in the variability between results in figure 7a,b confirms that 10 simulations is sufficient to explore the role of variability at the beginning of the experiment. It is also interesting to note that these simulation results are consistent with our previous observations. For example, the in vitro spheroids in figure 3 have r o ð0Þ ¼ 232 mm and we see that it takes until t = 6 days for a clear necrotic core to form in the equatorial cross section. By contrast, the spheroid in figure 6 is larger with r o ð0Þ ¼ 274 mm and we see a clear necrotic core at t ¼ 3 days. This highlights the importance of taking great care with measurements at the beginning of the experiment [17].

Quantitatively matching experimental and mathematical spheroids
Results in figure 8 compare the temporal evolution of r o ðtÞ, r a ðtÞ and r n ðtÞ, from our suite of experiments and simulations.
The data in figure 8 show the value in working with a stochastic model since the experimental measurements are quite variable, with estimates of r a ðtÞ and r n ðtÞ more variable than estimates of r o ðtÞ. This difference in variability is because we measure r o ðtÞ automatically every 6 h, whereas measurements of r a ðtÞ and r n ðtÞ require manual harvesting, fixing and imaging, and so we obtain daily measurements only. Similarly to §3.4, we compare experimental results of average data in simulations with and without explicit variability at t = 0. The experiment-IBM comparison in figure 8a corresponds to the case where we simulate 10 identically prepared realizations of the IBM, where each simulated spheroid has the same initial radius r o ð0Þ ¼ 245 mm, and we see that the average simulation results capture the average trends in the experimental measurements well, but the IBM simulations do not capture observed variability in the evolution of r a ðtÞ or r n ðtÞ. By contrast, the experiment-IBM comparison in figure 8b, where we deliberately mimic the experimental variability at t = 0, captures both the average experimental trends and variability in the experimental data quite well. Again, the difference between results in figure 8a,b suggests that accurately incorporating the initial variability in the experimental data is critical if we wish to capture the observed variability during the experiments with the IBM. Data in figure 8 show that r o ðtÞ increases approximately linearly, whereas the development of the internal structure is more complicated, with r a ðtÞ initially decreasing for the first day before growing at a similar rate as r o ðtÞ. The necrotic core does not form until approximately t = 4 days. While our IBM-experimental comparison in figure 8 suggests that the IBM can quantitatively capture experimental trends, we have obtained this match with a careful choice of parameters without undertaking a more rigorous parameter estimation exercise [38]. The sensitivity of the IBM predictions to the parameter values is briefly explored in the electronic supplementary material, S11.

Conclusion and future work
We developed an IBM that simulates 4D tumour spheroids with FUCCI. IBM simulations show we can successfully reproduce qualitative and quantitative patterns of spatial and temporal differences in cell cycle status observed experimentally. This heterogeneity is driven by spatial and temporal variations in nutrient availability, which we model using a reaction-diffusion equation. An important advantage of the IBM is our ability to extract and describe measurements that are difficult to obtain in vitro. In particular, we show how to visualize both the growing population within the spheroid together with the associated spatial patterns of nutrient concentration over time. Furthermore, the IBM makes it very simple to explore how various features contribute to the overall variability in spheroid development, and we find that relatively small variations in the initial size of the spheroid lead to relatively pronounced differences in spheroid size and composition at later times [17]. We conclude our investigation by showing that we can quantitatively match the spatial and temporal development of a series of in vitro 4D spheroids using the WM793B human primary melanoma cell line with a careful choice of parameters. We anticipate that tumour spheroids formed with different cell lines will be able to be simulated with our IBM, but will require different parameter values.
Overall, our modelling philosophy is always to work with the simplest possible mechanisms required to capture our experimental observations. Naturally, this means that there are many ways that the IBM can be extended. For example, here we make the simple assumption that spheroid growth is regulated by a single diffusible nutrient, which seems appropriate for our data. If, however, experiments show that it is important to consider multiple nutrients in unison, our IBM framework can be extended to deal with this. Similarly, we focus on symmetric spheroid growth to be consistent with our experiments, but it is straightforward to relax this assumption by specifying a different initial arrangement of agents, or by allowing asymmetric nutrient delivery by, for example, a blood vessel as we explore in electronic supplementary material, S12. Another point that could be revisited is that we implement the simplest possible cell migration mechanism where the direction of motion is random. While this assumption appears reasonable for our data, it is possible to bias the migration in response to either the nutrient concentration, the gradient of the nutrient concentration, or the density of agents by explicitly describing agent-to-agent interactions [30]. Each of these potential extensions could be incorporated into our IBM framework with the aim of potentially increasing the biological fidelity of the model. However, we caution against this simplistic approach since these mechanisms also increase the number of unknown model parameters required for simulation. To minimize issues with parameter identifiability, we prefer to work with a minimal model [38]. If, however, future experimental measurements indicate that our minimal assumptions need revising, our IBM framework is sufficiently flexible to incorporate such extensions. Another option for future refinement is to conduct a more thorough parameter estimation exercise [38]. Here, we carefully chose parameters so that data from the IBM matches our experimental data, but future analysis could include a more rigorous assessment of parameter estimation, and we leave this for future consideration.