Growth and adaptation mechanisms of tumour spheroids with time-dependent oxygen availability

Abstract

Tumours are subject to external environmental variability. However, in vitro tumour spheroid experiments, used to understand cancer progression and develop cancer therapies, have been routinely performed for the past fifty years in constant external environments. Furthermore, spheroids are typically grown in ambient atmospheric oxygen (normoxia), whereas most in vivo tumours exist in hypoxic environments. Therefore, there are clear discrepancies between in vitro and in vivo conditions. We explore these discrepancies by combining tools from experimental biology, mathematical modelling, and statistical uncertainty quantification. Focusing on oxygen variability to develop our framework, we reveal key biological mechanisms governing tumour spheroid growth. Growing spheroids in time-dependent conditions, we identify and quantify novel biological adaptation mechanisms, including unexpected necrotic core removal, and transient reversal of the tumour spheroid growth phases.

Introduction

In vivo tumours are subject to various types of environmental variability, for example due to fluctuating oxygen and nutrient availability [1–4]. To study cancer progression and develop cancer therapies, tumour spheroid experiments have been successfully and routinely performed for the past fifty years [2, 5–12]. However, tumour spheroid experiments are typically performed in constant environments and focus on the overall size of spheroids [2, 13–18]. By experimentally controlling oxygen availability and using mathematical modelling and statistical uncertainty quantification, we develop a new framework to study the impact of external environmental variability on the growth of tumour spheroids and their internal structure. Using our framework we identify and quantify novel biological adaptation mechanisms driven by environmental variability. This work begins to bridge the gap between in vitro and in vivo conditions, and lays the foundation for future experimental, mathematical, and statistical spheroid studies.

Oxygen availability is of particular importance since it is vital to the effectiveness of cancer therapies, such as chemotherapy and radiotherapy [1, 19, 20], and can be controlled in spheroid experiments. However, spheroid experiments are typically performed in ambient atmospheric conditions (21% oxygen), sometimes referred to as normoxia [13, 16]. In contrast, untreated tumours typically grow in variable hypoxic conditions (0.3–4.2% oxygen) [1, 20–22]. While many single-cell studies, and some spheroid studies, explore the role of environmental variability [2, 22–27], oxygen parameters critical to reproduce results are commonly not reported [28].

To visualise spheroid growth in normoxia, hypoxia, and time-dependent oxygen conditions we use fluorescent ubiquitination cell cycle indicator (FUCCI) transduced cell lines and hypoxia markers (Fig 1A–1E) [13, 14, 29, 30]: nuclei of cells in gap 1 (G1) phase fluoresce red, shown in magenta for clarity (Fig 1D); nuclei of cells in synthesis, gap 2, and mitotic (S/G2/M) phases fluoresce green (Fig 1D); and, regions of hypoxia are indicated by cyan (Fig 1B, 1C and 1E). Spheroids grown in constant normoxia experience three phases of growth (Fig 1A–1C and 1F) [10, 15, 16]. In phase (i) spheroids grow exponentially as all cells are able to proliferate, indicated by the presence of cells in the S/G2/M phases throughout the tumour spheroid shown by green (Fig 1A). In phase (ii) cells in the central region of the spheroid arrest in G1 phase while cells at the periphery continue to proliferate resulting in inhibited growth (Fig 1B). This arrested region is thought to arise due to spatial differences in nutrient availability, possibly oxygen, and/or a build up of metabolic waste from cells. In phase (iii) the spheroid is characterised by three regions: a central region composed of a necrotic core, 0 < r < Rn(t); an intermediate region of living but proliferation-inhibited cells, Rn(t) < r < Ri(t); and, a region at the periphery composed of living and proliferating cells, Ri(t) < r < Ro(t) (Fig 1C) [10, 13–16]. In comparison to spheroids grown in normoxia, spheroids grown in hypoxia form their necrotic core earlier, the distance from the edge of the spheroid to the hypoxic region and overall size are smaller (Fig 2).

[Figure omitted. See PDF.]

Fig 1. Impact of external environment on the structure of growing tumour spheroids: A focus on oxygen availability.

(a-c) Tumour spheroid growth in standard experimental protocols occurs in three phases. Experimental images shown for FUCCI-transduced human melanoma WM983b spheroids grown in normoxia. (a-b) Experimental images of the equatorial plane of spheroids on Day 2 and 3 after seeding. (c) 3D z-stack representation of half of a spheroid on Day 8 after seeding. Scale bars are 200μm. Colours in (a-c) correspond to cell cycle schematic shown in (d): cells in G1 phase (magenta) and cells in S/G2/M phase (green). Pimonidazole staining reveals the hypoxic regions of spheroid (cyan). (e) Schematic for spherically symmetric spheroid structure representing a quadrant of the equatorial plane of a spheroid. Spheroids in normoxia experience three phases of growth, resulting in a spheroid with three regions at later times: a central region composed of a necrotic core, 0 < r < Rn(t) (black); an intermediate region of living but proliferation-inhibited cells, Rn(t) < r < Ri(t) (magenta); and, a region at the periphery composed of living and proliferating cells, Ri(t) < r < Ro(t) (green). The hypoxic radius, Rp(t) (cyan) satisfies Rn(t) ≤ Rp(t) ≤ Ro(t). (f-j) Schematics for oxygen conditions and time evolution of spheroid structure and overall size in (f) normoxia, (g) hypoxia, (h) deoxygenation experiments, and (i-j) re-oxygenation experiments where spherical symmetry is maintained. Note in (i-j) spheroids transiently undergo the growth phases in reverse. Greyscale shading in (f-j) represent growth phases. (k-l) Schematics for re-oxygenation experiments where we do not observe spherically symmetry due to (k) necrotic core removal and (l) movement of the necrotic core without removal. For additional experimental images for each condition with the boundaries of the proliferating, inhibited and hypoxic regions overlayed see Text A.2 in S1 File.

https://doi.org/10.1371/journal.pcbi.1010833.g001

[Figure omitted. See PDF.]

Fig 2. Tumour spheroid growth in normoxia and hypoxia.

(a,d) Experimental images of the equatorial plane of FUCCI-transduced WM983b spheroids grown in (a) normoxia (21% oxygen) and (d) hypoxia (2% oxygen). (a) Images shown on Day 2, 3, 4, 6, and 8 after seeding. (d) Images shown on Day 2, 4, 6, and 8 after seeding. Scale bars are 200μm. (b,e) Schematics for oxygen partial pressure within spheroids for (b) normoxia and (e) hypoxia, where the solid blue/orange lines represent the oxygen partial pressure inside/outside of the hypoxic region and the hypoxic threshold is determined by the presence pimonidazole staining (blue-dashed). (c,f) Time evolution of Ro(t) (green), Ri(t) (magenta), Rn(t) (black), and Rp(t) (cyan) corresponding to schematic in Fig 1E for spheroids grown in (c) normoxia, and (f) hypoxia. Note that each spheroid measurement is an end-point measurement. For additional experimental images for each condition with the boundaries of the proliferating, inhibited and hypoxic regions overlayed see Text A.2 of S1 File.

https://doi.org/10.1371/journal.pcbi.1010833.g002

To investigate environmental variability we perform additional experiments in time-dependent oxygen conditions. In these experiments we observe various tumour spheroid adaptation mechanisms (Fig 1H–1L). For instance, in re-oxygenation experiments we discover a novel adaptation process where the necrotic core of the spheroid that has formed prior to re-oxygenation moves within the spheroid and in certain situations exits the spheroid as a single object (Fig 1K and S1 Movie). Further, for fifty years tumour spheroid growth has been described by three sequential growth phases but re-oxygenation experiments show that spheroids can transiently experience these phases in the reverse order (Fig 1I and 1J). Other observations from these experiments agree with intuitive expectations, but have not previously been explored nor quantified.

Many mathematical models describing tumour growth have been proposed and developed over the past fifty years, reviewed in [31–36]. Here we summarise studies that are of direct relevance to this study. First, we are interested in estimating the oxygen partial pressure profile within each individual spheroid at each time point. Grimes et al. [37] approach this problem using a mathematical model describing diffusion and consumption of oxygen. Using their oxygen diffusion model they estimate and validate the oxygen partial pressure within a spheroid, estimate the outer radius when the necrotic region forms, and the oxygen consumption rate. Gomes et al. [17] subsequently estimate the oxygen partial pressure within spheroids grown in normoxia (21% oxygen) and 5% oxygen which they refer to as physioxia. Their spheroid experiments, performed with HCT113 colon adenocarcinoma cells, show that spheroids grown in physioxia are smaller in overall size than spheroids grown in normoxia. Furthermore, Gomes et al. [17] measure the distance from the spheroid edge where 50% of the cells are proliferating and regions of hypoxia. They conclude that oxygen partial pressure is a rate-limiting parameter for cell proliferation. Here we briefly follow the approach of Grimes et al. [37] to analyse spheroid snapshots and then extend this analysis by exploring the temporal evolution of spheroid structure.

Throughout we quantitatively analyse the temporal evolution of the overall size and structure of tumour spheroids using mathematical modelling and statistical uncertainty quantification, including profile likelihood analysis and Bayesian inference approaches. We start with the seminal Greenspan mathematical model [10, 15, 16, 34]. Greenspan’s model describes the three phases of growth and is relatively simple in comparison to other models [31–36]. This simplicity is a great advantage. Each mechanism and parameter in Greenspan’s model has a biologically meaningful interpretation. Further, these parameters can be identified and estimated with the experimental data that we collect in this study [15, 16], which is not the case for other more complicated mathematical models [38–40]. In addition, here we show that we can extend Greenspan’s model to analyse growth in time-dependent external oxygen partial pressures while retaining physical and biologically insightful interpretations of results. Recent work by Lewin et al. [41] has extended Greenspan’s model in a comparable manner. However, in their theoretical study the external environmental conditions are constant throughout time and rapid changes in the oxygen partial pressure within the spheroid and to spheroid structure are driven by a rapid loss in living cell volume due to radiation induced cell death. With our new models we reveal biological adaptation mechanisms in de-oxygenation and re-oxygenation experiments and estimate key biological parameters.

In the following we first analyse spheroid experiments in normoxia and hypoxia. Such experiments demonstrate that Greenspan’s model describes the experimental data remarkably well. Further, that oxygen mechanisms accurately describe the growth and formation of the necrotic core, whereas other mechanisms, possibly waste mechanisms, likely result in the growth and formation of the inhibited region. We then extend the mathematical model to interpret deoxygenation and re-oxygenation experiments, providing quantitative insights to biological adaption mechanisms throughout. We conclude by describing unexpected behaviours observed in re-oxygenation experiments.

Results

