Serotonergic Axons as Fractional Brownian Motion Paths: Insights Into the Self-Organization of Regional Densities

Introduction

All cells in vertebrate brains are surrounded by a matrix of highly tortuous fibers that release serotonin (5-hydroxytryptamine, 5-HT), a major neurotransmitter. Presently the self-organization and dynamics of this matrix are not understood beyond neuroanatomical descriptions. Altered densities of serotonergic (serotonin-releasing) fibers have been associated with many mental disorders and conditions, including Autism Spectrum Disorder (Azmitia et al., 2011), Major Depressive Disorder (Numasawa et al., 2017), epilepsy (Maia et al., 2019), and exposure to 3,4-methylenedioxymethamphetamine (MDMA, Ecstasy) (Adori et al., 2011). Predictive computational models can significantly advance this research and support its biomedical applications. Motivated by this potential, we used a stochastic process framework to develop a model of serotonergic fibers and reproduced some key features of their density distribution in the mouse brain.

Serotonergic fibers are axons of neurons whose bodies are located in several brainstem clusters known as the raphe nuclei (Stuesse et al., 1991; Jacobs and Azmitia, 1992; Hornung, 2003). In mammals, these neurons mature early in development. They begin synthesizing serotonin around embryonic day 11–13 in the mouse and rat brains (Lidov and Molliver, 1982; Hendricks et al., 1999; Hawthorne et al., 2010) and around 5 weeks of gestation in the human brain (Sundstrom et al., 1993; Mai and Ashwell, 2004). In the adult mammalian brain, serotonergic axons are unusual in their ability to regenerate, with potential implications for the efforts to restore other axon systems after injury (Hawthorne et al., 2011; Jin et al., 2016; Kajstura et al., 2018). A recent study has shown that they may share this property with other axons in the ascending reticular activating system (Dougherty et al., 2019).

Recent studies have revealed great diversity of serotonergic neurons (Okaty et al., 2015, 2019), the functional significance of which is an active area of research (Ren et al., 2018). Paradoxically, transgenic mouse models with no serotonin synthesis in the brain during development have no gross neuroanatomical abnormalities and show only mild behavioral deficits (Mosienko et al., 2015; Pratelli and Pasqualetti, 2019). Serotonergic neurons can release other major neurotransmitters, such as glutamate (Okaty et al., 2019) and perhaps GABA (Stamp and Semba, 1995; Okaty et al., 2019). In the raphe nuclei, serotonergic neurons coexist with many other neurons (Cardozo Pinto et al., 2019; Schneeberger et al., 2019), some of which may participate in stereotypic synaptic arrangements (Soiza-Reilly et al., 2013).

The development of serotonergic fibers is currently conceptualized to proceed in two or three stages: the initial growth in well-defined fiber tracts, followed by extensive arborization and eventual dispersal in “terminal” fields (Lidov and Molliver, 1982; Carrera et al., 2008; Kiyasova and Gaspar, 2011; Jin et al., 2016; Maddaloni et al., 2017; Donovan et al., 2019). This orderly sequence is generally consistent with experimental observations at the level of local fiber populations, visualized with tract-tracing techniques or quantified with density measures. It is also theoretically appealing in that it mirrors the development of brain projections that connect two well-defined brain regions (e.g., the lateral geniculate nucleus and the primary visual cortex).

There is little doubt that the initial serotonergic projections form well-defined paths, which we have studied in our own research (Janušonis et al., 2004; Slaten et al., 2010). These paths are well-described because they are followed by axon bundles, thus facilitating their visualization in time and space. In contrast, the processes that lead to the formation of regional fiber densities remain poorly understood. Fundamentally, they require a rigorous description of the behavior of single fibers, each one of which has a unique, meandering trajectory. Since these processes are the main focus of the present study, we note several important challenges.

The first detailed morphological description of single serotonergic axons has become available only recently (Gagnon and Parent, 2014). This study has reconstructed a small set of axons originating in the dorsal raphe nucleus and has concluded that they travel through multiple brain regions, rarely branching in some of them and producing profuse arborizations in others. However, true branching points are difficult to distinguish from highly tortuous fiber segments that simply pass each other, especially in bright field microscopy (used in this study). A branching point can be unambiguously demonstrated only by examining an individual fiber at high resolution in all three-dimensions, within the physical section (Pratelli et al., 2017). Even when a confocal system with high-power objectives is used, a branching point can be difficult to distinguish from fibers crossing each other at sub-micrometer distances (Janušonis and Detering, 2019; Janušonis et al., 2019). The extent of local sprouting in regeneration also remains an open problem (Hawthorne et al., 2010; Jin et al., 2016), but current evidence suggests that it may not be significant (Jin et al., 2016). The rapidly developing methods of super-resolution microscopy and tissue expansion are well-positioned to provide definitive answers to some of these questions (Janušonis et al., 2019; Wassie et al., 2019).

