PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ploscompComputational BiologyView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS Comput Biol. 2010 April; 6(4): e1000749.
Published online 2010 April 22. doi:  10.1371/journal.pcbi.1000749
PMCID: PMC2858670

A Reaction-Diffusion Model of Human Brain Development

Karl J. Friston, Editor

Abstract

Cortical folding exhibits both reproducibility and variability in the geometry and topology of its patterns. These two properties are obviously the result of the brain development that goes through local cellular and molecular interactions which have important consequences on the global shape of the cortex. Hypotheses to explain the convoluted aspect of the brain are still intensively debated and do not focus necessarily on the variability of folds. Here we propose a phenomenological model based on reaction-diffusion mechanisms involving Turing morphogens that are responsible for the differential growth of two types of areas, sulci (bottom of folds) and gyri (top of folds). We use a finite element approach of our model that is able to compute the evolution of morphogens on any kind of surface and to deform it through an iterative process. Our model mimics the progressive folding of the cortical surface along foetal development. Moreover it reveals patterns of reproducibility when we look at several realizations of the model from a noisy initial condition. However this reproducibility must be tempered by the fact that a same fold engendered by the model can have different topological properties, in one or several parts. These two results on the reproducibility and variability of the model echo the sulcal roots theory that postulates the existence of anatomical entities around which the folding organizes itself. These sulcal roots would correspond to initial conditions in our model. Last but not least, the parameters of our model are able to produce different kinds of patterns that can be linked to developmental pathologies such as polymicrogyria and lissencephaly. The main significance of our model is that it proposes a first approach to the issue of reproducibility and variability of the cortical folding.

Author Summary

The anatomical variability of the human brain folds remains an unclear and challenging issue. However it is clear that this variability is the product of the brain development. Several hypotheses coexist for explaining the rapid development of cortical sulci and it is of the highest interest that understanding their variability would improve the comparison of anatomical and functional data across cohorts of subjects. In this article we propose to extend a model of cortical folding based on interactions between growth factors that shape the cortical surface. First the originality of our approach lies in the fact that the surface on which these mechanisms take place is deformed iteratively and engenders geometric patterns that can be linked to cortical sulci. Secondly we show that some statistical properties of our model can reflect the reproducibility and the variability of sulcal structures. At the end we compare different patterns produced by the model to different pathologies of brain development.

Introduction

The development of the human brain from the early gestational weeks to the buckling of the first folds at around 20 weeks follows a narrow pathway between determinism and pure randomness. On the one hand normal adult individuals offer quite similar - from a pure qualitative and descriptive point of view - folding structures: gyri and sulci. On the other hand we observe morphological variabilities between different brains [1]. This variability can reach extreme states in the case of rare abnormalities of the developing brain - such as lissencephaly, polymicrogyria or corpus callosum agenesis. The origin of variability remains an unclear and challenging issue [2] but it is, however, obvious that environmental factors have a deep impact on the sulcal and gyral pattern since even monozygotic twins exhibit important anatomical differences [3].

Folding or buckling are very general processes in nature and among living organisms. Especially one of the most studied step in the morphogenesis of metazoans is gastrulation which corresponds to a symmetry breaking of the spherical embryo and an invagination. The origin of this folding remains unknown even if mechanical factors are undoubtedly implied [4]. More disconcerting, it is shown in [4] that different mechanical actions (constriction, contraction, traction, gel swelling) can lead to similar shapes of the sea urchin primary gastrula.

In these conditions it raises the issue of realistic modeling of far more complex buckling processes such as the gyrification of mammal brains. In this regard it is important to inspect carrefully previous models of gyrification.

Le Gros Clark [5] raises first that the cortex grows by surface expansion rather than by increasing its thickness. He suggested that the expansion of the brain is constrained by the skull and basal ganglia and that compressive stresses cause sulcation. However experiments on sheeps whose large quantities of cortical and subcortical structures were ablated at the end of cellular migration revealed, at term, gyri and sulci of normal size and configuration [6]. This model and its refutation give us a way to categorize the hypotheses on the gyrification depending on whether they involve intra- or extra-cortical processes, in other terms intrinsic or extrinsic.