Here we focus on WM983b spheroids in normoxia, hypoxia and deoxygenation experiments. Similar results for WM793b and WM164 spheroids are discussed in Text E and Text F in S1 File. For re-oxygenation experiments, we compare results from WM983b and WM793b spheroids as we observe a range of behaviours. In Text F of S1 File we discuss additional WM164 re-oxygenation experiments.

Oxygen diffusion alone is insufficient to describe spheroid growth

We capture end-point equatorial plane images for spheroids grown in normoxia and hypoxia measuring Ro(t), Ri(t), Rn(t), and Rp(t) (Fig 2A, 2C, 2D and 2F) (Methods: Image processing). These measurements are remarkably consistent within each condition and time point (Fig 2C and 2F). Comparing spheroids grown in normoxia and hypoxia, we observe vastly different tumour growth dynamics and internal structure (ξn(t) = Rn(t)/Ro(t), ξi(t) = Ri(t)/Ro(t), ξp(t) = Rp(t)/Ro(t)), even when comparing spheroids of similar size (Fig 2C and 2F, and Fig U in S1 File).

For deeper mechanistic insight we use mathematical modelling and statistical uncertainty quantification to interpret our observations. Specifically, we use Greenspan’s mathematical model [10] and show that it accurately describes spheroid growth in both normoxia and hypoxia. Then using parameter estimation we identify biological mechanisms that differ between normoxia and hypoxia. Key model details are now discussed, for further details see Methods Greenspan’s mathematical model and Text C.1 of S1 File. The model assumes each spheroid is spherically symmetric and maintained by cell-cell adhesion or surface tension. The independent variables are time t [days], and radial position, r [μm]. Conservation of volume gives an equation describing the time evolution of the outer radius, Ro(t) [μm],(1)where s [day−1] is the rate at which cell volume is produced by mitosis per unit volume of living cells, and λ [day−1] is the proportionality constant describing the rate at which cell volume is lost from the necrotic core. In these experiments, Eq (1) simplifies as Ri(t) ≥ Rn(t). This restricts our attention to two interpretations of the model that differ with respect to how Ri(t) is defined. In the following discussion, we refer to these interpretations as hypotheses and show that hypothesis 2, where oxygen mechanisms drive Rn(t) and waste mechanisms drive Ri(t), is more consistent with the spheroids considered in this study.

In hypothesis 1 (Fig 3A), oxygen diffuses with diffusivity, k [m2 s−1], and is consumed by living cells at a constant rate, α [m3 kg−1 s−1]. The external oxygen partial pressure is p∞ [%]. Oxygen diffusion is fast relative to the growth of the spheroid, so that the oxygen partial pressure within the spheroid, p(r(t)) for 0 ≤ r ≤ Ro(t), is governed by(2)where H(⋅) is the Heaviside function and Ω [mmHg kg m−3] is a conversion constant from volume of oxygen per unit tumour mass to oxygen partial pressure [37]. The inhibited radius, Ri(t), is implicitly defined by p(r(t)) = pi [%] provided the oxygen partial pressure is sufficiently large (Fig 3A), and Ri(t) = 0 otherwise. Note we assume that pi [%] is a constant oxygen partial pressure independent of the external oxygen partial pressure.

[Figure omitted. See PDF.]

Fig 3. Mechanisms governing tumour spheroid growth in normoxia and hypoxia.

(a,b) Schematics for two hypotheses. (a) Hypothesis 1: oxygen diffusion and consumption drives formation and growth of necrotic and inhibited regions. Rn(t) [μm] and Ri(t) [μm] implicitly defined via oxygen partial pressure thresholds pn [%] and pi [%], respectively. (b) Hypothesis 2: oxygen diffusion and consumption drives formation and growth of necrotic region, waste production and diffusion drives formation and growth of inhibited region. Rn(t) [μm] implicitly defined via oxygen partial pressure thresholds pn [%]. Ri(t) [μm] implicitly defined via waste threshold βi [mol μm−3]. (c-h) Analysis of spheroid snapshots assuming hypothesis 1. Box charts representing estimates of (c) outer radius when necrotic region forms, Rc [μm], (d) oxygen consumption rate, α [m3 kg−1 s−1], and (f) pi [%] for spheroids grown in normoxia and hypoxia. (e) Comparison of measured and predicted Rp(t) [μm] (R2 = 0.751) suggesting the estimated oxygen partial pressure within each spheroid is reasonable. (g,h) Estimated oxygen partial pressure within (g) a spheroid grown in normoxia and (h) a spheroid grown in hypoxia. Insets in (g,h) show FUCCI signal and pimonidazole staining from Day 8 and 6, respectively, with scale bar is 200 μm. In (g,h) horizontal red and blue dashed lines correspond to estimates of pi [%] from spheroid grown in normoxia and hypoxia, respectively. The estimate of pi from normoxia in (g) suggests the spheroid growing in hypoxia in (i) should be entirely arrested under hypothesis 1. (i-o) Analysis of temporal evolution of spheroid structure using Greenspans mathematical model and hypothesis 2. (i-j) Comparison of experimental data (box charts) with Greenspan’s mathematical model (solid lines) simulated with the maximum likelihood estimates of parameters for (i) normoxia and (j) hypoxia. In (i,j) time evolution of outer radius, Ro(t) (green), inhibited radius, Ri(t) (magenta), and necrotic radius, Rn(t) (black). Data represent an average of twelve spheroids on days 2, 3, 4, 6, and 8 for spheroids grown in normoxia and an average of seven spheroids on days 2, 4, 6, and 8 for spheroids grown in hypoxia (Methods: Experimental methods). (k-o) Profile likelihoods for parameters of Greenspans mathematical model: (k) initial outer radius, Ro [μm], (l) proliferation rate, s [day−1], (m) outer radius when necrotic region first forms, Rc [μm], (n) dimensionless parameter relating proliferation rate and mass loss from necrotic core, γ [-], and (o) dimensionless parameter relating oxygen and waste mechanisms, Q [-] (Methods: Greenspan’s mathematical model).

https://doi.org/10.1371/journal.pcbi.1010833.g003

In hypothesis 2 (Fig 3B), diffusible metabolic waste is produced by living cells at a constant rate per unit volume, P [mol μm−3 day−1], and diffuses with diffusivity κ [μm2 day−1]. Waste diffusion is fast relative to the growth of the spheroid, so the waste concentration within the spheroid, β(r(t)) [mol μm−3], is governed by(3)The inhibited radius, Ri(t), is implicitly defined through β(r(t)) = βi [mol μm−3] provided the waste concentration is sufficiently large (Fig 3B), and Ri(t) = 0 otherwise.

Both hypothesis 1 and 2 assume that Rn(t) is implicitly defined by p(Rn(t)) = pn provided the oxygen partial pressure is sufficiently small, and Rn(t) = 0 otherwise. Informed by experimental results [37], we set pn = 0 [%].

Analysis of the model provides an analytical expression for the time when the inhibited region forms (Eq 7.1). For hypothesis 2, the inhibited region is predicted to form at the same time independent of oxygen dynamics, provided the spheroids grown in normoxia and hypoxia are initially the same size. This appears consistent with results in Fig 2A–2D where the inhibited region has formed on day 2 for both conditions. In contrast, with hypothesis 1 the time to form the inhibited region depends on oxygen mechanisms and specifically p∞, but without knowledge of additional parameters further insights are unclear.

By incorporating statistical uncertainty quantification to estimate parameters of the model we gain further mechanistic insights. A key assumption common to hypothesis 1 and 2 is that oxygen diffusion and consumption drives the time evolution of Rn(t). We test this assumption directly by analysing radial measurements of each spheroid at each time point independently [37]. Using measurements of Ro(t) and Rn(t) we estimate: the outer radius when the necrotic region first forms, Rc; α; and Rp(t) (Fig 3C and 3D, and Text D.1.1 in S1 File). Image processing to measure Rp(t) is more challenging than for Ro(t) due to gradients in the pimonidazole signal (Text B in S1 File). However, we find good agreement between experimentally measured and predicted values of Rp(t) > (R2 = 0.751, Fig 3E). This approach allows us to estimate the oxygen partial pressure within each spheroid at each time point (Fig 3G and 3H). From the results in Fig 3E, 3G and 3H, and Fig V in S1 File we conclude that oxygen of diffusion alone is a reasonable and sufficient mechanism to describe the formation of the necrotic core in WM983b spheroids [37]. Furthermore, these results suggest that the following assumptions are reasonable: pn = 0; a constant oxygen consumption rate within the spheroid; spherical symmetry; oxygen at the edge of the spheroid can be approximated with the oxygen settings on the incubator.

Given that Rn(t) is reasonably described by oxygen mechanisms we now examine hypothesis 1. Hypothesis 1 assumes that oxygen mechanisms alone drive the time evolution of Ri(t). We then estimate pi using the estimated oxygen partial pressure within each spheroid, p(r(t)) for 0 < r < Ro(t), measurements of Ri(t) and the definition p(Ri(t)) = pi (Fig 3A, 3G and 3H). Estimates of pi are consistently larger for spheroids grown in normoxia compared to spheroids grown in hypoxia (Fig 3F), consistent with results for other cell lines [17]. This is inconsistent with hypothesis 1. Specifically, results from normoxia suggest that spheroids grown in hypoxia should have larger inhibited regions than experimentally measured (Fig 3G and 3H). Similarly, results from spheroids grown in hypoxia suggest that spheroids grown in normoxia should have smaller inhibited regions than experimentally measured. Further, we do not restrict ourselves to the notion that oxygen must be directly responsible for formation and growth of the inhibited region as it is unclear why cells within spheroids grown in normoxia arrest at a higher oxygen partial pressure than cells within spheroids grown in hypoxia. These results provide strong evidence to suggest that oxygen alone is insufficient to describe the formation of the inhibited region across multiple oxygen conditions for this cell line.

To test hypothesis 2, which assumes that waste mechanisms drive the time evolution of Ri(t), we first analyse measurements of each spheroid at each time point independently. Experimentally measuring waste within spheroids is challenging, so we estimate the waste concentration within each spheroid and use measurements of Ri(t) to estimate the outer radius when the inhibited region first forms, (Fig W in S1 File). We observe that is larger for spheroids grown in normoxia than hypoxia, which may be due to changes in βi or P or κ (Fig W in S1 File). These results do not provide sufficient evidence to reject hypothesis 2.

To test whether hypothesis 2 is reasonable we estimate model parameters for spheroids grown in normoxia and hypoxia. Specifically, we estimate the five key parameters: Θg = (Ro (0), s, Rc, γ, Q), where γ = λ/s [-] and Q [-] are dimensionless quantities (Methods: Parameter estimation and identifiability analysis). Simulating the model with the maximum likelihood estimate (MLE), , for normoxia shows good agreement with the experimental data from spheroids grown in normoxia (Fig 3I). Similarly, simulating the model at for hypoxia shows good agreement with the experimental data from spheroids grown in hypoxia (Fig 3J). These results suggest the model accurately captures the dynamics of tumour spheroid growth.

