{'fraction_non_alphanumeric': 0.034601823708206685, 'fraction_numerical': 0.017458966565349544, 'mean_word_length': 4.506963042313872, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Surface tension controls the onset of gyrification in brain organoids

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


INTRODUCTION
The formation of folded structures in human and animal brains allows 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] apparently confirm 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, related to the LIS1 heterozygous (+/−) mutation, is correlated to nutritional disorders, alterations in muscle tone and severe psychomotor and mental retardation [Dobyns, 1993]. The mathematical description of the brain sulci embryogenesis can provide new insights to understand the mechanisms underlying this disease. Figure 1: Cells in the bulk (left) and on the free surface (right) of a cellular lattice. Adhesion forces generated by the surrounding cells are denoted by the red arrows. The sum of all these forces on an internal cell is zero, while it is non-zero and perpendicular to the boundary for a cell on the free surface.
A well developed framework to model growth induced mechanical instabilities 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 to elastically distort the body and generate residual stresses [Hoger, 1985, Rodriguez et al., 1994. Differential growth and residual stresses are involved in the morphogenesis of tissues [Ambrosi et al., 2017a] 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 . 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].
Such mechanical unbalance 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. The introduction of a small surface tension regularizes the problem, penalizing the short wavelength modes Ciarletta, 2010, Dervaux andBen Amar, 2011]. At small length scales, elasto-capillary forces can deform soft gel beams [Mora et al., 2010[Mora et al., , 2013 even inducing mechanical instabilities [Taffetani andCiarletta, 2015, Xuan andBiggins, 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.

GREGATES
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, and so on. From a macroscopic point of view, cellular aggregates 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., 2017b, 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., 2017b].
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, also solids possess surface tension [Style et al., 2017] and 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 µ be the shear modulus of the cellular agglomerate and γ the surface tension, then the capillary length l c is defined as [Mora et al., 2013, Style et al., 2017 Whenever this length-scale is of the same order of 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.

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 we assume that the spheroid is composed of an incompressible elastic material. Let T be the Cauchy stress tensor, the balance of the linear and angular momentum reads div T = 0.
If the spheroid is subjected to surface tension γ, the boundary condition reads [Style et al., 2017] where R denotes the radial position and K is twice the mean curvature. We prove that the undeformed reference configuration is in mechanical equilibrium. Assuming that residual stresses are absent [Lee et al., 2019], the Cauchy stress is given by [Ogden, 1997] T = −pI where 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). Assuming that a surface tension γ acts at the external boundary of the spheroid, the boundary condition reads of the surface and n is the external normal. From Eq. (2) we get which is nothing else that the Young-Laplace equation. We obtain an expression for γ, such as From the work of Lee et al. [2019], we get that the typical radius of a spheroid is ∼ 400 µm and the internal pressure p is about 500 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 γ 0.1 N/m.

ELASTIC MODEL OF BRAIN ORGANOIDS
In this Section, we illustrate a model of brain organoids, described as growing hyperelastic bodies subjected to surface tension

KINEMATICS
We denote by X the material position coordinate. The experiments of Karzbrun et al. [2018] suggest to adopt a two-dimensional domain equipped with polar coordinates. Let be the reference configuration of the organoid. We indicate with ϕ : Ω 0 → R 2 the deformation field, so that the actual configuration of the body Ω is given by ϕ(Ω 0 ). Then x = ϕ(X) be the actual position of the point X and the displacement vector is defined as u(X) = ϕ(X) − X.
Let F be the deformation gradient, i.e. F = Grad ϕ. We exploit a multiplicative decomposition of the deformation gradient (known as Kröner-Lee decomposition [Kröner, 1959, Lee, 1969) to model the growth, so that where the growth tensor G accounts for the local inelastic distortion due to the body growth while F e describe the elastic distortion necessary to maintain the geometrical compatibility of the body and to balance the external and internal forces [Rodriguez et al., 1994]. Since human cells are mainly composed of water, it is reasonable to model the organoid as an incompressible medium, namely we enforce that 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 i the radius of the lumen, we denote these two regions by Ω 0L and Ω 0R : and their images through ϕ are denoted by Ω L and Ω R respectively. We assume that the growth tensor G takes diagonal form where the scalar quantity g is the growth rate of the cortex with respect to the lumen. We now introduce some mechanical constitutive assumptions.

MECHANICAL CONSTITUTIVE ASSUMPTIONS AND FORCE BALANCE EQUATIONS
We assume that the organoids are composed of an incompressible hyperelastic material with strain energy density ψ. The first Piola-Kirchhoff stress P and the Cauchy stress tensors T are then given by where p is the Lagrange multiplier enforcing the incompressibility constraint (4). The balance of the linear and angular momentum reads in the material and actual reference frame, respectively. We assume that the center of the organoid is fixed, i.e.
while a constant surface tension γ acts at the external boundary of the organoid, so that [Style et al., 2017] where n is the outer normal in spatial coordinates and K is the oriented curvature of the boundary curve L parametrized clockwise, i.e.
The Lagrangian form of the boundary condition (8) is obtained performing a pull-back Finally, we enforce the continuity of the stress at the interface R = R i , so that We assume that the organoid is composed by a neo-Hookean elastic material, so that the strain energy density is given by where I 1 is trace of the right Cauchy-Green tensor C = F T F. The stress tensors P and T read 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.

EQUILIBRIUM RADIALLY-SYMMETRIC SOLUTION
Let (r, θ) be the actual radial and polar coordinates of a point. Let (E R , E Θ ) and (e r , e θ ) 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 The deformation gradient expressed in polar coordinates reads It is immediate to notice that In the cortex, from the incompressibility equation (4), we get r r = R.
Performing an integration and imposing that R i = r i , we get It remains to determine the pressure field p. First, we notice that, inverting and differentiating (15), we obtain respectively. The curvature of the boundary line is so that the boundary condition (8) reads Since the deformation depends only on the radial position r, the balance of the linear and angular momentum in polar coordinates reads where T ij are the components of the Cauchy stress tensor T in polar coordinates. From (13), (14) and (16), the Cauchy stress in the cortex reads We can integrate (17) from r o to r, obtaining We can find the pressure field in the cortex (i.e. for r i < r < r o by plugging Eq. (18) into Eq. (19), obtaining Finally, we impose the continuity of the stress (11) at r = r i = R i to get the pressure for r < r i . Since the lumen remains undeformed, the Cauchy stress reads T = −p L I, and, using Eq. (18) and Eq. (20), we can write the solution In the next section a linear stability analysis of the solution given by Eqs. (21)-(22) is perform.

INCREMENTAL EQUATIONS
In order to investigate the stability of the radially symmetric solution, we apply the theory of incremental deformations superposed on finite strains [Ogden, 1997]. Let δu be the incremental displacement field and let Γ = grad δu. We introduce the push-forward of the incremental Piola-Kirchhoff stress in the axisymmetric deformed configuration; such a tensor is given by where A 0 is the fourth order tensor of instantaneous elastic moduli, δp is the increment of the Lagrangian multiplier that imposes the incompressibility constraint. The two dots operator (:) denotes the double contraction of the indices where the convention of summation over repeated indices is used.
The components of the tensor A 0 for a neo-Hookean material, are given by where B e = F e F T e and δ ik is the Kronecker Delta. The incremental equilibrium equation and the linearised form of the incompressibility constraint are given by The linearised form of the kinematic constraint (7) and of the boundary condition (8) complement the incremental equations where δK is the increment of the curvature. Finally, we enforce the continuity of the incremental displacement of the stress at the interface, i.e.
In the following, we rewrite the incremental problem given by the Eqs. (24)-(26) into a more convenient form using the Stroh [1962] formulation.

STROH FORMULATION
We rewrite the incremental problem in non-dimensional form using the length scale R o and the shear modulus µ. The behaviour of the problem is governed by the non-dimensional parameters in addition to the growth parameter g.
For the sake of brevity, we introduce the multi-index W = {L, 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 W and v W the components of δu W while δP W rr and δP W rθ are the components of the incremental stress projected along the radial vector e r . We can reduce the system of partial differential equations (24) into a system of ordinary differential equations by assuming the following ansatz for the incremental displacement, pressure and stress: where m ∈ {n ∈ N | n ≥ 2} is the circumferential wavenumber and U W , V W , s W rr and s W rθ are scalar functions of r. By substituting (31) into (24), we obtain the following expression for Q L and Q C Q L (r) = −s L rr (r) + p L U L (r) where p L , p 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 where η W is the displacement-traction vector: The matrix N W ∈ R 4×4 is the Stroh matrix and it has the following sub-block form For the lumen (r < r i ), the sub-blocks read: In the cortex (r i < r < r o ), the sub-blocks are given by where β C = r 2 r 2 + (g 2 − 1)α 2 R . Interestingly, 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 i .

INCREMENTAL SOLUTION FOR THE LUMEN
We follow the procedure used in [Dervaux and Ben Amar, 2011] and in [Balbi et al., 2018]. Since Eq. (34) with W = L is a system of ODEs with constant coefficients, its solution can be rewritten in terms of eigenvalues and eigenvectors of N L .
The eigenvalues of N L are λ 1 = m − 1, λ 2 = m + 1, λ 3 = −m + 1 and λ 4 = −m − 1. The general integral of Eq. (34) is given by where w i are the eigenvectors of N L associated with the eigenvalues λ i , i = 1, . . . , 4. Since the incremental solution must satisfy the kinematic constraint (25), we immediately get that c 3 = c 4 = 0 while 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 i (i.e. by enforcing the equations (27)).

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 N C given by (36) are not constant. To overcome this difficulty, we implement a numerical code to solve based on the impedance matrix method [Biryukov, 1985, Biryukov et al., 1995. We first introduce the conditional impedance matrix, defined as where the name conditional derives from the fact that its expression depend on an auxiliary condition at r = r o , which in this case is provided by the boundary condition (26). In the following paragraphs we expose a procedure to construct the matrix Z C (r, r o ). The incremental boundary condition Eq. (26) reads where the incremental curvature δK is given by (see the Appendix A for the details of the computation): From the ansatz of variable separation in Eqs.(29)-(30) and using the incremental form of the incompressibility constraint (24), the boundary condition (40) is equivalent to We can now define the auxiliary impedance matrix [Norris and Shuvalov, 2010] as so that the boundary condition (41) is equivalent to the equation We introduce the matricant called conditional matrix, defined as the solution of the problem From the definition (44), the Stroh form of the incremental problem (34) and the Eq. (43), we get From (45), the conditional impedance matrix is given by From now on, we omit the dependence of Z C wherever convenient for sake of simplicity. By using Eq. (39), we can rewrite the Stroh problem in the cortex (34) into a differential Riccati equation. Indeed, from (34), with W = C, we obtain Substituting Eq. (46) into Eq. (47) we get the Riccati equation We integrate Eq. (48) from r o to r i , using as initial condition the auxiliary impedance matrix defined in (42) To construct a bifurcation criterion, we follow Balbi et al. [2018]. From the continuity of the displacement-traction vector η C (r i ) = η L (r i ) (27) we get

so that non-null solutions of the incremental problem exist if and only
where with w i defined in (38). For a fixed value of the control parameter g we integrate the Riccati equation (48) from r = r o up to r = r i making use of the the software MATHEMATICA 11.3 (Wolfram Research, Champaign, IL, USA). We iteratively increase the control parameter g until the bifurcation criterion (49) is satisfied.

DISCUSSION OF THE RESULTS
First, we need to identify an interval of interest for the adimensional parameter α γ . 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 N/m. According to Karzbrun et al. [2018], the shear modulus of the wild-type brain organoid is µ 900 Pa (Young modulus E 2.7 kPa) while µ 333 Pa (Young modulus E 1 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 o = 400 µm [Karzbrun et al., 2018], α γ ranges between 0.25 for the wildtype 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 α R and α γ , 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 cr as the minimum g m for m ≥ 2 and the critical mode m cr as the wavenumber corresponding to g cr .
In Fig. 2, we plot the critical values g cr and m cr versus α γ for two different values of α R , i.e. in the first one α R = 0.9, while in the other one α R = 0.95.
We observe that, for relatively small values of α γ , the marginal stability threshold g cr increases monotonously as α γ increases, while the critical wavenumber m cr decreases. There is a change in the behaviour of the instability when the parameter α γ is sufficiently large: the critical wavenumber m cr increases suddenly and the marginal stability threshold g cr remains nearly constant about g cr 2.5. The threshold at which this transition occurs strongly depends on α R , as shown in Fig. 2. Indeed, when α R = 0.9, g cr increases from 1.745 to 2.481 with α γ ∈ (0, 1.34), while when α R = 0.95, g cr Table 1: Solutions of the linearised incremental problem at different α R and different α γ . The amplitude of the incremental radial displacement has been set equal to 0.05R o for the sake of graphical clarity.
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 α γ and α R . For small values of α γ , we observe that the instability mainly releases elastic energy at the free boundary, displaying a wrinkling pattern: the wavenumber decreases and the critical mode displays a more rounded boundary as we increase α γ . Furthermore, for large values of α γ 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.
We also investigate the influence of α R on the instability by fixing α γ (see Fig. 3). In the absence of surface tension (i.e. α γ = 0), we observe that g cr decreases monotonously as α R increases (see Fig. 3a). The behaviour is the opposite in the presence of surface tension, where the marginal stability threshold g cr monotonously increases with α R . Moreover, as one can observe from Figs. 3c-3d, the transition from a surface instability to an interfacial instability is anticipated when α γ increases. For the range of parameters in which the critical mode displays an interfacial instability, we observe that g cr increases linearly with α R . As it concerns the critical wavenumber m 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 correctly captures the main features of organoid development. First, for small values of α γ , a morphological transition takes place at the free boundary, giving rise to a wrinkling pattern: as α γ increases, we observe a decrease of the critical wavenumber and an 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 observed that the elastic modulus of the cells is 2.7 times lower than the one of healthy cells, so that α γ is much larger. 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-3 and from the morphological diagram of Table 1, as we increase α γ the number of wrinkles decreases and the critical threshold increases, in accordance with Karzbrun et al. [2018]. Furthermore, for large values of α γ , 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 α γ .
Notwithstanding the good agreement with the experimental results of our model, it is to be reported that solid models of brain organoids have been criticized recently by Engstrom et al. [2018]. The authors interestingly 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 wrin-kles. 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.

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, 2018]: where (e 1 , e 2 ) represents the canonical vector basis. Given X n ∈ Ω n , the components represent the referential radial coordinate normalized with respect to the external radius, and the referential polar angle, respectively: X n 2 = Θ, as represented in Fig. 4. The function 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 n 1 = 0 up to X n 1 = α R . In the cortex we instead use a structured mesh (i.e. for α R < X n 1 < 1), see Fig. 5. We implement a numerical code based on the mixed finite element method to enforce the incompressibility constraint (4). We discretise the displacement field u using continuous, piecewise quadratic functions, while we approximate the pressure through piecewise constant functions. The corresponding mixed finite element is the 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 µ and R o as we did at the continuum level. We set the following boundary conditions We solve the discretised form of the equilibrium equation (6) in the Lagrangian form using a Newton method. The control parameter g o is incremented incremented of δg when the Newton method converges, the found numerical solution is used as initial guess for the new Newton cycle. The increment δg is automatically reduced near the theoretical marginal stability threshold and when the Newton method does not converge. The numerical simulation is stopped when δg < 10 −6 . To trigger the mechanical instability, a small perturbation of an amplitude of ∼ 10 −5 is applied at the free boundary of the mesh.
The numerical algorithm is implemented in Python through the open-source computing platform FEniCS (version 2018.1) [Logg et al., 2012]. We use PETSc [Balay et al., 2018] as linear algebra back-end and MUMPS [Amestoy et al., 2000] as linear solver.

RESULTS OF THE FINITE ELEMENT SIMULATIONS
In this section, we discuss the results of the numerical simulations for α R = 0.9. In Fig. 6, 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. 6b) the free boundary is smoother and rounded.
We now show a quantitative analysis of the numerical solution depicted in Fig. 6. We define ∆r as the amplitude of the pattern at the free surface, i.e.
where r denotes the actual radial position of the point with polar coordinates (R, Θ).
In Fig. 7, we show how the amplitude of the pattern ∆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.   Figure 7: Bifurcation diagrams for (a) α γ = 0 and (b) α γ = 0.5. The orange square denotes the theoretical marginal stability threshold computed as exposed in Section 4. The good agreement between the linear stability analysis and the finite element code outcomes validates the numerical algorithm.
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 be the total mechanical energy of the organoid. We compute the ratio of the energy E th of the base solution given by Eqs. (21)-(22) and the energy E num arising from the numerical simulations.
In Fig. 8 we plot the ratio E num /E th versus the control parameter g for both α γ = 0 and α γ = 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).
Interestingly, we observe that the thickness of the ridges is higher than the one of the furrows for α γ = 0 (see Fig. 9a), while the behaviour is the opposite in the case of α γ = 0.5 (see Fig. 9b).   Figure 8: Plot of the ratio E th /E num for (a) α γ = 0 and (b) α γ = 0.5. The orange square denotes the theoretical marginal stability threshold computed as exposed in Section 4.   Figure 9: Thickness of the cortex of ridges (blue) and furrows (green) for (a) α γ = 0 and (b) α γ = 0.5. The latter situation, in which the the thickness of ridges is higher that the one of furrows, corresponds to the "antiwrinkling" behaviour described in [Engstrom et al., 2018].
The latter case is tipical of brain organoids as observed in [Engstrom et al., 2018]. These numerical results prove 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 tipical size is of order of decimeters. In this case, keeping γ and µ fixed, α γ reduces of 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 α γ = 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. 7 show that the amplitude of the pattern increases very rapidly beyond the marginal stability threshold. As one can observe in Fig. 6a, 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 particularily 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.

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. The physical basis of the model is stated in Section 2, where we have estimated the magnitude of surface tension acting on a solid multicellular aggregate. Using the experimental data reported in [Lee et al., 2019], we have showed that the tissue surface tension acting on a solid multicellular aggregate is of the order of 10 −1 N/m.
We have then 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 cell 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 homogeneous growth rate of the cortex, α R and α γ , which are the lumen radius and the capillary length respectively, normalized with respect to the initial radius of the organoid.
We have computed a radially symmetric solution and we have studied its linear stability in Section 4. We have employed analytical and numerical tools, such as the theory of incremental deformations, the Stroh formulation and the impedance matrix method. In particular, we have proposed a novel extension of the last method to take into account surface tension. The outcomes are discussed in Section 4.5. Surface tension influences both the critical wavenumber m cr and the marginal stability threshold g cr . Experimental observations of Karzbrun et al. [2018] reported that, when the LIS1 +/− mutation is present, the shear modulus of the organoid is reduced and the critical wavenumber decreases compared with the healthy organoids. In our model the softening of cells, due to lissencephaly, corresponds to an increase of the parameter α γ , so that surface tension is predominant and the critical wavenumber decreases, in accordance with the experiments performed by Karzbrun et al. [2018]. Moreover, we have observed that, for larger α γ , 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 and in Tab. 1.
Finally, in Section 5, we have implemented a finite element code to approximate the fully non-linear problem. The numerical algorithm is based on a mixed variational formulation and a Newton method. In Fig. 6-9 the outcomes of the numerical simulations are reported. The results show that the 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. Moreover, we show that the "antiwrinkling" behaviour of the cortex (namely the thickness of the outer layer is larger in the furrows) is due to the presence of a surface energy. This strengthens the importance of considering surface tension in the modelling of cellular aggregates.
Summing up, in this paper we have developed a new model of brain organoids based on the theory of non-linear elasticity coupled with active phenomena, such as growth and the tissue surface tension. The outcomes of the model provide a novel explanation of the experimental results reported in [Karzbrun et al., 2018], showing that lissencephaly is caused by a competition between surface tension and bulk elastic energy. Indeed, organoids affected by lissencephaly are much softer than the healthy ones, so that the role of surface tension is predominant versus bulk elasticity.
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.