In the same extra-cortical point of view a recent and very popular model considers that the folding of the brain takes its origin in the mechanical tensions produced by the white matter fibers [7]. This model has been recently tested in [8] with a finite element model of cortical folding.

At the opposite there are numerous hypotheses arguing that the cortical folding has intrinsic origins. In [9] the differential growth of cortical layers causes sulcation and can explain anomalies of folding such as polymicrogyria and lissencephaly in terms of different mechanical properties of cortical tissues. Other models use mechanical hypotheses on the cortex such as elasticity or plasticity [10], [11]. In particular in [11] the authors suggest that the cortical folding is only a consequence of its growth modulated by anisotropies in mechanical or geometric properties. Even if this last model reproduces several characteristics of a growing cortex, it remains however implemented in 2D and does not explain completely where the anisotropies come from.

In the intrinsic origins of folding we encounter also purely morphogenetic hypotheses in which cortical convolutions are under genetic control [12]. In [13] it is proposed that the different cytoarchitectonic areas are provided by a protomap, that is a layer of predetermined neuronal units.

In the next part we will see in detail another hypothesis for the cortical folding which is based on Turing instabilities [14], [15]. We will show that these approaches of brain development can also be linked to the sulcal roots model proposed in [2]. These last authors offer indeed a descriptive model of the human sulcation based on the concept of sulcal roots that are elementary atoms around which the brain folding organizes itself. This concept has strong similarities with the sulcal pits one [16], [17]. Sulcal pits correspond to the deepest points of the sulci whose reproducibility has been demonstrated rigourously in [17]. We will see that initial conditions of the reaction diffusion process can have an interpretation in terms of sulcal roots or sulcal pits.

In this article we investigate the origin of anatomical variability from the early development and we propose a phenomenological model of the folding which is based on the putative existence of Turing morphogens. After recalling briefly some mathematical aspects of the model, we present the numerical schemes used for implementing the equations on a surface and for the deformation of the surface. We show some qualitative and quantitative results of the model. In particular we link sulcal pits maps to the average folding patterns across several realizations of a same noisy initial condition. And we study the variability of our model and demonstrate that it can lead to different modes of variability of one sulcus.

Model

Cartwright [14] suggested first a similarity between labyrinthine Turing patterns and brain gyri. He proposed to model the cortical anatomy as obtained by a reaction-diffusion mechanism and more especially by Van Der Pol-Fitz Hugh Nagumo equations. The mechanism will correspond to axonal pathfinding with diffusing chemicals that activate axon growth or inhibit it [18].

Very recently the Global Intermediate Progenitor (GIP) model has been proposed to explain the appearance of transversal or sectorial sulci following the Intermediate Progenitor hypothesis [19]. The GIP model is based on the BVM system of reaction-diffusion equations [20] that mimics the patterning of the subventricular zone.

We aim at extending the analogy first formulated by Cartwright and the GIP model using a system of reaction-diffusion equations that will modify the surface on which the equations take place. Namely the reaction diffusion system models the non-linear interaction of two morphogens An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e001.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e002.jpg or in other terms growth factors. These growth factors are characterized by the fact that they will tend to deform the surface on which they evolve. We can note that this approach has been previously used in the case of plant growth [21] [22] and that reaction-diffusion equations on surfaces have also been adopted to apply textures on meshes [23] [24].

In our model initial conditions of the reaction-diffusion equations have strong similarities with the sulcal roots described in [2], that are the initial seeds of the folding process.

Our choice of the reaction-diffusion equations differs from the one of Cartwright and the GIP model since we have adopted the Gray-Scott model:

equation image
(1)
equation image
(2)