Alongside the point estimates and , we are interested in forming approximate 95% confidence intervals for each of the parameters. To perform this analysis we employ profile likelihood analysis (Methods: Parameter estimation and identifiability analysis). All profile likelihoods computed here are narrow and well-formed around a single central peak, the MLE, indicating that parameters are identifiable and that a relatively narrow range of parameters give a similar match to the data as the MLE (Fig 3K–3O) [16]. Approximate 95% confidence intervals are obtained from these profile likelihoods for each parameter. The profile likelihoods for the initial outer radius, Ro(0), do not overlap and agree with observations that Ro(0) is smaller for spheroids grown in hypoxia than normoxia (Fig 3K). The profile likelihood for s interestingly estimates that the rate of cell proliferation per unit volume is faster in hypoxia than normoxia (Fig 3L). This result may seem surprising as a simplistic assumption would be that less oxygen results in less proliferation. However, our result is consistent with observations from other cell lines where an intermediate level of hypoxia encourages more proliferation than normoxia [42–44]. Profile likelihoods for Rc are consistent with estimates of Rc obtained by analysing spheroid measurements independently, a good consistency check for the two methods (Fig 3C and 3M). The other profile likelihoods for γ and Q overlap suggesting that these parameters are consistent across normoxic and hypoxic conditions (Fig 3N and 3O). Posterior densities and prediction intervals, estimated using Bayesian inference, also show good agreement with results here and the experimental data (Fig X in S1 File). Similar results hold for the other two cell lines (Text E and Text F in S1 File).

In the remainder of this study, we take the most fundamental approach and proceed with Greenspan’s mathematical model and interpret the governing mechanisms with hypothesis 2. We note that other biological mechanisms may also be relevant. However, as the model already appears to capture the key dynamics (Fig 3I and 3J) we will avoid overcomplicating the model at this stage. Furthermore, in the following analysis of deoxygenation and re-oxygenation experiments we necessarily extend the model.

Adaptation to deoxygenation

The mechanisms underlying how tumour spheroids adapt to time-dependent external environments is unclear. Here, we perform a series of deoxygenation experiments, where spheroids grown in normoxia are transferred to hypoxic conditions at t = ts [days]. Analysing spheroid snapshots reveals how spheroids, and in particular their internal structure, adapt. Extending Greenspan’s mathematical model and using parameter estimation, we identify and quantify key biological adaptation mechanisms.

In the deoxygenation experiments we set p∞ = 21 [%] for 0 < t ≤ ts [days], p∞ = 2 [%] for ts < t < 8 [days], and ts = 2 [days] (Fig 4A). At ts = 2 [days] all spheroids are in phase (i) of growth with proliferating cells throughout and no inhibited or necrotic region (Day 2 of Fig 4B). At ts + 1 [days] the FUCCI signal in the central region of the spheroid is blurred, relative to the signal at the periphery, indicating dying and dead cells (Day 3 of Fig 4B and Fig T in S1 File). Therefore, we identify this central region as the necrotic core (Text B in S1 File). Experimental images at later times show that Rn(t), Ri(t) and Rp(t) continue to increase but at a much slower rate (Days 4–8 of Fig 4B and 4C). Throughout the experiment Ro(t) remains approximately constant (Days 2–8 of Fig 4B and 4C) confirming that the most important changes involve the internal structure and not the overall spheroid size. Further, ξn(t), ξi(t), and ξp(t) approach values observed at late times for spheroids grown in hypoxia (Fig 4C and Fig Y in S1 File).

[Figure omitted. See PDF.]

Fig 4. Analysis of deoxygenation experiments reveals tumour spheroid adaptation mechanisms.

(a) Schematic for deoxygenation experiment where the external oxygen partial pressure switches from normoxia to hypoxia at ts = 2 [days]. (b) Experimental images of the equatorial plane of WM983b spheroids on Days 2, 3, 4, 6, and 8 after seeding. Scale bar is 200 μm. Colours in (b) correspond to cell cycle schematic shown in Fig 1D and hypoxic regions are shown by pimonidazole staining (cyan). Necrotic region forms on Day 3 (blurry central region) and grows (dark central region on Days 4, 6, and 8). (c) Comparison of experimental data (box charts) with deoxygenation mathematical model (solid line) simulated with the maximum likelihood estimates of parameters. Time evolution of outer radius Ro(t) (green), inhibited radius Ri(t) (magenta), hypoxic radius Rp(t) (cyan), and necrotic radius Rn(t) (black). Blue vertical dashed lines in (b-c) indicate ts. (d) Oxygen partial pressure within a spheroid estimated at ts under normoxia (red) and hypoxia (blue). (e-f) Schematics for spheroid structure in: (e) normoxia at t = ts; and, (f) hypoxia immediately after ts. In (f) the predicted necrotic radius, , is greater than the actual necrotic radius, Rn(ts), due to the reduced external oxygen partial pressure. Cells in region die due to lack of oxygen and increase the volume of the necrotic core. (g-i) Analysis of spheroid snapshots in deoxygenation experiment. (g) Estimates of oxygen consumption rate, α [m3 kg−1 s−1], from analysis of spheroid snapshots (orange squares) compared to prediction intervals generated from the deoxygenation mathematical model. (h) Estimates of the outer radius when the inhibited region forms, [μm], from analysis of spheroid snapshots (orange squares) compared to prediction intervals generated from the mathematical model. In (g-h) deoxygenation model simulated with the posterior means of the parameters (black), 50% posterior region of prediction interval (dark grey), and 97.5% posterior region of prediction interval (light grey). (i) Posterior density estimates for the timescales of adaptation of the oxygen consumption rate and outer radius when the inhibited region forms denoted by τα (purple) and (red).

https://doi.org/10.1371/journal.pcbi.1010833.g004

To interpret these deoxygenation experiments we extend Greenspan’s mathematical model [10]. We assume that the change in p∞ at ts is instantaneous, which is reasonable since the switch from normoxia to hypoxia requires only 1–2 minutes when the spheroids are transferred between incubators. This time is very short in comparison to the duration of the experiment and time interval between data points. Similarly, we assume that the oxygen partial pressure within the spheroid adapts to the change in p∞ instantaneously, which is reasonable since oxygen takes approximately 10 seconds to diffuse across a distance of 100 μm [10]. Then we estimate the oxygen partial pressure within the spheroid at ts under normoxic and hypoxic conditions (Methods: Mathematical model to interpret deoxygenation experiments, Fig 4D). Immediately after ts the predicted necrotic radius, denoted and implicitly defined by , is greater than the actual necrotic radius, Rn(t), specifically (Fig 4D–4F).

Before considering the region , recall that parameter estimates from spheroids grown in normoxia and hypoxia differ. Specifically, α (Fig 3D), λ = γs (Fig 3L and 3N), s (Fig 3L), and (Fig W in S1 File) are all different. Therefore, we expect that these parameter values will evolve in time after ts. To account for such changes we define the following, for t ≥ ts,(4.1) (4.2) (4.3)(4.4)where τα [days], τλ [days], τs [days], and [days] denote timescales of adaptation for α, λ, s, and , respectively. Further, the new constants in Eq (4) with subscripts n and h, for example αh and αn, represent parameter estimates from spheroids grown in normoxia and hypoxia, respectively. The other parameters (k, Ω, κ), are assumed to be constants. Hence, [μm2], [-], and γ(t) = λ(t)/s(t) [-] are functions of time.

In the region we assume cells die and increase the size of the necrotic core at a rate [day−1] per unit volume (Fig 4F). Conservation of volume for the necrotic core at time t, Vn(t), gives (Text C.2.1 in S1 File)(5)Volume is converted to radius for comparison with experimental data using . At later times the term involving dominates the right hand side of Eq (5) and Rn(t) tends to . Eq (1), obtained by conservation of volume arguments, remains valid by including the time dependence in s(t) and λ(t). At ts there is no immediate change in the waste concentration within the spheroid and so no immediate change to Ri(t). However, as changes with time (Eq (4.4)) the waste concentration within the spheroid and Ri(t) evolve over time in part directly due to deoxygenation.

