# Viscocapillary Instability in Cellular Spheroids

Matthieu Martin1 and Thomas Risler1 1Laboratoire Physico-Chimie Curie,
Institut Curie, PSL Research University, Sorbonne Université, CNRS, 26 rue
d’Ulm, 75005 Paris, France<EMAIL_ADDRESS>
###### Abstract

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.

††: New J. Phys.

##  1 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 [16, 15]
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
[22, 21], 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 [26, 18, 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 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 toward the centre of the spheroid [42, 34].
Under some circumstances, the spheroid reaches a steady-state size [43, 44,
34].

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 H.P. 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, 56, 57, 41, 58, 59], and an effective surface
tension at tissue boundaries has been shown to play an important role in their
shaping [38, 60, 61, 39, 40].

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.

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 toward 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.

The instability can be triggered by a change of internal properties such as
cell-cell adhesion strength, cell-renewal rate, or metabolite supply.

##  2 Description of the model

We consider a multicellular spheroid embedded in a culture medium with a
constant, physiological external pressure $P^{\rm ext}$, and which contains a
given concentration $\rho^{\rm 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 $\mathbf{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

$\nabla\cdot\mathbf{v}=k_{\rm p}\,.$ (1)

Here, $k_{\rm p}$ is the overall cell-production rate, considering cell
division and cell death, and $\nabla$ 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

$\nabla\cdot\boldsymbol{\sigma}=0\,,$ (2)

where $\boldsymbol{\sigma}$ 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 $-P\,\mathbf{1}$, where $P$ is the tissue
pressure and $\mathbf{1}$ the unity tensor. The dynamic part
$\boldsymbol{\sigma}^{\prime}$ 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 $\eta$ and $\zeta$, taking into account the long-
term effects of cell production [37, 56, 57]. We obtain

$\boldsymbol{\sigma}^{\prime}=\eta\left[\nabla\otimes\mathbf{v}+(\nabla\otimes\mathbf{v})^{\rm
T}-\frac{2}{3}(\nabla\cdot\mathbf{v})\,\mathbf{1}\right]+\zeta(\nabla\cdot\mathbf{v})\,\mathbf{1}\,,$
(3)

where $\otimes$ denotes the tensorial product and
$(\nabla\otimes\mathbf{v})^{\rm T}$ the transposed tensor of
$\nabla\otimes\mathbf{v}$. The remarkable absence of any compression modulus
is due to the nonconservation of cell number [64, 57].

The system of equations is closed by specifying the expression of the cell-
production rate $k_{\rm p}$. For simplicity, we assume that $k_{\rm p}$ is
independent of the stress $\boldsymbol{\sigma}$ and linearly dependent on the
metabolite concentration $\rho$. This leads to

$k_{\rm p}=\kappa\rho-k_{0}\,,$ (4)

where $\kappa$ 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 $\alpha$, as it has been observed in avascular tumour growth [43]:

$\partial_{t}\rho=D\Delta\rho-\alpha\rho\,.$ (5)

Here, $\partial_{t}$ denotes the partial time derivative and $\Delta$ 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 center of the spheroid $\mathbf{r}=\mathbf{0}$, the cell-velocity field
vanishes and the metabolite concentration remains positive. At the outer
surface, the metabolite concentration equals the external concentration
$\rho^{\rm ext}$ and the cell velocity equals that of the interface. Labelling
$R$ the spheroid stationary radius and $\delta R$ its perturbation, these
conditions read $\rho_{r=R+\delta R}=\rho^{\rm ext}$ and

$\partial_{t}(\delta R)=\mathbf{v}_{r=R+\delta R}\cdot\mathbf{e}_{r}\,,$ (6)

where $\mathbf{e}_{r}$ is the unit radial vector and $\mathbf{v}$ is evaluated
at the perturbed interface location. Finally, the outer surface is subjected
solely to the isotropic, external pressure. The tangential component
$\sigma_{\rm nt}$ of the stress tensor therefore vanishes and its normal
component is given by Laplace’s law with surface tension $\gamma$:
$(\sigma_{\rm nn})_{r=R+\delta R}=-P^{\rm ext}-\gamma H$, where $H$ is the
local curvature. Note that the first boundary condition mentioned here, in
defining a specific location for the center of the spheroid, breaks Galilean
invariance. We comment in the following on the signification and consequences
of this boundary condition.

##  3 Stationary solution

###  3.1 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:

$\frac{1}{r^{2}}\,\partial_{r}(r^{2}\,v_{r})=k_{\rm p}$ (7)

and

$\frac{D}{r^{2}}\,\partial_{r}(r^{2}\partial_{r}\rho)=\alpha\rho\,,$ (8)

together with the corresponding boundary conditions $(\rho^{\rm
stat})_{r=R}=\rho^{\rm ext}$ and $(v_{r}^{\rm stat})_{r=R}=0$. The stationary
force-balance condition and expressions of the stress-tensor components are
given in A.

###  3.2 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}=\sqrt{D/\alpha}$ and the reduced variables $\bar{r}=r/l_{D}$ and
$\bar{R}=R/l_{D}$, the metabolite concentration and the cell-production rate
read

$\rho^{\rm stat}=\rho^{\rm
ext}\frac{\bar{R}}{\sinh{\bar{R}}}\,\frac{\sinh{\bar{r}}}{\bar{r}}$ (9)

and

$k^{\rm stat}_{\rm p}=\kappa\rho^{\rm
ext}\,\frac{\bar{R}}{\sinh{\bar{R}}}\,\frac{\sinh{\bar{r}}}{\bar{r}}-k_{0}\,.$
(10)

We can see from Eq. 10 that a stationary solution with spherical symmetry
exists as long as $\kappa\rho^{\rm ext}>k_{0}$ to have a positive cell-
division rate at the outer rim of the spheroid. Finally, the cell-velocity
field reads

$v^{\rm stat}_{r}=\kappa\rho^{\rm
ext}\,\frac{R}{\sinh{\bar{R}}}\left[\frac{\cosh{\bar{r}}}{\bar{r}}-\frac{\sinh{\bar{r}}}{\bar{r}^{2}}\right]-\frac{1}{3}k_{0}r\,.$
(11)

The boundary condition $(v_{r}^{\rm stat})_{r=R}=0$ then leads to the
following equation for the stationary radius:

$\frac{1}{\bar{R}}\coth{\bar{R}}-\frac{1}{\bar{R}^{2}}=\frac{k_{0}}{3\kappa\rho^{\rm
ext}}\,.$ (12)

Note that the product $\kappa\rho^{\rm ext}$ is linked to $k_{0}$ and the
cell-production rate at the outer rim $k_{\rm p}^{\rm ext}$ by $k_{\rm p}^{\rm
ext}=\kappa\rho^{\rm ext}-k_{0}$. The complete expressions of all the other
quantities in this stationary state are given in B.

##  4 Mode computation

###  4.1 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 $\theta$, and azimuthal angle $\phi$, and
we study the perturbations with axial symmetry around $\theta=0$. The model
equations with axial symmetry read

$\frac{1}{r^{2}}\,\partial_{r}(r^{2}\,v_{r})+\frac{1}{r\sin{\theta}}\,\partial_{\theta}(\sin{\theta}\,v_{\theta})=k_{\rm
p}$ (13)

for the continuity equation 1 and

$\partial_{t}\rho=\frac{D}{r^{2}}\,\partial_{r}(r^{2}\partial_{r}\rho)+\frac{D}{r^{2}\sin{\theta}}\,\partial_{\theta}(\sin{\theta}\partial_{\theta}\rho)-\alpha\rho$
(14)

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 C.

###  4.2 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 ref. [65], the
perturbations $\delta R$, $\delta v_{r}$, $\delta P$, and $\delta\rho$ are
expanded onto the basis of Legendre polynomials
$\left(P_{n}(\cos{\theta})\right)_{n\in\mathbb{N}}$, and the perturbation
$\delta v_{\theta}$ is expanded onto the basis of the Gegenbauer polynomials
$\left(I_{n}(\cos{\theta})\right)_{n\in\mathbb{N}}$, where
$I_{n}=(P_{n-1}-P_{n+1})/(2n+1)$. The components $v_{r}^{(n)}$, $P^{(n)}$,
$\rho^{(n)}$, and $v_{\theta}^{(n)}$ of respectively $\delta v_{r}$, $\delta
P$, $\delta\rho$, and $\delta v_{\theta}$ under this expansion are functions
of the radial coordinate $r$, as the components $R^{(n)}$ of the expansion of
the interface location $\delta R$ are simple numbers. Explicitly, the $\delta
v_{r}$ expansion for example reads

$\delta
v_{r}(r,\theta)=\sum_{n=0}^{\infty}v_{r}^{(n)}\,P_{n}(\cos{\theta})\,,$ (15)

as the $\delta v_{\theta}$ expansion reads

$\delta
v_{\theta}(r,\theta)=\sum_{n=1}^{\infty}v_{\theta}^{(n)}\,\frac{I_{n+1}(\cos{\theta})}{\sin{\theta}}\,,$
(16)

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

$\partial_{r}^{2}\rho^{(n)}+\frac{2}{r}\,\partial_{r}\rho^{(n)}-\left[\frac{\alpha+\omega_{n}}{D}+\frac{n(n+1)}{r^{2}}\right]\rho^{(n)}=0\,,$
(17)

where $\rho^{(n)}$ is the component for the mode number $n$ of the nutrient
field $\rho$ and $\omega_{n}$ the corresponding growth rate. The components of
the perturbed velocity field are then determined by the perturbed continuity
equation

$\partial_{r}\,v_{r}^{(n)}+\frac{2}{r}\,v_{r}^{(n)}+\frac{1}{r}\,v_{\theta}^{(n)}=\kappa\,\rho^{(n)}$
(18)

as well as by the force-balance equations, further given in D.

###  4.3 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
$[(\alpha+\omega_{n})/D]\rho^{(n)}$ can be neglected in front of
$[n(n+1)/r^{2})]\rho^{(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\geq 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\geq 2$:

$\displaystyle\rho^{(n)}$ $\displaystyle=c_{n}\bar{r}^{n}$ $\displaystyle
v_{r}^{(n)}$ $\displaystyle=a_{n}\bar{r}^{n+1}+b_{n}\bar{r}^{n-1}$
$\displaystyle v_{\theta}^{(n)}$
$\displaystyle=\left[-(n+3)a_{n}+\kappa\,R\,c_{n}\right]\bar{r}^{n+1}-(n+1)b_{n}\bar{r}^{n-1}\,,$
(19)

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 E.

The integration contants are then determined using the boundary conditions.
Plugging these solutions into the kinematic equation 6, we finally get the
following mode growth rates:

$\omega_{n}=\frac{2n^{2}+5n+3}{2n^{2}+4n+3}\,k_{\rm p}^{\rm
ext}-\frac{1}{2n^{2}+4n+3}\left[\frac{n+1}{3}\,k_{0}\left(\frac{R}{l_{D}}\right)^{2}+(2n^{2}+5n+2)\frac{\gamma}{2\eta}\,\frac{n}{R}\right]$
(20)

for all $n\geq 2$. For the modes $n=0$ and $n=1$, we get separately
$\omega_{0}=k_{\rm p}^{\rm ext}-(1/9)k_{0}\bar{R}^{2}$ and
$\omega_{1}=(1/3)k_{\rm p}^{\rm ext}-(1/18)k_{0}\bar{R}^{2}$. These latter
expressions are independent of the surface tension $\gamma$. This is because,
for these two modes, perturbations in the curvature occur only to second order
or higher.

##  5 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 $\alpha/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 toward the spheroid core.

Considering the stationary condition equation 12, we can verify that the modes
$n=0$ and $n=1$ are always stable (see E, equation E). 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 $\gamma q_{n}/(2\eta)$ 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 $\lambda_{n}\sim 2\pi 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 $\omega_{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_{\rm max}$ with the following asymptotic
expression, linear in $R$: $n_{\rm max}\simeq\sqrt{\eta\kappa\rho^{\rm
ext}/(\gamma l_{D})}\,R$. This scaling indicates that the associated
wavelength converges toward 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 $\omega_{\rm max}\simeq\kappa\rho^{\rm
ext}-\sqrt{\gamma\kappa\rho^{\rm ext}/(\eta l_{D})}$.

In the generic case where metabolite diffusion is not necessarily fast
compared to metabolite consumption or perturbation growth, the solution for
$\rho^{(n)}$ in equation 17 is a function of the associated growth rate
$\omega_{n}$. Equation 6 then becomes an implicit equation for the growth rate
$\omega_{n}$, which cannot be solved analytically. We report the implicit
equation corresponding to the mode $n=0$ as an example in F, equation 39.

To compute the growth rates $\omega_{n}$ numerically, we now estimate the
different parameter values. The shear viscosity of cellular aggregates has
been estimated in different experiments, leading to $\eta\simeq 10^{4}-10^{5}$
Pa$\cdot$s [37, 56, 39]. Tissue surface tensions have been measured for
different tissue types. Measurements for $\gamma$ give values ranging from a
fraction up to several millinewton per meter [37, 60, 61, 39]. We further
assume a typical cell-division rate at the outer surface of the spheroid
$k_{\rm p}^{\rm ext}$ 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}\cdot$h-1 for growth factors [43]
to 1.5 10-3 cm${}^{2}\cdot$h-1 for glucose [68, 69, 43], or even larger values
for oxygen [70, 43]. In the following, we shall investigate the influence of a
variation of this particular parameter. Finally, to obtain radii of a few
hundred 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 $\alpha$ of several tenths per day for
$D\sim$10-6 cm${}^{2}\cdot$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.

Figure 2:  Stationary state of a multicellular spheroid. (a,b) The cell-
production rate $k_{\rm p}$ (a) and the radial cell-velocity field $v_{r}$ (b)
are displayed as functions of the radial distance $r$ for four different
values of the stationary radius $R$. Parameters common to the four sets of
curves are: $\eta=80$ kPa$\cdot$s, $\gamma=100$ $\mu$N$\cdot$m-1, and $k_{\rm
p}^{\rm ext}=k_{0}=1$ d-1 (d-1 stems for “per day”). The stationary state
depends on the diffusion coefficient of metabolites $D$ and metabolite-
absorption rate $\alpha$ through their ratio $D/\alpha=l_{D}^{2}$ only. We
display the resulting stationary state for four different values of this
ratio. The corresponding stationary radii and metabolite penetration lengths
are $R\simeq$ 803 $\mu$m and $l_{D}\simeq$ 170 $\mu$m (blue curves), $R\simeq$
622 $\mu$m and $l_{D}\simeq$ 131 $\mu$m (orange curves), $R\simeq$ 526 $\mu$m
and $l_{D}\simeq$ 111 $\mu$m (green curves), and $R\simeq$ 464 $\mu$m and
$l_{D}\simeq$ 98 $\mu$m (red curves). (c,d) To give a better feeling of the
spherical symmetry, we display in panel (c) the stationary metabolite-
concentration profile with respect to its outer-surface value $\rho/\rho^{\rm
ext}$ and in panel (d) the corresponding radial cell-velocity field $v_{r}$
(in $\mu$m$\cdot$d-1), for the parameter set corresponding to the blue curves
of panels (a) and (b). Plots are made in a plane of symmetry of the spheroid.

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 labeled 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 figure 2a and b, using the same colour code.

Figure 3:  Viscocapillary mode growth rates of a multicellular spheroid as
functions of the mode number $n$, obtained (a) in the analytic limit of
section 4.3, (b,c,d) numerically using the full model in both limits of slow
(b,c) and fast (d) diffusion. (a) Parameters are the same as those of figure 2
using the same colour code. The mode-growth rates depend on the diffusion
coefficient of metabolites $D$ and metabolite-absorption rate $\alpha$ through
their ratio $D/\alpha=l_{D}^{2}$ only. We have $l_{D}\simeq$ 170 $\mu$m (blue
circles), $l_{D}\simeq$ 131 $\mu$m (orange squares), $l_{D}\simeq$ 111 $\mu$m
(green diamonds), and $l_{D}\simeq$ 98 $\mu$m (red triangles). (b,c) Real (b)
and imaginary (c) parts of the mode-growth rates computed numerically with a
diffusion coefficient of metabolites equal to that of growth factors:
$D=5\cdot 10^{-10}$ cm${}^{2}\cdot$s-1 [43]. The metabolite-absorption rate is
varied accordingly to keep the same steady state as in figure 2: $\alpha=0.15$
d-1 (blue circles), $\alpha=0.25$ d-1 (orange squares), $\alpha=0.35$ d-1
(green diamonds), and $\alpha=0.45$ d-1 (red triangles). (d) Mode-growth rates
computed numerically with a diffusion coefficient of metabolite equal to that
of glucose: $D=5\cdot 10^{-7}$ cm${}^{2}\cdot$s-1 [68, 69, 43]. The
metabolite-absorption rate is scaled accordingly by a factor of a thousand to
perturb around the same steady state as previously: $\alpha=150$ d-1 (blue
circles), $\alpha=250$ d-1 (orange squares), $\alpha=350$ d-1 (green
diamonds), $\alpha=450$ d-1 (red triangles). The mode-growth rates in (a) and
(d) are real.

The system is unstable as soon as at least one perturbation mode displays a
growth rate $\omega$ 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 G,
figure 4. Figure 3a 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 $\alpha$ 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 3a, orange squares). This is associated with
a stationary radius $R_{5}\simeq$ 622 $\mu$m and corresponds to an unstable
wavelength $\lambda_{5}\simeq$ 782 $\mu$m.

In figure 3b 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\cdot 10^{-10}$
cm${}^{2}\cdot$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 $\omega$, corresponding to the
characteristic frequencies of the unstable oscillations (see figure 3c). 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 3d displays the mode growth rates obtained with a diffusion coefficient
$D$ a thousand times that of figure 3b,c, corresponding approximately to the
diffusion coefficient of small nutrient molecules such as glucose [68, 69,
43]. With this value of the diffusion coefficient, we have $(\ell_{D}\cdot
v)/D\sim 4.10^{-3}\ll 1$ when $v$ is estimated from the curves shown in figure
2b in the least favorable case, satisfying largely the required condition for
neglecting convective flows. We show in figure 3d that, perturbing around the
same steady state by rescaling the metabolite-absorbing rate $\alpha$, 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\geq 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 3a for small values of $n$. This stems
from the fact that, in the results of figure 3d, we scale the metabolite
consumption rate with the diffusion coefficient to keep the same steady-state
sizes as in figure 3b, as the analytic limit was obtained for small values of
$\alpha/D$ with respect to $1/R^{2}$. We show in H, figure 5, the mode
structure obtained with intermediate values of the diffusion coefficient $D$
between those of figures 3b,c and 3d, following the same rescaling procedure
of the metabolite-consumption rate $\alpha$.

As mentioned at the end of section 2, our boundary conditions specify that the
cell-velocity field vanishes at the center of the spheroid located at
$\mathbf{r}=\mathbf{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.

##  6 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, 60, 40].

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 4. 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 [80, 81, 76]. This limitation keeps the microtumour
to sizes typically smaller than 1–2 mm in diameter, until the angiogenic
switch is triggered [82, 83, 81].

Tumour-growth models have explored the effects of a wide variety of biological
processes [84, 85, 44, 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.

We thank F. Brochard-Wyart, D. Gonzalez-Rodriguez, K. Guevorkian, J.-F.
Joanny, and J. Prost for insightful discussions and useful comments on the
manuscript. This work received support from the LabEx Cell(n)Scale (former
CelTisPhyBio), grants ANR-11-LABX-0038 and ANR-10-IDEX-0001-02.

##  Appendix A Additional stationary mechanical equations

The stationary, non-trivial component of the force-balance equation 2 reads

$\partial_{r}\sigma_{rr}+\frac{2}{r}(\sigma_{rr}-\sigma_{\theta\theta})=0\,,$
(21)

and the non-zero components of the dynamic part of the stress tensor
$\sigma^{\prime}$ as given by equation 3 read

$\displaystyle\sigma^{\prime}_{rr}$
$\displaystyle=2\eta\,\partial_{r}v_{r}+\left(\zeta-\frac{2}{3}\eta\right)\nabla\cdot\mathbf{v}$
$\displaystyle\sigma^{\prime}_{\theta\theta}$
$\displaystyle=\sigma^{\prime}_{\phi\phi}=2\eta\,\frac{v_{r}}{r}+\left(\zeta-\frac{2}{3}\eta\right)\nabla\cdot\mathbf{v}\,,$
(22)

where

$\nabla\cdot\mathbf{v}=\partial_{r}v_{r}+\frac{2}{r}v_{r}\,.$ (23)

The corresponding boundary condition reads

$(\sigma_{rr}^{\rm stat})_{r=R}=-P^{\rm ext}-\frac{2\gamma}{R}\,.$ (24)

##  Appendix B Additional stationary expressions

The non-trivial components of the dynamic part of the stress tensor read

$\sigma^{\prime\,\rm stat}_{rr}=\kappa\rho^{\rm
ext}\,\frac{\bar{R}}{\sinh{\bar{R}}}\left[-4\eta\frac{\cosh{\bar{r}}}{\bar{r}^{2}}+\left(\zeta+\frac{4}{3}\eta+\frac{4\eta}{\bar{r}^{2}}\right)\frac{\sinh{\bar{r}}}{\bar{r}}\right]-\zeta
k_{0}$ (25)

and

$\sigma^{\prime\,\rm stat}_{\theta\theta}=\kappa\rho^{\rm
ext}\,\frac{\bar{R}}{\sinh{\bar{R}}}\left[2\eta\frac{\cosh{\bar{r}}}{\bar{r}^{2}}+\left(\zeta-\frac{2}{3}\eta-\frac{2\eta}{\bar{r}^{2}}\right)\frac{\sinh{\bar{r}}}{\bar{r}}\right]-\zeta
k_{0}\,.$ (26)

The pressure is given by

$P^{\,\rm stat}=P^{\,\rm stat}_{\bar{r}=0}+\kappa\rho^{\rm
ext}\,\left(\zeta+\frac{4}{3}\eta\right)\,\frac{\bar{R}}{\sinh(\bar{R})}\left(\frac{\sinh(\bar{r})}{\bar{r}}-1\right)\,,$
(27)

where

$P^{\,\rm stat}_{\bar{r}=0}=P^{\rm ext}+\frac{2\gamma}{R}-\zeta
k_{0}+\frac{\kappa\rho^{\rm
ext}}{\bar{R}^{2}}\left[4\eta\left(1-\bar{R}\coth{\bar{R}}\right)+\left(\zeta+\frac{4}{3}\eta\right)\frac{{\bar{R}}^{3}}{\sinh{\bar{R}}}\right]\,.$
(28)

##  Appendix C Additional perturbed axisymmetric equations

With axisymmetry, the non-trivial components of the force-balance equation 2
read

$\displaystyle\partial_{r}\sigma_{rr}+\frac{1}{r}\partial_{\theta}\sigma_{r\theta}+\frac{1}{r}\left[2\sigma_{rr}-\sigma_{\theta\theta}-\sigma_{\phi\phi}+\sigma_{r\theta}\cot{\theta}\right]$
$\displaystyle=0$
$\displaystyle\partial_{r}\sigma_{r\theta}+\frac{1}{r}\partial_{\theta}\sigma_{\theta\theta}+\frac{1}{r}\left[3\sigma_{r\theta}+(\sigma_{\theta\theta}-\sigma_{\phi\phi})\cot{\theta}\right]$
$\displaystyle=0\,,$ (29)

and the non-zero components of the dynamic part of the stress tensor, as given
by equation 3, are

$\displaystyle\sigma^{\prime}_{rr}$
$\displaystyle=2\eta\,\partial_{r}v_{r}+\left(\zeta-\frac{2}{3}\eta\right)\nabla\cdot\mathbf{v}$
$\displaystyle\sigma^{\prime}_{\theta\theta}$
$\displaystyle=\frac{2\eta}{r}(v_{r}+\partial_{\theta}v_{\theta})+\left(\zeta-\frac{2}{3}\eta\right)\nabla\cdot\mathbf{v}$
$\displaystyle\sigma^{\prime}_{\phi\phi}$
$\displaystyle=\frac{2\eta}{r}(v_{r}+v_{\theta}\cot{\theta})+\left(\zeta-\frac{2}{3}\eta\right)\nabla\cdot\mathbf{v}$
$\displaystyle\sigma^{\prime}_{r\theta}=\sigma^{\prime}_{\theta r}$
$\displaystyle=\frac{\eta}{r}\left(\partial_{\theta}v_{r}+r\partial_{r}v_{\theta}-v_{\theta}\right)\,,$
(30)

where

$\nabla\cdot\mathbf{v}=\partial_{r}v_{r}+\frac{2}{r}v_{r}+\frac{1}{r\sin{\theta}}\partial_{\theta}(\sin{\theta}\,v_{\theta})\,.$
(31)

##  Appendix D Additional perturbed axisymmetric mode decomposition

The components of the perturbed velocity field are determined by the perturbed
continuity equation 18 as well as by the force-balance equations, which reduce
to

$\displaystyle\left(2\partial_{r}^{2}+\frac{4}{r}\partial_{r}-\frac{4+n(n+1)}{r^{2}}\right)\,v_{r}^{(n)}+\left(-\frac{3}{r^{2}}+\frac{1}{r}\partial_{r}\right)\,v_{\theta}^{(n)}$
$\displaystyle=\frac{1}{\eta}\partial_{r}\bar{P}^{(n)}$
$\displaystyle\left(r\partial_{r}^{2}+2\partial_{r}-\frac{2n(n+1)}{r}\right)\,v_{\theta}^{(n)}-n(n+1)\left(\frac{4}{r}+\partial_{r}\right)\,v_{r}^{(n)}$
$\displaystyle=-\frac{n(n+1)}{\eta}\bar{P}^{(n)}\,,$ (32)

where $\bar{P}^{(n)}=P^{(n)}-\left(\zeta-(2/3)\eta\right)\kappa\rho^{(n)}$.

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 $\delta R$, it is given by

$H=\frac{2}{R}-[\partial_{\theta}^{2}+(\cot{\theta})\partial_{\theta}+2]\frac{\delta
R}{R^{2}}\,.$ (33)

The different quantities evaluated at $r=R+\delta R$ read, to first order in
perturbations:

$\displaystyle(\sigma_{\rm nn})_{r=R+\delta R}$
$\displaystyle=\sigma_{rr}^{\rm stat}+\frac{d\sigma_{rr}^{\rm
stat}}{dr}\,\delta R+\delta\sigma_{rr}-2\sigma_{r\theta}^{\rm
stat}\,\frac{\partial_{\theta}\delta R}{R}$ $\displaystyle(\sigma_{\rm
nt})_{r=R+\delta R}$ $\displaystyle=\sigma_{r\theta}^{\rm
stat}+\frac{d\sigma_{r\theta}^{\rm stat}}{dr}\,\delta
R-(\sigma_{\theta\theta}^{\rm stat}-\sigma_{rr}^{\rm
stat})\frac{\partial_{\theta}\delta R}{R}+\delta\sigma_{r\theta}$
$\displaystyle\rho_{r=R+\delta R}$ $\displaystyle=\rho^{\rm
stat}+\frac{d\rho^{\rm stat}}{dr}\,\delta R+\delta\rho$
$\displaystyle(v_{r})_{r=R+\delta R}$ $\displaystyle=v_{r}^{\rm
stat}+\frac{dv_{r}^{\rm stat}}{dr}\,\delta R+\delta v_{r}\,.$ (34)

In the right-hand sides, all the quantities that depend on $r$ are evaluated
at $r=R$. Note that $\sigma_{r\theta}^{\rm stat}$ is formally written here but
is actually zero in the spherical symmetric case of our stationary state.
Similarly, the stationary velocity $v_{r}^{\rm stat}$ is also zero at the
outer stationary boundary $r=R$.

##  Appendix E Additional solutions of the perturbed axisymmetric mode
decomposition in the limit of fast diffusion

The components of the pressure and total stress tensor read

$\displaystyle P^{(n)}$
$\displaystyle=\left[\frac{(4n+6)\eta}{n}\,\frac{a_{n}}{R}+\frac{(n-6)\eta+3n\zeta}{3n}\,\kappa\,R\,c_{n}\right]\,\bar{r}^{n}$
$\displaystyle\sigma^{(n)}_{rr}$
$\displaystyle=\frac{\eta}{R}\left[\left(\frac{2n^{2}-2n-6}{n}\,a_{n}-\frac{n-2}{n}\,\kappa\,R\,c_{n}\right)\bar{r}^{n}+2(n-1)\,b_{n}\bar{r}^{n-2}\right]$
$\displaystyle\sigma^{(n)}_{r\theta}$
$\displaystyle=\frac{\eta}{R}\left[\left(-2n(n+2)\,a_{n}+n\,\kappa\,R\,c_{n}\right)\bar{r}^{n}-2(n^{2}-1)\,b_{n}\bar{r}^{n-2}\right]\,,$
(35)

where the perturbed pressure $\delta P$ and the perturbed component
$\delta\sigma_{rr}$ of the stress tensor are decomposed on the basis of
Legendre polynomials, and the perturbed component $\delta\sigma_{r\theta}$ of
the stress tensor is decomposed on the basis of Gegenbauer polynomials111Note
that the other non-trivial components of the stress tensor,
$\delta\sigma_{\theta\theta}$ and $\delta\sigma_{\phi\phi}$, have expansion
coefficients both on the Legendre and Gegenbauer polynomials and are not
reported here.. The three integration constants $a_{n}$, $b_{n}$, and $c_{n}$
are further determined by the three boundary conditions at the outer surface
of the multicellular spheroid, corresponding to the continuity of the normal
and tangential components of the stress tensor, as well as that of the
metabolite concentration.

Solving the boundary-condition equations D for the metabolite field with its
expression given by equation 4.3, we get

$c_{n}=-\frac{1}{3}\frac{k_{0}R}{\kappa\,l_{D}^{2}}\,R^{(n)}\,.$ (36)

Doing the same for the stress tensor gives an equation for $a_{n}$ and
$b_{n}$:

$\begin{pmatrix}\frac{n^{2}-n-3}{n}&n-1\\\
-\frac{n+2}{n+1}&-\frac{n-1}{n}\end{pmatrix}\cdot\begin{pmatrix}a_{n}\\\
b_{n}\end{pmatrix}=\begin{pmatrix}\left[1-\frac{n(n+1)}{2}\right]\frac{\gamma}{\eta\,R}-\frac{n-2}{6n}\,k_{0}\left(\frac{R}{l_{D}}\right)^{2}+2(\kappa\rho^{\rm
ext}-k_{0})\\\
\frac{1}{6(n+1)}\,k_{0}\left(\frac{R}{l_{D}}\right)^{2}+\kappa\rho^{\rm
ext}-k_{0}\end{pmatrix}\,R^{(n)}\,.$ (37)

Using further the kinematic condition equation 6 leads to the mode growth
rates given by equation 20.

For the modes $n=0$ and $n=1$, the expressions reported in section 4.3 can
further be expressed as functions of either $k_{0}$ and $\bar{R}=R/l_{D}$ only
or $\kappa\rho^{\rm ext}$ and $\bar{R}$ only, thanks to equation 12. This
leads to

$\displaystyle\omega_{0}$
$\displaystyle=k_{0}\left[\frac{\bar{R}^{2}}{3(\bar{R}\coth{\bar{R}}-1)}-1-\frac{1}{9}\bar{R}^{2}\right]=\kappa\rho^{\rm
ext}\left[\frac{4}{3}+\frac{3}{\bar{R}^{2}}-\frac{(9+\bar{R}^{2})\coth{\bar{R}}}{3\bar{R}}\right]$
$\displaystyle\omega_{1}$
$\displaystyle=k_{0}\left[\frac{\bar{R}^{2}}{9(\bar{R}\coth{\bar{R}}-1)}-\frac{1}{3}-\frac{1}{18}\bar{R}^{2}\right]=\kappa\rho^{\rm
ext}\left[\frac{1}{2}+\frac{2}{\bar{R}^{2}}-\frac{(6+\bar{R}^{2})\coth{\bar{R}}}{6\bar{R}}\right]\,.$
(38)

With these expressions, we can easily verify that these two modes are always
stable.

##  Appendix F Implicit equation for the mode $n=0$ in the generic case

Considering the generic case of equation 17, the growth rates $\omega_{n}$ are
given by implicit equations, which cannot be solved analytically. We report
here as an example the implicit equation giving the mode $n=0$, which
corresponds to the simplest one:

$\omega_{0}=k_{\rm p}^{\rm
ext}-\frac{1}{3}k_{0}\frac{\bar{R}^{2}}{\sinh{{\bar{R}_{0}}}}\left[\frac{\cosh{\bar{R}_{0}}}{\bar{R}_{0}}-\frac{\sinh{\bar{R}_{0}}}{{\bar{R}_{0}}^{2}}\right]\,,$
(39)

where $\bar{R}=R/l_{D}$ with $l_{D}=\sqrt{D/\alpha}$ as before, and
$\bar{R}_{0}=R/l_{0}$ with $l_{0}=\sqrt{D/(\alpha+\omega_{0})}$.

##  Appendix G Shapes of the lowest-order instability modes

We illustrate in figure 4 the shapes of the lowest-order instability modes.

Figure 4:  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 $\phi$, and all modes depicted here
are invariant under this transformation.

##  Appendix H Mode structure for intermediate values of the diffusion
coefficient of metabolites

We illustrate in figure 5 the mode structure for two values of the diffusion
coefficient of metabolites $D$, intermediate between those of figures 3b,c and
3d.

Figure 5:  Mode growth rates as functions of the mode number $n$ for
intermediate values of the diffusion coefficient of metabolites. Parameters
are the same as those used in figure 3b,c, using the same colour code, except
for the diffusion coefficient of metabolites $D$, which here takes the
respective values (a) $D=10\cdot 10^{-10}$ cm${}^{2}\cdot$s-1 and (b)
$D=20\cdot 10^{-10}$ cm${}^{2}\cdot$s-1, and the metabolite consumption rate
$\alpha$, which is scaled accordingly to keep the same metabolite penetration
lengths $l_{D}$ and associated steady states.

## References

## References

  * [1] Saffman P G, Taylor S G and S F R 1958 Proc. R. Soc. Lond. A 245 312–329 ISSN 0080-4630, 2053-9169
  * [2] Bensimon D, Kadanoff L P, Liang S, Shraiman B I and Tang C 1986 Rev. Mod. Phys. 58 977
  * [3] Mullins W W and Sekerka R F 1964 J. Appl. Phys. 35 444–451 ISSN 0021-8979
  * [4] Langer J S 1980 Rev. Mod. Phys. 52 1–28
  * [5] Rayleigh L 1900 Scientific Papers vol II (Cambridge University Press, Cambridge, England) p 200
  * [6] Taylor G I 1950 Proc. R. Soc. Lond. Ser. Math. Phys. Sci. 201 192–196
  * [7] Risler T, Peilloux A and Prost J 2015 Phys. Rev. Lett. 115 258104
  * [8] Kalziqi A, Yanni D, Thomas J, Ng S L, Vivek S, Hammer B K and Yunker P J 2018 Phys. Rev. Lett. 120 018101
  * [9] Williamson J J and Salbreux G 2018 Phys. Rev. Lett. 121 238102
  * [10] Alert R, Blanch-Mercader C and Casademunt J 2019 Phys. Rev. Lett. 122 088104
  * [11] Ben-Jacob E, Schochet O, Tenenbaum A, Cohen I, Czirók A and Vicsek T 1994 Nature 368 46–49 ISSN 0028-0836
  * [12] Ben-Jacob E, Cohen I and Levine H 2000 Adv. Phys. 49 395–554 ISSN 0001-8732, 1460-6976
  * [13] Sander L M 1986 Nature 322 789–793 ISSN 0028-0836
  * [14] Fujikawa H and Matsushita M 1989 J. Phys. Soc. Jpn. 58 3875–3878 ISSN 0031-9015
  * [15] Wang Q and Zhao X 2015 Sci. Rep. 5 8887 ISSN 2045-2322
  * [16] Li B, Cao Y P, Feng X Q and Gao H 2012 Soft Matter 8 5728 ISSN 1744-683X, 1744-6848
  * [17] Dervaux J and Amar M B 2008 Phys. Rev. Lett. 101 068101
  * [18] Liang H and Mahadevan L 2011 Proc. Natl. Acad. Sci. 108 5516–5521 ISSN 0027-8424, 1091-6490
  * [19] Savin T, Kurpios N A, Shyer A E, Florescu P, Liang H, Mahadevan L and Tabin C J 2011 Nature 476 57–62 ISSN 0028-0836
  * [20] Ben Amar M and Jia F 2013 Proc. Natl. Acad. Sci. 110 10525–10530 ISSN 0027-8424, 1091-6490
  * [21] Durel J F and Nerurkar N L 2020 Curr. Opin. Genet. Dev. 63 45–52 ISSN 0959-437X
  * [22] Shyer A E, Tallinen T, Nerurkar N L, Wei Z, Gil E S, Kaplan D L, Tabin C J and Mahadevan L 2013 Science 342 212–218 ISSN 0036-8075
  * [23] Tallinen T, Chung J Y, Rousseau F, Girard N, Lefèvre J and Mahadevan L 2016 Nat. Phys. 12 588–593 ISSN 1745-2481
  * [24] Drasdo D 2000 Phys. Rev. Lett. 84 4244
  * [25] Hannezo E, Prost J and Joanny J F 2012 Phys. Rev. Lett. 109 018101
  * [26] Kücken M and Newell A C 2005 J. Theor. Biol. 1 71–83 ISSN 0022-5193
  * [27] Tozluoǧlu M, Duda M, Kirkland N J, Barrientos R, Burden J J, Muñoz J J and Mao Y 2019 Dev. Cell 51 299–312.e4 ISSN 1534-5807
  * [28] Eng R C and Sampathkumar A 2020 Curr. Biol. 30 R402–R405 ISSN 0960-9822
  * [29] Sutherland R M 1988 Science 240 177–184 ISSN 0036-8075, 1095-9203
  * [30] Hirschhaeuser F, Menne H, Dittfeld C, West J, Mueller-Klieser W and Kunz-Schughart L A 2010 J. Biotechnol. 148 3–15 ISSN 0168-1656
  * [31] Costa E C, Moreira A F, de Melo-Diogo D, Gaspar V M, Carvalho M P and Correia I J 2016 Biotechnol. Adv. 34 1427–1441 ISSN 0734-9750
  * [32] Helmlinger G, Netti P A, Lichtenbeld H C, Melder R J and Jain R K 1997 Nat. Biotechnol. 15 778–783 ISSN 1087-0156
  * [33] Montel F, Delarue M, Elgeti J, Malaquin L, Basan M, Risler T, Cabane B, Vignjevic D, Prost J, Cappello G and Joanny J F 2011 Phys. Rev. Lett. 107 188102
  * [34] Delarue M, Montel F, Caen O, Elgeti J, Siaugue J M, Vignjevic D, Prost J, Joanny J F and Cappello G 2013 Phys. Rev. Lett. 110 138103
  * [35] Khalifat N, Beaune G, Nagarajan U, Winnik F M and Brochard-Wyart F 2016 Jpn. J. Appl. Phys. 55 1102A8 ISSN 1347-4065
  * [36] Dolega M E, Delarue M, Ingremeau F, Prost J, Delon A and Cappello G 2017 Nat. Commun. 8 14056 ISSN 2041-1723
  * [37] Forgacs G, Foty R A, Shafrir Y and Steinberg M S 1998 Biophysical Journal 74 2227–2234 ISSN 0006-3495
  * [38] Lecuit T and Lenne P F 2007 Nat. Rev. Mol. Cell Biol. 8 633–644 ISSN 1471-0072
  * [39] Guevorkian K, Colbert M J, Durth M, Dufour S and Brochard-Wyart F 2010 Phys. Rev. Lett. 104 218101
  * [40] Manning M L, Foty R A, Steinberg M S and Schötz E M 2010 Proc. Natl. Acad. Sci. 107 12517–12522 ISSN 0027-8424, 1091-6490
  * [41] Gonzalez-Rodriguez D, Guevorkian K, Douezan S and Brochard-Wyart F 2012 Science 338 910–917 ISSN 0036-8075, 1095-9203
  * [42] Dorie M J, Kallman R F, Rapacchietta D F, Van Antwerp D and Huang Y R 1982 Exp. Cell Res. 141 201–209 ISSN 0014-4827
  * [43] Jiang Y, Pjesivac-Grbovic J, Cantrell C and Freyer J P 2005 Biophys. J. 89 3884–3894 ISSN 0006-3495
  * [44] Byrne H M 2010 Nat. Rev. Cancer 10 221–230 ISSN 1474-175X
  * [45] Greenspan H P 1976 J. Theor. Biol. 56 229–242 ISSN 0022-5193
  * [46] Ferreira S C, Martins M L and Vilela M J 2002 Phys. Rev. E 65 021907
  * [47] Cristini V, Lowengrub J and Nie Q 2003 J. Math. Biol. 46 191–224 ISSN 0303-6812, 1432-1416
  * [48] Macklin P and Lowengrub J 2007 J. Theor. Biol. 4 677–704 ISSN 0022-5193
  * [49] Chaplain M A J, Ganesh M and Graham I G 2001 J. Math. Biol. 42 387–423 ISSN 0303-6812, 1432-1416
  * [50] Khain E and Sander L M 2006 Phys. Rev. Lett. 96 188103
  * [51] Byrne H M and Chaplain M A J 1996 Math. Comput. Model. 24 1–17 ISSN 0895-7177
  * [52] Frieboes H B, Zheng X, Sun C H, Tromberg B, Gatenby R and Cristini V 2006 Cancer Res. 66 1597–1604 ISSN 0008-5472, 1538-7445
  * [53] Chatelain C, Balois T, Ciarletta P and Ben Amar M 2011 New J. Phys. 13 115013 ISSN 1367-2630
  * [54] Drasdo D and Hoehme S 2012 New J. Phys. 14 055025 ISSN 1367-2630
  * [55] Ciarletta P 2013 Phys. Rev. Lett. 110 158102
  * [56] Marmottant P, Mgharbel A, Käfer J, Audren B, Rieu J P, Vial J C, van der Sanden B, Marée A F M, Graner F and Delanoë-Ayari H 2009 Proc. Natl. Acad. Sci. 106 17271–17275 ISSN 0027-8424, 1091-6490
  * [57] Ranft J, Basan M, Elgeti J, Joanny J F, Prost J and Jülicher F 2010 Proc. Natl. Acad. Sci. 107 20863–20868 ISSN 0027-8424, 1091-6490
  * [58] Popović M, Nandi A, Merkel M, Etournay R, Eaton S, Jülicher F and Salbreux G 2017 New J. Phys. 19 033006 ISSN 1367-2630
  * [59] Petridou N I, Grigolon S, Salbreux G, Hannezo E and Heisenberg C P 2019 Nat. Cell Biol. 21 169–178 ISSN 1476-4679
  * [60] Schötz E M, Burdine R D, Jülicher F, Steinberg M S, Heisenberg C P and Foty R A 2008 HFSP J. 2 42–56 ISSN 1955-2068
  * [61] Mgharbel A, Delanoë-Ayari H and Rieu J P 2009 HFSP J. 3 213–221 ISSN 1955-2068
  * [62] Montel F, Delarue M, Elgeti J, Vignjevic D, Cappello G and Prost J 2012 New J. Phys. 14 055008 ISSN 1367-2630
  * [63] Delarue M, Joanny J F, Jülicher F and Prost J 2014 Interface Focus 4 20140033 ISSN 2042-8898, 2042-8901
  * [64] Basan M, Risler T, Joanny J F, Sastre-Garau X and Prost J 2009 HFSP J. 3 265–272 ISSN 1955-2068
  * [65] Happel J and Brenner H 2012 Low Reynolds Number Hydrodynamics: With Special Applications to Particulate Media (Springer Science & Business Media) ISBN 978-94-009-8352-6
  * [66] Basan M, Joanny J F, Prost J and Risler T 2011 Phys. Rev. Lett. 106 158101
  * [67] Risler T and Basan M 2013 New J. Phys. 15 065011 ISSN 1367-2630
  * [68] Freyer J P and Sutherland R M 1983 Determination of Diffusion Constants for Metabolites in Multicell Tumor Spheroids Oxygen Transport to Tissue—IV Advances in Experimental Medicine and Biology ed Bicher H I and Bruley D F (Boston, MA: Springer US) pp 463–475 ISBN 978-1-4684-7790-0
  * [69] Casciari J J, Sotirchos S V and Sutherland R M 1988 Cancer Res. 48 3905–3909 ISSN 0008-5472, 1538-7445
  * [70] Mueller-Klieser W 1984 Biophys. J. 46 343–348 ISSN 0006-3495
  * [71] Salbreux G, Joanny J F, Prost J and Pullarkat P 2007 Phys. Biol. 4 268 ISSN 1478-3975
  * [72] Callan-Jones A C, Joanny J F and Prost J 2008 Phys. Rev. Lett. 100 258106
  * [73] Blanch-Mercader C and Casademunt J 2013 Phys. Rev. Lett. 110 078102
  * [74] Mongera A, Rowghanian P, Gustafson H J, Shelton E, Kealhofer D A, Carn E K, Serwane F, Lucio A A, Giammona J and Campàs O 2018 Nature 561 401–405 ISSN 1476-4687
  * [75] Hashmi A, Tlili S, Perrin P, Martinez-Arias A and Lenne P F 2020 bioRxiv 2020.05.21.105551
  * [76] Aguirre-Ghiso J A 2007 Nat. Rev. Cancer 7 834–846 ISSN 1474-175X
  * [77] Uhr J W, Scheuermann R H, Street N E and Vitetta E S 1997 Nat. Med. 3 505–509
  * [78] Karrison T G, Ferguson D J and Meier P 1999 J. Natl. Cancer Inst. 91 80–85 ISSN 0027-8874
  * [79] Black W C and Welch H G 1993 N. Engl. J. Med. 328 1237–1243 ISSN 0028-4793
  * [80] Folkman J and Kalluri R 2004 Nature 427 787–787 ISSN 0028-0836
  * [81] Naumov G N, Akslen L A and Folkman J 2006 Cell Cycle 5 1779–1787 ISSN 1538-4101
  * [82] Hanahan D and Folkman J 1996 Cell 86 353–364
  * [83] Semenza G L 2003 Annu. Rev. Med. 54 17–28 ISSN 0066-4219
  * [84] Preziosi L 2003 Cancer Modelling and Simulation (CRC Press) ISBN 978-0-203-49489-9
  * [85] Tracqui P 2009 Rep. Prog. Phys. 72 056701 ISSN 0034-4885
  * [86] Risler T 2015 New J. Phys. 17 055011 ISSN 1367-2630
  * [87] Scott D and Weston K 2015 Converg. Sci. Phys. Oncol. 1 026002 ISSN 2057-1739
  * [88] Porta C L and Zapperi S 2017 The Physics of Cancer (Cambridge University Press) ISBN 978-1-108-15033-0
  * [89] Lamprecht S, Schmidt E M, Blaj C, Hermeking H, Jung A, Kirchner T and Horst D 2017 Nat. Commun. 8 1406 ISSN 2041-1723
  * [90] Lenos K J, Miedema D M, Lodestijn S C, Nijman L E, van den Bosch T, Ros X R, Lourenço F C, Lecca M C, van der Heijden M, van Neerven S M, van Oort A, Leveille N, Adam R S, Melo F d S E, Otten J, Veerman P, Hypolite G, Koens L, Lyons S K, Stassi G, Winton D J, Medema J P, Morrissey E, Bijlsma M F and Vermeulen L 2018 Nat. Cell Biol. 20 1193–1202 ISSN 1476-4679
  * [91] van der Heijden M, Miedema D M, Waclaw B, Veenstra V L, Lecca M C, Nijman L E, van Dijk E, van Neerven S M, Lodestijn S C, Lenos K J, de Groot N E, Prasetyanti P R, Varea A A, Winton D J, Medema J P, Morrissey E, Ylstra B, Nowak M A, Bijlsma M F and Vermeulen L 2019 Proc. Natl. Acad. Sci. 116 6140–6145 ISSN 0027-8424, 1091-6490