Minimal Morphoelastic Models of Solid Tumour Spheroids: A Tutorial

Tumour spheroids have been the focus of a variety of mathematical models, ranging from Greenspan’s classical study of the 1970 s through to contemporary agent-based models. Of the many factors that regulate spheroid growth, mechanical effects are perhaps some of the least studied, both theoretically and experimentally, though experimental enquiry has established their significance to tumour growth dynamics. In this tutorial, we formulate a hierarchy of mathematical models of increasing complexity to explore the role of mechanics in spheroid growth, all the while seeking to retain desirable simplicity and analytical tractability. Beginning with the theory of morphoelasticity, which combines solid mechanics and growth, we successively refine our assumptions to develop a somewhat minimal model of mechanically regulated spheroid growth that is free from many unphysical and undesirable behaviours. In doing so, we will see how iterating upon simple models can provide rigorous guarantees of emergent behaviour, which are often precluded by existing, more complex modelling approaches. Perhaps surprisingly, we also demonstrate that the final model considered in this tutorial agrees favourably with classical experimental results, highlighting the potential for simple models to provide mechanistic insight whilst also serving as mathematical examples.


Introduction
Cancer is a disease that impacts the lives of tens of millions of people worldwide each year and represents a leading cause of death (Sung et al. 2021). The growing prevalence and severity of the disease have driven rapid advancements in our understanding of the biology that underpins tumour growth, as highlighted by the evolving characterisation of the Hallmarks of Cancer in the renowned works of Weinberg (2000, 2011). A subset of this vast body of research has considered the broad range of stimuli that are known to affect the behaviour of biological cells and tissues, including cancer cells. These stimuli include, but are not limited to, the availability of nutrients for growth, mechanical forces acting on tissues, and electric fields (Vaupel et al. 1989;Pavlova and Thompson 2016;Northcott et al. 2018;Sengupta and Balla 2018;Kolosnjaj-Tabi et al. 2019). Whilst the study of a single stimulus is often difficult or intractable in vivo, experimental assays have provided a means to focus on one or two stimuli at a time. For instance, a common approach for studying the early stages of avascular tumour growth is to consider three-dimensional collections of cancer cells known as tumour spheroids (Hirschhaeuser et al. 2010), which are thought to better emulate in vivo environments than alternative two-dimensional assays whilst still enabling the targeted study of tumour growth stimuli. Spheroid assays have been used to study the impact on tumour growth of multiple stimuli, such as how nutrient availability affects cancer development (Kunz-Schughart et al. 1998;Murphy et al. 2022), and to gain insight into the mechanical inhibition of growth, as in the now-classical work of Helmlinger et al. (1997), which we consider in detail below.
In addition to the range of experimental investigations that have involved the use of tumour spheroids, a multitude of mathematical models have been developed to study spheroids. These efforts, which are part of the emergent field of mathematical oncology, range from simple single-compartment ordinary differential equation (ODE) models to complex, multiscale schemes and hybridised partial differential equation (PDE) and agent-based methods, which differ in complexity, spatial resolution, and scale. One of the earliest and best known models is that of Greenspan (1972), which considers how the composition of a tumour spheroid evolves as the growing tumour limits the availability of diffusing nutrients (in this case oxygen) to the central core of the spheroid. Greenspan's PDE approach, in which the behaviour of cells is driven by the local nutrient concentration, has since been adapted by many authors and adds to the breadth of mathematical methods that have been employed in the study of cancer. The reviews of Araujo and McElwain (2004), Roose et al. (2007), and Bull and Byrne (2022) provide a comprehensive summary of these theoretical approaches.
Since Greenspan's early work, many mathematical models have focussed on exploring how nutrient availability and spatial constraints limit tumour growth (Ward and King 1997;Sherratt and Chaplain 2001;Murphy et al. 2022). In contrast, however, the notion of mechanical feedback remains relatively unexplored in theoretical works, despite mechanical effects being increasingly appreciated as significant in many biological settings. For instance, the work of Helmlinger et al. (1997) demonstrated that mechanical resistance to growth can markedly limit the growth of tumour spheroids, with resistance in this case being imparted via an agarose gel that surrounds the spheroids. More recent experimental studies add weight to Helmlinger et al.'s conclusions, such as that of Cheng et al. (2009), which considered the effects of externally imposed stresses on tumours and measured the impacts of mechanical stress on cell proliferation and apoptosis. Notwithstanding these experimental results, there is no consensus about how mechanical cues alter growth dynamics on the tissue and cell scales. This uncertainty has spawned a range of phenomenological continuum models of mechanically influenced tumour growth (Chen et al. 2001;Roose et al. 2003;Byrne and Drasdo 2009;Mollica 2004, 2002;Ambrosi et al. 2017;Ambrosi and Preziosi 2009;Byrne 2003), which have successfully reproduced both tumour growth curves and profiles of accumulated solid stress (Nia et al. 2018). These theoretical studies have made different modelling choices, most notably in the posited constitutive couplings between mechanics and growth. A key challenge is finding a coupling that represents the least complex relation needed to generate experimentally observed profiles of growth and stress. Such simplicity is often desirable in mathematical models when detailed understanding of the biological mechanisms is lacking, and such 'minimal ingredients' models generally facilitate both ready interpretation and analytical study, the latter of which can provide rigorous characterisation of model dynamics and behaviours. Such characterisations are largely absent from existing solid mechanical models of spheroid growth. In particular, it remains to be established whether agreement between numerical solutions of existing mechanical tumour models and experimental data depends strongly on the particular parameter regimes employed, or whether they reproduce the observed phenomena more generally.
Motivated by these observations, the scientific aim of this tutorial is to develop a minimal model of tumour spheroid growth that reproduces observed growth dynamics, under varying external conditions, and permits rigorous characterisation of model behaviours. In pursuit of this goal, we will adopt an iterative and expository approach to model development, beginning with a simple, established foundation and successively posing extensions and modifications in order to realise a number of desirable properties. To facilitate the development of such a simple mathematical model, each of our models will be based on the solid mechanical framework of morphoelasticity, introduced by Rodriguez et al. (1994) and reviewed by Ambrosi et al. (2011) andKuhl (2014). The theory of morphoelasticity has been applied broadly to problems of biological growth and often leads to models that are analytically tractable and numerically straightforward, such as a recent model of the human eye (Kimpton et al. 2021), as described by Goriely (2017).
Throughout our exploration of tumour growth models, we will seek a model that exhibits a number of properties, each of which will focus on robustly reproducing features of experimentally observed tumour dynamics and having behaviour consistent Mechanically influenced tumour spheroid growth. Experimentally observed growth dynamics of tumour spheroids in various extracellular media are indicated by empty circles, as reported by Helmlinger et al. (1997), digitised by Yan et al. (2021). Tumours growing in agarose gels of higher concentration experience reduced growth compared to evolution in free suspension (shown black). Classical models of spheroid growth are capable of capturing growth in free suspension to excellent accuracy, as highlighted by the least-squares fit of a Greenspan-inspired model to the free-suspension data, shown in black. Details of the fitting are found in Appendix A with our understanding of tumour growth. Before we describe these properties, it will be helpful to first illustrate a known impact of mechanical factors on tumour growth dynamics. To this end, in Fig. 1 we showcase a selection of the experimental data reported by Helmlinger et al. (1997), as digitised by Yan et al. (2021); various growth curves correspond to differing levels of mechanical resistance exerted on growing spheroids embedded in agarose gels of various concentrations and, hence, stiffnesses. From these datapoints alone, it is clear that the mechanical properties of the external medium can significantly impact tumour growth dynamics in this system, with increasing stiffness reducing a spheroid's capability to grow.
The simplest models of tumour growth neglect mechanical effects and assume that growth depends only on the availability of diffusing nutrients, such as oxygen. Such models have proved successful in reproducing the growth of tumour spheroids in free suspension. We illustrate this by fitting a model inspired by Greenspan's seminal model (Greenspan 1972) of nutrient-limited growth to the free-suspension growth curve in Fig. 1; full details of the employed model are provided in Sect. 3.1, and the fitting process is summarised in Appendix A. Whilst the mathematical model exhibits excellent agreement with the experimental observations of Helmlinger et al. (1997) for tumour growth in free suspension, it is not able to simultaneously fit all three datasets shown, as mechanical effects are neglected. Motivated by this, a further aim of this tutorial is to formulate Greenspan's classical approach within a solid mechanical framework. In particular, we will seek to capture the phenomena reported by Helmlinger et al., with the specific goal of fitting a mechanical model to the full range of dynamics shown in Fig. 1.
With this additional goal in mind, we specify a basic requirement of any model that we will develop: it must give rise to growth curves that are qualitatively similar to those of Fig. 1. More specifically, growth curves must be monotonic and should capture a profile of development that is approximately of exponential, linear, then saturating character, as observed in Helmlinger's et al. experiments and as is canonical of tumour spheroid growth (Bull and Byrne 2022). Here, we define saturating dynamics to be those whose growth rate tends monotonically to zero, noting that this does not guarantee that the tumour size itself is bounded 1 . However, a stable nonzero steady state of tumour size is often associated with these growth curves, the existence of which we will also look to guarantee analytically.
As solid mechanics will necessarily play a key role, we will also seek certain properties that relate to the mechanical stress within the spheroid. A minimal such requirement is that the solid stresses in the tumour should be bounded, so that a model will not predict that a tumour experiences arbitrarily large internal stresses as time increases. Whilst this intuitive property might seem to be elementary to realise, we will show that it does not necessarily hold in even simple cases and therefore, requires appropriate consideration. Finally, we will also seek out the ability of a model to reproduce profiles of solid stress that are similar to those that have been estimated from experimental data. In particular, we will aim to qualitatively match profiles of residual stress, the term given to solid stresses accumulated in a tumour during growth that remain present when the tissue experiences no external mechanical load. These mechanical properties, along with the features introduced above, are summarised in Table 1.
In summary, in this tutorial we will seek to develop a continuum model of avascular spheroid growth, one in which growth is regulated by mechanical effects and nutrient availability. In attempting to realise the noted desirable properties, we will describe and explore a number of intermediate models, highlighting how an iterative, minimalistic approach to model construction can provide insight into emergent behaviours and facilitate the development of mathematical models that are simple, interpretable, and robust. We will begin by incorporating the established foundation of Greenspan (1972) within the framework of morphoelasticity, striving throughout for simplicity in order to enable exploratory and analytical study.

