Surface tension controls the onset of gyrification in brain organoids

TL;DR

This study models brain organoids as elastic bodies influenced by surface tension, revealing how cortex growth and cell stiffness affect gyrification and the development of brain folds, with implications for understanding lissencephaly.

Contribution

It introduces a morpho-elasticity-based mathematical model incorporating surface tension to explain gyrification mechanics in brain organoids, highlighting the role of cell stiffness in lissencephaly.

Findings

Surface tension significantly influences organoid mechanics.

Gyrification is triggered by cortex growth and surface tension.

Reduced cell stiffness delays or inhibits folding onset.

Abstract

Understanding the mechanics of brain embryogenesis can provide insights on pathologies related to brain development, such as lissencephaly, a genetic disease which cause a reduction of the number of cerebral sulci. Recent experiments on brain organoids have confirmed that gyrification, i.e. the formation of the folded structures of the brain, is triggered by the inhomo-geneous growth of the peripheral region. However, the rheology of these cellular aggregates and the mechanics of lissencephaly are still matter of debate. In this work, we develop a mathematical model of brain organoids based on the theory of morpho-elasticity. We describe them as non-linear elastic bodies, composed of a disk surrounded by a growing layer called cortex. The external boundary is subjected to a tissue surface tension due the intercellular adhesion forces. We show that the resulting surface energy is relevant at the small length scales of brain organoids and significantly affects the mechanics of cellular aggregates. We perform a linear stability analysis of the radially symmetric configuration and we study the post-buckling behaviour through finite element simulations. We find that the process of gyrification is triggered by the cortex growth and modulated by the competition between two length scales: the radius of the organoid and the capillary length due to surface tension. We show that a solid model can reproduce the results of the in-vitro experiments. Furthermore, we prove that the lack of brain sulci in lissencephaly is caused by a reduction of the cell stiffness: the softening of the organoid strengthens the role of surface tension, delaying or even inhibiting the onset of a mechanical instability at the free boundary.

Figures