which exhibits a high number of patterns for differents values of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e005.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e006.jpg [25].

The mathematical analysis of the model has been previously conducted in [26]. As in [15] we have adopted a phenomenological approach and it is interesting to note that our model and the BVM one differ only by a quadratic term. Linear and cubic terms in the kinetic reactions are present in the two models.

Surface deformation

Following ideas of [21] and [22], we suppose that the evolution of the studied surface An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e007.jpg is driven by the morphogens An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e008.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e009.jpg. An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e010.jpg is the inhibitor and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e011.jpg the activator. In mathematical terms we have that

equation image
(3)

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e013.jpg is the normal to the surface and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e014.jpg a function of the two morphogens. The simplest case for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e015.jpg that we have adopted in the following is a linear function of one morphogen:

equation image

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e017.jpg is a parameter in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e018.jpg.

Since the surface on which evolve the morphogens is modified with time, we have to adapt the equations (1) and (2) to take into account the geometric changes. The problem of reaction-diffusion on growing domains has been well-studied in the past years. It leads generally to add convective and dilution terms to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e019.jpg (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e020.jpg respectively) that can be combined in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e021.jpg where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e022.jpg represents the flow velocity of the growing surface [27]. However this result does not directly apply to surfaces and we have to refer to [28] to see the influence on the curvature changes on the reaction-diffusion equations.

The model proposed in [28] consists in adding a term reflecting the modification of the surface metric along time. If the surface An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e023.jpg is parameterized by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e024.jpg then equations (1) and (2) read:

equation image
(4)
equation image
(5)

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e027.jpg stands for the Laplace-Beltrami operator. An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e028.jpg is the determinant of the metric An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e029.jpg associated to the surface, that is:

equation image

In the following we will use An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e031.jpg instead of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e032.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e033.jpg instead of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e034.jpg for simplicity reasons but one has to remember that the surface on which the equations are defined is changing along time.

Numerical implementation

Since we work on discrete meshes we have used a finite element method to discretize the linear terms in the equations 1 and 2.

First we derive a weak formulation of the system on

equation image

with non-linear terms included in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e036.jpg using a test function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e037.jpg:

equation image

Then, integrating by part the Laplacian term, as in [29], we get:

equation image

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e040.jpg is the metric associated to the Riemannian manifold An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e041.jpg.

Next we work on a discrete tessellation An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e042.jpg of the surface An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e043.jpg composed of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e044.jpg vertices. We define An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e045.jpg functions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e046.jpg which are continuous piecewise affine, with the property to be equal to 1 at node An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e047.jpg and 0 at all other triangle nodes. They are the basis functions for the approximation on the functional space of finite dimension An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e048.jpg. So any function An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e049.jpg continuous piecewise affine reads: An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e050.jpg.

The weak formulation becomes:

equation image

It is possible to treat the non-linear term with the following approximation as in [30]:

equation image

then discretizing implicitly and explicitly in time between An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e053.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e054.jpg and writing with matricial expressions:

equation image

with

equation image
equation image

and by definition

equation image

So we can deduce:

equation image

On each triangle we have

equation image
equation image

where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e062.jpg is the area of triangle An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e063.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e064.jpg is the height of triangle An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e065.jpg from vertex An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e066.jpg.

At last we need to compute An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e067.jpg. This can be performed on each triangle An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e068.jpg from the expression of the metric tensor An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e069.jpg [31]:

equation image

Remarks

To be exact we should take into account the fact that the finite elements vary along time or in other terms that the mass matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e071.jpg depends on the evolving surface. This point is addressed in [32] but it is not possible to apply this framework in our case since we should know the geometry of the mesh at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e072.jpg to compute the mass matrix An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e073.jpg which is not the case because the mesh at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e074.jpg is deduced from An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e075.jpg. However the authors in [32] propose an equivalent weak formulation in which the velocity term is present. In our case the dilution term replaces the velocity term. Moreover the expression of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e076.jpg should involve An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e077.jpg but for the same reason as previously we have not the knowledge of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e078.jpg at time n.