If serotonergic fiber densities are determined by local arborization, it is unclear how fibers can be restricted to specific “terminal” regions, such as cerebral cortical layers (Linley et al., 2013). Because of their high degree of tortuosity, they are likely to cross over to adjacent areas, suggesting a subtle balance between region-specific branching and a diffusion-like process. Furthermore, the concept of “terminal region” is ambiguous for the serotonergic axons that typically do not form conventional synapses and can release serotonin at virtually any segment of their trajectory, based on in vivo and in vitro observations of axon varicosities (Benzekhroufa et al., 2009; Gagnon and Parent, 2014; Quentin et al., 2018). Serotonergic neurons also can release serotonin from the soma, dendrites, and growth cones, effectively making their entire membrane surface active (Ivgy-May et al., 1994; Quentin et al., 2018). It should be noted that serotonergic axons may also form conventional synapses (Papadopoulos et al., 1987), but the extent of this “wiring” transmission (Agnati and Fuxe, 2014) is currently unknown and continues to be debated.

In this study, we model individual serotonergic fibers as paths of a stochastic process that reflects their physical properties and show that regional arborization or other local control is not necessary to arrive at a good approximation of the observed fiber densities. Instead, these densities may strongly depend on the geometry of the brain. Our novel approach may offer insights into the self-organization of serotonergic densities in development and may also explain their stability in adulthood, without assuming the permanence of individual fiber trajectories. Consistent with this hypothesis, a recent study has shown that in the adult mouse brain regenerating serotonergic fibers do not follow the pathways left by degenerated fibers but can still restore the layer-specific densities after cortical injury (Jin et al., 2016).

We focus on Fractional Brownian Motion (FBM), a process first described under this name by Mandelbrot and Van Ness in 1968 (Mandelbrot and Van Ness, 1968). FBM and related stochastic processes have emerged as flexible and theoretically rich models in a variety of physical and biological systems. In particular, they have been used to understand the behavior of filamentous objects, such as biopolymer chains and chromosomes (Polovnikov et al., 2018, 2019).

The application of Brownian random walks to polymer chains dates back to the seminal works of Paul Flory, a founder of modern polymer physics and the winner of the 1974 Nobel Prize in Chemistry (Flory, 1986). In this theoretical approach, a fully flexible polymer in a poor solvent can be represented by a discrete random walk, in which each step corresponds to a monomer. However, these models are insufficient to reproduce polymer behavior in good solvents. Self-avoiding walks (SAWs) on lattices may be useful in these systems, but they are difficult to treat analytically and in computer simulations. FBM offers a more powerful theoretical solution, in which the fractal dimension of the path (the polymer chain in space) can be reproduced by a scaling exponent that controls the decay of the long-range correlation of the FBM steps (Chakravarty and Sebastian, 1997; Qian et al., 1998). This fractal dimension is a measure of how often the chain intersects itself (Hu and Nualart, 2005; Polovnikov et al., 2019). The same mathematical framework also allows modeling active motion, for instance, in living biological cells (Reverey et al., 2015). In the polymer language, these trajectories correspond to long chains that can show persistence in their spatial direction.

Theoretically, FBM extends the normal Brownian motion (BM), which for over a century has served as a standard model to describe simple diffusion and other similar processes (e.g., simple polymer dynamics and stock markets). While BM assumes independence between non-overlapping increments, FBM expands this model by allowing non-zero correlations. The sign and strength of the increment correlation is determined by the Hurst index (H), which defines two fundamentally different FBM regimes. If 0 < H < ½, two neighboring increments are negatively correlated, which produces highly jittery, “anti-persistent” trajectories (also known as “rough paths”). If ½ < H <1, two neighboring increments are positively correlated, which produces “persistent” trajectories that tend to maintain their current direction. In precise terms, the correlation between two neighboring increments is given by 2^2H−1 − 1 and the mean-square displacement follows the power law 〈x² 〉 ~ n^2H, where n is the number of performed steps. In this framework, BM becomes a special case of FBM, represented by H = ½. According to the scaling of the mean-square displacement, which can be (in the number or steps) sublinear (for 0 < H < ½) or superlinear (for ½ < H < 1), FBM can be classified as subdiffusion or superdiffusion, respectively. Since the trajectories of serotonergic fibers are considerably less jittery than BM and have the tendency to maintain their current direction, one can expect to capture their behavior with superdiffusive FBM. In the polymer language, H < ½ would increase the tendency to create coiled configurations, whereas H > ½ would lead to more stretched configurations. The ballistic limit H = 1 would correspond to a fully stretched chain.