d in (Lee et al., [2019](#bib.bib46))). From these data and from Eq. ([3](#S2.E3)), we estimate that the surface tension acting on the spheroid is $\gamma\simeq 0.1\,\mathrm{N}/\mathrm{m}$.

Tables

Table 1. Table 1 : Solutions of the linearised incremental problem at different α R subscript 𝛼 𝑅 \alpha_{R} and different α γ subscript 𝛼 𝛾 \alpha_{\gamma} . The amplitude of the incremental radial displacement has been set equal to 0.05 ​ R o 0.05 subscript 𝑅 o 0.05R_{\text{o}} for the sake of graphical clarity.

	${\alpha }_{\gamma }=0$	${\alpha }_{\gamma }=0.5$	${\alpha }_{\gamma }=1$	${\alpha }_{\gamma }=2$
${\alpha }_R=0.9$
${\alpha }_R=0.95$

Equations

$${\ell }_{\text{c}}=\frac{\gamma }{\mu }.$$

$${\Omega }_{\text{s}}=\left\{\mathbf{X}\in {\mathbb{R}}^3||\mathbf{X}|<R_{\text{o}}\right\},$$

$$\mathrm{div}\mathsf{T}=0.$$

$$\mathsf{T}\mathbf{n}=\gamma \mathcal{K}\mathbf{n},\text{at }R=R_{\text{o}}$$

$$\mathsf{T}=-p\mathsf{I}$$

$$p=\frac{2\gamma }{R_{\text{o}}}$$

$$\gamma =\frac{p{R}_{\text{o}}}{2}.$$

$${\Omega }_0=\left\{\mathbf{X}=[R\cos \Theta ,R\sin \Theta ]\in {\mathbb{R}}^2|0\le R<R_{\text{o}}\text{ and }0\le \Theta <2\pi \right\}$$

$$\mathsf{F}={\mathsf{F}}_{\text{e}}\mathsf{G}$$

$$\det {\mathsf{F}}_{\text{e}}=1.$$

$${\Omega }_{0\text{L}}=\{\mathbf{X}\in \Omega |R<R_{\text{i}}\},{\Omega }_{0\text{C}}=\{\mathbf{X}\in \Omega |R_{\text{i}}<R<R_{\text{o}}\},$$

$$\mathsf{G}=\{\begin{array}{rlrl}& \mathsf{I}& & \text{if }R<R_{\text{i}},\\ & g\mathsf{I}& & \text{if }{R}_{\text{i}}<R<R_{\text{o}},\end{array}$$

$$E_{\mathrm{interface}}={\int }_{\Omega }[{\gamma }_{\mathrm{membrane}}+{\gamma }_{\mathrm{coverslip}}]dv=({\gamma }_{\mathrm{membrane}}+{\gamma }_{\mathrm{coverslip}})|{\Omega }_0|\det \mathsf{F}.$$

$$\mathsf{P}=(\det \mathsf{G}){\mathsf{G}}^{-1}\frac{\partial {\psi }_0({\mathsf{F}}_{\text{e}})}{\partial {\mathsf{F}}_{\text{e}}}-p{\mathsf{F}}^{-1}\mathsf{T}=\frac{1}{\det \mathsf{F}}\mathsf{FP}$$

$$\mathrm{Div}\mathsf{P}=0\text{ in }{\Omega }_{0\text{L}},{\Omega }_{0\text{C}},\text{or}\mathrm{div}\mathsf{T}=0\text{ in }{\Omega }_{\text{L}},{\Omega }_{\text{C}}$$

$$\mathbf{u}(0)=0$$

$$\mathsf{T}\mathbf{n}=\gamma \mathcal{K}\mathbf{n}$$

$$\mathcal{L}(\Theta )=\mathbf{\phi }\left([R_o\cos (\Theta ),R_o\sin (-\Theta )]\right).$$

$${\mathsf{P}}^T\mathbf{N}=(\det \mathsf{F})\gamma \mathcal{K}{\mathsf{F}}^{-T}\mathbf{N}.$$

$$\underset{R\to R_{\text{i}}^{-}}{\lim }{\mathsf{P}}^T\mathbf{N}=\underset{R\to R_{\text{i}}^{+}}{\lim }{\mathsf{P}}^T\mathbf{N}.$$

$$\psi (\mathsf{F})=\det (\mathsf{G}){\psi }_0({\mathsf{F}}_{\text{e}})=(\det \mathsf{G})\frac{\mu }{2}\left(\mathrm{tr}({\mathsf{F}}_{\text{e}}^T{\mathsf{F}}_{\text{e}})-2\right).$$

$$\{\begin{array}{rl}& \mathsf{P}=\mu (\det \mathsf{G}){\mathsf{G}}^{-1}{\mathsf{G}}^{-T}{\mathsf{F}}^{-T}-p{\mathsf{F}}^{-1},\\ & \mathsf{T}=\mu \mathsf{F}{\mathsf{G}}^{-1}{\mathsf{G}}^{-T}{\mathsf{F}}^T-p\mathsf{I}.\end{array}$$

$$\mathbf{\phi }(\mathbf{X})=r(R){\mathbf{e}}_r.$$

$$\mathsf{F}=\mathrm{diag}\left(r^{\prime },\frac{r}{R}\right).$$

$$r(R)=R\text{for }R<R_{\text{i}},$$

$$r^{\prime }r=g^{2}R.$$

$$r(R)=g\sqrt{R^2+\left(\frac{1}{g^2}-1\right)R_{\text{i}}^2}.$$

$$\{\begin{array}{rl}& R=\frac{1}{g}\sqrt{r^2+(g^2-1)R_{\text{i}}^2},\\ & r^{\prime }=g\frac{\sqrt{r^2+(g^2-1)R_{\text{i}}^2}}{r},\end{array}\text{in }{R}_{\text{i}}<R<R_{\text{o}}$$

$$r_{\text{o}}=r(R_{\text{o}})=g\sqrt{R_{\text{o}}^2+\left(\frac{1}{g^2}-1\right)R_{\text{i}}^2}.$$

$$\mathsf{T}{\mathbf{e}}_r=-\frac{\gamma }{r_{\text{o}}}{\mathbf{e}}_r.$$

Full text

In this work, we develop a mathematical model of brain organoids based on the theory of morpho-elasticity. We describe them as non-linear elastic bodies, composed of a disk surrounded by a growing layer called cortex. The external boundary is subjected to a tissue surface tension due the intercellular adhesion forces. We show that the resulting surface energy is relevant at the small length scales of brain organoids and affects the mechanics of cellular aggregates. We perform a linear stability analysis of the radially symmetric configuration and we study the post-buckling behaviour through finite element simulations.

We find that the process of gyrification is triggered by the cortex growth and modulated by the competition between two length scales: the radius of the organoid and the capillary length generated by surface tension. We show that a solid model can reproduce the results of the in-vitro experiments. Furthermore, we prove that the lack of brain sulci in lissencephaly is caused by a reduction of the cell stiffness: the softening of the organoid strengthens the role of surface tension, delaying or even inhibiting the onset of a mechanical instability at the free boundary.

1 Introduction

The formation of folded structures in human and animal brains makes it possible to increase the extension of the cerebral cortex, packing a larger number of neurons in a limited space. The creation of these furrows and ridges called sulci and gyri, respectively, is fundamental for a healthy development of the brain in embryogenesis. The mechanics underlying this morphogenetic phenomenon is not still completely understood.

Recent experiments performed on human brain organoids (Karzbrun et al., 2018), i.e. cell agglomerates cultured in vitro that reproduce the morphogenesis of organs, apparently confirm the hypothesis that sulci are generated by brain cortex buckling triggered by growth (Ronan et al., 2013; Budday et al., 2014; Bayly et al., 2014). In (Karzbrun et al., 2018), the authors observed an increased growth of the cortex with respect to the underlying lumen. In some pathological situations such as lissencephaly, a genetic mutation, the physiological generation of brain sulci is inhibited or even suppressed. This disease, caused by the LIS1 heterozygous ($+/-$) mutation, is correlated to nutritional disorders, alterations in muscle tone, severe psychomotor and mental retardation (Dobyns, 1993). The mathematical description of brain sulci embryogenesis can provide new insights to understand the mechanisms underlying this disease.

A well developed framework to model mechanical instabilities induced by growth is the theory of morpho-elasticity (Goriely, 2017), where living tissues are treated as growing elastic materials. A spatially inhomogeneous growth generates microstructural misfits, leading to a geometrically incompatible relaxed configuration. The restoration of the compatibility requires elastically distorting the body and generating residual stresses (Hoger, 1985; Rodriguez et al., 1994).

Differential growth and residual stresses are involved in the morphogenesis of tissues such as intestinal villi (Balbi and Ciarletta, 2013; Ben Amar and Jia, 2013; Ciarletta et al., 2014) and they enhance the mechanical strength of several biological structures, such as arteries (Chuong and Fung, 1986).

A first model of the experiments on brain organoids (Karzbrun et al., 2018) has been developed by Balbi et al. (2018). The authors model the organoid as a non-linear elastic material, where gyrification is triggered by a remodelling of the cortex and the contraction of the lumen. In their model, the selection of the critical wavelength is dictated by different mechanical properties of the lumen and the cortex. Despite the good agreement with experimental results, Engstrom et al. (2018) noticed that brain organoids exhibit an unconventional behaviour: the cortex is thinner in correspondence of sulci and thicker in correspondence of gyri, in contrast with the morphology predicted by elastic models.

In this paper, we revisit the model proposed by Balbi et al. (2018) to overcome the limitations remarked by Engstrom et al. (2018), proposing a different explanation of lissencephaly.

In cellular aggregates, cohesion among cells is due to adhesion forces induced by adhesion molecules (Turlier and Maître, 2015; Maître et al., 2015). Internal cells are surrounded by other cells, so that the sum of all these forces is zero and each cell is in mechanical equilibrium. Conversely, cells at the boundary of the agglomerate possess a portion of their membrane which is not in contact with other cells: the total adhesion force acting on such cells is non null and it is perpendicular to the free surface of the cellular agglomerate (see Fig. 1 for a graphical representation). These forces generate deformation and the appearance of a boundary layer at the periphery that can be treated as a surface effect called tissue surface tension (Steinberg, 1963).

This phenomenon is reminiscent of the mechanics of surface tension in fluids and soft gels (Style et al., 2017). Since organoids and embryos are characterized by small length scales, surface effects arising from cohesion forces cannot be neglected. The presence of tissue surface tension has been used in fluid models of cellular agglomerates (Foty et al., 1996; Davis et al., 1997; Forgacs et al., 1998) but it is usually overlooked in solid models.

Surface tension has been used within the framework of the theory of non-linear elasticity to regularize the Biot (1963) instability, i.e. the buckling instability of an elastic half-space under compression. Without taking into account surface tension, all the wrinkling modes become unstable at the same compression rate (Ben Amar and Ciarletta, 2010). Non-linear elasticity and particularly models based on elastic bilayers (Li et al., 2012; Cao et al., 2012; Jia et al., 2015; Jin et al., 2019) have been widely employed to model growth phenomena in living tissues, such as the morphogenesis of intestinal villi (Ben Amar and Jia, 2013; Ciarletta et al., 2014) and of brain sulci (Tallinen et al., 2014; Holland et al., 2018). The majority of the literature (see for instance (Jia et al., 2018; Balbi et al., 2019) for solid cylindrical bi-layers) focuses their attention on growth induced instability without considering surface tension. A remarkable exception is provided by the work of Dervaux and Ben Amar (2011), where surface tension is introduced in a cylindrical bi-layer as a regularization of the Biot instability: when the substrate is much more rigid than the coating, the critical wavelength goes to zero. A possible method to regularize this instability is to introduce a small surface tension that introduces a new length scale (the capillary length), leading to a finite critical wavelength (see Ben Amar and Ciarletta (2010)). At small length scales, elasto-capillary forces can deform soft gel beams (Mora et al., 2010, 2013) even inducing mechanical instabilities (Taffetani and Ciarletta, 2015; Xuan and Biggins, 2016). Surface tension can also enhance the resistance to fracture in soft solids (Liu et al., 2014; Hui et al., 2016).

The work is organized as follows: in Section 2 we justify the assumption of a solid model for brain organoids and we compute an estimation of the surface tension acting on a solid cellular aggregate. In Section 3 we develop the elastic model of brain organoids. In Section 4 we perform a linear stability analysis of the radially symmetric configuration and in Section 5 we implement a finite element code to study the post-buckling behaviour. Finally, in Section 6 we discuss the outcomes of our model together with some concluding remarks.

2 Intercellular adhesion generates surface tension in cellular aggregates

At the micro-scale, cellular aggregates are composed by several constituents that, in bulk, have a solid or a fluid-like mechanical nature, like cells, the extracellular matrix, the interstitial fluid. From a macroscopic point of view, these agglomerates can be treated as continuum media but their rheology is still a matter of debate. In fact, cellular aggregates are frequently modelled as fluids (Foty et al., 1994; Manning et al., 2010) that can bear external loading thanks to the tissue surface tension. Another point of view is that cellular aggregates behave as active viscoelastic solids (Kuznetsova et al., 2007; Ambrosi et al., 2017; Karzbrun et al., 2018).

Some important features of biological tissues mechanics cannot be reproduced by fluid-like models. First, there are experimental evidences that cell mitosis and apoptosis (i.e. the cell division and death respectively) are regulated by mechanical stress (Montel et al., 2012b, a). In particular, in (Cheng et al., 2009), the authors report an increased cell duplication in the regions where the compressive stress exerted by the surrounding material on the tumour spheroid is minimum. This spatially inhomogeneous growth can be explained only by using a solid description of the cellular agglomerate: the stress tensor of a fluid at rest corresponds to an hydrostatic pressure which is independent on the spatial position (Ambrosi et al., 2017).

Furthermore, contrarily to fluids, solids can store mechanical stress even in the absence of external loads. These stresses are called residual (Hoger, 1985) and they are created when differential growth in solid matter develops microstructural misfits. These geometrical incompatibilities are restored by elastic distortions of the body, generating mechanical stress (Rodriguez et al., 1994; Goriely, 2017; Riccobelli et al., 2019).

Contrary to the fluid approach, the tissue surface tension of cellular aggregates is frequently neglected in solid models. Nonetheless, solids possess surface tension too (Style et al., 2017): it can play an important role, when the aggregate is very soft or has a small size. In fact, surface tension introduces a new length scale in the problem: let $\mu$ be the shear modulus of the cellular agglomerate and $\gamma$ the surface tension, then the capillary length $\ell_{\text{c}}$ is defined as (Mora et al., 2013; Style et al., 2017)

[TABLE]

Whenever this length scale is of the same order as the characteristic length of the body (e.g. the radius of a spheroid) surface tension cannot be neglected: it can produce a non-negligible deformation (Mora et al., 2013; Style et al., 2013; Mora and Pomeau, 2015) and it can even induce mechanical instabilities (Mora et al., 2010; Taffetani and Ciarletta, 2015; Xuan and Biggins, 2016). Cellular aggregates are very soft and the effect of the surface tension can be highly relevant, as we show in the following.

2.1 Estimation of the surface tension acting on a multicellular spheroid

There are experimental evidences that a tensile skin, having the thickness of a couple of cells, generates an isotropic compression inside multicellular aggregates (Lee et al., 2019). This phenomenon can be explained as a manifestation of tissue surface tension induced by intercellular adhesion: the tensile skin is indeed very thin and such boundary layer can be treated as a surface tension.

Modelling the unloaded multicellular spheroid as a ball occupying the domain

[TABLE]

we assume that the spheroid is composed of an incompressible elastic material. Let $\mathsf{T}$ be the Cauchy stress tensor, the balance of the linear and angular momentum reads

[TABLE]

If the spheroid is subjected to a tissue surface tension $\gamma$, the boundary condition reads (Style et al., 2017)

[TABLE]

where $R$ denotes the radial position, $\mathcal{K}$ is twice the mean curvature, and $\bm{n}$ is the external normal.

We now show that the undeformed reference configuration is in mechanical equilibrium. Experimental observations (Lee et al., 2019) lead us to state that residual stresses are absent. The Cauchy stress is then given by (Ogden, 1997)

[TABLE]

where $\mathsf{I}$ is the identity and $p$ is the pressure field that enforces the incompressibility constraint. A constant pressure field $p$ satisfies the equilibrium equation (1). By imposing Eq. (2), we get

[TABLE]

which is nothing but the Young-Laplace equation. We invert the previous equation with respect to $\gamma$, obtaining

[TABLE]

From the work of Lee et al. (2019), we get that the typical radius of a spheroid is $\sim 400\,\mu\mathrm{m}$ and the internal pressure $p$ is about $500\,\text{Pa}$ (Fig. 5d in (Lee et al., 2019)). From these data and from Eq. (3), we estimate that the surface tension acting on the spheroid is $\gamma\simeq 0.1\,\mathrm{N}/\mathrm{m}$.

Remark 1**.**

In this Section, we have estimated the tissue surface tension acting on a spheroid. Unfortunately, up to our knowledge, there are not similar experimental data that can be used to perform the same computation for brain organoids. There are some measurements of the surface tension acting on different embryonic tissues when they are treated as fluids (Schötz et al., 2008). In their work, all the measured values have the same order of magnitude. From these data, we can assume that $\gamma\simeq 0.1\,\mathrm{N}/\mathrm{m}$ is a qualitative estimate of tissue surface tension acting on a brain organoid.

3 Elastic model of brain organoids

In this Section, we illustrate a model of brain organoids, described as growing hyperelastic bodies subjected to surface tension.

3.1 Kinematics

We denote by $\bm{X}$ the material position coordinate. Since brain organoids are cell aggregates confined in a narrow space by a coverslip and a membrane (Karzbrun et al., 2018), we model them as two dimensional objects. Let

[TABLE]

be the reference configuration of the organoid. We indicate with $\bm{\varphi}:\Omega_{0}\rightarrow\mathbb{R}^{2}$ the deformation field, so that the actual configuration of the body $\Omega$ is given by $\bm{\varphi}(\Omega_{0})$. Let $\bm{x}=\bm{\varphi}(\bm{X})$ be the actual position of the point $\bm{X}$ and the displacement vector is defined as $\bm{u}(\bm{X})=\bm{\varphi}(\bm{X})-\bm{X}$. Let $\mathsf{F}$ be the deformation gradient, i.e. $\mathsf{F}=\operatorname{Grad}\bm{\varphi}$. We exploit a multiplicative decomposition of the deformation gradient (known as Kröner-Lee decomposition (Kröner, 1959; Lee, 1969)) to model the growth of the organoid, so that

[TABLE]

where the growth tensor $\mathsf{G}$ accounts for the local inelastic distortion due to the body growth, while $\mathsf{F}_{\text{e}}$ describes the elastic distortion necessary to maintain the geometrical compatibility of the body and to balance the external and internal forces (Rodriguez et al., 1994). As commonly done for biological tissues, it is reasonable to model organoids as incompressible media, namely we enforce that

[TABLE]

In (Karzbrun et al., 2018), the authors identified two distinct regions in brain organoids: an internal lumen and an external ring, called cortex, the latter being characterized by a faster growth. Indicating with $R_{\text{i}}$ the radius of the lumen, we denote these two regions by $\Omega_{0\text{L}}$ and $\Omega_{0\text{R}}$:

[TABLE]

and their images through $\bm{\varphi}$ are denoted by $\Omega_{\text{L}}$ and $\Omega_{\text{R}}$ respectively.

We assume that the growth tensor $\mathsf{G}$ takes diagonal form

[TABLE]

where the scalar quantity $g$ is the growth rate of the cortex with respect to the lumen. Tissue surface tension acts at both the interfaces organoid-coverslip and organoid-membrane. Let us call $\gamma_{\mathrm{membrane}}$ and $\gamma_{\mathrm{coverslip}}$ the tension acting on the organoid at the interface with these two surfaces. Their contribution to the total energy is given by

[TABLE]

Since the deformation is isochoric, $\det\mathsf{F}=\det\mathsf{G}$. The contribution of the surface tension acting on both the membrane and the coverslip, given by $E_{\mathrm{interface}}$, does not change when the organoid undergoes an elastic deformation and it does not influence either the stability or the morphology of the organoid. This agrees with the observations of Karzbrun et al. (2018): the authors tried different coverslips and membranes, each one characterized by a different adhesion to the organoid and they did not notice remarkable differences (see Fig. S9 of the supplementary material of Karzbrun et al. (2018)). From now on, since $E_{\text{interface}}$ does not influence the deformation of the body, it is not explicitly included in the total energy.

We now introduce some mechanical constitutive assumptions.

3.2 Mechanical constitutive assumptions and force balance equations

We assume that the organoids are composed of a homogenous hyperelastic material, having strain energy density $\psi$. The first Piola–Kirchhoff stress $\mathsf{P}$ and the Cauchy stress tensors $\mathsf{T}$ are then given by

[TABLE]

where $p$ is the Lagrange multiplier enforcing the incompressibility constraint (4). The balance of the linear and angular momentum reads

[TABLE]

in the material and actual reference frame, respectively.

We assume that the center of the organoid is fixed, i.e.

[TABLE]

while a constant surface tension $\gamma$ acts at the external boundary of the organoid, so that (Style et al., 2017)

[TABLE]

where $\bm{n}$ is the outer normal in spatial coordinates and $\mathcal{K}$ is the oriented curvature of the boundary curve $\bm{\mathcal{L}}$ parametrized clockwise, i.e.

[TABLE]

The Lagrangian form of the boundary condition (8) is obtained performing a pull-back

[TABLE]

Finally, we enforce the continuity of the stress at the interface $R=R_{\text{i}}$, so that

[TABLE]

We can assume that brain organoids behave as isotropic media, since they are composed of neural progenitors and not fully developed neurons (Karzbrun et al., 2018). Moreover, in the fully developed brain, both the gray and the white matter are essentially isotropic (Budday et al., 2017a, b). As it concerns the internal structure, brain organoids are complex cell agglomerates composed of different biological materials, such as the extra-cellular matrix, cells, interstitial fluid. According to the observations reported in Karzbrun et al. (2018), the extra cellular matrix and the cellular cytoskeleton account for the elasticity of the organoid. We have assumed that the two elastic phases in the organoids generate a homogeneous hyperelastic material. There are several works which show a very good fit of the experimental data by assuming a Mooney–Rivlin model (especially when the stretches are moderate (Destrade et al., 2015; Balbi et al., 2019)). Such an energy reduces to the neo-Hookean one whenever the deformation is planar. Hence, the strain energy density is given by

[TABLE]

The first Piola–Kirchhoff and Cauchy stress tensors read, respectively,

[TABLE]

Summing up, Eqs. (4) and (6), together with the kinematic constraint Eq. (7) and the boundary condition Eq. (10) define the non linear elastic problem. In the next Section we look for a radially symmetric solution.

3.3 Equilibrium radially-symmetric solution

Let $(r,\,\theta)$ be the actual radial and polar coordinates of a point. Let $(\bm{E}_{R},\,\bm{E}_{\Theta})$ and $(\bm{e}_{r},\,\bm{e}_{\theta})$ be the local vector basis in polar coordinates in the Lagrangian and Eulerian reference frame, respectively. We look for a radially-symmetric solution of the form

[TABLE]

The deformation gradient expressed in polar coordinates reads

[TABLE]

It is immediate to notice that

[TABLE]

where $r_{\text{i}}=r(R_{\text{i}})=R_{\text{i}}$. In the cortex, from the incompressibility constraint given by Eq. (4), we get

[TABLE]

Performing an integration and imposing that $R_{\text{i}}=r_{\text{i}}$, we get

[TABLE]

It remains to determine the pressure field $p$. First, we notice that, inverting and differentiating (15), we obtain

[TABLE]

respectively. The curvature of the boundary line is $-r_{\text{o}}^{-1}$, where

[TABLE]

The boundary condition Eq. (8) reads

[TABLE]

Since the deformation depends only on the radial position $r$, the balance of the linear and angular momentum in polar coordinates reads

[TABLE]

where $T_{ij}$ are the components of the Cauchy stress tensor $\mathsf{T}$ in polar coordinates. From Eqs. (13), (14) and (16), the Cauchy stress in the cortex reads

[TABLE]

We can integrate Eq. (17) from $r$ to $r_{\text{o}}$, obtaining

[TABLE]

We can find the pressure field in the cortex (i.e. for $r_{\text{i}}<r<r_{\text{o}}$) by plugging Eq. (18) into Eq. (19), obtaining

[TABLE]

Finally, we impose the continuity of the stress Eq. (11) at $r=r_{\text{i}}=R_{\text{i}}$ to get the pressure for $r<r_{\text{i}}$. Since the lumen remains undeformed, the Cauchy stress reads

[TABLE]

Using Eq. (18) and Eq. (20), we can write the solution

[TABLE]

In the next section a linear stability analysis of the solution given by Eqs. (21)-(22) is performed, but first we remark some important aspects regarding the boundary conditions.

Remark 2**.**

In this model, we have assumed that the only force acting on the boundary of the organoid is due to tissue surface tension. We remark that, in the experiments reported in (Karzbrun et al., 2018), organoids are embedded in Matrigel. Karzbrun and co-authors tried different concentration of the gel, without registering remarkable changes to the morphology of the wrinkling (see Fig. S9 (Karzbrun et al., 2018)), so that we are lead to conjecture that Matrigel has a negligible effect on the onset of the instability.

Matrigel is characterized by an elastic modulus of the same order as the organoid one ($500$ Pa Soofi et al. (2009)). In each experiment, eleven organoids are placed in a dish having a diameter of 6 cm, so that the distance among organoids is of the order of centimeters. The distance among organoids is much larger than the displacement induced by the organoid growth, so that the response of the Matrigel is linear and can be modelled as the action of linear springs acting at the boundary, as done for similar motivations in (Riccobelli and Ciarletta, 2018b). We have verified through this approach that the role of Matrigel is negligible. For the readers interested in the modeling of the Matrigel embedment, we expose in Appendix B some details on the linear stability analysis in the presence of linear springs at the boundary.

4 Linear stability analysis

4.1 Incremental equations

We apply the theory of incremental deformations superposed on finite strains (Ogden, 1997) to investigate the stability of the radially symmetric solution. Let $\delta\bm{u}$ be the incremental displacement field and let $\mathsf{\Gamma}=\operatorname{grad}\delta\bm{u}$. We introduce the push-forward of the incremental Piola-Kirchhoff stress in the axisymmetric deformed configuration, given by

[TABLE]

where $\mathcal{A}_{0}$ is the fourth order tensor of instantaneous elastic moduli, $\delta p$ is the increment of the Lagrangian multiplier that imposes the incompressibility constraint. The two dots operator $(:)$ denotes the double contraction of the indices

[TABLE]

where the convention of summation over repeated indices is used. The components of the tensor $\mathcal{A}_{0}$ for a neo-Hookean material are given by

[TABLE]

where $\mathsf{B}_{\text{e}}=\mathsf{F}_{\text{e}}\mathsf{F}_{\text{e}}^{T}$ and $\delta_{ik}$ is the Kronecker delta. The incremental equilibrium equation and the linearised form of the incompressibility constraint read respectively

[TABLE]

The linearised form of the kinematic constraint (7) and of the boundary condition (8) complement the incremental equations

[TABLE]

where $\delta\mathcal{K}$ is the increment of the curvature. Finally, we enforce the continuity of the incremental displacement of the stress at the interface

[TABLE]

In the following, we rewrite the incremental problem given by the Eqs. (LABEL:eq:deltaP)-(26) into a more convenient form using the Stroh (1962) formulation.

4.2 Stroh formulation

We rewrite the incremental problem in non-dimensional form adopting as the characteristic length scale and shear modulus $R_{\text{o}}$ and $\mu$ respectively. The behaviour of the problem is governed by the non-dimensional parameters

[TABLE]

in addition to the growth parameter $g$.

For the sake of brevity, we introduce the multi-index ${\text{W}}=\{{\text{L}},{\text{C}}\}$. The quantities with subscript L are computed in the lumen, while the ones in the cortex have the subscript C.

We denote with $u_{\text{W}}$ and $v_{\text{W}}$ the components of $\delta\bm{u}_{\text{W}}$ while $\delta P^{\text{W}}_{rr}$ and $\delta P^{\text{W}}_{r\theta}$ are the components of the incremental stress along the radial direction $\bm{e}_{r}$. We can reduce the system of partial differential equations (LABEL:eq:deltaP) into a system of ordinary differential equations by assuming the following ansatz for the incremental displacement, pressure and stress:

[TABLE]

where $m\in\{n\in\mathbb{N}\;|\;n\geq 2\}$ is the circumferential wavenumber. By substituting Eqs. (29)-(31) and Eq. (33) into Eq. (LABEL:eq:deltaP), we obtain the following expression for $Q^{\text{L}}$ and $Q^{\text{C}}$

[TABLE]

where $p_{\text{L}}$, $p_{\text{C}}$ are defined in Eq. (22).

By the choices (29)-(32) and using a well established procedure, the incremental problem can be rewritten in the Stroh (1962) form

[TABLE]

where $\bm{\eta}_{\text{W}}$ is the displacement-traction vector defined as

[TABLE]

The matrix $\mathsf{N}_{\text{W}}\in\mathbb{R}^{4\times 4}$ is the Stroh matrix and it has the following sub-block form

[TABLE]

For the lumen ($r<r_{\text{i}}$), the sub-blocks read:

[TABLE]

In the cortex ($r_{\text{i}}<r<r_{\text{o}}$), the sub-blocks are given by

[TABLE]

where

[TABLE]

We remark that the coefficient of the Stroh matrix are constant in the lumen (see Eq. (35)). This allows us to solve analytically the incremental problem for $r<r_{\text{i}}$.

4.3 Incremental solution for the lumen

We follow the procedure proposed in (Dervaux and Ben Amar, 2011) and in (Balbi et al., 2018). Since Eq. (34) with $\text{W}=\text{L}$ is a system of ODEs with constant coefficients, its solution can be rewritten in terms of eigenvalues and eigenvectors of $\mathsf{N}_{\text{L}}$. The eigenvalues of $\mathsf{N}_{\text{L}}$ are $\lambda_{1}=m-1$, $\lambda_{2}=m+1$, $\lambda_{3}=-m+1$ and $\lambda_{4}=-m-1$. The general integral of Eq. (34) is given by

[TABLE]

where $\bm{w}_{i}$ are the eigenvectors of $\mathsf{N}_{\text{L}}$ associated with the eigenvalues $\lambda_{i}$, $i=1,\dots,\,4$. Since the incremental solution must satisfy the kinematic constraint Eq. (25), we immediately get that $c_{3}=c_{4}=0$ while

[TABLE]

The two constants $c_{1}$ and $c_{2}$ will be fixed by imposing the continuity of the displacement and of the stress at $r=r_{\text{i}}$ (i.e. by enforcing Eq. (27)).

4.4 Numerical procedure for the solution in the cortex

The incremental problem in the cortex cannot be solved analytically since the coefficient of the Stroh matrix $\mathsf{N}_{\text{C}}$ given by Eq. (36) are not constant. To overcome this difficulty, we implement a numerical code based on the impedance matrix method (Biryukov, 1985; Biryukov et al., 1995).

We first introduce the conditional impedance matrix, defined as

[TABLE]

Such a matrix is called conditional since its expression depends on an auxiliary condition at $r=r_{\text{o}}$, in this case the boundary condition Eq. (26). In the following paragraphs we expose a procedure to construct the matrix $\mathsf{Z}_{\text{C}}(r,\,r_{\text{o}})$. The incremental boundary condition Eq. (26) reads

[TABLE]

where the incremental curvature $\delta\mathcal{K}$ is given by (see the Appendix A for the details of the computation):

[TABLE]

From the ansatz of variable separation in Eqs.(29)-(30) and using the incremental form of the incompressibility constraint Eq. (LABEL:eq:deltaP), the boundary condition Eq. (40) is equivalent to

[TABLE]

We can now define the auxiliary impedance matrix (Norris and Shuvalov, 2010) as

[TABLE]

so that the boundary condition Eq. (41) is equivalent to the equation

[TABLE]

We introduce the matricant

[TABLE]

called conditional matrix, defined as the solution of the problem

[TABLE]

From Eq. (44), the Stroh form of the incremental problem given by Eq. (34) and Eq. (43), we get

[TABLE]

so that

[TABLE]

From Eq. (45), the conditional impedance matrix is given by

[TABLE]

From now on, we omit the dependence of $\mathsf{Z}_{\text{C}}$ wherever convenient for sake of simplicity. By using Eq. (39), we can rewrite the Stroh problem given by Eq. (34) into a differential Riccati equation. Indeed, from Eq. (34), we get

[TABLE]

Substituting Eq. (46) into Eq. (47) we get the following differential Riccati equation

[TABLE]

We integrate Eq. (48) from $r_{\text{o}}$ to $r_{\text{i}}$, using as initial condition the auxiliary impedance matrix defined in Eq. (42), i.e.

[TABLE]

To construct a bifurcation criterion, we follow Balbi et al. (2018): from the continuity of the displacement-traction vector $\bm{\eta}_{\rm C}(r_{\text{i}})=\bm{\eta}_{\rm L}(r_{\text{i}})$ Eq. (27) we get

[TABLE]

so that non-null solutions of the incremental problem exist if and only if

[TABLE]

where

[TABLE]

with $\bm{w}_{i}$ defined in Eq. (38). For a fixed value of the control parameter $g$ we integrate the Riccati Eq. (48) from $r=r_{\text{o}}$ up to $r=r_{\text{i}}$ making use of the software Mathematica 11.3 (Wolfram Research, Champaign, IL, USA). We iteratively increase the control parameter $g$ until the bifurcation criterion Eq. (49) is satisfied.

4.5 Discussion of the results

First, we need to identify an interval of interest for the adimensional parameter $\alpha_{\gamma}$. We have estimated the surface tension of cellular aggregates in Section 2.1. From the stress profiles reported by Lee et al. (2019), we obtain a surface tension of the order of $10^{-1}\mathrm{N/m}$. According to Karzbrun et al. (2018), the shear modulus of the wild-type brain organoid is $\mu\simeq 900\,\mathrm{Pa}$ (Young modulus $E\simeq 2.7\,\mathrm{kPa}$) while $\mu\simeq 333\,\mathrm{Pa}$ (Young modulus $E\simeq 1\,\mathrm{kPa}$) for the unhealthy ones, afflicted by lissencephaly (i.e. when the mutation LIS1 $+/-$ is present). Provided that the typical radius of the organoid is about $R_{\text{o}}=400\,\mu\mathrm{m}$ (Karzbrun et al., 2018), $\alpha_{\gamma}$ ranges between $0.25$ for the wild-type organoids and $0.75$ for the ones affected by lissencephaly.

Let us now discuss the results of the linear stability analysis. For fixed values of the adimensional parameters $\alpha_{R}$ and $\alpha_{\gamma}$, we denote by $g_{m}$ the first value of $g$ such that the bifurcation criterion (49) is satisfied for the wavenumber $m$. We define the critical threshold $g_{\text{cr}}$ as the minimum $g_{m}$ for $m\geq 2$ and the critical mode $m_{\text{cr}}$ as the wavenumber corresponding to $g_{\text{cr}}$.

In Fig. 2, we plot the critical values $g_{\text{c}r}$ and $m_{\text{c}r}$ versus $\alpha_{\gamma}$ for two different values of $\alpha_{R}$, i.e. in the first one $\alpha_{R}=0.9$, while in the other one $\alpha_{R}=0.95$.

We observe that, for relatively small values of $\alpha_{\gamma}$, the marginal stability threshold $g_{\text{c}r}$ increases monotonously as $\alpha_{\gamma}$ increases, while the critical wavenumber $m_{\text{cr}}$ decreases. There is a change in the behaviour of the instability when the parameter $\alpha_{\gamma}$ is sufficiently large: the critical wavenumber $m_{\text{cr}}$ increases suddenly and the marginal stability threshold $g_{cr}$ remains nearly constant about $g_{\text{cr}}\simeq 2.5$.

The threshold at which this transition occurs strongly depends on $\alpha_{R}$, as shown in Fig. 2. Indeed, when $\alpha_{R}=0.9$, $g_{\text{c}r}$ increases from $1.745$ to $2.481$ with $\alpha_{\gamma}\in(0,1.34)$, while when $\alpha_{R}=0.95$, $g_{\text{c}r}$ increases in approximately the same range as for $\alpha_{R}=0.9$, i.e. $(1.671,2.510)$, but $\alpha_{\gamma}$ varies in a smaller interval, i.e. $\alpha_{\gamma}\in(0,0.86)$.

To study the morphology of the critical mode, we have integrated Eq. (46) to compute the incremental displacement field, as described in (Destrade et al., 2009). We depict in Tab. 1 a morphological diagram where we show the solution of the incremental problem for different values of $\alpha_{\gamma}$ and $\alpha_{R}$. For small values of $\alpha_{\gamma}$, we observe that the instability mainly releases elastic energy at the free boundary, displaying a wrinkling pattern: as we increase $\alpha_{\gamma}$ the wavenumber decreases and the critical mode displays a more rounded boundary. Furthermore, for larger values of $\alpha_{\gamma}$ there is a drastic change in the features of the instability: the morphological transition localizes at the interface between the cortex and the lumen with a high critical wavenumber. Let us call $\alpha_{\gamma}^{\rm tr}$ the smallest value at which, fixing $\alpha_{R}$, the instability localizes at the interface between the cortex and the lumen. We plot in Figure 3 $\alpha_{\gamma}^{\rm tr}$ versus $\alpha_{R}$. Interestingly, the transition from a surface to an interfacial instability can take place only if surface tension is relatively large. It is to be remarked that, whenever a surface instability takes place, the stretch induced by growth is quite large at the boundary between cortex and lumen. It is possible that, instead of wrinkles with short wavelength, interfacial creases appear. However, the nucleation of sharp creases is the result of a nonlinear instability, not detectable by a linear analysis (Ciarletta and Truskinovsky, 2019).

We also investigate the influence of $\alpha_{R}$ on the instability by fixing $\alpha_{\gamma}$ (see Fig. 4). In the absence of surface tension (i.e. $\alpha_{\gamma}=0$), we observe that $g_{\rm cr}$ decreases monotonously as $\alpha_{R}$ increases (see Fig. 4(a)). The behaviour is the opposite in the presence of surface tension, where the marginal stability threshold $g_{\text{cr}}$ monotonously increases with $\alpha_{R}$. For the range of parameters in which the critical mode displays an interfacial instability, we observe that $g_{\text{cr}}$ increases linearly with $\alpha_{R}$.

As regards the critical wavenumber $m_{\rm cr}$, we can see that it increases in all the cases, both in the presence and in the absence of surface tension.

We observe that our model captures the main features of organoid development. First, for small values of $\alpha_{\gamma}$, a morphological transition takes place at the free boundary, giving rise to a wrinkling pattern: as $\alpha_{\gamma}$ increases, we observe a decrease of the critical wavenumber and a higher marginal stability threshold. This is in agreement with the experiments of Karzbrun et al. (2018): if the cells have the LIS1 $+/-$ mutation, the authors measured that the elastic modulus of the cells is $2.7$ times lower than the one of healthy cells. In our model this reduction of the stiffness is equivalent to an increase of $\alpha_{\gamma}$. They have also reported a reduction of the number of folds in organoids affected by lissencephaly. As one can observe from the plots of Figs. 2-4 and from the morphological diagram of Table 1, as we increase $\alpha_{\gamma}$ the number of wrinkles decreases and the critical threshold increases, in accordance with Karzbrun et al. (2018)’s experiments (these results are shown in Karzbrun et al. (2018), Figure 5). Furthermore, for large values of $\alpha_{\gamma}$, wrinkles at the free surface are completely absent, as happens in the most serious case of lissencephaly.

Compared to the model proposed by Balbi et al. (2018), in which the authors do not take into account the tissue surface tension, our theoretical description presents some advantages. In fact, we do not introduce different shear moduli for the cortex and the lumen to modulate the critical wavenumber and the critical growth threshold. This choice is motivated by the experimental results of Karzbrun et al. (2018): the authors reported a unique value of elastic modulus for the organoid and they did not experimentally measure a change in the shear moduli between the cortex and the lumen.

In our model the selection of the critical wavenumber is controlled by the competition between surface capillary energy and bulk elasticity. Furthermore, we are able to justify the complete absence of surface wrinkles in the most severe cases of lissencephaly, which corresponds to the case of large $\alpha_{\gamma}$: large values of $\alpha_{\gamma}$ corresponds to both a high value of surface tension and a very soft material, the shear modulus is decreasing.

Notwithstanding the good agreement with the experimental results of our model, it is to be reported that models based on solid mechanics of brain organoids have been recently criticized by Engstrom et al. (2018). The authors observe that the folds of the cortex display an “antiwrinkling” behaviour: the cortex is thicker in correspondence of furrows and thinner at the ridges of wrinkles. The authors claim that solid models do not show this feature and, thus, they are inadequate to model multicellular aggregates. In the next section we implement a numerical code to approximate the fully non-linear problem and we show that the “antiwrinkling” phenomenon is provoked by tissue surface tension.

5 Post-buckling analysis

5.1 Description of the numerical method

In this section we show the results of the numerical approximation of the non-linear problem given by Eqs. (6)-(7)-(10) to investigate the post-buckling behaviour of the organoid. We use as computational domain the rectangle obtained through the conformal mapping corresponding to the polar coordinate transformation, as done in (Riccobelli and Ciarletta, 2018a): let

[TABLE]

Given $\bm{X}_{\text{n}}\in\Omega_{\text{n}}$, the components represent the referential radial coordinate normalized with respect to the external radius, and the referential polar angle, respectively:

[TABLE]

as represented in Fig. 5. The function

[TABLE]

maps the computational domain to a half circle, which represent half of the reference configuration. The full domain can be obtained thanks to the axial symmetry of the problem.

We discretise the computational domain through the software GMSH (Geuzaine and Remacle, 2009). We use a triangular grid, with a progressive refinement of the elements from $X_{1}^{\text{n}}=0$ up to $X_{1}^{\text{n}}=\alpha_{R}$. In the cortex we instead use a structured mesh (i.e. for $\alpha_{R}<X_{1}^{\text{n}}<1$), see Fig. 6.

We implement a numerical code based on the mixed finite element method to enforce the incompressibility constraint (4). We discretise the displacement field $\bm{u}$ using continuous, piecewise quadratic functions, while we approximate the pressure through piecewise constant functions. The corresponding mixed finite element is the $\bm{P}_{2}-P_{0}$ element, which is numerically stable for the incompressible hyperelastic problem (Boffi et al., 2013). We use an index $h$ when we refer to the discretised counterpart of the mathematical quantities. We adimensionalise the system of equation with respect to $\mu$ and $R_{\text{o}}$ as we did at the continuum level. We set the following discrete boundary conditions

[TABLE]

where $(\bm{e}_{1},\,\bm{e}_{2})$ represents the canonical vector basis.

We solve the discretised form of the equilibrium equation (6) in the Lagrangian form using a Newton method. The control parameter $g$ is incremented of $\delta g$ when the Newton method converges, the numerical solution is used as initial guess for the following Newton cycle. The increment $\delta g$ is automatically reduced near the theoretical marginal stability threshold and when the Newton method does not converge. The numerical simulation is stopped when $\delta g<10^{-6}$. To trigger the mechanical instability, a small perturbation of an amplitude of $\sim 10^{-5}$, having the shape of the critical mode computed in Section 4, is applied at the free boundary of the mesh. We have numerically verified that the wavelength of the buckled pattern is not sensitive to the applied imperfection and the only effect is a slight anticipation of the instability threshold.

The numerical algorithm is implemented in Python through the open-source computing platform FEniCS (version 2018.1) (Logg et al., 2012). The computation of the weak form and of the Jacobian necessary to solve each step of the Newton method are computed from the total energy through the library UFL (Alnæs et al., 2014). Surface tension is introduced in the numerical algorithm following the implementation proposed in (Mora et al., 2013). We use PETSc (Balay et al., 2018) as linear algebra back-end and MUMPS (Amestoy et al., 2000) as linear solver.

5.2 Results of the finite element simulations

In this section, we discuss the results of the numerical simulations for $\alpha_{R}=0.9$.

In Fig. 7, we plot the buckled configuration of the organoid. As predicted by the critical modes of the linear stability analysis plotted in Table 1, in presence of surface tension (Fig. 7(b)) the free boundary is smoother and rounded.

We define $\Delta r$ as the amplitude of the pattern at the free surface

[TABLE]

where $r$ denotes the actual radial position of the point with polar coordinates $(R,\,\Theta)$.

In Fig. 8, we show how the amplitude of the pattern $\Delta r$ evolves with respect to the growth rate $g$. We observe that there is an excellent agreement with the marginal stability thresholds computed in the previous section, proving the consistence of the numerical code with respect to the theoretical predictions.

Both the bifurcation diagrams exhibit a continuous transition from the unbuckled to the buckled configuration, displaying the typical behaviour of a supercritical pitchfork bifurcation. Let

[TABLE]

be the total mechanical energy of the organoid. We compute the ratio of the energy $E_{\text{th}}$ of the base solution given by Eqs. (21)-(22) and the energy $E_{\text{num}}$ arising from the numerical simulations.

In Fig. 9 we plot the ratio $E_{\text{num}}/E_{\text{th}}$ versus the control parameter $g$ for both $\alpha_{\gamma}=0$ and $\alpha_{\gamma}=0.5$. In both cases, the buckled configuration exhibits a total lower mechanical energy with respect to the unbuckled state. Furthermore we observe that the energy lowers continuously, confirming that the bifurcation is supercritical.

Finally, we compute the thickness of the cortex at the ridges and at the furrows of the buckled configuration (which correspond to the gyri and the sulci of the fully developed brain respectively). We observe that the thickness of the ridges is higher than the one of the furrows for $\alpha_{\gamma}=0$ (see Fig. 10(a)), while the behaviour is the opposite in the case of $\alpha_{\gamma}=0.5$ (see Fig. 10(b)). The latter case is typical of brain organoids, as observed in (Engstrom et al., 2018).

These numerical results suggest that the “antiwrinkling” behaviour, as named by Engstrom et al. (2018), is due to the presence of surface tension, which is highly relevant due to the small radius of a brain organoid. If we consider a fully developed brain, its typical size is of order of decimeters. In this case, keeping $\gamma$ and $\mu$ fixed, $\alpha_{\gamma}$ is reduced by five order of magnitudes with respect to the case of the organoid and the contribution of surface tension to the total energy becomes negligible. Experimental results show that gyri are thicker than sulci (Holland et al., 2018), differently from what happens at the small length scales of the brain organoid. This is in agreement with the outcomes of our model in the case $\alpha_{\gamma}=0$, confirming that the “antiwrinkling” behaviour is caused by the competition of bulk elastic energy and the surface tension.

However, our numerical algorithm has some limitations. The Newton method does not converge anymore slightly past the marginal stability threshold. In fact, the bifurcation diagrams of Fig. 8 show that the amplitude of the pattern increases very rapidly beyond the marginal stability threshold. As one can observe in Fig. 7(a), the deformation tends to localize near the furrows, highlighting a possible subcritical transition to a folded state, which would be in agreement with the experimental observations of Karzbrun et al. (2018). The numerical approximation of a layer growing on a substrate with similar mechanical properties is particularly tricky, even numerically (Fu and Ciarletta, 2015). Future efforts will be devoted to the improvement of the numerical scheme, implementing an arclength continuation method to study the wrinkle-to-fold transition in brain organoids.

6 Discussion and concluding remarks

In this work, we have developed a model of brain organoids to describe the formation of cerebral sulci and to investigate the influence of surface tension on such a morphogenetic process. In Section 2 we have computed that the tissue surface tension acting on a solid multicellular spheroid is $\gamma=10^{-1}\,\mathrm{N}/\mathrm{m}$, using the data reported in Lee et al. (2019). We exploit this measure as a qualitative estimate of tissue surface tension of a generic cellular aggregate. In fact, measures performed on different embryonic tissues modelled as fluids all have the same order of magnitude (Schötz et al., 2008). Thus, we do not expect that the tissue surface tension acting on brain organoids dramatically differs from the one acting on multicellular spheroids.

Then, we have built in Section 3 an elastic model of brain organoids. They are described as disks surrounded by a growing rim and subjected to a surface tension generated by intercellular adhesion forces (Manning et al., 2010). We have assumed that the two regions (disk and outer rim) are composed of the same incompressible neo-Hookean material. The system is governed by the dimensionless parameters $g$, i.e. the growth rate of the cortex with respect to the lumen, $\alpha_{R}$ and $\alpha_{\gamma}$, which are the lumen radius and the capillary length, normalized with respect to the initial radius of the organoid, respectively.

We have computed a radially symmetric solution and we have studied its linear stability in Section 4 using the theory of incremental deformations. We have rewritten the linear stability analysis into an optimal Hamiltonian system using the Stroh formulation (Fu, 2007). The impedance matrix method is adapted to take into account the boundary contribution of surface tension.

The outcomes are discussed in Section 4.5. The introduction of a tissue surface tension of the same order as the estimate computed in Section 2.1 strongly influence the stability of the base solution and alters the critical wavelength. This result suggests that the interplay between elasticity and tissue surface tension plays a crucial role in controlling pattern selection in embryo morphogenesis. The predicted critical wavelength is always finite (see Fig. 4(a)). In our model the softening of cells, due to lissencephaly, corresponds to an increase of the parameter $\alpha_{\gamma}$, strengthening the role of tissue surface tension and leading to a lower critical wavenumber. Experimental observations of Karzbrun et al. (2018) report a similar behaviour: when the LIS1 $+/-$ mutation is present, the shear modulus of the organoid is reduced and the wavenumber decreases compared with the healthy organoids. This suggest that the reduction of brain sulci due to lissencephaly is due to the competition between elastic and surface energies in early embryogenesis. Moreover, we have observed that, for larger $\alpha_{\gamma}$, a transition from a surface to an interfacial instability occurs: buckling localizes at the interface between the cortex and the lumen. In this case the cortex remains smooth as one can observe in the most severe cases of lissencephaly. The results are reported in Figs. 2-3-4 and in Tab. 1.

Finally, in Section 5, we have implemented a finite element code to approximate the fully non-linear problem. The algorithm is based on a mixed variational formulation and the Newton method. The outcomes of the numerical simulations are reported in Fig. 7-10. These results show that tissue surface tension rounds the external boundary. Both in the presence and in the absence of surface energy, the bifurcation is supercritical, displaying a continuous transition from the unbuckled to the buckled state. Contrarily to the bilayer model without surface tension, our model reproduces “antiwrinkling” behaviour of the cortex (namely the thickness of the outer layer is larger in the furrows). This suggests that this unconventional variation of cortex thickness is due to tissue surface energy. This strengthens the importance of considering surface tension in the modelling of cellular aggregates.

Summing up, our model suggests a possible purely mechanical explanation of a number of open questions:

• •

brain sulci in organoids are generated by a mechanical instability triggered by cortex growth;

• •

the “anti-wrinkling” behaviour is induced by the presence of surface tension which tends to reduce the perimeter of the organoid as proved by the numerical simulations;

• •

the reduction of stiffness of the cell decreases the role of elasticity and enhances the effect of surface tension (leading to an increased $\alpha_{\gamma}$ while keeping $\gamma$ constant). A reduction of stiffness of the organoid has different consequences: the onset of bifurcation is delayed, the critical wavelength increases, the cortex is thicker and, in the most severe cases, the surface instability is absent.

This work opens the path to other studies, such as the analysis of the effect of surface tension on the growth of embryos and tumour spheroid. Our paper proposes a theoretical explanation of purely mechanical nature of known experiments; of course further investigations will be needed to validate it. What emerges from our model is that tissue surface tension cannot be disregarded whenever a living tissue is characterized by relatively small length scales or by small elastic moduli.

Future efforts will include the improvement of the numerical scheme in order to capture possible secondary bifurcations. Possible extensions include the implementation of an arclength continuation method to improve the numerical convergence in presence of turning points. From an experimental point of view, it would be important to quantitatively measure the surface tension acting on the organoid cortex. Another possible line of research is the study of the influence of surface tension on the growth of cellular aggregates.

Acknowledgments

The authors wish to thank D. Ambrosi, P. Ciarletta, A. Marzocchi and G. Noselli for helpful suggestions and fruitful discussions. DR gratefully acknowledge the support of the European Research Council (AdG-340685 MicroMotility).

Appendix A Computation of the incremental curvature

The external boundary in the actual incremental configuration is given by

[TABLE]

since the parametrization of $\bm{\alpha}$ is counter-clockwise (9). By differential geometry, the oriented curvature of $\bm{\alpha}$ is given by

[TABLE]

where with ′ we denote the derivative with respect to $\theta$. By combining Eq. (50) and Eq. (51) we obtain

[TABLE]

We can linearise the relation (52) with respect to $u$, $v$, and their derivatives to get the following expression of the incremental curvature

[TABLE]

Appendix B Role of the Matrigel embedment

Let $k_{s}$ be the elastic constant of the springs per unit length. Let $\mu_{M}\simeq\mu$ be the shear modulus of the Matrigel. The boundary condition Eq. (10) modifies into

[TABLE]

The radially symmetric solution given by Eqs. (21)-(22) is still a solution of the modified problem, the only difference is that the pressure increases of a constant $k_{\mathrm{s}}R_{o}(1-R_{o}/r_{o})$. Denoting by $R_{M}$ the characteristic distance among the organoids, from a dimensional evaluation we obtain

[TABLE]

We adimensionalize $k_{s}$ with respect to $\mu$ and $R_{\mathrm{o}}$, obtaining

[TABLE]

From this computation it is likely that the force exerted by the Matrigel on the organoid is negligible. To verify it, we have modified the linear stability of exposed in Section 4 to account for the Matrigel embedment.

Following the procedure exposed in (Riccobelli and Ciarletta, 2018b), we can compute the marginal stability threshold using the same algorithm exposed in the manuscript by just modifying Eq. (43) into

[TABLE]

The marginal stability threshold are only slightly influenced by the presence of the linear springs at the boundary, as shown by Fig. 11. From these results, it is clear that the presence of Matrigel has very little influence on the behaviour of the organoid, as observed in (Karzbrun et al., 2018).