For the surface deformation step, we translate equation (3) by simply modifying each vertex An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e079.jpg:

equation image
(6)

Although this incremental procedure can rapidly lead to abnormal deformation of the original mesh, we have used triangle refinements in order to prevent this issue. When the area of a triangle exceeds two times the averaged area of the triangles of the original mesh we simply divide the triangle in four triangles constructed from the three midpoints of each side Fig. 1.

Figure 1
Mesh refinement procedure.

In our implementation we do not prohibit self-intersection which would increase considerably the computation time. However we can say that we escape this issue by not solving for too long time but also taking a parameter An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e081.jpg for the deformation along the normal not too large. By this we mimic also the real brain expansion that is not confronted to the problem of gyri collision.

Results

Labyrinthine patterns

First we can model the growth of a normal brain with the value An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e082.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e083.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e084.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e085.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e086.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e087.jpg. The initialization corresponds to a slight perturbation of the stable equilibrium An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e088.jpg in a position of a sphere composed of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e089.jpg vertices. This perturbation consists of a broad line with An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e090.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e091.jpg where An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e092.jpg is white noise of amplitude An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e093.jpg (see first picture on Fig. 2).

Figure 2
Evolution of the Gray-Scott model coupled to surface deformation.

Note that the evolution specified by the coupled reaction-diffusion equations and the surface deformation leads to a progressive folding of the initial sphere on Fig. 2. It is possible to extract an order parameter at each time step which consists of the number of folds or sulci (see Fig. 3 left). This index is defined from the curvature map of the surface. At each vertex of the mesh we compute the mean curvature following [33]. Once this curvature map has been obtained we compute automatically the number of sulci, that is, the number of connected components whose curvature is inferior to 0. For this we use a region growing algorithm: we start from a vertex whose curvature is inferior to 0 and build a connected region of vertices whose curvature is inferior to 0. We repeat this procedure until there are no more initial seeds.

Figure 3
Evolution of the number of folds along time in 50 simulations.

We can observe on Fig. 3 that the number of sulci is equal to 0 on the interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e094.jpg then increases quasi linearly on An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e095.jpg and reaches a sort of plateau on An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e096.jpg.

Moreover we propose a simple way to characterize the spatial stability of the folds along time. In other terms we demonstrate that the position of folds formed at different time instants remains relatively stable. We extract a map of the curvature An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e097.jpg at each time instant An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e098.jpg. And we define a thresholded map An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e099.jpg by An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e100.jpg. This map depends on the mesh on which it is defined so we interpolate it on the final mesh (which has the largest number of vertices). We note that we use a smoothed version of the meshes in order to avoid problems of interpolation. The smoothing of the folded meshes has been performed by using an iterative process that consists, at each iteration, to replace a node of the mesh by the mean of its neighbors. So the maps An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e101.jpg are defined on the same domain and we can compute an average map

equation image

Intuitively this quantity represents the proportion of the temporal interval during which a fold is present at each position An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e103.jpg. This yields the map shown on Fig. 4.

Figure 4
Stability of folds in time.

We can see on this figure that the average map An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e108.jpg is not uniform but has patterns. In other terms we can observe a certain stability of the folds along time. In particular there are parts of the initial sphere that never belong to a fold. We note that the maximum of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e109.jpg is not An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e110.jpg since there are no folds during the temporal interval An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e111.jpg which represent An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e112.jpg of the full temporal interval on which the simulation has been performed.

Influence of noise

In this part we investigate the influence of noise in the spatial position of the folds. In particular we aim at demonstrating that the reaction diffusion mechanism is able to produce reproducible folds at certain specific locations but can also engender variability at other locations. For this we simulate An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e113.jpg realizations of the folding process from different noisy initial conditions An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e114.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e115.jpg. We consider the curvature maps An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e116.jpg and the thresholded maps An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e117.jpg at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e118.jpg and interpolate them on the same smoothed mesh. Then we sum the binary thresholded maps in order to see areas of reproducibility:

