Designing and interpreting 4D tumour spheroid experiments

Tumour spheroid experiments are routinely used to study cancer progression and treatment. Various and inconsistent experimental designs are used, leading to challenges in interpretation and reproducibility. Using multiple experimental designs, live-dead cell staining, and real-time cell cycle imaging, we measure necrotic and proliferation-inhibited regions in over 1000 4D tumour spheroids (3D space plus cell cycle status). By intentionally varying the initial spheroid size and temporal sampling frequencies across multiple cell lines, we collect an abundance of measurements of internal spheroid structure. These data are difficult to compare and interpret. However, using an objective mathematical modelling framework and statistical identifiability analysis we quantitatively compare experimental designs and identify design choices that produce reliable biological insight. Measurements of internal spheroid structure provide the most insight, whereas varying initial spheroid size and temporal measurement frequency is less important. Our general framework applies to spheroids grown in different conditions and with different cell types.


A Mathematical model A.1 Nutrient only
Here we recall Greenspan's mathematical model governing the evolution of R o (t), R i (t), and R n (t) (Figures 1f-g, S1a-c). We consider conservation of mass to govern the evolution of R o (t). Assuming: i) all living cells are identical and an incompressible mass of constant volume; ii) cell division occurs instantaneously relative to the growth time of the tumour, and each daughter cell occupies the same volume as any other cell; iii) the proliferation rate is a constant, s, for cells which have sufficient nutrient; and, iv) the mass density of living cells is constant and equal to density of necrotic debris; then conservation of mass is equivalent to conservation of volume, giving, where A is the total volume of living cells at any time, t; B is the initial volume of living cells at time t = 0; C is the total volume of cells produced in t ≥ 0; D is the total volume of necrotic debris at time t; and, E is the total volume lost in the necrotic core in t ≥ 0. Greenspan's mathematical model describes the three phases of growth and time evolution of the outer radius, R o (t) (green), inhibited radius, R i (t) (magenta), and necrotic radius, R n (t) (black). (b) Schematic for Greenspan mathematical model. Nutrient diffuses within the tumour spheroid and is consumed by living cells. (c) Snapshot of nutrient concentration, c(r(t)) for 0 < r < R o (t), for a tumour spheroid in phase (iii). External nutrient concentration is c ∞ . Inhibited radius, R i (t), and necrotic radius, R n (t), are defined as the radii where the nutrient concentration first reaches the thresholds c i and c n , respectively. of a sphere with radius r is 4πr 3 /3 and 4πr 2 , respectively, gives proliferation of living cells where H(·) is the Heaviside step function.
To determine the full evolution of the system we solve equations (S.4) and (S.5) together with the nutrient threshold for inhibition, c i , and the nutrient threshold for necrosis, c n , which implicitly define the radius of the inhibited region, R i (t), and the radius of the necrotic region, R n (t), respectively, if the nutrient concentration inside the spheroid is sufficiently small otherwise R i (t) = 0 or R n (t) = 0.
Note that the equation (S.4) for nutrient does not involve any temporal derivative so the only initial condition required to solve the full system of equations (S.4) and (S.5) is the initial outer radius, The solution of equation (S.5) is, where c ∞ − c n = α 3k The necrotic region first forms when the nutrient concentration reaches c n at the centre, which occurs when R n (t) = 0 and r = 0 in equation (S.8), which gives a critical outer radius, Also recall that R i (t) corresponds to c(R i (t), t) = c i which we can substitute into equation (S. 7) to give, . (S.10) Since the inhibited region first forms when the nutrient concentration reaches c i at the centre of the spheroid and the necrotic region forms after the inhibited region, setting R n (t) = 0 and r = 0 on right-hand-side of equation (S.10) gives the outer radius when the inhibited region first forms We can then define a useful dimensionless quantity, Q 2 = R 2 d /R 2 c = (c ∞ − c i ) / (c ∞ − c n ), which is related to the time when phase (ii) begins.
Equations (S.4), (S.5), and (S.6) can now be solved in each of phase (i), (ii), and (iii). To provide valuable insights into the structure of the solutions to the Greenspan model it helps to consider the non-dimensional form of the equations and their solutions. To non-dimensionalise we rescale time with s to give τ = st and rescale lengths with R c via ξ o (t) = R o (t)/R c , ξ i (t) = R i (t)/R c , and ξ n (t) = R n (t)/R c . Then phase (ii) starts when ξ o (t) = Q and phase (iii) starts when ξ o (t) = 1. We now consider each phase in turn.

