Quantitative analysis of tumour spheroid structure

Abstract
Editor's evaluation
Introduction
Materials and methods
Results
Discussion
Appendix 1
Appendix 2
Appendix 3
Data availability
References
Article and author information
Metrics

Abstract

Tumour spheroids are common in vitro experimental models of avascular tumour growth. Compared with traditional two-dimensional culture, tumour spheroids more closely mimic the avascular tumour microenvironment where spatial differences in nutrient availability strongly influence growth. We show that spheroids initiated using significantly different numbers of cells grow to similar limiting sizes, suggesting that avascular tumours have a limiting structure; in agreement with untested predictions of classical mathematical models of tumour spheroids. We develop a novel mathematical and statistical framework to study the structure of tumour spheroids seeded from cells transduced with fluorescent cell cycle indicators, enabling us to discriminate between arrested and cycling cells and identify an arrested region. Our analysis shows that transient spheroid structure is independent of initial spheroid size, and the limiting structure can be independent of seeding density. Standard experimental protocols compare spheroid size as a function of time; however, our analysis suggests that comparing spheroid structure as a function of overall size produces results that are relatively insensitive to variability in spheroid size. Our experimental observations are made using two melanoma cell lines, but our modelling framework applies across a wide range of spheroid culture conditions and cell lines.

Editor's evaluation

In this work, the authors test the hypothesis that tumor spheroids initiated with different numbers of cells grow to similar limiting sizes. The authors use a combination of experimental and mathematical techniques to examine this hypothesis with two melanoma cell lines. The authors find that spheroid structure and size are relatively insensitive to variations in initial number of cells, and suggest this finding may generalize to other cell lines.

https://doi.org/10.7554/eLife.73020.sa0

Introduction

Three-dimensional tumour spheroids provide an accessible and biologically realistic in vitro model of early avascular tumour growth (Hirschhaeuser et al., 2010; Cui et al., 2017). Spheroids play a vital role in cancer therapy development, where the effect of a putative drug on spheroid growth is an indicator of efficacy (Smalley et al., 2008; Santiago-Walker et al., 2009; Alexander and Friedl, 2012; LaBarbera et al., 2012; Loessner et al., 2013; Beaumont et al., 2015; Langhans, 2018; Theard et al., 2020). In this context, reproducibility and uniformity in spheroid sizes is paramount (Ivascu and Kubbies, 2006; Friedrich et al., 2009; Eilenberger et al., 2021), yet variability in the initial and final spheroid size is rarely accounted for, meaning subtle differences go undetected. We address this by developing a mathematical and statistical framework to study spheroid structure as a function of size, allowing us to ascertain whether initial spheroid size significantly affects growth dynamics.

Compared with traditional two-dimensional cell culture, spheroids closely mimic an avascular tumour microenvironment where spatial differences in the availability of nutrients strongly influence growth (Mark et al., 2020). We observe that spheroids grow to a limiting size that is independent of the number of cells used to initiate the experiment (Figure 1a–f), leading us to hypothesise that spheroids have a limiting structure (Folkman and Hochberg, 1973). This behaviour is consistent with untested predictions of mathematical models of tumour progression (Greenspan, 1972; Adam and Maggelakis, 1990; Groebe and Mueller-Klieser, 1996; Byrne and Chaplain, 1998; Ward and King, 1999; Araujo and McElwain, 2004; Wallace and Guo, 2013; Sarapata and de Pillis, 2014; Flegg and Nataraj, 2019; Murphy et al., 2021; Figure 1g). Many mathematical models assume that spheroid growth eventually ceases due to a balance between growth at the spheroid periphery and mass loss at the spheroid centre, driven by the spatial distribution of nutrients and metabolites (Figure 1h; Greenspan, 1972; Gomes et al., 2016). We analyse highly detailed experimental data from a large number of spheroids to answer fundamental biological and theoretical questions. Firstly, we study the effect of initial spheroid size on the transient and limiting spheroid structure. The initial size of spheroids is often highly variable (Mark et al., 2020), yet is rarely accounted for in statistical analysis. Secondly, we study the relationship between spheroid size and structure using a mathematical model that describes growth inhibition due to the spatial distribution of nutrients and metabolites.

Figure 1

Download asset Open asset

Experimental data and mathematical model.

(**a–f**) Growth of WM983b and WM793b spheroids over three weeks, initiated using approximately 2500, 5000 and 10,000 cells. The solid curve represents average outer radius and the coloured region corresponds to a 95% prediction interval (mean ± 1.96 std). (**c–f**) Size distribution of WM983b spheroids at days 10, 14, 18, and 21 for each initial seeding density. (**g–h**) Dynamics of the Greenspan, 1972 model, which describes three phases of growth and the development of a stable spheroid structure under assumptions of nutrient and waste diffusion. We denote by $R$ the spheroid radius, $ϕ$ the relative radius of the arrested region and $η$ the relative radius of the necrotic core. (i) Optical sections showing three phases of growth in the experimental data (WM983b spheroids initiated with 2500 cells at days 3, 7, and 14). Colouring indicates cell nuclei positive for mKO2 (magenta), which indicates cells in gap 1; cell nuclei positive for mAG (green), which indicates cells in gap 2; and cell nuclei stained with DRAQ7 (blue), which indicates necrosis.

We study spheroids grown at three seeding densities from human melanoma cells (Herlyn et al., 1985; Herlyn, 1990) transduced with the fluorescent ubiquitination cell cycle indicator (FUCCI) (Sakaue-Sawano et al., 2008; Haass et al., 2014; Kienzle et al., 2017; Spoerri et al., 2021). FUCCI technology discriminates between cells in different stages of the cell cycle, namely gap 1 (before synthesis) and gap 2 (after DNA replication), allowing us to identify regions containing actively cycling cells, and regions where the majority of the cells are viable but in cell cycle arrest. We grow spheroids for up to 24 days to allow sufficient time to observe growth inhibition. We summarise experimental images using three measurements of spheroid structure: (1) the overall size of each spheroid; (2) the size of the inhibited region (which we define as the region where the majority of cells are in gap 1); and, (3) the size of the necrotic core.

It is widely accepted that the eventual inhibition of spheroid growth arises through three phases (Figure 1g and i; Wallace and Guo, 2013; Spoerri et al., 2017; Flegg and Nataraj, 2019). During phase 1, for spheroids that are sufficiently small, we observe cycling cells throughout. In phase 2, spheroids develop to a size where cells in the spheroid centre remain viable but enter cell cycle arrest, potentially due to a higher concentration of metabolites in the spheroid centre (Weiswald et al., 2015; Masuda et al., 2016). Finally, during phase 3 the spheroid develops a necrotic core. Eventually, the loss of cells within the spheroid balances growth at the spheroid periphery, stalling net overall growth.

Whether spheroids reach the size required for necrosis to develop relates to experimental design choices such as the experimental duration and initial seeding density, among many other factors. Our hypothesis is that, provided the availability of nutrients is maintained in the cell culture, the structure of a spheroid is eventually a function of spheroid size, independent of the initial seeding density. This presents us with a technical challenge and a biological opportunity for protocol refinement. For example, we find that the initial aggregation of cells into spheroids occurs over several days (Spoerri et al., 2017), a timescale similar to that of cell proliferation. Therefore, the growth of spheroids over a short experimental duration may be significantly influenced by differences in initial seeding density, potentially confounding differences due to variations in cell behaviour between experimental conditions and limiting the reproducibility of experiments. Our analysis of late-time spheroid structure circumvents this by studying structure as function of overall size instead of time. The primary benefit of this approach is that inferences are insensitive to variations in the initial seeding density.

We take a likelihood-based approach to estimating parameters (Lehmann et al., 1998) employ profile likelihood analysis to produce approximate confidence intervals (Raue et al., 2009; Pawitan, 2013; Browning et al., 2021a) and develop a likelihood-ratio-based hypothesis test to assess consistency in results between seeding densities. Firstly, we work solely with a statistical model that describes the average sizes of the spheroid, inhibited region and necrotic core at each observation time. Secondly, we apply a simple mechanistic model that describes spheroid progression due to a balance between growth at the spheroid periphery and mass loss due to necrosis in the spheroid centre. Following the seminal work of Greenspan, 1972, we assume that nutrients and wastes from living cells are at diffusive equilibrium, leading to a functional relationship between spheroid size and inner structure. Comparing model predictions to experimental observations allows us to assess whether the underlying assumptions of the Greenspan model are appropriate, providing valuable information for model refinement. As we are primarily interested in spheroid structure and model validation, we focus our analysis on comparing the structure at different observation times and seeding densities rather than a more typical approach that calibrates the mathematical to all data simultaneously (Murphy et al., 2021).

We are motivated to work with a simple mathematical model instead of a more complex (and potentially more biologically realistic) alternative (Ward and King, 1997; Ward and King, 1999; Roose et al., 2007; Byrne, 2010; Bull et al., 2020) for two reasons. Firstly, complex models are often highly parameterised (Gutenkunst et al., 2007; Gábor and Banga, 2015; Raman et al., 2017). Given the practical difficulties in extracting detailed measurements from spheroids, we do not expect to be able to reliably estimate parameters in many complex models; that is, we expect parameters to be practically non-identifiable (Raue et al., 2009). Working with a simple model avoids over-parameterisation allowing for a better comparison between experimental conditions. Secondly, Greenspan’s model encapsulates our central hypothesis that spheroid structure is purely a function of spheroid size, and captures the key features of spheroid growth seen in the experimental data with a low-dimensional, interpretable, parameter space.

Materials and methods

Experimental methods

Request a detailed protocol