equation image

On Fig. 5 left, we can clearly observe that the sum of the binary maps representing the averaged pattern of folding have a spatial structure and do not organize randomly. In particular we notice a big longitudinal fold that comes across the surface and seems to be very reproducible among the 50 simulations. Moreover we note that other smaller reproducible folds are positioned along the main fold on both sides. This figure echoes the average cortical surface of 222 hemispheres that we have displayed on Fig. 5 right. This surface is the one described in [34] and has been visualized with anatomist [35]. The white arrow represents the cingulate sulcus that is comparable to our main fold in the 50 simulations while the three white stars show secondary folds that are parts of the para-cingulate sulcus and that we can link to the smaller reproducible folds of our model.

Figure 5
Reproducibility of the model compared to real data.

On Fig. 6 we illustrate however the variability of the main fold through three different scenarios of buckling. On the first line, left figure, we can see a mode that follows the main distribution previously described on Fig. 5 that is, in which the main fold is in one part. On the right is shown a left hemisphere of a real brain displayed with anatomist [35] on which the superior temporal sulcus (pink) is in one part. On the two following lines we represent two other modes for the main fold, in two and three parts respectively, and their correspondence on real anatomies with a superior temporal sulcus (pink) in two and three parts respectively.

Figure 6
Three modes of variability and their correspondence on real anatomies.

More generally Fig. 7 shows the different modes of variability of the main fold among the 50 simulations. In An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e121.jpg of cases it is composed of one segment, in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e122.jpg of two segments, in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e123.jpg of three segments and in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e124.jpg of four and five segments. The determination of connected components has been done by visual inspection rather than by automatic methods that tend to increase artificially the number of segments in the main fold.

Figure 7
Different modes of the main fold.

Phase diagram and pathologies of gyrification

It is possible to represent directly the influence of one or several parameters of the model (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e125.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e126.jpg in our case) on the qualitative properties of the patterns. We vary the parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e127.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e128.jpg linearly over a spatial domain from An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e129.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e130.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e131.jpg to An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e132.jpg respectively with steps of An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e133.jpg which yields 28 couples An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e134.jpg. On Fig. 8 we display the surfaces obtained at time An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e135.jpg. The color represents the curvature of these surfaces (red: positive curvature, blue: negative curvature). Note that the star-like patterns obtained with values An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e136.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e137.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e138.jpg, An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e139.jpg are just an artifact corresponding to the structure of the spherical mesh.

Figure 8
Phase diagram of the Gray-Scott model.

As suggested by [14] it is possible to link the qualitative nature of the obtained patterns to different modes of brain development, i.e in particular to pathologies or anomalies (see Fig. 9).

Figure 9
Patterns of folding and pathologies of development.

So we can see that for An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e143.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e144.jpg no patterns emerge. This state appears similar to lissencephaly, a pathology in which the brain is smooth and offers no gyri or sulci. The values An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e145.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e146.jpg might correspond to a normal brain with stripe-like patterns of gyrification. At last the range An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e147.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e148.jpg show spot-like patterns which make one think of polymicrogyria.

Discussion

