1 Introduction

The surface of the human brain is characterised by a complex pattern of folds (gyri) and troughs (sulci), allowing for a high area-to-volume ratio [1,2,3]. This intricate structure has increasingly been linked to intellectual ability, marking reduced cortical folding as indicative of cerebral impairments such as “smooth brain" (lissencephaly). Indeed, afflicted human brains exhibit a markedly lower degree of gyrification, resulting in reduced life expectancy and intellectual disability [4, 5]. Recent advancements in the field of stem-cell research provide a controlled, in vitro model system allowing the study of gyrification in the form of brain organoids: cultured, three-dimensional arrangements of pluripotent stem cells replicating some of the key features of human brain development [6, 7]. Although the brain organoid model has been widely used in furthering understanding of a wealth of diseases [8,9,10], the underlying physical mechanism governing gyrification has yet to be pinned down decisively.

In this context, recent work by Karzbrun et al. broke new ground by experimentally probing brain organoids for the onset of an interface instability that results in the formation of folds on the organoid [11] (see Figure 1a). In particular, they show the development of brain organoids over the course of several days, observing the self-organisation of a concentric shell of cells around a spherical cavity (lumen). Due to their active motility, the cell nuclei continually move radially inward and outward, dividing at the inner surface and eventually accumulating at the outer surface. They find prior to the onset of mechanical instability an increase in the density of the nuclei, as well as the aspect ratio, which indicates compression. Taking inspiration from polymer gel models [12], they subsequently argue that these instabilities emerge from the differential swelling of the inner and outer cortex, inducing compressive stress [1, 3].

Building on these observations, Balbi et al. [13] showed that an interplay between the lumen compression and the remodelling of the cortex determines the interface instability of the organoids. Riccobelli and Bevilasqua further showed [14] that surface tension generated by intercellular adhesion in cellular aggregates also contributes to determining the onset of the interface instability. Recently, Engstrom et al. [15] argued that such elastic instabilities based on differential growth sketch an incomplete picture of the folding of brain organoids. Instead, they introduce a system-spanning fibrous model of organoids with an elastic core ensnared by a growing, fluid-like film, suggesting that the details of the microstructure play an important role in the emergence and structure of the wrinkles.

Notwithstanding these important contributions, here we consider a hitherto overlooked aspect of the microstructural complexity of brain organoids in active stress generation by the cells, showing that activity can induce folding at the surface of model organoids in the form of hydrodynamic instabilities. Active stress between cells has been shown to induce surface instabilities in epithelial cell layers [16, 17] and surface deformations in membranes [18]. The model presented here expands on this list and is motivated by experiments in which the brain organoids are treated with cytoskeletal-inhibiting drugs [19], where a marked decrease in the number of folds exhibited and the sharpness of the folds are observed [11]. Since the cytoskeletal filaments inside the cells continuously generate active stress, we conjecture that in addition to a purely elastic wrinkling phenomenon following differential swelling, the activity-driven instability can also contribute to the observed folding behaviour.

The concept of active stress generation by cells is well established in the context of active matter, which describes a material class that exists far from thermodynamic equilibrium by virtue of a local supply of energy. This energy is then converted into work by the constituents, resulting in active stress [20]. Examples include biological systems such as actin filaments [21], microtubule bundles [22] and eukaryotic epithelial cells [23, 24]. Furthermore, recent discrete numerical modelling of epithelial shells [25] has demonstrated the role of activity in inducing morphological changes of the interface by considering the cell–cell tension in the apical and basal surface planes of the cells.

Here, we present a generic, continuum, two-phase framework, in which we describe the brain organoid as an active gel, to study the spatiotemporal dynamics of active, self-deforming surfaces. This generic approach allows us to effectively describe organisation of cells on long time and length scales in terms of the hydrodynamics of the active gel, as well as include the orientational order and contractile activity of the cell cytoskeleton, without reference to the microscopic details [26]. The hydrodynamic nature of the model that accounts for the flow of cells, together with the explicit role of the contractile active stress generation by the cell cytoskeleton, is what sets it apart from those mentioned above.