Geometry and Set-Up
Throughout this tutorial, we will model a tumour spheroid as a morphoelastic solid (Goriely 2017), considering its deformation due to the processes of growth and elastic relaxation whilst assuming strict spherical symmetry, as illustrated in Fig. 2. Deformations are captured by the time-dependent gradient tensor F, which encodes the mapping from the Lagrangian spherical coordinates (R, , ) of the initial spheroid, assumed to be stress free, to the Eulerian spherical coordinates (r , θ, φ) that represent the deformed configuration of the spheroid. With our assumption of spherical Table 1 Desirable properties of a minimal model of tumour spheroid growth. The Greenspan-inspired model used in Fig. 1 where r = r (R, t) is a function of time t ≥ 0 and the Lagrangian radial coordinate R ∈ [0, B], with r (R, 0) = R and where B is the initial radius of the spheroid. Throughout, we assume that the mapping between the initial and deformed configurations is such that det F(R, t) > 0 for all t ≥ 0 and all R ∈ [0, B], so that the deformation preserves the orientation of the material and is locally injective for all t. We also assume that the spheroid undergoes no topological changes, so that r (0, t) = 0, and we denote the outer radius of the deformed tumour by b(t) := r (B, t).

Morphoelasticity
Following Goriely (2017), the central assumption of morphoelasticity is that the tensor F can be multiplicatively decomposed into two components: one that represents the growth of the material and another that captures the elastic response of the grown material: where A and G are second-order tensors that represent the effects of elasticity and growth, respectively. Appealing to spherical symmetry, this relation can be explicitly written as ⎡ where α r , α θ , α φ and γ r , γ θ , γ φ are non-negative scalar elastic stretches and growth stretches, respectively. In particular, γ r captures material growth in the radial direction, whilst γ θ ≡ γ φ encodes circumferential growth. The equivalence between γ θ and γ φ is a consequence of spherical symmetry and, analogously, we have α θ ≡ α φ . Here, a growth stretch larger than one corresponds to the addition of material, whilst resorption occurs when the growth stretch is less than 1 (but larger than 0). Taking the determinant of Eq. (3) yields the following quasistatic partial differential equation for the radial coordinate where we have replaced α φ and γ φ by α θ and γ θ , respectively. Henceforth, for simplicity and consistent with our goal of pursuing a minimal model of spheroid growth, we will assume that the material is incompressible, so that elastic deformations do not alter the volume. With this being equivalent to imposing the condition det A = α r α 2 θ = 1, Eq. (4) now takes on the simplified form and elasticity no longer plays an explicit role in the instantaneous configuration of the spheroid. In particular, given the instantaneous growth stretches, one may readily integrate this equation, with the condition r (0, t) = 0, to determine r (R, t), without further consideration of mechanical effects. From this, one might be tempted to neglect mechanics altogether; however, we will see that incorporating mechanics in an evolution law for the growth stretches is needed in order to reproduce the properties listed in Table 1. Writing α r = α −2 θ and defining α := α θ , the decomposition of Eq. (3) then yields which will later enable us to eliminate the elastic stretch α from calculations and to transform between Lagrangian (R) and Eulerian (r ) coordinates.