Our model extends the initial proposal of Cartwright in [14] where no geometrical deformation of the cortical surface was proposed. We have demonstrated that it was possible to combine a reaction diffusion mechanism to a surface deformation in order to produce a model of the gyrification process. This approach is not new since it has been applied to model plant growth [21] but it seems to be the first to tackle the very old and controversial problem of brain folding in terms of reaction diffusion coupled to surface deformation. However the question about the origin of the morphogens used in our model remains open. In [14] the activation/inhibition process is supposed to model the mechanical tensions due to white matter fibers so the morphogenetic approach becomes indirect and extrinsic. On the contrary we prefer to view the folding process as the result of an intrinsic phenomenon, promoted by morphogens that decide the cytoarchitechtony. Different cytoarchitechtonic areas would correspond to different gyri and the limits between areas to sulci. This idea, suggested one century ago by Broadmann, has been recently pointed out in [36]. Moreover in [15], the GIP model supposes that the morphogens responsible for the patterning of subventricular zone could be some specific genes such as Pax6, Ngn2, Id4. Our model supports this hypothesis since mutations in the Pax6 gene for instance can be responsible for polymicrogyria [37], so the parameters An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e149.jpg and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e150.jpg of the model could reflect different gene expression of Pax6. We can also mention an alternative scenario for pattern formation that has been recently exposed in [38] and does not necessarily require the interaction of a long range inhibitor (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e151.jpg) and a short range activator (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e152.jpg) as in our case. In [38] an activator-activator mechanism combined with domain growth can also lead to pattern formation.

In our model we investigate also the variability of folding along the development of one individual and across several individuals - that is several realizations of the model. First we can see that for an unique development the position of the sulci remains stable along time. This result may seem trivial but is required for our model to produce definite patterns of gyrification that can be compared between different realizations of the model. Secondly the study of folding variability among 50 random realizations of the model reveals two important characteristics. The folding does not organize randomly even if we add noise to the initial condition of the reaction diffusion process. We have shown on one example that a main structure emerges that is strongly reproducible among several simulations. We can find a direct analogy between this main fold and the primary folds described in the literature [1], [39]. Primary folds are indeed characterized by their early time of appearance and their reproducibility across subjects. If we follow the comparison we can link the smaller structures found on Fig. 5 to the secondary or tertiary folds that are more posterior and variable than the primary ones. Our average map on Fig. 5 left echoes the average cortical surface of 222 hemispheres displayed on Fig. 5 right. This average surface has been computed in [34] and used also to represent an average map of sulcal pits density in [17]. In particular the main fold found on our simulations can be compared to the cingulate sulcus while the smaller folds around evoke the small pits of the paracingulate sulcus.

Moreover we have shown that in spite of its strong reproducibility the main fold could be broken in two separate parts by a gyrus. This result echoes previous studies [39], [40] where it is shown that some primary sulci reveal variability in their topology. For instance in [40] the superior temporal sulcus is continuous in one third of the cases (An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e153.jpg on the left and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e154.jpg on the right), in two segments in An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e155.jpg on the right and An external file that holds a picture, illustration, etc.
Object name is pcbi.1000749.e156.jpg on the left. The gyri that separate our main structure in two or three parts could also be interpreted in terms of ‘pli de passage’, which is a fold that can divide a sulcus in two sulci or just be buried at the bottom of a sulcus [2].

On a more theoretical point of view, our results on the reproducibility of the folds seem to confirm the impact of growth domain on the robust selection of patterns as it has been previously shown in [41] [42]. In our study there remains however some points that will require some theoretical developments, in particular about the existence of Turing instabilities that occur in the simulations. Some results have been obtained recently for isotropic domain growth or specific growth function [43] [38].

In conclusion we have proposed an extended framework for modelling the cortical folding. It is based on a system of coupled reaction-diffusion equations defined on a surface that evolves through the action of morphogens. We show that for some parameters the model gives rise to geometric patterns that can be related to cortical sulci. We also demonstrate that under the effect of noise the system yields morphological variability in these cortical structures. Moreover changing slightly the values of the parameters of the model can have an important influence on the nature of the created patterns which suggest a link toward pathologies of the brain development such as lissencephaly or polymicrogyria. In future developments we plan to investigate the difficult issue of estimating good values of parameters with respect to a given sequence of cortical surfaces across development.

Acknowledgments

We would like to thank Oliver Lyttelton, Maxime Boucher, Steven Robbins, and Alan Evans for having given us the possibility to use the average surface of 222 hemispheres.