We perform coarse-grained, continuum simulations in two dimensions, showcasing that the number of folds is dependent on active stress generation by the contractile actomyosin machinery of the cells. An activity threshold needs to be surpassed for folding to occur, similar to experiments where the organoid needs to surpass a critical nuclear cell density for folding to occur. In this vein, our model naturally incorporates the effect of cytoskeletal-inhibiting drugs on the folding behaviour of the organoid, whereas alternative approaches, e.g. based on differential swelling, generally require the assumption of an actively contracted organoid core [13]. Additionally, we present a linear stability analysis of the governing equations in order to characterise the onset of the activity-driven instabilities. Taken together, our results suggest that active stress generation provides a currently overlooked mechanism for cortical folding that is complementary to existing models.

Fig. 1
figure 1

Depiction of brain organoids observed in experiments and the active nematic model we use to illustrate the role of contractile stress in folding. (a) Fluorescence image during organoid growth. Actin filaments found in the organoid cells are coloured green, and the cell nuclei are coloured red. Figure adapted from [11] and the scale bar is \(200\mu m\). (b) Schematic representation of the theoretical model. Filled red ellipsoids represent cell alignment, and the corresponding black arrows indicate contractile stress. The model assumes uniform active stress, which is sufficient to capture the folding of the interfaces. The remainder of the organoid is coloured green to represent the organoid cells. Solid black lines indicate undeformed organoid (low activity), and dashed black lines indicate a deformed organoid (high activity) due to the active stress

Full size image

2 Model

We model the brain organoid (see Figure 1) as a ring of an active gel (with inner radius a and outer radius b) [27], representing the cortex, surrounding a passive isotropic cavity (within the inner interface at radius a), which represents the lumen. Here, we differentiate between the cortex and lumen region by introducing a binary order parameter \(\phi \), with \(\phi \sim 1\) in the cortex and \(\phi \sim 0\) within the lumen region. Within this cortex region, the cells and their associated cytoskeletal filaments are extended radially with clear orientational order (see Figure 1a). To account for the orientational order associated with this microstructural feature of the cortex, we introduce a nematic order parameter \({\mathsf {Q}}\) that is a symmetric, traceless tensor \(Q_{\alpha \beta }=S\left( q_\alpha q_\beta -\frac{1}{2}\delta _{\alpha \beta }\right) \) [28], where S represents the magnitude of the orientational order and q indicates direction [29]. The use of this nematic order parameter is well established in the study of cellular and subcellular systems with elongated constituents such as the cell cytoskeleton [30, 31] and confluent tissues of epithelial or fibroblast cells [32, 33]. By employing this nematic order parameter, our formulation allows for the dynamics of cell alignment within the cortex to be explicitly accounted for. Furthermore, this mapping allows us to describe the cortex as an active nematic system, for which the continuum equations are well established in the literature [34].

Within the active nematic framework, we evolve the binary order parameter \(\phi \) and the orientational order parameter \({\mathsf {Q}}\) through the Cahn-Hilliard [35] and the Beris–Edwards equations [36, 37], respectively:

$$\begin{aligned} \partial _t \phi + \partial _k \left( v_k\phi \right)&= \mu , \end{aligned}$$
(1)
$$\begin{aligned} \left( \partial _t+v_k\partial _k\right) Q_{ij}-\tilde{S}_{ij}&= H_{ij}, \end{aligned}$$
(2)

where \(\varvec{v}\) is the velocity that advects order parameter fields, and \({\mathsf {\tilde{S}}}\) represents the co-rotation terms describing the response of elongated particles to velocity gradients. The latter is given by

$$\begin{aligned} \begin{aligned} \tilde{S}_{ij}=\xi E_{ij}+Q_{ik}\varOmega _{kj}-\varOmega _{ik}Q_{kj}, \end{aligned} \end{aligned}$$
(3)

with the strain rate tensor \(E_{ij}=\left( \partial _iv_j+\partial _jv_i\right) /2\) and the vorticity tensor \(\varOmega _{ij}=\left( \partial _iv_j-\partial _jv_i\right) /2\) describing the symmetric and asymmetric parts of the velocity gradient tensor, respectively. In addition, \(\xi \) denotes a flow-alignment parameter, which determines the collective response of the orientation field to gradients in the velocity field. These relaxation dynamics are governed by minimising the free energy of the system with respect to \(\phi \) and \({\mathsf {Q}}\) through the normalised chemical potential