The human melanoma cell lines WM793b (Herlyn et al., 1985) and WM983b (Herlyn, 1990) were genotypically characterised (Hoek et al., 2006; Smalley et al., 2007a; Smalley et al., 2007b), grown as described in Spoerri et al., 2017 supplemented with 1% penicillin-streptomycin (ThermoFisher, Massachusetts, United States), and authenticated by short tandem repeat fingerprinting (QIMR Berghofer Medical Research Institute, Herston, Australia). All cell lines were transduced with fluorescent ubiquitination-based cell cycle indicator (FUCCI) constructs as described in Haass et al., 2014; Spoerri et al., 2017. Wells within a flat-bottomed 96-well plate were prepared with 50 µL non-adherent 1.5% agarose to prevent cell-to-substrate attachment and promote the formation of a single centrally located spheroid (Spoerri et al., 2021). Cells were seeded into each well at a density of approximately 2500, 5000, and 10,000 cells in 200 µL of medium. A medium change was performed every 2–4 days.

Spheroids were harvested and fixed with 4% paraformaldehyde at day 3, 4, 5, 7, 10, 12, 14, 16, 18, 21, and 24; mounted in 2% low melting agarose; placed in a refractive-index-matched clearing solution Spoerri et al., 2021; and imaged using fluorescent confocal microscopy to obtain high-resolution images at the equator of each spheroid (Olympus FV3000, Olympus, Tokyo, Japan). To minimise variability due to the vertical position of each image, spheroids are fixed in place using an agarose gel, and equatorial images are defined as the cross-section with the largest cross-sectional area. To obtain the result in Figure 1i, we selectively stain spheroids with DRAQ7 (ThermoFisher, Massachusetts, United States), which indicates necrosis (Kienzle et al., 2017; Spoerri et al., 2021). Staining, fixation, and microscopy are repeated to obtain at least 20 WM983b spheroids at day 18 (spheroids initially seeded with 5000 and 10,000 cells) and day 21 (spheroids seeded with 2500 cells); and at least 10 spheroids for all other conditions. Data are then randomly subsampled to obtain exactly 10 and 20 spheroids for each initial condition and observation day where possible. Time-lapse phase-contrast and fluorescent channel images are obtained at 6 hr intervals for up to 24 spheroids for each initial condition using an Incucyte S3 (Sartorius, Goettingen, Germany).

Data processing

Request a detailed protocol

We apply a semi-automated data processing algorithm to summarise experimental images with three measurements (Figure 1h; Browning and Murphy, 2021b). Firstly, we calculate the outer radius, $R$ , based on a sphere with the same cross-sectional area as the image obtained. Secondly, the radius of the inhibited region, $R_{i}$ . We calculate the radius of the inhibited region by determining the average distance from the spheroid periphery where the signal from mAG (FUCCI green), which indicates cells in gap 2, falls below a threshold value, taken to be 20% of the maximum area-averaged green signal. We find this choice leads to accurate results (Figure 2). Finally, the radius of the necrotic core, $R_{n}$ , which is identified using texture recognition (stdfilt, Mathworks, 2021). The regions identified using the algorithm are shown in Figure 2. Full details of the image processing algorithms are available in Browning and Murphy, 2021b and additional images are available as supplementary material.

Figure 2

Download asset Open asset

Late-time progression of WM983b spheroids, randomly sampled from the 10 spheroids imaged from each condition (additional images in Supplementary file 2).

Overlaid are the three boundaries identified by the image processing algorithm: the entire spheroid, the inhibited region and the necrotic region. Each image shows a 800 × 800 µm field of view. Colouring indicates cell nuclei positive for mKO2 (magenta), which indicates cells in gap 1; and cell nuclei positive for mAG (green), which indicates cells in gap 2.

Mathematical model

Request a detailed protocol

Following Greenspan, 1972, we make two minimal assumptions regarding growth inhibition and necrosis (Figure 1h). Firstly, that growth inhibition, or cell cycle arrest, is a result of a chemical inhibitor that originates from the metabolic waste of living cells (Laurent et al., 2013). This inhibitor is produced by living cells at rate $β_{prod} [m o l / d]$ and diffuses with diffusivity $β_{diff} [μ m^{2} / d]$ . At the outer boundary of the spheroid, we assume that the concentration of inhibitor is zero. Cells enter arrest in regions where the inhibitor concentration is greater than $β_{crit} [m o l / μ m^{3}]$ . Secondly, cycling cells require nutrients that are plentifully available in the surrounding medium at concentration $ω_{\infty} [m o l / μ m^{3}]$ . The nutrient is consumed by cycling cells at a constant rate $ω_{cons} [m o l / d]$ and diffuses with diffusivity $ω_{diff} [μ m^{2} / d]$ . Cells die in regions where the nutrient concentration is less than $ω_{crit} [m o l / μ m^{3}]$ .

In regions where the nutrient concentration is sufficiently high and the inhibitor concentration sufficiently low, we assume that cells proliferate exponentially at the per-volume rate $s [/ d]$ . Furthermore, we assume that cell debris is lost from the necrotic core at the per-volume rate $λ [/ d]$ .

It is convenient to define two non-dimensional parameters

Q^{2} = \frac{ω_{cons}}{ω_{diff} (ω_{\infty} - ω_{crit})} \times \frac{β_{crit} β_{diff}}{β_{prod}} < 1,

and

γ = \frac{λ}{s} > 0.

The parameter $Q$ quantifies the balance between nutrient and inhibitor concentration and $γ$ quantifies the balance between cell growth and the loss due to necrosis. The restriction $Q < 1$ arises since we observe an inhibited region form before the necrotic region (Greenspan, 1972). Since the resultant equations depend only on $Q$ and $γ$ , the constituents of $Q$ , namely $β_{prod}$ , $β_{diff}$ , $β_{crit}$ , $ω_{cons}$ , $ω_{diff}$ , $ω_{\infty}$ , and $ω_{crit}$ , cannot be uniquely identified unless prior knowledge from other experiments is considered (Murphy et al., 2017), perhaps in a Bayesian framework (Browning et al., 2019). In contrast, the constituents of $γ$ , namely $λ$ and $s$ , can be identified if information relating to the per-volume cell proliferation rate $s$ is available, perhaps from phase one spheroid growth data.

We take the standard approach and model the spheroid as a single spherical mass (Greenspan, 1972; Araujo and McElwain, 2004). We denote by $R$ the radius of the spheroid, $ϕ = R_{i} / R$ the relative radius of the inhibited region, and $η = R_{n} / R$ the relative radius of the necrotic core (Figure 1h). We note that $R > 0$ and $0 \leq η \leq ϕ < 1$ . Noting that nutrient and inhibitor diffusion occurs much faster than cell proliferation, we assume that the chemical species are in diffusive equilibrium, leading to

\frac{d R}{d t} = \underset{G r o w t h i n c y c l i n g r e g i o n}{\underset{⏟}{\frac{s}{3} (1 - ϕ^{3}) R}} - \underset{M a s s l o s s f r o m n e c r o t i c c o r e}{\underset{⏟}{\frac{s}{3} γ η^{3} R .}}

A distinguishing feature of Greenspan’s model is that the inner structure of the spheroid, quantified by $(ϕ, η)$ , is determined solely by the spheroid radius, and not by time. We denote

0 = f_{s} (ϕ, η; R, Q, R_{c}),

as a function describing this relationship, and refer to the relationship between the spheroid radius, $R$ , and the inner structure, $(ϕ, η)$ , as the structural model. Here, we define $R_{c}$ as the radius at which necrosis first occurs. For $R < R_{c}$ , nutrient is available throughout the spheroid above the critical concentration $ω_{crit}$ .

During phases 1 and 2, there is no necrotic core ( $η = 0$ ) and the solution to Equation 4 is given by

ϕ^{2} = max (0, 1 - \frac{Q^{2} R_{c}^{2}}{R^{2}}), R < R_{c} .

During phase 3, $R > R_{c}$ and $f_{s}$ is given by

f_{s} (ϕ, η; R, Q, R_{c}) = (\begin{matrix} 2 R^{2} η^{3} - 3 R^{2} η^{2} + R^{2} - R_{c}^{2}, \\ R^{2} ϕ^{3} + (Q^{2} R_{c}^{2} - R^{2} (1 + 2 η^{3})) ϕ + 2 η^{3} R^{2} \end{matrix}) .

To investigate the limiting structure of spheroids, we consider the solution to the mathematical model where the outer radius is no longer increasing: the dynamics have reached a steady-state. Experimental observations suggest that this occurs during phase 3. We denote $\bar{R} = lim_{t \to \infty} R (t)$ the limiting radius and $(\bar{ϕ}, \bar{η})$ the associated limiting structure. The steady-state model is the solution of

