Functional optimality of the sulcus pattern of the human brain

We develop a mathematical model of information transmission across the biological neural network of the human brain. The overall function of the brain consists of the emergent processes resulting from the spread of information through the neural network. The capacity of the brain is therefore related to the rate at which it can transmit information through the neural network. The particular transmission model under consideration allows for information to be transmitted along multiple paths between points of the cortex. The resulting transmission rates are governed by potential theory. According to this theory, the brain has preferred and quantized transmission modes that correspond to eigenfunctions of the classical Steklov eigenvalue problem, with the reciprocal eigenvalues quantifying the corresponding transmission rates. We take the model as a basis for testing the hypothesis that the sulcus pattern of the human brain has evolved to maximize the rate of transmission of information between points in the cerebral cortex. We show that the introduction of sulci, or cuts, in an otherwise smooth domain indeed increases the overall transmission rate. We demonstrate this result by means of numerical experiments concerned with a spherical domain with a varying number of slits on its surface.


Introduction
The complex convoluted structure of the cerebral cortex is a hallmark of many mammalian brains and, in particular, of the human brain. Its specific surface morphology and the associated mechanisms underlying its growth remain the subject of active topical research. Open challenges range from identifying the triggers for cortical folding processes to understanding the correlation between brain structure and brain function.
Mammalian brains are composed of an outer cortical layer, which is referred to as the grey matter and which primarily consists of neuronal and glial cell bodies. In contrast, the inner subcortical core acts as a transmission structure and comprises axons connecting the neurons on the cerebral cortex. Such axons are coated with the electrically insulating substance myelin, which is white in appearance and lends the subcortical core the name white matter. Axons typically show complex tree structures with myriads of branch points, and they make connections with other cells at junctions called synapses. The profuse branching of axons results in a total of approximately 0.15 quadrillion synapses, an exceedingly large number when compared to an estimated 86 billion neurons in the average adult male human brain (Hormuzdi et al., 2004). Due to the geometrical constraint imposed by the cranium, within which the brain must be embedded, the formation of cortical folds is the only way of increasing the cerebral surface area, which in turn results in an increase in the total number of neurons and a simultaneous minimization of their relative separation distance (Zilles et al., 2013). Since the total number of neurons and axons as well as the length of connections directly correlates to signaling speed and total information transmission rate, the cerebral surface area is viewed as a strong indicator of intelligence (Roth and Dicke, 2005). Furthermore, several studies spanning diverse neurological disorders have revealed that the degree of cortical folding can serve as an indicator of abberations in brain development. For instance, reduced gyrification is found to be an inherent feature of the brain pathology in schizophrenia (Nesvåg et al., 2014), whereas increased frontal cortical folding is observed in the context of autism (Hardan et al., 2004).
The particular mechanisms that drive the folding process during growth are as yet not fully understood. Van Essen (Van Essen, 1997) first introduced the axonal tension hypothesis which assumes that tension along axons in the human connectome is the primary driving force for cortical folding. Based on this theory, axonal tension acts against hydrostatic pressure that is generated internally by the cerebrospinal fluid. Since long-distance connections enter and leave the cortex exclusively through the subcortical white matter, strongly interconnected regions will be pulled together during brain development due to tension along axons, whereas weakly connected regions are allowed to drift apart. Neighboring areas that are weakly interconnected are thus separated by outward folds, or gyri, whereas strongly interconnected areas are separated by inward folds, or sulci. As an alternative theory, the grey matter hypothesis postulates growth processes during cortical development to be governed by the driving forces of differential growth (Richman et al., 1975). Based on this model, folding processes are predicted when the growth of the outer (supragranular) cortical layer exceeds that of the inner (infragranular) layer (Zilles et al., 2013). Differential growth is thus introduced into the theory as a mechanism to release residual mechanical stresses by allowing for surface buckling. More recent studies combine the hypotheses of tension-mediated and differential-growth-induced cortical folding by modeling the cortex as a morphogenetically growing outer layer in combination with the subcortex as a strain-driven growing inner core (Budday et al., 2014). Both theories of differential growth and axonal tension, however, disagree with experimental findings. While dissection experiments did not reveal significant tangential tension within developing gyri as postulated by the axonal tension hypothesis (Xu et al., 2010), the differential growth hypothesis relies on unrealistic differences in stiffness parameters between the cortex and the underlying subcortical layers (Bayly et al., 2014).
Rather than how the human brain cortex folds during growth, in this work we address the question of why it is folded in the first place. We suppose that the overall function of the brain consists of the emergent processes resulting from the spread of information through the neural network. The capacity of the brain is therefore related to the rate at which it can transmit information through the neural network. The hypothesis underlying the present work is that the sulcus pattern of the human brain has evolved to maximize the rate of transmission of information between points in the cerebral cortex.
In order to render the problem in tractable mathematical terms, we formulate a simple graph model that relates the rate of transmission between points of the cortex to a problem in potential theory (cf. (Avena-Koenigsberger et al., 2014;Misic et al., 2015) for the biological basis of this model) and the classical Steklov eigenvalue problem (Steklov, 1902). The particular transmission model under consideration allows for information to be transmitted along multiple paths between points of the cortex. We regard the human brain as a biological neural network, i. e., a collection of interconnected neurons. The interface between neurons consists of several axon terminals connected via synapses to dendrites. At any time, a neuron in the network can have an activation in the form of an action potential spike. The activation then spreads to all other connected neurons, which in turn become activated. From an information-theoretical point of view, each neuronal activation may be regarded as a bit of information. The overall function of the brain then consists of the emergent processes resulting from the spread of information through the neural network. The capacity of the brain is therefore related to the rate at which it can transmit information through the neural network.
A remarkable outcome of this model is that the transmission of information within the brain is quantized. Thus, we find that the brain has preferred transmission modes that correspond to eigenfunctions of the classical Steklov eigenvalue problem (Steklov, 1902), with the reciprocal eigenvalues quantifying the corresponding transmission rates. The Steklov spectrum of the brain thus collects all the preferred modes of transmission, or eigenfunctions, of the brain. The Steklov eigenvalue problem originally arose in connection with hydrodynamics and has been extensively studied (cf. (Girouard and Polterovich, 2017) for a recent review). The Steklov spectrum is discrete and the reciprocal Steklov eigenvalues, or transmission rates, have an accumulation point at zero (Moiseev, 1964;Kopachevskii and Krein, 2001;Brock, 2001).
The model may be taken as a basis for testing the hypothesis that the sulcus pattern of the human brain has evolved to maximize the rate of transmission of information between points on the cerebral cortex. Specifically, the question is whether the sulcus pattern may be understood as the shape that minimizes the Steklov eigenvalues among all domains contained within a fixed confining set (the cranium). This type of shape optimization differs somewhat from the classical Steklov shape optimization problem which is concerned with competitor domains of fixed measure (cf., e. g., (Bogosel et al., 2017;Girouard and Polterovich, 2017)). While the full shape optimization problem is beyond the scope of this paper, we take some steps in that direction. In particular, we show that the introduction of sulci, or cuts, in an otherwise smooth domain indeed increases the overall transmission rate. We additionally demonstrate this result by means of numerical experiments concerned with spherical domains with a varying number of slits on its surface.