Phase (i)
In phase (i), all cells are free to proliferate and the nutrient concentration is sufficiently high, i.e. c(r, t) > c i for 0 ≤ r ≤ R o (t), such that there is no inhibited or necrotic region (Figure 1(a)). Phase (i) ends when the nutrient concentration at the centre of the spheroid equals the inhibited threshold, when c(0, t) = c i and R o (t) = R d .
Since R i (t) = 0 and R n (t) = 0, equation (S.4) becomes 12) giving, Non-dimensionalising gives, 14) Given the solution in equation (S.14) we determine R o (t) by reintroducing s and R c , Note that R i (t) = 0 and R n (t) = 0 throughout phase (i). Hence, we have obtained R o (t), R i (t), R n (t) throughout this phase.

Phase (ii)
In phase (ii) the spheroid experiences inhibited growth due to a core of inhibited cells and outer region of freely proliferating cells (Figure 1(b)). Phase (ii) ends when the necrotic core forms. Since with initial condition y i (τ ) = 0 at τ = τ 1 and terminating condition y i (τ ) 2 = 1 − Q 2 . The constraint used to derive equation (S.18) and the termination condition for phase (ii) are obtained with the 6 following argument. In phase (ii) equation (S.10) is Non-dimensionalising equations (S.8) and (S.19), using definitions of ξ o (τ ), ξ i (τ ), Q, and combining the resulting expressions gives , which gives the constraint used to derive equation (S.18). Using the fact that phase (ii) ends when ξ o (τ ) = 1 and rearranging gives the termination condition for y i (τ ).
Numerically solving equation (S.18), using MATLAB's in-built ode15s differential equation solver [1] with absolute and relative tolerances set to 1 × 10 −8 , we obtain y i (τ ) for phase (ii). To obtain Recall R n (t) = 0 throughout phase (ii). Hence, we have obtained R o (t), R i (t), R n (t) throughout this phase.

Phase (iii)
In phase (iii) the spheroid experiences inhibited growth due to an outer proliferating region, an intermediate region of inhibited cells, and a necrotic core (Figure 1(a)). At steady state there is a balance between the number of cells that are proliferating in the outer region and mass lost from the necrotic core.
Given the solution for y n (τ ) and y i (τ ) we obtain R o (t), R i (t) and R n (t) using the following approach.
Greenspan's mathematical model assumes that tumour spheroids experience three phases of growth [2]. While we find experimental evidence confirming that many tumour spheroids experience three phases of growth (Figures 1, S10, S11, S12, S14, and S15), we also find experimental evidence suggesting tumour spheroids seeded with a higher number of cells may form in phase (ii) (Figures S13 and S16 is the unknown variable. We determine the three solutions of this polynomial using the MATLAB function roots [3] and define R n (0) as the only physically realistic real-valued solution which satisfies 0 < R n (0) < R o (0). Similarly, to obtain R i (0), we rewrite equation (S.10) as the following cubic is the unknown variable and R o (0) and R n (0) are known. We then define R i (0) as the only physically realistic real-valued solution which satisfies R n (0) < R i (0) < R o (0). For statistical identifiability analysis we assume spheroids may form in phase (i), phase (ii), or phase (iii).
The approach taken in phases (i), (ii), and (iii) means that we do not require knowledge of the values of the parameters c ∞ , c n , c i , k and α but instead only the value of This reduces the number of parameters describing the evolution of the spheroid from eight to five.
The three pieces of information no longer consider regard the nutrient concentration which we do not directly measure in this study and has been explored in other studies [4,5]. Furthermore, equations (S.8)-(S.10) show that there are constraints on the relationships between R o (t), R i (t), R n (t) which can be explored further.