{\begin{aligned} 0 & = 1 - {\bar{ϕ}}^{3} - γ {\bar{η}}^{3}, \\ 0 & = f_{s} (\bar{ϕ}, \bar{η}; \bar{R}, Q, R_{c}), \end{aligned}

subject to $R > R_{c}$ . By defining $ρ = \bar{η} / \bar{ϕ} \in (0, 1)$ , we find a semi-analytical solution to the steady-state model (Appendix 1).

The behaviour in the steady-state model is characterised by three parameters, $θ = (Q, R_{c}, γ)$ . We denote the solution to Equation 7 (i.e. the steady-state model) as

m (θ) : (Q, R_{c}, γ) \to (\bar{R}, \bar{ϕ}, \bar{η}) .

Equation 8 can be thought of as a map from the parameter space to the limiting structure of the spheroid. This demonstrates that the parameters are identifiable only when all three variables, $(\bar{R}, \bar{ϕ}, \bar{η})$ , are observed, since the two-dimensional observation space $(R, η)$ cannot uniquely map to the entire three-dimensional parameter space $(Q, R_{c}, γ)$ . As a consequence, the model parameters cannot be uniquely identified from steady-state information unless phase 3 information that includes measurements of the inhibited region—using FUCCI or another marker of cell cycle arrest—is considered alongside measurements of necrotic core and overall spheroid size.

Statistical model

Request a detailed protocol

While the mathematical model is deterministic, experimental observations of spheroid structure can be highly variable. To account for this, we take the standard approach and assume that the mathematical model describes the expected behaviour and experimental observations are multivariate normally distributed (Lehmann et al., 1998). Aside from accounting for biological variability, the observation process captures variability introduced during imaging and image processing.

Denoting $x_{i} = (R_{i}, ϕ_{i}, η_{i})$ as experimental observation i of the spheroid size and structure, we assume that

x_{i} \sim f (x; μ, Σ) = N (μ, Σ),

where $μ = (R, ϕ, η)$ is the mean of each component of $x$ , $N (μ, Σ)$ denotes a multivariate normal distribution with mean μ and covariance $Σ$ . To account for increased variability at later time points (Figure 1a–b), we estimate $Σ$ as the sample covariance associated with experimental observations of $x_{i}$ at each time, $t$ . For steady-state analysis, we calculate the covariance using the pooled sample covariance from all seeding densities.

We refer to Equation 9 as the statistical model. To connect experimental observations to the mathematical model, we substitute $μ = m (θ)$ in Equation 9.

Inference

We take a likelihood-based approach to parameter inference and sensitivity analysis. Given a set of observations $X = {x_{i}}_{i = 1}^{n}$ , the log-likelihood function is

ℓ (θ; X) = \sum_{i} \log f (x_{i}; m (θ), Σ),

where $f (x; μ, Σ)$ is the multivariate normal probability density function (Equation 9). Although we take a purely likelihood-based approach to inference, we note that our implementation is equivalent to a Bayesian approach where uniform priors encode existing knowledge about parameters, a common choice (Hines et al., 2014; Simpson et al., 2020).

We apply maximum likelihood estimation to obtain point estimates of the parameters for a given set of experimental observations. The maximum likelihood estimate (MLE) is given by

\hat{θ} = \underset{θ}{a r g m a x} ℓ (θ; X) .

We solve Equation 11 numerically to within machine precision using a local optimisation routine (Powell, 2009; Johnson, 2021). In Figure 3, we show point estimates obtained for a bivariate problem using maximum likelihood estimation.

Figure 3

Download asset Open asset

We calculate approximate confidence intervals (CI) using profile likelihood and confidence regions (CR) using contours of the normalised likelihood function.

Results demonstrate estimates of $Q$ and $R_{c}$ using the structural model, Equation 6, and data from WM983b spheroids at day 14 initiated using 5000 cells. Point estimates are calculated using the maximum likelihood estimate (white marker). The boundaries of regions are defined as contours of the log-likelihood function. Univariate confidence intervals are constructed by profiling the log-likelihood and using a threshold of approximately −1.92 for a 95% confidence interval.

Confidence regions and hypothesis tests

Request a detailed protocol

We take a log-likelihood based approach to compute confidence regions and marginal univariate confidence intervals for model parameters (Pawitan, 2013). In a large sample limit, Wilks’ Theorem provides a limiting distribution for the log-likelihood ratio statistic, such that

2 [ℓ (\hat{θ}) - ℓ (θ)] \sim χ^{2} (ν)

where $ν = d i m (θ)$ and $χ^{2} (ν)$ Is the $χ^{2}$ Distribution with $ν$ degrees of freedom. Therefore, an approximate $α$ level confidence region is given by

θ : ℓ (θ) \geq ℓ (\hat{θ}) - \frac{Δ_{ν, α}}{2},

where $Δ_{ν, α}$ is the $α$ level quantile of the $χ^{2} (ν)$ distribution.

Hypothesis tests

Request a detailed protocol

To compare parameters between initial conditions, we perform likelihood-ratio-based hypothesis test based on the distribution provided in Equation 13 (Pawitan, 2013). We denote by ${\hat{θ}}_{*}$ the MLE computed using data from all initial seeding densities, $X_{*}$ , simultaneously. Similarly, to compare parameter estimates from spheroids initially seeded with 2500 and 5000 cells, we denote by ${\hat{θ}}_{N}$ the MLE using a subset of data from spheroids seeded using $N \in {2500, 5000}$ cells. The test statistic is given by

T = 2 (- ℓ ({\hat{θ}}_{*}) + \sum_{N} ℓ ({\hat{θ}}_{N})) \sim χ^{2} (ν)

where $ν$ is number of additional parameters in the case where a different parameter combination is used to describe each initial condition. An approximate p-value is therefore given by $1 - F_{χ^{2} (ν)} (T)$ , where $F_{χ^{2} (ν)}$ is the cumulative distribution function for the $χ^{2} (ν)$ distribution.

Marginal confidence intervals

Request a detailed protocol

The profile likelihood method (Raue et al., 2009; Boiger, 2016) allows for the construction of univariate confidence interval of each parameter. Firstly, we partition the parameter space such that $θ = (ψ, λ)$ where $ψ$ is the parameter of interest and $λ$ is a vector containing the remaining parameters. Taking the supremum of the log-likelihood function over $λ$ and normalising using the MLE gives the normalised profile log-likelihood

{\hat{ℓ}}_{p} (ψ; X) = sup_{λ} ℓ (ψ, λ; X) - ℓ (\hat{θ}; X), {\hat{ℓ}}_{p} \leq 0.

An approximate 95% confidence interval is given by Equation 13 as the region where ${\hat{ℓ}}_{p} (ψ; X) \geq - Δ_{1, 0.95} / 2 \approx - 1.92$ (Pawitan, 2013). We compute the profile log-likelihood numerically using a local optimisation routine (Powell, 2009) with either the MLE, or the nearest profiled point (Boiger, 2016) as an initial guess. In Figure 3, we show profile likelihoods for a bivariate problem.

Confidence regions

Request a detailed protocol

We construct two-dimensional confidence regions using Equation 13 (we construct three-dimensional confidence regions using a sequence of two-dimensional slices). First, we find a point on the boundary of the region, denoted $θ_{0}$ such that $ℓ (θ_{0}) = ℓ (\hat{θ}) - Δ_{ν, α} / 2$ , using bisection to machine precision. Next, we integrate along the likelihood annihilating field; that is, we move in a direction perpendicular to the gradient of the likelihood to obtain a set of points on the level set $ℓ (θ) = ℓ (θ_{0})$ , given by

\begin{array}{l} \frac{d θ}{d t} = (\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}) \nabla_{θ} ℓ (θ), θ (0) = θ_{0} . \end{array}

This calculation is demonstrated for a bivariate problem in Figure 3.

The gradient for the statistical model, $\nabla_{μ} ℓ (μ)$ , can be calculated to within machine precision using automatic differentiation (Revels et al., 2016). For the mathematical model, we apply the identity

\nabla_{θ} ℓ (θ) = J_{m} (θ) \nabla_{μ} ℓ (m (θ)),

where $J_{m} (θ)$ is calculated analytically (Appendix 1).

Results

To assess the limiting structure of spheroids and the effect of initial seeding density, we analyse confocal sections of a large number of spheroids across three seeding densities using the WM983b cell line. We show a subset of these images in Figure 2 and summarise images with three concentric annular measurements: the spheroid radius, $R$ ; the relative radius of the inhibited region, $ϕ$ ; and the relative radius of the necrotic core, $η$ (Figure 1h). In addition to spheroids from different initial conditions tending towards a similar overall size (as seen from time-lapse data in Figure 1a–f), these results show that spheroids develop similar structures by day 21.

First, we fit the statistical model to the experimental data by estimating the mean of each measurement, denoted $μ = (R, ϕ, η)$ . We obtain a maximum likelihood estimate and an approximate 95% confidence interval for each initial condition at observation days 12–21 (Figure 4a–c). On average, spheroids of all seeding densities increase in size from day 12 to day 18. In agreement with earlier observations from time-lapse data in Figure 1e–h, we see that spheroids initiated at different seeding densities tend toward similar limiting sizes. Between days 18 and 21, spheroids seeded with 5000 and 10,000 cells decrease in average size, potentially indicating a period of decay after a limiting size is reached.

Figure 4

Download asset Open asset

Estimates of parameters using the structural model with data from various time points.

In (**a–c**), parameters are the mean of each observation: $(R, ϕ, η)$ . In (**d–e**), parameters are those in the structural model: $(R, Q, R_{c})$ . In (f), estimates of $γ$ are obtained by calibrating observations to the steady-state model. As estimates $Q$ and $R_{c}$ can be derived from the structural model (Equation 6), which applies at any time during phase 3, we expect to see similar parameter estimates across observation times. As estimates of $γ$ can only be obtained from the steady-state model, which assumes the outer radius is no longer increasing, we do not expect to see similar parameter estimates across observation times. Bars indicate an approximate 95% confidence interval.

Figure 4b and c show estimates relating to the sizes of the inhibited region, $ϕ$ , and necrotic core, $η$ . We see remarkable consistency in $ϕ$ across seeding densities, tending toward a value of 90% in all cases: this corresponds to an actively cycling region with volume approximately 27% of the total spheroid volume. The necrotic core increases significantly in size from days 12 to 21, and late time estimates of $η$ are quantitively similar between seeding densities.

Next, we calibrate the mathematical model to identify any mechanistic differences between seeding densities. Parameters $Q$ and $R_{c}$ can be estimated using the structural model (Equation 6) at any time point. To estimate $γ$ we must invoke the steady-state model (Equation 7), which assumes that the overall growth of the spheroid has ceased. Therefore, we expect to see consistency in estimates of $Q$ and $R_{c}$ between observation days but do not expect the same for estimates of $γ$ .