In this new mathematical model (Eqs (8.1)–(8.11)) there are fifteen parameters Θd = (Ro(0), αn, αh, τα, , , , sn, sh, τs, λn,λh, τλ, , . However, using Bayesian inference for parameter estimation the biological adaptation mechanisms become clearer (Methods: Parameter estimation and identifiability analysis). We identify fast adaptation to deoxygenation in α(t) with mean τα = 0.23 [days] and for with mean [days] (Fig 4G–4J). Furthermore, simulating the new deoxygenation mathematical model we find good agreement with the experimental measurements of Ro(t), Rn(t), and Ri(t) (Fig 4C and Fig Z in S1 File). Therefore, our new mathematical model provides a mechanistic description to the observed growth dynamics in the variable external environment and appears to capture key adaptation mechanisms.

Adaptation to re-oxygenation

We also perform re-oxygenation experiments, where spheroids grown in normoxia are transferred to hypoxic conditions at time ts. These re-oxygenation experiments exhibit a range of unexpected biological adaptation mechanisms for each cell line that appear to depend on: ts; spheroid size at re-oxygenation, Ro(ts); and necrotic core fraction at re-oxygenation, ξn(ts) = Rn(ts)/Ro(ts).

First we focus on slower growing WM793b spheroids and ts = 2 [days] (Fig 5A and 5C). In hypoxic conditions prior to re-oxygenation, experimental images show a large hypoxic region (Day 2 in Fig 5A, 5C and 5G). However, after re-oxygenation at ts + 1 [days] there is no hypoxic region (Day 3 in Fig 5C and 5G). Spheroid growth after deoxygenation appears to progress similar to spheroids that are grown in normoxia throughout (Days 3–8 in Fig 5C and 5G).

[Figure omitted. See PDF.]

Fig 5. Tumour spheroids exhibit a range of adaption mechanisms in response to re-oxygenation.

(a,b) Schematics for re-oxygenation experiments, where the external oxygen environment switches from hypoxia to normoxia at (a) ts = 2 [days] and (b) ts = 4 [days]. (c,g) Hypoxia recovery and growth for WM793b cell line with ts = 2 [days] with (c) experimental images and (g) radial measurements. (d,h) Necrotic core infilling and growth for WM793b cell line with ts = 4 [days] with (d) experimental images and (h) radial measurements. Green dashed line in (c,d,g,h,i,j) indicate ts. Colours in (g,h) Ro(t) (green), Ri(t) (magenta), Rn(t) (black), and Rp(t) (cyan). (e,f) Schematics for tumour spheroid structure due to re-oxygenation at ts. (i-l) Removal of necrotic core and growth for WM983b cell line. (i) Confocal images for experiment with ts = 2 [days]. (j) Confocal images for experiment with ts = 4 [days]. (k) Bright-field images for experiment with ts = 2.5 [days]. (l) Schematic for removal of necrotic core and growth in WM983b cell line. (m) Bright-field images for experiment with ts = 5.5 [days]. (n) Schematic for movement of necrotic core and growth in WM983b cell line. Colours in (c,d,i,j) correspond to cell cycle schematic shown in Fig 1D and hypoxic regions are shown by pimonidazole staining (cyan). Scale bars in (c,d,i,j,k,m) are 200 μm.

https://doi.org/10.1371/journal.pcbi.1010833.g005

For WM793b spheroids and ts = 4 [days], a necrotic core forms before ts (Day 4 in Fig 5D and 5H). However, after deoxygenation, at ts + 2 [days] there is no necrotic core (Day 6 in Fig 5D and 5H). At ts + 4 [days] spheroids are either in phase (i) or phase (ii) (Day 8 in Fig 5D and 5H). While traditional tumour spheroid experiments progress through phase (i), (ii), and (iii) sequentially, these experiments show that spheroids can transition transiently through the growth phases in reverse order before subsequently growing in the usual order.

To interpret the WM793b re-oxygenation experiments we proceed analogously to the deoxygenation experiments. We extend Greenspan’s mathematical model to account for differences in parameter estimates between normoxia and hypoxia. Estimating the oxygen partial pressure within the spheroid at ts, the predicted necrotic core, denoted , is smaller than the actual necrotic core, Rn(ts), provided Rn(ts) > 0. In the region where the necrotic core is now supplied with oxygen, we assume that size of the necrotic core decreases at a rate [day−1] per unit volume. We assume that a fraction, 0 ≤ ν ≤ 1, of the volume lost from the necrotic core recovers from the harsh oxygen conditions and increases the population of living cells, and the remaining volume lost from the necrotic core diffuses out of the spheroid and does not influence Ri(t). Using Bayesian inference we estimate the seventeen model parameters, Θr = (Ro(0), Rn(0), αn, αh, τα, , , , sn, sh, τs, λn, λh, τλ, , , ν) and simulate the mathematical model. We observe good agreement with the experimental data suggesting our new re-oxygenation mathematical model captures key mechanisms underlying adaptation and growth (Fig 5G and 5H, and Figs BB and CC in S1 File).

Results for WM983b spheroids are unexpected. We hypothesised, based on the exploration of the mathematical model, that spheroid growth dynamics may occur in reverse as observed for WM793b spheroids. However, we did not anticipate that WM983b spheroids would lose their symmetrical internal structure and necrotic core due to re-oxygenation (Fig 5I–5N, S1 and S2 Movies). The WM983b experiments are performed at four different re-oxygenation times ts = 2, 2.5, 4, and 5.5 [days] (Fig 5I, 5J, 5K and 5M). For all experiments, spheroids at ts prior to re-oxygenation are in phase (iii). For ts = 2, confocal microscopy reveals that at ts + 1 [days] there is a necrotic region but it is not at the centre of the spheroid (Day 3 of Fig 5I). At later times there is no necrotic region and growth proceeds analogous to spheroids in normoxia. Similarly, for ts = 4 (Fig 5J).

To explore this unusual behaviour for ts = 2.5 [days] and ts = 5.5 [days] we perform experiments using the IncuCyte S3 live cell imaging system (Sartorius, Goettingen, Germany) and obtain hourly bright-field images after re-oxygenation. In Fig 5K and S1 Movie with ts = 2.5 [days], initially the necrotic core of the spheroid is visible as a dark central region. At later times the necrotic core is located closer to the edge of the spheroid and the radially symmetric internal structure is lost (12, 24 and 36 hours after re-oxygenation in Fig 5K). At ts + 2 [days] the necrotic core appears to have exited the spheroid as a single object (48 hours after re-oxygenation in Fig 5K). Tracking the position of the necrotic core relative to the spheroid suggests the necrotic core moves randomly (Text D.3.1 in S1 File).

We observe similar behaviour for WM983b spheroids with ts = 5.5 [days] (Fig 5M and S2 Movie). However, likely due to the larger Ro(ts) and ξn(ts) here, the necrotic core is close to the edge of the spheroid but does not exit as a single object (1.5 and 3 days after ts in Fig 5M). Instead as the spheroid grows necrotic matter forms at the centre of the spheroid and appears to merge with the necrotic matter located closer to the periphery (3, 4.5, and 6 days after ts in Fig 5M). As the WM983b spheroids do not maintain spherical symmetry we do not interpret these experimental data with the re-oxygenation mathematical model. Schematics describing the behaviours are presented in Fig 5L and 5N.

Confocal images for WM164 spheroids suggest that symmetry is maintained (Figs P and Q in S1 File), and that WM164 spheroids transition transiently through the growth phases in reverse order before subsequently growing in the usual order. In contrast to WM793b spheroids, WM164 spheroids can resume growth while maintaining a necrotic core (Fig CC in S1 File). These different behaviours appear to be dependent on ts, Ro(ts) and ξn(ts). To interpret WM164 re-oxygenation experiments we use the re-oxygenation model (Text F in S1 File). However, additional care should be exercised interpreting WM164 re-oxygenation results. Brightfield time-lapse images show that mass from the necrotic core can move to the periphery and exit the spheroid (Fig R in S1 File and S3 Movie).

Discussion

Tumours grow in complicated fluctuating external environments. However, spheroid experiments used to study tumours are typically performed in constant external environments and with oxygen partial pressures that are much greater than in vivo. To explore this gap between in vitro and in vivo conditions, we analyse a series of tumour spheroid experiments using mathematical modelling and statistical uncertainty quantification. Growing spheroids in time-dependent external oxygen conditions reveals a range of behaviours not observed in standard experimental protocols. For fifty years, tumour spheroid growth has been characterised by three sequential growth phases: phase (i) exponential growth; phase (ii) reduced exponential growth; and phase (iii) saturation. However, here in re-oxygenation experiments, spheroids can transiently undergo these growth phases in reverse. Furthermore, spheroids can lose their spherically symmetric structure and necrotic core. Deoxygenation experiments also show that large changes to the internal structure of spheroids can occur while the overall size remains constant. Overall, our results suggest that oxygen and internal structure play pivotal roles in spheroid growth and should be taken into account when interpreting spheroid experiments. This is important as many studies do not provide sufficient information to replicate oxygen conditions and do not measure spheroid internal structures.

Tumour spheroid growth is a complex process involving multiple mechanisms. However, the contribution of each mechanism to growth is unclear using experimentation alone. To quantitatively explore which mechanisms contribute to spheroid growth we use the seminal Greenspan mathematical model [10]. The model describes the growth of spheroids in normoxia and hypoxia remarkably well. Moreover, our analysis suggests that growth and formation of the necrotic core is reasonably described by oxygen diffusion and consumption, whereas the growth and formation of the inhibited region is more accurately described by waste production and diffusion. Using statistical uncertainty quantification we show that the rates at which different biological processes occur differ between normoxia and hypoxia. Therefore, external environmental conditions should be taken into account when interpreting tumour spheroid experiments. Previous studies analysing previously available experimental data with Greenspan’s model [15, 16] have not been able to distinguish between the two mechanisms. Further, by considering the temporal evolution of spheroid structure, our results build on studies analysing spheroid snapshots only [17, 37]. As the model captures the key dynamics, we do not include other biological mechanisms that may be relevant in future studies, for example glucose [4]. Introducing additional mechanisms prior to developing the deoxygenation and re-oxygenation models would complicate the analysis, and likely result in parameters being non-identifiable and not physically interpretable. Both of which we aim to avoid. Many mathematical models have been developed with additional mechanisms, but they have not been quantitatively tested with experimental data [25, 31–34]. The experimental data and framework that we provide here are suitable to test such models.

Deoxygenation and re-oxygenation experiments reveal how spheroid overall size and internal structure adapt to changes in the external environment. However, without a mathematical modelling and statistical uncertainty quantification framework like we use here, the mechanisms underlying adaptation are challenging to identify and interpret. Extending Greenspan’s mathematical model allows us to interpret, analyse, and describe these deoxygenation and re-oxygenation experiments remarkably well. Parameter estimation and identifiability analysis, using profile likelihood analysis and Bayesian inference, allows us to identify and quantify the contributions of key biological mechanisms to adaptation and growth. For both models, the analysis identifies a narrow range of behaviours that differ in terms of the rates of adaptation in ts < t < ts + 1 [days] and long term dynamics (Text D.2.1 in S1 File). These modelling predictions raise interesting questions (Text D.2.1 in S1 File). In comparison to standard experimental protocols, we collect more data per time point and over a longer duration. Further, we build on previous studies [13, 15, 16] to improve on standard experimental designs by measuring the internal structure of spheroids and hypoxic regions in addition to the overall size. Even these improvements to standard protocols are insufficient to identify all adaptation processes. As with all studies, additional data would be beneficial. In particular, frequent measurements of oxygen and internal structure at early times would be useful but are challenging to obtain experimentally (Text D.2 in S1 File).

Re-oxygenation experiments reveal unexpected necrotic core removal in WM983b spheroids. To the best of our knowledge this behaviour has not been previously described. The exact mechanisms underlying this behaviour are unclear. We hypothesise that small asymmetry at re-oxygenation, possibly in the distribution of proliferating cells, in combination with changes to cell-cell adhesion and physical interactions contribute. We also observe that mass from the necrotic core of WM164 spheroids can move to the periphery and exit the spheroid. It is unclear whether these observations are relevant to other cell lines and to in situ tumours. Interesting future work is to explore these unusual behaviours in greater detail experimentally. For example, can this phenomenon be induced by other external environmental changes, drug treatments, and in vivo. As the necrotic core removal destroys the spherical symmetry of the spheroid we do not use Greenspan’s mathematical model to analyse these observations. It would be interesting future work to ascertain the physical and biological mechanisms that could explain these observations, possibly by extending a multiphase model [45, 46], studies on cellular blebbing [47], compound droplet models [48], or models with different treatments of necrotic material [46].

This work lays the foundation for further studies bridging the gap between clinical conditions and standard experimental protocols. Here, we consider normoxia and hypoxia and switches between normoxia and hypoxia. The impact of other oxygen conditions on spheroid growth and cell cycle progression within growing spheroids are also worth consideration, for example to mimic in vivo oxygen gradients [49], in vivo vascularisation, disrupted oxygen supplies, or cyclic hypoxia [1, 22]. Microfluidic devices may be one useful approach [25], but challenges visualising the internal structure of spheroids throughout such experiments must be overcome. Intracellular responses to oxygen changes are also of interest [50, 51]. Results in this study also depend upon Greenspan’s mathematical model, and therefore assumptions of spherical symmetry. Interesting future work would be to explore relaxing spherical symmetry assumptions. Further interesting work would be to repeat the experiments using the EF5 hypoxia marker and compare results, repeat the experiments with different cell lines, explore similarities and differences with in situ tumours, and explore alternative mechanisms that might explain why estimates of pi differ across oxygen conditions. Alternative mathematical models to explain these experiments and provide further useful insights and predictions may also be of interest. For example, one could explore replacing the time-dependent adaptation mechanisms with functional forms incorporating additional biological mechanisms. When developing such alternative models one should take into account challenges with parameter identifiability that can arise as model complexity increases [38]. We also note that our framework is well suited to explore the role of other changing external conditions on spheroid growth, for example nutrient availability and mechanical confinement [4, 11], and can be extended to explore different treatment strategies, for example radiotherapy and chemotherapy [34, 41].

Materials and methods

Mathematical modelling

Greenspan’s mathematical model.

Key elements of Greenspan’s mathematical model for fixed p∞ are included in the main text. Further details of the model derivation are included in Text C.1 in S1 File. Recall that this model has two interpretations, that we refer to as hypotheses 1 and 2.

* Hypothesis 1 assumes that the necrotic and inhibited regions are both driven by oxygen diffusion and consumption.

* Hypothesis 2 assumes that the necrotic region is driven by oxygen diffusion and consumption whereas the inhibited region is driven by waste production and diffusion.

Here, we present the governing differential-algebraic system of equations for the outer radius, Ro(t), necrotic radius, Rn(t), and inhibited radius, Ri(t), for both hypotheses 1 and 2,(6.1)(6.2)(6.3)In these experiments, Eq (6.1) simplifies as Ri(t) ≥ Rn(t), consistent with our parameter choices (Methods: Parameter estimation and identifiability analysis) [16]. For both hypothesis 1 and 2: the outer radius when the necrotic region forms is Rc = [6kp∞/(αΩ)]1/2; Eq (6.1) arises from conservation of volume; and, Eq (6.2) is obtained by evaluating the oxygen partial pressure within the spheroid at the necrotic threshold. For hypothesis 1, Eq (6.3) is obtained by evaluating the oxygen partial pressure within the spheroid at the oxygen inhibited threshold, pi, and Q2 = (p∞ − pi)/p∞ so the left hand side of Eq (6.3) is . In contrast, for hypothesis 2, Eq (6.3) arises by evaluating the waste concentration within the spheroid at the waste inhibited threshold, βi, and Q2 = βiκαΩ/(Pkp∞) so the left hand side of Eq (6.3) is . We solve the system of Eqs (6.1)–(6.3) numerically (Text C.1.2 in S1 File).

Analysing the model, the inhibited region forms at [10](7.1)For hypothesis 1, Eq (7.1) is t = (3/s) log ([6k(p∞ − pi)/α]1/2/Ro(0)). For hypothesis, 2 Eq (7.1) is t = (3/s) log ([6βiκ/P]1/2/Ro(0)).

Mathematical model to interpret deoxygenation experiments.

Key elements of the mathematical model derived to interpret deoxygenation experiments are included in the main text. Further details of the model derivation are included in Text C.2.2 in S1 File. Here, we present the governing equations for 0 < t < ts and t > ts.

For 0 < t < ts we solve Greenspan’s mathematical model [10] in normoxia interpreting the governing mechanisms with hypothesis 2. The differential-algebraic system of Eqs (6.1)–(6.3) are solved to determine Ro(t), Rn(t), and Ri(t).

After deoxygenation, t > ts, we extend Greenspan’s mathematical model to account for adaptation to hypoxia. Rewriting Eqs (4), (5) and (6.1), and solving Eq (6.2) for instead of Rn(t), gives the governing differential-algebraic system of equations(8.1)(8.2) (8.3) (8.4) (8.5) (8.6) (8.7)(8.8) (8.9) (8.10) (8.11)Note that in the long time limit t → ∞, we recover Greenspan’s mathematical model for normoxia (Eqs (6.1)–(6.3)). Specifically, α(t) → αh, λ(t) → λh, s(t) → sh and as t → ∞. Further, the term involving dominates the right hand side of Eq (8.2) as t → ∞, so as t → ∞. We solve the system of Eqs (8.1)–(8.11) numerically (Text C.2.2 in S1 File).

Mathematical model to interpret re-oxygenation experiments.

Key details of the mathematical model to interpret re-oxygenation experiments are included in Text C.3 in S1 File.

Parameter estimation and identifiability analysis

Parameter estimation and identifiability analysis are performed using profile likelihood analysis [15, 16, 52, 53], Bayesian inference [54–56], and global optimisation techniques, as now detailed. Throughout we exclude outliers in the experimental data. To detect outliers for each experiment we analyse each measurement type, Ro(t), Rn(t), and Ri(t) at each time point independently using MATLABs isoutlier function with method quartiles.

Greenspan’s model.

Parameter estimation and identifiability analysis for Greenspan’s model are first performed using profile likelihood analysis (Fig 3), see [16]. The Bayesian inference approach to estimate parameters of Greenspan’s model is discussed in the following.

Deoxygenation experiments.

Parameter estimation for the deoxygenation model (Eqs (8.1)–(8.11)) is performed using global optimisation and Bayesian inference techniques. We now explain how we estimate the fifteen model parameters, Θd = (Ro(0), αn, αh, τα, , , , sn, sh, τs, λn, λh, τλ, , . Informed by experimental measurements we set Rn(ts) = 0 and Ri(ts) = 0, but they could be included as additional parameters in future work.

We revisit Greenspan’s mathematical model for spheroids grown in normoxia and hypoxia for first estimates of sn, sh, λn and λh. Starting with experimental data from spheroids grown in normoxia, we fit a normal distribution to the initial outer radius measurements using the MATLAB fitdist function. Then for each spheroid we estimate Rc using Eq (6.2) given measurements of Ro(t) and Rn(t) and fit a normal distribution using the MATLAB fitdist function. Similarly, we estimate using Eq (6.3) and measurements of Ro(t), Rn(t) and Ri(t) and fit a normal distribution using the MATLAB fitdist function. Next, we seek to estimate the five parameters of Greenspan’s model, , using the MATLAB package MCMCstat developed by Marko Laine [57, 58]. Detailed information on the MCMCstat package is available on the GitHub repository (https://mjlaine.github.io/mcmcstat/).

Before we can use the MCMCstat package we require good first estimates of Θn and the mean squared error, mse. To provide a good first estimate for Θn, we perform global optimisation using the MATLAB GlobalSearch function with settings: fmincon sqp algorithm; MaxTime = 60 [minutes]; NumTrialPoints = 5000; and lowerbounds and upperbounds informed by fitted normal distributions for Ro(0), Rc, and and previous results [16]. To estimate the mse we use experimental measurements as observations and simulate the deoxygenation model with the estimate of Θn from global optimisation to obtain predicted values. Next we use the MCMCstat package to generate four MCMC chains with 250,000 samples and enable automatic sampling and estimation of the error standard deviation. Next we discard the first 50,000 samples of each of the four chains as burn-in, resulting in each chain containing 200,000 samples. For other MCMCstat package options we use the default settings. To assess convergence of the MCMC chains we compute the potential scale reduction factor, , [56] for each parameter and observe convergence, where convergence corresponds to (Table B in S1 File). Performing posterior checks, using 50,000 samples from the chains to generate prediction intervals and comparing with the experimental data, suggests the parameter estimates are reasonable. This process is repeated for the hypoxia data to estimate Θh (Table C in S1 File). The posteriors generated here are for sn, sh, λn and λh.

Next, we analyse the deoxygenation experimental data at each spheroid and time point independently. To estimate , , from Eq (8.8) we use global minimisation and the MCMCstat package. For estimates of for each spheroid and at each time point independently we use Eq (6.2) and measurements of Ro(t) and Rn(t). Similarly, to estimate αn, αh, τα from Eq (8.5) we use global minimisation and the MCMCstat package. For estimates of α we use Eqs (8.11) and (6.3) and measurements of Ro(t), Rn(t) and Ri(t). Note that to estimate Rc we assume that . To estimate Ro(0) we fit a normal distribution using the MATLAB fitdist function to measurements of Ro(t) at the first time point.

Next we perform a global minimisation to estimate Θd, using the estimates of Ro(0), αn, αh, τα, , , , sn, sh, λn, λh to inform lowerbounds and upperbounds. To estimate Θd we use global minimisation and the MCMCstat package (Table D In S1 File). Performing posterior checks, using 50,000 samples from the chain to generate prediction intervals and comparing with the experimental data, suggests the parameter estimates are reasonable (Figs Z, BB, and CC in S1 File).

Re-oxygenation experiments.

Parameter estimation for the re-oxygenation model formed by Eqs (S.36.1)-(S.36.11) in Text C.3 of S1 File, with parameters Θr = (Ro(0), Rn(0), αn, αh, τα, , , , sn, sh, τs, λn, λh, τλ, , , ν), is analogous to the approach used for the deoxygenation model. Summary statistics for the MCMC chains and MCMC diagnostics are shown in Tables E and F in S1 File. Comparing prediction intervals and with the experimental data suggests the parameter estimates are reasonable (Figs BB and CC in S1 File).

Experimental methods

Cell culture. The human melanoma cell lines established from primary (WM793b) and metastatic cancer sites (WM983b, WM164) were provided by Prof. Meenhard Herlyn, The Wistar Institute, Philadelphia, PA, [59]. All cell lines were previously transduced with fluorescent ubiquitination-based cell cycle indicator (FUCCI) constructs [13, 14]. Cell lines were previously genotypically characterised [13, 60–62], and authenticated by short tandem repeat fingerprinting (QIMR Berghofer Medical Research Institute, Herston, Australia). The cells were cultured in melanoma cell medium (“Tu4% medium”): 80% MCDB-153 medium (Sigma-Aldrich, M7403), 20% L-15 medium (Sigma-Aldrich, L1518), 4% fetal bovine serum (ThermoFisher Scientific, 25080–094), 5 mg ml−1 insulin (Sigma-Aldrich, I0516), 1.68 mM CaCl2 (Sigma-Aldrich, 5670) [14]. Cell lines were checked routinely for mycoplasma and tested negative using the MycoAlert MycoPlasma Detection Kit (Lonza) and polymerase chain reaction [63].

Spheroid generation, culture, and experiments. Spheroids were generated in 96-well cell culture flat-bottomed plates (3599, Corning), with 5000 total cells/well, using 50 μl total/well non-adherent 1.5% agarose to promote formation of a single spheroid per well [30]. From previous work we expect that different seeding densities, in the range 1250–10000 total cells/well, will provide similar results [15, 16]. For all experiments spheroids formed after 2 days for WM793b, WM164 and WM983b. On day 3 and 7 of each experiment 50% of the medium in each well was replaced with fresh medium (200 μl total/well). Each experiment was performed for 8 days. This choice of experimental duration is informed by previous experiments with these cell lines so that necrotic and inhibited region form prior to the end of the experiments and so that parameters of Greenspan’s model are practically identifiable and can be estimated [16].

For normoxia experiments, cells and spheroids were grown and formed in an incubator with standard settings: 37°C, 5% CO2 [13, 14]; referred to as the normoxia incubator. For hypoxia experiments, cells were cultured in the normoxia incubator and then spheroids were grown in a hypoxia incubator with settings: 37°C, 5% CO2, 2% O2. For deoxygenation experiments, cells and spheroids were formed and grown in the normoxia incubator. At the time of deoxygenation, the relevant plate(s) of spheroids were manually transferred from the normoxia incubator to the hypoxia incubator. Similarly, for the re-oxygenation experiments spheroids were grown in the hypoxia incubator then at the time of re-oxygenation the relevant plate(s) of spheroids were manually transferred to the normoxia incubator. The time to move plates was 1–2 minutes.

To estimate the outer, necrotic, inhibited, and hypoxic radii, we use a high-throughput method of mounting, clearing and imaging [64]. Results using this optical sectioning method are consistent with different methods such as cryosectioning [64]. Spheroids maintained in the relevant incubator were harvested, fixed with 4% paraformaldehyde (PFA), and stored in phosphate buffered saline solution (PBS), sodium azide (0.02%), Tween-20 (0.1%), and DAPI (1:2500) at 4°C, on days 2, 3, 4, 6 and 8 after seeding. For all hypoxia measurements, spheroids were stained with 100 mM pimonidazole for three hours, prior to fixation. Spheroids were then permeabilized with 0.5% Triton X-100 in PBS for one hour; blocked in antibody dilution buffer (Abdil) [65] for 24 hours; stained with a 1:50 anti-pimonidazole mouse IgG1 monoclonal antibody (Hypoxyprobe-1 MAb1) in Abdil for 48 hours; washed in PBS with 0.1% Tween-20 for 6 hours; placed in a 1:100 solution of goat anti-mouse Alexa Flour 647 in Abdil for 48 hours; and, finally washed for 6 hours in PBS. Then for imaging, fixed spheroids were set in place using low melting 2% agarose and optically cleared in 500 μl total/well high refractive index mounting solution (Quadrol 9% wt/wt, Urea 22% wt/wt, Sucrose 44% wt/wt, Triton X-100 0.1% wt/wt, water) for 2 days in a 24-well glass bottom plate (Cellvis, P24–1.5H-N) before imaging to ensure accurate measurements [64, 66, 67]. Images were then captured using an Olympus FV3000 confocal microscope with the 10 × objective focused on the equatorial plane of each spheroid.

As the unexpected necrotic core movement for WM983b cell line was observed in the re-oxygenation experiments, the re-oxygenation experiments were repeated for all cell lines alongside a control normoxia condition. Spheroids were cultured into three 96-well plates (3599, Corning): plate (i) control for normoxia; plate (ii) re-oxygenation 2.5 days after seeding; and, plate (iii) re-oxygenation 5.5 days after seeding. Each plate consisted of 32 spheroids of each cell line. The plates were placed inside the IncuCyte S3 live cell imaging system (Sartorius, Goettingen, Germany) incubator (37°C, 5% CO2). IncuCyte S3 settings were chosen to image with the 10 × objective. For plate (i) images were captured every 2 hours for the first three days and then every 4 hours for the remainder of the experiment. For plate (ii) and (iii) images were captured every hour for three and seven days, respectively.

Image processing. Confocal microscopy images were converted to TIFF files in ImageJ and then processed with custom MATLAB scripts that use standard MATLAB image processing toolbox functions. Area was converted to an equivalent radius (r2 = A/π). These scripts are freely available on Zenodo with DOI:10.5281/zenodo.5121093 [68], with modifications to account for pimonidazole staining and blurred central regions due to hypoxia and deoxygenation discussed in Text B in S1 File. Images captured with the IncuCyte S3 were processed with custom MATLAB scripts that use standard MATLAB image processing toolbox functions (Text D.3.1 in S1 File). Note that we identify the boundary between the inhibited and proliferating region to be the radial position at which 20% of cells are cycling.

Statistics and Reproducibility. Details of practical parameter identifiability analysis and the Bayesian inference are presented in Methods: Parameter estimation and identifiability analysis. Each radial measurements is represented as an individual data point in relevant figures, with non-filled circles representing outliers (Methods: Parameter estimation and identifiability analysis), and are summarised using box charts. Table A in S1 File details the number of measurements at each time point for each cell line and experimental data analysed during the study are available on a GitHub repository (https://github.com/ryanmurphy42/Murphy2022SpheroidOxygenAdaptation). We note that 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 late times. Data for these spheroids was excluded. In a previous study we assess experimental designs [16] and use this to inform that our sample size is sufficient in this study. Randomisation and blinding was not possible.

Supporting information

S1 File. Additional data and results.

This file includes: a summary of experimental data; experimental images; additional details of image processing; additional details of mathematical modelling; additional results for WM983b spheroids; additional results for WM793b cell line; additional results for WM164 cell line; and summary statistics of MCMC chains and MCMC diagnostics.

Fig A: Experimental images of WM983b tumour spheroids in Experiment 1 - Normoxia. Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig B: Experimental images of WM983b tumour spheroids in Experiment 2 - hypoxia. Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig C: Experimental images of WM983b tumour spheroids in Experiment 3 - deoxygenation at ts = 2 [days] (blue dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig D: Experimental images of WM983b tumour spheroids in Experiment 6 - Re-oxygenation at ts = 2.5 [days]. Scale bars are 400μm.

Fig E: Experimental images of WM983b tumour spheroids in Experiment 7 - Re-oxygenation at ts = 5.5 [days]. Scale bars are 400μm.

Fig F: Experimental images of WM793b tumour spheroids in Experiment 1 - Normoxia. Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig G: Experimental images of WM793b tumour spheroids in Experiment 2 - hypoxia. Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig H: Experimental images of WM793b tumour spheroids in Experiment 3 - deoxygenation at ts = 2 [days] (blue dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig I: Experimental images of WM793b tumour spheroids in Experiment 4 - Re-oxygenation at ts = 2 [days] (green dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig J: Experimental images of WM793b tumour spheroids in Experiment 5 - Re-oxygenation at ts = 4 [days] (green dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig K: Experimental images of WM793b tumour spheroids in Experiment 6 - Re-oxygenation at ts = 2.5 [days]. Scale bars are 400μm.

Fig L: Experimental images of WM793b tumour spheroids in Experiment 7 - Re-oxygenation at ts = 5.5 [days]. Scale bars are 400μm.

Fig M: Experimental images of WM164 tumour spheroids in Experiment 1 - Normoxia. Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig N: Experimental images of WM164 tumour spheroids in Experiment 2 - hypoxia. Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig O: Experimental images of WM164 tumour spheroids in Experiment 3 - deoxygenation at ts = 2 [days] (blue dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig P: Experimental images of WM164 tumour spheroids in Experiment 4 - Re-oxygenation at ts = 2 [days] (green dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig Q: Experimental images of WM164 tumour spheroids in Experiment 5 - Re-oxygenation at ts = 4 [days] (green dashed line). Top set of images shows spheroids with FUCCI signal only. Bottom set of images show spheroids with FUCCI signal and pimonidazole staining. Scale bars are 400μm.

Fig R: Experimental images of WM164 tumour spheroids in Experiment 6 - Re-oxygenation at ts = 2.5 [days]. Scale bars are 400μm.

Fig S: Experimental images of WM164 tumour spheroids in Experiment 7 - Re-oxygenation at ts = 5.5 [days]. Scale bars are 400μm.

Fig T: Image processing to estimate Rn(t) when the FUCCI signal in the central region of the spheroid is blurred. Images show one WM983b spheroid one day after deoxygenation. Scale bars are 200 μm. (a) Spheroid with FUCCI signals (green and magenta). (b) Spheroid with FUCCI signals and pimonidazole staining (cyan). (c-e) Blurred region identified using the polygon section tool in ImageJ. Spheroid with: (c) FUCCI-green signal only (green), (d) FUCCI-red signal only (magenta), and (e) pimonidazole staining only (cyan). (f-g) Images in (a) and (b) with blurred region removed.

Fig U: Additional experimental measurements comparing spheroids grown in normoxia and hypoxia. (a-b) Estimates of ξn(t) = Rn(t)/Ro(t). (a-b) Estimates of ξn(t) = Rn(t)/Ro(t) with time for (a) normoxia and (b) hypoxia. (c-d) Estimates of ξi(t) = Ri(t)/Ro(t) with time for (c) normoxia and (d) hypoxia. (e-f) Estimates of ξp(t) = Rp(t)/Ro(t) with time for (e) hypoxia and (f) hypoxia. (g-h) Estimates of (g) ξn(t), (h) ξi(t), and (i) ξp(t) with against outer radius, Ro(t).

Fig V: Oxygen diffusion and consumption is sufficient to describe formation and growth of necrotic core but not sufficient to describe formation and growth of inhibited region. Estimates of Rc (a) per spheroid; (b) with time; (c) against Ro(t). Estimates of α (d) per spheroid; (e) with time; (f) against Ro(t). Estimates of pi (g) per spheroid; (h) with time; (i) against Ro(t). In (a-i) each data point represents a single spheroid. Data points are only included for pi if the spheroid is in phase (ii) or phase (iii), consistent with when equations used to estimate pi are valid.

Fig W: Analysing whether production and diffusion of waste describes formation and growth of inhibited region. Estimates of the outer radius when the inhibited region first forms, (a) box chart; (b) per spheroid; (c) with time; (d) against Ro(t). In (b-d) each data point represents a single spheroid. Data points are only included only if the spheroid is in phase (ii) or phase (iii), consistent with when the equations used to estimate are valid.

Fig X: Bayesian inference to estimate parameters of Greenspan’s mathematical model for spheroids grown in normoxia and hypoxia. (a-e) Posterior densities for Greenspan model parameters: (a) Ro(0), (b) Rc, (c) , (d) s, (e) λ. Prediction intervals for (f) normoxia. (g) hypoxia. In (f-g) colour bands, in decreasing darkness, represent 50%, 75%, 95%, 97.5%, and 99.5% prediction intervals.

Fig Y: Additional results from deoxygenation experiments for WM983b spheroids. Measurements of (a) ξn(t) = Rn(t)/Ro(t), (b) ξi(t) = Ri(t)/Ro(t), and (c) ξp(t) = Rp(t)/Ro(t).

Fig Z: Prediction intervals for (a) deoxygenation experimental data and (b) deoxygenation experimental data with additional synthetic data points at t = 2.5 [days]. In (a-b) colour bands, in decreasing darkness, represent 50%, 75%, 95%, 97.5%, and 99.5% prediction intervals. Additional synthetic data at t = 2.5 not shown.

Fig AA: Direction of movement of necrotic core to edge of spheroid appears random in WM983b re-oxygenation experiments. (a) Schematic for re-oxygenation experiment, with ts = 2.5 [days]. (b) Exemplar experimental brightfield images of a single spheroid with image processing to detect spheroid boundary and centroid (red) and necrotic core boundary and centroid (cyan). (c) Trajectory of centroid of spheroid for spheroid imaged in (b). (d) Trajectory of centroid of necrotic core for spheroid imaged in (b). (e) Five exemplar trajectories of the centroid of the necrotic core relative to the centroid of the spheroid. In (e) the thick blue trajectory corresponds to spheroid imaged in (b).

Fig BB: Additional results for WM793b spheroids. (a-d) Mechanisms governing tumour spheroid growth in normoxia and hypoxia. (a) Box chart for estimated outer radius when necrotic region forms, Rc [μm]. (b) Box chart for estimated oxygen consumption rate, α [m3 kg−1 s−1]. (c) Comparison of measured and predicted Rp(t) when pimonidazole staining is present. Note this does not include images where pimonidazole staining is present but does not surround the necrotic core, for example Day 8 of Fig G. (d) Box chart for estimated oxygen partial pressure defining inhibited region from hypothesis 1, pi [%]. (e-i) Experimental data and prediction intervals for (e) normoxia, (f) hypoxia, (g) deoxygenation, (h) re-oxygenation with ts = 2 [days], and (i) re-oxygenation with ts = 4 [days].

Fig CC: Additional results for WM164 spheroids. (a-d) Mechanisms governing tumour spheroid growth in normoxia and hypoxia. (a) Box chart for estimated outer radius when necrotic region forms, Rc [μm]. (b) Box chart for estimated oxygen consumption rate, α [m3 kg−1 s−1]. (c) Comparison of measured and predicted Rp(t) when pimonidazole staining is present. (d) Box chart for estimated oxygen partial pressure defining inhibited region from hypothesis 1, pi [%]. (e-i) Experimental data and prediction intervals for (e) normoxia, (f) hypoxia, (g) deoxygenation, (h) re-oxygenation with ts = 2 [days], and (i) re-oxygenation with ts = 4 [days].

Table A: Number of spheroids imaged with confocal microscopy for the WM983b, WM793b, and WM164 cell lines. For WM983b re-oxygenation experiments, denoted by *, we focus on brightfield images.

Table B: Greenspan’s model parameters and MCMC diagnostics for normoxia experiments. Summary statistics of the MCMC chains include: mean; standard deviation, σ; and, 25%, 50%, and 75% quartiles, Q25%, Q50%, and Q75%, respectively. To assess convergence of the MCMC chains we compute the potential scale reduction factor, , [56] where convergence corresponds to .

Table C: Greenspan’s model parameters and MCMC diagnostics for hypoxia experiments. Summary statistics of the MCMC chains include: mean; standard deviation, σ; and, 25%, 50%, and 75% quartiles, Q25%, Q50%, and Q75%, respectively. To assess convergence of the MCMC chains we compute the potential scale reduction factor, , [56] where convergence corresponds to .

Table D: Deoxygenation model parameters and MCMC diagnostics for deoxygenation experiments. Summary statistics of the MCMC chains include: mean; standard deviation, σ; and, 25%, 50%, and 75% quartiles, Q25%, Q50%, and Q75%, respectively. To assess convergence of the MCMC chains we compute the potential scale reduction factor, , [56] where convergence corresponds to .

Table E: Re-oxygenation model parameters and MCMC diagnostics for re-oxygenation experiments with ts = 2 [days]. Summary statistics of the MCMC chains include: mean; standard deviation, σ; and, 25%, 50%, and 75% quartiles, Q25%, Q50%, and Q75%, respectively. To assess convergence of the MCMC chains we compute the potential scale reduction factor, , [56] where convergence corresponds to .

Table F: Re-oxygenation model parameters and MCMC diagnostics for re-oxygenation experiments with ts = 4 [days]. Summary statistics of the MCMC chains include: mean; standard deviation, σ; and, 25%, 50%, and 75% quartiles, Q25%, Q50%, and Q75%, respectively. To assess convergence of the MCMC chains we compute the potential scale reduction factor, , [56] where convergence corresponds to .

https://doi.org/10.1371/journal.pcbi.1010833.s001

(PDF)

S1 Movie. Necrotic core movement and removal in WM983b spheroids in response to re-oxygenation.

Time-lapse brightfield microscopy movie for three days following re-oxygenation at ts = 2.5 [days]. The necrotic core of the spheroid is initially visible as a dark central region. As time progresses the necrotic core is located closer to the edge of the spheroid and the symmetric internal structure is lost. At later times the necrotic core appears to exit the spheroid as a single object.

https://doi.org/10.1371/journal.pcbi.1010833.s002

(MP4)

S2 Movie. Necrotic core movement in WM983b spheroids in response to re-oxygenation.

Time-lapse brightfield microscopy movie for seven days following re-oxygenation at ts = 5.5 [days]. The necrotic core of the spheroid is initially visible as a dark central region. At later times the necrotic core is close to the edge of the spheroid but does not exit as a single object. As the spheroid grows necrotic matter forms at the centre of the spheroid and appears to merge with the necrotic matter located closer to the periphery.

https://doi.org/10.1371/journal.pcbi.1010833.s003

(MP4)

S3 Movie. Loss of the necrotic core in WM164 re-oxygenation experiments.

Time-lapse brightfield microscopy movie for three days following re-oxygenation at ts = 2.5 [days]. The necrotic core of the spheroid is initially visible as a dark central region. As time progresses mass from the necrotic core at the centre of the spheroid moves towards the periphery and exits the spheroid. The spheroid then appears to resume growth.

https://doi.org/10.1371/journal.pcbi.1010833.s004

(MP4)

Acknowledgments

We thank Dr Alexander P. Browning and Dr Patrick B. Thomas for helpful discussions, and John Blake for guidance using IncuCyte. 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, The Wistar Institute, Philadelphia, PA, for providing the cell lines.

Citation: Murphy RJ, Gunasingh G, Haass NK, Simpson MJ (2023) Growth and adaptation mechanisms of tumour spheroids with time-dependent oxygen availability. PLoS Comput Biol 19(1): e1010833. https://doi.org/10.1371/journal.pcbi.1010833

References

1. Bader SB, Dewhirst MW, Hammond EM. Cyclic hypoxia: An update on its characteristics, methods to measure it and biological implications in cancer. Cancers. 2021;13:23.

2. Kumar B, Adebayo AK, Prasad M, Capitano ML, Wang R, Bhat-Nakshatri P, et al. Tumor collection/processing under physioxia uncovers highly relevant signaling networks and drug sensitivity. Science Advances. 2022;8:eabh3375. pmid:35020422

3. Vaupel P, Kallinowski F, Okunieff P. Blood flow, oxygen and nutrient supply and metabolic environment of human tumours: A review. Cancer Research. 1989;49:6449–6465.

4. Mueller-Klieser W, Freyer J, Sutherland R. Influence of glucose and oxygen supply conditions on the oxygenation of multicellular spheroids. British Journal of Cancer. 1986;53:345–353. pmid:3964538

5. Hirschhaeuser F, Menne H, Dittfled C, West J, Mueller-Klieser W, Kunz-Schughart LA. Multicellular tumor spheroids: An underestimated tool is catching up again. Journal of Biotechnology. 2010;148:3–15. pmid:20097238

6. Costa EC, Moreira AF, de Melo-Diogo D, Gaspar VM, Carvalho MP, Correia IJ. 3D tumor spheroids: an overview on the tools and techniques used for their analysis. Biotechnology Advances. 2016;34:1427–1441. pmid:27845258

7. Nath S, Devi GR. Three-dimensional culture systems in cancer research: Focus on tumor spheroid model. Pharmacology & Therapeutics. 2016;163:94–108. pmid:27063403

8. Mehta G, Hsiao AY, Ingram M, Luker GD, Takayama S. Opportunities and challenges for use of tumor spheroids as models to test drug delivery and efficacy. Journal of Controlled Release. 2012;10:192–204. pmid:22613880

9. Weiswald L, Bellet D, Dangles-Marie V. Spherical cancer models in tumor biology. Neoplasia. 2015;17:1–15. pmid:25622895

10. Greenspan HP. Models for the growth of a solid tumor by diffusion. Studies in Applied Mathematics. 1972;51:317–340.

11. Le Maout V, Alessandri K, Gurchenkov B, Bertin H, Nassoy P, Sciumè G. Role of mechanical cues and hypoxia on the growth of tumor cells in strong and weak confinement a dual in vitro in silico approach. Science Advances. 2020;6:eaaz7130. pmid:32232163

12. Beaumont KA, Mohana-Kumaran N, Haass NK. Modeling melanoma in vitro and in vivo. Healthcare. 2014;2:27–46.

13. Haass NK, Beaumont KA, Hill DS, Anfosso A, Mrass P, Munoz MA, et al. Real-time cell cycle imaging during melanoma growth, invasion, and drug response. Pigment Cell & Melanoma Research. 2014;27:764–776. pmid:24902993

14. Spoerri L, Beaumont KA, Anfosso A, Haass NK. Real-time cell cycle imaging in a 3D cell culture model of melanoma. Methods in Molecular Biology. 2017;1612:401–416. pmid:28634959

15. Browning AP, Sharp JA, Murphy RJ, Gunasingh G, Lawson B, Burrage K, et al. Quantitative analysis of tumour spheroid structure. eLife. 2021;10:e73020. pmid:34842141

16. Murphy RJ, Browning AP, Gunasingh G, Haass NK, Simpson MJ. Designing and interpreting 4D tumour spheroid experiments. Communications Biology. 2022;5:91. pmid:35075254

17. Gomes A, Guillaume L, Grimes DR, Fehrenbach J, Lobjois V, Ducommun B. Oxygen partial pressure is a rate-limiting parameter for cell proliferation in 3D spheroids grown in physioxic culture condition. PLoS One. 2016;11:e0161239. pmid:27575790

18. Grimes DR, Kannan P, McIntyre A, Kavanagh A, Siddiky A, Wigfield S, et al. The role of oxygen in avascular tumor growth. PLoS One. 2016;11:e0153692. pmid:27088720

19. Muz B, de la Puente P, Azab F, Azab AB. The role of hypoxia in cancer progression, angiogenesis, metastasis, and resistance to therapy. Hypoxia. 2015;3:83–92. pmid:27774485

20. McKeown SR. Defining normoxia, physoxia and hypoxia in tumours-implications for treatment response. British Journal of Radiology. 2014;87:20130676. pmid:24588669

21. Carreau A, El Hafny-Rahbi B, Matejuk A, Grillon C, Kieda C. Why is the partial oxygen pressure of human tissues a crucial parameter? Small molecules and hypoxia. Journal of Cellular and Molecular Medicine. 2011;15:1239–1253. pmid:21251211

22. Celora GL, Bader SB, Hammond EM, Maini PK, Pitt-Francis J, Byrne HM. A DNA-structured mathematical model of cell-cycle progression in cyclic hypoxia. Journal of Theoretical Biology. 2022;545:111104. pmid:35337794

23. Lee P, Chandel NS, Simon MC. Cellular adaptation to hypoxia through hypoxia inducible factors and beyond. Nature Reviews Molecular Cell Biology. 2020;21:268–283. pmid:32144406

24. Goto T, Kaida A, Miura M. Visualizing cell-cycle kinetics after hypoxia/reoxygenation in HeLa cells expressing fluorescent ubiquitination-based cell cycle indicator (Fucci). Experimental Cell Research. 2015;339:389–396. pmid:26500111

25. Grist SM, Nasseri SS, Laplatine L, Schmok JC, Yao D, Hua J, et al. Long-term monitoring in a microfluidic system to study tumour spheroid response to chronic and cycling hypoxia. Scientific Reports. 2019;9:17782. pmid:31780697

26. Fridman IB, Ugolini GS, VanDelinder V, Cohen S, Konry T. High throughput microfluidic system with multiple oxygen levels for the study of hypoxia in tumor spheroids. Biofabrication. 2021;13:035037.

27. Riffle S, Hegde RS. Modeling tumor cell adaptations to hypoxia in multicellular tumor spheroids. Journal of Experimental & Clinical Cancer Research. 2017;36:102. pmid:28774341

28. Al-Ani A, Toms D, Kondro D, Thundathil J, Yu Y, Ungrin M. Oxygenation in cell culture: Critical parameters for reproducibility are routinely not reported. PLoS One. 2018;13:e0204269. pmid:30325922

29. Sakaue-Sawano A, Kurokawa H, Morimura T, Hanyu A, Hama H, Osawa H, et al. Visualizing spatiotemporal dynamics of multicellular cell-cycle progression. Cell. 2008;132:487–498. pmid:18267078

30. Spoerri L, Gunasingh G, Haass NK. Fluorescence-based quantitative and spatial analysis of tumour spheroids: a proposed tool to predict patient-specific therapy response. Frontiers in Digital Health. 2021;3:668390. pmid:34713141

31. Araujo RP, McElwain DLS. A history of the study of solid tumour growth: the contribution of mathematical modelling. Bulletin of Mathematical Biology. 2004;66:1039. pmid:15294418

32. Byrne HM. Dissecting cancer through mathematics: from the cell to the animal model. Nature Reviews Cancer. 2010;10:221–230. pmid:20179714

33. Roose T, Chapman SJ, Maini PK. Mathematical models of avascular tumor growth. SIAM Review. 2007;49:179–208.

34. Bull JA, Byrne HM. The hallmarks of mathematical oncology. Proceedings of the IEEE. 2022;110:523–540.

35. Bull JA, Mech F, Quasier T, Waters SL, Byrne HM. Mathematical modelling reveals cellular dynamics within tumour spheroids. PLoS Computational Biology. 2020;16:e1007961. pmid:32810174

36. Jagiella N, Müller B, Müller M, Vignon-Clementel IE, Drasdo D. Inferring growth control mechanisms in growing multi-cellular spheroids of NSCLC cells from spatial-temporal image data. PLoS Computational Biology. 2016;12:e1004412. pmid:26866479

37. Grimes DR, Kelly C, Bloch K, Partridge M. A method for estimating the oxygen consumption rate in multicellular tumour spheroids. Journal of the Royal Society Interface. 2014;11:20131124. pmid:24430128

38. Browning AP, Simpson MJ. Geometric analysis enables biological insight from complex non-identifiable models using simple surrogates. arXiv. 2022.

39. Jin W, Spoerri L, Haass NK, Simpson MJ. Mathematical model of tumour spheroid experiments with real-time cell cycle imaging. Bulletin of Mathematical Biology. 2021;83:1–23. pmid:33743088

40. Klowss JJ, Browning AP, Murphy RJ, Carr EJ, Plank MJ, Gunasingh G, et al. A stochastic mathematical model of 4D tumour spheroids with real-time fluorescent cell cycle labelling. Journal of the Royal Society Interface. 2022;19:20210903. pmid:35382573

41. Lewin TD, Maini PK, Moros EG, Enderling H, Byrne HM. The evolution of tumour composition during fractionated radiotherapy implications for outcome. Bulletin of Mathematical Biology. 2018;80:1207–1235. pmid:29488054

42. Santilli G, Lamorte G, Carlessi L, Ferrari D, Nodari LR, Binda E, et al. Mild hypoxia enhances proliferation and multipotency of human neural stem cells. PLoS One. 2010;5:e8575. pmid:20052410

43. Grayson WL, Zhao F, Bunnell B, Ma T. Hypoxia enhances proliferation and tissue formation of human mesenchymal stem cells. Biochemical and Biophysical Research Communications.2007;358:948–953. pmid:17521616

44. Danet GH, Pan Y, Luongo JL, Bonnet DA, Simon MC. Expansion of human SCID-repopulating cells under hypoxic conditions. The Journal of Clinical Investigation. 2003;112:126–135. pmid:12840067

45. Ward JP, King JR. Mathematical modelling of avascular-tumour growth. II: Modelling growth saturation. IMA Journal of Mathematics Applied in Medicine and Biology. 1999;16:171–211. pmid:10399312

46. Lewin TD, Maini PK, Moros EG, Enderling H, Byrne HM. A three phase model to investigate the effects of dead material on the growth of avascular tumours. Mathematical Modelling of Natural Phenomena. 2020;15:22.

47. Woolley TE, Gaffney EA, Oliver JM, Baker RE, Waters SL, Goriely A. Cellular blebs: pressure-driven, axisymmetric, membrane protrusions. Biomechanics and Modeling in Mechanobiology. 2014;13:463–476. pmid:23857038

48. Landman KA. On the crenation of a compound liquid droplet. Studies in Applied Mathematics. 1983;69:51–73.

49. Leedale J, Colley HE, Gaskell H, Williams DP, Bearon RN, Chadwick AE, et al. In silico-guided optimisation of oxygen gradients in hepatic spheroids. Computational Toxicology. 2019;12:100093.

50. Cavadas MAS, K NL, Cheong A. Hypoxia-inducible factor (HIF) network insights from mathematical models. Cell Communication and Signaling. 2013;11:42. pmid:23758895

51. Leedale J, Herrmann A, Bagnall J, Fercher A, Papkovsky D, Sée V, et al. Modeling the dynamics of hypoxia inducible factor-1α (HIF-1α) within single cells and 3D cell culture systems. Mathematical Biosciences. 2014;258:33–43. pmid:25245610

52. Pawitan Y. In All Likelihood: Statistical Modelling And Inference Using Likelihood. Oxford, UK: Oxford University Press; 2001.

53. Simpson MJ, Baker RE, Vittadello ST, Maclaren OJ. Practical parameter identifiability for spatio-temporal models of cell invasion. Journal of the Royal Society Interface. 2020;17:2020055. pmid:32126193

54. Browning AP, Haridas P, Simpson MJ. A Bayesian sequential learning framework to parameterise continuum models of melanoma invasion into human skin. Bulletin of Mathematical Biology. 2019;81:676–698. pmid:30443704

55. Collis J, Connor AJ, Paczkowski M, Kannan P, Pitt-Francis J, Byrne HM, et al. Bayesian calibration, validation and uncertainty quantification for predictive modelling of tumour growth: a tutorial. Bulletin of Mathematical Biology. 2017;79:939–974. pmid:28290010

56. Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. Bayesian Data Analysis. 3rd ed. New York: Chapman and Hall/CRC; 2013.

57. Haario H, Saksman E, Tamminen J. An adaptive Metropolis algorithm. Bernoulli. 2001;7:223–242.

58. Haario H, Laine M, Mira A, Saksman E. DRAM: Efficient adaptive MCMC. Statistics and Computing. 2006;16:339–354.

59. Hsu MY, Elder DE, Herylyn M. Melanoma: the Wistar melanoma (WM) cell lines. Human Cell Culture. 2002;1:259–274.

60. Hoek KS, Schlegel NC, Brafford P, Sucker A, Ugurel S, Kumar R, et al. Metastatic potential of melanomas defined by specific gene expression profiles with no BRAF signature. Pigment Cell Research. 2006;19:290–302. pmid:16827748

61. Smalley KSM, Contractor R, Haass NK, Kulp AN, Atilla-Gokcumen GE, Williams DS, et al. An organometallic protein kinase inhibitor pharmacologically activates p53 and induces apoptosis in human melanoma cells. Cancer Research. 2007;67:209–217. pmid:17210701

62. Smalley KSM, Contractor R, Haass NK, Lee JT, Nathanson KL, Medina CA, et al. Ki67 expression levels are a better marker of reduced melanoma growth following MEK inhibitor treatment than phospho-ERK levels. British Journal of Cancer. 2007;96:445–449. pmid:17245336

63. Uphoff CC, Drexler HG. Detecting mycoplasma contamination in cell cultures by polymerase chain reaction. Methods in Molecular Biology. 2011;731:93–103. pmid:21516400

64. Gunasingh G, Browning AP, Haass NK. Rapid optical clearing for semi-high-throughput analysis of tumour spheroids. Journal of Visualized Experiments. 2022;186:e64103

65. Cold Spring Harbor Laboratory Press. Antibody Dilution Buffer (Abdil); 2018.

66. Costa EC, Silva DB, Moreira AF, Correia IJ. Optical clearing methods: an overview of the techniques used for the imaging of 3D spheroids. Biotechnology & Bioengineering. 2019;116:2742–2763. pmid:31282993

67. Susaki EA, Shimizu C, Kuno A, Tainaka K, Li X, Nishi K, et al. Versatile whole-organ/body staining and imaging based on electrolyte-gel properties of biological tissues. Nature Communications. 2020;11:1982. pmid:32341345

68. Browning AP, Murphy RJ. Image processing algorithm to identify structure of tumour spheroids with cell cycle labelling. Zenodo. 2021;