A.2 Nutrient and waste
The mathematical model presented in the methods section of the main manuscript and supplementary discussion A.1 is a special case of Greenspan's model [2]. The general Greenspan model proposes the inhibited region is a result of a build up of waste produced from live or dead cells and the necrotic region forms due to a lack of nutrient. Here, we consider the alternative case where waste is produced from live cells only and show that, for the measurements we obtain, it is equivalent to the nutrient only case we consider in the main manuscript ( Figure S2). We do not consider waste produced only from dead cells in this study since with that model the necrotic core must form before the inhibited region which is not what we observe in these experiments (Figure 1(b)).
In comparison to the nutrient only model in supplementary discussion A.1, the model with nutrient and waste requires an additional equation for the evolution of waste concentration, β(r, t).
The full system of governing equations are, where β i κ/P are the parameters related to nutrient and α/k(c ∞ − c n ) are the parameters related to waste. This new definition of Q provides a different interpretation of the data since Q now represents a combination of waste and nutrient parameters. Importantly, with this new definition of Q there are two cases to consider: i) Q ≤ 1, and ii) Q > 1. Previously, we only considered case (i). In case (ii) where Q > 1 the necrotic core forms before the inhibited region. We do not observe this scenario in the experiments that we perform and therefore we restrict the attention of this study to case (i) Lack of nutrient forms necrotic region Build up of waste produced from live cells forms inhibited region Parameter identifiability using statistical profile likelihood analysis is outlined in the methods section of the main manuscript. We now provide further details.
We partition the full set of observations y o into sets of observations y o o , y o n , and y o i corresponding to experimental measurements of R o (t), R n (t), and R i (t). For computational accuracy, we perform calculations using the log-likelihood which is, assuming data independence, and, σ 2 o , σ 2 n , and σ 2 i correspond to pooled variances of the three measurement types R o (t), R n (t), and R i (t), respectively [7,8]. We approximate σ 2 o ≈ s 2 o , σ 2 n ≈ s 2 n , and σ 2 i ≈ s 2 i , where s 2 o , s 2 n , s 2 i are pooled sample variances of the outer, necrotic, and inhibited radius measurements, respectively [9]. The pooled sample variance for the outer radius is defined as log f (y o n,j ; y n,j (ψ, ϕ), σ 2 n )