$$\begin{aligned} \mu&= \varGamma _\phi \left( \frac{\delta \mathcal {F}}{\delta \phi } - \partial _k \left( \frac{\delta \mathcal {F}}{\partial _k \delta \phi } \right) \right) \end{aligned}$$
(4)

and the normalised molecular field

$$\begin{aligned} H_{ij}=- \varGamma \frac{\delta \mathcal {F}}{\delta Q_{ij}}, \end{aligned}$$
(5)

respectively. Here, \(\varGamma _\phi \) is a mobility coefficient and \(\varGamma \) indicates the rotational diffusion coefficient, both of which set the rate of relaxation towards the minimum of the free energy

$$\begin{aligned} \mathcal {F}\left[ {\mathsf {Q}},\phi \right] =\int d^2\varvec{r} \left( f_Q + f_{\nabla Q} + f_{\phi } + f_{\nabla \phi } \right) . \end{aligned}$$
(6)

It can be seen that this free energy includes both bulk and gradient contributions in terms of the tensor order parameter \({\mathsf {Q}}\) and the binary order parameter \(\phi \), which read

$$\begin{aligned} f_{Q}= & {} \frac{1}{2}\mathcal {C}\left( \phi S_n-2Q_{ij}Q_{ij}\right) ^2, \end{aligned}$$
(7)
$$\begin{aligned} f_{\nabla Q}= & {} \frac{1}{2}L\partial _k Q_{ij}\partial _k Q_{ij}, \end{aligned}$$
(8)
$$\begin{aligned} f_\phi= & {} \frac{1}{2}\mathcal {A}\phi ^2\left( 1-\phi \right) ^2, \end{aligned}$$
(9)
$$\begin{aligned} f_{\nabla \phi }= & {} \frac{1}{2}K\partial _k\phi \partial _k\phi , \end{aligned}$$
(10)

where \(S_n\) denotes the equilibrium value of the orientational order parameter S and \(\mathcal {C}, L, \mathcal {A}, K\) are model parameters. For a more detailed discussion of the model, as well as the used parameter values, see the SI.

The evolution of the binary and orientational order parameters is subsequently completed by means of a coupling to the evolution of the velocity field \(\varvec{v}\) that is described by the generalised, incompressible Navier–Stokes equations

$$\begin{aligned} \begin{aligned} \partial _iv_i&=0,\\ \rho \left( \partial _t+v_k\partial _k\right) v_i&=\partial _j\varPi _{ij}, \end{aligned} \end{aligned}$$
(11)

where \(\varPi _{ij}\) describes the full stress tensor that includes pressure and is given in the SI. The dominant terms in the stress tensor are the viscous and capillary stress, which depend on the viscosity of the fluid \(\eta \) and the surface tension \(\sigma \) of the binary-order parameter \(\phi \) [38]. We point out that, in our simulations, the surface tension is no independent model parameter, but rather it is accessed in terms of existing model parameters through the expression

$$\begin{aligned} \sigma =\frac{1}{6}\sqrt{\left( K+\frac{1}{2}L\right) \left( \mathcal {A}_{\text {binary}}+\mathcal {C}_{\text {LQ}}\right) }, \end{aligned}$$
(12)

which can be derived by considering the free-energy cost of the interface (see the SI for details). Similar surface tensions have been calculated for solely \(\phi \)-dependent free energies [35], but here we have expanded the description to account for the \({\mathsf {Q}}\) dependence of the free energy explicitly. For simplicity, we assume surface properties are equal between the outer and inner surfaces.

Fig. 2
figure 2

Temporal evolution of the surface folding. The modelled cortex region starts as a circle (green denotes the organoid region \(\phi =1\)), but due to the onset of an active nematic instability, the surfaces fold over time for high activity (depicted here is \(\alpha =0.0008\)). Displayed timesteps are (0,100000 150000,175000) LB times from left to right. Red lines illustrate coarse-grained orientation field of the cytoskeleton of the cells q

Full size image

