An experimental platform for stochastic analyses of single serotonergic fibers in the mouse brain

The self-organization of the serotonergic matrix, a massive axon meshwork in all vertebrate brains, is driven by the structural and dynamical properties of its constitutive elements. Each of these elements, a single serotonergic axon (fiber), has a unique trajectory and can be supported by a soma that executes one of the many available transcriptional programs. This “individuality” of serotonergic neurons necessitates the development of specialized methods for single-fiber analyses, both at the experimental and theoretical levels. We developed an integrated platform that facilitates experimental isolation of single serotonergic fibers in brain tissue, including regions with high fiber densities, and demonstrated the potential of their quantitative analyses based on stochastic modeling. Single fibers were visualized using two transgenic mouse models, one of which is the first implementation of the Brainbow toolbox in this system. The trajectories of serotonergic fibers were automatically traced in the three spatial dimensions with a novel algorithm, and their properties were captured with a single parameter associated with the directional von Mises-Fisher probability distribution. The system represents an end-to-end workflow that can be imported into various studies, including those investigating serotonergic dysfunction in brain disorders. It also supports new research directions inspired by single-fiber analyses in the serotonergic matrix, including supercomputing simulations and modeling in physics.

1. Introduction

The structural and dynamical properties of single serotonergic fibers determine their population-scale behavior and ultimately lead to the emergence of the serotonergic matrix in all vertebrate brains. This matrix, composed of a very large number of meandering fibers and their branches, is invariably present in mammals (Wilson and Molliver, 1991; Awasthi et al., 2021), birds (Sako et al., 1986; Challet et al., 1996), reptiles (Bennis et al., 1990; Challet et al., 1991), amphibians (Ueda et al., 1984; Bhat and Ganesh, 2023), bony fishes (Lillesaar, 2011; López and González, 2014), and cartilaginous fishes (Stuesse et al., 1995; Carrera et al., 2008). The fundamental role of the serotonergic matrix remains enigmatic; in adulthood, it may support neuroplasticity (Lesch and Waider, 2012; Teissier et al., 2017; Huang et al., 2021; Daws et al., 2022; Morgan et al., 2023), with potential implications for artificial neural networks (Lee et al., 2022).

The serotonergic fibers are long axons of neurons that express the tryptophan hydroxylase 2 gene (Tph2) and can synthesize the neurotransmitter serotonin (5-hydroxytryptamine, 5-HT) from the amino acid tryptophan. In mammals, their somata are restricted to the raphe nuclei of the brainstem. The rostral raphe nuclei (the dorsal raphe (DR), the median raphe, and the caudal linear nucleus) project anteriorly to supply fibers to virtually all midbrain and forebrain brain regions (Jacobs and Azmitia, 1992; Hornung, 2003).

Understanding the deeper structure of this system may require novel conceptual frameworks. Serotonin is not the only neuroactive compound that can be produced from tryptophan [e.g., tryptophan can also be metabolized into kynurenic acid, quinolinic acid, or N, N-dimethyltryptamine (DMT)] (Dean, 2018; Dean et al., 2019; Yabut et al., 2019; Paul et al., 2022; Vargas et al., 2023), many serotonergic neurons also release glutamate and other neurotransmitters (and therefore can be considered “other-ergic” neurons that co-release serotonin) (Okaty et al., 2019), and the population of the serotonergic neurons is now known to consist of transcriptionally diverse, nested cell clusters (Okaty et al., 2019; Ren et al., 2019; Okaty et al., 2020), some with distinct neuroanatomical profiles (Ren et al., 2018; Okaty et al., 2020).

The transcriptomics studies have highlighted the importance of the understanding of serotonergic neurons at the single-cell level (Okaty et al., 2020). Individual neurons can choose among several “adjacent” transcriptional networks, depending on their developmental path and activation state (Okaty et al., 2019), and are also affected by stochastic noise (Elowitz et al., 2002; Raj and van Oudenaarden, 2008). These networks are likely to represent dynamical attractors and are therefore computationally constrained (i.e., not every set of gene activity is equally stable).

Further, individual serotonergic fibers acquire unique identities because their trajectories are strongly stochastic (Janušonis and Detering, 2019). It does not mean that any trajectory is possible in practice. Individual trajectories may represent “sample paths” (realizations) of rigorously-definable spatial stochastic processes, which again implies computational constrains. We have demonstrated the potential of this conceptualization in recent studies in which serotonergic axons were modeled as paths of reflected fractional Brownian motion (FBM), a stochastic continuous-time process (Janušonis et al., 2020, 2023). In these supercomputing simulations, performed in geometric shapes based on the mouse brain, the resultant fiber densities approximated the actual serotonergic fiber densities, with no other biological information. These models require improvements to reflect neural tissue heterogeneities (Wang et al., 2023), but they can eventually bridge the stochasticity on the microscopic scale (the uniqueness of fiber trajectories in each brain) and the determinism on the macroscopic scale [the predictability of regional fiber densities (Awasthi et al., 2021)].

The serotonergic neurons are among the first neurons to acquire a mature phenotype in the developing brain (Lidov and Molliver, 1982; Hendricks et al., 1999; Hawthorne et al., 2010). Somewhat paradoxically, they also appear to remain highly dynamic in the adult brain. In particular, serotonergic fibers can robustly regenerate after a traumatic injury (Jin et al., 2016; Kajstura et al., 2018; Cooke et al., 2022), with new paths (Jin et al., 2016). This potential and the challenge of maintaining the physical continuity of extremely long, unfasciculated fibers suggests that a healthy adult brain may always contain degenerating and regenerating serotonergic fibers. This in turn implies that the fiber trajectories of given serotonergic neurons may change during an individual’s life span. Because of technical challenges, this prediction currently cannot be directly verified in experimental studies. However, it is known that the maintenance of the serotonergic fiber matrix depends on the adult expression of Lmx1b and Pet1, two genes essential for the embryonic development of serotonergic neurons (Kitt et al., 2022), and that the stochastic serotonergic trajectories acquire interesting computational properties if they do not stay constant (Lee et al., 2022).

The large-scale properties of the serotonergic matrix can be theoretically predicted from the properties of single serotonergic fibers—but not the other way around. Specifically, the same patterns of regional fiber densities can be produced by different models, and major alterations in fiber densities, such as in autism (Azmitia et al., 2011) or epilepsy (Maia et al., 2019), does not imply specific changes in the properties of single fibers. Importantly, these properties include not only trajectories, but the pattern and probability distribution of branching events, the spatial and temporal features of varicosities (fiber “swellings”), and other features. The branching of serotonergic fibers remains poorly understood because their paths often cross at distances that fall below the limit of optical resolution, making the definitive identification of branching points difficult even in three-dimensional, high-resolution images (Janušonis et al., 2019). Branching serotonergic fibers are more easily accessible in primary brainstem cultures (Hingorani et al., 2022), but the artificial culture environments limit the applicability of these observations in natural neural tissue. The field also lacks a fundamental understanding of the serotonergic varicosities. Their distribution may reflect the developmental age of the fiber (Maddaloni et al., 2017), but the same serotonergic fiber may also contain varicosities that differ in size, shape, and spacing—both in brain tissue (Baizer, 2001; Gagnon and Parent, 2014; Maddaloni et al., 2017) and in vitro (Hingorani et al., 2022). This variability does not necessarily imply different varicosity types and may instead reflect a dynamic process frozen at different stages along the fiber trajectory (in fixed preparations). Serotonergic varicosities are likely to be fluid structures, as suggested by experimental studies of mice and Drosophila with elevated serotonin levels (Daubert et al., 2010), a rat model of epilepsy (Maia et al., 2019), and analyses of varicosities in other neurotransmitter systems. Specifically, the appearance of varicosities may reflect the current state of the fiber segment, the effects of its local microenvironment, and purely stochastic fluctuations (Hellwig et al., 1994; Hatada et al., 1999; Shepherd et al., 2002; Gu, 2021; Ma et al., 2022). Notably, the dynamics of varicosities on the same fiber are likely to be correlated, necessitating unambiguous discrimination among individual fibers that often cross and intertwine in the same spatial location. An added complexity is that different serotonergic fibers in the same location may execute different transcriptional programs. Therefore, disregarding the identity of fibers is equivalent to a gross distortion of the covariance structure among these various elements, potentially leading to incorrect statistical calculations and conclusions.

At present, reliable identification and analyses of extended segments of single serotonergic fibers present serious challenges, especially in brain tissue. In particular, it requires knowledge about the continuity of a given fiber. Only a handful of recent studies have focused on descriptions of single serotonergic fibers at this level of precision (Gagnon and Parent, 2014; Maddaloni et al., 2017; Janušonis et al., 2019; Hingorani et al., 2022). Furthermore, specialized quantitative methods for their analyses are virtually absent. For example, the tortuosity index, often used to describe meandering fibers (Jin et al., 2016; Pratelli et al., 2017), is a relatively unstable metric in that it depends on the measured fiber length and imaging dimensionality and also can vary drastically in different realizations of the same spatial stochastic process.

This study aimed to (i) develop flexible and accessible experimental mouse models to isolate single serotonergic fibers in brain tissue, independently of their serotonin content or the expression of the serotonin transporter (SERT), (ii) build an algorithm for high-precision tracing of single fibers in noisy three-dimensional images, and (iii) develop a statistical model to achieve a succinct description of fiber trajectories that can vary strongly in their individual appearance. The three parts constitute a workflow that can be imported into various experimental setups.

2. Materials and methods

2.1. Transgenic mice

Three transgenic mouse lines were used in the study: Tph2-iCreER (Tg (Tph2-icre/ERT2)6Gloss/J; RRID: IMSR_JAX:016584), ROSA^mT/mG (B6.129(Cg)-Gt(ROSA)26Sor^{tm4(ACTB-tdTomato,-EGFP)Luo}/J; RRID: IMSR_JAX:007676), and Brainbow 3.2 (line 7) (Tg(Thy1-Brainbow3.2)7Jrs/J; RRID: IMSR_JAX:021227). The Tph2-iCreER mice express a tamoxifen-inducible Cre recombinase under the control of the promoter of the Tph2 gene that encodes a key enzyme in the serotonin synthesis pathway. The ROSA^mT/mG (Cre reporter) mice ubiquitously express a cell membrane-targeting tdTomato, which after Cre-recombination is replaced with a cell membrane-targeting enhanced green fluorescent protein (EGFP) (Muzumdar et al., 2007). The Brainbow mice contain the Brainbow 3.2 construct under neuron-specific regulatory elements from the Thy1 gene (Cai et al., 2013). This construct includes sequences for three cell membrane-targeting (farnesylated) fluorescent proteins (mOrange2, EGFP, and mKate2) that after Cre-recombination are expressed in random intensity combinations in individual cells. The construct also includes a fourth sequence (in the stop cassette) for a mutated (non-fluorescent and nucleus-targeting), Phialidium-derived yellow fluorescent protein [PhiYFP(Y65A)] that can be used to assess the transcriptional activity of the locus in specific brain regions, before recombination (Cai et al., 2013).