Formulation of the problem
In order to formulate a mathematical model of information transmission in the brain, we consider a simple neural network in the form of a cubic lattice. 1 The nodes of the lattice represent the neurons and the bonds the synapses. Other lattices, including random networks, can be treated likewise without essential change in the outcome. Let aZ 3 denote the lattice of points x = (al 1 , al 2 , al 3 ), with a > 0 the lattice parameter and l 1 , l 2 , l 3 ∈ Z, where here and subsequently Z denotes the set of integer numbers. Two points x, y ∈ aZ 3 are nearest neighbors if |x − y| = a. A path in aZ 3 is a sequence of points such that every consecutive pair is also a pair of neighboring points.
Next, we turn to the propagation of signals through the lattice. Consider a signal that starts at x and subsequently traverses a path in the lattice. At every point along the path, the signal has the choice of moving to one of the 6 neighbors of the point with probability 1/6 (cf. (Avena-Koenigsberger et al., 2014;Misic et al., 2015) for the biological basis of this model). By these set of rules, the paths available to the lattice for the transmission of information are random walks. If we define the one-step transition probability as (1) p(x, y) = 1/6, if x, y neighbors, 0, otherwise.
Then, by the Markov property of random walks, is the probability that the signal traverse a path Γ = {x, z 1 , z 2 , . . . , z k−1 , y} joining x to y. A frequentist interpretation of these probabilities is that signals have a choice of paths to travel between points of the neural lattice and that (2) gives the frequency with which a signal originating at x traverses a particular path Γ to reach another point y.
There is a well-known connection between random walks and harmonic functions (cf., e. g., (Treves, 1970)). We recall that a function u : We follow the notation of (Treves, 1970), which may be consulted for background on the connection between Brownian motion and potential theory.
i. e., if its value at every point x of the lattice equals its average over the neighbors of x. The one-step shift or averaging operator u → P u is defined as and the discrete Laplacian as where I is the identity. From these definitions, it follows that a function u is harmonic if and only if i. e., if it is a solution of Laplace's equation. By the maximum principle of harmonic functions, it follows that any bounded harmonic function over the entire lattice aZ 3 is necessarily constant. Consider now the discrete Poisson equation where f is a distribution of sources over aZ 3 . Since, as already noted, the kernel of the discrete Laplacian consists of constant functions, for solutions to exist the Fredholm alternative requires f to sum to zero, in which case solutions are determined up to an additive constant. Inserting (5) into (7) and solving for u we obtain is the free-space discrete Green's function. By translation-invariance, we have G(x, y) = F (x−y), where F is the fundamental solution of the discrete Laplacian.
The sought connection between potential theory and transmission of information along lattice paths can now be forged as follows. Expanding (9) in Neumann series gives With P as in (4), it is readily shown that this series converges uniformly. Evaluating the powers P k explicitly, we obtain where is the discrete Dirac function. For x = y, it follows from (1) and (4) that the term P (x, y) in (11) contributes to the sum only if y is a neighbor of x. Likewise, the product p(x, z 1 )p(z 1 , z 2 ) · p(z k−1 , y) is non-zero only if the points {x, z 1 , z 2 , . . . , z k−1 , y} define a path Γ joining x to y. Thus, eq. (11) can be recast in the revealing form where p(Γ) is defined in (2) and P(x, y) denotes the set of all paths joining x to y in the lattice. The path-sum representation (13) shows that the Green's function G(x, y) of the discrete Laplacian is the sum of contributions p(Γ) arising from all paths joining x to y.
In the neural network representation of the brain, we may regard p(Γ), eq. (2), as the probability that a signal originating at x reach y through the path Γ ∈ P(x, y). From representation (13), it then follows that G(x, y) is the total rate at which information injected into the network at x is transmitted to y. This identification establishes the sought connection between the transmission of information through the neural network and discrete potential theory. Since the size of the neurons is much smaller than the overall size of the brain, we may further expect the continuum limit a → 0 to supply a good approximation. This continuum limit can indeed be effected rigorously (Treves, 1970). However the analysis is technical and beyond the scope of the present work. Therefore, in the sequel we proceed formally and simply replace the preceding discrete framework by its formally equivalent continuum counterpart.