Importantly, we also take into account the effect of active stress generated within the cortex. While the activity of the cells is clearly manifested in the brain organoid experiments in the form of active contraction of cells within the cortex [11, 39], rather surprisingly, to our current knowledge, no prior work has explored the impact of this apparent activity in the dynamics and morphology of the organoids. To account for this contractile activity, we introduce an active stress in the form of coarse-grained stresslets that represent the contractile force dipoles that are generated by the actomyosin machinery of the cell cytoskeleton [20, 40, 41] (see Figure 1b). For simplicity, and to show the generic effect of including active stress we assume a uniform spread of dipoles throughout the organoid, resulting in uniform active stress throughout the active ring. Coarse graining over the dipolar force fields leads to an additional, active contribution to the stress tensor [40, 41]

$$\begin{aligned} \varPi _{ij}^{\text {active}}=\alpha Q_{ij}, \end{aligned}$$
(13)

with \(\alpha \) a proportionality constant scaling with activity. The sign of \(\alpha \) determines the nature of the active stress, with \(\alpha > 0\) for contractile and \(\alpha < 0\) for extensile active stress. Due to the contractility of the stem cells, we use \(\alpha > 0\) in the case of the model brain organoid considered here. However, the activity-induced instability is a generic mechanism that is also expected for extensile activities corresponding, for example, to stress generation due to cell division events [42, 43].

The simulations start with zero velocity and the director field oriented along the radial direction (Figure 2a), corresponding to the radially elongated cells between the lumen and the outer interface of the organoid observed in the experiments [11]. Unless otherwise specified, the active region is initialised with the inner interface at radius \(a=45\) and outer interface at radius \(b=170\).

Fig. 3
figure 3

Folding transition through an increase in contractility. The folding number N for different activities and surface tensions found from simulations (a) and the linear stability analysis in a linear (b) and log (c) scale. Simulation snapshots of the organoid for surface tension \(\sigma =0.0014\) and activities \(\alpha =0.000025\) (d), \(\alpha =0.0003\) (e) and \(\alpha =0.0008\) (f). (e-f) are taken just before formation of first topological defects

Full size image

3 Simulation results

We investigate the onset of the interface instability, which resembles gyrification (Fig. 2), by determining the number of folds that form on the initially circular interface of the active ring. We measure the number of folds (also referred to as the folding or wrinkling number N [11]) by first measuring the radial distance of the points on the outer interface from the centre of the ring. The number of folds is then determined from the number of maxima in the radial distance signal as a function of azimuthal coordinates at the onset of the instability. While small perturbations dominate this initial amplitude signal, we notice that a well-defined number of peaks are established relatively quickly. This number remains constant until secondary, nonlinear effects begin to dominate the bulk active system, leading to the nucleation of topological defects and the emergence of active turbulence [34, 44]. To establish the role of activity in the interface instability and determine the number of folds, we focus only on the time span before the bulk instability and the creation of topological defects.

We begin by characterising the number of folds for varying activity and surface tension. The results are represented in the stability diagram (Figure  3a), which clearly demonstrates the competition between the destabilising effect of active stress and the stabilising impact of the surface tension. Increasing activity results in a larger number of folds on the interface, while larger surface tension suppresses the instability and leads to a smaller number of folds. The latter phenomenology is in line with predictions of a purely elastic model of a brain organoid [14]. Interestingly, no folding is observed for sufficiently small activities, indicating that there exists a threshold for the active stress exerted by the cells in order to create folds on the interface. This distinguishes the interface instability in the ring geometry from the well-known thresholdless hydrodynamic instability of unconfined active nematics [40]. More importantly, since we expect the activity to increase with cell nuclear density in the cortex, we conjecture that this observation explains the experimental results of Ref. [11], in which no folding was observed for cell nuclear densities below a critical threshold. Together, these results show that active stress alone can result in hydrodynamic instability of a model organoid, suggesting that activity-induced instabilities can provide a previously overlooked, generic and complementary mechanism to the differential growth mechanism to govern the emergence of folds on brain organoids.

Fig. 4
figure 4

Quantification of the role of activity and organoid thickness. (a) The wrinkling index as a function of activity. For low activity, the wrinkling index is 1 (circular) while it increases linearly after a contractile stress threshold. Measurements are performed at a fixed simulation time \(t=300000\). (b) Wavelength (defined as the inverse of the folding number times the outer radius) for different organoid thicknesses, where the inner radius a is varied while retaining the same outer radius b