Mechanics and Constitutive Assumptions
We will assume that the spheroid is composed of an isotropic hyperelastic material, such that it is characterised by a strain energy function W . With W thereby a function of the principal stretches, which are the eigenvalues of the right Cauchy-Green tensor AA T , we may generically write W = W (α r , α θ , α φ ). The Cauchy stress tensor σ is then given by where the pressure p is the Lagrange multiplier required to enforce incompressibility and the (i, j) th entry of ∂ W /∂ A is defined to be the derivative of W with respect to the ( j, i) th entry of A. In our spherically symmetric system, we write leading to the scalar relations for the radial stress σ r (R, t) and the hoop stress σ θ (R, t). Eliminating the Lagrange multiplier and definingW such that W =W (α), we obtain the single equation relating the stresses in the circumferential and radial directions. In the absence of body forces, conservation of linear momentum reads where the divergence operator is with respect to the Eulerian coordinate system. Using the assumed spherical symmetry and the relation of Eqs. (10), (11) gives where we are treating σ r as a function of r , recalling that r = r (R, t). Seeking simplicity, we suppose that the tumour is a neo-Hookean material, so that where μ > 0 is the material-dependent shear modulus and α −4 + 2α 2 is the first invariant (trace) of the right Cauchy-Green tensor in our simplified setting. Hence, the conservation equation for the radial stress becomes Transforming to Lagrangian coordinates using Eq. (6a) and eliminating α via Eq. (6b), we obtain a quasistatic Lagrangian PDE for the radial stress, To model resistance to expansion due to external material, we prescribe a compressive radial stress on the boundary of the spheroid that opposes growth. Specifically, we impose where κ ≥ 0 encodes the stiffness of the surrounding medium, though we note that generalisations of this relation are straightforward to accommodate. With this boundary condition, which reduces to growth in free suspension when κ = 0, we can integrate Eq. (15) inwards from the boundary to yield σ r , with integration being performed numerically in practice. We can then construct σ θ via or, equivalently,

Growth Dynamics
It remains to specify how the growth stretches evolve with time from their common initial value of unity. Even in the case of isotropic growth, which we will assume hereafter, it is unclear how the growth of tissue should depend on the state of the tumour. Indeed, it is this dependence that we will explore and vary throughout this tutorial, seeking a phenomenological growth law that gives rise to the canonical growth dynamics defined earlier (exponential, linear, saturating) whilst being free from unphysical behaviours. Each growth law will take the same basic form, which we state generically in terms of γ := γ r = γ θ as 1 γ for a fixed rate constant k and initial condition of unity, where the function f encodes the dependence of growth on the stress tensor σ (r , t) and a generic diffusible nutrient with concentration c(r , t). We interpret this nutrient as oxygen, as in Greenspan (1972). Seeking to model the growth of an avascular tumour, so that the only source of nutrients is via diffusion from the outer boundary of the spheroid, we suppose that the non-negative nutrient concentration c(r , t) is governed by the PDE wherever c is nonzero, where D is the diffusion coefficient of the nutrient in the tumour medium and λ is the constant consumption rate of nutrient by the tissue. We impose the simple boundary condition c(b(t), t) = c ∞ at the surface of the tumour, which allows us to define the diffusive lengthscale L := √ Dc ∞ /λ and characteristic timescale T = 1/kc ∞ of the spheroid problem. Seeking a quasistatic solution of Eq. (20) that is a function purely of the Eulerian coordinate r , noting that the timescales of diffusion are typically much shorter than those of biological growth, Eq. (20) reduces to At the centre of the tumour, symmetry considerations imply a no-flux condition ∂c/∂r (0, t) = 0, which leads to the solution However, as the nutrient concentration must be non-negative, this solution is not valid if b(t) is sufficiently large, with c predicted to be negative at the core of the spheroid. There are at least two resolutions to this problem. One route simply modifies the consumption term in the differential equation to include a dependence on the nutrient concentration itself, with the simplest specifying that the consumption is proportional to the concentration. This gives rise to a non-negative solution that can be written in terms of hyperbolic functions. Alternatively, if we suppose that Eq. (21) holds whenever c > 0, one can obtain an appropriate piecewise solution that only differs from Eq.
, which we explore in more detail in Appendix B. However, both of these options give rise to tumour dynamics that are essentially indistinct from those of Eq. (22). Hence, we will assume that Eq. (22) applies without further consideration, seeking simplicity in the analysis that follows and implicitly considering spheroids that are sufficiently small so as to justify this assumption. It is straightforward, but notationally cumbersome, to pursue our analysis with either of the alternative nutrient profiles and relax this assumption of smallness.

Governing Equations
For completeness, we now state the full system of equations governing the evolution of the solid tumour, incorporating each of the assumptions detailed above: along with the boundary and initial conditions Integrating Eq. (23a) in space, taking a Lagrangian time derivative, and changing integration variable yield the degenerate partial differential equation for r = r (R, t).

A Minimal Nutrient-Limited Growth Model
Our first model for tumour growth draws inspiration from the classical model of Greenspan for nutrient-limited growth. In Greenspan's model, the rate of growth is determined by the local nutrient concentration, and thresholds of nutrient concentration determine whether a tissue is classified as proliferating, quiescent, or necrotic. Adopting this principle leads to the minimal growth law whereĉ ∈ (0, c ∞ ) is a fixed nutrient threshold, below which tissues reduce in size due to lack of nutrient availability. This simple form captures the notion that, given greater nutrient availability, growth will be accelerated, whilst a lack of nutrient results in cell death and decay. In particular, we will define necrotic tissue via the nutrient threshold condition c <ĉ, whilst tissues with c ≥ĉ will be referred to as proliferating or growing. Of note, we have simplified Greenspan's original model by omitting a threshold for quiescence, instead distinguishing only growing and necrotic regimes.

Steady States of Growth
Inserting this growth law into our governing equations modifies Eq. (25) to the simple relation With c(r , t) given by Eq. (23d), this integral can be evaluated explicitly, yielding the temporal evolution equation for material points. In particular, taking R = B gives an explicit ODE for the outer radius b(t) of the spheroid, We may solve this equation to give the cumbersome but elementary explicit form valid for a positive initial radius B andĉ < c ∞ , illustrated in Fig. 3a. Alternatively, a direct analysis of Eq. (29) yields the steady states of the dynamics, at which db/dt = 0. These steady solutions are readily seen to be b = 0 and and it is straightforward to show that these states are linearly unstable and stable, respectively. Hence, so long as B = b(0) > 0, the tumour evolves to the state b = b * N , with growth limited by nutrient availability. The nonzero steady state features a necrotic core at the centre of the tumour, surrounded by a proliferative rim of tissue. The steady radius of the necrotic core, denoted by rĉ, can be computed as so that rĉ/b * N = 3/ √ 15 and the necrotic region occupies approximately 46% of the tumour volume, independent of the model parameters. In line with this analysis, the distribution of nutrient in the tumour at steady state is shown in Fig. 3b.

Material Turnover
Though we have described the spheroid as being at steady state when b = b * N , the tissue inside the tumour is far from idle. In particular, assuming that b(t) = b * N , the spatial ODE of Eq. (28) is trivially modified to which captures the non-steady dynamics of material within the spheroid. When viewed as an ordinary differential equation in time for r (R, t) for fixed R, this equation admits the same steady states as that for b(t), so that the steady states are simply r = 0 and r = b * N . However, for R ∈ [0, B), the linear stability of these stationary solutions is reversed, with r = b * N being unstable whilst r = 0 is stable. Hence, material points in the interior of the spheroid move away from the outer proliferating rim and towards the central necrotic region. This behaviour, which might be expected of nutrient-limited growth, is illustrated in Fig. 3a and motivates a feature of the numerical implementation later used to simulate the governing equations, as described in Appendix A.