Connection with the Steklov eigenvalue problem
A connection between the preceding graph model and the classical Steklov eigenvalue problem can be forged as follows. We consider a neural network occupying a domain Ω. We further consider a distribution h ∈ H −1/2 (∂Ω) of signals exchanged between points of the cortex ∂Ω, with In the continuum limit, the corresponding transmission rate is We suppose that the brain has preferred transmission modes that maximize R(h) locally subject to the zero-sum condition (14) and the normalization constraint We can combine the maximum transmission rate objective and the constraints into the Lagrangian where λ and µ are Lagrange multipliers, to be maximized with respect to u ∈ H 1 (Ω) and h ∈ H −1/2 (∂Ω). The stationarity of F demands that From the zero-sum condition (14), we find where |∂Ω| denotes the area of ∂Ω. Inserting into (19) and solving for h(x), we obtain We thus conclude that the transmission mode h(x) admits the representation From the normalization constraint (17), we additionally find whereupon (22) becomes This expression yields a general representation of the transmission modes h(x) in terms of potentials v(x) with zero mean trace (23) over the boundary. Inserting (25) into (18) and using (23), we obtain the reduced functional to be maximized locally with respect to v ∈ H 1 (Ω) subject to the zero mean trace constraint (23). The corresponding Euler-Lagrange equations are where n is the outward unit normal and we write .

