Viscocapillary Instability in Cellular Spheroids
VViscocapillary Instability in Cellular Spheroids
Matthieu Martin and Thomas Risler Laboratoire Physico-Chimie Curie, Institut Curie, PSL Research University,Sorbonne Universit´e, CNRS, 26 rue d’Ulm, 75005 Paris, FranceE-mail: [email protected]
Abstract.
We describe a viscocapillary instability that can perturb the sphericalsymmetry of cellular aggregates in culture, also called multicellular spheroids. In thecondition where the cells constituting the spheroid get their necessary metabolitesfrom the immediate, outer microenvironment, a permanent cell flow exists withinthe 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 thetissue that can destabilise the aggregate. The proposed instability is viscocapillary innature and does not rely on external heterogeneities, such as a pre-existing patternof blood vessels or the presence of a substrate on which the cells can exert pullingforces. It arises for sufficiently large cell-cell adhesion strengths, cell-renewal rates,and metabolite supplies, as described by our model parameters. Since multicellularspheroids in culture are good model systems of small, avascular tumours, mimicking themetabolite concentration gradients found in vivo , we can speculate that our descriptionapplies to microtumour instabilities in cancer progression.
Submitted to:
New J. Phys. a r X i v : . [ phy s i c s . b i o - ph ] F e b iscocapillary Instability in Cellular Spheroids
1. Introduction
Interface instabilities in systems driven far from equilibrium have been extensivelystudied 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-Shawcell [1, 2], the Mullins-Sekerka instability, which stems from the diffusive transport ofthe latent heat of solidification in unidirectional solidification [3, 4], and the Rayleigh-Taylor instability, corresponding to the fingering of an interface between two immisciblefluids 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 theseclassical condensed-matter instabilities have been identified in living systems. In tissuesor 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 fingeringhave been observed [11, 12]. Such a coupling can lead to fractal branching patterns viathe process of diffusion-limited aggregation [13, 14] or other types of branching patterns,depending on the bacterial morphotype [12]. In tissues, mechanical instabilities havebeen recognised to play a potential role in different morphological processes and patternsexhibited by growing cell populations [15]. Examples are the wrinkling patterns ofgrowing, soft surfaces [16, 15] such as those of leaves and flowers [17, 18], the large-scalelooping 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 abuckling phenomenon [24, 25], potentially driven by the differential growth of adjacenttissue layers [26, 18, 27, 28], or by other curling and crumpling instabilities due toanisotropic growth [17].Here, we focus on the stability of the spherical growth of cellular aggregates. Suchexperimental systems are used to study the growth dynamics and cellular structure ofmicroscopic tumours with realistic metabolite concentration gradients [29]. They havealso been used as anti-cancer therapy test platforms, mirroring the three dimensionalcellular context and therapeutically relevant pathophysiological gradients of in-vivo tumours [30, 31]. Cellular aggregates also permit the study of the effects of differentperturbations on the growth or cellular-duplication dynamics, such as a change in theexternal mechanical constraints [32, 33, 34, 35, 36]. Most often, an effective surfacetension exists between the aggregate and its direct environment, which makes it look likea spheroid, justifying the denomination of multicellular spheroid [37, 38, 39, 40, 41]. Thesupply of metabolites from the microenvironment is responsible for an inhomogeneousdistribution of cell divisions within the spheroid, with an increased cell proliferation atits periphery and an increased cell death in its core [30]. As a consequence, cells flowfrom 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 sphericalshape can be unstable without changing the average spheroid size, due to the shear iscocapillary Instability in Cellular Spheroids
2. 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 achemical substance necessary for cellular proliferation. This substance can be oxygen,growth factors, glucose or other nutrients, and will be referred to with the generic termof ‘metabolites’ in the following. The tissue within the spheroid is characterised in acontinuum theory by a cell-number density and a cell-velocity field v . Analyses of thestress distribution in multicellular spheroids have shown that, at steady state underphysiological osmotic conditions, the cell density is essentially homogeneous throughoutone aggregate [62, 63]. We therefore consider an incompressible tissue, for which thecell-number density is constant and the continuity equation reads ∇ · v = k p . (1) iscocapillary Instability in Cellular Spheroids AB
Schematic drawing of the proposed viscocapillary instability in amulticellular spheroid. The original spherical shape is indicated by the external blue-dashed circle. The scaled metabolite concentration is depicted as a scaled colourgradient. Dividing cells (red) are located mostly close to the outer surface. Closeto the centre, cells are deprived of metabolites and die (yellow dashed). An inner bluecircle indicates cells equidistant from the spheroid centre, and the tissue-flow field isqualitatively depicted by white (resp. black) arrows for inward (resp. outward) cellvelocities. At point A , cells are pushed inward faster than at point B , because morelayers of dividing cells push the cells in A toward the spheroid core. This createsviscous shear stresses that, in reaction, push the cells located above A outwards.Capillary effects due to tissue surface tension at the outer boundary, however, favourstability. Depending on parameter values, the original spherical shape can or notbecome unstable by this viscocapillary mechanism. 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 togive its divergence. Neglecting inertia and in the absence of bulk external forces, forcebalance reads ∇ · σ = 0 , (2)where σ denotes the total stress tensor. We further split the stress tensor into adynamic part and a velocity-independent part. For an isotropic tissue, the latter reads − P , where P is the tissue pressure and the unity tensor. The dynamic part σ (cid:48) encodes the rheological properties of the tissue. The timescale of interest here is largecompared to those of individual cellular processes, such as cell-cell rearrangementsand cell renewal. We can therefore model the tissue as a viscous fluid with effectiveshear and bulk viscosities η and ζ , taking into account the long-term effects of cellproduction [37, 56, 57]. We obtain σ (cid:48) = η (cid:20) ∇ ⊗ v + ( ∇ ⊗ v ) T −
23 ( ∇ · v ) (cid:21) + ζ ( ∇ · v ) , (3) iscocapillary Instability in Cellular Spheroids ⊗ denotes the tensorial product and ( ∇ ⊗ v ) T the transposed tensor of ∇ ⊗ v .The remarkable absence of any compression modulus is due to the nonconservation ofcell number [64, 57].The system of equations is closed by specifying the expression of the cell-productionrate k p . For simplicity, we assume that k p is independent of the stress σ and linearlydependent on the metabolite concentration ρ . This leads to k p = κρ − k , (4)where κ and k 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 metaboliteconcentration with a constant absorption rate α , as it has been observed in avasculartumour growth [43]: ∂ t ρ = D ∆ ρ − αρ . (5)Here, ∂ t denotes the partial time derivative and ∆ the Laplacian operator. Note that theconvective term has been ignored in this equation, which is justified if nutrient diffusionis fast compared to its convective transport by the cell flow.In the center of the spheroid r = , the cell-velocity field vanishes and the metaboliteconcentration remains positive. At the outer surface, the metabolite concentrationequals the external concentration ρ ext and the cell velocity equals that of the interface.Labelling R the spheroid stationary radius and δR its perturbation, these conditionsread ρ r = R + δR = ρ ext and ∂ t ( δR ) = v r = R + δR · e r , (6)where e r is the unit radial vector and v is evaluated at the perturbed interfacelocation. Finally, the outer surface is subjected solely to the isotropic, external pressure.The tangential component σ nt of the stress tensor therefore vanishes and its normalcomponent is given by Laplace’s law with surface tension γ : ( σ nn ) r = R + δR = − P ext − γ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 boundarycondition.
3. Stationary solution
The stationary equations are characterised by spherical symmetry with radial coordinate r . The stationary continuity equation 1 and metabolite-diffusion equation 5 reduce to:1 r ∂ r ( r v r ) = k p (7)and Dr ∂ r ( r ∂ r ρ ) = αρ , (8) iscocapillary Instability in Cellular Spheroids ρ stat ) r = R = ρ ext and ( v stat r ) r = R =0. The stationary force-balance condition and expressions of the stress-tensorcomponents are given in Appendix A. To integrate this system of equations, we start by integrating the metabolite diffusionequation 8. With the characteristic metabolite-penetration length l D = (cid:112) D/α andthe reduced variables ¯ r = r/l D and ¯ R = R/l D , the metabolite concentration and thecell-production rate read ρ stat = ρ ext ¯ R sinh ¯ R sinh ¯ r ¯ r (9)and k statp = κρ ext ¯ R sinh ¯ R sinh ¯ r ¯ r − k . (10)We can see from Eq. 10 that a stationary solution with spherical symmetry exists aslong as κρ ext > k to have a positive cell-division rate at the outer rim of the spheroid.Finally, the cell-velocity field reads v stat r = κρ ext R sinh ¯ R (cid:20) cosh ¯ r ¯ r − sinh ¯ r ¯ r (cid:21) − k r . (11)The boundary condition ( v stat r ) r = R = 0 then leads to the following equation for thestationary radius: 1¯ R coth ¯ R − R = k κρ ext . (12)Note that the product κρ ext is linked to k and the cell-production rate at the outer rim k extp by k extp = κρ ext − k . The complete expressions of all the other quantities in thisstationary state are given in Appendix B.
4. Mode computation
We now investigate the linear stability of this stationary solution to axisymmetricperturbations. We choose a system of coordinates composed of radial coordinate r , polarangle θ , and azimuthal angle φ , and we study the perturbations with axial symmetryaround θ = 0. The model equations with axial symmetry read1 r ∂ r ( r v r ) + 1 r sin θ ∂ θ (sin θ v θ ) = k p (13)for the continuity equation 1 and ∂ t ρ = Dr ∂ r ( r ∂ r ρ ) + Dr sin θ ∂ θ (sin θ∂ θ ρ ) − αρ (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 Appendix C. iscocapillary Instability in Cellular Spheroids To integrate this system of equations to linear order in perturbations, we expand theangular dependence of the different perturbative fields onto the basis of axisymmetric,spherical harmonics. Following ref. [65], the perturbations δR , δv r , δP , and δρ areexpanded 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 − − P n +1 ) / (2 n + 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 , asthe components R ( n ) of the expansion of the interface location δR are simple numbers.Explicitly, the δv r expansion for example reads δv r ( r, θ ) = ∞ (cid:88) n =0 v ( n ) r P n (cos θ ) , (15)as the δv θ expansion reads δv θ ( r, θ ) = ∞ (cid:88) n =1 v ( n ) θ I n +1 (cos θ )sin θ , (16)where the sum starts at n = 1 since the mode n = 0 is purely radial. Using theseexpansions, the metabolite diffusion equation 5 leads to ∂ r ρ ( n ) + 2 r ∂ r ρ ( n ) − (cid:20) α + ω n D + n ( n + 1) r (cid:21) ρ ( n ) = 0 , (17)where ρ ( n ) is the component for the mode number n of the nutrient field ρ and ω n thecorresponding growth rate. The components of the perturbed velocity field are thendetermined by the perturbed continuity equation ∂ r v ( n ) r + 2 r v ( n ) r + 1 r v ( n ) θ = κ ρ ( n ) (18)as well as by the force-balance equations, further given in Appendix D. To solve equation 17, we first consider the regime where the relaxation or growth ofthe perturbation modes as well as metabolite consumption are slow compared withmetabolite diffusion over the characteristic lengths involved, at most equal to thespheroid radius. In that limit, the term [( α + ω n ) /D ] ρ ( n ) can be neglected in frontof [ n ( n + 1) /r )] ρ ( 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 wellfor n = 1. For n ≥
2, the solution of equation 17 in this approximation is a simplepower law, and we can further obtain the other perturbed quantities analytically. We iscocapillary Instability in Cellular Spheroids n ≥ ρ ( n ) = c n ¯ r n v ( n ) r = a n ¯ r n +1 + b n ¯ r n − v ( n ) θ = [ − ( n + 3) a n + κ R c n ] ¯ r n +1 − ( n + 1) b n ¯ r n − , (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 linearcombinations of these three integration constants. The other obtained expressions aregiven in Appendix E.The integration contants are then determined using the boundary conditions.Plugging these solutions into the kinematic equation 6, we finally get the followingmode growth rates: ω n = 2 n + 5 n + 32 n + 4 n + 3 k extp − n + 4 n + 3 (cid:34) n + 13 k (cid:18) Rl D (cid:19) + (2 n + 5 n + 2) γ η nR (cid:35) (20)for all n ≥
2. For the modes n = 0 and n = 1, we get separately ω = k extp − (1 / k ¯ R and ω = (1 / k extp − (1 / k ¯ R . These latter expressions are independent of the surfacetension γ . This is because, for these two modes, perturbations in the curvature occuronly to second order or higher.
5. Results
We can now discuss the instability in this fast-diffusion regime. The first term inequation 20 is destabilising and proportional to the cell-production rate at the outersurface. The second term is stabilising and results for one part from cell-death processescontrolled by the parameters α/D (via l D ) and k , and from the other part by surfacetension. Increasing the viscosity lowers the contribution of surface tension, destabilisingthe spheroid. This underlines the mechanical origin of the instability, which relieson internal viscous stresses within the tissue, generated by differential cell flows (seefigure 1). Therefore, this instability can only exist around a kinematic steady state withnonzero 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 Appendix E, equation E.4). For the othermodes, in the absence of surface tension, the first term in equation 20 is dominant atlarge n , meaning that, without this contribution, the spheroid is always unstable witha rate asymptotically equal to that of cell division at its outer rim. Surface tensionhowever stabilises the spheroid, since it contributes by a term scaling as γq n / (2 η ) atlarge n , where q n = n/R . Depending on the values of the different parameters, wetherefore expect a potential instability to develop at a finite value of n , correspondingto a finite wavelength λ n ∼ πR/n . iscocapillary Instability in Cellular Spheroids R . Since the equation characterising this mode is in generalfourth order in n , we investigate separately the limits of small and large spheroids.In the limit of small radii, the least stable mode is ω , since curvature is large andsurface tension strongly stabilises all modes for larger values of n . In the limit oflarge radii, the most unstable mode occurs at n = n max with the following asymptoticexpression, linear in R : n max (cid:39) (cid:112) ηκρ ext / ( γl D ) R . This scaling indicates that theassociated 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 whathas been proposed for epithelial tissues [66, 67]. The corresponding growth rate reads ω max (cid:39) κρ ext − (cid:112) γκρ ext / ( ηl D ).In the generic case where metabolite diffusion is not necessarily fast compared tometabolite consumption or perturbation growth, the solution for ρ ( n ) in equation 17is a function of the associated growth rate ω n . Equation 6 then becomes an implicitequation for the growth rate ω n , which cannot be solved analytically. We report theimplicit 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 differentparameter values. The shear viscosity of cellular aggregates has been estimated indifferent experiments, leading to η (cid:39) − Pa · s [37, 56, 39]. Tissue surface tensionshave been measured for different tissue types. Measurements for γ give values rangingfrom a fraction up to several millinewton per meter [37, 60, 61, 39]. We further assumea typical cell-division rate at the outer surface of the spheroid k extp and a cell-death ratein the absence of metabolites k of one per day. Cellular growth within the spheroidcan be limited by different types of metabolites. Depending typically on the molecularsize of a given metabolite, its diffusion coefficient can take different values. Estimates ofthe diffusion coefficient of growth factors and glucose in avascular tumours range from10 − cm · h − for growth factors [43] to 1.5 10 − cm · h − for glucose [68, 69, 43], or evenlarger values for oxygen [70, 43]. In the following, we shall investigate the influenceof a variation of this particular parameter. Finally, to obtain radii of a few hundredmicrometres, we choose to have comparable values for the characteristic penetrationlength of metabolites l D . This leads to values for the metabolite-consumption rate α of several tenths per day for D ∼ − cm · h − , scaled accordingly when D is varied tokeep 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 outersurface and disappear in the centre (panel (a)), due to the lack of metabolites penetratingthe 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 ofcellular flows in multicellular spheroids using fluorescently labeled particles [34].We illustrate in figure 3 the central result of our study, that is the obtained modestructure of the instability, for the four different steady states presented in figure 2a iscocapillary Instability in Cellular Spheroids
250 r ( μ m ) v r ( μ m / d a y ) Figure 2.
Stationary state of a multicellular spheroid. (a,b) The cell-production rate k p (a) and the radial cell-velocity field v r (b) are displayed as functions of the radialdistance r for four different values of the stationary radius R . Parameters common tothe four sets of curves are: η = 80 kPa · s, γ = 100 µ N · m − , and k extp = k = 1 d − (d − stems for “per day”). The stationary state depends on the diffusion coefficient ofmetabolites D and metabolite-absorption rate α through their ratio D/α = l D only.We display the resulting stationary state for four different values of this ratio. Thecorresponding stationary radii and metabolite penetration lengths are R (cid:39) µ mand l D (cid:39) µ m (blue curves), R (cid:39) µ m and l D (cid:39) µ m (orange curves), R (cid:39) µ m and l D (cid:39) µ m (green curves), and R (cid:39) µ m and l D (cid:39) µ m (redcurves). (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 ρ/ρ ext and in panel (d) the corresponding radial cell-velocity field v r (in µ m · d − ), forthe parameter set corresponding to the blue curves of panels (a) and (b). Plots aremade in a plane of symmetry of the spheroid. and b, using the same colour code. The system is unstable as soon as at least oneperturbation mode displays a growth rate ω with a positive real part. When this isthe case, we expect the fastest growing modes to develop first, corresponding to themaximum of each series of points presented in these plots. A visual display of theshapes associated with the deformation modes n = 0 to 5 is shown in Appendix G,figure G1. Figure 3a shows the mode structure in the approximation of the analyticsolution of section 4.3. There is a transition at finite wavelength as a function of themetabolite-consumption rate α from a stable to an unstable regime, with a range ofunstable modes. These appear for example when nutrient consumption decreases atfixed external cell-division rate, surface tension and viscosity. With the parameterschosen here, the first mode to become unstable is n = 5 (figure 3a, orange squares). iscocapillary Instability in Cellular Spheroids - - - R e [ ω ] ( / d a y ) d)
Viscocapillary mode growth rates of a multicellular spheroid as functions ofthe mode number n , obtained (a) in the analytic limit of section 4.3, (b,c,d) numericallyusing the full model in both limits of slow (b,c) and fast (d) diffusion. (a) Parametersare the same as those of figure 2 using the same colour code. The mode-growthrates depend on the diffusion coefficient of metabolites D and metabolite-absorptionrate α through their ratio D/α = l D only. We have l D (cid:39) µ m (blue circles), l D (cid:39) µ m (orange squares), l D (cid:39) µ m (green diamonds), and l D (cid:39) µ m (redtriangles). (b,c) Real (b) and imaginary (c) parts of the mode-growth rates computednumerically with a diffusion coefficient of metabolites equal to that of growth factors: D = 5 · − cm · s − [43]. The metabolite-absorption rate is varied accordingly to keepthe same steady state as in figure 2: α = 0 .
15 d − (blue circles), α = 0 .
25 d − (orangesquares), α = 0 .
35 d − (green diamonds), and α = 0 .
45 d − (red triangles). (d) Mode-growth rates computed numerically with a diffusion coefficient of metabolite equal tothat of glucose: D = 5 · − cm · s − [68, 69, 43]. The metabolite-absorption rate isscaled accordingly by a factor of a thousand to perturb around the same steady stateas previously: α = 150 d − (blue circles), α = 250 d − (orange squares), α = 350 d − (green diamonds), α = 450 d − (red triangles). The mode-growth rates in (a) and (d)are real. This is associated with a stationary radius R (cid:39) µ m and corresponds to an unstablewavelength λ (cid:39) µ m.In figure 3b and c, we display respectively the real and imaginary parts of the modegrowth rates obtained in the full model solved numerically, with the parameter sets offigure 2 and a diffusion coefficient D = 5 · − cm · s − , corresponding approximatelyto that estimated for growth factors [43]. The system also displays first an instability atfinite wavelength (for the mode n = 3, green diamonds), but the modes n = 0 and n = 1can now be unstable. In addition, these instabilities are oscillatory, as characterised bynon-zero imaginary parts of ω , corresponding to the characteristic frequencies of the iscocapillary Instability in Cellular Spheroids D/l D be larger than the amplitude v of the cellular flow within the aggregate. Figure 3d displays the mode growthrates obtained with a diffusion coefficient D a thousand times that of figure 3b,c,corresponding approximately to the diffusion coefficient of small nutrient moleculessuch as glucose [68, 69, 43]. With this value of the diffusion coefficient, we have( (cid:96) D · v ) /D ∼ . − (cid:28) v is estimated from the curves shown in figure 2b in theleast favorable case, satisfying largely the required condition for neglecting convectiveflows. We show in figure 3d that, perturbing around the same steady state by rescalingthe metabolite-absorbing rate α , the instability occurs at a similar radius and finitewavelength as those reported above. In addition, the growth rates of the high-ordermodes ( n ≥ n = 0 and n = 1. This is the signature of the fact that here metabolite diffusionis sufficiently fast to allow for an almost instantaneous response of the inner cells toperturbations of the outer surface. We however do not recover the analytic resultsof figure 3a for small values of n . This stems from the fact that, in the results offigure 3d, we scale the metabolite consumption rate with the diffusion coefficient tokeep the same steady-state sizes as in figure 3b, as the analytic limit was obtained forsmall values of α/D with respect to 1 /R . We show in Appendix H, figure H1, themode structure obtained with intermediate values of the diffusion coefficient D betweenthose of figures 3b,c and 3d, following the same rescaling procedure of the metabolite-consumption rate α .As mentioned at the end of section 2, our boundary conditions specify that thecell-velocity field vanishes at the center of the spheroid located at r = , whichbreaks Galilean invariance. Doing so, the mode n = 1 here corresponds to an actualdeformation of the flow pattern within the spheroid, with an outer boundary displacedwith respect to the point where the flow pattern converges. As a consequence, this modeis not necessarily marginal, contrary to many standard spherical-harmonic perturbationanalyses. A similar interpretation of the mode n = 1 is found, e.g., in the context of thedeformations of the actin cortex of a spherical cell [71], where it corresponds to a relativetranslation of the inner part of the cell cortex with respect to the outer boundary ofthe cell, rather than to a global translation of the whole system. Therefore, even ifthis mode corresponds to a global translation of the inner boundary of the cortex withno deformation, the relative positions of the inner and outer boundaries of the cortexdo 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 displacementsof the different cell layers with each other. Each cell layer is purely translated in themode 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 iscocapillary Instability in Cellular Spheroids
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 metabolitesdiffusing from their microenvironment. The present instability relies neither on cellmotility nor on the presence of external forces, but rather stems from the presenceof viscous shear stresses generated by the spatial organisation of cell renewal withinthe tissue. We have shown that the instability develops at a finite wavenumber,which reflects a balance of viscous shear stresses with those stemming from surfacetension. We propose that this mechanism could be observed in multicellular spheroids inculture, which would be an ideal system for testing the influence of different parameterscontrolling the instability, such as tissue viscosity, surface tension, and metabolitesupply. The former two could be changed, e.g., by varying the expression of proteinsimplicated in cell-surface adhesion or actin-cortex contractility [38, 60, 40].The proposed instability here evokes other already proposed instabilities in thecontext of the cell cytoskeleton, driven by actin-polymerisation dynamics [72, 73]. Inthese studies of the stability of cell fragments on a substrate, an inward flow is drivenby actin polymerisation at the outer edge and actin depolymerisation in the bulk. Asa result, an originally circular cell fragment can become unstable and spontaneouslyacquire a polarisation. Important differences between our current study and theseprevious works however exist. In the stability analysis of circular cell fragments, thegeneration 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 continuouslywithin the spheroid. The second difference is the rheology, which here corresponds aStoke flow with no contact with an external substrate, as in the case of cell fragmentsthere 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 worksare two dimensional. As a result, our flow pattern has no anchoring to the external worldand requires a minimal thickness over which cells divide to become unstable.In addition to be applicable to multicellular aggregates, one can wonder if similar iscocapillary Instability in Cellular Spheroids in vivo . During development, transition from solid-like to fluid-like tissue properties, tissue surface tension, and flow patterns have been shown toplay a crucial role (see, e.g., [74]). Recently, three-dimensional aggregates of mouseembryonic stem cells have been shown to undergo a first morphological transformationfrom a spherical into an oblong shape during gastrulation, associated with a reducedlevel of E-cadherin expression at the developing tip [75]. Such shapes resemble asuperposition of instability modes such as n = 2 and 3, and potentially higher, asillustrated in figure G1. Interestingly with respect to our current study, the polarisationof E-cadherin expression precedes the onset of tip formation, and when the level ofE-cadherin expression is maintained high, the aggregate remains generally devoid ofany pole [75]. These observations suggest a role for a reduced surface tension in thedevelopment of the protrusion, similar to what we are proposing here. However, inthese examples, and to our knowledge more generally during development, it seems thatan original inhomogeneity in the tissue rheological parameters—such as surface tensionand viscosity—is at the origin of the shape formation. Such inhomogeneities are howeverabsent in our proposed mechanism, which relies solely on the presence of a permanentflow of duplicating cells.A domain to which we can speculate that the present mechanism applies is theevolution of microtumours after a long period of dormancy. Small primary tumoursor early metastases often enter a state where their sizes remain steady, before theyresume growth or disappear [76]. Such a dormant state can last for a long timeand is at the origin of late cancer reappearance, sometimes years after the originaltreatment [77, 78]. It is recognised that microscopic, clinically occult tumours arevery common in the population and that only a tiny fraction of them ever becomesclinically relevant [79, 80, 81]. Understanding the factors that can destabilise a dormanttumour is therefore of crucial importance. The main mechanisms at the origin ofsuch steady states are angiogenic dormancy, cellular dormancy (G0-G1 arrest) andimmunosurveillance [76]. In angiogenic dormancy, the tumour is limited in its growthby the lack of metabolites, which are brought by blood vessels that do not penetrate thetumour [80, 81, 76]. This limitation keeps the microtumour to sizes typically smallerthan 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 biologicalprocesses [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 theirsurfaces [89, 90], which corresponds to the steady-state patterns of cell duplicationson which our current study relies. In addition, clonal expansion largely depends on thelocation of a clone within the tumour [91], suggesting that differences in geometricalor physical properties within the tumour are major contributors to heterogeneousclonal expansion. This latter observation leads us to speculate that, after a firstinstability such as the one proposed here, the resulting irregular shape creates differentmicroenvironments for different parts of the tumour, further driving different epigeneticand maybe even later genetic transformations by diverse selection processes. iscocapillary Instability in Cellular Spheroids Acknowledgments
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 workreceived 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 ∂ r σ rr + 2 r ( σ rr − σ θθ ) = 0 , (A.1)and the non-zero components of the dynamic part of the stress tensor σ (cid:48) as given byequation 3 read σ (cid:48) rr = 2 η ∂ r v r + (cid:18) ζ − η (cid:19) ∇ · v σ (cid:48) θθ = σ (cid:48) φφ = 2 η v r r + (cid:18) ζ − η (cid:19) ∇ · v , (A.2)where ∇ · v = ∂ r v r + 2 r v r . (A.3)The corresponding boundary condition reads( σ stat rr ) r = R = − P ext − γR . (A.4) Appendix B. Additional stationary expressions
The non-trivial components of the dynamic part of the stress tensor read σ (cid:48) stat rr = κρ ext ¯ R sinh ¯ R (cid:20) − η cosh ¯ r ¯ r + (cid:18) ζ + 43 η + 4 η ¯ r (cid:19) sinh ¯ r ¯ r (cid:21) − ζk (B.1) iscocapillary Instability in Cellular Spheroids σ (cid:48) stat θθ = κρ ext ¯ R sinh ¯ R (cid:20) η cosh ¯ r ¯ r + (cid:18) ζ − η − η ¯ r (cid:19) sinh ¯ r ¯ r (cid:21) − ζk . (B.2)The pressure is given by P stat = P stat¯ r =0 + κρ ext (cid:18) ζ + 43 η (cid:19) ¯ R sinh( ¯ R ) (cid:18) sinh(¯ r )¯ r − (cid:19) , (B.3)where P stat¯ r =0 = P ext + 2 γR − ζk + κρ ext ¯ R (cid:20) η (cid:0) − ¯ R coth ¯ R (cid:1) + (cid:18) ζ + 43 η (cid:19) ¯ R sinh ¯ R (cid:21) . (B.4) Appendix C. Additional perturbed axisymmetric equations
With axisymmetry, the non-trivial components of the force-balance equation 2 read ∂ r σ rr + 1 r ∂ θ σ rθ + 1 r [2 σ rr − σ θθ − σ φφ + σ rθ cot θ ] = 0 ∂ r σ rθ + 1 r ∂ θ σ θθ + 1 r [3 σ rθ + ( σ θθ − σ φφ ) cot θ ] = 0 , (C.1)and the non-zero components of the dynamic part of the stress tensor, as given byequation 3, are σ (cid:48) rr = 2 η ∂ r v r + (cid:18) ζ − η (cid:19) ∇ · v σ (cid:48) θθ = 2 ηr ( v r + ∂ θ v θ ) + (cid:18) ζ − η (cid:19) ∇ · v σ (cid:48) φφ = 2 ηr ( v r + v θ cot θ ) + (cid:18) ζ − η (cid:19) ∇ · v σ (cid:48) rθ = σ (cid:48) θr = ηr ( ∂ θ v r + r∂ r v θ − v θ ) , (C.2)where ∇ · v = ∂ r v r + 2 r v r + 1 r sin θ ∂ θ (sin θ v θ ) . (C.3) Appendix D. Additional perturbed axisymmetric mode decomposition
The components of the perturbed velocity field are determined by the perturbedcontinuity equation 18 as well as by the force-balance equations, which reduce to (cid:18) ∂ r + 4 r ∂ r − n ( n + 1) r (cid:19) v ( n ) r + (cid:18) − r + 1 r ∂ r (cid:19) v ( n ) θ = 1 η ∂ r ¯ P ( n ) (cid:18) r∂ r + 2 ∂ r − n ( n + 1) r (cid:19) v ( n ) θ − n ( n + 1) (cid:18) r + ∂ r (cid:19) v ( n ) r = − n ( n + 1) η ¯ P ( n ) , (D.1)where ¯ P ( n ) = P ( n ) − ( ζ − (2 / η ) κρ ( n ) . iscocapillary Instability in Cellular Spheroids H at the outer surface of the spheroid. To first order in perturbations in δR ,it is given by H = 2 R − [ ∂ θ + (cot θ ) ∂ θ + 2] δRR . (D.2)The different quantities evaluated at r = R + δR read, to first order in perturbations:( σ nn ) r = R + δR = σ stat rr + dσ stat rr dr δR + δσ rr − σ stat rθ ∂ θ δRR ( σ nt ) r = R + δR = σ stat rθ + dσ stat rθ dr δR − ( σ stat θθ − σ stat rr ) ∂ θ δRR + δσ rθ ρ r = R + δR = ρ stat + dρ stat dr δR + δρ ( v r ) r = R + δR = v stat r + dv stat r dr δR + δv r . (D.3)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 symmetriccase of our stationary state. Similarly, the stationary velocity v stat r is also zero at theouter stationary boundary r = R . Appendix E. Additional solutions of the perturbed axisymmetric modedecomposition in the limit of fast diffusion
The components of the pressure and total stress tensor read P ( n ) = (cid:20) (4 n + 6) ηn a n R + ( n − η + 3 nζ n κ R c n (cid:21) ¯ r n σ ( n ) rr = ηR (cid:20)(cid:18) n − n − n a n − n − n κ R c n (cid:19) ¯ r n + 2( n − b n ¯ r n − (cid:21) σ ( n ) rθ = ηR (cid:2) ( − n ( n + 2) a n + n κ R c n ) ¯ r n − n − b n ¯ r n − (cid:3) , (E.1)where the perturbed pressure δP and the perturbed component δσ rr of the stress tensorare decomposed on the basis of Legendre polynomials, and the perturbed component δσ rθ of the stress tensor is decomposed on the basis of Gegenbauer polynomials ‡ . Thethree integration constants a n , b n , and c n are further determined by the three boundaryconditions at the outer surface of the multicellular spheroid, corresponding to thecontinuity of the normal and tangential components of the stress tensor, as well asthat of the metabolite concentration.Solving the boundary-condition equations D.3 for the metabolite field with itsexpression given by equation 19, we get c n = − k Rκ l D R ( n ) . (E.2) ‡ Note that the other non-trivial components of the stress tensor, δσ θθ and δσ φφ , have expansioncoefficients both on the Legendre and Gegenbauer polynomials and are not reported here. iscocapillary Instability in Cellular Spheroids a n and b n : (cid:32) n − n − n n − − n +2 n +1 − n − n (cid:33) · (cid:32) a n b n (cid:33) = (cid:104) − n ( n +1)2 (cid:105) γη R − n − n k (cid:16) Rl D (cid:17) + 2( κρ ext − k ) n +1) k (cid:16) Rl D (cid:17) + κρ ext − k R ( n ) . (E.3)Using further the kinematic condition equation 6 leads to the mode growth rates givenby equation 20.For the modes n = 0 and n = 1, the expressions reported in section 4.3 can furtherbe expressed as functions of either k and ¯ R = R/l D only or κρ ext and ¯ R only, thanksto equation 12. This leads to ω = k (cid:20) ¯ R
3( ¯ R coth ¯ R − − −
19 ¯ R (cid:21) = κρ ext (cid:20)
43 + 3¯ R − (9 + ¯ R ) coth ¯ R R (cid:21) ω = k (cid:20) ¯ R
9( ¯ R coth ¯ R − − −
118 ¯ R (cid:21) = κρ ext (cid:20)
12 + 2¯ R − (6 + ¯ R ) coth ¯ R R (cid:21) . (E.4)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 ω n are given by implicitequations, which cannot be solved analytically. We report here as an example theimplicit equation giving the mode n = 0, which corresponds to the simplest one: ω = k extp − k ¯ R sinh ¯ R (cid:20) cosh ¯ R ¯ R − sinh ¯ R ¯ R (cid:21) , (F.1)where ¯ R = R/l D with l D = (cid:112) D/α as before, and ¯ R = R/l with l = (cid:112) D/ ( α + ω ). Appendix G. Shapes of the lowest-order instability modes
We illustrate in figure G1 the shapes of the lowest-order instability modes.
Appendix H. Mode structure for intermediate values of the diffusioncoefficient of metabolites
We illustrate in figure H1 the mode structure for two values of the diffusion coefficientof metabolites D , intermediate between those of figures 3b,c and 3d. References [1] Saffman P G, Taylor S G and S F R 1958
Proc. R. Soc. Lond. A
Rev. Mod. Phys. iscocapillary Instability in Cellular Spheroids n = 0
Schematic illustration of the lowest-order, axisymmetric, perturbativemodes. The modes n = 0 to n = 5 are represented with an arbitrary amplitude. Ineach schematic representation, the unperturbed spherical shape is represented with adashed line. In the first, n = 0 mode, the axisymmetry is indicated as a rotationalinvariance in φ , and all modes depicted here are invariant under this transformation. - - - ω ( / d a y ) - - - R e [ ω ] ( / d a y ) a)
Mode growth rates as functions of the mode number n for intermediatevalues of the diffusion coefficient of metabolites. Parameters are the same as thoseused in figure 3b,c, using the same colour code, except for the diffusion coefficient ofmetabolites D , which here takes the respective values (a) D = 10 · − cm · s − and(b) D = 20 · − cm · s − , and the metabolite consumption rate α , which is scaledaccordingly to keep the same metabolite penetration lengths l D and associated steadystates.[3] Mullins W W and Sekerka R F 1964 J. Appl. Phys. Rev. Mod. Phys. Proc. R. Soc. Lond. Ser. Math. Phys. Sci.
Phys. Rev. Lett.
Phys. Rev.Lett. iscocapillary Instability in Cellular Spheroids [9] Williamson J J and Salbreux G 2018 Phys. Rev. Lett.
Phys. Rev. Lett.
Nature
Adv. Phys. Nature
J. Phys. Soc. Jpn. Sci. Rep. Soft Matter Phys. Rev. Lett.
Proc. Natl. Acad. Sci.
Nature
Proc. Natl. Acad. Sci.
Curr. Opin. Genet. Dev. Science
Nat. Phys. Phys. Rev. Lett. Phys. Rev. Lett.
J. Theor. Biol. Dev.Cell Curr. Biol. R402–R405 ISSN 0960-9822[29] Sutherland R M 1988
Science
J. Biotechnol.
Biotechnol. Adv. Nat. Biotechnol. Phys. Rev. Lett.
Phys. Rev. Lett.
Jpn. J. Appl. Phys. Nat. Commun. Biophysical Journal Nat. Rev. Mol. Cell Biol. Phys. Rev. Lett.
Proc. Natl. Acad. Sci.
Science
Exp. Cell Res. iscocapillary Instability in Cellular Spheroids [43] Jiang Y, Pjesivac-Grbovic J, Cantrell C and Freyer J P 2005 Biophys. J. Nat. Rev. Cancer J. Theor. Biol. Phys. Rev. E J. Math. Biol. J. Theor. Biol. J. Math. Biol. Phys. Rev. Lett. Math. Comput. Model. Cancer Res. New J. Phys. New J. Phys. Phys. Rev. Lett.
Proc. Natl. Acad. Sci.
Proc. Natl. Acad. Sci.
New J.Phys. Nat. Cell Biol. HFSPJ. HFSP J. New J. Phys. Interface Focus HFSP J. Low Reynolds Number Hydrodynamics: With Special Applicationsto 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.
New J. Phys. Oxygen Transport to Tissue—IV
Advances in Experimental Medicineand 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. Biophys. J. Phys. Biol.
268 ISSN 1478-3975[72] Callan-Jones A C, Joanny J F and Prost J 2008
Phys. Rev. Lett.
Phys. Rev. Lett.
Nature bioRxiv iscocapillary Instability in Cellular Spheroids [76] Aguirre-Ghiso J A 2007 Nat. Rev. Cancer Nat. Med. J. Natl. Cancer Inst. N. Engl. J. Med.
Nature
Cell Cycle Cell Annu. Rev. Med. Cancer Modelling and Simulation (CRC Press) ISBN 978-0-203-49489-9[85] Tracqui P 2009
Rep. Prog. Phys. New J. Phys. Converg. Sci. Phys. Oncol. 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. Nat. Cell Biol. Proc. Natl. Acad. Sci.116