The Growth of Solid Stress
Whilst the simple growth law of Eq. (26), which does not incorporate mechanics, produces a plausible growth curve, we now consider whether it predicts realistic mechanical stress. Numerical solution of the governing equations of solid stress is straightforward, and the details of our implementation are discussed in Appendix A. In Fig. 4 , we show an illustrative stress profile at an instant during growth. Hoop stresses are compressive at the proliferating boundary and become tensile towards the centre of the tumour, whilst the radial stress is tensile throughout, with κ = 0 in this figure.
Turning our attention to the evolution of the solid stress, we recall that, even at a steady state of b(t), there is a constant turnover of material within the tumour. In particular, the tissue at the outer boundary of the spheroid, which remains there throughout the dynamics, grows continuously, with c = c ∞ on r = b(t). Noting the simplicity that this boundary condition affords, we can explicitly write the evolution of the growth stretch at the boundary as This then yields the elastic stretch Fig. 4 Accumulated stresses in a nutrient-driven spheroid. The radial and hoop stresses σ r and σ θ at an instant during growth are plotted as a function of r , shown as solid and dashed curves, respectively. Tensile radial stress increases in magnitude towards the centre of the spheroid, whilst the hoop stress is compressive at the outer boundary of the tumour and becomes tensile at the core.
Here, we have adopted the parameters of Fig. 3, taken κ = 0, and sampled at time t/T = 10 which in turn allows us to evaluate the derivative of σ r via Eq. (14) as for positive constants C 1 and C 2 . Perhaps surprisingly, the derivative of radial stress eventually increases in magnitude exponentially with time, driven by the exponential growth of the spheroid at the boundary, even at a steady state of the spheroid radius. Further, noting from Eq. (16) that σ r is constant on the boundary at a steady state of b(t) and that σ θ is a linear combination of σ r and ∂σ r /∂r , we can also conclude that σ θ is eventually compressive and grows exponentially in time. Specifically, as t → ∞. Hence, this model predicts that the solid stress in the tumour is unbounded, growing exponentially. Numerically, we observe similar behaviour in the radial stress.
The existence of this behaviour, which is confirmed to be present across parameter regimes via numerical simulations, lends itself to two distinct interpretations. One sees this prediction of ever-accumulating stress as capturing the phenomenon of spheroid shedding, whereby material is seen to break off from a tumour accompanying the destabilisation of the spherical structure (Giverso and Ciarletta 2016). In this context, one might interpret the unlimited stresses as indicative of a symmetrybreaking or topology-changing instability, though such events are beyond the reach of our framework.
Alternatively, the growing stresses might be seen instead as an unphysical consequence of our modelling assumptions. We adopt such a viewpoint, seeking not to overreach in the interpretation of this model and its emergent dynamics. Hence, we will explore alternative growth laws that replace the minimal form of Eq. (26), with our goal being to preclude the generation of unbounded stress.