Monotonicity with respect to cutting
The Steklov eigenvalues have a number of domain monotonicity properties that bear directly on the present discussion (cf., e. g., Kulczycki and Kuznetsov (2009)). We recall that the Steklov eigenvalues admit the variational characterization (cf., e. g., Lamberti and Provenzano (2013); Bogosel et al. (2017); Girouard and Polterovich (2017)) (30) σ n = inf where the infimum is taken over all n-dimensional subspaces V n of H 1 (Ω)/R = {u ∈ H 1 (Ω),´∂ Ω u(x) dS(x) = 0}. The corresponding eigenfunctions v n represent transmission modes of the brain and the inverse eigenvalues 1/σ n give the corresponding transmission rates. The beneficial effect of sulci to brain function can be deduced from (30) as follows. Let Ω be the domain of the cranium and let {σ n } be its Steklov spectrum. Let Γ be a collection of cuts (sulci) performed on the boundary Ω, set Ω = Ω\Γ, and let {σ n } be the Steklov spectrum of the slit domain Ω . From the variational characterization (30) of the Steklov eigenvalues, we have since H 1 (Ω) ⊂ H 1 (Ω ). We thus conclude that the introduction of sulci increases the transmission rate of all the transmission modes of the brain.
The same argument shows that the transmission rate also increases when existing cuts are made deeper. These monotonicity properties explain the functional benefit of sulci to the function of the brain.