Results in Figure 4d show remarkable consistency in estimates of $Q$ across seeding densities until day 18, suggesting that the balance between nutrient availability and waste concentration (Equation 1) is maintained throughout the experiment and is similar between seeding densities. Between days 18 and 21, estimates of $Q$ for spheroids initially seeded with 5000 and 10,000 cells increase significantly, suggesting a behavioural change during this time; we attribute this to a final period of decay. Estimates of $R_{c}$ do not show the consistency between observation days we might expect if $f_{s}$ (Equation 6) holds for the experimental data. Rather, estimates of $R_{c}$ decrease between days 12 and 21, indicating $f_{s}$ may be misspecified. Results in Figure 4f show that estimates of $γ$ decrease with time to a similar value for all seeding densities. We interpret this asymptotic decrease as an indication that spheroids approach a limiting structure since estimates of $γ$ are strictly only valid when growth has ceased. Closer inspection of results in Figure 4f show a delay in estimates of $γ$ between spheroids seeded with 2500 cells and the other seeding densities. Whereas the larger spheroids reach a limiting size by day 18, the smaller spheroids are still growing. It is not until day 21 that estimates of $γ$ are comparable across all seeding densities.

Next, we analyse the limiting structure of spheroids across each initial seeding density. As spheroids initially seeded with 5000 and 10,000 cells decrease in average size from day 18 to day 21, we compare day 18 data from these high densities to day 21 data from spheroids initially seeded with 2500 cells. Results in Figure 5 show profile log-likelihoods for each parameter in the mathematical model. In Figure 5, we show 3D confidence regions for parameters in the statistical and mathematical models, respectively. We see that both profile log-likelihoods and 3D confidence regions overlap, indicating that parameter estimates are consistent between seeding densities.

Figure 5

Download asset Open asset

Comparison of WM983b spheroids between each initial seeding density at day 18 (spheroids seeded with 5000 or 10000 cells) and day 21 (2500).

(**a–c**) Profile likelihoods for each parameter, which are used to compute approximate confidence intervals (Table 1). (d) 95% confidence region for the full parameter space. 95% confidence regions for (d) the mean of each observation at steady state $(\bar{R}, \bar{φ}, \bar{η})$ and (e) the model parameters $(Q, R_{c}, γ)$ .

To compare quantitatively parameter estimates between seeding densities, we tabulate maximum likelihood estimates, approximate 95% confidence intervals, and results of a likelihood-ratio-based hypothesis test for both models in Table 1. The late-time sizes of spheroids initiated with 5,000 and 10,000 cells are statistically consistent (p=0.62), as is their structure (p=0.69). We find evidence to suggest that spheroids seeded with 2500 cells, even at day 21, are smaller (p=0.04); however, the overall size and structure of the spheroids seeded with 2500 and 5000 cells are statistically consistent (p=0.20). We find no significant differences in model parameters between seeding densities and note that the conclusion of overall statistical consistency between seeding densities is identical for the mathematical model.

Table 1

Parameter estimates and approximate confidence intervals for each initial conditions.

Also shown are p-values for likelihood-ratio-based hypothesis tests for parameter equivalence between seeding densities.

Parameter	θ₂₅₀₀	θ₅₀₀₀	θ₁₀₀₀₀	p _2500,5000	p _5000,10000
$R$	340.0 (331.0, 349.0)	353.0 (344.0, 361.0)	356.0 (347.0, 365.0)	0.0420	0.617
$ϕ$	0.899 (0.889, 0.908)	0.895 (0.886, 0.905)	0.901 (0.891, 0.911)	0.617	0.406
$η$	0.719 (0.674, 0.764)	0.716 (0.671, 0.761)	0.742 (0.696, 0.788)	0.940	0.438
μ				0.202	0.687
$Q$	0.75 (0.696, 0.811)	0.758 (0.704, 0.818)	0.771 (0.711, 0.838)	0.854	0.767
$R_{c}$	149.0 (127.0, 171.0)	156.0 (133.0, 178.0)	145.0 (121.0, 168.0)	0.672	0.503
$γ$	0.737 (0.598, 0.916)	0.768 (0.624, 0.953)	0.657 (0.532, 0.816)	0.792	0.308
$θ$				0.202	0.687

Next, we investigate the relationship between spheroid structure and spheroid size from day 3 to day 21 (Figure 6). We again see evidence of a period of eventual decay that occurs after a limiting size has been reached in our experiments. To validate the structural relationship suggested by Greenspan’s model, we plot the solution to the structural model (Equation 6) using parameters estimated using the steady-state model (Table 1). The overall trend throughout all three phases of growth in the mathematical model—made only using information from days 18 and 21—is remarkably consistent with experimental measurements Figure 6. We find an explanation for the inconsistent estimates of $R_{c}$ observed in Figure 4e. During phase 3, the mathematical model predicts a non-linear relationship between $R$ , $ϕ$ and $η$ (Equation 6). In contrast, the trend in the data is close to linear. We confirm this in Figure 6 by calibrating a linear model of the form

(R (τ), ϕ (τ), η (τ)) = (R_{c}, ϕ_{c}, 0) + τ \hat{q},

Figure 6

Download asset Open asset

Data from days 3 to 21 (WM983b) and days 4 to 24 (WM793b) for all initial conditions.

Solid curves in (a) show the solution to the mathematical model (Equation 6) using the maximum likelihood estimate calculated using the steady-state data (Table 1). Solid curves in (c) show the solution to the mathematical model (Equation 6) using using the maximum likelihood estimate calculated using day 24 data. In (b) and (d), we fit a linear model to phase 3 data (indicated by coloured markers). The p value corresponds to a hypothesis test where the linear model parameters are the equivalent for all initial conditions. Shown in black is the best fit linear model.

to phase 3 data using a total least squares approach that accounts for uncertainty in the independent variable $τ$ (Appendix 2). Here, $τ = 0$ at the start of phase 3. Performing an approximate likelihood-ratio-based hypothesis test confirms that the behaviour in spheroids of all initial conditions is statistically consistent (p=0.56). That is, the spheroid structure where necrosis first occurs (at $τ = 0$ ), $(R_{c}, ϕ_{c}, 0)$ , and the direction in which it develops, $\hat{q}$ , do not appear to depend on the initial seeding density.

In Figure 6, we perform a similar analysis on spheroids grown from WM793b cells. Whereas WM983b spheroids approach a limiting size by the conclusion of the experiment (Figure 1a), spheroids grown from the WM793b do not (Figure 1b). Results in Appendix 3 examine parameter estimates from the mathematical and statistical models through time for the WM793b spheroids, demonstrating that the outer radius increases monotonically until day 24 for all initial conditions. These results also suggest consistency in estimates of $Q$ across observation days and seeding densities. Performing a likelihood-ratio-based hypothesis test indicates that phase 3 is independent of the initial seeding density (p=0.36).

Discussion

Time-lapse measurements of WM983b spheroids over a 21-day experiment show a cessation in overall growth as the spheroids reach a limiting size. Consistent with largely untested predictions of classical mathematical models (Greenspan, 1972; Ward and King, 1999; Araujo and McElwain, 2004), these limiting sizes appear to be independent of the initial seeding density. Motivated by these observations, we develop a quantitative framework to study spheroid structure as a function of overall size. We aim to answer two fundamental questions: Do these spheroids have a limiting structure? Is the late-time behaviour independent of the initial seeding density?

We find compelling evidence that WM983b spheroids have a limiting structure that is independent of the initial seeding density. This assumption is routinely invoked in mathematical models of tumour structure but is yet to be experimentally verified. Given that we observe spheroids to eventually reduce in size, we compare structural measurements at days when the average outer radius for each initial seeding density is largest. First, we establish that spheroids seeded with 5000 and 10,000 cells have similar limiting sizes (353 µm and 356 µm, respectively; p=0.62) and that spheroids seeded with 2500 cells are slightly smaller at late time (340 µm). This result highlights one of the challenges in determining the limiting structure of spheroids: it is unclear whether there is a difference or whether the smaller spheroids would continue to grow given additional time. Despite this discrepancy, we find a statistically consistent limiting structure, with a necrotic core of 73% of the outer radius and an inhibited region of 90% of the outer radius, indicating a proliferative periphery approximately 35 µm (two to three cell diameters) thick.

By examining spheroid structure throughout the entire experiment (Figure 6), we establish a relationship between spheroid structure and size that is independent of initial seeding density. This result is significant as it suggests that variability in size and structure may be primarily attributable to time. For example, spheroids that are smaller than average on a given observation day may have been seeded at a lower density. Statistical techniques, such as ODE-constrained mixed effects models, can be applied to elucidate sources of intrinsic variability, such as variability in the initial seeding density (Wang et al., 2012; Hasenauer et al., 2014). It is common in the literature to compare spheroids with and without a putative drug after a fixed number of days (Friedrich et al., 2009). However, our analysis suggests that comparing the structure of spheroids of a fixed size may be more insightful; this approach obviates variability due to initial seeding density, increasing the sensitivity of statistical tests to small effects. A corollary is that since inferences relating to spheroid structure are independent of spheroid size, experiments can be initiated with a larger number of cells to decrease the time until spheroids reach phase 3.

Given our observations of WM983b spheroids across seeding densities, an apparent conclusion of our analysis is that statistically consistent phase 3 behaviour implies a statistically consistent limiting structure. If true, this suggests that an experimentalist only has to investigate phase 3 behaviour to reach a conclusion relating to the limiting structure. Analysis of both the mathematical model and experimental results for WM793b spheroids indicate that this is not the case. In the mathematical model, $f_{s}$ (Equation 6) characterises the structure solely in terms of parameters $Q$ and $R_{c}$ , whereas $γ$ —which relates to the ratio of cell proliferation and loss due to necrosis (Equation 2)—determines the steady-state. We see this for WM793b spheroids, as phase 3 behaviour is independent of the initial seeding density (Figure 6), but time-lapse data of the overall growth (Figure 1b) gives no indication that spheroids of different densities will tend toward the same limiting size.