(S.29)
Given the five-dimensional parameter space that we are searching to find the maximum likelihood estimate and the four-dimensional parameter space we search to find profile likelihoods, we sequentially determine the maximum likelihood estimate (MLE) and profile likelihoods. All subsequent minimisation optimisations are performed using functions in MATLABs global optimisation toolbox. Specifically, we use the GlobalSearch function [10] where we create the following optimisation problem structure. We set the local solver to be the fmincon function using the sequential quadratic programming (sqp) algorithm, MaxIterations = 2500 and MaxFunctionEvaluations = 5000. The objective function is defined as the argument of the minimisation of the right hand side of equation 3. Thirdly, we calculate an estimate of the confidence intervals using a profile likelihood threshold value of 0.15, which can be approximately calibrated via simulation or the χ 2 −distribution [11].
Specifically, we start at either end of the simple parameter bounds until we determine the first grid point where the normalised profile likelihood, L p ψ; y 0 = exp l p ψ; y 0 − max θ l θ; y 0 , is greater than 0.15. We then set new lower and upper bounds as being two grid points to the left or right of that location, respectively. Note that a more sophisticated approach to determine the approximate 95% confidence intervals is applied in step seven to compute the results shown in Table S3, which is not required here.
4. Fourth, we repeat the search for the maximum likelihood estimate using the new lower and upper bounds with the same settings as we first used.
5. Fifth, we repeat the calculations for the profile likelihoods using the new lower and upper bounds.
6. Sixth, we determine the maximum likelihood estimate to be the value across all calculations which maximises the likelihood. We form the final profile likelihood from steps two and four and present the normalised likelihood function, L p ψ; y 0 = exp l p ψ; y 0 − max θ l θ; y 0 , in figures.
7. To compute approximate 95% confidence intervals for each parameter, as shown in Table S3, we form the profile likelihood from steps two and four. Next, we start at either end of the simple parmater bounds until we determine the first grid point where the normalised profile likelihood, L p ψ; y 0 = exp l p ψ; y 0 − max θ l θ; y 0 , is greater than profile likelihood threshold value of 0.15 [11]. These two grid points are a first approximation of the lower and upper 95% confidence interval boundaries. Finally, to obtain a more accurate estimate of the approximate 95% confidence interval boundaries, we consider each of these two grid points in turn as the FirstGuess for the MATLAB function fsolve [12], and use linear interpolation, with the MATLAB function interp1 [13].
Key software, including precise details of the implementation of this section, is freely available on a GitHub repository (https://github.com/ryanmurphy42/4DSpheroids Murphy2021).

B.2 Parameter bounds
To interpret s we consider the evolution of the tumour spheroid in phase (i). Equation (S.13) can be written in terms of volume V (t), recalling that the volume of a sphere is 4πR 3 o (t)/3, as

B.3 Pooled sample variances
To identify parameters we assume distinct pooled sample variances for the outer, necrotic, and inhibited radius measurements, as opposed to a single pooled variance for all measurements. In Figure S3 we plot the pooled sample variances for different experimental designs which justify the use of a variance for each measurement type.

C Experimental data C.1 Outer radius experimental measurements and images
The IncuCyte S3 live cell imaging system is a useful tool that we use to obtain many outer radius measurements. The other outer radius measurements are obtained from confocal microscopy.
We start with 24 spheroids in the IncuCyte S3 live cell imaging system for each cell line and initial spheroid size and image every 6 hours for the duration of the experiment. However, some measurements could not be obtained primarily due to blurring of the automated imaging, spheroids not forming properly, or spheroids losing their structural integrity at very late time. In Table S1 we show the total number of measurements obtained at 24 hour intervals starting from Day 0 which corresponds to the time that we determined as when spheroid formation ends and growth begins

C.3 Confocal microscopy supplementary experimental images
Here we present confocal microscopy images of spheroids formed with the WM793b, WM983b, and WM164 cell lines. In the images we outline each spheroids outer boundary, inhibited region, and necrotic region.

C.3.4 3D rendering
Here we present a 3D rendering of a confocal microscopy image z-stack of half of a FUCCI-melanoma WM793b spheroid 17 days after formation with 5000 cells.
Figure S18: 3D rendering of half of a FUCCI-melanoma WM793b spheroid 17 days after formation with 5000 cells. Scale bar 200 µm.

D.1 Results in tables
In all figures with profile likelihoods we include a red-dashed horizontal line at 0.15 indicating the 95% confidence interval threshold value [11]. Here, in Table S3 Figure 1 shows that varying the temporal resolution in Design 1 is not sufficient to predict necrotic and inhibited radii. Here, in Figures S19 and S20, we show that varying the temporal resolution using Designs 2 and 3, respectively, gives consistent results across temporal resolutions A, B, and C.

D.2 Measurement times and experimental duration
Next we consider four additional experimental designs that use different temporal measurements Temporal Resolution D: the first 4 days (Day 1, 2, 3

D.3 Profile likelihoods for R o (0)
To perform statistical identifiability analysis we we treat the initial outer radius, R o (0), as a parameter. Here, in Figure S23 we show that profile likelihoods for R o (0) are consistent across temporal resolutions and experimental designs.

E Synthetic data: WM793b
To confirm that profile likelihood analysis works as expected, we generate synthetic data from Greenspan's mathematical model using known parameters. We then explore when these known parameters are recovered using the varying experimental designs considered in the main manuscript: with these known parameters. Next, to obtain one noisy synthetic outer radius measurement we record the outer radius from Greenspan's model generated from the known parameters at one time point. Next, we sample a normal distribution with zero mean and variance given by experimentally obtained outer radius pooled sample variance s 2 o = 9.35. We add this sampled noise to the recorded outer radius measurement. We repeat this process to obtain additional outer radius measurements.
Similarly, we repeat this process to obtain necrotic and inhibited radius measurements, using experimentally obtained pooled sample variances s 2 n = 15.89, and s 2 i = 33.12, respectively. We generate 10 measurements, or 48 when exploring the role of additional measurements in supplementary discussion E.4, of the outer radius, inhibited radius and necrotic radius every half day from day 0 to day 20.

E.1 Outer radius measurements are not sufficient to predict inhibited and necrotic radii
Similarly to Figure 2, we observe in Figure S30  However, inspecting the profile likelihoods in Figures S30g-j shows that, while the known parameters are captured, the profiles are wide suggesting that parameters are non-identifiable. This means that many parameter values give a similar match to the outer radius experimental data and these 49 parameters do not necessarily agree with the inhibited and necrotic radii measurements.   Figure S31).

E.3 Role of initial spheroid size and experiment duration
In Greenspan's model a change in the initial radius, R o (0), corresponds to a shift in time ( Figure   S32a). We now consider the role of initial spheroid size and experiment duration. As before, we use the MLE from spheroids formed with 5000 cells per spheroid to generate synthetic data. To generate synthetic data for spheroids formed with 1250, 2500, and 10000 cells per spheroid we use the MLE obtained from spheroids formed with 5000 cells per spheroid and only update R o (0). To update R o (0) we use the MLE from Design 3 applied to experimental data obtained from WM793b spheroids formed with 1250, 2500, and 10000 cells per spheroid, respectively.
We assume that each experiment is performed to Day 6 after formation, and use Design 3 with

E.4 Increasing number of measurements
In biological experiments it is time consuming and expensive to increase the number of measurements obtained. However, by generating synthetic data we can easily simulate additional measurements.
We generate 48 measurements of the outer radius, inhibited radius, and necrotic radius every half day from Day 0 to Day 20. We choose 48 measurements since this corresponds to half of a 96-well plate and is extremely large in comparison to typical experiments. These results show that many measurements of Design 2 may provide good insight in this extreme scenario but that Design 3 still provides most insight.

F Parameter identifiability analysis for WM983b
The main manuscript focuses on results for the human melanoma WM793b cell line. Here, we present the corresponding results for the human melanoma WM983b spheroids formed with 2500, 5000, and 10000 cells.
All key observations made in reference to the WM793b cell line hold for the WM983b cell line.
Specifically, in Figures S34, S35 and S36, for spheroids formed with 2500, 5000, and 10000 cells, respectively, we show that varying the temporal resolution using only Design 1 is insufficient to determine necrotic and inhibited radii. In Figures S37, S38, and S39 for spheroids formed with 2500, 5000, and 10000 cells, respectively, we show that Design 3 provides most insight. In Figure S40 we show that information gained across experiments with different initial spheroid sizes is consistent.

G Parameter identifiability analysis for WM164
The main manuscript focuses on results for the human melanoma WM793b cell line. Here, we present analogous results for the human melanoma WM164 spheroids formed with 1250, 2500, 5000, and 10000 cells per spheroid. These spheroids are more challenging to interpret as we will now explain.
In experiments WM164 spheroids formed after 3 days. These spheroids were larger in size than other spheroids considered in this work, with the initial radius of WM164 spheroids formed with 1250 cells per spheroid larger than and similar size to WM983b and WM793b spheroids formed with 10000 cells per spheroid, respectively. The WM164 spheroids had relatively poor spherical symmetry [14], grew rapidly and many spheroids lost their structural integrity nine days after seeding formed with 1250, 2500, and 5000 cells per spheroid, and seven days after seeding for spheroids formed with 10000 cells spheroid. In addition, confocal microscopy could not be performed on day 7 after seeding for spheroids formed with 5000 and 10000 cells per spheroid due loss of structural integrity during harvesting. Identification of the necrotic region using image processing was more challenging, than for other cell lines, as a well-defined necrotic region did not form prior to the termination of the experiment. Therefore, necrotic region measurements for these spheroids are more subjective and uncertain. Spheroid boundaries were less smooth, so settings on the IncuCyte S3 live cell imaging system were updated to measure the largest brightfield object area with max eccentricity to 0.75 and sensitivity 20. These outer radius measurements were then manually reviewed to confirm accuracy.
We perform analysis for WM164 spheroids using Days 1, 2, 3, 4 and 5 after formation, where measurements could be obtained. This means that we do not include the last day of outer radius measurements for spheroids formed with 1250, 2500, and 5000 cells per spheroid. This allows us to compare the final outer radius measurement to Greenspan's model simulated with the MLE as a predictive test. However, for spheroids formed with 10000 cells per spheroid we include all data points so cannot form a predictive test, but this is because we seek to obtain as much information as possible in the shorter experimental duration. While the experimental duration for WM164 spheroids is relatively short in comparison to the WM793b and WM983b experiments, these experiments are still performed for multiple days longer than previous WM164 spheroid experiments [15].
To perform the analysis we update initial parameter bounds, used for practical parameter iden- insight, here we compare results obtained from spheroids with different initial sizes using Design 3.
In Figure S41a, we observe four distinct narrow peaks for R o (0) corresponding to the four initial spheroid sizes, which is expected. For s, we observe that profile likelihoods overlap showing information obtained for s is relatively consistent for different initial spheroid sizes ( Figure S41b).
Interestingly and importantly, we observe four distinct peaks for R c ( Figure S41c). This lack of consistency is different to the other two cell lines considered and strongly suggests information gained across initial spheroid sizes is not consistent. This is supported by direct observation of the experimental data where spheroids formed with 2500 cells have necrotic cores on day 4, whereas similar sized spheroids on day 2 formed with 5000 and 10000 cells per spheroid do not. Profile likelihoods for γ are wide for spheroids formed with 5000 and 10000 cells per spheroid, and narrow and overlapping for spheroids formed with 1250 and 2500 cell per spheroid, showing that γ requires more necrotic core measurements to be identified ( Figure S41d). Profile likelihoods for Q suggest that Q decreases as the initial spheroid size increases ( Figure S41e). This result for Q is less consistent and in constrast to results from other cell lines, where the profiles for Q overlapped for all spheroid sizes. Overall, we conclude that, possibly due to their lack of spherical symmetry, WM164 spheroids are more challenging to interpret and information gained using spheroids of different sizes is not consistent.
To support these results, we show along the diagonal of Figure S41f the solution of the mathematical model evaulated at the MLE associated with each initial spheroid size compared to the experimental measurements. In doing so we demonstrate that we accurately predict the last outer radius measurement using previous days measurements for spheroids formed with 1250, 2500, and 5000 cells per spheroid. However, on the off-diagonals of Figure 4f, we compare how the Greenspan model simulated with the MLE from one initial spheroid size predicts data from different initial spheroid sizes by only changing the initial radius. These off-diagonal results show that using information from one spheroid size to predict the behaviour of a different spheroid size is not always accurate. For example, using information gained from spheroids formed with 10000 cell per spheroid poorly predicts the behaviour of spheroids formed with 1250 cell per spheroid, as the time evolution of the outer radius is not accurately predicted at late time and the inhibited and necrotic regions form much earlier than predicted (top right of Figure S41f).