Numerical experiments
We further illustrate the connection between sulci and rate of transmission by means of selected numerical experiments. Specifically, we compare the Steklov eigenvalues of a sphere with those of slit-spheroidal domains and show that the introduction of slits does indeed reduce the Steklov eigenvalues. We carry out all calculations be recourse to the finite element method. The accuracy of the method is assessed and controlled with the aid of the known Steklov spectrum of the spherical domain.
5.1. Spherical domain. In order to compute the Steklov spectrum it is convenient to reformulate the problem as an eigenvalue problem for the Laplace operator with Neumann boundary conditions, with mass density ρ concentrated on the boundary ∂Ω. The eigenvalues of problem (32) are known to coincide with those of the classical Steklov eigenvalue problem (Lamberti and Provenzano, 2013;Arrieta et al., 2008).
eigenvalues have multiplicity We assess the accuracy of the finite element discretization by means of the unit ball test case just described.We discretize Ω by means of ≈ 11, 000 linear tetrahedral elements and tile the boundary Ω by means of ≈ 2, 000 membrane elements of uniform areal mass density. The interior degrees of freedom are subsequently eliminated by means of static condensation, resulting in an eigenvalue problem involving the surface degrees of freedom only. The 9 lowest eigenvalues of the computed Steklov spectrum are compared in Table 1 to the exact analytical values. As can be seen from the table, the lowest non-zero eigenvalue is computed with a relative error of 0.3 %. In general, the numerical accuracy of eigenvalues σ j is known to decrease with increasing j. For instance the eigenvalue σ 25 = 5 is obtained with a relative error of 1.0 %. We further note that numerical eigenvalues have the expected multiplicities. 5.2. Slit-spheroidal domains. Next we proceed to compute the Steklov spectrum of slit-spheroidal domains. To this end, we start with a spherical domain of radius 6.7 cm, which matches the average male human brain volume of 1260 cm 3 (Cosgrove et al., 2007). Subsequently, slit-spheroidal domains are constructed by successively introducing slits of varying depth, intended to represent idealized sulci. Geometry definition as well as meshing are performed using the Siemens PLM Software Femap with NX Nastran. Representative geometries and finite element meshes are shown in Fig. 1. Fig. 2 shows the lowest 25 Steklov eigenvalues σ j of slit-spheroidal geometries with up to 5 slits. The finite element meshes used in calculations contain ≈ 7, 000 linear tetrahedral elements and ≈ 5, 000 membrane elements. As may be seen from the figure, the eigenvalues are found to decrease for an increasing number of slits.
To further elucidate the influence of the slit depth, Fig. 3 shows the Steklov spectra of the lowest 25 eigenvalues σ j of a one-slit-spheroidal geometry for an increasing depth of the slit while keeping the width of the slit constant. Meshes with a slit depth of 0.67 cm as well as 2.68 cm are generated and the calculated eigenvalues are depicted in comparison to corresponding values for the smooth domain. As may be seen from the figure, all eigenvalues decrease with increasing slit depth.
These results are in keeping with the monotonicity arguments of Section 4, which show that the Steklov eigenvalues (respec., transmission rates) decrease (respec., increase) when cuts are introduced in a domain and when the cuts are made deeper.