As $f_{s}$ determines the relationship between spheroid size and structure at any time point, we expect estimates of $Q$ and $R_{c}$ to be similar when calibrated to data from different days. This is the case for estimates of $Q$ (Figure 4d), but estimates of $R_{c}$ decrease with time (Figure 4e). While the mathematical model captures the same overall behaviour observed in the experiments, it is evident from the discrepancy observed during phase 3 (Figure 6) that $f_{s}$ is misspecified. Our assumptions of nutrient and waste at diffusive equilibrium and a hard threshold for growth inhibition and necrosis give rise to $f_{s}$ that is cubic in $ϕ$ and $η$ . Since the empirical relationship for the cell lines we investigate is approximately linear, the model underestimates the radius at which phase 3 begins, $R_{c}$ . At the loss of mechanistic insight, one approach to rectify this discrepancy is to construct a purely phenomenological relationship where $f_{s}$ is piecewise linear. A second approach is to revisit fundamental modelling assumptions to develop a mechanistic description of the relationship between spheroid structure and overall size that is consistent with our experimental observations for these cell lines.

Our observations for WM983b and WM793b melanoma cell lines do not preclude a form of $f_{s}$ that is cubic for other cell lines or experimental conditions. In our framework, the behaviour of spheroids is characterised by the empirical relationship between spheroid size and structure. Therefore, despite misspecification in parameter estimates of $R_{c}$ , we can compare spheroids grown with WM793b and WM983b cell lines by comparing the structural relationship observed in the experimental data (Figure 6). In this case, we observe that radius at which the necrotic core develops is much smaller in WM983b spheroids than for WM793b spheroids. While we cannot elucidate the biological factors that lead to this difference from our analysis, we postulate that differences in the diffusion or consumption of nutrients by cells of each cell line may contribute.

We have restricted our analysis of spheroid structure to three measurements that quantify the sizes of the spheroid, inhibited region and necrotic core. While the spheroid and necrotic core sizes are objective measurements, the boundary of the inhibited region is not. Our approach is to identify the distance from the spheroid periphery where the density of cells in gap 2 falls below 20% of the maximum. We find this semi-automated approach produces excellent results and enables high-throughput analysis of hundreds of spheroids; however, it does not take advantage of all the information available in the experimental images. Mathematical models that explicitly include variation in cell density through space (Ward and King, 1999; Jin et al., 2021) may be appropriate, however are typically heavily parameterised, limiting the insight obtainable from typical experimental data. The mass-balance model coupled to a model describing the relationship between spheroid size and structure avoids these issues and, despite model simplicity, we are still able to gain useful biological insight.

Conclusion

Reproducibility and size uniformity are paramount in practical applications of spheroid models. Yet, the effect of intentional or unintentional variability in spheroid size on the inner structure that develops is not well understood. We present a quantitative framework to analyse spheroid structure as a function of overall size, finding that the outer radius characterises the inner structure of spheroids grown from two melanoma cell lines. Further, we find that the initial seeding density has little effect on the structure that develops. These results attest to the reproducibility of spheroids as an in vitro research tool. While we analyse data from two melanoma cell lines, our focus on commonly reported spheroid measurements allows our framework to be applied more generally to a other cell lines and culture conditions. It is routine to compare spheroid size and structure of spheroids at a pre-determined time, our results suggest a refined protocol that compares the structure of spheroids at a pre-determined overall size.

Given the prominence of spheroids in experimental research, there is a surprising scarcity of experimentally validated mathematical models that can be applied to interpret data from these experiments. We find that one of the earliest and simplest models of tumour progression—the seminal model of Greenspan, 1972—can give valuable insights with a parameter space that matches the level of detail available from spheroid structure data. Given that we establish an empirical relationship between spheroid size and structure independent of both time and the initial spheroid size, we suggest future theoretical work to identify mechanisms that give rise to this relationship, perhaps through equation learning (Lagergren et al., 2020). To aid in validating theoretical models of spheroid growth, we make our highly detailed experimental data freely available.

Data availability

Code, data, and interactive figures are available as a Julia module on GitHub at github.com/ap-browning/Spheroids (Browning, 2021c; copy archived at swh:1:rev:27f9e32bb702cb56a62bacaae1e49746a3c4342d). Code used to process the experimental images is available on Zenodo (Browning and Murphy, 2021b).

Appendix 1

Steady-state model solution

The steady-state, denoted $(\bar{R}, \bar{ϕ}, \bar{η})$ is given by setting $d R / d t = 0$ (Equation 3 in the main text) yielding the non-linear system of equations

\begin{array}{ll} 0 = 1 - {\bar{ϕ}}^{3} - γ {\bar{η}}^{3}, \\ 0 = 2 {\bar{R}}^{2} {\bar{η}}^{3} - 3 {\bar{R}}^{2} {\bar{η}}^{2} + {\bar{R}}^{2} - {\bar{R}}_{c}^{2}, \\ 0 = {\bar{R}}^{2} {\bar{ϕ}}^{3} + (Q^{2} {\bar{R}}_{c}^{2} - {\bar{R}}^{2} (1 + 2 {\bar{η}}^{3})) \bar{ϕ} + 2 {\bar{η}}^{3} {\bar{R}}^{2} . \end{array}

Applying the substitution $ρ = \bar{η} / \bar{ϕ}$ , where 0 ≤ ρ ≤ 1, and algebraic manipulation allows the solution to Equations 19 to be expressed as the root of $f (ρ; Q, γ)$ , where

f (ρ; Q, γ) = \sum_{m = 0}^{12} c_{m} ρ^{m},

and where

\begin{array}{ll} c_{0} = 3 Q^{2} - 3 Q^{4} + Q^{6}, \\ c_{1} = 0, \\ c_{2} = - 9 Q^{2}, \\ c_{3} = 18 Q^{2} - 18 Q^{4} + 6 Q^{6} - 2 γ + 9 Q^{2} γ - 9 Q^{4} γ + 3 Q^{6} γ, \\ c_{4} = 27 Q^{4}, \\ c_{5} = - 36 Q^{2} - 9 Q^{2} γ, \\ c_{6} = 36 Q^{2} - 36 Q^{4} - 15 Q^{6} - 6 γ + 36 Q^{2} γ - 36 Q^{4} γ \\ + 12 Q^{6} γ - 3 γ^{2} + 9 Q^{2} γ^{2} - 9 Q^{4} γ^{2} + 3 Q^{6} γ^{2}, \\ c_{7} = 54 Q^{4} + 27 Q^{4} γ, \\ c_{8} = - 36 Q^{2} - 36 Q^{2} γ, \\ c_{9} = 24 Q^{2} - 24 Q^{4} + 8 Q^{6} + 36 Q^{2} γ - 36 Q^{4} γ - 15 Q^{6} γ - 6 γ^{2} + 18 Q^{2} γ^{2} \\ - 18 Q^{4} γ^{2} + 6 Q^{6} γ^{2} - γ^{3} + 3 Q^{2} γ^{3} - 3 Q^{4} γ^{3} + Q^{6} γ^{3}, \\ c_{10} = 54 Q^{4} γ, \\ c_{11} = - 36 Q^{2} γ, \\ c_{12} = 8 γ . \end{array}

Since $ρ$ is subject to the constraint $0 \leq ρ \leq 1$ , we solve $0 = f (ρ; Q, γ)$ using bisection (Implemented to within machine precision using Roots.jl), which is guaranteed to converge provided there exists only one root in the interval $0 \leq ρ \leq 1$ . In Figure 1a, we demonstrate that in the parameter region of interest ( $0 < Q < 1$ , $γ > 0$ ) there exists only a single solution to Equation 20. We do this by finding all 12 roots of Equation 20 (Implemented by finding the eigenvalues of the characteristic matrix using Polynomials.jl) and counting the number of real roots where ( $0 \leq ρ \leq 1$ ).

The solution to Equation 19 is then given by

\bar{R} = f_{R} (ρ, \bar{ϕ}, θ) = \frac{R_{c}}{(1 - ρ ϕ) \sqrt{1 + 2 ρ ϕ}},

\bar{ϕ} = f_{ϕ} (ρ, θ) = \frac{1}{(1 + γ ρ^{3})^{1 / 3}},

\bar{η} = f_{η} (ρ, \bar{ϕ}) = ρ ϕ,

where $θ = (Q, R_{c}, γ)$ .

In Appendix 1—figure 1b, we compare a numerical solution to the transient model to the semi-analytical solution for the steady state showing an excellent match. All algorithms used to produce the results relating to the mathematical model are available on Github in Module/Greenspan.jl.

Appendix 1—figure 1

Download asset Open asset

Number of solutions of Equation 20.

(a) Number of solutions to Equation 20 subject to the constraint 0 ≤ ρ ≤ 1. Dashed line indicates the region of interest, where γ > 0 and 0 < Q < 1. (b) Comparison between a long-term solution to the transient model and the semi-analytical solution to the steady state, where Q = 0.8, γ = 1, R_c = 150, s = 1 and R₀ = 100.

Jacobian of the steady-state model

In the main document, we denote the solution to Equation 19 as $m (θ)$ . Here, we demonstrate how given a value $(\bar{R}, \bar{ϕ}, \bar{η}) = m (θ)$ , we can obtain an analytical expression for the model Jacobian,

J_{m} (θ) = \frac{\partial m}{\partial θ} .

Given $ρ$ , we can form an analytical expression for Equation 22. Noting that the coefficients of Equation 20 are functions of $θ$ , we consider