We also would like to thank Roberto Toro for his comments on the manuscript.

Footnotes

The authors have declared that no competing interests exist.

The authors received no specific funding for this article.

References

1. Welker W. Why does cerebral cortex fissure and fold? A review of determinants of gyri and sulci. Cerebral Cortex. 1990;8:3–136.
2. Régis J, Mangin J, Ochiai T, Frouin V, Riviére D, et al. “Sulcal Root” Generic Model: a Hypothesis to Overcome the Variability of the Human Cortex Folding Patterns. Neurologia medico-chirurgica. 2005;45:1–17. [PubMed]
3. Biondi A. Are the brains of monozygotic twins similar? A three-dimensional MR study. American Journal of Neuroradiology. 1998;19:1361–1367. [PubMed]
4. Davidson L. How do sea urchins invaginate? Using biomechanics to distinguish between mechanisms of primary invagination. Development. 1995;121:2005–2018. [PubMed]
5. Le Gros Clark W. Deformation patterns in the cerebral cortex. Essays on Growth and Form. 1945:1–22.
6. Barron D. An experimental analysis of some factors involved in the development of the fissure pattern of the cerebral cortex. J Exp Zool. 1950;113:553–581.
7. Van Essen D. A tension-based theory of morphogenesis and compact wiring in the central nervous system. Nature. 1997;385:313–8. [PubMed]
8. Geng G, Johnston L, Yan E, Britto J, Smith D, et al. Biomechanisms for modelling cerebral cortical folding. Medical Image Analysis. 2009;13:920–930. [PubMed]
9. Richman D, Stewart R, Hutchinson J, Caviness V., Jr Mechanical Model of Brain Convolutional Development. Science. 1975;189:18–21. [PubMed]
10. Raghavan R, Lawton W, Ranjan S, Viswanathan R. A continuum mechanics-based model for cortical growth. Journal of Theoretical Biology. 1997;187:285–296.
11. Toro R, Burnod Y. A Morphogenetic Model for the Development of Cortical Convolutions. Cerebral Cortex. 2005;15:1900–1913. [PubMed]
12. Rakic P. Neuroscience: genetic control of cortical convolutions. Science Signaling. 2004;303:1983. [PubMed]
13. Rakic P. Specification of cerebral cortical areas. Science. 1988;241:170–176. [PubMed]
14. Cartwright J. Labyrinthine Turing Pattern Formation in the Cerebral Cortex. Journal of Theoretical Biology. 2002;217:97–103. [PubMed]
15. Striegel D, Hurdal M. Chemically Based Mathematical Model for Development of Cerebral Cortical Folding Patterns. PLoS Comput Biol. 2009;5 [PMC free article] [PubMed]
16. Lohmann G, von Cramon D, Colchester A. Deep sulcal landmarks provide an organizing framework for human cortical folding. Cerebral Cortex. 2008;18:1415–1420. [PubMed]
17. Im K, Jo H, Mangin J, Evans A, Kim S, et al. Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface. Cerebral Cortex 2009 [PubMed]
18. Tessier-Lavigne M, Goodman C. The molecular biology of axon guidance. Science. 1996;274:1123. [PubMed]
19. Kriegstein A, Noctor S, Martnez-Cerdeño V. Patterns of neural stem and progenitor cell division may underlie evolutionary cortical expansion. Nature Reviews Neuroscience. 2006;7:883–890. [PubMed]
20. Barrio R, Varea C, Aragón J, Maini P. A two-dimensional numerical study of spatial pattern formation in interacting Turing systems. Bulletin of mathematical biology. 1999;61:483–505. [PubMed]
21. Harrison L, Wehner S, Holloway D. Complex morphogenesis of surfaces: theory and experiment on coupling of reaction-diffusion patterning to growth. Faraday Discussions. 2002;120:277–293. [PubMed]
22. Holloway D, Harrison L. Pattern selection in plants: Coupling chemical dynamics to surface growth in three dimensions. Annals of Botany. 2008;101:361. [PMC free article] [PubMed]
23. Turk G. Generating textures on arbitrary surfaces using reaction-diffusion. Proceedings of the 18th annual conference on Computer graphics and interactive techniques. 1991:289–298.
24. Witkin A, Kass M. Reaction-diffusion textures. ACM SIGGRAPH Computer Graphics. 1991;25:299–308.
25. Pearson J. Complex Patterns in a Simple System. Science. 1993;261:189–192. [PubMed]
26. McGough J, Riley K. Pattern formation in the gray-scott model. Nonlinear Analysis: Real world Applications. 2004;5:105–121.
27. Crampin E, Hackborn W, Maini P. Pattern formation in reaction-diffusion models with nonuniform domain growth. Bulletin of Mathematical Biology. 2002;64:747–769. [PubMed]
28. Plaza R, Sánchez-Garduño F, Padilla P, Barrio R, Maini P. The Effect of Growth and Curvature on Pattern Formation. Journal of Dynamics and Differential Equations. 2004;16:1093–1121.
29. Chung M, Taylor J. Diffusion smoothing on brain surface via finite element method. 2004. In: Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI)
30. Pennacchio M, Simoncini V. Efficient algebraic solution of reaction-diffusion systems for the cardiac excitation process. Journal of Computational and Applied Mathematics. 2002;145:49–70.
31. Delaunoy A, Fundana K, Prados E, Heyden A. Convex multi-region segmentation on manifolds. 2009. Proceedings of the 12th International Conference on Computer Vision, Kyoto, Japan.
32. Dziuk G, Elliott C. Finite elements on evolving surfaces. IMA Journal of Numerical Analysis. 2007;27:262.
33. Cachia A, Mangin J, Riviere D, Kherif F, Boddaert N, et al. A primal sketch of the cortex mean curvature: a morphogenesis based approach to study the variability of the folding patterns. IEEE transactions on medical imaging. 2003;22:754–765. [PubMed]
34. Lyttelton O, Boucher M, Robbins S, Evans A. An unbiased iterative group registration template for cortical surface analysis. Neuroimage. 2007;34:1535–1544. [PubMed]
36. Fischl B, Rajendran N, Busa E, Augustinack J, Hinds O, et al. Cortical folding patterns and predicting cytoarchitecture. Cerebral Cortex. 2007;18:1973–80. [PubMed]
37. Mitchell T, Free S, Williamson K, Stevens J, Churchill A, et al. Polymicrogyria and absence of pineal gland due to PAX6 mutation. Annals of neurology. 2003;53:658–662. [PubMed]
38. Madzvamuse A, Gaffney E, Maini P. Stability analysis of non-autonomous reaction-diffusion systems: the effects of growing domains. Journal of Mathematical Biology. 2009:1–32. [PubMed]
39. Ono M, Kubik S, Abernathey C. Atlas of the cerebral sulci. Georg Thieme Verlag 1990
40. Ochiai T, Grimault S, Scavarda D, Roch G, Hori T, et al. Sulcal pattern and morphology of the superior temporal sulcus. Neuroimage. 2004;22:706–719. [PubMed]
41. Maini P, Crampin E, Madzvamuse A, Wathen A, Thomas R. Implications of domain growth in morphogenesis. Mathematical modelling & computing in biology and medicine. 2003:67.
42. Barrass I, Crampin E, Maini P. Mode Transitions in a Model Reaction–Diffusion System Driven by Domain Growth and Noise. Bulletin of mathematical biology. 2006;68:981–995. [PubMed]
43. Madzvamuse A. Stability analysis of Reaction-Diffusion Systems with constant coefficients on growing domains. International Journal of Dynamical Systems and Differential Equations. 2008;1:250–262.

Articles from PLoS Computational Biology are provided here courtesy of Public Library of Science