Full size image

To draw more parallels between the experimental results on brain organoids and the active nematic ring model presented here, we reproduce two measurements of the experimental paper [11]. First, we measure the wrinkling index L of the active nematic ring for different activities at a fixed time (see Fig. 4a). The wrinkling index is a measurement of the curvature of the interface, defined as the ratio of the contour length normalised to the length of the maximally-protruding outer convex contour. With this definition, \(L=1\) denotes a perfect circle without any folding. We find that for low activities, the wrinkling index remains 1, indicating the active ring remains a circle and no wrinkles are detected. If we increase the activity further, the wrinkling index starts to increase linearly, similar to the linear increase with nuclear cell density found in experiments.

We also measure the effect of changing the active nematic ring thickness t. The experiments of Karzbrun et al. [11] demonstrated that treating the organoid with blebbistatin resulted in a change in organoid thickness and that the wavelength between different convective wrinkles depends on the thickness. We mimic this set-up by varying the inner radius of the isotropic cavity (lumen), while keeping the outer surface radius constant, and define the wavelength \(\lambda \) as the contour length over the folding number N. In line with the experimental results, the wavelength \(\lambda \) increases with the thickness t up until a certain thickness threshold. After this, the wavelength becomes independent of the thickness (see Fig. 4b). Together, these numerical results show that modelling brain organoids as a rings of active nematics can replicate several experimental observations. Next, in order to gain more insight into the nature of the activity-induced instability of the interface, we present a linear stability analysis of the governing equations of the model system.

4 Linear stability analysis

The simulation results point to a possible role of activity in driving the interface instability. To provide a better understanding of the possible activity-induced instability, we perform a linear stability analysis on the surfaces of the model organoid at inner radius a and outer radius b. For simplicity, we consider sharp interfaces, and retain identical surface tensions \(\sigma _a=\sigma _b=\sigma \). Furthermore, in order to allow for further analytical treatment we perform the stability analysis in the limit of overdamped friction. The active forces that put the organoid surfaces under contractile stress originate at the perturbed surfaces, as—in line with experimental observations [11]—we assume the orientation field of the cells retains perpendicular alignment to the interfaces. This corresponds to the limit of strong active anchoring [45]. The perturbations of the interfaces thus directly induce nematic distortions, which in turn give rise to active forces. Through this mechanism, we probe the model organoid surfaces for instabilities that resemble gyrification by applying infinitesimal sinusoidal perturbations of the inner and outer interface of the ring of the form:

$$\begin{aligned} a(t)&=a_0+\delta _a(t)e^{in\theta }, \end{aligned}$$
(14)
$$\begin{aligned} b(t)&=b_0+\delta _b(t)e^{in\theta }, \end{aligned}$$
(15)

where \(\delta (t)\) denotes the infinitesimal perturbation amplitude, n is an integer wavenumber and \(\theta \) denotes the azimuthal coordinate. The wavenumber with the fastest-growing instability is the folding number N that corresponds to the number of folds observed numerically from random perturbations (Figure 3b-c).

We take the expansion around the quiescent state \(\overline{v}_r=0\), \(\overline{v}_\theta =0\) and \(\overline{p}(r)=\int \mathop {}\!\mathrm {d}r\,\alpha S/r\). The perturbations to the surface then result in perturbative corrections to the velocity of the form \(v_r(r,t)=\overline{v}_r+R(r,t)e^{in\theta }\) and \(v_\theta (r,t)=\overline{v}_\theta +\varTheta (r,t) e^{in\theta }\), as well as a perturbative correction to the pressure \(p(r,t)=\overline{p}(r)+P(r,t)e^{in\theta }\). The perturbations evolve according to the Navier–Stokes equations, Eq. (11), where we include only the isotropic and active contributions to the stress tensor. Discarding terms of higher-order than linear in perturbation, we find