\begin{array}{ll} \frac{\partial}{\partial c_{i}} (0) = 0 = \sum_{m = 0}^{12} \frac{\partial}{\partial c_{i}} (c_{m} ρ^{m}) = \frac{\partial}{\partial c_{i}} (c_{i} ρ^{i}) + \sum_{\begin{matrix} m = 0 \\ m \neq i \end{matrix}}^{12} \frac{\partial}{\partial c_{i}} (c_{m} ρ^{m}), \\ = ρ^{i} + c_{i} i ρ^{i - 1} \frac{\partial ρ}{\partial c_{i}} + \sum_{\begin{matrix} m = 0 \\ m \neq i \end{matrix}}^{12} c_{m} m ρ^{m - 1} \frac{\partial ρ}{\partial c_{i}}, \\ = ρ^{i} + \frac{\partial ρ}{\partial c_{i}} \sum_{m = 0}^{12} c_{m} m ρ^{m - 1}, \end{array}

which yields

\frac{\partial ρ}{\partial c_{i}} = \frac{- ρ^{i}}{\sum_{m = 0}^{12} m c_{m} ρ^{m - 1}} = - ρ^{i} {(\frac{\partial f}{\partial ρ})}^{- 1} .

Therefore,

\frac{\partial ρ}{d θ} = \frac{\partial ρ}{\partial c} \frac{\partial c}{\partial θ},

where $c = (c_{0}, c_{1}, . . ., c_{12}); \partial ρ / \partial c = (\partial ρ / \partial c_{0}, . . ., \partial ρ / \partial c_{12})$ and $\partial c / \partial θ$ is the Jacobian of c with respect to θ.

Therefore, we have that

\frac{d \bar{ϕ}}{d θ} = \frac{\partial f_{ϕ}}{\partial \bar{ϕ}} + \frac{\partial f_{ϕ}}{\partial ρ} \frac{d ρ}{d θ},

and it follows that

\frac{d \bar{η}}{d θ} = \frac{\partial f_{η}}{\partial \bar{ϕ}} \frac{\partial \bar{ϕ}}{d θ} + \frac{\partial f_{η}}{\partial ρ} \frac{d ρ}{d θ},

\frac{d \bar{R}}{d θ} = \frac{\partial f_{R}}{\partial \bar{ϕ}} \frac{d \bar{ϕ}}{d θ} + \frac{\partial f_{R}}{\partial ρ} \frac{d ρ}{d θ} + \frac{\partial f_{R}}{\partial θ} .

Therefore, an analytical expression for $J_{m} (θ)$ (Equation 22) is given by

J_{m} (θ) = (\frac{d \bar{R}}{d θ}, \frac{d \bar{ϕ}}{d θ}, \frac{d \bar{η}}{d θ}) .

Appendix 2

Total squares regression

In typical least-squares estimation we fit a model of the form

y_{i} = a + b x_{i} + ε_{y, i},

where $ε_{y, i} \sim N (0, σ_{y})$ is assumed to be a normally distributed error component in $y$ component (Markovsky and Van Huffel, 2007), and $(a, b)$ are model parameters. Least-squares and maximum likelihood estimates $(\hat{a}, \hat{b})$ can then be found by minimising the sum-square error

(\hat{a}, \hat{b}) = \underset{(a, b)}{a r g m i n} \sum_{i} (y_{i} - (a + b x_{i}))^{2} .

We demonstrate this in Appendix 2—figure 1. In typical least squares estimation, we minimise the vertical distance between the data points and the regression line (blue dashed).

In the main document, we fit a linear model to data of the form $(R, ϕ, η)$ , where each component contains an error term. In two-dimensions, this is akin to a model of the form

y_{i} = a + b x_{i} + ε_{y, i} + b ε_{x, i} .

where we have included an additional error term $ε_{x, i} \sim N (0, σ_{x})$ , assumed to be a normally distributed error component in x_i. In this case, the least squares estimate is given by minimising the total perpendicular distance between the data points and the regression line (Appendix 2—figure 1, blue solid) (Markovsky and Van Huffel, 2007).

In the main paper, we fit a linear model of the form

(R (τ), ϕ (τ), η (τ)) = (R_{c}, ϕ_{c}, 0) + τ \hat{q},

parameterised by $R_{c}, ϕ_{c}$ and a unit vector $\hat{q}$ .

If we denote $y_{0} = (R_{c}, ϕ_{c}, 0)$ and $y_{1} = (R_{c}, ϕ_{c}, 0) + \hat{q}$ , then the shortest distance between observation $x_{i} = (R_{i}, ϕ_{i}, η_{i})$ is given by

d (x_{i}; R_{c}, ϕ_{c}, \hat{q}) = \frac{‖ (x_{i} - y_{0}) \times (x_{i} - y_{1}) ‖}{‖ y_{0} - y_{1} ‖},

where $∥ \cdot ∥$ denotes the Frobenius norm, and × denotes the vector cross product.

Therefore, least-squares estimates of the parameters can then be found by minimising the sum-square error

\underset{(R_{c}, ϕ_{c}, \hat{q})}{m i n} \sum_{i} d (X_{i}; R_{c}, ϕ_{c}, \hat{q}) .

Approximating the likelihood

To implement a log-likelihood-ratio based hypothesis test, we must approximate the likelihood at the parameter estimates. To do this, we note that the total square error, denoted $ε_{i}^{2}$ , is of the form

ε_{i}^{2} = c_{1} ε_{x, i}^{2} + c_{2} ε_{y, i}^{2} + c_{3} ε_{z, i}^{2},

where $ε_{x, i}$ , $ε_{y, i}$ , and $ε_{z, i}$ are normally distributed with variances $σ_{x}^{2}$ , $σ_{y}^{2}$ and $σ_{z}^{2}$ , respectively. If $σ_{x}^{2} = σ_{y}^{2} = σ_{z}^{2}$ , $ε_{i}^{2}$ would have an approximate chi-squared distribution by the Welch-Satterthwaite equation (Welch, 1947), a special case of the gamma distribution. Therefore, we approximate the distribution of $ε_{i}^{2}$ by fitting a gamma-distribution to the observed square error when a total squares estimate is fit to the combined data (Figure 1b).

Therefore, the approximate log-likelihood is given by

ℓ (R_{c}, ϕ_{c}, \hat{q}) = \sum_{i} \log f_{Γ} (d^{2} (x_{i}; R_{c}, ϕ_{c}, \hat{q})),

where $f_{Γ} (\cdot)$ is the probability density function of the fitted gamma function.

Appendix 2—figure 1

Download asset Open asset

Fitting experimental data to linear model.

(a) Comparison between typical least-squares error (blue dashed), and total-least-squares error (blue solid). (b) Square error observed in the data and fitted gamma distribution.

Log-likelihood-ratio based test

We denote ${\hat{θ}}_{0} = (R_{c}, ϕ_{c}, \hat{q})$ the maximum likelihood estimate when the data from all initial conditions is pooled, and ${\hat{θ}}_{N} = (R_{c}, ϕ_{c}, \hat{q})$ the estimates from initial condition $N \in {2500, 5000, 10000}$ . As the models must be nested for the likelihood-ratio test, we estimate the noise model, $f_{Γ} (\cdot)$ , using the estimates from the pooled data.

The test-statistic is given by

λ = ℓ ({\hat{θ}}_{2500}) + ℓ ({\hat{θ}}_{5000}) + ℓ ({\hat{θ}}_{10000}) - ℓ ({\hat{θ}}_{0}),

where $λ \sim χ_{ν}^{2}$ , and

ν = \dim ({\hat{θ}}_{2500}) + \dim ({\hat{θ}}_{5000}) + \dim ({\hat{θ}}_{10000}) - \dim ({\hat{θ}}_{0}) = 8.

Our implementation of this test is provided on GitHub in Module/Inference in the function lm_orthogonal_test.

Appendix 3

Results for WM793b

Appendix 3—figure 1

Download asset Open asset

Data availability

