Viscocapillary Instability in Cellular Spheroids

We describe a viscocapillary instability that can perturb the spherical symmetry of cellular aggregates in culture, also called multicellular spheroids. In the condition where the cells constituting the spheroid get their necessary metabolites from the immediate, outer microenvironment, a permanent cell flow exists within the spheroid from its outer rim where cells divide toward its core where they die. A perturbation of the spherical symmetry induces viscous shear stresses within the tissue that can destabilise the aggregate. The proposed instability is viscocapillary in nature and does not rely on external heterogeneities, such as a pre-existing pattern of blood vessels or the presence of a substrate on which the cells can exert pulling forces. It arises for sufficiently large cell-cell adhesion strengths, cell-renewal rates, and metabolite supplies, as described by our model parameters. Since multicellular spheroids in culture are good model systems of small, avascular tumours, mimicking the metabolite concentration gradients found in vivo, we can speculate that our description applies to microtumour instabilities in cancer progression.


Introduction
Interface instabilities in systems driven far from equilibrium have been extensively studied in solid-state physics. Classical examples are the Saffman-Taylor instability, which occurs when a fluid of lower viscosity displaces a more viscous one in a Hele-Shaw cell [1,2], the Mullins-Sekerka instability, which stems from the diffusive transport of the latent heat of solidification in unidirectional solidification [3,4], and the Rayleigh-Taylor instability, corresponding to the fingering of an interface between two immiscible fluids of different densities when the heavier fluid is placed on top of the lighter [5,6].
Instabilities originating from similar coupling terms as those responsible for these classical condensed-matter instabilities have been identified in living systems. In tissues or bacterial colonies, growth and cell divisions may give rise to similar or new out-of-equilibrium phenomena [7][8][9][10]. For example, in the case of bacterial-colony growths, patterns similar to those associated with aggregation phenomena and viscous fingering have been observed [11,12]. Such a coupling can lead to fractal branching patterns via the process of diffusion-limited aggregation [13,14] or other types of branching patterns, depending on the bacterial morphotype [12]. In tissues, mechanical instabilities have been recognised to play a potential role in different morphological processes and patterns exhibited by growing cell populations [15]. Examples are the wrinkling patterns of growing, soft surfaces [15,16] such as those of leaves and flowers [17,18], the large-scale looping morphology of the gut [19][20][21] and the generation of its surface villi [21,22], or the formation of cortical convolutions [23]. Such instabilities may emerge from a buckling phenomenon [24,25], potentially driven by the differential growth of adjacent tissue layers [18,[26][27][28], or by other curling and crumpling instabilities due to anisotropic growth [17].
Here, we focus on the stability of the spherical growth of cellular aggregates. Such experimental systems are used to study the growth dynamics and cellular structure of microscopic tumours with realistic metabolite concentration gradients [29]. They have also been used as anti-cancer therapy test platforms, mirroring the three dimensional cellular context and therapeutically relevant pathophysiological gradients Figure 1. Schematic drawing of the proposed viscocapillary instability in a multicellular spheroid. The original spherical shape is indicated by the external blue-dashed circle. The scaled metabolite concentration is depicted as a scaled colour gradient. Dividing cells (red) are located mostly close to the outer surface. Close to the centre, cells are deprived of metabolites and die (yellow dashed). An inner blue circle indicates cells equidistant from the spheroid centre, and the tissue-flow field is qualitatively depicted by white (resp. black) arrows for inward (resp. outward) cell velocities. At point A, cells are pushed inward faster than at point B, because more layers of dividing cells push the cells in A towards the spheroid core. This creates viscous shear stresses that, in reaction, push the cells located above A outwards. Capillary effects due to tissue surface tension at the outer boundary, however, favour stability. Depending on parameter values, the original spherical shape can or not become unstable by this viscocapillary mechanism.
of in-vivo tumours [30,31]. Cellular aggregates also permit the study of the effects of different perturbations on the growth or cellular-duplication dynamics, such as a change in the external mechanical constraints [32][33][34][35][36]. Most often, an effective surface tension exists between the aggregate and its direct environment, which makes it look like a spheroid, justifying the denomination of multicellular spheroid [37][38][39][40][41]. The supply of metabolites from the microenvironment is responsible for an inhomogeneous distribution of cell divisions within the spheroid, with an increased cell proliferation at its periphery and an increased cell death in its core [30]. As a consequence, cells flow from the outer rim towards the centre of the spheroid [34,42]. Under some circumstances, the spheroid reaches a steady-state size [34,43,44].
In the present paper, we investigate theoretically whether this steady-state spherical shape can be unstable without changing the average spheroid size, due to the shear stresses created by differential inward-directed cell flows. Similar shape instabilities have been proposed in the context of the growth of small, spherical tumours, starting with the seminal work of Greenspan [45]. These instabilities can be described in terms of reaction-diffusion processes, driven by the transport of growth-affecting factors. Some factors such as nutrients, acidity, or chemotherapeutic agents diffuse from external sources to the tumour cells [46][47][48], while others are produced by the tumour cells themselves in a positive or negative feedback loop [49,50]. Other descriptions are based on external adhesion cues or elastic heterogeneities. Cell-cell adhesions can participate in maintaining tumour compactness and radially symmetric geometry [51,52] but may drive phase separations between different sub-populations of cells [53]. Differential growth processes in heterogeneous elastic media or differential pulling forces can also trigger shape instabilities in two or three dimensions [54,55]. Most of the aforementioned instabilities however rely on cell migration, elastic or poroelastic tissue models, or numerical simulations. For a multicellular spheroid grown in a culture medium, however, there is no external medium on which the cells can pull, and only internal stresses can contribute to a potential shape instability. Also, experiments and modelling suggest that, on long timescales, cell-cell rearrangements lead to an effective viscous rheology of the tissue [37,41,[56][57][58][59], and an effective surface tension at tissue boundaries has been shown to play an important role in their shaping [38-40, 60, 61].
In the present work, we rely on a viscous description of cellular tissues on long timescales to establish analytically a new instability of multicellular spheroids that is viscocapillary in nature. This new instability does not require external heterogeneities or the presence of a substrate on which cells could exert pulling forces. Rather, it is powered by viscous shear stresses that build up within the spheroid due to permanent cell renewal, as illustrated in figure 1. The instability can be triggered by a change of internal properties such as cell-cell adhesion strength, cell-renewal rate, or metabolite supply.

Description of the model
We consider a multicellular spheroid embedded in a culture medium with a constant, physiological external pressure P ext , and which contains a given concentration ρ ext of a chemical substance necessary for cellular proliferation. This substance can be oxygen, growth factors, glucose or other nutrients, and will be referred to with the generic term of 'metabolites' in the following. The tissue within the spheroid is characterised in a continuum theory by a cell-number density and a cell-velocity field v. Analyses of the stress distribution in multicellular spheroids have shown that, at steady state under physiological osmotic conditions, the cell density is essentially homogeneous throughout one aggregate [62,63]. We therefore consider an incompressible tissue, for which the cell-number density is constant and the continuity equation reads Here, k p is the overall cell-production rate, considering cell division and cell death, and ∇ denotes the spatial derivative operator, contracted with the cell-velocity field to give its divergence. Neglecting inertia and in the absence of bulk external forces, force balance reads where σ denotes the total stress tensor. We further split the stress tensor into a dynamic part and a velocity-independent part. For an isotropic tissue, the latter reads −P1, where P is the tissue pressure and 1 the unity tensor. The dynamic part σ encodes the rheological properties of the tissue. The timescale of interest here is large compared to those of individual cellular processes, such as cell-cell rearrangements and cell renewal. We can therefore model the tissue as a viscous fluid with effective shear and bulk viscosities η and ζ, taking into account the long-term effects of cell production [37,56,57]. We obtain where ⊗ denotes the tensorial product and (∇ ⊗ v) T the transposed tensor of ∇ ⊗ v. The remarkable absence of any compression modulus is due to the nonconservation of cell number [57,64].
The system of equations is closed by specifying the expression of the cell-production rate k p . For simplicity, we assume that k p is independent of the stress σ and linearly dependent on the metabolite concentration ρ. This leads to where κ and k 0 are two positive phenomenological constants. Within the spheroid, metabolites diffuse with a coefficient D and are consumed or absorbed by the cells. Similarly, we assume a linear dependence of this absorption term in the metabolite concentration with a constant absorption rate α, as it has been observed in avascular tumour growth [43]: Here, ∂ t denotes the partial time derivative and Δ the Laplacian operator. Note that the convective term has been ignored in this equation, which is justified if nutrient diffusion is fast compared to its convective transport by the cell flow.
In the centre of the spheroid r = 0, the cell-velocity field vanishes and the metabolite concentration remains positive. At the outer surface, the metabolite concentration equals the external concentration ρ ext and the cell velocity equals that of the interface. Labelling R the spheroid stationary radius and δR its perturbation, these conditions read ρ r=R+δR = ρ ext and where e r is the unit radial vector and v is evaluated at the perturbed interface location. Finally, the outer surface is subjected solely to the isotropic, external pressure. The tangential component σ nt of the stress tensor therefore vanishes and its normal component is given by Laplace's law with surface tension γ: where H is the local curvature. Note that the first boundary condition mentioned here, in defining a specific location for the centre of the spheroid, breaks Galilean invariance. We comment in the following on the signification and consequences of this boundary condition.

Stationary equations
The stationary equations are characterised by spherical symmetry with radial coordinate r. The stationary continuity equation (1) and metabolite-diffusion equation (5) reduce to: and together with the corresponding boundary conditions (ρ stat ) r=R = ρ ext and (v stat r ) r=R = 0. The stationary force-balance condition and expressions of the stress-tensor components are given in appendix A.

Stationary solutions
To integrate this system of equations, we start by integrating the metabolite diffusion equation (8). With the characteristic metabolite-penetration length l D = D/α and the reduced variablesr = r/l D andR = R/l D , the metabolite concentration and the cell-production rate read and We can see from equation (10) that a stationary solution with spherical symmetry exists as long as κρ ext > k 0 to have a positive cell-division rate at the outer rim of the spheroid. Finally, the cell-velocity field reads v stat The boundary condition (v stat r ) r=R = 0 then leads to the following equation for the stationary radius: Note that the product κρ ext is linked to k 0 and the cell-production rate at the outer rim k ext p by k ext p = κρ ext − k 0 . The complete expressions of all the other quantities in this stationary state are given in appendix B.

Perturbed axisymmetric equations
We now investigate the linear stability of this stationary solution to axisymmetric perturbations. We choose a system of coordinates composed of radial coordinate r, polar angle θ, and azimuthal angle φ, and we study the perturbations with axial symmetry around θ = 0. The model equations with axial symmetry read for the continuity equation (1) and for the metabolite diffusion equation (5). Contrary to the stationary system of equations, the cell-velocity field cannot be solved independently of the force-balance condition (2). The other, coupled equations, are given in appendix C.

Linear decomposition
To integrate this system of equations to linear order in perturbations, we expand the angular dependence of the different perturbative fields onto the basis of axisymmetric, spherical harmonics. Following reference [65], the perturbations δR, δv r , δP, and δρ are expanded onto the basis of Legendre polynomials (P n (cos θ)) n∈N , and the perturbation δv θ is expanded onto the basis of the Gegenbauer polynomials (I n (cos θ)) n∈N , where I n = (P n−1 − P n+1 )/(2n + 1). The components v (n) r , P (n) , ρ (n) , and v (n) θ of respectively δv r , δP, δρ, and δv θ under this expansion are functions of the radial coordinate r, as the components R (n) of the expansion of the interface location δR are simple numbers. Explicitly, the δv r expansion for example reads as the δv θ expansion reads where the sum starts at n = 1 since the mode n = 0 is purely radial. Using these expansions, the metabolite diffusion equation (5) leads to where ρ (n) is the component for the mode number n of the nutrient field ρ and ω n the corresponding growth rate. The components of the perturbed velocity field are then determined by the perturbed continuity equation as well as by the force-balance equations, further given in appendix D.

Explicit solution in the limit of fast diffusion
To solve equation (17), we first consider the regime where the relaxation or growth of the perturbation modes as well as metabolite consumption are slow compared with metabolite diffusion over the characteristic lengths involved, at most equal to the spheroid radius. In that limit, the term [(α + ω n )/D]ρ (n) can be neglected in front of [n(n + 1)/r 2 )]ρ (n) . This approximation is certainly not valid for the mode n = 0, which needs to be computed separately. The calculation happens to be singular as well for n = 1. For n 2, the solution of equation (17) in this approximation is a simple power law, and we can further obtain the other perturbed quantities analytically. We have, for all n 2: where a n , b n , and c n are three integration constants, all proportional to the amplitude R (n) of the perturbed radius. The other quantities can be further expressed as linear combinations of these three integration constants. The other obtained expressions are given in appendix E. The integration constants are then determined using the boundary conditions. Plugging these solutions into the kinematic equation (6), we finally get the following mode growth rates: for all n 2. For the modes n = 0 and n = 1, we get separately ω 0 = k ext p − (1/9)k 0R 2 and These latter expressions are independent of the surface tension γ. This is because, for these two modes, perturbations in the curvature occur only to second order or higher.

Results
We can now discuss the instability in this fast-diffusion regime. The first term in equation (20) is destabilising and proportional to the cell-production rate at the outer surface. The second term is stabilising and results for one part from cell-death processes controlled by the parameters α/D (via l D ) and k 0 , and from the other part by surface tension. Increasing the viscosity lowers the contribution of surface tension, destabilising the spheroid. This underlines the mechanical origin of the instability, which relies on internal viscous stresses within the tissue, generated by differential cell flows (see figure 1). Therefore, this instability can only exist around a kinematic steady state with nonzero permanent cell flows, here from the outer surface towards the spheroid core.
Considering the stationary condition equation (12), we can verify that the modes n = 0 and n = 1 are always stable (see appendix E, equation (E.4)). For the other modes, in the absence of surface tension, the first term in equation (20) is dominant at large n, meaning that, without this contribution, the spheroid is always unstable with a rate asymptotically equal to that of cell division at its outer rim. Surface tension however stabilises the spheroid, since it contributes by a term scaling as γq n /(2η) at large n, where q n = n/R. Depending on the values of the different parameters, we therefore expect a potential instability to develop at a finite value of n, corresponding to a finite wavelength λ n ∼ 2πR/n.
It is interesting to investigate the behaviour of the most unstable mode as a function of the stationary radius R. Since the equation characterising this mode is in general fourth order in n, we investigate separately the limits of small and large spheroids. In the limit of small radii, the least stable mode is ω 0 , since curvature is large and surface tension strongly stabilises all modes for larger values of n. In the limit of large radii, the most unstable mode occurs at n = n max with the following asymptotic expression, linear in R: n max ηκρ ext /(γl D )R. This scaling indicates that the associated wavelength converges towards a finite value at large radius R. In this limit, we asymptotically reach the case of a flat surface, with an instability that evokes what has been proposed for epithelial tissues [66,67]. The corresponding growth rate reads ω max κρ ext − γκρ ext /(ηl D ).
In the generic case where metabolite diffusion is not necessarily fast compared to metabolite consumption or perturbation growth, the solution for ρ (n) in equation (17) is a function of the associated growth rate ω n . Equation (6) then becomes an implicit equation for the growth rate ω n , which cannot be solved analytically. We report the implicit equation corresponding to the mode n = 0 as an example in appendix F, equation (F.1).
To compute the growth rates ω n numerically, we now estimate the different parameter values. The shear viscosity of cellular aggregates has been estimated in different experiments, leading to η 10 4 − 10 5 Pa s [37,39,56]. Tissue surface tensions have been measured for different tissue types. Measurements for γ give values ranging from a fraction up to several millinewton per metre [37,39,60,61]. We further assume a typical cell-division rate at the outer surface of the spheroid k ext p and a cell-death rate in the absence of metabolites k 0 of one per day. Cellular growth within the spheroid can be limited by different types of metabolites. Depending typically on the molecular size of a given metabolite, its diffusion coefficient can take different values. Estimates of the diffusion coefficient of growth factors and glucose in avascular tumours range from 10 −6 cm 2 h −1 for growth factors [43]   even larger values for oxygen [43,70]. In the following, we shall investigate the influence of a variation of this particular parameter. Finally, to obtain radii of a few 100 micrometres, we choose to have comparable values for the characteristic penetration length of metabolites l D . This leads to values for the metabolite-consumption rate α of several tenths per day for D ∼ 10 −6 cm 2 h −1 , scaled accordingly when D is varied to keep l D constant.
We illustrate in figure 2 the obtained results for the stationary solution of the model. The stationary-state profiles show that cells divide preferentially close to the outer surface and disappear in the centre (panel (a)), due to the lack of metabolites penetrating the tissue (panel (c)). As a result, cells flow inwards from the outer surface to the centre (panels (b) and (d)). This result is in agreement with experimental measurements of cellular flows in multicellular spheroids using fluorescently labelled particles [34].
We illustrate in figure 3 the central result of our study, that is the obtained mode structure of the instability, for the four different steady states presented in figures 2(a) and (b), using the same colour code. The system is unstable as soon as at least one perturbation mode displays a growth rate ω with a positive real part. When this is the case, we expect the fastest growing modes to develop first, corresponding to the maximum of each series of points presented in these plots. A visual display of the shapes associated with the deformation modes n = 0 to 5 is shown in appendix G, figure G1. Figure 3(a) shows the mode structure in the approximation of the analytic solution of section 4.3. There is a transition at finite wavelength as a function of the metabolite-consumption rate α from a stable to an unstable regime, with a range of unstable modes. These appear for example when nutrient consumption decreases at fixed external cell-division rate, surface tension and viscosity. With the parameters chosen here, the first mode to become unstable is n = 5 ( figure 3(a), orange squares). This is associated with a stationary radius R 5 622 μm and corresponds to an unstable wavelength λ 5 782 μm.
In figures 3(b) and (c), we display respectively the real and imaginary parts of the mode growth rates obtained in the full model solved numerically, with the parameter sets of figure 2 and a diffusion coefficient D = 5 × 10 −10 cm 2 s −1 , corresponding approximately to that estimated for growth factors [43]. The system also displays first an instability at finite wavelength (for the mode n = 3, green diamonds), but the modes n = 0 and n = 1 can now be unstable. In addition, these instabilities are oscillatory, as characterised by non-zero imaginary parts of ω, corresponding to the characteristic frequencies of the unstable oscillations (see figure 3(c)). This result stems from the existence of a delay in the response of the cells located in the bulk of the aggregate to a given perturbation of the outer surface, due to the finite kinetics of metabolite penetration into the spheroid.
Our model ignores the convective cellular flow in the metabolite diffusion equation (5). For this approximation to be valid, one requires that D/l D be larger than the amplitude v of the cellular flow within the aggregate. Figure 3(d) displays the mode growth rates obtained with a diffusion coefficient D a thousand times that of figures 3(b) and (c), corresponding approximately to the diffusion coefficient of small nutrient molecules such as glucose [43,68,69]. With this value of the diffusion coefficient, we have 1 when v is estimated from the curves shown in figure 2(b) in the least favorable case, satisfying largely the required condition for neglecting convective flows. We show in figure 3(d) that, perturbing around the same steady state by rescaling the metabolite-absorbing rate α, the instability occurs at a similar radius and finite wavelength as those reported above. In addition, the growth rates of the high-order modes (n 5-6) are largely unchanged. Interestingly, the oscillatory instability is lost for n = 0 and n = 1. This is the signature of the fact that here metabolite diffusion is sufficiently fast to allow for an almost instantaneous response of the inner cells to perturbations of the outer surface. We however do not recover the analytic results of figure 3(a) for small values of n. This stems from the fact that, in the results of figure 3(d), we scale the metabolite consumption rate with the diffusion coefficient to keep the same steady-state sizes as in figure 3(b), as the analytic limit was obtained for small values of α/D with respect to 1/R 2 . We show in appendix H, figure H1, the mode structure obtained with intermediate values of the diffusion coefficient D between those of figures 3(b)-(d), following the same rescaling procedure of the metabolite-consumption rate α.
As mentioned at the end of section 2, our boundary conditions specify that the cell-velocity field vanishes at the centre of the spheroid located at r = 0, which breaks Galilean invariance. Doing so, the mode n = 1 here corresponds to an actual deformation of the flow pattern within the spheroid, with an outer boundary displaced with respect to the point where the flow pattern converges. As a consequence, this mode is not necessarily marginal, contrary to many standard spherical-harmonic perturbation analyses. A similar interpretation of the mode n = 1 is found, e.g., in the context of the deformations of the actin cortex of a spherical cell [71], where it corresponds to a relative translation of the inner part of the cell cortex with respect to the outer boundary of the cell, rather than to a global translation of the whole system. Therefore, even if this mode corresponds to a global translation of the inner boundary of the cortex with no deformation, the relative positions of the inner and outer boundaries of the cortex do vary, which corresponds to a spatial variation of its overall thickness. As a result, the mode n = 1 is not marginal. Similarly, here, we investigate relative displacements of the different cell layers with each other. Each cell layer is purely translated in the mode n = 1, but the relative positions of the different layers varies within the spheroid. Therefore, even though the overall external shape of the spheroid is unchanged, the cell-flow pattern is modified by this perturbation.
At threshold, the unstable modes are expected to grow exponentially in time, albeit potentially on long timescales. If this scenario is expected at sufficiently small amplitudes, where the linear regime of small deformations around the spherical shape is valid, large-amplitude deformations are expected to follow another, more complicated dynamics. In particular, we expect that when perturbation amplitudes reach a finite fraction of the spheroid radius, geometric nonlinearities will induce more complex flow patterns because of an asymmetry between inward and outward deformations. We expect eventually that outward protrusions would at late stages outgrow inward protrusions, the latter being limited by the original spheroid size. We can speculate that such an unbalance would lead to a global growth of the overall spheroid mass.

Discussion
In this work, we have shown the potential existence of an instability in spherical tissues, which can develop from steady states that are limited by the supply of metabolites diffusing from their microenvironment. The present instability relies neither on cell motility nor on the presence of external forces, but rather stems from the presence of viscous shear stresses generated by the spatial organisation of cell renewal within the tissue. We have shown that the instability develops at a finite wavenumber, which reflects a balance of viscous shear stresses with those stemming from surface tension. We propose that this mechanism could be observed in multicellular spheroids in culture, which would be an ideal system for testing the influence of different parameters controlling the instability, such as tissue viscosity, surface tension, and metabolite supply. The former two could be changed, e.g., by varying the expression of proteins implicated in cell-surface adhesion or actin-cortex contractility [38,40,60].
The proposed instability here evokes other already proposed instabilities in the context of the cell cytoskeleton, driven by actin-polymerisation dynamics [72,73]. In these studies of the stability of cell fragments on a substrate, an inward flow is driven by actin polymerisation at the outer edge and actin depolymerisation in the bulk. As a result, an originally circular cell fragment can become unstable and spontaneously acquire a polarisation. Important differences between our current study and these previous works however exist. In the stability analysis of circular cell fragments, the generation of new material occurs only at the outer surface, where actin polymerises. In our current model, cell production is a global, bulk effect, which varies continuously within the spheroid. The second difference is the rheology, which here corresponds a Stoke flow with no contact with an external substrate, as in the case of cell fragments there is a Darcy flow, rendering the two types of instability different. Associated to that, the third difference is that our current study is three dimensional as these previous works are two dimensional. As a result, our flow pattern has no anchoring to the external world and requires a minimal thickness over which cells divide to become unstable.
In addition to be applicable to multicellular aggregates, one can wonder if similar instabilities could arise in vivo. During development, transition from solid-like to fluid-like tissue properties, tissue surface tension, and flow patterns have been shown to play a crucial role (see, e.g., [74]). Recently, three-dimensional aggregates of mouse embryonic stem cells have been shown to undergo a first morphological transformation from a spherical into an oblong shape during gastrulation, associated with a reduced level of E-cadherin expression at the developing tip [75]. Such shapes resemble a superposition of instability modes such as n = 2 and 3, and potentially higher, as illustrated in figure G1. Interestingly with respect to our current study, the polarisation of E-cadherin expression precedes the onset of tip formation, and when the level of E-cadherin expression is maintained high, the aggregate remains generally devoid of any pole [75]. These observations suggest a role for a reduced surface tension in the development of the protrusion, similar to what we are proposing here. However, in these examples, and to our knowledge more generally during development, it seems that an original inhomogeneity in the tissue rheological parameters-such as surface tension and viscosity-is at the origin of the shape formation. Such inhomogeneities are however absent in our proposed mechanism, which relies solely on the presence of a permanent flow of duplicating cells.
A domain to which we can speculate that the present mechanism applies is the evolution of microtumours after a long period of dormancy. Small primary tumours or early metastases often enter a state where their sizes remain steady, before they resume growth or disappear [76]. Such a dormant state can last for a long time and is at the origin of late cancer reappearance, sometimes years after the original treatment [77,78]. It is recognised that microscopic, clinically occult tumours are very common in the population and that only a tiny fraction of them ever becomes clinically relevant [79][80][81]. Understanding the factors that can destabilise a dormant tumour is therefore of crucial importance. The main mechanisms at the origin of such steady states are angiogenic dormancy, cellular dormancy (G0-G1 arrest) and immunosurveillance [76]. In angiogenic dormancy, the tumour is limited in its growth by the lack of metabolites, which are brought by blood vessels that do not penetrate the tumour [76,80,81]. This limitation keeps the microtumour to sizes typically smaller than 1-2 mm in diameter, until the angiogenic switch is triggered [81][82][83].
Tumour-growth models have explored the effects of a wide variety of biological processes [44,[84][85][86][87][88]. In support of the current surface mechanism, it has recently been shown that colon cancer xenografts grow primarily from their surfaces [89,90], which corresponds to the steady-state patterns of cell duplications on which our current study relies. In addition, clonal expansion largely depends on the location of a clone within the tumour [91], suggesting that differences in geometrical or physical properties within the tumour are major contributors to heterogeneous clonal expansion. This latter observation leads us to speculate that, after a first instability such as the one proposed here, the resulting irregular shape creates different microenvironments for different parts of the tumour, further driving different epigenetic and maybe even later genetic transformations by diverse selection processes.
The proposed instability can be triggered by a change of internal properties such as cell-cell adhesion strength, cell-renewal rate, metabolite supply, or other parameters affecting the overall spheroid size. In a microtumour, such changes might be multifactorial, e.g. ageing, a change in the person's metabolism, in the immune system's activity, or in drug delivery or efficiency. While our model does not address the long-time evolution of these parameters, a multicellular aggregate or a microtumour can still be unstable by the mechanism proposed here. Our study could therefore be of importance for determining which parameters control the spherical stability. Using multicellular spheroids as model systems, it could on the long run participate in guiding which aspects of tumour development should be targeted by medication. 2η The pressure is given by P stat = P stat r=0 + κρ ext ζ + To express the boundary conditions with axisymmetry, we need to express the local curvature H at the outer surface of the spheroid. To first order in perturbations in δR, it is given by The different quantities evaluated at r = R + δR read, to first order in perturbations: (σ nn ) r=R+δR = σ stat rr + In the right-hand sides, all the quantities that depend on r are evaluated at r = R. Note that σ stat rθ is formally written here but is actually zero in the spherical symmetric case of our stationary state. Similarly, the stationary velocity v stat r is also zero at the outer stationary boundary r = R. Figure G1. Schematic illustration of the lowest-order, axisymmetric, perturbative modes. The modes n = 0 to n = 5 are represented with an arbitrary amplitude. In each schematic representation, the unperturbed spherical shape is represented with a dashed line. In the first, n = 0 mode, the axisymmetry is indicated as a rotational invariance in φ, and all modes depicted here are invariant under this transformation.