The breeders were purchased from The Jackson Laboratory (JAX), and the colonies were kept in a vivarium on a 12:12 light–dark cycle with free access to food and water. The ROSA^mT/mG mice were kept and bred in the homozygous state, as suggested by the supplier. Mouse litters, including the offspring of crosses, were PCR-genotyped by using DNA from toe biopsies collected at postnatal day 7. Genomic DNA was isolated with the QIAamp Fast DNA Tissue Kit (Qiagen #51404), and the transgene sequences were amplified with an Eppendorf Mastercycler pro S using iTaq DNA polymerase (Bio-Rad #1708870), a dNTP mix (Bio-Rad #1708874), and primers with sequences provided by JAX (using the cycling protocol suggested by the supplier of the polymerase). The amplicons were analyzed with a 2200 TapeStation (Agilent Technologies). All animal procedures have been approved by the UCSB Institutional Animal Care and Use Committee.

2.2. EGFP labeling of serotonergic fibers

Tph2-iCreER and ROSA^mT/mG mice were crossed, and their offspring with the transgenes were allowed to reach adulthood (at least 10 weeks of age). The induction of Cre-recombination was achieved by administering tamoxifen (Millipore-Sigma #T5648) dissolved in corn oil (Millipore-Sigma #C8267). The stock concentration (20 mg/mL) was stored at −20°C, and mice were injected intraperitoneally with around 0.1 mL of the solution (at 75 mg/kg) for 5 consecutive days. Following the tamoxifen treatment, the mice were allowed to survive for 1 week to 1 month. The brains of eight mice were processed with immunohistochemistry and examined with microscopy.

2.3. Stochastic multicolor labeling of serotonergic fibers with Brainbow AAVs

The Brainbow adeno-associated viruses (AAVs) (AAV-EF1a-BbTagBY [Addgene #45185-AAV9] and AAV-EF1a-BbChT [Addgene #45186-AAV9]) (Cai et al., 2013) were stored at −75°C. Before use, they were diluted to 1.5 × 10¹² vg/mL each in a tube containing sterile, alcohol-free saline. Adult (at least 10 weeks of age) Tph2-iCreER mice were anesthetized with an intraperitoneal injection of a mixture of ketamine (100 mg/kg) and xylazine (10 mg/kg) and placed in a small-animal stereotaxic frame. Further anesthesia was maintained with inhaled isoflurane, and the animal’s core temperature was maintained at 37°C using a TCAT-2DF temperature controller (a closed loop rectal probe-heating pad system; Physitemp Instruments). The incision area (with hair removed) was given a subcutaneous injection of lidocaine and disinfected with Betadine. A rostrocaudal skin incision was made with a scalpel, and a small hole was drilled in the skull directly over the DR. The needle of a 10 μL-Hamilton microsyringe was lowered into the dorsal raphe, and 1–2 μL of the AAV mixture was slowly pressure-injected at the lambda, 3.5 mm ventral to the dura. The needle remained in the DR for 5–10 min and was slowly withdrawn. The skin incision was closed with sterile wound clips, and the animal was monitored until it fully recovered.

One week after the surgery, the induction of Cre-recombination was achieved as described for the EGFP labeling. Following the tamoxifen treatment, the mice were allowed to survive for 1–3 months to ensure fluorophore transport to distal axon segments.

In this study, nine mice received Brainbow-AAV injections (after a pilot series to fine-adjust the stereotaxic coordinates). Of these cases, two mice had no detectable fibers in the DR, three mice had sparsely labeled fibers in the DR and in some other brain regions, and four mice had densely labeled fibers in the DR, with many fibers in other brain regions (all after signal enhancement with immunohistochemistry). In one mouse in the latter group, the densities of Brainbow-labeled fibers in the entire brain were comparable to those of the recently published monochrome map (Awasthi et al., 2021). The actual success rate of the procedure is difficult to assess because of gradual technical improvements and unexpected delays in tissue processing caused by the COVID-19 pandemic (in 2020–2021). Based on our experience, we recommend long transport times (around 3 months) and storing sections in a cryoprotectant at −20° (if they do not get immunostained in the next several days).

We also attempted to label serotonergic fibers by crossing Tph2-iCreER and Brainbow mice. Their offspring with the transgenes were allowed to reach adulthood (at least 10 weeks of age), and the induction of Cre-recombination was performed as described for the EGFP labeling. Following the tamoxifen treatment, the mice were allowed to survive for around 1 month.

2.4. Immunohistochemistry

Mice were euthanized with CO₂ and their brains were dissected into 0.1 M phosphate-buffered saline (PBS, pH 7.2). They were immediately immersion-fixed in 4% paraformaldehyde overnight at 4°C, cryoprotected in 30% sucrose for 2 days at 4°C, and sectioned coronally at 40 μm thickness on a freezing microtome. The sections were stored in PBS and processed further with immunohistochemistry (IHC).

2.4.1. Cre IHC

Sections were rinsed in PBS and treated with a heat-induced epitope retrieval (HIER) procedure (untreated sections produced no signal). The HIER was performed in a basic Tris–HCl buffer solution (10 mM, pH 9.0), with 30 s-microwave pulses at 30–50% power for 4–6 min (to prevent section disintegration). The sections were rinsed in PBS, blocked in 2% normal donkey serum (NDS) for 30 min at room temperature, and incubated on a shaker in rabbit anti-Cre IgG (1:1000; Abcam #ab216262) with 2% NDS and 0.5% Triton X-100 (TX) in PBS for 3 days at 4°C. They were rinsed in PBS three times (10 min each) and incubated in Cy3-conjugated donkey anti-rabbit IgG (1:200; Jackson ImmunoResearch #711-165-152) with 2% NDS and 0.3% TX for 90 min at room temperature. The sections were rinsed in PBS three times (10 min each), mounted onto gelatin/chromium-subbed glass slides, allowed to air-dry, and coverslipped with ProLong Gold Antifade Mountant with no DAPI (ThermoFisher Scientific, #P36930).

2.4.2. PhiYFP(Y65A) IHC

Sections were rinsed in PBS, blocked in 2% normal donkey serum (NDS) for 30 min at room temperature, and incubated on a shaker in rabbit anti-PhiYFP(Y65A) IgG (1:500; kindly provided by Dr. Dawen Cai) with 2% NDS and 0.5% TX in PBS for 3 days at 4°C. They were rinsed in PBS three times (10 min each) and incubated in Cy3-conjugated donkey anti-rabbit IgG (1:400; Jackson ImmunoResearch #711-165-152) with 2% NDS for 90 min at room temperature. The sections were rinsed in PBS three times (10 min each), mounted onto gelatin/chromium-subbed glass slides, allowed to air-dry, and coverslipped with ProLong Gold Antifade Mountant with no DAPI.

2.4.3. EGFP IHC

Sections were rinsed in PBS, blocked in 2% normal donkey serum (NDS) for 30 min at room temperature, and incubated on a shaker in rabbit anti-GFP IgG (1:500; Abcam #ab6556) with 2% NDS and 0.3% TX in PBS for 2–3 days at 4°C. They were rinsed in PBS three times (10 min each) and incubated in AlexaFluor 488-conjugated donkey anti-rabbit IgG (1:1000; ThermoFisher #A21206) with 2% NDS for 90 min at room temperature. The sections were rinsed in PBS three times (10 min each), mounted onto gelatin/chromium-subbed glass slides, allowed to air-dry, and coverslipped with ProLong Gold Antifade Mountant with DAPI (ThermoFisher Scientific, #P36931).

2.4.4. Brainbow IHC

Sections were rinsed in PBS, blocked in 2% normal donkey or goat serum (NDS or NGS, matching the host of the secondary antibodies) for 30 min at room temperature, and incubated on a shaker in rabbit anti-GFP IgG (1:500; Abcam #ab6556), rat anti-mTFP IgG (1:500; Ximbio #155264), and guinea pig anti-TagRFP IgG (1:500; Ximbio #155267), with 2% NDS or NGS and 0.3% TX in PBS for 3 days at 4°C. They were rinsed in PBS three times (10 min each) and incubated in either a donkey or goat set of secondary antibodies diluted in PBS with 2% NDS or NGS, respectively, for 90 min at room temperature. The donkey set consisted of AlexaFluor 488-conjugated donkey anti-rabbit IgG (1:1000; ThermoFisher #A21206), AlexaFluor 594-conjugated donkey anti-rat IgG (1:500; Jackson ImmunoResearch #712-585-150), and AlexaFluor 647-conjugated donkey anti-guinea pig IgG (1:500; Jackson ImmunoResearch #706-605-148). The goat set consisted of AlexaFluor 488-conjugated goat anti-rabbit IgG (1:1000; ThermoFisher #A11034), AlexaFluor 546-conjugated goat anti-rat IgG (1:1000; ThermoFisher #A11081), and AlexaFluor 647-conjugated goat anti-guinea pig IgG (1:1000; ThermoFisher #A21450). The sections were rinsed in PBS three times (10 min each) and stored in PBS at 4°C or immediately mounted onto gelatin/chromium-subbed glass slides, allowed to air-dry, and coverslipped with ProLong Gold Antifade Mountant with no DAPI. Initially, the primary antibody set also included chicken anti-mCherry IgY (1:500; Ximbio #155259) and the secondary antibody set included either AlexaFluor 594-conjugated donkey anti-chicken IgY (1:500; Jackson ImmunoResearch #703-585-155) or AlexaFluor 546-conjugated goat anti-chicken IgY (1:1000; ThermoFisher #A11040). However, this fourth signal was eventually eliminated because of its relative weakness and overlap with another imaging channel. A similar reduction of fluorophore channels from four to three was used by the authors of the Brainbow system (Cai et al., 2013). The two secondary antibody sets produced comparable results.

2.5. Microscopy imaging

Epifluorescence imaging was performed in one or two of three channels (Cy3, GFP, and DAPI) on a Zeiss AxioVision Z1 using a 5× objective (NA 0.13) or a 10× objective (NA 0.45).

Confocal imaging was performed in two channels (AlexaFluor 488 and DAPI) or three (Brainbow) channels (AlexaFluor 488, AlexaFluor 594 or 546, and AlexaFluor 647) on a Leica SP8 resonant scanning confocal system. High-power images were obtained using a 63 × oil objective (NA 1.40) with the XY-resolution of 60–70 nm/pixel and the Z-resolution of 300 nm/optical section. Typical z-stacks consisted of around 100 optical sections. The figures show maximum-intensity projections.

2.6. Computational analyses

In simulations, sample unit-vectors drawn from the von Mises-Fisher (vMF) directional probability distribution were generated using the Wood algorithm (Wood, 1994; Hoff, 2009). The vMF distribution is parametrized by the concentration parameter (κ > 0), which is analogous to the inverse of the variance of the normal distribution (Janušonis and Detering, 2019). Briefly, the algorithm first generates a sample vector that has the required distribution around the fixed mean direction of ${\mu }^0=\left(0,0,1\right)$ and then rotates the vector to any other mean direction $\mu =\left({\mu }_1,{\mu }_2,{\mu }_3\right)$ (where μ is again a unit-vector). In the first step, (i) a 2D-sample unit-vector $v$ is generated from the uniform distribution on the circle, (ii) a sample scalar w is generated from $\left(-1,1\right)$ from the probability distribution that is proportional to $e^{\kappa w}$ (e.g., using rejection sampling), and (iii) the 3D-sample vector $u^0$ is then constructed as $u^0=\left(v\sqrt{1-w^2},w\right)$ . In the second step, the obtained vector is rotated to the required population mean with ${\rm M}{u}^0$ , where the orthogonal 3 × 3 matrix ${\rm M}$ can be given by

${\rm M}=\left(\begin{array}{ccc}\hfill {\mu }_3+\frac{{{\mu }_2}^2}{1+{\mu }_3}\hfill & \hfill -\frac{{\mu }_1{\mu }_2}{1+{\mu }_3}\hfill & \hfill {\mu }_1\hfill \\ \hfill -\frac{{\mu }_1{\mu }_2}{1+{\mu }_3}\hfill & \hfill 1-\frac{{{\mu }_2}^2}{1+{\mu }_3}\hfill & \hfill {\mu }_2\hfill \\ \hfill -{\mu }_1\hfill & \hfill -{\mu }_2\hfill & \hfill {\mu }_3\hfill \end{array}\right)\text{.}$

Fibers can then be simulated with repeated sampling, in which the direction of each step has the vMF distribution with a constant κ and a mean direction iteratively defined as the direction of the previous step.

Conversely, if the directions of all steps are already known (e.g., from experimental data), the fiber can be iteratively rotated such that the current step direction always coincides with ${\mu }^0=\left(0,0,1\right)$ , and then the direction of the next step can be determined with respect to this standardized direction. If the original next direction is given by $u$ , the standardized next direction is given by ${{\rm M}}^{-1}u$ , where the elements of $M$ are constructed from the original current direction $\mu$ . In particular, if the next direction is the same as the current direction, one obtains ${{\rm M}}^{-1}\mu ={\mu }^0$ , as expected. After all directions have been standardized to the corresponding previous directions, this unit-vector sample is used to estimate κ. Several methods are available that are based on the length of the mean vector of the sample ( $L$ , where $0\le L\le 1$ ). If all steps maintain the same direction (i.e., the fiber is “infinitely” rigid), $L=1$ . If, conversely, all directions are equally possible (including the unrealistic return to the previous point), $L=0$ (in the limit). Therefore, higher $L$ values should correspond to higher κ values. A precise estimate of κ for any $L$ can be obtained by using a method based on a fixed-point iteration (Tanabe et al., 2007). More conveniently, a relatively accurate κ estimate can also be obtained with the formula (Tanabe et al., 2007)

$\frac{L\left(3-L^2\right)}{1-L^2}\text{.}$

The reported κ values were calculated with this method. It should be noted that, if $L\ge 0.9$ , an even simpler formula can be used, $1/\left(1-L\right)$ (Mardia and Jupp, 2000). It follows that this formula can be used to estimate κ values that are greater than 10 (which makes it highly applicable in studies of serotonergic fibers). Both formulas assume that $L$ has been calculated in the three spatial dimensions; however, L-based estimators for lower and higher dimensions are available (Mardia and Jupp, 2000; Tanabe et al., 2007). The toolbox of directional statistics is rapidly expanding and includes goodness-of-fit and similarity tests, regression models, and other methods analogous to but mathematically different from those in classical statistics (Pewsey and Garcia-Portugués, 2021).

All simulation and analysis scripts were written in Wolfram Mathematica 13.

3. Results

3.1. The isolation of single fibers with EGFP labeling

The induction of Cre-recombination is an inherently stochastic event at the single-cell level, but the efficiency of this process can be partially controlled with the dose and administration regimen of the inducing agent (tamoxifen) (Zhang et al., 2005; Yoshinobu et al., 2021). Once the production of the fluorophore has been initiated, it can visualize all cell processes, especially if the fluorophore targets the cell membrane (Muzumdar et al., 2007; Cai et al., 2013). These properties can support the visualization of sparse, theoretically random subsets of neurons and their axons (Economo et al., 2016), with each axon fully labeled (given a sufficient transport time).

Following the induction of the Tph2-dependent Cre-recombination, strong EGFP expression was observed in the rostral raphe nuclei. EGFP-positive somata were detectable in unstained sections (Figure 1A); the signal was further verified and enhanced with GFP-immunohistochemistry (Figures 1B,C). After a relatively short post-induction time (less than 2 weeks), uninterrupted single fibers could be visualized in various brain regions (Figures 1D–K). Some of these fibers could be traced over considerable distances (Figure 1D) and contained rich statistical information about the underlying stochastic process. These fibers were used in the random-walk modeling described further. Some fibers had growth cone-like swellings (Figures 1E,F) and resembled developing serotonergic axons in the embryonic mouse brain (Hingorani et al., 2022). However, these structures require further analyses because their functional identity and terminal location cannot be confirmed in this preparation (e.g., the fiber may continue in an adjacent section). Some fibers in several regions (e.g., the inferior colliculus and the dentate gyrus) contained fibers with multiple large swellings (Figures 1G,H). They resembled spheroids, abnormal varicosities associated with dystrophic serotonergic fibers (Daubert et al., 2010). These segments may be in the process of natural degeneration, also because their general morphology resembles that of collapsing axons disconnected from their soma (Shaw and Bray, 1977; Baas and Heidemann, 1986), but they may also represent fibers making extension or branching decisions. True branching events were also easily detected in the tissue (Figures 1I–K). Some branching patterns were nearly indistinguishable from those observed with high resolution in primary serotonergic neuron cultures (Hingorani et al., 2022). Overall, the transgenic model achieved reliable isolation of single serotonergic fibers in a number of brain regions.

FIGURE 1

Figure 1. Serotonergic somata and fibers in two mice with serotonergic neuron-specific EGFP expression (Tph2-iCreER; ROSA^mT/mG). These mice were allowed to survive for 10 days after the tamoxifen treatment and were 5 months old at the time of the tissue collection. (A) An epifluorescence image of unstained EGFP-positive somata (green) in the dorsal raphe (DR). (B,C) Epifluorescence images of GFP-positive somata (green) visualized with GFP-immunohistochemistry in the rostral (B) and caudal (C) DR. Cell nuclei are blue (DAPI). (D–K) Confocal images of GFP-positive fibers (green) visualized with GFP-immunohistochemistry in the superior colliculus (D,E), the inferior colliculus (G,I), and the polymorphic layer of the dentate gyrus of the hippocampus (F,H,J,K). Cell nuclei are blue (DAPI). The images show an isolated fiber (D; arrows), growth cone-like dilated structures (E,F; asterisks), varicosity-like structures (G,H), and branching (I–K; asterisks). The 3D-connectivity of the branching points was verified in the corresponding confocal z-stacks. Some elements (arrowheads) are shown enlarged in insets. Aq, aqueduct; c, capillary; mlf, medial longitudinal fasciculus. Scale bars = 200 μm in (A,B), 100 μm in (C), and 10 μm in (D) [applies to (D–K)].

3.2. The isolation of single fibers with Brainbow labeling

In high-throughput applications, reliable isolation of single fibers should be combined with the visualization of many fibers in the same volume. Simultaneous visualization of many fibers is also important in studies of interactions among different fibers, including the extent of their spatial intermingling (e.g., a small tissue volume may contain fiber segments from many individual fibers or may be dominated by segments from the same few, tortuous fibers). A technical solution is offered by the Brainbow 3.2 toolbox (Cai et al., 2013) which is uniquely well suited for analyses of serotonergic fibers. To our knowledge, it has not been used in this system; an early version of Brainbow (Brainbow 1.0) has been used to visualize serotonergic somata in raphe nuclei (Weber et al., 2009), but it does not support neurite labeling.

In the initial approach, we induced the Tph2-dependent Cre-recombination in mice carrying the Brainbow 3.2 transgene (genotyped offspring of Tph2-iCreER and Brainbow crosses) and used immunohistochemistry to visualize the expression of the Brainbow fluorophores in the rostral raphe nuclei. We observed no Brainbow signal, which was likely caused by a low or negligible expression of Thy1 in the raphe region of this transgenic line, perhaps due to a position effect (Dr. Dawen Cai, personal communication). The Thy1 gene supports strong transgene expression in many, but not all, neuron types (Cai et al., 2013), and another study has reported a weak or undetectable brainstem expression of a different transgene driven by the Thy1 promoter (Dana et al., 2018).

To directly test this explanation, we performed an immunohistochemical staining for PhiYFP, a reporter of the transcriptional activity of the Brainbow transgene region before Cre-recombination (Cai et al., 2013). Strongly PhiYFP-positive cells were found in the neocortex and the hippocampus, but no immunoreactivity was detected in the rostral raphe nuclei and adjacent midbrain structures (Figure 2). This finding can inform other studies attempting Brainbow 3.2 labeling in this neuroanatomical region.

FIGURE 2

Figure 2. Epifluorescence images of PhiYFP-immunoreactivity in (A) the somatosensory cortex, (B) the hippocampus, and (C) the dorsal raphe (DR) of a mouse with the Brainbow 3.2 transgene. The PhiYFP-signal is absent from the DR. I and VI, cortical layers I and VI; Aq, aqueduct; CA1, field CA1 of the hippocampus; DG, dentate gyrus; DRD, dorsal DR; DRV, ventral DR; mlf, medial longitudinal fasciculus. Scale bar = 100 μm.

We next attempted an alternative approach, by injecting Brainbow AAVs (Cai et al., 2013) directly into the dorsal raphe of Tph2-iCreER mice (with the induction of the Tph2-dependent Cre-recombination 1 week later). An immunohistochemical staining for three Brainbow fluorophores (EGFP, mTFP, and TagRFP) revealed labeled somata that were restricted to the raphe complex and showed varying fluorophore-intensity combinations across individual cells, as expected (Figure 3). Importantly, strong and dense labeling of serotonergic fibers was found in the entire brain of some mice, including the diencephalon and telencephalon (Figure 4). Generally, the Brainbow-labeling of fibers showed a considerable variability across animal cases, likely because of the uncontrolled differences in the stereotaxic injections, the efficiency of the AAV-transfection, the efficiency of the tamoxifen-induced Cre-recombination, and other experimental and physiological factors. Our best case was an adult male with a post-induction (transport) time of 3 months. The reliability of this approach may be further improved; at this time, we recommend maximizing the probability of a successful outcome by using batches of 5–10 animals. It should be noted that a single successfully labeled brain can provide an enormous amount of information about individual serotonergic fibers in all brain regions, including their spatial interplay.

FIGURE 3

Figure 3. Confocal images of the immunoreactivity of the three Brainbow fluorophores (with the red, green, and blue channels merged) in the dorsal raphe (DR) of a Tph2-iCreER mouse (male) that received an intracranial injection of the Brainbow-AAVs into the DR. This mouse was allowed to survive for around 3 months after the tamoxifen treatment and was around 1 year old at the time of the tissue collection. (A) A low-power image of the DR. (B) A high-power image of the lateral DR. (C) A high-power image of the ventral DR. Aq, aqueduct; DRD, dorsal DR; DRL, lateral DR; DRV, ventral DR; mlf, medial longitudinal fasciculus. Scale bars = 150 μm (A), 30 μm (B,C).

FIGURE 4

Figure 4. Confocal images of fibers immunoreactive for the three Brainbow fluorophores (with the red, green, and blue channels merged) in the diencephalon and telencephalon of a Tph2-iCreER mouse (male) that received an intracranial injection of Brainbow-AAVs in the dorsal raphe (DR). This mouse was allowed to survive for around 3 months after the tamoxifen treatment and was around 1 year old at the time of the tissue collection. (A,B) The left and right habenulas in the same section. In this geometrically constrained space, individual fibers appear to produce local high-density islands (pink and green in the images). (C) The basolateral amygdala. (D–G) Four enlarged parts of (C) that show highly tortuous but readily separable trajectories (D,E) and unambiguous branching points (F,G). Scale bars = 20 μm (A–C) and 5 μm (D–G).

Some individual fibers appeared to claim considerable territories in the habenula, a nucleus unusual in its highly constrained geometry (Figures 4A,B). It is possible that fibers can become trapped in this region, producing highly tortuous trajectories as they interact with the physical borders. In contrast, fibers appeared to be well intermixed in the basolateral amygdala (Figure 4C), a region with one of the highest densities of serotonergic fibers in the entire brain (Awasthi et al., 2021). Despite this high density, the Brainbow labeling revealed clearly identifiable individual paths (Figure 4D), including full turns (Figure 4E), as well as unambiguous branching events (Figure 4F), some followed by an immediate intermingling of the branches with fibers of different identities (Figure 4G).

3.3. The quantification of fiber trajectories: theoretical considerations

This study focuses on one essential characteristic of serotonergic fibers, their trajectories in the three-dimensional space. Considering the individual uniqueness of these trajectories, their rigorous quantification requires some theoretical assumptions.

We started with a top-down approach, by producing a series of computer simulations of serotonergic fibers modeled as step-wise random walks based on the von Mises-Fisher (vMF) directional probability distribution (Figure 5). This distribution is parametrized by the unitless concentration parameter (κ) that determines the stiffness of the fiber (low κ values make the fiber more flexible, high κ values make it stiffer).

FIGURE 5

Figure 5. Simulated fibers produced by a step-wise random walk, in which the direction of each step was drawn from the von Mises-Fisher probability distribution with the concentration parameter κ. In (A–C), the diameter of the fibers is 1 μm and each step is 1.5 μm in length. The tortuosity index (t) of each fiber, calculated as the fiber length divided by the Euclidean distance between the fiber ends (Jin et al., 2016; Pratelli et al., 2017), is also shown. (A) Six realizations of a walk with κ = 20 and 30 steps (the fiber length of 45 μm), approximately corresponding to the mean segment length in a 40 μm-thick section (Janušonis et al., 2019). Each fiber is shown in three non-perpendicular orientations that demonstrate caveats of two-dimensional interpretations (the same fiber may appear very different, depending on its orientation with respect to the imaging plane). (B) Six realizations of a walk with κ = 20 and 300 steps (the fiber length of 450 μm). The tortuosity index can be suboptimal in that it can strongly vary in fibers produced by the same stochastic process. (C) Single realizations of walks with six κ values (from 2 to 64) and 300 steps (the fiber length of 450 μm). (D) A simulation of Brainbow-labeled fibers. A set of 300 fibers, each with κ = 15 and 300 steps of 1.5 μm, was computed in a volume matching the z-stack in Figure 4C (185 μm × 185 μm × 21 μm). Each fiber was seeded randomly in the volume (using the uniform distribution in the three dimensions). The segments of the fibers that escaped the volume were considered to be outside the “section” and are not shown. Each fiber was assigned a random color. The fiber radius at each step was independently drawn from the normal distribution with the mean of 0.20 μm and the standard deviation of 0.15 μm. (E) A magnified part of (D).

The selected model was based on our previous theoretical investigations (Janušonis et al., 2019; Janušonis and Detering, 2019). The shown simulations are more physically realistic in that they include an accurate scaling of the fiber diameter (set to 1 μm) with respect to its length and also use a physically-informed step (set to 1.5 μm). The simulations show the typical appearance of fiber segments in 40 μm-thick sections (Figure 5A), computational extrapolations of fibers over longer distances, with the corresponding tortuosity indices (Figure 5B), and the dependence of the fiber appearance on the concentration parameter, again with the corresponding tortuosity indices (Figure 5C). The model performs well in capturing the geometry and inherent variability of single serotonergic fibers. In addition, the simulations demonstrate the limitations of the tortuosity index which can vary strongly in different realizations of the same process (Figure 5B). The model can also produce realistic simulations of Brainbow-labeled fibers (Figures 5D,E).

Since fibers are assumed to advance in steps that maintain the same length (but not the same direction), the step length should be ideally optimized with an unbiased and physically-informed procedure. It is easy to see at the qualitative level that an excessively small step can make a fiber appear stiffer that it actually is, and that an excessively large step may lead to unstable computational estimates, because the step can skip over important direction changes (Figure 6A). This problem can become more severe with noise, especially if the step and the noise signal are comparable in magnitude (Figures 6B,C). One solution is to always indicate the step at which the concentration parameter was estimated [e.g., κ(2 μm) = 30], but such arbitrary decisions can complicate quantitative comparisons among different studies. We used two calculations to investigate the dependence of the concentration parameter on the step length, with the goal of finding its optimal (“natural”) value.

FIGURE 6

Figure 6. Qualitative effects of the sampling step along a fiber trajectory. (A) Fiber trajectories can be described with a natural step that reflects the physical structure of the fibers (left). With no noise, smaller steps (oversampling) will make a fiber appear more rigid and produce higher concentration parameter (κ) values (middle). In contrast, larger steps (undersampling) can make a fiber appear more flexible and produce lower κ values; but it can also result in computationally unstable directional distributions (right). (B) With noise, unnaturally small steps can also become unstable because their order of magnitude may become comparable to that of the noise. (C) This phenomenon is demonstrated with the contour of mountains, where an unnaturally small sampling step may capture the contour of relatively small, irrelevant objects and grossly reduce the estimate of κ (i.e., the contour may appear less rigid than it actually is).

First, we simulated a long fiber, with a known value of the concentration parameter (κ = 20), and computed the estimates of this parameter with the original step (set to one unit), as well as with fractional steps and integer-multiples of the original step (Figure 7A). The fractional and integer-multiple steps modeled excessively small and large steps, respectively. As expected, the real parameter value was recovered at step = 1. However, the parameter was grossly overestimated with the smaller steps and slight underestimated with the larger steps. It should be noted that we did not include extremely large steps (e.g., jumping from the beginning of the fiber to its mid-point), which would make the estimates unstable. Importantly, the plot produced an “inflection” point, corresponding to the correct step length.

FIGURE 7

Figure 7. A quantitative method to obtain an optimal step for the sampling of fiber trajectories. (A) The estimated κ values of a simulated fiber (shown in the inset) with the actual κ = 20. The fiber length is 1,000 steps (one length-unit each). The smaller steps were produced by dividing each original step into 2, 3, … 9, and 10 steps (producing steps of 0.1–0.5 units), and the larger steps were produced by multiplying the original step by factors of 2–10 (producing steps of 2–10 units). Even if the original step were unknown, it could be recovered at the sharp “inflection” of the plot. (B) The estimated κ values of 213 serotonergic fibers in the mouse somatosensory cortex that have been visualized with 5-HT-immunohistochemistry [using a published protocol (Janušonis, 2018)], imaged with confocal microscopy, and manually traced through three-dimensional z-stacks with Simple Neurite Tracer, an ImageJ plugin (the inset shows traced fibers in one z-stack). The plugin automatically selects a step of approximately 0.15 μm (Janušonis and Detering, 2019), which is considerably smaller than a typical fiber diameter. A gradual increase in the step length shows that the original step strongly overestimates κ. (C) The mean estimated κ values of the fiber set in (B), which reveals an “inflection region” at the factor values of around 4–10. Considering the diameter of a typical fiber (around 1 μm), these factors (corresponding to steps of approximately 0.6–1.5 μm) is a well-balanced choice.

Second, we attempted the same approach in a set of actual serotonergic fibers that have been manually traced by several trained individuals blind to the model (Figure 7B). In this case, the “correct” step length was not known but could be recovered in a similar “inflection region.” The plot of the population means (Figure 7C) put it at around 0.6–1.5 μm, which is in register with the typical diameter of serotonergic fibers. We set the step at 0.75–1.50 μm in the following analyses.

3.4. High-precision tracing of single fibers with a novel algorithm

After a fiber has been isolated with transgenic methods and imaged with high resolution in the three spatial dimensions, it has to be converted to an array of XYZ-coordinates (i.e., an N × 3 matrix, where N is the number of points). These coordinates (more precisely, the resultant step-vectors) are sufficient to produce a numerical estimate of the concentration parameter. This estimate becomes more accurate with larger step numbers, with one important caveat: all steps should reliably belong to the same fiber (i.e., the trace cannot accidentally jump to another, adjacent fiber). Identity errors cannot be easily mitigated with larger samples; a single long, accurate trace can produce a more reliable estimate than a set of unreliable traces.

The conversion of visualized serotonergic fibers to trajectory arrays is a non-trivial problem. Manual tracing is highly unreliable because human tracers typically have to move up and down z-stacks (a significant perceptual and memory challenge in noisy images), have to make subjective decisions in fibers whose varicosities alternate with signal interruptions, and cannot maintain a steady level of performance due to fatigue. General-purpose systems, such as Imaris (Bitplane) have no built-in knowledge of the physical properties of serotonergic fibers, which (in our hands) results in frequent incorrect decisions. Deep learning-based segmentation systems can separate serotonergic fiber segments from the background, also in whole-brain images, but they currently cannot reliably trace the individual, uninterrupted trajectories of fibers over useful distances (Friedmann et al., 2020).

We developed an algorithm that overcomes some of these difficulties. It prioritizes accuracy over volume and speed and is designed for high-resolution, three-dimensional confocal images of serotonergic fibers labeled with a single fluorophore. In particular, it is highly appropriate for EGFP-labeled serotonergic fibers in the described transgenic model. In the near future, it can be extended to Brainbow-labeled fibers (with an incorporation of the color dimension). The key components of this algorithm are shown in Figure 8. We next describe it in detail.

FIGURE 8

Figure 8. (A) A seed point (yellow) and the initial search circle around it (with a radius of around 2 μm; left), with a corresponding polar plot of the brightness values on the circle (where larger radii correspond to stronger brightness; right). The polar plot indicates the position of the second point, which initiates the automated tracing procedure. (B) At each current fiber point (an XY-coordinate in a specific optical section), the “flashlight” tracing algorithm scans the intensity of the pixel values in an arc that is centered at this point, with a radius R and the subtended angle [−α, α]. The middle point of the arc (0) corresponds to the “expected” direction of the fiber (marked with the yellow segment), which is based on the current point and the point before it. The values of R and α can be adjusted to reflect the general tortuosity of the fiber. For clarity, the first segment is shown straight and the value of R is greatly exaggerated (it is typically on the order of the fiber width). A specified number of optical sections above and below the current optical sections are also scanned in the same way. (C) If the algorithm fails to detect any intensity peaks in the given subset of optical sections, the step can be automatically extended (the dark blue region). Since this extension increases the risk of false-positives (e.g., it may detect a segment of another fiber in the vicinity), the search angle can be automatically decreased. The extended R is a multiple of the original R (0.5R, 1R, 1.5R, 2R, …); in the figure, it is shown smaller for compactness. (D) An example of a decision made by the tracing algorithm in one step (taken from the tracing sequence of an actual EGFP-labeled fiber). The center plot shows the pixel intensity values in the current optical section (z = 0) in the arc centered at the current XY-position with R = 2 μm and α = 0.6π radians (the current direction is given by α = 0). The plots to the left and right of the middle plot show the corresponding pixel values at the same current XY-position and direction in the three optical sections below and above the current optical section (typically, more adjacent sections are used). In this example, the trace turns slightly to the right and three optical sections down. The trace then continues iteratively with the new XY-position, optical section, and direction.

The z-stack of one fluorophore channel is imported as a series of high-resolution, grayscale TIFF images (with the XY-scaling of around 60 nm per pixel and the distance between two adjacent optical sections of around 300 nm). The images are auto-contrasted (with the brightness range of 0–1) and Gaussian-blurred using a one-pixel radius. The initial seed point is selected on a fiber in a specified (“current”) optical section, manually or automatically, and the second point is automatically detected as the brightest point on the circle whose center is at the seed point and the radius is set at a fixed value (R). The second point becomes the “current point,” and the “current direction” is determined by the vector connecting the two points. This initiates the tracing sequence (Figure 8A).

In the current optical section, the brightness of pixels is examined on an arc that is centered at the current point. The radius (tracing step) and angle of the arc are given by R and a fixed interval [−α, α], respectively, where the angle of 0° corresponds to the current direction (Figure 8B). The number of increments (2n) from −α to α should provide a sufficiently dense coverage of the arc; in particular, the arc increment (α/n) should be less than one-half of the expected fiber width. By sampling the arc, brightness peaks within a fixed brightness interval [B_min, B_max] are detected, but the peaks whose width (spatial extent) at B_min falls below a set threshold (W_min) are not kept. The latter procedure, along with the described Gaussian blur, efficiently eliminates pixels whose strong brightness is noise-related and does not extend to adjacent pixels. In the remaining peaks, the peak closest to the current direction is selected (if there are more than one), and its point with the brightest value is recorded as the next candidate point.

In addition to the current optical section, this procedure is repeated in N sections above and N sections below the current optical section (if they exist), using the same arc. Among the next candidate points across the entire subset of optical sections, the brightest point is selected. It becomes the next current point of the fiber. The next current direction is determined by the vector connecting the XY-coordinates of the two most recent points, and the procedure continues iteratively. It is automatically terminated when brightness peaks are no longer detected or when an edge of the stack is reached. It can also be terminated manually, after a fixed number of steps.

We note that the algorithm is “anisotropic” in that within each optical section it prefers a bright value that is closest to the current direction, but among the optical sections it prefers the brightest value in the set (which may not be closest to the current direction in all sections). This choice was based on tests in actual z-stacks and it reflects an inherent property of confocal imaging: confocal images are sets of 2D-images, not true 3D-images with isotropic resolution.

As described, the algorithm would not be able to handle transient interruptions in signal brightness that accidentally fall on the search arc. To prevent this, the algorithm can automatically extend its search radius to multiples of R (0.5R, 1R, 1.5R, 2R, …), if no brightness peaks are detected across the optical section subset (Figure 8C). The set includes one radius that is smaller than the original radius. However, the increasing extensions (numbered i = 1, 2, 3, 4, …, i_max, where i = 2 corresponds to the original R) increase the risk that the trace will jump to another fiber located close to the fiber that is being traced. Therefore, these extensions can be optionally accompanied by an automatically narrowing search angle, given by $\alpha \left(i\right)={\alpha }_o{e}^{-k\left(i-2\right)}$ , where α₀ is the original α and k is a constant (the “narrowing strength”). The maximal number of extensions should be restricted (i_max) to a small value, especially in dense areas. If the next point is detected with the extensions, in the next step the extended R and the narrowed α revert to their original values.

Since the arc is set in a 2D-plane, the actual distances between adjacent points in the 3D-trace may slightly differ from the original R. Extended R values also contribute to these deviations. However, final traces can be automatically fine-tuned with a 3D-step of a fixed length, by simply sliding along the original trace and re-marking (“nudging”) all points such that they are separated by the same exact Euclidean distance. Since the original steps are small, the recalculated trajectory is typically visually indistinguishable from the original one. However, the constant step improves the reliability of concentration parameter estimates.

In some instances, the algorithm may fail to correctly process a step. It can then be manually forced to bridge the gap, with an automatic continuation. With an optimal parameter set, such instances should be rare. In some cases, these adjustments are necessary to make the trace follow a specific branch, when two options are available. The question of whether either of the branches can be considered a continuation of the original fiber (with regard to its prior history) poses interesting questions in stochastic modeling.

Figure 8D shows an example of the automatic detection of the next trace point in an actual EGFP-labeled serotonergic fiber. The typical parameter values are given in Table 1. The algorithm is currently implemented in Wolfram Mathematica, a leading computer language in functional programming, but it can be moved to other languages such as Python or Java.

TABLE 1

Table 1. Typical values of the tracing algorithm.

3.5. The quantification of fiber trajectories: an experimental example

After the image of a fiber has been converted into an array of XYZ-coordinates, with a constant step in the 3D-space, the trace can be iteratively rotated such that the current step always coincides with the unit-vector (0, 0, 1). The direction of the next step is then recorded as a unit-vector, and the next step becomes the current step (Figure 9). The length of the mean of the collected vector sample (shown in red in Figure 9) is used to calculate an estimate of the vMF concentration parameter (κ), as described in detail in the Materials and Methods.

FIGURE 9

Figure 9. Assuming the fiber is composed of constant-length segments, each of which has a direction drawn from the von Mises-Fisher distribution (with the expected direction of each step given by the direction of the previous segment), the value of the concentration parameter (κ) can be estimated by sequentially rotating each current segment (a unit-vector) to the vertical position (shown just below the plane), recording the direction of the next segment (represented by the unit vector just above the plane), and calculating the length (L) of the mean of all direction vectors (shown in red). If $L>0.9$ , $\kappa =1/\left(1-L\right)$ is a good approximation. In particular, if all direction vectors point in the same direction, L approaches 1 and κ approaches infinity, as expected for a perfectly rigid, straight fiber.

We computed the concentration parameters in a sample of EGFP-labeled fibers in the inferior colliculus of mice with the ROSA^mT/mG transgene, after the Tph2-dependent Cre-recombination (Figure 10). Most of the values clustered around κ = 20 (for the used step of 1.5 μm), but two of them were around κ = 40. Larger samples are required to achieve a better understanding of the distribution of this parameter, which may reveal differences among fibers belonging to transcriptionally different neurons or located in different brain regions. For example, a recent study has found two functionally distinct types of serotonergic fibers in the hippocampus (Luchetti et al., 2020). Such analyses fall outside the scope of this methodological study, but we are currently pursuing some of these directions.

FIGURE 10

Figure 10. (A,B) Confocal images of GFP-immunoreactive fibers in the inferior colliculus of two mice with serotonergic neuron-specific EGFP expression (Tph2-iCreER; ROSA^mT/mG). (C) An automatically traced fiber in the z-stack corresponding to (A). The fiber is located in an area with a relatively high fiber density, but the algorithm can still follow its trajectory because fiber intersections can be disambiguated in the 3D-volume. (D) A set of automatically traced fibers in a relatively noisy z-stack corresponding to (B). In (C,D), each yellow dot marks a detected trajectory point. The tracing quality was verified manually across all optical sections. The estimated κ values (at the step of 1.5 μm) are shown in bold. (E) The X-, Y, and Z-coordinates (as a function of the step number) of the traced fiber in (C). (F) The X-, Y, and Z-coordinates (as a function of the step number) of one of the traced fibers in (D) (marked with an asterisk). In (E,F), the Z-coordinate is constrained within a range of 40 μm (because of the thickness of the physical sections). Scale bar = 20 μm.

3.6. An assessment of the reliability of traces

In order to assess the reliability of traces calculated by the algorithm, we digitally dissected a sample volume from a Brainbow z-stack obtained in the basolateral amygdala (Figure 4C). The volume had a smaller XY-area but the same Z-span (Figure 11A). In each optical section, the three fluorophore channels were merged and converted to grayscale to make “color” information unavailable to the algorithm. Nearly all (20) serotonergic fibers present in the volume were automatically traced, and the continuity of the traces was manually checked across all optical sections, point-by-point. If the algorithm made an incorrect decision, the point was manually corrected and the algorithm was allowed to automatically proceed. In total, 18 of the 637 points required a correction (2.8% points). The final traces are shown superimposed on the original image (Figure 11B; the image is in color for better readability) and without it (Figure 11C).

FIGURE 11

Figure 11. An assessment of the accuracy of the tracing algorithm. (A) The maximum-intensity projection of a subset volume of a z-stack of the basolateral amygdala containing Brainbow-labeled fibers (Figure 4C). Scale bar = 10 μm. (B) Twenty automatically traced fibers (with the step of 0.75 μm and no “color” information). Nine traces required no manual correction of any of the points, six traces required one point to be corrected, three traces required two points to be corrected, and two traces required three points to be corrected. In total, 18 of 637 points (2.8%) had to be corrected. Many of the corrections were due to two fibers crossing in the same optical section (in grayscale images). Single-point errors are often easily detectable in maximum-intensity projections (e.g., the trace may jump to another fiber) and can be corrected before the dataset is used in stochastic analyses. (C) The traces shown in isolation. The mean of the estimated κ values was 26 (at the step of 0.75 μm). (D–F) The maximum-intensity projections of the “red,” “green,” and “blue” channels, respectively. The asterisks (one, two, and three) show three “color-isolated” fibers and their corresponding traces.

The algorithm constructs fiber traces based on grayscale pixel intensities in space. In contrast, Brainbow-labeling assigns each fiber a “color” signature and therefore allows for isolation that is independent of spatial information. It suggests that fiber traces obtained in grayscale-converted 3D-images can be verified against “color-filtered” images of the same brain region. A filter can be produced by restricting each fluorophore channel to a narrow range of pixel values or linking these values proportionally. Some fibers can be isolated with the simple restriction of the signal to one fluorophore channel (Figures 11D−F). Future algorithm versions can incorporate this “color dimension” (in addition to the space dimensions) and further reduce the error rate in point detection.

4. Discussion

We developed an integrated system for analyses of single serotonergic fibers in brain tissue. It includes two transgenic mouse models for the visualization of individual fibers, a mathematical model of fiber trajectories, a high-precision algorithm for tracing fibers in 3D-images, and procedures to estimate a single theoretical parameter from complex fiber paths. This system can be used in many experimental setups, including those that use serotonergic perturbations, and can be readily deployed in any mouse brain regions. The study focused on the rostral raphe nuclei, but the caudal raphe nuclei also express Tph2, making this system extendable to the lower brainstem (with serotonergic fibers descending to the spinal cord). Since both mouse models are based on Tph2-dependent Cre-recombination, the single-fluorophore model may express EGFP in the gut because the Tph2 gene is also active in enteric neurons (Gershon, 2009; Neal et al., 2009; Yabut et al., 2019) and the ROSA^mT/mG construct is transcribed in the gut (Muzumdar et al., 2007). However, the overwhelming amount of peripheral serotonin is synthesized by the gut enterochromaffin cells that express a different gene, Tph1 (Gershon and Tack, 2007). Because the Brainbow AAVs are injected intracranially, the labeled somata should be restricted to the raphe complex in this model.

It should be noted that a number of other Cre models have been produced for studies of serotonergic neurons. Some of them may not label all serotonergic neurons and some may label additional neuron populations (Cardozo Pinto et al., 2019). They can be used in combination with various Cre reporters (Navabpour et al., 2020) and may work well with the Brainbow AAVs.

Our fiber-tracing algorithm seeks to harvest all available spatial information in the 3D-image and reconstructs fiber trajectories in step-by-step walks. It makes a minimal set of physical assumptions (e.g., the fiber diameter cannot exceed a certain value and the fiber cannot make 180° turns) and would perform nearly perfectly in high-resolution images with no noise and no signal interruptions along the fiber. The algorithm is useful because this set of assumptions is often sufficient to obtain accurate traces in images of real neural tissue where these ideal conditions can never be met.

When the imaging resolution is low, the fiber has significant signal interruptions, and/or the noise is high, the algorithm may fail. Importantly, in these cases the image itself may contain insufficient information for a correct decision at a given tracing step. This limitation cannot be overcome by any computational method without additional information or further theoretical assumptions (which are never cost-free). In the current implementation, the user can simply override individual step decisions. Care should be exercised in these instances because human visual perception is heavily assumption-based.

Important additional information can be provided by Brainbow-labeling. For example, a signal gap (Maddaloni et al., 2017; Hingorani et al., 2022) is likely to be followed by several candidate fiber-segments, each of which can be reasonably connected to the current end; however, only one of these segments may have the same “color.” It is also possible that none of the segments has the expected “color,” thus indicating that the fiber simply left the section and there is no actual gap.

If further assumptions are introduced, they should capture a deeper understanding of what serotonergic fibers can and cannot do, at a “zoomed-out” level. Convolutional neural networks (CNNs) are particularly well positioned to accomplish this task. A recently introduced CNN, based on the U-Net architecture, can isolate large-diameter serotonergic fiber segments from the background in light-sheet microscopy images of the entire brain (Friedmann et al., 2020). However, this system operates well below the resolution necessary to capture many serotonergic fibers and is not trained to extract information about fiber continuity. Despite these challenges, future CNNs may be able to automatically reject traces that are accurately built on local pixel information (as implemented in our algorithm) but are still unrealistic when viewed as a whole. It should be noted that CNN-based decisions carry an inherent cost in that they can also reject unusual (e.g., other species’ or pathology-induced) trajectories that the network has never seen. Our algorithm is essentially free of these risks.

In general, the pixel-analytical and artificial neural network (ANN) approaches are complementary in that the former excels at “microvision” (with some nearsightedness) and the latter at “macrovision” (with some farsightedness). This distinction is not strict, and future ANNs may be able to capture the continuous trajectories of individual serotonergic fibers in 3D-images. Success in this area depends on the available corpus of annotated training samples. With regard to serotonergic fibers, such high-quality annotations cannot be produced manually. Brainbow labeling provides an accurate, high-throughput method to support this effort: it effectively sorts image pixels by “color” thus “auto-annotates” individual fibers with no human intervention. Training sets can be further augmented with simulated fibers, computer-generated at various densities in unlimited quantities (e.g., using our model). Once the ANN has acquired knowledge about the normal morphology and behavior of serotonergic fibers, it can be used to automatically analyze tissue samples labeled with a single fluorophore, including human brain biopsies and postmortem samples processed with immunohistochemistry.

We note that currently two mathematical models of serotonergic fibers exist. The first model is based on a step-wise walk guided by the vMF-distribution (Janušonis and Detering, 2019) and was used in this study. The second model is based on fractional Brownian motion (FBM) (Janušonis et al., 2020, 2023). Both models appear highly promising but are not mathematically compatible. One essential difference, important in practical applications, is that the vMF-model is purely spatial and has no time term; in contrast, FBM-paths develop in space as a function of continuous time (notably, with no definable “velocity” at any time point). Since the real-time dynamics are typically unavailable in experimental fiber analyses (in most cases, fibers are studied in fixed preparations), the vMF-model is superior in this case. However, the FBM-model promises to eventually capture the “deep” structure of serotonergic fibers, as they evolve in both space and time. Live-imaging recordings of growing serotonergic fibers are needed to advance these studies; efforts have already been made in this direction (Jin et al., 2016; Hingorani et al., 2022). Interestingly, normal Brownian motion (a special case of FBM) and the vMF-distribution are related (Mardia and Jupp, 2000; Gatto, 2013). It should also be noted that the field of “active Brownian particles” may offer other theoretical tools (Romanczuk et al., 2012). We are actively investigating the applicability of these frameworks in our larger research program.

Directional statistics deals with probability distributions on manifolds (e.g., circles, spheres) rather than familiar Euclidean spaces (e.g., lines, planes) (Mardia and Jupp, 2000). A trivial example is given by probability distributions on the unit circle, where the same point can be referred to as 0 or 360 (in degree units). Therefore, caution should be exercised in applications of standard statistical tests to concentration parameter estimates. Directional methods exist to test unit-vector samples for their consistency with the von Mishes-Fisher distribution (against other probability distributions), to compare the concentration parameter estimates in two and more samples (which is relevant to fiber comparisons), to perform ANOVA-and regression-like tests, and to carry out other procedures analogous to those is classical statistics (Mardia and Jupp, 2000; Pewsey and Garcia-Portugués, 2021). A number of R packages are available for these calculations (Pewsey and Garcia-Portugués, 2021), of which Directional (Tsagris et al., 2023) is particularly useful for three-dimensional analyses. We note that each estimate of the concentration parameter is based on a relatively large sample of step-directions and therefore comes with reliability information (e.g., a value based on 20 steps is less reliable than that based on 100 steps). Several methods have been proposed to compute the confidence intervals of κ estimates (Mardia and Jupp, 2000; Matheson et al., 2023), which can be obtained for each fiber.

Long and accurate fiber traces capture the entire distribution of step-directions and thus accurately estimate the population value of the concentration parameter. This value may in turn control the regional density of serotonergic fibers, with developmental or clinical implications (Janušonis and Detering, 2019). In Autism Spectrum Disorder, both altered serotonergic fiber densities and abnormalities of single serotonergic fibers have been observed (Azmitia et al., 2011). In a mouse model, fluoxetine exposure has had a regional effect both on the morphology of individual serotonergic fibers and their density (Nazzi et al., 2019). Generally, the links between the properties of individual serotonergic fibers and their resultant densities remain poorly understood; our experimental platform directly supports this effort.

The proposed framework also supports the rapidly growing interest in single serotonergic neurons and fibers, far beyond their classical (population-based) conceptualizations. Recent studies include such diverse approaches as single-cell RNA-seq (Okaty et al., 2020), high-resolution analyses in primary brainstem cultures (Hingorani et al., 2022), supercomputing-based modeling (Janušonis et al., 2020, 2023), potential applications in artificial neural networks (Lee et al., 2022), and the development of novel, broadly-applicable models in physics, directly inspired by the properties of the serotonergic axons (Vojta et al., 2020; Saito, 2023; Wang et al., 2023). It is conceivable that these new frameworks will produce new fundamental insights into the brain as a self-organizing system.

Data availability statement

The datasets presented in this study can be found in the Dryad online repository at: https://datadryad.org/stash/dataset/doi:10.25349/D95W36.

Ethics statement

The animal study was approved by UCSB Institutional Animal Care and Use Committee. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

KCM developed all experimental methods and contributed to all theoretical discussions. JHH made major contributions in assisting KCM in all experimental procedures and in the maintenance of transgenic mouse colonies. SJ initiated and supervised the project, performed all computational analyses, and wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the National Science Foundation (grants #1822517 and #2112862 to SJ), the National Institute of Mental Health (#MH117488 to SJ), and the California NanoSystems Institute (Challenge-Program Development grants to SJ). We acknowledge the use of the NRI-MCDB Microscopy Facility and the Leica SP8 Resonant Scanning Confocal Microscope (supported by National Science Foundation MRI grant #1625770).

Acknowledgments

We acknowledge the technical contributions of Kirra Thomas, George Grant, and Marissa Millwater (members of the Janušonis laboratory). This work benefitted from manual fiber tracing performed at the beginning of the project by Kimia Boldaji, Natalie Mohr, Jenna Sanfilippo, and Stephanie Valdez. We thank Benjamin Lopez (UCSB Neuroscience Research Institute) and Jennifer Smith (UCSB California NanoSystems Institute) for their technical support in microscopy imaging and genotyping. We also thank Dawen Cai (University of Michigan) for his advice in the search for solutions to technical challenges in the Brainbow labeling.

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Awasthi, J. R., Tamada, K., Overton, E. T. N., and Takumi, T. (2021). Comprehensive topographical map of the serotonergic fibers in the male mouse brain. J. Comp. Neurol. 529, 1391–1429. doi: 10.1002/cne.25027

PubMed Abstract | CrossRef Full Text | Google Scholar

Azmitia, E. C., Singh, J. S., and Whitaker-Azmitia, P. M. (2011). Increased serotonin axons (immunoreactive to 5-HT transporter) in postmortem brains from young autism donors. Neuropharmacology 60, 1347–1354. doi: 10.1016/j.neuropharm.2011.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Baas, P. W., and Heidemann, S. R. (1986). Microtubule reassembly from nucleating fragments during the regrowth of amputated neurites. J. Cell Biol. 103, 917–927. doi: 10.1083/jcb.103.3.917

PubMed Abstract | CrossRef Full Text | Google Scholar

Baizer, J. S. (2001). Serotonergic innervation of the primate claustrum. Brain Res. Bull. 55, 431–434. doi: 10.1016/s0361-9230(01)00535-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Bennis, M., Gamrani, H., Geffard, M., Calas, A., and Kah, O. (1990). The distribution of 5-HT immunoreactive systems in the brain of a saurian, the chameleon. J. Hirnforsch. 31, 563–574.

PubMed Abstract | Google Scholar

Bhat, S. K., and Ganesh, C. B. (2023). Organization of serotonergic system in Sphaerotheca breviceps (Dicroglossidae) tadpole brain. Cell Tissue Res. 391, 67–86. doi: 10.1007/s00441-022-03709-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, D., Cohen, K. B., Luo, T., Lichtman, J. W., and Sanes, J. R. (2013). Improved tools for the Brainbow toolbox. Nat. Methods 10, 540–547. doi: 10.1038/nmeth.2450

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardozo Pinto, D. F., Yang, H., Pollak Dorocic, I., de Jong, J. W., Han, V. J., Peck, J. R., et al. (2019). Characterization of transgenic mouse models targeting neuromodulatory systems reveals organizational principles of the dorsal raphe. Nat. Commun. 10:4633. doi: 10.1038/s41467-019-12392-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrera, I., Molist, P., Anadon, R., and Rodriguez-Moldes, I. (2008). Development of the serotoninergic system in the central nervous system of a shark, the lesser spotted dogfish Scyliorhinus canicula. J. Comp. Neurol. 511, 804–831. doi: 10.1002/cne.21857

PubMed Abstract | CrossRef Full Text | Google Scholar

Challet, E., Miceli, D., Pierre, J., Repérant, J., Masicotte, G., Herbin, M., et al. (1996). Distribution of serotonin-immunoreactivity in the brain of the pigeon (Columba livia). Anat Embryol (Berl) 193, 209–227. doi: 10.1007/bf00198325

CrossRef Full Text | Google Scholar

Challet, E., Pierre, J., Repérant, J., Ward, R., and Miceli, D. (1991). The serotoninergic system of the brain of the viper, Vipera aspis. An immunohistochemical study. J. Chem. Neuroanat. 4, 233–248. doi: 10.1016/0891-0618(91)90015-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Cooke, P., Janowitz, H., and Dougherty, S. E. (2022). Neuronal redevelopment and the regeneration of neuromodulatory axons in the adult mammalian central nervous system. Front. Cell. Neurosci. 16:872501. doi: 10.3389/fncel.2022.872501

PubMed Abstract | CrossRef Full Text | Google Scholar

Dana, H., Novak, O., Guardado-Montesino, M., Fransen, J. W., Hu, A., Borghuis, B. G., et al. (2018). Thy1 transgenic mice expressing the red fluorescent calcium indicator jRGECO1a for neuronal population imaging in vivo. PLoS One 13:e0205444. doi: 10.1371/journal.pone.0205444

PubMed Abstract | CrossRef Full Text | Google Scholar

Daubert, E. A., Heffron, D. S., Mandell, J. W., and Condron, B. G. (2010). Serotonergic dystrophy induced by excess serotonin. Mol. Cell. Neurosci. 44, 297–306. doi: 10.1016/j.mcn.2010.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Daws, R. E., Timmermann, C., Giribaldi, B., Sexton, J. D., Wall, M. B., Erritzoe, D., et al. (2022). Increased global integration in the brain after psilocybin therapy for depression. Nat. Med. 28, 844–851. doi: 10.1038/s41591-022-01744-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Dean, J. G. (2018). Indolethylamine-N-methyltransferase polymorphisms: genetic and biochemical approaches for study of endogenous N, N-dimethyltryptamine. Front. Neurosci. 12:232. doi: 10.3389/fnins.2018.00232

PubMed Abstract | CrossRef Full Text | Google Scholar

Dean, J. G., Liu, T., Huff, S., Sheler, B., Barker, S. A., Strassman, R. J., et al. (2019). Biosynthesis and extracellular concentrations of N, N-dimethyltryptamine (DMT) in mammalian brain. Sci. Rep. 9:9333. doi: 10.1038/s41598-019-45812-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Economo, M. N., Clack, N. G., Lavis, L. D., Gerfen, C. R., Svoboda, K., Myers, E. W., et al. (2016). A platform for brain-wide imaging and reconstruction of individual neurons. elife 5:e10566. doi: 10.7554/eLife.10566

PubMed Abstract | CrossRef Full Text | Google Scholar

Elowitz, M. B., Levine, A. J., Siggia, E. D., and Swain, P. S. (2002). Stochastic gene expression in a single cell. Science 297, 1183–1186. doi: 10.1126/science.1070919

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedmann, D., Pun, A., Adams, E. L., Lui, J. H., Kebschull, J. M., Grutzner, S. M., et al. (2020). Mapping mesoscale axonal projections in the mouse brain using a 3D convolutional network. Proc. Natl. Acad. Sci. U. S. A. 117, 11068–11075. doi: 10.1073/pnas.1918465117

PubMed Abstract | CrossRef Full Text | Google Scholar

Gagnon, D., and Parent, M. (2014). Distribution of VGLUT3 in highly collateralized axons from the rat dorsal raphe nucleus as revealed by single-neuron reconstructions. PLoS One 9:e87709. doi: 10.1371/journal.pone.0087709

CrossRef Full Text | Google Scholar

Gatto, R. (2013). The von Mises-Fisher distribution of the first exit point from the hypersphere of the drifted Brownian motion and the density of the first exit time. Stat Prob Lett 83, 1669–1676. doi: 10.1016/j.spl.2013.03.010

CrossRef Full Text | Google Scholar

Gershon, M. D. (2009). Enteric serotonergic neurones … Finally! J. Physiol. 587:507. doi: 10.1113/jphysiol.2008.167676

PubMed Abstract | CrossRef Full Text | Google Scholar

Gershon, M. D., and Tack, J. (2007). The serotonin signaling system: from basic understanding to drug development for functional GI disorders. Gastroenterology 132, 397–414. doi: 10.1053/j.gastro.2006.11.002

CrossRef Full Text | Google Scholar

Gu, C. (2021). Rapid and reversible development of axonal varicosities: a new form of neural plasticity. Front. Mol. Neurosci. 14:610857. doi: 10.3389/fnmol.2021.610857

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatada, Y., Wu, F., Silverman, R., Schacher, S., and Goldberg, D. J. (1999). En passant synaptic varicosities form directly from growth cones by transient cessation of growth cone advance but not of actin-based motility. J. Neurobiol. 41, 242–251. doi: 10.1002/(SICI)1097-4695(19991105)41:2<242::AID-NEU7>3.0.CO;2-#

PubMed Abstract | CrossRef Full Text | Google Scholar

Hawthorne, A. L., Wylie, C. J., Landmesser, L. T., Deneris, E. S., and Silver, J. (2010). Serotonergic neurons migrate radially through the neuroepithelium by dynamin-mediated somal translocation. J. Neurosci. 30, 420–430. doi: 10.1523/jneurosci.2333-09.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Hellwig, B., Schüz, A., and Aertsen, A. (1994). Synapses on axon collaterals of pyramidal cells are spaced at random intervals: a Golgi study in the mouse cerebral cortex. Biol. Cybern. 71, 1–12. doi: 10.1007/bf00198906

PubMed Abstract | CrossRef Full Text | Google Scholar

Hendricks, T., Francis, N., Fyodorov, D., and Deneris, E. S. (1999). The ETS domain factor Pet-1 is an early and precise marker of central serotonin neurons and interacts with a conserved element in serotonergic genes. J. Neurosci. 19, 10348–10356. doi: 10.1523/JNEUROSCI.19-23-10348.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

Hingorani, M., Viviani, A. M. L., Sanfilippo, J. E., and Janušonis, S. (2022). High-resolution spatiotemporal analysis of single serotonergic axons in an in vitro system. Front. Neurosci. 16:994735. doi: 10.3389/fnins.2022.994735

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoff, P. D. (2009). Simulation of the matrix Bingham-von Mises-Fisher distribution, with applications to multivariate and relational data. J. Comput. Graph. Stat. 18, 438–456. doi: 10.1198/jcgs.2009.07177

CrossRef Full Text | Google Scholar

Hornung, J. P. (2003). The human raphe nuclei and the serotonergic system. J. Chem. Neuroanat. 26, 331–343. doi: 10.1016/j.jchemneu.2003.10.002

CrossRef Full Text | Google Scholar

Huang, C. X., Zhao, Y., Mao, J., Wang, Z., Xu, L., Cheng, J., et al. (2021). An injury-induced serotonergic neuron subpopulation contributes to axon regrowth and function restoration after spinal cord injury in zebrafish. Nat. Commun. 12:7093. doi: 10.1038/s41467-021-27419-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, B. L., and Azmitia, E. C. (1992). Structure and function of the brain serotonin system. Physiol. Rev. 72, 165–229. doi: 10.1152/physrev.1992.72.1.165

CrossRef Full Text | Google Scholar

Janušonis, S. (2018). Some galeomorph sharks express a mammalian microglia-specific protein in radial ependymoglia of the telencephalon. Brain Behav. Evol. 91, 17–30. doi: 10.1159/000484196

PubMed Abstract | CrossRef Full Text | Google Scholar

Janušonis, S., and Detering, N. (2019). A stochastic approach to serotonergic fibers in mental disorders. Biochimie 161, 15–22. doi: 10.1016/j.biochi.2018.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Janušonis, S., Detering, N., Metzler, R., and Vojta, T. (2020). Serotonergic axons as fractional Brownian motion paths: insights into the self-organization of regional densities. Front. Comput. Neurosci. 14:56. doi: 10.3389/fncom.2020.00056

PubMed Abstract | CrossRef Full Text | Google Scholar

Janušonis, S., Haiman, J. H., Metzler, R., and Vojta, T. (2023). Predicting the distribution of serotonergic axons: a supercomputing simulation of reflected fractional Brownian motion in a 3D-mouse brain model. Front. Comput. Neurosci. 17:1189853. doi: 10.3389/fncom.2023.1189853

PubMed Abstract | CrossRef Full Text | Google Scholar

Janušonis, S., Mays, K. C., and Hingorani, M. T. (2019). Serotonergic axons as 3D-walks. ACS Chem. Neurosci. 10, 3064–3067. doi: 10.1021/acschemneuro.8b00667

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., Dougherty, S. E., Wood, K., Sun, L., Cudmore, R. H., Abdalla, A., et al. (2016). Regrowth of serotonin axons in the adult mouse brain following injury. Neuron 91, 748–762. doi: 10.1016/j.neuron.2016.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Kajstura, T. J., Dougherty, S. E., and Linden, D. J. (2018). Serotonin axons in the neocortex of the adult female mouse regrow after traumatic brain injury. J. Neurosci. Res. 96, 512–526. doi: 10.1002/jnr.24059

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitt, M. M., Tabuchi, N., Spencer, W. C., Robinson, H. L., Zhang, X. L., Eastman, B. A., et al. (2022). An adult-stage transcriptional program for survival of serotonergic connectivity. Cell Rep. 39:110711. doi: 10.1016/j.celrep.2022.110711

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, C., Zhang, Z., and Janušonis, S. (2022). Brain serotonergic fibers suggest anomalous diffusion-based dropout in artificial neural networks. Front. Neurosci. 16:949934. doi: 10.3389/fnins.2022.949934

PubMed Abstract | CrossRef Full Text | Google Scholar

Lesch, K. P., and Waider, J. (2012). Serotonin in the modulation of neural plasticity and networks: implications for neurodevelopmental disorders. Neuron 76, 175–191. doi: 10.1016/j.neuron.2012.09.013

CrossRef Full Text | Google Scholar

Lidov, H. G., and Molliver, M. E. (1982). Immunohistochemical study of the development of serotonergic neurons in the rat CNS. Brain Res. Bull. 9, 559–604. doi: 10.1016/0361-9230(82)90164-2

CrossRef Full Text | Google Scholar

Lillesaar, C. (2011). The serotonergic system in fish. J. Chem. Neuroanat. 41, 294–308. doi: 10.1016/j.jchemneu.2011.05.009

CrossRef Full Text | Google Scholar

López, J. M., and González, A. (2014). Organization of the serotonergic system in the central nervous system of two basal actinopterygian fishes: the Cladistians Polypterus senegalus and Erpetoichthys calabaricus. Brain Behav. Evol. 83, 54–76. doi: 10.1159/000358266

PubMed Abstract | CrossRef Full Text | Google Scholar

Luchetti, A., Bota, A., Weitemier, A., Mizuta, K., Sato, M., Islam, T., et al. (2020). Two functionally distinct serotonergic projections into hippocampus. J. Neurosci. 40, 4936–4944. doi: 10.1523/jneurosci.2724-19.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, D., Deng, B., Sun, C., McComb, D. W., and Gu, C. (2022). The mechanical microenvironment regulates axon diameters visualized by cryo-electron tomography. Cells 11:2533. doi: 10.3390/cells11162533

PubMed Abstract | CrossRef Full Text | Google Scholar

Maddaloni, G., Bertero, A., Pratelli, M., Barsotti, N., Boonstra, A., Giorgi, A., et al. (2017). Development of serotonergic fibers in the post-natal mouse brain. Front. Cell. Neurosci. 11:202. doi: 10.3389/fncel.2017.00202

PubMed Abstract | CrossRef Full Text | Google Scholar

Maia, G. H., Soares, J. I., Almeida, S. G., Leite, J. M., Baptista, H. X., Lukoyanova, A. N., et al. (2019). Altered serotonin innervation in the rat epileptic brain. Brain Res. Bull. 152, 95–106. doi: 10.1016/j.brainresbull.2019.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardia, K. V., and Jupp, P. E. (2000). Directional statistics. New York: John Wiley & Sons.

Google Scholar

Matheson, I. C., Malhotra, R., and Keane, J. T. (2023). A von Mises–Fisher distribution for the orbital poles of the plutinos. Mon. Not. R. Astron. Soc. 522, 3298–3307. doi: 10.1093/mnras/stad1208

CrossRef Full Text | Google Scholar

Morgan, A. A., Alves, N. D., Stevens, G. S., Yeasmin, T. T., Mackay, A., Power, S., et al. (2023). Medial prefrontal cortex serotonin input regulates cognitive flexibility in mice. bio Rxiv. doi: 10.1101/2023.03.30.534775

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzumdar, M. D., Tasic, B., Miyamichi, K., Li, L., and Luo, L. (2007). A global double-fluorescent Cre reporter mouse. Genesis 45, 593–605. doi: 10.1002/dvg.20335

PubMed Abstract | CrossRef Full Text | Google Scholar

Navabpour, S., Kwapis, J. L., and Jarome, T. J. (2020). A neuroscientist's guide to transgenic mice and other genetic tools. Neurosci. Biobehav. Rev. 108, 732–748. doi: 10.1016/j.neubiorev.2019.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Nazzi, S., Maddaloni, G., Pratelli, M., and Pasqualetti, M. (2019). Fluoxetine induces morphological rearrangements of serotonergic fibers in the hippocampus. ACS Chem. Neurosci. 10, 3218–3224. doi: 10.1021/acschemneuro.8b00655

PubMed Abstract | CrossRef Full Text | Google Scholar

Neal, K. B., Parry, L. J., and Bornstein, J. C. (2009). Strain-specific genetics, anatomy and function of enteric neural serotonergic pathways in inbred mice. J. Physiol. 587, 567–586. doi: 10.1113/jphysiol.2008.160416

PubMed Abstract | CrossRef Full Text | Google Scholar

Okaty, B. W., Commons, K. G., and Dymecki, S. M. (2019). Embracing diversity in the 5-HT neuronal system. Nat. Rev. Neurosci. 20, 397–424. doi: 10.1038/s41583-019-0151-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Okaty, B. W., Sturrock, N., Escobedo Lozoya, Y., Chang, Y., Senft, R. A., Lyon, K. A., et al. (2020). A single-cell transcriptomic and anatomic atlas of mouse dorsal raphe Pet1 neurons. elife 9:e55523. doi: 10.7554/eLife.55523

PubMed Abstract | CrossRef Full Text | Google Scholar

Paul, E. R., Schwieler, L., Erhardt, S., Boda, S., Trepci, A., Kämpe, R., et al. (2022). Peripheral and central kynurenine pathway abnormalities in major depression. Brain Behav. Immun. 101, 136–145. doi: 10.1016/j.bbi.2022.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Pewsey, A., and Garcia-Portugués, E. (2021). Recet advances in directional statistics. TEST 30, 1–58. doi: 10.1007/s11749-021-00759-x

CrossRef Full Text | Google Scholar

Pratelli, M., Migliarini, S., Pelosi, B., Napolitano, F., Usiello, A., and Pasqualetti, M. (2017). Perturbation of serotonin homeostasis during adulthood affects serotonergic neuronal circuitry. eNeuro 4, ENEURO.0376–ENEU16.2017. doi: 10.1523/eneuro.0376-16.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Raj, A., and van Oudenaarden, A. (2008). Nature, nurture, or chance: stochastic gene expression and its consequences. Cells 135, 216–226. doi: 10.1016/j.cell.2008.09.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, J., Friedmann, D., Xiong, J., Liu, C. D., Ferguson, B. R., Weerakkody, T., et al. (2018). Anatomically defined and functionally distinct dorsal raphe serotonin sub-systems. Cells 175, 472–487.e20. doi: 10.1016/j.cell.2018.07.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, J., Isakova, A., Friedmann, D., Zeng, J., Grutzner, S. M., Pun, A., et al. (2019). Single-cell transcriptomes and whole-brain projections of serotonin neurons in the mouse dorsal and median raphe nuclei. elife 8:e49424. doi: 10.7554/eLife.49424

PubMed Abstract | CrossRef Full Text | Google Scholar

Romanczuk, P., Bär, M., Ebeling, W., and Lindner, B. (2012). Active Brownian particles. Eur. Phys. J. Special Topics 202, 1–162. doi: 10.1140/epjst/e2012-01529-y

CrossRef Full Text | Google Scholar

Saito, T. (2023). Langevin analogy between particle trajectories and polymer configurations. Phys. Rev. E 107:034502. doi: 10.1103/PhysRevE.107.034502

PubMed Abstract | CrossRef Full Text | Google Scholar

Sako, H., Kojima, T., and Okado, N. (1986). Immunohistochemical study on the development of serotoninergic neurons in the chick: I. distribution of cell bodies and fibers in the brain. J. Comp. Neurol. 253, 61–78. doi: 10.1002/cne.902530106

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaw, G., and Bray, D. (1977). Movement and extension of isolated growth cones. Exp. Cell Res. 104, 55–62. doi: 10.1016/0014-4827(77)90068-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Shepherd, G. M., Raastad, M., and Andersen, P. (2002). General and variable features of varicosity spacing along unmyelinated axons in the hippocampus and cerebellum. Proc. Natl. Acad. Sci. U. S. A. 99, 6340–6345. doi: 10.1073/pnas.052151299

PubMed Abstract | CrossRef Full Text | Google Scholar

Stuesse, S. L., Stuesse, D. C., and Cruce, W. L. (1995). Raphe nuclei in three cartilaginous fishes, Hydrolagus colliei, Heterodontus francisci, and Squalus acanthias. J. Comp. Neurol. 358, 414–427. doi: 10.1002/cne.903580308

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanabe, A., Fukumizu, K., Oba, S., Takenouchi, T., and Ishii, S. (2007). Parameter estimation for von Mises-Fisher distributions. Comput. Stat. 22, 145–157. doi: 10.1007/s00180-007-0030-7

CrossRef Full Text | Google Scholar

Teissier, A., Soiza-Reilly, M., and Gaspar, P. (2017). Refining the role of 5-HT in postnatal development of brain circuits. Front. Cell. Neurosci. 11:139. doi: 10.3389/fncel.2017.00139

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsagris, M., Athineou, G., Adam, C., Sajib, A., Amson, E., and Waldstein, M. J. (2023). Package ‘directional’: A collection of functions for directional data analysis. CRAN Available at: cran.r-project.org/web/packages/Directional

Google Scholar

Ueda, S., Nojyo, Y., and Sano, Y. (1984). Immunohistochemical demonstration of the serotonin neuron system in the central nervous system of the bullfrog, Rana catesbeiana. Anat. Embryol. (Berl) 169, 219–229. doi: 10.1007/bf00315627

PubMed Abstract | CrossRef Full Text | Google Scholar

Vargas, M. V., Dunlap, L. E., Dong, C., Carter, S. J., Tombari, R. J., Jami, S. A., et al. (2023). Psychedelics promote neuroplasticity through the activation of intracellular 5-HT_2A receptors. Science 379, 700–706. doi: 10.1126/science.adf0435

PubMed Abstract | CrossRef Full Text | Google Scholar

Vojta, T., Halladay, S., Skinner, S., Janušonis, S., Guggenberger, T., and Metzler, R. (2020). Reflected fractional Brownian motion in one and higher dimensions. Phys. Rev. E 102:032108. doi: 10.1103/PhysRevE.102.032108

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Balcerek, M., Burnecki, K., Chechkin, A. V., Janušonis, S., Ślęzak, J., et al. (2023). Memory-multi-fractional Brownian motion with continuous correlations. Phys. Rev. Res. 5:L032025. doi: 10.1103/PhysRevResearch.5.L032025

CrossRef Full Text | Google Scholar

Weber, T., Böhm, G., Hermann, E., Schütz, G., Schönig, K., and Bartsch, D. (2009). Inducible gene manipulations in serotonergic neurons. Front. Mol. Neurosci. 2:24. doi: 10.3389/neuro.02.024.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, M. A., and Molliver, M. E. (1991). The organization of serotonergic projections to cerebral cortex in primates: regional distribution of axon terminals. Neuroscience 44, 537–553. doi: 10.1016/0306-4522(91)90076-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wood, A. T. A. (1994). Simulation of the von Mises-Fisher distribution. Commun. Stat. Simula 23, 157–164.

Google Scholar

Yabut, J. M., Crane, J. D., Green, A. E., Keating, D. J., Khan, W. I., and Steinberg, G. R. (2019). Emerging roles for serotonin in regulating metabolism: new implications for an ancient molecule. Endocr. Rev. 40, 1092–1107. doi: 10.1210/er.2018-00283

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshinobu, K., Araki, M., Morita, A., Araki, M., Kokuba, S., Nakagata, N., et al. (2021). Tamoxifen feeding method is suitable for efficient conditional knockout. Exp. Anim. 70, 91–100. doi: 10.1538/expanim.19-0138

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Fujitani, Y., Wright, C. V., and Gannon, M. (2005). Efficient recombination in pancreatic islets by a tamoxifen-inducible Cre-recombinase. Genesis 42, 210–217. doi: 10.1002/gene.20137

PubMed Abstract | CrossRef Full Text | Google Scholar