Code, data, and interactive figures are available as a Julia module on GitHub (https://github.com/ap-browning/Spheroids copy archived at https://archive.softwareheritage.org/swh:1:rev:27f9e32bb702cb56a62bacaae1e49746a3c4342d). Code used to process the experimental images is available on Zenodo (https://doi.org/10.5281/zenodo.5121093).

The following data sets were generated

1. Browning AP
(2021) Github
ID v.0.6.2. Quantitative analysis of tumour spheroid structure.

https://github.com/ap-browning/Spheroids
1. Browning AP
2. Murphy RJ
(2021) Zenodo
Image processing algorithm to identify structure of tumour spheroids with cell cycle labelling.

https://doi.org/10.5281/zenodo.5121093

References

1. Adam JA
2. Maggelakis SA
(1990) Diffusion regulated growth characteristics of a spherical prevascular carcinoma
Bulletin of Mathematical Biology 52:549–582.

https://doi.org/10.1007/BF02462267
- PubMed
- Google Scholar
1. Alexander S
2. Friedl P
(2012) Cancer invasion and resistance: interconnected processes of disease progression and therapy failure
Trends in Molecular Medicine 18:13–26.

https://doi.org/10.1016/j.molmed.2011.11.003
- PubMed
- Google Scholar
1. Araujo RP
2. McElwain DLS
(2004) A history of the study of solid tumour growth: the contribution of mathematical modelling
Bulletin of Mathematical Biology 66:1039–1091.

https://doi.org/10.1016/j.bulm.2003.11.002
- PubMed
- Google Scholar
1. Beaumont KA
2. Anfosso A
3. Ahmed F
4. Weninger W
5. Haass NK
(2015) Imaging- and Flow Cytometry-based Analysis of Cell Position and the Cell Cycle in 3D Melanoma Spheroids
Journal of Visualized Experiments 1:e53486.

https://doi.org/10.3791/53486
- PubMed
- Google Scholar
1. Boiger R
(2016) Integration based profile likelihood calculation for PDE constrained parameter estimation problems
Inverse Problems 32:125009.

https://doi.org/10.1088/0266-5611/32/12/125009
- Google Scholar
(2019) A Bayesian Sequential Learning Framework to Parameterise Continuum Models of Melanoma Invasion into Human Skin
Bulletin of Mathematical Biology 81:676–698.

https://doi.org/10.1007/s11538-018-0532-1
- PubMed
- Google Scholar
(2021a) Model-based data analysis of tissue growth in thin 3D printed scaffolds
Journal of Theoretical Biology 528:110852.

https://doi.org/10.1016/j.jtbi.2021.110852
- PubMed
- Google Scholar
Software
1. Browning AP
2. Murphy RJ
(2021b) Image processing algorithm to identify structure of tumour spheroids with cell cycle labelling
Zenodo.

https://doi.org/10.5281/zenodo.5121093
Software
1. Browning AP
(2021c) Spheriods, version swh:1:rev:27f9e32bb702cb56a62bacaae1e49746a3c4342d
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:c6768d440d317ed4739b1ac9cae0299c5f86edf1;origin=https://github.com/ap-browning/Spheroids;visit=swh:1:snp:cd949bd0004234322ca6971c60cd14dbdecf1ee1;anchor=swh:1:rev:27f9e32bb702cb56a62bacaae1e49746a3c4342d
1. Bull JA
2. Mech F
3. Quaiser T
4. Waters SL
5. Byrne HM
(2020) Mathematical modelling reveals cellular dynamics within tumour spheroids
PLOS Computational Biology 16:e1007961.

https://doi.org/10.1371/journal.pcbi.1007961
- PubMed
- Google Scholar
1. Byrne HM
2. Chaplain MAJ
(1998) Necrosis and Apoptosis: Distinct Cell Loss Mechanisms in a Mathematical Model of Avascular Tumour Growth
Journal of Theoretical Medicine 1:223–235.

https://doi.org/10.1080/10273669808833021
- Google Scholar
1. Byrne HM
(2010) Dissecting cancer through mathematics: from the cell to the animal model
Nature Reviews. Cancer 10:221–230.

https://doi.org/10.1038/nrc2808
- PubMed
- Google Scholar
1. Cui X
2. Hartanto Y
3. Zhang H
(2017) Advances in multicellular spheroids formation
Journal of the Royal Society, Interface 14:20160877.

https://doi.org/10.1098/rsif.2016.0877
- PubMed
- Google Scholar
1. Eilenberger C
2. Rothbauer M
3. Selinger F
4. Gerhartl A
5. Jordan C
6. Harasek M
7. Schädl B
8. Grillari J
9. Weghuber J
10. Neuhaus W
11. Küpcü S
12. Ertl P
(2021) A Microfluidic Multisize Spheroid Array for Multiparametric Screening of Anticancer Drugs and Blood-Brain Barrier Transport Properties
Advanced Science 8:e2004856.

https://doi.org/10.1002/advs.202004856
- PubMed
- Google Scholar
1. Flegg JA
2. Nataraj N
(2019) Mathematical Modelling and Avascular Tumour Growth
Resonance 24:313–325.

https://doi.org/10.1007/s12045-019-0782-8
- Google Scholar
1. Folkman J
2. Hochberg M
(1973) Self-regulation of growth in three dimensions
The Journal of Experimental Medicine 138:745–753.

https://doi.org/10.1084/jem.138.4.745
- PubMed
- Google Scholar
(2009) Spheroid-based drug screen: considerations and practical approach
Nature Protocols 4:309–324.

https://doi.org/10.1038/nprot.2008.226
- PubMed
- Google Scholar
1. Gábor A
2. Banga JR
(2015) Robust and efficient parameter estimation in dynamic models of biological systems
BMC Systems Biology 9:74.

https://doi.org/10.1186/s12918-015-0219-2
- PubMed
- Google Scholar
(2016) Oxygen Partial Pressure Is a Rate-Limiting Parameter for Cell Proliferation in 3D Spheroids Grown in Physioxic Culture Condition
PLOS ONE 11:e0161239.

https://doi.org/10.1371/journal.pone.0161239
- PubMed
- Google Scholar
1. Greenspan HP
(1972) Models for the Growth of a Solid Tumor by Diffusion
Studies in Applied Mathematics 51:317–340.

https://doi.org/10.1002/sapm1972514317
- Google Scholar
1. Groebe K
2. Mueller-Klieser W
(1996) On the relation between size of necrosis and diameter of tumor spheroids
International Journal of Radiation Oncology, Biology, Physics 34:395–401.

https://doi.org/10.1016/0360-3016(95)02065-9
- PubMed
- Google Scholar
(2007) Universally sloppy parameter sensitivities in systems biology models
PLOS Computational Biology 3:0030189.

https://doi.org/10.1371/journal.pcbi.0030189
- PubMed
- Google Scholar
1. Haass NK
2. Beaumont KA
3. Hill DS
4. Anfosso A
5. Mrass P
6. Munoz MA
7. Kinjyo I
8. Weninger W
(2014) Real-time cell cycle imaging during melanoma growth, invasion, and drug response
Pigment Cell & Melanoma Research 27:764–776.

https://doi.org/10.1111/pcmr.12274
- PubMed
- Google Scholar
(2014) ODE constrained mixture modelling: a method for unraveling subpopulation structures and dynamics
PLOS Computational Biology 10:e1003686.

https://doi.org/10.1371/journal.pcbi.1003686
- PubMed
- Google Scholar
1. Herlyn M
2. Thurin J
3. Balaban G
4. Bennicelli JL
5. Herlyn D
6. Elder DE
7. Bondi E
8. Guerry D
9. Nowell P
10. Clark WH
(1985)
Characteristics of cultured human melanocytes isolated from different stages of tumor progression

Cancer Research 45:5670–5676.
- PubMed
- Google Scholar
1. Herlyn M
(1990) Human melanoma: development and progression
Cancer Metastasis Reviews 9:101–112.

https://doi.org/10.1007/BF00046337
- PubMed
- Google Scholar
(2014) Determination of parameter identifiability in nonlinear biophysical models: A Bayesian approach
The Journal of General Physiology 143:401–416.

https://doi.org/10.1085/jgp.201311116
- PubMed
- Google Scholar
(2010) Multicellular tumor spheroids: an underestimated tool is catching up again
Journal of Biotechnology 148:3–15.

https://doi.org/10.1016/j.jbiotec.2010.01.012
- PubMed
- Google Scholar
1. Hoek KS
2. Schlegel NC
3. Brafford P
4. Sucker A
5. Ugurel S
6. Kumar R
7. Weber BL
8. Nathanson KL
9. Phillips DJ
10. Herlyn M
11. Schadendorf D
12. Dummer R
(2006) Metastatic potential of melanomas defined by specific gene expression profiles with no BRAF signature
Pigment Cell Research 19:290–302.

https://doi.org/10.1111/j.1600-0749.2006.00322.x
- PubMed
- Google Scholar
1. Ivascu A
2. Kubbies M
(2006) Rapid generation of single-tumor spheroids for high-throughput cell function and toxicity analysis
Journal of Biomolecular Screening 11:922–932.

https://doi.org/10.1177/1087057106292763
- PubMed
- Google Scholar
1. Jin W
2. Spoerri L
3. Haass NK
4. Simpson MJ
(2021) Mathematical model of tumour spheroid experiments with real-time cell cycle imaging
Bulletin of Mathematical Biology 83:44.

https://doi.org/10.1007/s11538-021-00878-4
- PubMed
- Google Scholar
Software
1. Johnson SG
(2021) The NLopt module for Julia, version v0.6.2
Github.

https://github.com/JuliaOpt/NLopt.jl
1. Kienzle A
2. Kurch S
3. Schlöder J
4. Berges C
5. Ose R
6. Schupp J
7. Tuettenberg A
8. Weiss H
9. Schultze J
10. Winzen S
11. Schinnerer M
12. Koynov K
13. Mezger M
14. Haass NK
15. Tremel W
16. Jonuleit H
(2017) Dendritic Mesoporous Silica Nanoparticles for pH-Stimuli-Responsive Drug Delivery of TNF-Alpha
Advanced Healthcare Materials 6:201700012.

https://doi.org/10.1002/adhm.201700012
- PubMed
- Google Scholar
(2012) The multicellular tumor spheroid model for high-throughput cancer drug discovery
Expert Opinion on Drug Discovery 7:819–830.

https://doi.org/10.1517/17460441.2012.708334
- PubMed
- Google Scholar
(2020) Biologically-informed neural networks guide mechanistic modeling from sparse experimental data
PLOS Computational Biology 16:e1008462.

https://doi.org/10.1371/journal.pcbi.1008462
- PubMed
- Google Scholar
1. Langhans SA
(2018) Three-Dimensional in Vitro Cell Culture Models in Drug Discovery and Drug Repositioning
Frontiers in Pharmacology 9:6.

https://doi.org/10.3389/fphar.2018.00006
- PubMed
- Google Scholar
(2013) Multicellular tumor spheroid models to explore cell cycle checkpoints in 3D
BMC Cancer 13:73.

https://doi.org/10.1186/1471-2407-13-73
- PubMed
- Google Scholar
Book
(1998) Theory of Point Estimation
Springer.

https://doi.org/10.1007/b98854
- Google Scholar
(2013) Growth of confined cancer spheroids: a combined experimental and mathematical modelling approach
Integrative Biology 5:597–605.

https://doi.org/10.1039/c3ib20252f
- PubMed
- Google Scholar
1. Mark C
2. Grundy TJ
3. Strissel PL
4. Böhringer D
5. Grummel N
6. Gerum R
7. Steinwachs J
8. Hack CC
9. Beckmann MW
10. Eckstein M
11. Strick R
12. O’Neill GM
13. Fabry B
(2020) Collective forces of tumor spheroids in three-dimensional biopolymer networks
eLife 9:e51912.

https://doi.org/10.7554/eLife.51912
- PubMed
- Google Scholar
1. Markovsky I
2. Van Huffel S
(2007) Overview of total least-squares methods
Signal Processing 87:2283–2302.

https://doi.org/10.1016/j.sigpro.2007.04.004
- Google Scholar
1. Masuda F
2. Ishii M
3. Mori A
4. Uehara L
5. Yanagida M
6. Takeda K
7. Saitoh S
(2016) Glucose restriction induces transient G2 cell cycle arrest extending cellular chronological lifespan
Scientific Reports 6:19629.

https://doi.org/10.1038/srep19629
- PubMed
- Google Scholar
Website
1. Mathworks
(2021) Texture analysis
Accessed April 4, 2021.

https://au.mathworks.com/help/images/texture-analysis-1.html
1. Murphy KC
2. Hung BP
3. Browne-Bourne S
4. Zhou D
5. Yeung J
6. Genetos DC
7. Leach JK
(2017) Measurement of oxygen tension within mesenchymal stem cell spheroids
Journal of the Royal Society, Interface 14:20160851.

https://doi.org/10.1098/rsif.2016.0851
- PubMed
- Google Scholar
Preprint
(2021) Designing and Interpreting 4D Tumour Spheroid Experiments
bioRxiv.

https://doi.org/10.1101/2021.08.18.456910
- Google Scholar
Book
1. Pawitan Y
(2013)
In All Likelihood: Statistical Modelling and Inference Using Likelihood

Oxford University Press.
- Google Scholar
Book
1. Powell MJD
(2009)
The BOBYQA Algorithm for Bound Constrained Optimization without Derivatives

Department of Applied Mathematics and Theoretical Physics.
- Google Scholar
(2017) Delineating parameter unidentifiabilities in complex models
Physical Review. E 95:032314.

https://doi.org/10.1103/PhysRevE.95.032314
- PubMed
- Google Scholar
1. Raue A
2. Kreutz C
3. Maiwald T
4. Bachmann J
5. Schilling M
6. Klingmüller U
7. Timmer J
(2009) Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood
Bioinformatics 25:1923–1929.

https://doi.org/10.1093/bioinformatics/btp358
- PubMed
- Google Scholar
Preprint
(2016) Forward-Mode Automatic Differentiation in Julia
arXiv.

https://arxiv.org/abs/1607.07892
- Google Scholar
(2007) Mathematical Models of Avascular Tumor Growth
SIAM Review 49:179–208.

https://doi.org/10.1137/S0036144504446291
- Google Scholar
1. Sakaue-Sawano A
2. Kurokawa H
3. Morimura T
4. Hanyu A
5. Hama H
6. Osawa H
7. Kashiwagi S
8. Fukami K
9. Miyata T
10. Miyoshi H
11. Imamura T
12. Ogawa M
13. Masai H
14. Miyawaki A
(2008) Visualizing spatiotemporal dynamics of multicellular cell-cycle progression
Cell 132:487–498.

https://doi.org/10.1016/j.cell.2007.12.033
- PubMed
- Google Scholar
(2009) Melanocytes: from morphology to application
Skin Pharmacology and Physiology 22:114–121.

https://doi.org/10.1159/000178870
- PubMed
- Google Scholar
1. Sarapata EA
2. de Pillis LG
(2014) A comparison and catalog of intrinsic tumor growth models
Bulletin of Mathematical Biology 76:2010–2024.

https://doi.org/10.1007/s11538-014-9986-y
- PubMed
- Google Scholar
(2020) Practical parameter identifiability for spatio-temporal models of cell invasion
Journal of the Royal Society, Interface 17:20200055.

https://doi.org/10.1098/rsif.2020.0055
- PubMed
- Google Scholar
(2007a) An organometallic protein kinase inhibitor pharmacologically activates p53 and induces apoptosis in human melanoma cells
Cancer Research 67:209–217.

https://doi.org/10.1158/0008-5472.CAN-06-1538
- PubMed
- Google Scholar
1. Smalley KSM
2. Contractor R
3. Haass NK
4. Lee JT
5. Nathanson KL
6. Medina CA
7. Flaherty KT
8. Herlyn M
(2007b) Ki67 expression levels are a better marker of reduced melanoma growth following MEK inhibitor treatment than phospho-ERK levels
British Journal of Cancer 96:445–449.

https://doi.org/10.1038/sj.bjc.6603596
- PubMed
- Google Scholar
1. Smalley KS
2. Lioni M
3. Noma K
4. Haass NK
5. Herlyn M
(2008) In vitro three-dimensional tumor microenvironment models for anticancer drug discovery
Expert Opinion on Drug Discovery 3:1–10.

https://doi.org/10.1517/17460441.3.1.1
- PubMed
- Google Scholar
(2017) Real-Time Cell Cycle Imaging in a 3D Cell Culture Model of Melanoma
Methods in Molecular Biology 1612:401–416.

https://doi.org/10.1007/978-1-4939-7021-6_29
- PubMed
- Google Scholar
(2021) Fluorescence-Based Quantitative and Spatial Analysis of Tumour Spheroids: A Proposed Tool to Predict Patient-Specific Therapy Response
Frontiers in Digital Health 3:668390.

https://doi.org/10.3389/fdgth.2021.668390
- PubMed
- Google Scholar
(2020) Marked synergy by vertical inhibition of EGFR signaling in NSCLC spheroids shows SOS1 is a therapeutic target in EGFR-mutated cancer
eLife 9:e58204.

https://doi.org/10.7554/eLife.58204
- PubMed
- Google Scholar
1. Wallace DI
2. Guo X
(2013) Properties of tumor spheroid growth exhibited by simple mathematical models
Frontiers in Oncology 3:51.

https://doi.org/10.3389/fonc.2013.00051
- PubMed
- Google Scholar
1. Wang L
2. Cao J
3. Ramsay JO
4. Burger DM
5. Laporte CJL
6. Rockstroh JK
(2012) Estimating mixed-effects differential equation models
Statistics and Computing 24:111–121.

https://doi.org/10.1007/s11222-012-9357-1
- PubMed
- Google Scholar
1. Ward JP
2. King JR
(1997) Mathematical modelling of avascular-tumour growth
IMA Journal of Mathematics Applied in Medicine and Biology 14:39–69.

https://doi.org/10.1093/imammb/14.1.39
- PubMed
- Google Scholar
1. Ward JP
2. King JR
(1999) Mathematical modelling of avascular-tumour growth. II: Modelling growth saturation
IMA Journal of Mathematics Applied in Medicine and Biology 16:171–211.

https://doi.org/10.1093/imammb/16.2.171
- PubMed
- Google Scholar
(2015) Spherical cancer models in tumor biology
Neoplasia 17:1–15.

https://doi.org/10.1016/j.neo.2014.12.004
- PubMed
- Google Scholar
1. Welch BL
(1947) The generalisation of student’s problems when several different population variances are involved
Biometrika 34:28–35.

https://doi.org/10.1093/biomet/34.1-2.28
- PubMed
- Google Scholar

Article and author information

Author details

Alexander P Browning
1. School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia
2. ARC Centre of Excellence for Mathematical and Statistical Frontiers, Queensland University of Technology, Melbourne, Australia
Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing

Contributed equally with
Jesse A Sharp

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8753-1538
Jesse A Sharp
1. School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia
2. ARC Centre of Excellence for Mathematical and Statistical Frontiers, Queensland University of Technology, Melbourne, Australia
Contribution
Conceptualization, Formal analysis, Investigation, Writing – review and editing

Contributed equally with
Alexander P Browning

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-2865-4853
Ryan J Murphy

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia

Contribution
Data curation, Investigation, Methodology, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-9844-6734
Gency Gunasingh

The University of Queensland Diamantina Institute, The University of Queensland, Brisbane, Australia

Contribution
Data curation, Investigation, Writing – review and editing

Competing interests
No competing interests declared
Brodie Lawson
1. School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia
2. ARC Centre of Excellence for Mathematical and Statistical Frontiers, Queensland University of Technology, Melbourne, Australia
Contribution
Methodology, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1317-5988
Kevin Burrage
1. School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia
2. ARC Centre of Excellence for Mathematical and Statistical Frontiers, Queensland University of Technology, Melbourne, Australia
3. Department of Computer Science, University of Oxford, Oxford, United Kingdom
Contribution
Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8111-1137
Nikolas K Haass

The University of Queensland Diamantina Institute, The University of Queensland, Brisbane, Australia

Contribution
Funding acquisition, Investigation, Project administration, Writing – review and editing

Contributed equally with
Matthew Simpson

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-3928-5360
Matthew Simpson

School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia

Contribution
Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing – review and editing

Contributed equally with
Nikolas K Haass

For correspondence
matthew.simpson@qut.edu.au

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6254-313X

Funding

Australian Research Council (DP200100177)

Nikolas K Haass
Matthew Simpson

ARC Centre of Excellence for Mathematical and Statistical Frontiers (CE140100049)

Alexander P Browning
Jesse A Sharp

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Patrick Thomas for helpful comments and discussions and John Blake for guidance using the Incucyte S3. We thank Jennifer Flegg and one anonymous referee for helpful comments. NKH and MJS are supported by the Australian Research Council (DP200100177). APB and JAS are supported by the ARC Centre of Excellence for Mathematical and Statistical Frontiers (CE140100049). This research was carried out at the Translational Research Institute (TRI), Woolloongabba, QLD. TRI is supported by a grant from the Australian Government. We thank the staff in the microscopy core facility at TRI for their technical support. We thank Prof. Atsushi Miyawaki, RIKEN, Wako-city, Japan, for providing the FUCCI constructs, Prof. Meenhard Herlyn and Ms. Patricia Brafford, The Wistar Institute, Philadelphia, PA, for providing the cell lines.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.