$$\begin{aligned} {\left\{ \begin{array}{ll} \begin{aligned} &{}\partial _tR=-\frac{1}{\rho }P'+\frac{\alpha S}{\rho }\frac{n^2-1}{b-a}\left\{ \frac{b-r}{a^2}\delta _a +\frac{r-a}{b^2}\delta _b\right\} \\ &{} \; \; \; \; \; \; \; \; \; \, -\chi R\\ &{}\partial _t\varTheta =-\frac{in}{\rho }\frac{P}{r}-\chi \varTheta \\ &{}R'+\frac{R}{r}+in\frac{\varTheta }{r}=0, \end{aligned} \end{array}\right. } \end{aligned}$$
(16)

where \(\rho \) denotes the density of the active nematic, we introduced the overdamped friction coefficient \(\chi \) to model viscous effects using an argument similar to Stokes’ law [46] (see SI), and primes denote derivatives with respect to the radial coordinate r. The last term on the first line represents the active forces induced by the perturbed surfaces, which we assume decay linearly away from the surface.

Next, we search for exponential solutions of the form \(\propto e^{\omega t}\), where we tacitly assume both surfaces exhibit the same growth rate \(\omega \) [47]. This presupposes a hydrodynamic coupling between the organoid surfaces, in accordance with our numerical simulations, and rationalises the dependence of the wrinkling wavelength \(\lambda \) on the organoid thickness (see Figure 4b). In the SI, we make this coupling explicit.

Upon inserting this ansatz, Eq. (16) reduces to a set of ODEs for the perturbative velocities and pressure, which we solve sequentially. Subsequently, demanding continuity and stress balance at the surfaces of the model organoid (see SI) yields a dispersion relation for the growth rate \(\omega \):

$$\begin{aligned} \omega ^4+\frac{2\chi }{\rho }\omega ^3+A\omega ^2+\chi B\omega +C=0, \end{aligned}$$
(17)

where the coefficients A, B, C embed dependency on the wave number and model parameters (see the SI for explicit form of these terms).

To determine the stability of a perturbation with wavenumber n, we inspect the real part of the corresponding growth rate \(\omega \) and associate the wrinkling number N with the fastest-growing mode. Interestingly, numerical checks of our analytical results show that all unstable surface perturbations are of the same sign, indicating an undulation mode in accordance with experiments [11]. Furthermore, we recover a clear folding transition, as—for a given surface tension \(\sigma \)—nonzero wrinkling numbers are only recovered past a critical activity-to-surface tension ratio (see Figure 3b-c and the SI). This is in agreement with the numerical simulations found in Figure 3a, as well as with the experiments of Ref. [11], where no wrinkling instability is observed for low cell nuclear densities.

5 Conclusion

We present an active matter model system for describing the experimentally observed folds on brain organoids. By treating the organoid as a contractile, active, biphasic annulus, we demonstrate a generic instability mechanism for the interface deformation due solely to active stress generation by the cytoskeleton of the cells within the active cortex that surrounds a passive lumen without introducing growth, swelling, or any other planar surface forces. Combining numerical simulations with linear stability analyses, we show that the activity-induced instability mechanism occurs over a well-defined activity threshold that depends on the mechanical properties of the system, such as surface tension and the thickness of the model organoid. Since the active stress generation is expected to increase with increasing cell nuclear density within the cortex, these results are in line with the experimental measurements that establish the emergence of folds on the organoid surface above a critical nuclear density [11].

Furthermore, the increase in folding number N upon increasing activity \(\alpha \) (see Figure 3) lends credence to the proposed model in light of molecular perturbation experiments on blebbistatin-treated brain organoids that show a marked decrease in the number of folds the organoids exhibit upon treatment with the contractility-inhibiting drug [11].

The activity-induced mechanism proposed here complements previously suggested mechanisms of differential swelling between the inner and outer cortex [48], adhesion-induced surface tension [14], and the fluid-elastic fibrous model [15] for the generation of organoid surface deformations.

It is important to note that the model presented here provides a generic framework based on accounting for (i) microstructural complexity of the organoid cortex based on the orientational order of the cells, (ii) the biphasic nature of the organoid based on an active, viscoelastic cortex surrounding a fluid-like lumen region, and (iii) active stress generation by the cells. As such, the model allows for isolating the effects of active stress generation from potential impacts of growth-induced mechanisms on the formation of surface instabilities. A more complete picture of the relative importance of different mechanisms leading to brain organoid deformation should allow for activity-induced, growth-induced, and differential adhesion-induced mechanisms to be accounted for in one unifying framework, which should be the focus of future studies.