Concluding remarks
This paper establishes a connection, through a neural network model, between information transmission between points in the brain cortex and the classical Steklov eigenvalue problem. The implications of this connection are manifold. Firstly, according to the model the modes of transmission of the  Figure 3. Steklov spectra of the lowest 25 eigenvalues σ j of a one-slit-spheroidal geometry for an increasing depth of the slit.
brain are described spatially by the Steklov eigenfunctions and, therefore, are quantized, with the corresponding eigenvalues supplying the transmission rates of the modes. In particular, the lowest Steklov eigenvalue gives the largest possible rate of transmission of the brain. However, the transmission rates 1/σ j have an accumulation point at the origin, eq. (29), with the result that slow transmission rates are densely distributed and quantization is ostensibly lost in that range. The physiological significance of the fast transmission modes, corresponding to the lowest Steklov eigenvalues, is yet to be elucidated. We have also shown that the transmission efficiency of a domain increases when cuts are introduced. Specifically, the Steklov eigenvalues decrease upon the introduction of the cuts and decrease further as the cuts are made deeper. This type of monotonicity provides an explanation for the functional benefit of the sulci, regarded as cuts in an otherwise smooth brain. However, since the transmission efficiency increases monotonically with the depth of the cuts, no regular brain shape can be expected to maximize transmission efficiency or, equivalently, minimize the Steklov eigenvalues. A well-posed optimization problem can be obtained by introducing additional constraints, e. g., on the cortical surface area. Indeed, the classical shape optimization problems for the Steklov eigenvalues consider competing domains of fixed boundary area (cf., e. g., Bogosel et al. (2017); Girouard and Polterovich (2017)). Such a strong surface area constraint is not appropriate to the case at hand, in which the brain shapes of interest fill the cranium and differ only in the sulcus, or cut, pattern. An alternative is to penalize, instead of constraining, the cortical surface area. An objective function of this type is, for instance, for some constants A and α and with |∂Ω| the surface area. The additional term A |∂Ω| α represents the physiological cost of cortical surface area. The problem is then to minimize F over all brain domains Ω contained in the cranium K. Our conjecture is that the optimal shape is obtained by cutting K into an optimal sulcus pattern, but a rigorous analysis of this conjecture is beyond the scope of the present work.
Objective functions of the form (35) can also potentially explain scaling relations and size effects that have been uncovered through systematic analyses of the variation in cortical folding across large samples of mammalian species. For instance, Mota and Herculano-Houzel (2017) have shown that the degree of cortical folding scales uniformly across species, across individuals, and within individual cortices. Specifically, for all noncetacean gyrencephalic species they find the scaling relation Remarkably, this power law is significantly superlinear, which implies that, as the total surface area increases, the brain becomes increasingly folded. We recall that the ratio |∂Ω|/|∂K| of total cortical area to exposed cortical area is related to the classical gyrification index (Zilles et al., 1988). The scaling law (36) shows that the gyrification index is a function of brain size, with larger brains having a more folded cortex. Scaling relations such as (36) can be rationalized within the present framework as follows. Suppose that the fundamental transmission rate 1/σ 1 (Ω) obeys an optimal scaling law of the form for some constants C > c > 0 and exponents and δ. We recall that scaling laws are said to be optimal if they entail power-law lower and upper bounds with matching exponents, in this case and δ. Optimal scaling laws were developed in mathematics in connection with energy-minimizing branched microstructures Müller, 1992, 1994;Choksi et al., 1999;Conti, 2000). From dimensional considerations, we must have Suppose, in addition, that the optimal brain shape satisfies equipartition between the two terms of the objective function (35). Then, we have Inserting into (37) and assuming that the upper and lower bounds are tight, we have Comparing with (36) we finally find which gives the exponent β in terms of the optimal-scaling exponents and δ and the cost exponent α. If, for instance, we assume that the cortical surface area cost is proportional to the cortical surface areas, α = 1, then (42) specializes to where we have used (38). Identity (42) gives the gyrification exponent β in terms of the optimal scaling exponents for the fundamental Steklov eigenvalue.
Beyond the question of gyrification, we remark that the Steklov eigenfunctions provide a convenient orthogonal basis for the spatial representation of brain activity. Thus, for a given patient-specific brain geometry Ω the Steklov eigenfunctions and eigenvalues can be computed numerically, e. g., by recourse to the finite-element method as in Section 5 or by other means (Akhmetgaliyev et al., 2017). We also recall that, according to the theory put forth in Section 3, the brain activity is described by a potential u over Ω. It thus follows that any pattern of brain activity u can be decomposed into, and then represented as a sum of, Steklov modal components.
We also note that the gradient of the potential, J = ∇u, measures the electrical current density associated with the activity of the brain. By Ampère's law, ∇ × B = µ 0 J, this current density induces a magnetic field B = ∇ × A, where A is a vector potential satisfying the gage condition ∇ · A = 0. This vector potential, and the corresponding magnetic field, are given by the Biot-Savart law and extend outside the brain and the cranium. Though weak, the exterior fields can be detected and measured, e. g., with the aid of sensitive detectors known as superconducting quantum interference devices (SQUIDs). Such measurements are known as magnetoencephalograms, or MEGs (cf., e. g., Cohen and Halgren (2004)).
A compelling alternative consists of reversing the process and inducing activity patterns in the brain through the application of external magnetic fields, a process known as transcranial magnetic stimulation (TMS) (cf., e. g., George and Belmaker (2000)). During a TMS procedure, a magnetic field generator in the form of a coil is placed near the head of the patient. Evidence suggests that TMS is effective against neuropathic pain and treatment-resistant major depressive disorder. However, the magnetic field generators, or coils, used at present are not tailored to patient-specific geometries and are relatively delocalized. The detailed knowledge, through patient-specific calculations, of the Steklov spectra of an individual brain opens up the way for stimulating specific transmission modes, including the fundamental modes at which the brain performs at its greatest capacity. Such mode-specific stimulation could be achieved by externally applying to the cranium shaped electromagnetic fields corresponding to the fundamental Steklov modes. The therapeutic benefits of such tailored TMS procedures are yet to be ascertained.