A Modified Growth Law
Motivated by experimental evidence of stress-mediated regulation of cell proliferation (Delarue et al. 2014;Helmlinger et al. 1997), and building on previous works in which stress has been incorporated into the regulation of spheroid growth (Ambrosi et al. 2012(Ambrosi et al. , 2017Ciarletta et al. 2013;Ambrosi and Mollica 2004), we now couple the growth dynamics to stress. Explicitly, we pose where n is an as-yet-undefined non-negative function of σ that encodes the effects of stress on growth. This form captures the intuitive principle that stress may modify the dynamics of growing tissues, whilst the decay of necrotic material is unaltered and consistent with Greenspan's assumptions. It is unclear how the stress modifier n should depend on the stresses experienced by the tumour: should it be a function of the radial stress, the hoop stress, or some other measure? Here, recalling that we impose a condition on the radial stress at the boundary and as a first exploration, we will suppose that n = n(σ r ), so that growth is affected by the local radial stress, though we later explore an alternative.
In specifying the functional form of n, we note that, unless the infimum of n is zero, the stress accumulation argument of the previous section would hold with minor modification, with this growth law therefore also leading to unbounded stresses at steady state. Seeking to avoid such an unphysical phenomenon, we consider whereσ ≤ 0 is a threshold parameter. This piecewise linear function, which is weakly increasing in σ r , prohibits local growth when σ r is sufficiently negative but has no effect when σ r is non-negative. This amounts to the hypothesis that compressive stresses restrict growth, with tensile stresses having no similar effect, qualitatively

Impact of Stress-Limited Growth
By design, this law prohibits the growth of material that is under too much radial compression. At first glance, this might appear to solve the problem of unbounded stresses, with our argument for exponentially growing stress at the boundary no longer applying. However, with solid stress being an inherently non-local quantity, we will see that the locality of our growth law still allows for stress accumulation past any given threshold. Indeed, we exemplify such stress-limited dynamics in Fig. 5, takinĝ σ /μ = −100 and incorporating the compressive effects of the surrounding medium The growth curve of a tumour where growth is restricted by the local radial stress, with the outer radius b(t) and the boundary of the necrotic core shown as solid and dotted curves, respectively. The tumour radius initially appears to saturate, though experiences an increase in growth rate around t/T = 100, qualitatively distinct from the nutrient-driven model of Fig. 3. b The radial stress σ r at the centre of the tumour is shown as a function of time, from which an approximately exponential accumulation of stress is apparent, despite our stress-limited growth law. Inset is the dynamics at early times, with the stress initially compressive. c The spheroid composition at t/T = 200, shaded by growth rate, highlighting a large quiescent rim of tissue whose growth has been arrested by accumulated radial stress. Beneath this rim is a region of proliferation, shaded red, with a decaying necrotic core being present inside the grey dashed curve. Here,ĉ/c ∞ = 4/5, κ/μ = 316.2,σ /μ = −100, and B = L by taking κ/μ = 100. The configuration of the spheroid at the final simulated time is shown in Fig. 5c, shaded by growth rate, from which we note the presence of a quiescent outer rim of tissue whose growth has been arrested by compressive radial stress, in line with our posited growth law. However, Fig. 5b highlights the rapid accumulation of solid radial stress in the interior of the spheroid, with the now-internal proliferative rim driving material turnover into the necrotic core. Hence, in spite of our stress-limited framework, growing stresses remain a realisable and undesirable behaviour. This example also showcases another unrealistic consequence of this growth law. In particular, focusing on the growth curve of Fig. 5a, we note that the dynamics near saturation are not reminiscent of typical growth profiles, with the rate of tumour growth visibly increasing, rather than decreasing, around t/T = 100. These dynamics are due to the relaxation of accumulated stress in the proliferating rim of the tumour, the latter driven by the decay of the necrotic core and causing the corresponding increase in the stress-modulated growth rate. This leads us to question our treatment of necrosis.

A Modified, History-Dependent Growth Law
A basic assumption of our growth laws thus far has been the decay of necrotic tissue, drawing inspiration from Greenspan's classical work. However, with no precise notion of material clearance present in our model, it is not clear that having such a sink of tumour mass is appropriate. The shrinkage of material also appears to be partly driving stress accumulation in the spheroid, along with the generation of atypical growth curves, as we saw in Sect. 3.2. Hence, recalling our assumption of incompressibility, we will adopt an alternative approach to necrosis, supposing instead that necrotic material simply ceases to proliferate, a condition that is permanent. This model of Fig. 6 Growth of a stress-limited spheroid with persistent necrotic tissue. a The tumour growth curve and the necrotic radius r (R N (t), t). b The stress-driven saturation of radial stress at the centre of the spheroid. c The composition of the tumour and the distribution of radial stress as it approaches steady state. The spheroid consists of a central necrotic region surrounded by a wide rim of stress-arrested tissue The dashed circle marks the current boundary of the necrotic region, with radius r (R N ) at steady state. The approximately uniform radial stress distribution entails that the hoop stress is approximately equal to the radial stress, with the two quantities being indistinguishable at the resolution of these plots. Here,ĉ/c ∞ = 9/10, κ/μ = 1, σ /μ = −100, and B = L/100 nutrient-starved tissues also overcomes a potential limitation of Greenspan's approach, which allows necrotic tissue to return to a proliferating state. Concretely, we pose where R N (t) denotes the radius of the necrotic core of the spheroid. This law prevents the growth or decay of necrotic tissue, whilst leaving the dynamics of perfused tissue unaltered from the stress-dependent law of Sect. 3.2.1. Formally, the Lagrangian quantity R N is defined via the somewhat cumbersome expression with the intuitive interpretation that R N (t) is non-decreasing and bounded below by Greenspan's nutrient-determined necrotic radius at time t. This weakly increasing quantity captures the desired permanence of necrosis, whilst incorporating the principle of Greenspan's threshold-based definition. The growth law of Eq. (40) has an immediate consequence: the growth rate is always non-negative. Hence, the growth rate must vanish everywhere in order for the tumour to attain a steady state. This contrasts our previous explorations, which had the potential to admit a steady state of tumour radius where proliferation was balanced by the supposed decay of necrotic material. The absence of such a behaviour in our modified model results directly from the persistence of necrotic tissue, with no mechanism of material clearance being present in this simple model. However, a steady state is still attainable, resulting from the total arrest of nutrient-rich tissues by the imposed external stress. Sample growth dynamics corresponding to this model are presented in Fig. 6 and exemplify the desired exponential-linear-saturating growth curve. Further, the stress evolution shown in Fig. 6b demonstrates the saturation of Fig. 7 Stress-modulated tumour growth without a steady state. a A growth curve corresponding to a model tumour whose growth, whilst limited in theory by compressive radial stress, is accelerating after a period of apparent saturation. b The stress profile at t/T = 450, highlighting a region of rapid transition between high-magnitude stresses in the core and compressive radial stresses at the boundary. c The composition of the tumour at t/T = 450, with the narrow but persistent region of proliferation shown red, shaded by growth rate. The inner grey region corresponds to the necrotic core, whilst the outermost rim of the tumour has been arrested due to the compressive stress from the external medium. Here,ĉ/c ∞ = 4/5, κ/μ = 102.4, σ /μ = −100, and B = L mechanical stress within the spheroid as it approaches steady state, with the entire spheroid becoming quiescent as t → ∞, as illustrated in Fig. 6c.

The Potential for Unbounded Evolution
Though our refined model appears promising, the lack of a decay-driven steady state raises a new issue: it is not guaranteed that there is a non-negative steady state in all parameter regimes. Indeed, we realise such a case in Fig. 7, where the evolution of the tumour radius deviates from the saturating profile that is typical of tumour spheroids, instead growing indefinitely. In particular, after what appears to be the onset of saturation, the tumour experiences an increase in growth rate, with a thin proliferating band of tissue growing within an outer quiescent rim. This appears to be the result of tensile radial stresses building up inside the tumour at a faster rate than can be compensated for by the linearly compressive boundary condition, yielding a narrow non-vanishing region between the necrotic core and the stress-arrested tissue where the growth rate is not identically zero. Such a persisting region of mild solid stress is visible in Fig. 7c, with the corresponding radial and hoop stresses being illustrated in Fig. 7b. To remove this possibility and guarantee the existence of a nonzero steady state in all parameter regimes, as per the desired properties set out in Table 1, we further modify our growth law, focusing on the consequences of locality.

Recovering a Robust Steady State
With the previous stress-dependent spheroid model having forgone any guarantee of a steady state in favour of stress-dependent growth and the persistence of necrotic tissue, it remains to recover the desired growth curve in all parameter regimes where external resistance is applied. In modifying our growth law appropriately, it is instructive to further consider the stress distribution in the ever-growing tumour of Fig. 7, as illustrated in Fig. 7b. In particular, we note that the compressive radial stress arrests growth only in the outer portion of the tumour, where σ r <σ , with the rest of spheroid able to grow if sufficiently perfused with nutrient. It is this locality of the stress modulation that allows for this varied composition. Hence, we will modify the local nature of our stress dependence, instead modulating the growth rate by a non-local measure of stress. With reference to the biological tissues that we seek to model, such a non-local effect may be interpreted as the result of inter-cell signalling (Maia et al. 2018;Aasen et al. 2016). In particular, noting that the lengthscale of the proliferating region of tissue is dictated by the diffusion of nutrients, the diffusion of signalling molecules represents a plausible mechanism for cell-cell communication within the well-perfused regions of the spheroid.
With this interpretation in mind, we propose the following non-local growth law: Here, in place of the local radial stress, we have adopted a measure that is both global, in that it accounts for the stress throughout the tumour, and directionally unbiased, with minR{σ r (R, t), σ θ (R, t)} being the largest magnitude compressive stress experienced by the tissue in any direction. This latter property arises as σ r and σ θ are the eigenvalues of the (diagonal) stress tensor. In practice, with reference to the composition of the spheroid shown in Fig. 7, this growth law prevents the formation of narrow bands of proliferating tissue within an arrested outer rim. Equation (42) is identical in structure to that of Sect. 3.3 and, hence, inherits the desirable properties of each of our considered models. In particular, it maintains a nutrient dependence reminiscent of Greenspan's seminal work, modified to consider a persistent core of necrotic tissue in the absence of material clearance. This law also enables solid stress to regulate growth within the tumour, whilst also guaranteeing the existence of a steady state. This latter assertion warrants justification. As the rate of growth under our new law is non-negative, r (R, t) is weakly increasing in time for all material points, which follows immediately from Eq. (25). In particular, the outer radius of the tumour is weakly increasing in time, with its rate of change being zero precisely at a steady state. Supposing that the tumour grows in a resistive external medium, so that κ > 0, the compressive boundary condition of Eq. (16) implies that σ r at the boundary decreases as the spheroid radius increases. Hence, assuming that the tumour does not attain a steady state by another means, the radial stress at the boundary necessarily approaches the thresholdσ . Thus, owing to the now-global dependence of the growth rate on the radial solid stress, tumour growth arrests throughout the entire spheroid, so that a steady state necessarily exists. This argument also places an upper bound on the tumour radius that at which σ r (B, t) =σ , as specified by the compressive boundary condition of Eq. (16). Of note, this reasoning applies to any form of radial stress boundary condition, subject to the assumption that it is strictly decreasing and unbounded below in the  curves (a, d, g), the evolution of radial stress at the centre of the tumour (b, e, h), and the distribution of radial stress at the final simulated timepoint (c, f, i). In all cases, we observe qualitatively plausible saturating growth curves, accompanied by saturating radial and hoop stresses (not shown) tumour radius. In the absence of a compressive external stress, i.e. κ = 0, the above argument does not apply, though an alternative argument involving the hoop stress allows us to partially recover the guarantee of a steady state in this case, discussed briefly in Appendix C.
In order to exemplify the robustness of our proposed model and that it generates qualitatively plausible growth dynamics, we present a number of simulated spheroids in Fig. 8, adopting the parameters of Figs. 5, 6 and 7 in turn, two of which previously gave rise to undesirable growth curves or ceaselessly accumulating stresses.

Generating Plausible Residual Stresses
The above model satisfies all but one of the criteria set out in the Introduction, giving rise to plausible growth curves in all parameter regimes and being free from singular or inadmissible behaviours. However, we have not yet considered the plausibility of the generated stress profiles, except for guaranteeing that they remain bounded. In particular, we now seek to augment Eq. (42) in order to qualitatively match experimentally determined profiles of residual stress, i.e. stress in the absence of external loads. In particular, the radial component σ R r of the residual Cauchy stress tensor σ R satisfies whilst the corresponding residual hoop stress σ R θ satisfies In a theoretical study inspired by experimental works, it has been predicted that the residual hoop stresses in tumours may be tensile (σ R θ > 0) at the boundary of the spheroid and compressive (σ R θ < 0) further inside the tumour (Stylianopoulos et al. 2012). When such a spheroid is cut radially, the relaxation of these stresses leads to an opening of the tumour, measurements of which lead to estimates for the residual hoop stresses (Stylianopoulos et al. 2012;Ambrosi et al. 2017;Guillaume et al. 2019). As the final goal of this tutorial, we seek to replicate key features of this qualitative profile of residual hoop stress, specifically that the residual hoop stress can be compressive within the tumour whilst also satisfying σ R θ > 0 at the surface. To make progress, we note that the growth stretch associated with this model increases monotonically in R for all times t, a property that it inherits from the nutrient profile, so that ∂γ /∂ R ≥ 0 everywhere. This immediately implies r ≤ γ R, which may be deduced via the integral formula The first equality in Eq. (45) results from the integration of Eq. (23a) in R, whilst the inequality stems from the monotonicity of γ in R. Hence, Eq. (43) gives ∂σ R r /∂ R ≤ 0, so that ∂σ R r /∂r ≤ 0 also. Separately, noting Eq. (44) and the stress-free boundary condition on σ R r at R = B, we can identify σ R θ = r (∂σ R r /∂r )/2 at R = B. Hence, the residual hoop stress σ R θ has the same sign as the spatial derivative of σ R r . Thus, we have that σ R θ ≤ 0 at the boundary, so that the residual hoop stress at the surface of the spheroid is never tensile in the model of Sect. 3.4.
From this analysis, it is clear that the monotonicity of γ poses a barrier to generating a realistic profile of residual stress. However, this guarantee of monotonicity was not present in some of our earlier models, in particular those of Sects. 3.2 and 3.3, which each employed a coupling of the growth rate to the local stress. Motivated by the rapidly varying profiles of hoop stress shown in Fig. 7b, we now reintroduce a local stress coupling to our growth law.
In order to retain the many desirable properties of our previous model, we augment the robust growth law of Sect. 3.4 to Fig. 9 Accumulated residual stresses and growth stretches in a tumour with both local and non-local stress dependence, following the model of Sect. 3.5. a The residual stress profile at steady state. At and near the boundary of the spheroid, the residual hoop stress is tensile, in line with experimental observations, and becomes compressive further inside the tumour. Accordingly, we see that the residual radial stress is increasing towards the boundary. b The growth stretches at steady state, whose non-monotonicity near the boundary is associated with the tensile residual hoop stresses shown in a. Here,ĉ/c ∞ = 1/4, κ/μ = 0.1, σ /μ = −1, β = 6.25, and B = L whereñ couples the growth rate to the locally experienced stress, analogous to n of Sect. 3.2. This growth law thereby captures both local and non-local stress responses, in addition to enabling nutrient availability to regulate growth. Further, and by design, this growth law retains the robustness properties of the previous model whilst the local stress term allows for nuanced and, notably, non-monotonic growth stretches, as we demonstrate below. As with n, there are countless choices of the functionñ, including its argument. Here, we have takenñ = n(βσ r ) for a parameter β > 0 that determines the relative sensitivity of the local stress response compared to the global term, so that the local response depends on the experienced radial stress, emphasising that the local and global stress responses may be wholly distinct in form. Even in this minimal case, the introduction of local stress dependence can give rise to profiles of residual stress with the desired features, though we note that such a profile is not guaranteed by this form of growth law. As an example, in Fig. 9 a we illustrate the residual stress profiles in a simulated tumour at steady state, which displays tensile hoop stresses on the boundary and a compressive region within the tumour. Commensurate with the stress at the boundary is the non-monotonic profile of the growth stretch, as shown in Fig. 9b.

Summary and Evaluation
Having constructed a model that satisfies the criteria set out in the Introduction, in Table 2 we summarise the iterative process that led to this final spheroid model. In this table, we identify each model by the subsection of Sect. 3 in which it was introduced, and highlight the relevant features.

Table 2
Summarising model properties. Summary of the models in this study and their characteristics; a dash denotes that a property has not been explored. Models are identified by the subsection of Sect. 3 in which they were introduced. Entries annotated with a are conditional: models 4 and 5 have a guaranteed steady state except in a subcase of free suspension, whilst plausible stress profiles are realisable, but not guaranteed, for  Helmlinger et al. (1997). The experimentally observed tumour growth dynamics described in Fig. 1 are shown in terms of the spheroid radius alongside fitted dynamics of the model of Sect. 3.4, with the latter displayed as solid curves. The good agreement between the experimental data and the model fits highlights the ability of the constructed model to capture the phenomenon of mechanically influenced growth. In particular, we have varied only the initial conditions and external stiffness parameters between datasets during fitting the model, so that differences in dynamics can be directly attributed to these factors As a final evaluation of the model of Sect. 3.5, we return to the experimental data of Helmlinger et al. (1997) illustrated in Fig. 1. Recalling the inability of the minimal classical model of Greenspan to capture the evidenced mechanical influences on tumour growth, we now fit our final model to the three displayed growth curves simultaneously, as described in Appendix A. During this fitting process, we vary only the initial condition B and external stiffness parameter κ between the three growth curves, with all other model parameters shared. The fitted dynamics displayed in Fig. 10 show good agreement between the model and experimental data. This agreement highlights the ability of our final model to capture the range of mechanically driven behaviours observed by Helmlinger et al. (1997), with these behaviours being directly linked to the mechanical parameters of our framework.

Discussion
In this tutorial, we have explored a sequence of models of solid spheroid growth, ranging from simple nutrient-limited growth dynamics to more complex stress-dependent growth laws. Inspired by the classical model of Greenspan (1972), we began our exploration of tumour development by assuming a minimal threshold-based approach to determine the growth rates of tissues within the spheroid, distinguishing between proliferating and necrotic cells via the local concentration of an abstracted nutrient. The resulting growth law, when coupled to the solid framework of morphoelasticity, enabled a thorough theoretical analysis of the emerging dynamics, with simplicity affording significant tractability in this case. However, this analysis uncovered an unavoidable and unphysical behaviour in the tumour dynamics, with solid stress increasing unboundedly, even at steady states of spheroid size. Thus, we conclude that the principles of Greenspan's approach are not well-suited to the modelling of solid spheroids, at least in the context of the considered framework and without appropriate refinement.
Following on from our analysis of a Greenspan-inspired growth law, we considered a range of refinements and modifications, introducing a dependence on both local and non-local solid stress and modifying the treatment of necrotic tissue. We demonstrated several unexpected consequences of our modifications, such as the limitless accumulation of stress in a model where stress limits and even halts tissue growth. This example particularly highlights how the non-locality of mechanical stress can complicate model analysis, with stress building up in non-growing regions of the spheroid due to growth elsewhere in the tumour. Nevertheless, despite the presence of this mechanical complication, we have been able to concretely reason about the behaviours of our final two models of tumour growth, concluding in both cases that the emergent dynamics necessarily reach a steady state of tumour size and solid stress, except for a subset of dynamics for growth in free suspension. This analysis highlights the benefits of employing minimal models, especially in the context of solid mechanics, with this reasoning rendered tractable by the simplicity of the morphoelastic framework in a spherical geometry.
A further benefit of the modelling framework explored in this tutorial is the readiness with which we have been able to explore and experiment with different growth laws. In particular, even when analysis has not been tractable, numerical solution has been straightforward, with only the solid stresses necessitating care in order to integrate a removable singularity. In turn, this has enabled a thorough consideration of the consequences of employing solid stress as a regulator of tumour growth, from which we have seen that stress appears feasible as a factor that affects spheroid development, in support of previous experimental and theoretical works (Helmlinger et al. 1997;Delarue et al. 2014;Ambrosi and Mollica 2004). In particular, we have seen that a combination of local and non-local stress dependencies can give rise to robustly reasonable profiles of growth and stress, in qualitative agreement with observed tumour dynamics and estimated residual stresses. Hence, our exploration supports the broad hypothesis that tumour-environment interactions, in the particular form of mechanical stress, can be an influential and nuanced determinant of spheroid growth.
Complimentary to the analysis and exploration presented here, there are numerous directions for future modifications to the modelling framework and for the assessment of the role of the environment in tumour growth dynamics. One such avenue includes the relaxation of the strict assumption of spherical symmetry, which may be of pertinence to the stability of model spheroids that are highly stressed, such as those encountered in Sect. 3.1. Alternatively, there is extensive scope for the refinement of our treatment of the tumour as a single solid phase, with the potential to broaden to a poroelastic framework, as in Ambrosi et al. (2017), or to a more general multiphase model. In particular, such an extension might include an alternative treatment of necrotic material and the inclusion of material clearance, with an explicit mechanism absent from our growth-focused exploration. Further, there is also the prospect of establishing quantitative agreement between the described models and additional experimental datasets.
In summary, we have posed and explored a hierarchy of mechanical tumour spheroid models, exploring and iterating upon simple growth laws that draw from classical study and experimental observations. In doing so, we have seen how even simple models can give rise to unexpected and unphysical behaviours when cast in the context of solid spheroids. Seeking to preclude such eventualities, we have explored how solid stress can be used to regulate growth in phenomenological models, considering various couplings of growth to stress and evidencing the plausibility of the mechanical environment as a driver of spheroid development. This sequence of modifications has culminated in a simple yet robust model of tumour growth that shows favourable agreement with experimental data, highlighting the potential for simple models to be more than just toy mathematical examples.

A.1: Adaptive Discretisation
The analysis of Sect. 3.1 and Appendix B predicts that material points will move towards the centre of nutrient-driven spheroids with a decaying necrotic core, with the exception of those precisely on the surface of the tumour. In the context of a computational implementation, this entails that any fixed discretisation of the Lagrangian domain will eventually become unsuited to capturing the entire Eulerian domain with sufficient resolution. This is illustrated in both Figs. 3a and 12 a, with r (R, t) → 0 as t → ∞ for R ∈ [0, B). To overcome this numerical issue, we rediscretise the Lagrangian domain at each timestep so that it corresponds to a uniform discretisation of the Eulerian domain, noting that r (R, t) is a bijective mapping between the two domains at any time t.
However, this remeshing procedure results in the discrete points of the Lagrangian domain being clustered close to R = B, to the point that these quantities rapidly become indistinguishable when using double precision arithmetic. We mitigate this numerical limitation by casting the governing equations of Eq. (23) in terms of a shifted Lagrangian variable,R := R − B, with values close to zero being better distinguished in fixed-precision arithmetic than those near B > 0.

A.2: Computing Radial Stress
Equation (23b), which gives the Lagrangian derivative of σ r , has a removable singularity at r = R = 0 that presents a significant numerical challenge. To circumvent this when integrating Eq. (23b) numerically via the trapezium rule, we replace this expression with its two-term Taylor expansion about R = 0 when R/B < 1/20. This expansion is calculated by first computing higher-order expansions of r (R, t) and γ (R, t), noting once more that from Eq. (23a). Writing the expansion of γ with respect to R-derivatives of γ , after significant calculation this leads to 2μγ as R → 0. Here, a prime denotes a derivative with respect to R, with γ and its derivatives being evaluated at R = 0. Derivatives of γ are estimated numerically with a second-order finite-difference scheme. The validity of this expansion and the resulting numerical scheme is confirmed by comparison with numerical solutions computed with high-resolution discretisations. to steady state is shown as a solid black curve, alongside the paths of material points, which are shown as solid grey curves and correspond to R ∈ [0, (1 − 10 −5 )B]. As in Sect. 3.1, after initial growth, material points move away from the steady outer edge of the spheroid and towards the necrotic core. The analytically predicted steady state, which here is only a leading-order approximation, is shown as a dashed black line, from which we note remarkable agreement between the numerical and asymptotic solutions even when b(t) is not large. The radius at which c(r , t) =ĉ is shown as a thin dotted curve. b A slice through the centre of the spheroid at steady state, shaded according to the nutrient concentration c, which highlights a key difference between these partially perfused dynamics and the nutrient-rich evolution of Fig. 3. The threshold for necrosis,ĉ, is shown as a dashed white curve, inside which the tissue shrinks in response to lack of nutrient. Here, we have taken parameters such thatĉ/c ∞ = 2/5 and B = L. For these parameter choices, the leading-order approximation to the steady state is b * N ∼ 5L/ √ 2 where c I denotes the solution of Eq. (22) in the main text. This elementary but notationally cumbersome solution is illustrated in Fig. 11.

B.1: Steady States of Growth
In the partially perfused regime, adopting the minimal nutrient-driven growth law of Sect. 3.1, Eq. (27) is modified to For r <r (t), recalling the expression for c II of Eq. (B3), it is clear that this spatial ODE results in movement towards the centre of the spheroid, whilst, for r >r (t), there is the possibility of a nonzero steady state, with growth balancing out necrotic decay. For given system parameters, this steady state is straightforward to compute numerically, though presents a challenge for analytical solution due to the presence of the threshold radiusr (t), which depends on b(t). However, we can make simple progress towards a steady-state solution for b(t) under the assumption that the tumour is large, by which we mean that b(t) is larger than the other lengthscales in the problem.
To see this, we first note that b −r represents the distance beyond which the nutrient at the boundary can no longer penetrate the tissue, a quantity that should be approximately independent of the size of the tumour when the spheroid is large. Hence, posing the asymptotic ansatzr (t) ∼ b(t) −r 0 + o(1) as b(t) → ∞ for some as-yet-undetermined constantr 0 , solving for the roots of Eq. (B4) gives to leading order, which simply givesr 0 = √ 2Dc ∞ /λ = √ 2L. Hence, for large spheroids, we have thatr (t) ∼ b(t) − √ 2L. Inserting this approximate expression into Eq. (B3), setting r = B(t), and evaluating the integral of Eq. (B6) yield which has leading-order steady states b = 0 and redefining the nutrient-limited steady state in this context for notational convenience. It is clear that only the nonzero steady state can be consistent with the asymptotic analysis, so that the only admissible steady state scales with c ∞ /ĉ, which may readily become large for sufficiently low thresholds for tissue necrosis. Further, inspecting the form of the approximate dynamics highlights that this steady state is linearly stable, consistent with the negative feedback present between tumour size and nutrient availability within the spheroid, the balance of which gives rise to this nutrient-limited steady state. This steady-state solution also gives the approximate steady-state value ofr (t), which we see is related to the system parameters simply asr ∼ √ 2L(c ∞ /ĉ − 1). Approximating the radius of the necrotic core byr , valid in the large-b * N limit, the fraction of the tumour volume occupied by the necrotic core is approximately Perhaps counter-intuitively, a decrease in the necrosis threshold results in an increase in the proportion of necrotic tissue at steady state. However, this is easily reconciled upon noting that the tumour radius at steady state is greatly increased by the change in this parameter, whilst the perfused rim of the spheroid remains at an approximately constant thickness. In order to illustrate the validity of our approximate solution for the steady dynamics, particularly the tumour radius at steady state, we numerically solve the exact ODEs governing the evolution of material points in this partially perfused tumour, showcasing a sample exploration in Fig. 12. Remarkably, even for the b(t) < 4L attained in this example, we see very good agreement between the asymptotically predicted steady state b * N and the numerically computed evolution, with agreement naturally improving in parameter regimes with larger steady states (i.e. those with reduced c/c ∞ ). From Fig. 12b, we note the distinct composition of the spheroid at steady state when compared to the nutrient-rich tumour of Fig. 3, with a large region devoid of Fig. 13 Nonzero steady states of nutrient-driven spheroid growth. Normalised by the diffusive lengthscale L, the numerically computed steady state b * N of spheroid growth in the absence of mechanical effects is shown as a dashed black curve as a function ofĉ/c ∞ . The analytically predicted steady states are shown as solid grey curves, with the exact result of the analysis of Sect. 3.1 shown forĉ/c ∞ ≥ 3/5 and the asymptotic result of Appendix B shown forĉ/c ∞ < 3/5. Though derived for large steady states, the asymptotic result displays a remarkable accuracy even for states of moderate magnitude. Spheroid dynamics were simulated until a time t/T = 100 with B = L, at which point the dynamics had converged to the steady state at the resolution of this Figure nutrient present in this case and the necrotic region occupying a larger proportion of the spheroid.
Additionally, with the identified steady states having been functions of only the ratiô c/c ∞ and the diffusive lengthscale L, we illustrate the normalised nonzero steady states across the admissible range ofĉ/c ∞ in Fig. 13, including the asymptotic results of this section, the exact results of Sect. 3.1, and numerically computed steady states. From this, we observe precise agreement between the numerical results and the expression of Sect. 3.1, valid for b ≤b, whilst the asymptotic results of this section are seen to retain substantial accuracy even for moderately sized steady states, which lie outside the regime of asymptotic validity.

B.2: Material Turnover
Figure 12a also highlights the motion of material points within the spheroid, which approach the centre of the necrotic core of the tumour as the spheroid approaches steady state, as in the case of Sect. 3.1. Here, we can make simple progress in determining the evolution of much of the material in the tumour at steady state, noting that the majority of the tissue lies in the central region where c = 0, with r (R, t) ≤r (t). In this region, the evolution of material points is governed by the simple spatial ordinary differential equation whose solution is exponential decay towards r = 0 at a rate kĉ. For material points lying above the thresholdr , the analysis is once again complicated by the presence of ther term implicitly defined by Eq. (B6). However, for large r , one can once again make asymptotic progress, leading to the conclusion that r (R, t) also has a steady solution r = b * N , though this state is linearly unstable for R ∈ [0, B), which mirrors the analysis of the simpler, nutrient-rich case of Sect. 3.1.

Appendix C: Recovering a Steady State Without External Compression
Consider the model of Sect. 3.4 when κ = 0, so that no external compression is applied. Suppose that no steady state is reached and that n of Eq. (42) attains a strictly positive minimum during the dynamics, so that n ≥ n m > 0 for some fixed n m and σ r and σ θ are bounded below. Then, bounding and integrating the growth law of Eq. (42) evaluated at R = B gives the inequality where we are mirroring the analysis of Sect. 3.1.3. Equation (14) then allows us to bound the derivative of radial stress as where D 1 and D 2 are positive constants. Substituting this inequality into Eq. (18) and noting that σ r | R=B = 0 when κ = 0, we have As the growth of large spheroids is limited by the diffusive lengthscale L, it is clear that b(t) grows sub-exponentially in time for large spheroids. Hence, recalling that we are assuming that growth is unbounded, the hoop stress at the boundary is dominated by the growing exponential in the second term of Eq. (C15), so that σ θ | R=B → −∞ and, consequently, n becomes zero at finite time. This contradicts our supposition, so that the tumour must evolve to a steady state, subject to the caveat that our argument does not apply to dynamics where n → 0 but n = 0. This argument can be trivially modified to apply to any compact interval of time, which guarantees the existence of n m > 0 if n is assumed to be continuous. By evaluating the bound of Eq. (C15) for given parameters, noting that the stress must lie aboveσ in order to avoid a steady state, this exponentially restricts the values of n m needed to generate perpetually unsteady dynamics. Whilst it may be possible to construct a functional form forσ so as to satisfy these bounds and generate dynamics without a steady state, we do not expect this to occur regularly in practice due to the exponential time dependence of Eq. (C15). Indeed, we invariably obtain steady states as the numerical solutions to the model of Sect. 3.4 when κ = 0 across a wide range of parameter values, in support of this claim.