https://orcid.org/0000-0002-6502-7066
Abstract
Association fibers connect different areas of the cerebral cortex over long distances and integrate information to achieve higher brain functions, particularly in humans. Prototyped association fibers are developed to the respective tangential direction throughout the cerebral hemispheres along the deepest border of the subplate during the fetal period. However, how guidance to remote areas is achieved is not known. Because the subplate is located below the cortical surface, the tangential direction of the fibers may be biased by the curved surface geometry due to Sylvian fissure and cortical poles. The fiber length can be minimized if the tracts follow the shortest paths (geodesics) of the curved surface. Here, we propose and examine a theory that geodesics guide the tangential direction of long association fibers by analyzing how geodesics are spatially distributed on the fetal human brains. We found that the geodesics were dense on the saddle-shaped surface of the perisylvian region and sparse on the dome-shaped cortical poles. The geodesics corresponded with the arrangement of five typical association fibers, supporting the theory. Thus, the geodesic theory provides directional guidance information for wiring remote areas and suggests that long association fibers emerge from minimizing their tangential length in fetal brains.
Introduction
Association fibers achieve higher brain functions in humans and connect long-distance special anatomical areas of the cerebral cortex to integrate their functions (Mori et al. 2005; Philippi et al. 2009; Schmahmann et al. 2009; Catani et al. 2012; Egger et al. 2015). There are five major association fibers, which connect Broca’s area to Wernicke’s area (the arcuate fasciculus), the cingulate gyrus to the entorhinal cortex (the cingulum (cg)), the frontal lobe to the temporal lobe (the uncinate fasciculus (unc)), the frontal lobe to the occipital lobe (the inferior fronto-occipital fasciculus (ifo)), and the occipital lobe to the temporal lobe (the inferior longitudinal fasciculus (ilf)), based on Mori’s atlas (Mori et al. 2005). Association fibers begin to form early in the fetal period at around 13 gestational weeks (GW), and all major fibers emerge by 20 GW (Huang et al. 2009; Takahashi et al. 2012; Žunić Išasegi et al. 2018; Horgos et al. 2020). Although the wiring of interhemispheric commissural fiber on the coronal section is well understood at the molecular level, including the axon guidance (Stoeckli 2018; Comer et al. 2019), a method to wire the long association fibers in a tangential direction within the hemisphere has not yet been discovered, likely due to experimental difficulty in cutting out the respective brain sections.
Association fibers are initially formed at the deepest border of the subplate layer, which appears only during the fetal period directly below the cortical layer (Kostović et al. 2002; Kanold and,Luhmann 2010; Duque et al. 2016; Hadders-Algra 2018; Molnár et al. 2019; Vasung et al. 2019). The subplate layer is a tangentially continuous tissue throughout the cerebral hemispheres (Kostović et al. 2002) and is well-developed in primates, especially in humans (Duque et al. 2016). Subplate neurons are sparsely distributed in this layer and are surrounded by glial cells and extracellular matrix, making the subplate layer uniform and elastic (Kostović et al. 2019). Subplate neurons are locally connected to each other through gap junctions (Moore et al. 2014) and chemical synapses (Kristt and Molliver 1976; Luhmann and Khazipov 2018; Luhmann et al. 2022). The locally connected subplate neurons show the membrane potential activities spontaneously or stimulus-responsively (Friauf and Shatz 1991; Moore et al. 2011; Ohtaka-Maruyama et al. 2018; Molnár et al. 2020; Ohtaka-Maruyama 2020), which can be propagated tangentially throughout the layer (Garaschuk et al. 2000; Sun and Luhmann 2007; Kilb et al. 2011). Moreover, despite the immaturity of the cortical plate at approximately 24 GW, electroencephalogram (EEG) and functional magnetic resonance imaging (fMRI) findings have indicated the presence of global neuronal activities throughout the cerebral hemispheres, likely derived from the subplate (Keunen et al. 2017; Thomason et al. 2017). Therefore, the subplate nexus hypothesis has been proposed to explain the developmental mechanism of association fibers. According to this hypothesis, the subplate layer, as a uniform continuum tangentially stretched across the cerebral hemispheres, creates the long-distance, “weak” connectivity (precursor of association fibers) by locally connecting neuronal fibers during the fetal period (Kostović 2020).
Because the subplate in the fetal period is located just below the cortical surface (Kostović et al. 2002), fiber formation may be biased by the smoothly curved landscape in the form of domes at the brain poles with positive Gaussian curvature and saddles on the perisylvian region (appearing around 11–29 GW) with negative Gaussian curvature (Fig. 1, left panel) (Van Essen 2020). Fiber formation between a pair of distant positions of the cortex has a higher cost, as the structural materials of the neuronal fibers (neurofilaments, microtubules, etc.) must cover a longer tract (Kostović and Rakic 1990; Kanold and Luhmann 2010; Bullmore and Sporns 2012; Oldham et al. 2022). Therefore, wiring likely follows the shortest tract to minimize the costs (Fig. 1, middle panel). We further assume that the association fibers are initially formed to take the shortest distance on a 2D surface determined by a specific laminar of the cerebrum. This is based on the histological studies showing that the developing cerebrum during the midgestation has a multilaminated structure of cells and fibers (Kostović et al. 1989) and the association fibers grow at the interface between the subplate and the external sagittal stratum (Žunić Išasegi et al. 2018). These assumptions allow us to derive a hypothetical principle for the long-distance fiber arrangement, hereafter referred to as the geodesic theory, in which the shortest path (i.e. geodesic (Kobayashi and Nomizu 1963) of a curved surface landscape provides the tangential direction of the prototyped association fibers. In the present study, we used open-source surface data at the midthickness of the cortex as an approximation for those at the deepest border of the subplate, where the association fibers are likely to grow. The closed surfaces of different lamina should be geometrically similar when the cerebrum has the smooth multilaminated structure without numerous gyri. Thus, this should be a valid approximation for the analysis of the geometric property of the surface landscape, which is a constraint on the generation of association fibers.
Representative methods used in this study. At 27 GW, the cortical surface geometry was represented by triangular meshes (left), and a geodesic between a pair of randomly selected positions was calculated (white line: a single geodesic, center). Overwriting |$1,500$| geodesics covered the entire brain surface (right).
The direction of the geodesics generally depends on the surface curvature (Kobayashi and Nomizu 1963). Geodesics, each connecting two arbitrary positions on the surface, are uniformly distributed on a sphere with spatially uniform curvature and can be non-uniformly distributed on surfaces with non-uniform curvature (e.g. sparse at domes with positive curvature and dense at saddles with negative curvature). Similarly, the brain poles and fissures/sulci are dome- and saddle-shaped, respectively. However, the spatial distribution of geodesics on cortical surfaces lacks evidence (Fig. 1, right panel). Assuming that a single geodesic corresponds to a single neuronal fiber, a bundle of multiple geodesics may account for the association fibers. Comparing the spatial distribution of geodesics with those of long-association fibers could validate the geodesic theory.
An appropriate stage for this comparison is around 27 GW, because at this time, the subplate layer comprises a large volume and all major association fiber prototypes are accumulated below the subplate layer (Vasung et al. 2010, 2017; Kostović 2020). At this age, the Sylvian fissure impact the overall landscape of the smooth surface of the subplate layer and the cortical layer, prior to the gyrification of the primary and secondary sulci (Van Essen 2020). Hence, in this study, we analyzed the cortical surface geometry at 27 GW to validate the geodesic theory.
Materials and methods
Brain samples
We used the left brain surface data, generated at the midthickness of the cortex using the surface reconstruction method (Hill et al. 2010), of a preterm infant with 27 GW published by Garcia, Robinson, et al. (2018). The corpus callosum was cut out, and the cross-section was masked. The brain surface data in gii format downloaded from BALSA (https://balsa.wustl.edu/) were converted to obj format files using the GIfTI library of Matlab (MathWorks Inc., USA) to create triangular mesh data. To modify the edge length of the triangular mesh to a mean of 1.25 mm, we sequentially applied a re-meshing function (uniform method, iterations = 3, smoothing = 0.1) and a smoothing function (curvature-dominant method) in Houdini (version 18.5; SideFX Inc., Canada) to the mesh. Then, the mesh of the corpus callosum was removed using Houdini’s Blast geometry node, resulting in 9,998 triangular mesh faces.
Geometrical analysis
We calculated the shortest paths (geodesics) on the brain surface data using the eikonal equation:
where |$T$| is the distance function and |$F(x,y)$| is the velocity of movement of the boundary. In this study, we set |$F(x,y) = 1$|. To obtain a finite element approximation of the eikonal equation, we used the fast marching method (Sethian 1996; Kimmel and Sethian 1998), which provides a distance function on the cortical surface between the start and end positions, which were sampled randomly among the nodes of the triangular meshes. Using |$T$|, the geodesic to the initial position was obtained by solving the following differential equation (Eq. 2, Fig. 1, central panel).
where |$X(s)$| represents the geodesic and |$ds = 0.625$| mm (distance of half of the edge length). From the end position, we sequentially computed the curve on the surface following the gradient of the distance function |$T$| while completing the termination condition as follows. The distance from the starting position was less than |$ds = 0.625$| mm or the inner product of the movement vector was negative, i.e., the direction of movement was reversed. Thus, the geodesic was given by a sequence of points and a complementary curve between them.
We calculated 50,000 geodesics on the cerebral cortex and confirmed that the ensemble of the geodesics covered almost the entire surface area. As an indicator of the local density of the geodesics, we computed the passage frequency of the geodesics at each node of the mesh. If the position of points representing the geodesic was contained within a radius of |$1.25$| mm of each node, the geodesic was judged to pass through that node. The Gaussian curvature on the mesh was calculated using the Gaussian function in Houdini.
We visualized the correspondence of the geodesics to the association fibers by applying “OR” and “AND” operations. “OR” took the logical sum of the geodesics passing at the two separated regions, each of which included a population of nodes, while “AND” took the logical product.
We focused on bundles of geodesics connecting distant spots around a node with a local peak of geodesic density to make comparisons with the densely packed fibers. We also focused on bundles running between distant regions where low density nodes widely spread to make comparisons with the divergent and convergent properties of fibers.
Results
Geodesics were dense nearby endpoint of fissures and sulci and sparse at cortical poles
To investigate the spatial distribution of the geodesics on the cortex (at 27 GW), we first observed the passage density of the geodesics at each position (see Materials and methods). We found several regions exhibiting 100% higher density than the average over the cortical surface (Fig. 2A, red), while other regions exhibited 30% lower density than the average (Fig. 2A, blue). In the lateral view of the brain, high local density was observed in the perisylvian regions (red and yellow in Fig. 2A1), while low density was observed around the frontal, temporal, and occipital poles (blue in Fig. 2A1). Additionally, in the medial view, high local density emerged in the cingulate (red region around white-masked region in Fig. 2A2) and hippocampal gyri in the medial view (red in Fig. 2A2) and in the insula in the frontal view (red in Fig. 2A3). Hereafter, the narrow regions surrounding a high density node and the wide regions spreading with low density nodes are referred to as hot spots (HSs; Fig. 2A, red) and cold regions (CRs; Fig. 2A, blue), respectively. The HSs were surrounded by moderately high-density regions (yellow in Fig. 2A1). In summary, the distribution of geodesics on the fetal cortex depended on the anatomical regions.
Association of geodesic distribution with cortical surface geometry. Heat map for normalized density of the geodesics (A1-3) and Gaussian curvature (B1-3) on the brain surface at 27 GW, where the inset of each panel denotes the respective value. Reference orientations: anterior (A), posterior (P), dorsal (D), ventral (V), left (L), and right (R) are placed in the right upper corner of each panel. Scale bar = 10 mm.
We found five HSs (hereafter, referred to as HS1, 2, 3, 4, and 5), each with a specific spatial distribution (Fig. 2A1 and 2, HS1-5). In the lateral view, HS1 was located around the inferior part of the frontal cortex (Fig. 2A1, HS1). HS2 was located around the anteroventral endpoint of the central sulcus and expanded posteriorly on the parietal cortex and anteriorly on the frontal cortex (Fig. 2A1, HS2). HS3 was located around the endpoint of the Sylvian fissure (Fig. 2A1) and distributed anteriorly toward HS2 on the parietal cortex and ventrally on the temporal cortex (Fig. 2A1, HS3). In the lateral and frontal views, HS4 was located around the anterior part of the large intracortical insula (Fig. 2A1 and 3, HS4). In the medial view, HS5 was located around the parieto-occipital sulcus (Fig. 2A2, HS5). In addition, geodesics were densely distributed in a ring shape (hereafter referred to as Hot Ring, HR) in the cingulate and hippocampal gyri (Fig. 2A2, HR).
Moreover, we observed three CRs (hereafter referred to as CR1, 2, and 3) corresponding to regions surrounding the brain poles. CR1 was located around the frontal pole and expanded dorsally on the frontal cortex (Fig. 2A, CR1). CR2 was located around the temporal pole and expanded posteriorly on the temporal pole (Fig. 2A, CR2). CR3 was located at the occipital pole and surrounded the calcarine sulcus in the medial view (Fig. 2A, CR3).
To understand how the geodesic density relates to the local geometry, we measured the Gaussian curvature. The curvature on all five HSs and part of the HR tended to be negative, indicating a saddle shape at the endpoint of the gyri and sulci (Fig. 2B in blue or green). In contrast, all three CRs mostly had positive Gaussian curvature, exhibiting a dome-shape (Fig. 2B in red or yellow). These findings indicate that the typical geometrical landscape of the cortical surface biases the geodesic density.
Dense geodesics corresponded to the arrangement of two association fibers
The correspondence of the HSs and HR to the anatomical regions (Fig. 2A) prompted us to examine whether the geodesics further corresponded to major association tracts, such as the superior longitudinal fasciculus (slf)/arcuate fasciculus and cingulum (cg), in reference to the MRI atlas of the white matter of the human brain (Oishi et al. 2011; Catani et al. 2012).
First, the geodesics passing through the HSs on the lateral surface were delineated to compare their distribution with the slf or arcuate fasciculus, a lateral tract composed of long and short fibers connecting the perisylvian cortex to the frontal, parietal, and temporal lobes. When selecting a bundle of geodesics connecting HS1, HS2, and HS3, we chose all nodes within a |$2.5$| mm radius from the point with the highest geodesic density for each of the HSs. We separately selected the geodesics passing through both HS1 (left gray cubes, Fig. 3A) and HS2 (central gray cubes, Fig. 3A), as well as the geodesics passing through HS2 and HS3 (right gray cubes, Fig. 3A) using the “AND” operation. By further applying the “OR” operation to the geodesics passing through both HS1 and HS2 and those passing through both HS2 and HS3, we obtained a bundle of geodesics passing through HS1, 2, and 3 (see Materials and methods). We found that this geodesic bundle traveled along the anterior–posterior direction over the moderately high-density region (Fig. 2A, yellow) of the frontal and parietal regions (white lines, Fig. 3A). Moreover, many geodesics passing through HS3 curved towards the temporal lobe, forming a large arc on the perisylvian cortex as shown in Fig. 3A. Thus, the shape and trajectory of the obtained bundle correspond to those of the slf or arcuate fasciculus.
Dense geodesics. White line: geodesic. (A) The geodesic bundle obtained by applying the “OR” operation to two bundles selected by the “AND” operation passing through HS1 and HS2 and through HS2 and HS3. (B) The geodesic bundle obtained by applying the “OR” operation to three bundles selected by the “AND” operation passing through the anterior-dorsal part, dorsal-posterior part, and posterior-ventral part of the HR (anterior part: anterior cingulate gyrus, dorsal part: isthmus cingulate gyrus, posterior part: posterior cingulate gyrus, ventral part: parahippocampal gyrus). The colors of the heat map, reference orientations and scale are the same as those in Fig. 2A.
Second, the geodesics passing through the HR on the medial surface were delineated to compare their distribution with the cg, a medial tract running immediately dorsal to the corpus callosum and along the temporal lobe, forming a large, C-shaped trajectory (Oishi et al. 2011). For each of the 4 selected parts on the HR, nodes within a |$2.5$| mm radius were chosen. We selected four parts of the HR: the anterior part on the anterior cingulate gyrus (right gray cubes, Fig. 3B), the dorsal part on the isthmus cingulate gyrus (upper gray cubes, Fig. 3B), the posterior part on the posterior cingulate gyrus (left gray cubes, Fig. 3B), and the ventral part on the parahippocampal gyrus (lower gray cubes, Fig. 3B). We applied the “OR” operation to three bundles of geodesics: those passing the anterior & dorsal parts, those passing the dorsal & posterior parts, and those passing the posterior & ventral parts, which were selected by the “AND” operation. As a result, a circular bundle of geodesics passing through the anterior cingulate, isthmus cingulate, posterior cingulate, and parahippocampal gyrus was obtained as shown in Fig. 3B. Moreover, a part of the bundle extended toward HS5 along the parieto-occipital fissure and a hot region along the calcarine fissure (Figs. 3B and 2A2). Thus, the shape of the obtained bundle of geodesics corresponded to that of the cg. In summary, dense geodesics passing through HSs and HR correspond to the arrangement of slf and cg association fibers (Table 1).
Geodesic bundles passing through HSs and CRs corresponded to major association fibers.
| Geodesic bundles | Geodesic characteristics | Association fiber |
|---|---|---|
| HS1 - HS2 - HS3 | High density | Superior longitudinal fasciculus (slf) / Arcuate fasciculus |
| HR - HS5 | High density | Cingulum (cg) |
| CR1 - HS4 - CR2 | Inter-poles | Uncinate fasciculus (unc) |
| CR1 - HS4 - CR3 | Inter-poles | Inferior fronto-occipital fasciculus (ifo) |
| CR2 - HR - CR3 | Inter-poles | Inferior longitudinal fasciculus (ilf) |
| Geodesic bundles | Geodesic characteristics | Association fiber |
|---|---|---|
| HS1 - HS2 - HS3 | High density | Superior longitudinal fasciculus (slf) / Arcuate fasciculus |
| HR - HS5 | High density | Cingulum (cg) |
| CR1 - HS4 - CR2 | Inter-poles | Uncinate fasciculus (unc) |
| CR1 - HS4 - CR3 | Inter-poles | Inferior fronto-occipital fasciculus (ifo) |
| CR2 - HR - CR3 | Inter-poles | Inferior longitudinal fasciculus (ilf) |
Geodesic bundles passing through HSs and CRs corresponded to major association fibers.
| Geodesic bundles | Geodesic characteristics | Association fiber |
|---|---|---|
| HS1 - HS2 - HS3 | High density | Superior longitudinal fasciculus (slf) / Arcuate fasciculus |
| HR - HS5 | High density | Cingulum (cg) |
| CR1 - HS4 - CR2 | Inter-poles | Uncinate fasciculus (unc) |
| CR1 - HS4 - CR3 | Inter-poles | Inferior fronto-occipital fasciculus (ifo) |
| CR2 - HR - CR3 | Inter-poles | Inferior longitudinal fasciculus (ilf) |
| Geodesic bundles | Geodesic characteristics | Association fiber |
|---|---|---|
| HS1 - HS2 - HS3 | High density | Superior longitudinal fasciculus (slf) / Arcuate fasciculus |
| HR - HS5 | High density | Cingulum (cg) |
| CR1 - HS4 - CR2 | Inter-poles | Uncinate fasciculus (unc) |
| CR1 - HS4 - CR3 | Inter-poles | Inferior fronto-occipital fasciculus (ifo) |
| CR2 - HR - CR3 | Inter-poles | Inferior longitudinal fasciculus (ilf) |
Interconnected geodesics between broad regions around brain poles corresponded to the other association fibers
Since the widely spread regions with low geodesics density (CRs) were observed around each of the brain poles, we examined if bundles of geodesics connecting the lobes (hereafter referred to as inter-pole geodesics, CR1-CR2, CR1-CR3, and CR2-CR3, Fig. 4) show divergent and convergent pathways. To demarcate these regions, we chose nodes within a |$15$| mm radius from each of the poles: the frontal and temporal poles were defined as the most anterior position of the anterior-posterior axis coordinate and the occipital pole as the most posterior position.
Inter-pole geodesics connecting brain poles. The bundle of geodesics with start and end positions located within CR1 and CR2 (A); CR1 and CR3 (B); and CR2 and CR3 (C). The white arrow points to the position where the bundle passes through HS4. The colors of the heat map, reference orientations and scale are the same as those in Fig. 2A.
The inter-pole geodesics connecting the inferior part of the frontal cortex (CR1) and the temporal lobe (CR2) exhibited a converged passage through a narrow part of HS4 (Fig. 4A3). This bundle of geodesics corresponded to the uncinate fasciculus (unc), a tract that connects the medial and lateral orbitofrontal cortex with the anterior temporal lobe and passes through the narrow region of the external capsule corresponding to partial HS4 (white arrow, Fig. 4A3).
The inter-pole geodesics connecting CR1 and the inferior and medial occipital cortex (CR3) exhibited a converged passage through a narrow part of HS4 and the hot regions on the medial occipito-temporal cortex (Fig. 4B). This bundle corresponded to the inferior fronto-occipital fasciculus (ifo), a tract that connects the orbitofrontal cortex and the inferior and medial occipital lobes passing through the narrow region of the external capsule corresponding to partial HS4 (Fig. 4B2 and white arrow, 4B3). Furthermore, the tangential relative positions of the geodesic bundles (compare white arrow between Fig. 4B3 and 4A3) were consistent with the anatomical property in which the fibers of the ifo enter the external capsule dorsally to the fibers of the unc. Note that some of the bundles that densely passed through the medial parietal lobe (Fig. 4B2) may correspond to the superior fronto-occipital fasciculus, for which periventricular fasciculus frontooccipitalis was reported in fetuses (Vasung et al. 2011) but existence in adults is subject to debate (Schmahmann et al. 2009; Oishi et al. 2011; Catani et al. 2012).
The inter-pole geodesics between CR2 and CR3 passed through the ventro-medial HR on the occipito-temporal cortex, as shown in Fig. 4C in yellow. This bundle corresponded to the ilf, a tract connecting the occipital and temporal lobes. Moreover, the tangential relative position of this bundle with respect to the bundle between CR1 and CR3 (compare Fig. 4C1 with 4B1) was consistent with the anatomical property in which the ilf runs parallel and lateral to the inferior longitudinal fasciculi (ifo). In summary, the inter-pole geodesics for fronto-temporal, fronto-occipital, and temporo-occipital connectivity exhibit converged pathways corresponding to the position of inter-pole association fibers (Table 1).
Discussion
Summary of results
As a hypothetical landmark for tangential long-distance fiber arrangement, we calculated the spatial distribution of the geodesics on the fetal cortical surface, where association fiber prototypes have already emerged (Huang et al. 2009; Takahashi et al. 2012; Horgos et al. 2020). We found that geodesics were dense at the saddle-shaped anatomical regions (HSs and HR in Fig. 2) and sparse at the dome-shaped cortical poles (CRs in Fig. 2). The positions of the multiple geodesics passing through the HSs and HR corresponded to those of the two major association fibers (Fig. 3 and Table 1), while those connecting the CRs at the cortical poles via the HSs or HR corresponded to the 3 other major association fibers (Figs. 3 and 4, and Table 1). These correspondent spatial distributions between association fibers and geodesics support that geodesics guide the tangential position of the long-distance association fibers according to the curved cortical surface landscape. This geodesic theory expands upon the subplate nexus hypothesis to provide a better understanding of the direction information for wiring long-distance association fibers.
A hypothetical mechanism for the development of association fibers
Neuronal activity serves as a potential physiological mechanism for the geodesic-mediated directional guidance. Subplate neurons exhibit spontaneous activity forming neural circuits via gap junctions (Moore et al. 2014) and chemical synapses (Kristt and Molliver 1976; Luhmann and Khazipov 2018; Luhmann et al. 2022). These neuronal circuits have been shown to play an important role in tangential neuronal fibers (Kostović 2020), as well as in radial fibers (Bystron et al. 2008; Kanold and Luhmann 2010; Hadders-Algra 2018; Luhmann and Khazipov 2018; Kanold 2019; Molnár et al. 2019, 2020; Ohtaka-Maruyama 2020). In the tangential direction, neuronal activity can be propagated to the subplate neural circuits throughout the entire brain (Garaschuk et al. 2000; Sun and Luhmann 2007; Kilb et al. 2011; Keunen et al. 2017; Thomason et al. 2017). A theoretical study has shown that such propagation of neuronal excitatory activity follows a geodesics on the curved surface of propagated media (e.g. brain surface) (Brazhnik et al. 1988; Davydov et al. 2000; Horibe et al. 2019). This suggests that neuronal activity can be densely propagated at the HSs (Figs. 2 and 4), thereby guiding neuronal fibers in a similar way to Hebbian learning (Hebb 1949). The network of subplate neurons is connected to the cortical neurons, which migrate in the radial direction towards the surface (Rakic 1988; Ohtaka-Maruyama et al. 2018). Cortical neurons grow their tangential network, using the network of subplate neurons as the prototype. When the subplate disappears via apoptosis during the fetal period (Kostović et al. 2002; Vasung et al. 2016), the network of subplate neurons is replaced by corticocortical networks (i.e. association fibers networks) (Kostović 2020). Together, the previous and present studies suggest that the association fibers network is formed by tangential propagation of neuronal activities following geodesics on the smooth surface of the fetal brain.
The geodesic theory for association fibers is applicable to other developmental stages and species
While the geodesic theory primarily focuses on association fibers as its main subject, there are three different fiber systems (i.e. association, projection, and commissural fibers) for long-distance fibers. According to the classical statement, long, earlier-developing, fiber systems are closer to the ventricles, and later-developing short pathways are closer to the cortex (Žunić Išasegi et al. 2018). Before the onset of gyrification, the spatial patterns of major growing fibers are characterized by predominantly tangential organization and prominent crossroads of commissural, projection and association fibers (Judaš et al. 2005; Žunić Išasegi et al. 2018). During the midfetal period, the corpus callosum fibers are located closest to the ventricle; thalamic and basal forebrain fibers situated in the internal sagittal stratum and external sagittal stratum, respectively; and prospective association fibers develop along the deep border with the subplate (Žunić Išasegi et al. 2018). This tangential organization at different surfaces between different kinds of fibers would justify the geodesic theory that explains the formation of the association fibers on a particular two-dimensional surface. On the other hand, the fact that the association, projection and commissural fibers intersect in narrow crossroad areas (Judaš et al. 2005; Vasung et al. 2011, 2017; Milos et al. 2020) suggests that the three fiber systems interact while running in different direction in 3D space, which is beyond the scope of the geodesic theory. Elucidating the mechanism by which the entire fiber systems are formed in an orderly fashion is a task for future theoretical studies.
The geodesic theory for association fibers can be extended to other development stages. Smoothly curved landscapes, such as the saddle shape around the Sylvian fissure (Fig. 2A3, around HS4) and the dome shape at the brain poles (Fig. 2A1, around CR1-3) are maintained between 9 and 27 GW, despite the volumetric growth of the cortex (Huang et al. 2009; Horgos et al. 2020). For example, at 9 GW, primary folding has already begun to develop, resulting in a saddle-shaped Sylvian fissure, while frontal and temporal poles exhibit a dome shape (Horgos et al. 2020). This consistency in the surface landscape suggests that the HSs of geodesics on the saddles (Fig. 2A, HS4) and the CRs on the domes interconnected by the geodesics (Fig. 4A, CR1-3) exist even at 9 GW. Importantly, the unc appears to connect the brain poles around 9 GW, which is the earliest connection among the association fibers (Ouyang et al. 2019; Horgos et al. 2020). Thus, the position of association fiber formation could correspond to the geodesics earlier than 27 GW. At later stages, short association fibers (e.g. U-shaped fiber) emerge along with a shorter-distance landscape of the secondary sulci than that of primary sulci (Takahashi et al. 2012; Das and Takahashi 2018). Thus, whether the position of short-distance fibers corresponds to the geodesics of the shorter-distance landscape should also be evaluated. In summary, examining the geodesic theory in early and late fetal developmental stages would provide directional guidance of long and short association fibers, respectively.
The geodesic theory is further applicable to other species. Primates share the primary sulci; however, species evolution modifies their brain surface landscape such as the depth and length of the sulci (Amiez et al. 2019; Miller et al. 2020; Friedrich et al. 2021). This morphological evolution of shared anatomical characteristics enables the modification of the respective distribution of geodesics, potentially affecting that of the association fibers. For example, the arcuate fasciculus is well-developed (strongly connected) in humans, playing an important role in language function (Dick and Tremblay 2012); however, it is immaturely and weakly connected in other primates (Rilling et al. 2008). In the human fetal cortex, the spatial distribution of geodesics corresponded to those of the arcuate fasciculus (HS1, HS2, HS3; Figs. 2A and 3A). Future studies should examine whether and how the density and direction of geodesics correspond to those of the arcuate fasciculus and other association fibers in other primate cortices, which would provide further understanding of the evolution of association fibers constrained by the cortical landscape.
Limitation of the study
The present theory does not incorporate the radial thickness heterogeneity of the subplate depending on the growth stages and cortical regions (Huang et al. 2009; Garcia, Kroenke, et al. 2018a; Terashima et al. 2021). Nevertheless, the geodesics reproduced the relative tangential position of the fibers. For example, the inter-pole geodesics between fronto-occipital poles corresponding to the ifo were located more dorsally than those between fronto-temporal poles corresponding to the unc (Fig. 4A3 and B3 white arrows), consistent with the fact that the ifo passes through a more dorsal area than the unc around the external capsule (Catani et al. 2012). Furthermore, the geodesics corresponding to the ilf and ifo reproduce this tangential relative position, in which the ilf runs parallel to the outer part of the ifo (Catani et al. 2012) (Fig. 4B2 and C2). The effect of the thickness heterogeneity on the geodesic distribution can be evaluated if the sampling probability of the start and end positions of the geodesics is set to depend on the respective thickness; however, in this study, this probability was set to be homogeneous.
Additionally, the geodesic theory includes the cortical surface geometry but not the radial direction. Here we used data from the surface generated at the midthickness of the cortex (expected cortical layer 4) (Hill et al. 2010). This should be geometrically similar to the surface at the deepest border of the subplate, where association fibers start to grow (Žunić Išasegi et al. 2018). Although identifying the deepest border of surface in MRI data is challenging, methods to obtain a boundary between subplate and intermediate zone (Vasung et al. 2016) or between superficial and deep layers of the subplate (Pogledic et al. 2020) would be useful. Using these surface data may improve the results of the geodesic analysis in future.
Consistency of the geodesic theory with the tangential fiber distribution, presented here rather qualitatively, could be quantitatively evaluated by comparing the geodesics of cortical landscape data with diffusion tensor imaging (DTI) tract data. A research article containing both cortical landscape imaging (MRI) and DTI tract for the early developmental stages of the human brain was published (Vasung et al. 2017). The consistency between the DTI tracts and the cortical surface geodesics, both represented by 3D curved bundles, can be quantitatively measured by using the proximity of 2 different 3D curves (e.g., the Wasserstein distance (Villani 2009)), thereby enabling to statistically verify the geodesic theory.
Conclusion
Our geometrical analysis revealed that the distribution of geodesics on the cortical surface of human fetal brains follows that of the 5 major association fibers. This finding supports the geodesic theory, in which the tangential direction of the association fibers follows geodesics, suggesting a unified understanding on the development and evolution of association fibers constrained by the landscape of the cortical surface.
Acknowledgments
We thank K. Matsushita for discussion and J. Horikawa for technical support in Houdini.
Funding
This work was supported in part by funds from Japan Society for Promotion of Science Grants-in-Aid for Scientific Research (26220004, 19KK0247, and 23H05425) and Japan Agency for Medical Research and Development (JP21gm1310012) to GT and the Japan Science and Technology Agency (JPMJCR2121) to KF.
Conflict of interest statement: None declared.
Author contributions
Kazuya Horibe (Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing-original draft, Writing-review & editing), Gentaro Taga (Conceptualization, Funding acquisition, Investigation, Methodology, Validation, Writing-original draft, Writing-review & editing), Koichi Fujimoto (Conceptualization, Funding acquisition, Investigation, Methodology, Resources, Supervision, Validation, Writing-original draft, Writing-review & editing)
Data availability
The codes and mesh data related to this article are available at https://github.com/KazuyaHoribe/GeodesicTheory. The MRI data related to this article are available on BALSA (https://balsa.wustl.edu/).
References