In addition, FBM has four convenient properties that make it a natural choice in this context. First, it is a continuous process, which is consistent with the time-continuity of axon growth. Second, it has stationary increments, which informally means that its statistical properties do not change as the process evolves (this assumption is reasonable from the biological perspective). Third, it is a self-similar process, which ensures that the estimation of H does not depend on the discretization grid of experimental observations. Since time-dependent information is difficult to obtain in growing fibers [e.g., with time-lapse imaging in live animals Jin et al., 2016], this property ensures robustness. Fourth, its increments are normally distributed. Assuming that randomness in the fiber trajectory arises from collision-like events in its microenvironment and that each of these events has a small effect on the trajectory, the total effect of these collisions inevitably leads to a normal distribution (by the Central Limit Theorem). Importantly, FBM is the only stochastic process with all of these properties (assuming mean-zero increments).

Although FBM was introduced over a half-century ago, it poses major challenges in theoretical analyses. This is due to the fact that FBM is neither a Markovian process nor a semimartingale (in contrast to BM). A particularly important problem for biological sciences is the behavior of FBM in bounded domains (e.g., in two- or three-dimensional shapes). Reflected BM is well-understood (Ito and McKean, 1965), but it is not until very recently that the first description of the properties of reflected FBM (rFBM) have become available, in one-dimensional domains (Wada and Vojta, 2018; Guggenberger et al., 2019; Wada et al., 2019). The present study is the first application of this theoretical framework to serotonergic fiber distributions, on the whole-brain scale. We determine the steady state distributions of superdiffusive rFBM in constrained brain-like domains and show that they approximate neuroanatomical observations.

Materials and Methods

Immunohistochemistry and Imaging

Two adult male mice (C57BL/6J, 8 months of age, The Jackson Laboratory) were deeply anesthetized with a mixture of ketamine (200 mg/kg) and xylazine (20 mg/kg) and perfused transcardially with saline, followed by 4% paraformaldehyde. Their brains were dissected, postfixed in 4% paraformaldehyde overnight at 4°C, cryoprotected in 30% sucrose overnight at 4°C, and sectioned coronally at 40 μm thickness on a freezing microtome. The sections were rinsed in 0.1 M phosphate-buffered saline (PBS, pH 7.2), incubated in 0.3% H₂O₂ in PBS for 20 min to suppress endogenous peroxidase activity, rinsed in PBS (3 times, 5 min each), blocked in 5% normal donkey serum (NDS) and 0.3% Triton X-100 (TX), and incubated in rabbit anti-serotonin transporter (SERT) IgG (1:5000; ImmunoStar, #24330) with 5% NDS and 0.3% TX in PBS for 2.5 days at 4°C on a shaker. The sections were rinsed in PBS (3 times, 10 min each), incubated in biotinylated donkey anti-rabbit IgG (1:1000; Jackson ImmunoResearch, #711-065-152) with 2% NDS and 0.3% TX in PBS, rinsed in PBS (3 times, 10 min each), incubated in the avidin-biotin-peroxidase complex (1:100; Vector Laboratories, #PK-6100), rinsed in PBS (3 times, 10 min each), developed with 3,3'-diaminobenzidine and H₂O₂ using nickel intensification for 5 min (Vector Laboratories, #SK-4100), rinsed in PBS (3 times, 5 min each), mounted onto gelatin/chromium-subbed slides, allowed to air-dry, and coverslipped with Permount. They were imaged on the Zeiss Axio Imager Z1 system with the following objectives: 1 × (NA = 0.025; used for whole sections), 10 × (NA = 0.45), 20 × (NA = 0.80), and 40 × (NA = 1.30, oil). All procedures have been approved by the UCSB Institutional Animal Care and Use Committee.

Preparation of 2D-Shapes

The outer border, ventricular spaces, and major white matter tracts of imaged sections were hand-traced in Adobe Illustrator CC by an expert trained in neuroanatomy. To reduce the number of contours, major tracts at the edge of the section were left out of the border contour. The contours were imported into Wolfram Mathematica 12. The outer border contour was split along the median (sagittal) symmetry line. The right side of the contour was smoothed with a moving average and reflected on the left side. The same procedure was used for internal contours symmetric with respect to the median line (e.g., the cerebral aqueduct). Inner contours away from the median line (e.g., the fornix) were smoothed on the right side with a moving average and reflected on the left side. Therefore, the final digitized contours (rational-valued arrays of X- and Y-coordinates) were perfectly bilaterally symmetric, compensating for minor sectioning plane deviations and real (minor) brain asymmetries.

The obtained contours were next reformatted for FBM simulations. They were transformed into N × 2 matrices, the rows of which represented consecutive, integer-valued Y-coordinates and the two columns of which represented the leftmost and rightmost X-coordinates of the contour (also integer-valued). To arrive at this format, the original contour coordinates were divided by an integer factor that after rounding produced at least four X-values for each consecutive Y-value, and the minimal and maximal X-values were chosen. Since this procedure effectively reduced the size of the contour, it was enlarged back to its original size by multiplying the integer coordinates by the same factor and filling in the new, empty rows with X-values obtained by linear interpolation between the nearest available X-coordinates. Because this format cannot encode concavities oriented along the Y-axis (e.g., the third ventricle), such concavities were stored as separate inner contours. In the study, all such concavities were centered on the median line, which allowed their easy capture with the maximal X-value left to the line and the minimal X-value right to the line. For the purpose of this study, all inner contours were treated as impenetrable obstacles, irrespective of their physical nature (ventricular spaces, outer border concavities, white matter tracts).

The computer simulations described in detail below were performed on the Frontera supercomputing system (NSF, Texas Advanced Computing Center).

Discrete Reflected FBM

In the simulations, each individual serotonergic fiber was represented as the trajectory of a discrete two-dimensional FBM (Qian, 2003). Consider a random walker starting at position r₀ and moving according to the recursion relation r_n+1 = r_n + ξ_n, where ξ_n is a two-component fractional Gaussian noise. This means that the x and y components of ξ_n are Gaussian random numbers with zero average and variance σ², and that each component features long-range correlations between the steps (but the x and y motion are independent of each other). The corresponding covariance function is given by $\langle {\xi }_{x,m}{\xi }_{x,n+m}\rangle =\langle {\xi }_{y,m}{\xi }_{y,n+m}\rangle =\frac{1}{2}{\sigma }^2[|n+1{|}^{2H}-2|n{|}^{2H}+|n-1{|}^{2H}]$, where H is the Hurst index. The Fourier-filtering method (Makse et al., 1996) was employed to generate these long-range correlated random numbers on the computer (Wada and Vojta, 2018).

If the random walker encounters a boundary (i.e., an outer or inner contour), it is reflected. This reflection can be implemented in several different ways, including a repulsive potential force or a simple reflection condition that either prevents the trajectory from entering the forbidden region or mirrors it at the reflecting wall. Extensive test calculations (Vojta et al., 2020) have demonstrated that the choice of the reflection condition does not affect the resulting probability density for rFBM outside a narrow region (just a few steps wide) at the boundary. The following simulations thus employed the condition that a step that would lead the random walker into the forbidden region was simply not carried out.

Reflected FBM-Paths in Brain-Like Shapes

The goal of the FBM simulations in brain-like shapes was to model the two-dimensional distributions of serotonergic fiber densities at different rostro-caudal levels. The brain is a three-dimensional structure, but two-dimensional simulations are sufficient to approximate the densities, provided that the brain geometry does not change too rapidly in the direction perpendicular to the selected section. Serotonergic fibers do not form fascicles and show no preferred orientation (i.e., they are spatially isotropic) in terminal fields (Janušonis et al., 2019). As the three spatial coordinates in FBM are independent of each other, this implies that the motion perpendicular to the selected section can be neglected [for further information, see Vojta et al. (2020)].

Most simulations were performed with H = 0.8, but values between H = 0.3 and 0.9 were also studied for comparison. Each dataset consisted of 960 individual fibers (FBM trajectories) starting from random positions inside the shape. Each trajectory had 2²² ≈ 4 million steps of size (standard deviation) σ = 1.29 μm. Sample trajectories are shown in Figure 1. The local density (d_s) was determined by counting the total number of random walk segments inside each cell (of linear size 12.9 μm) of a square grid covering the shape. The trajectories were sufficiently long for the relative densities to reach a steady state (i.e., they did not change if the trajectory lengths were increased further).

For comparison between the simulated fiber densities and the densities observed in the actual (immunostained) sections, the simulated densities were scaled to “optical densities” by the (Beer-Lambert law-like) transformation d_o = 1 − exp(−kd_s), where the attenuation parameter k was chosen such that the mean pixel value in the simulated section matched the mean pixel value in the image of the immunostained section. This transformation removed the dependence on the arbitrarily chosen number and length of the FBM trajectories, and it also realistically constrained density values to a finite interval.

Reflected FMB-Paths in a Ring With a Varying Curvature

The boundary curvature may strongly affect the accumulation of fibers, which has direct implications for brain neuroanatomy. For example, in the coronal plane the curvature of the cerebral cortex can be high near the sagittal plane, as the outer border turns toward the corpus callosum, and low in the lateral convexities. The cortical gyri and sulci of the human brain create additional variability of local curvature, which can further modulate the accumulation of serotonergic fibers. Therefore, a qualitative relationship between boundary curvatures and fiber densities can allow theoretical predictions in neuroanatomy and may obviate the need for detailed computer simulations for each individual shape. However, caution should be exercised in shapes with complex geometry because superdiffusive rFBM is not a local process, due to its long-range correlations. Therefore, its properties can be potentially affected by the geometry of the entire domain. A systematic quantitative analysis of rFBM in geometrically idealized spatial domains (in two and three dimensions) falls outside the scope of the present analysis and is presented elsewhere (Vojta et al., 2020).

In order to gain an insight into how the density of FBM-trajectories varies as a function of H and the contour curvature, especially in neuroanatomically-relevant shapes, we investigated a ring-like shape bounded by the contours $R_{outer}=100(1+0.5{cos}^2(4\phi ))$ and R_inner = 50 (defined in polar coordinates, where φ ∈ [0, 2π)). This shape can represent a cross section through an abstracted vertebrate brain. Vertebrate brains develop from the neural tube and remain topologically tube-like in adulthood, where the “hollow” inside of the tube is the ventricular space filled with the cerebrospinal fluid (CSF).

These simulations were performed with H = 0.3, 0.5, and 0.8. Each data set consisted of 100 individual fibers (FBM trajectories), starting from random positions inside the shape. Each trajectory had 2²⁰ ≈ 1 million steps of size (standard deviation) σ = 0.4. In the subdiffusive case (H = 0.3), it was necessary to increase the step size to 1 to ensure that the simulations reached the stationary regime. As before, the local fiber density (d_s) was determined by counting the total number of random walk segments inside each cell (of linear size 1) of a square grid covering the shape.

Simulation of Reflected FBM-Paths in a 2D-Disk With Crowding

The brain tissue is a highly crowded space, but little quantitative information is available about the extent and geometry of the extracellular space in different brain regions (Hrabetova et al., 2018). In order to assess the sensitivity of our results to crowding effects, we performed FBM simulations in a large disk of radius R = 100, filled with 1,013 small disk-shaped obstacles of radius R_obs = 2. The obstacles were located within a circle of radius 90, leaving an outer ring of width 10 unoccupied. This shape was intended to represent a highly abstracted neocortical region, with the empty outer rim representing cortical layer I. In adulthood, this layer is nearly devoid of neuron somata (but contains many dendrites and axons that for simplicity were ignored in the simulation).

The simulations were performed with H = 0.8. Each dataset consisted of 192 individual fibers (FBM trajectories) starting from random positions inside the shape. Each trajectory had 2²⁵ ≈ 34 million steps of size (standard deviation) σ = 0.4. As before, the local fiber density (d_s) was determined by counting the total number of random walk segments inside each cell (of linear size 0.5) of a square grid covering the shape.

Results

Serotonergic fiber densities are typically described with regard to specific neuroanatomical brain regions, with the assumption that they reflect the region's functional demands and are supported by local biological factors. Our immunolabeling results (Figure 2) are consistent with the previously reported density maps in the rat and hamster brains (Steinbusch, 1981; Morin and Meyer-Bernstein, 1999). In particular, they show high fiber densities in the superficial layers of the cerebral cortex, consistent with observations in the rat and ferret (Voigt and de Lima, 1991a; Linley et al., 2013).

In contrast to previous interpretations, we suggest that these data indicate a general tendency of serotonergic fibers to accumulate near the borders of neural tissue (at the pial or ependymal surfaces) (Figures 2A–E). This observation has been previously made in some brain regions [e.g., in the hamster thalamus Morin and Meyer-Bernstein, 1999] but to our knowledge has never been extended to the entire brain. We next show that this tendency is consistent with the behavior of rFBM in the superdiffusion regime.

Supercomputing FBM simulations (with H = 0.8) were performed in four two-dimensional shapes that closely approximated the shape of actual coronal sections. Since major white matter tracts may act as forbidden regions (Figure 2F), we included them as “obstacles” in the simulations. Qualitatively, FBM sample paths closely resembled the trajectories of individual serotonergic fibers (Figure 1).

The comparisons between the simulated and actual fiber densities are shown in Figures 3–6. These densities showed similar increases at the outer brain border, irrespective of the rostro-caudal level. A similar trend was observed around the ventricles. Interestingly, high simulated fiber densities were obtained in some neuroanatomically-defined brain regions, such as the lateral geniculate nucleus and the hypothalamus (Figures 4, 5), even though these regions were not specifically modeled and the increase was induced purely by the contour geometry. It should be noted that both of these regions have a convex border with a relatively high curvature. We investigated this potential association using a simple shape and a range of H values (Figure 7), which further supported this conjecture. It leads to verifiable predictions of how serotonergic fiber densities can vary across brain regions and may also support comparative neuroanatomy, where differences in the fiber densities across mammalian species may be caused, at least in part, by differences in the brain shapes (despite the highly similar neuroanatomical plans).

At the outer border, a considerable mismatch was observed between some gradients of densities (e.g., Figure 4E). Considering the neuroanatomical simplicity of the simulated shape (e.g., it contained no “cells”), this result is not surprising. Also, the simulated gradients depend on the value of H and the attenuation parameter of “optical transformation” (e.g., they can be made less steep, with an effect on the overall density intensity). Matching the simulated and actual gradients precisely is difficult because the true fiber density cannot be determined in immunostained sections without tracing every fiber. Also, many non-linear effects can take place between the section and the image sensor, even at optimal illumination settings.

The strongest discrepancy was found between the strong density spikes around white matter tracts in the simulations (where the tracts were modeled as impenetrable “obstacles”) and the virtual absence of such spikes in the actual (immunostained) sections. This suggests that these tracts cannot be modeled as “hard” obstacles and that other FBM-reflection models may reflect their properties more accurately.

Since a fixed H was used in the simulations, we investigated the sensitivity of the obtained results to a range of H values (Figures 6, 7). The density distribution patterns varied dramatically across the three diffusion regimes, as anticipated (Wada and Vojta, 2018). However, they were robust within the superdiffusion regime, suggesting that the results can be safely generalized to other H values, beyond the one that was used in the simulations (0.8).

Finally, we examined the sensitivity of the results to cell packing. Heterogeneous crowding is an essential property of neural tissue (Hrabetova et al., 2018), but currently little is known about the variability of the crowding density in brain regions. In the present study, we used a relatively simple model and tested whether the accumulation of fibers at borders could be reversed by many cell-like obstacles that can potentially retard or trap fibers in the interior region of the bounded domain (Figure 8). Despite the presence of many obstacle surfaces, the simulation produced only thin, high-density layers around the obstacles, with no major effect on the overall density distribution. This result supports the robustness of our findings and is generally consistent with experimental observations in the cerebral cortex, where the high density of serotonergic fibers may not be restricted to cortical layer I (which is virtually devoid of neuron somata) and may also include layers II-III (which contain densely packed neurons) (Voigt and de Lima, 1991b). However, this finding should not be overgeneralized without future studies of the interactions among different shape geometries, heterogeneous crowding, and the physical constraints of the serotonergic fibers (e.g., the simulated trajectories in our study had no physical width).

Discussion

We introduce a novel approach to the self-organization of serotonergic fibers in the brain and demonstrate that rFBM can replicate the behavior of these fibers at “hard” borders (such as the pial and ependymal surfaces of the brain). In contrast, major white matter tracts do not appear to significantly “reflect” serotonergic fibers, even though they constrain their trajectories as obstacles. This phenomenon may be due to the fact that many individual fibers can penetrate (and perhaps traverse) these tracts (Figure 2F). The presence of serotonergic fibers in some major pathways, such as the fasciculus retroflexus, has been noted in early neuroanatomical studies (Lidov and Molliver, 1982). These phenomena can be included in our model by treating the boundaries of white matter tracts differently from other boundaries, such as modeling them as soft repulsive potentials (Vojta et al., 2020).

Our model does not include the rich brain architecture and biological signals that may act on serotonergic fibers. While the simulated fiber densities are a good approximation of the actual fiber densities, some of these factors may be important for more accurate predictions. We briefly mention some of them.

We used two-dimensional shapes, even though the brain is a three-dimensional object. As explained in the section Reflected FBM-paths in Brain-like Shapes, the reduction to two dimensions unlikely produced major distortions. However, some discrepancies between the actual and simulated densities may be due to the geometry of coronal sections just rostral or caudal to the selected section, with some fibers “spilling into” or “leaking from” the examined section. This problem may be particularly significant where the geometry of ventricular spaces rapidly changes in the rostro-caudal direction (e.g., ventricular spaces “fuse” or “separate” in two-dimensional projections). Conceptually, the computational approach can be easily extended to three dimensions, provided a fully three-dimensional model of the brain geometry is available. However, the numerical effort for such simulations would be significantly higher. Specifically, it would require the parametrization of a two-dimensional manifold (the boundary), leading to increased computational costs close to the boundary (where the decision has to be made whether the walker is in the allowed region or not). One potential solution is using two-dimensional brain shapes obtained in two or three perpendicular planes (e.g., coronal, sagittal, axial), which may allow a more accurate reconstruction of densities in all three dimensions.

The brain can be viewed as a highly heterogeneous material. It is densely packed with cells, their processes, microvasculature, and other elements that serotonergic fibers cannot penetrate. Axons can be assembled into major white matter tracts which may require detours; for example, the number of axons in the primate corpus callosum can be on the order of 10⁷-10⁸ (Doty, 2007). The fine structure of the extracellular space (ECS) in the brain requires state-of-the-art experimental methods and currently is an active area of research (Nicholson and Hrabetova, 2017; Hrabetova et al., 2018). The stochastic geometry of the ECS may impose a particular covariance structure on traveling fibers (or a “memory effect”) (Morgado et al., 2002; Bénichou et al., 2013), possibly in a region-specific manner. Some of this geometry can be modeled with sphere packing (to represent cell bodies) (Picka, 2012). Similar problems arise in intracellular environments (Smith et al., 2017). In addition to the “hard” geometry, the stochastic properties of fibers may be affected by the viscoelastic properties of their environment (Cherstvy et al., 2019). In the case of FBM, it may lead to different H values in different brain regions, with implications for local fiber densities. Conversely, the estimated H values of fibers can be potentially used to obtain information about the structure and integrity of the ECS, with possible biomedical applications.

Since individual serotonergic fibers can be visualized, estimates of H can be obtained from experimental data. However, it requires overcoming a few technical challenges, which are a focus of our larger research program (in the present study, we assumed H = 0.8 and demonstrated the robustness of the main results). For example, tracing a single fiber in immunostained brain sections (e.g., 40 μm in thickness) is difficult because most fibers exit the section in the Z-direction, before advancing substantially in the X-Y plane (Janušonis et al., 2019). Modern tissue clearing methods and light-sheet microscopy allow direct 3D-imaging with no sectioning (Mano et al., 2018; Hillman et al., 2019), but light-sheet microscopy is only now approaching the sub-micrometer resolution needed for the imaging of individual serotonergic fibers (Chakraborty et al., 2019). The high fiber densities in most brain regions presents another problem, where single-fiber tracing has to be performed in the presence of interfering signals, such as other fibers and potential branching points. Advances in computer image analysis (Kayasandik et al., 2018; Falk et al., 2019), combined with transgenic technologies that allow labeling individual neurons and their processes with unique combinations of fluorophores [such as Brainbow 3.2 Cai et al., 2013], are well positioned to advance these efforts.

Our model assumes no self-avoidance or biological feedback signals that depend on fiber density. Several factors have been reported to control the growth of distribution of serotonergic fibers, but this research has been heavily influenced by the notion that orderly distribution cannot be achieved without tight biological control [see, e.g., Chen et al. (2017)]. The simulation results provide evidence that a considerable degree of self-organization can be achieved with simple assumptions, but it does not rule out these factors. Early studies have found that S100β, a biologically active protein, may promote the development of serotonergic fibers. Interestingly, S100β can be released from astrocytes, in response of activation of serotonin 5-HT_1A receptors, suggesting a positive feedback loop (Whitaker-Azmitia, 2001). Also, the absence of brain serotonin synthesis alters the development of normal fiber densities, but this effect appears to be strongly region-dependent (Migliarini et al., 2013). It may be mediated by the brain-derived neurotrophic factor (BDNF) which has long been implicated in the growth of serotonergic fibers (Mamounas et al., 1995; Migliarini et al., 2013). It remains unclear how serotonin affects fiber densities under physiological conditions because a different genetic model (with a less severe reduction of brain serotonin levels) has failed to reproduce these effects (Donovan et al., 2019). Recently, protocadherin-αc2 has been strongly implicated in the distribution of serotonergic fibers, through homophilic interaction between individual axons (Katori et al., 2009, 2017; Chen et al., 2017). Intriguingly, protocadherin-α mutants show pronounced increases in the fiber densities in layer I of the primary motor cortex and in the lacunosum-moleculare layer of the hippocampus, both of which are at the border of their respective brain regions. This study also has noted that in some other regions the “distribution of serotonin axonal terminals was […] dense at the periphery of each region but sparse in the center” (Katori et al., 2009). Our results suggest that this experimental result may reflect either a more pronounced rFBM behavior (in the absence of fiber interaction) or an rFBM with a higher H (induced by the genetic mutation).

It is important to note that FBM is not the only mathematical model that allows superdiffusion (or anomalous diffusion, more generally). Superdiffusion can also be modeled with Lévy flights or with continuous-time random walks (CTRWs) that have a heavy-tailed displacement probability density (Codling et al., 2008; Schulz et al., 2013; Metzler et al., 2014). These two processes allow large, instantaneous (spatially discontinuous) jumps from one location to another, but CTRWs are additionally parametrized by the waiting-time probability density (therefore, a Lévy flight can be viewed as a special CTRW). They are highly appropriate for some physical processes [e.g., the dynamics of molecular complexes jumping from one segment of a polymer to another, facilitated by folding-induced physical proximity Lomholt et al., 2005], but they do not well-represent axon growth where instantaneous jumps are not biologically realistic. Also, the trajectories of growing axons are likely to show long-range temporal correlations that are an inherent property of FBM (in addition to other useful properties reviewed in the Introduction; here we assume H ≠ 1/2). It should be noted that the long-range correlations in FBM extend to arbitrarily large distances (Biagini et al., 2010), which may exceed biological reality, but a possible theoretical refinement may be provided by stochastic processes in which the long-range correlations are cut off at a large but finite distance (Molina-Garcia et al., 2018). In addition, branching FBM-like processes may offer insights into how the bifurcation or arborization of serotonergic fibers can affect their steady state distribution. Direct simulations of macromolecular dynamics in confined domains can further enrich these studies; for example, in some simulations particles near the wall tend to stay near the wall (Chow and Skolnick, 2015), which may explain the tendency of serotonergic fibers to orient parallel to the edge immediately below the pia (for depths up to 25–50 μm; Figure 2). Finally, it has been recently demonstrated (Vojta et al., 2019) that the increased density close to a boundary arises from the non-equilibrium nature of FBM. A similar anomalous diffusion process in thermal equilibrium, modeled by the fractional Langevin equation, does not lead to accumulation at the boundary. This property of FBM is consistent with the active growth of the serotonergic fibers.

In summary, the present study demonstrates that FBM offers a promising theoretical framework for the modeling of serotonergic fibers. Since serotonin-releasing fibers are a part of the larger ascending reticular activating system, which releases other major neurotransmitters and has widespread projections, this framework may also be useful in advancing the understanding of other stochastic axon systems in vertebrate brains.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Ethics Statement

The animal study was reviewed and approved by UCSB Institutional Animal Care and Use Committee.

Author Contributions

SJ proposed the hypothesis of serotonergic fibers as paths of a stochastic process, performed the immunostaining, prepared the 2D-shape matrices for simulations, and wrote the first draft of the manuscript. TV guided the computational analyses of reflected FBM and performed all supercomputing simulations. ND suggested FBM as a potential model that allows scale-invariance and made other theoretical contributions. RM led the development of the model within the theoretical framework of anomalous diffusion processes. All authors are the Principal Investigators of their respective research programs.

Funding

This research was supported by the National Science Foundation (grants #1822517 and #1921515 to SJ and ND), the National Institute of Mental Health (grant #MH117488 to SJ and ND), the California NanoSystems Institute (Challenge grants to SJ and ND), the Research Corporation for Science Advancement (a Cottrell SEED Award to TV), and the German Research Foundation (DFG grant #ME 1535/7-1 to RM). RM acknowledges support from the Foundation of Polish Science through an Alexander von Humboldt Polish Honorary Research Scholarship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Kasie Mays (UCSB) for the assistance with tissue preparation.

References