Displaying the autonomic processing network in humans – a global tractography approach

Regulation of the internal homeostasis is modulated by the central autonomic system. So far, the view of this system is determined by animal and human research focusing on cortical and subcortical grey substance regions. To provide an overview based on white matter architecture, we used a global tractography approach to reconstruct a network of tracts interconnecting brain regions that are known to be involved in autonomic processing. Diffusion weighted imaging data were obtained from subjects of the human connectome project (HCP) database. Resulting tracts are in good agreement with previous studies assuming a division of the central autonomic system into a cortical (CAN) and a subcortical network (SAN): the CAN consist of three subsystems that encompass all cerebral lobes and overlap within the insular cortex: a parieto-anterior-temporal pathway (PATP), an occipito-posterior-temporo-frontal pathway (OPTFP) and a limbic pathway. The SAN on the other hand connects the hypothalamus to the periaqueductal grey and locus coeruleus, before it branches into a dorsal and a lateral part that target autonomic nuclei in the rostral medulla oblongata. Our approach furthermore reveals how the CAN and SAN are interconnected: the hypothalamus can be considered as the interface-structure of the SAN, whereas the insula is the central hub of the CAN. The hypothalamus receives input from prefrontal cortical fields but is also connected to the ventral apex of the insular cortex. Thus, a holistic view of the central autonomic system could be created that may promote the understanding of autonomic signaling under physiological and pathophysiological conditions.


Introduction
Breathing, strength of heart contraction or blood pressure -each aspect of autonomic physiology is tightly controlled by the autonomic nervous system ( Benarroch, 1993 ;Loewy, 1981 ). This is a network consisting of afferent and efferent fibers that modulate the internal state in a specialized but synergistic manner ( Jänig, 2006 ). Whereas postganglionic sympathetic fibers release norepinephrine that activates adrenergic receptors, postganglionic parasympathetic fibers release acetylcholine, a neurotransmitter that acts via muscarinic cholinergic receptors ( Gordan et al., 2015 ). Afferent fibers innervating baroand chemoreceptors provide sensory feedback of relevant actuating variables ( Klabunde, 2012 ). To adapt autonomic physiology to situational needs and external requirements, a central autonomic control is required to orchestrate the outflow of the autonomic nervous system ( Critchley, 2005 ). This is provided by the central autonomic system (CAS) that consists of the following two main components: 1. A cortical autonomic network (CAN) that modulates physiology according to "higher ", i.e. motor, cognitive or emotional demands nerve activity ( Beissner et al., 2013 ;Sklerov et al., 2019 ). Furthermore, using voxel-based morphometry (VBM), a correlation between the aerobic capacity and grey matter density in healthy subjects has been reported ( Peters et al., 2014 ).
Here, we present a white matter architecture-based approach that was devised to gain a holistic view of the entire central autonomic system (CAS). Using known grey matter regions of interest (ROIs) involved in autonomic signaling processes, we performed global tractography based on diffusion-weighted imaging (DWI). DWI sequences were extracted from a cohort of 100 subjects derived from the Human Connectome Project (HCP) database. The resulting processing network covers the cortical (CAN) as well as the subcortical autonomic network (SAN). The insular cortex and the hypothalamus turned out to be central interfaces that link the two subsystems. Tracked autonomic fiber bundles are integrated into current anatomical concepts and their impact on further research strategies becomes discussed.

Research strategy
To define fiber tracts belonging to central autonomic system (CAS), we used a research strategy that warrants an unbiased and at least semiquantitative visualization of particular fiber tracts in a human population ( Hosp et al., 2020 ). In brief: 1. Key anatomical structures involved in autonomic processing were identified via a literature search, and the relevant regions of interest (ROIs) were defined in MNI space. 2. Based on a normative connectome in MNI space ( n = 100 subjects; Human Connectome Project database, www.humanconnectome.org ), ROIs were applied in all possible combinations for fiber selection. Combinations of ROIs that were not connected (that is, those that did not yield any interlinking streamline) were excluded from all further investigations. By definition, a streamline connects a combination of ROIs if any of its supporting points lies within each ROI. 3. Global tractography was performed individually for each subject in an unconstrained manner as described in previous work ( Coenen et al., 2018 ;Hosp et al., 2019 ). As opposed to a local walker-based tractography, global fiber tracking has several advantages: it creates a fiber configuration that delivers the best explanation for the acquired diffusion-weighted MRI data, it is usually less sensitive to noise and fiber density is directly related to the measured data itself ( Fillard et al., 2011 ;Reisert et al., 2012Reisert et al., , 2011Schumacher et al., 2018aSchumacher et al., , 2018b. The above-defined ROIs were then warped from MNI to native (i.e. single subject) space using the deformation fields provided by CAT12 to select streamlines based on each subject's individual connectome. Thus, the tracking procedure itself occurred independently from ROIs, thereby abolishing the bias of predetermined starting and stopping points. 4. Fiber density maps and directional fiber density maps were computed for each of the bundles selected. After rendition in native space, density maps were normalized to template space and aggregated to yield a group density representation of the bundles in MNI standard space. 5. Deterministic bundle-specific tractography based on the tensor fields was performed for the visualization of particular fiber tracts.

Subjects and magnetic resonance imaging
We used diffusion data derived from 100 subjects from the Human Connectome Project (HCP; Q1, S3) data corpus (resolution 1.25 mm isotropic, three b-shells with 1000, 2000, 3000; for more details on the protocol and pre-processing, see ( Glasser et al., 2013 ), mean age 29 ± 3.7; 64 females, 36 males). For mapping into MNI template space, all of the accompanying T1 weighted images were subjected to CAT12 segmentation/normalization. For normalization to MNI space, we used the deformation fields provided by CAT12.

Definition of regions of interest (ROIs)
The key structures involved in central autonomic processing were defined by a literature search using the NIHS Public Archive for the Refereed Literature (PUBMED; https://www.nvbi.nlm.nih.gov ). The following search criteria were used: "central autonomic processing " and "central autonomic network ". Review articles published within the last 10 years were taken into account. At the 19.10.2020, the term "central autonomic processing " yielded 263 results, whereas the term "central autonomic network " yielded 195 partially overlapping results. Resulting papers were screened thoroughly to select altogether 21 papers that crucially focused on the anatomy of central networks related to autonomic function (summarized in Supplementary Table1). Whereas an involvement of some "higher " brain regions (e.g. prefrontal cortex) in autonomic regulations may be unique to humans ( Riganello et al., 2019 ), the basic anatomy and organization of subcortical autonomic networks are considered to be conserved among mammal species ( Benarroch, 1993 ;Silvani et al., 2016 ). Thus, research from animal models were also included into the literature selection. Based on the literature research, we identified multiple ROIs for streamline selection. With respect to cortical ROIs, masks were adapted from the Desikan-Killiany atlas (DKA; Kanaan et al., 2016 ). In detail, the following structures were selected: amygdala (AMG), insular cortex (IC), lingual gyrus (LG), medial prefrontal cortex (mPFC; i.e. medial and lateral orbitofrontal cortex according to the DKA nomenclature), midcingulate cortex (MCC, i.e. posteriorcingulate cortex according to the DKA nomenclature), precuneus (Pc), rostral anteriorcingulate cortex (rACC) and supramarginal gyrus (SMG). For the subcortical part, the following ROIs were selected: the hypothalamus (HT; mask adapted from Ilinsky et al. 2018), the periaqueductal grey (PAG; mask adapted from Edlow et al., 2012 ), the locus coeruleus (LC; mask adapted from Edlow et al., 2012 ). With respect to the autonomic centers of the medulla oblongata, two seeds were defined (Supplementary Figure 1): a dorsal medullar seed (DMS; ± 4, − 42, − 51; radius 3 mm) covers the dorsal nucleus of the vagus nerve (DNX), the nucleus ambiguus (NAmb) and the nucleus of the tractus solitarius (NTS; Paxinos and Huang 1995 ). DNX, NAmb and NTS are involved in the processing of cardiac and vascular afferent baroceptor feedback ( Kaufman and Jones, 2018 ). Furthermore, a lateral medullar seed (LMS; ± 6, − 38, − 51; radius 3 mm) covers the ventrolateral medulla (VLM; Paxinos and Huang 1995 ). The VLM harbors efferent sympathoexcitatory premotor neurons that trigger the systemic sympathetic response via activation of the intermediolateral nucleus of the thoracic spinal cord ( Mai and Paxinos, 2012 ).

Global tractography algorithm
Global tractography was performed individually for each subject as described ( Coenen et al., 2018 ;Hosp et al., 2019 ) and evaluated ( Fillard et al., 2011 ) previously. White-matter probability maps obtained from CAT12 were thresholded at a probability of 0.5 to determine the area of fiber reconstruction (Supplementary Figure 2). For tractography analysis, we followed the global approach described by Reisert et al. (2011 ; see https://www.uniklinik-freiburg.de/mren/research-groups/diffperf/fibertools.html ). As opposed to local walker-based tractography, global fiber tracking aims to find a fiber configuration that delivers the best explanation for the acquired diffusion-weighted MRI data. The optimization process is essentially similar to a polymerization process: the streamlines are initially short and fuzzy, whereas during the optimization phase, the connections proliferate, and fibers become increasingly congruent with the data. The algorithm is based on so-called "simulated annealing ". The system is initially set at a relatively high temperature, which is slowly decreased during iterations to progressively obtain more accurate results. Global fiber tractography is usually found to be less sensitive to noise, and fiber density is directly related to the measured data itself. From the two parameter sets provided by the toolbox Fig. 1. The parieto-anterior-temporal pathway (PATP): A. Top, ROIs used for streamline selection are framed by a solid line. Dashed lines indicate regions that are targeted by the tract but were not used for initial streamline selection. Please note that SMG and Pc were fused to a combined ROI to preserve the dichotomic ramification of streamlines at the level of posterior insula. Bottom, the scheme indicates the three-dimensional distribution of ROIs within the brain. B. An overview of the PATP produced by deterministic bundle-specific tractography: streamlines are attached to the precuneus (Pc) or the supramarginal gyrus (SMG) and project to the insular cortex (IC) and amygdala (AMG). From the bottleneck between IC and AMG, streamlines project to hypothalamus (HT), cross the midline via the anterior commissure or project to the temporal pole (TP) and the anterior part of the superior temporal gyrus (STG). C. To show the distribution of fibers at the amygdalo-insular junction, the group density representation map of the PATP in MNI space is superimposed onto a coronal T1w template. Color-coding in red (left) and blue (right) indicates the probability of occurrence of fiber streamlines in the entire group (in percentage). Brown shading denotes the amygdala, ocher shading denotes the insular cortex. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) , we chose the 'dense' preset for our analyses. In addition, to increase reproducibility, we increased the number of fibers using the following accumulation strategy: after the cooling-down phase, the temperature was again increased to 0.1, and the state further iterated for 10 7 iterations. This procedure was iterated over 5 rounds and the tracts resulting from each round were accumulated to obtain one final tractogramm that was five times larger than the initial one. This approach was proposed by ( Schumacher et al., 2018a ) and showed much higher re-test reliability.

Fiber density generation and aggregation
Fiber densities were computed by means of trilinear interpolation on an isotropic matrix with 1.5 mm resolution. Before averaging, the densities were thresholded at a cut-off value of 1 mm streamline length per voxel. Group averages of the streamline indicator images were built, and a color-coded scheme was used to indicate the voxel-wise probability (in percentage) of fiber streamline occurrence in the entire group. To understand the true extension of the tracts, the structure was overlaid onto a T1W template in MNI space (Supplementary Figure 3A). Directional fiber density maps were then obtained by rendering the rank-1 tensor formed by the tangent of the fibers. The tensor field representation allowed the calculation of mean values in the common additive manner, as for the scalar densities. The directional density maps were normalized as explained above. However, since the tensorial nature of the field had to be taken into account for normalization to MNI standard space, we used the Jacobian matrix of the associated template warp to map the tensor from subject space to MNI standard space. The resulting tensor fields in MNI standard space were then used for deterministic bundle-specific tractography (Supplementary Figure 3B). This was obtained by randomly placing seeds in high density regions (threshold > 10ˆ − 1), with a very loose stopping criterion (threshold > 10ˆ − 8) to allow inclusion of cortical areas. This streamlining algorithm is similar to the commonly-used FACT algorithm ( Mori et al., 1999 ;Mori and van Zijl, 2002 ). The medical imaging platform NORA was used for visualization and bundle specific tractography ( www.nora-imaging.org ). The resulting tracts varied in terms of robustness and reproducibility as indicated by the probability of fiber occurrence. A low probability of fiber occurrence can be explained by either high degree of inter-individual variability in tract anatomy, a low absolute number of streamlines, or a combination of both factors. To exclude tracts that were especially weak or variable, only those that included voxels with at least ≤ 40% probability of streamline occurrence were taken into account.

Results
Based on the literature, the central autonomic system can be categorized into a cortical autonomic network (CAN) that mediates modulation by higher order processes and a subcortical autonomic network (SAN) that connects medullar autonomic centers and brainstem nuclei to hypothalamus or insula ( Beissner et al., 2013 ;Benarroch, 1993 ;Dampney, 2016 ). Although the selection of ROI-combinations for the tracking process followed an unconstrained and random approach, resulting tracts fitted well into this arrangement. Thus, we decided to adopt this classification for the presentation of our results: 1 The cortical autonomic network (CAN): the CAN is a composite of three distinct pathways: the parieto-anterior-temporal, the occipitoposterior-temporo-frontal and the limbic pathway.
The parieto-anterior-temporal pathway (PATP) is a complex mesh that connects precuneus (Pc) and the supramarginal cortex (SMC) with the anterior insular cortex (IC) and the amygdala (AMG). From here, fibers project to the temporal pole (TP) and the anterior part of the superior temporal gyrus (STG) -but also to the hypothalamus (HT) and over the anterior commissure. To preserve the dichotomic ramification of streamlines between precuneus Dashed lines indicate regions that are targeted by the tract but were not used for streamline selection. Bottom, the scheme indicates the three-dimensional distribution of ROIs within the brain. B. An overview of the OPTFP produced by deterministic bundle-specific tractography: streamlines are attached to the lingual gyrus (LG) and project through the insular cortex (IC) and amygdala (AMG) to the medial prefrontal cortex (mPFC). At the level of the amygdalo-insular junction, streamlines project to posterior middle-and inferior temporal gyrus (MTG/ITG). The left-sided tract is more pronounced than the right one. C. To display the temporal projections at the amygdaloinsular junction, the group density representation map of the OPTFP in MNI space is superimposed onto a coronal T1w template. Colorcoding in red (left) and blue (right) indicates the probability of occurrence of fiber streamlines in the entire group (in percentage). Brown shading denotes the amygdala, ocher shading denotes the insular cortex. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) and supramarginal gyrus, both ROIs were integrated into a combined one. An overview of the entire pathway obtained from deterministic bundle-specific tractography is displayed in Fig. 1 A and B. The bottleneck region between anterior IC and AMG (i.e. the amygdalo-insular junction) is characterized by a fan-shaped distribution of fibers that is shown in Fig. 1 C. Group density representation maps underlying the bundle-specific tractography are presented in the Supplementary Figure 4. The occipito-posterior-temporo-frontal pathway (OPTFP) connects the lingual gyrus (LG) to the medial prefrontal gyrus (mPFC) via insular cortex (IC) and amygdala (AMG). An overview of the entire pathway obtained from deterministic bundle-specific tractography is displayed in Fig. 2 A and B. The left-sided tract is more pronounced than the right one. At the level of the amygdaloinsular junction, streamlines project to the posterior middle-and inferior temporal gyrus (MTG/ITG) as shown in Fig. 2 C. Group density representation maps underlying the bundle-specific tractography are presented in Supplementary Figure 5. The limbic pathway can be divided into a rostral-anterior and a midcingulate subsystem. Similar to the OPTFP, the rostral-anterior part connects the rostral anterior-cingulate cortex (rACC) to the insula and the lingual gyrus (LG). An overview of the entire pathway obtained from deterministic bundle-specific tractography is displayed in Fig. 3 A and B. Group density representation maps underlying the bundle-specific tractography are presented in Supplementary Figure 6. The mid-cingulate limbic pathway is a complex looped fiber-tract that resembles to the circuit of Papez ( Bhattacharyya, 2017 ): projections from the posterior mid-cingular cortex (MCC) target the inferior insular cortex (IC) and amygdala (AMG). From the amygdalo-insular junction, fibers target the superior temporal gyrus (STG) and project over the anterior commissure. Fibers further form a loop in parallel to the cingulate cortex and finally target the parahippocampal gyrus (GPH), the hippocampal formation (HF) and the AMG. An overview of the entire pathway obtained from deterministic bundle-specific tractography is shown in Fig. 4 A, B and C. Group Dashed lines indicate regions that are targeted by the tract but were not used for streamline selection. Bottom, the scheme indicates the three-dimensional distribution of ROIs within the brain. B -D. An overview of the mid-cingulate limbic pathway produced by deterministic bundle-specific tractography: projections from the posterior mid-cingular cortex (MCC) target the inferior insular cortex (IC) and amygdala (AMG). From the amygdalo-insular junction, fibers target the superior temporal gyrus (STG) and project over the anterior commissure (AC). Furthermore, fibers form a loop in parallel to the cingulate cortex including the rostral anteriorcingulate cortex (rACC) and finally target the parahippocampal gyrus (GPH), the hippocampal formation (HF) and the AMG. Interestingly, the left-sided tract is more pronounced than the right one. E. To display the temporal projections at the amygdalo-insular junction, the group density representation map in MNI space is superimposed onto a coronal T1w template. Color-coding in red (left) and blue (right) indicates the probability of occurrence of fiber streamlines in the entire group (in percentage). Brown shading denotes the amygdala, ocher shading denotes the insular cortex. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) density representation maps underlying the bundle-specific tractography are presented in Supplementary Figure 7.
1 The subcortical autonomic network (SAN): the subcortical autonomic network couples the autonomic centers of the medullar reticular formation to brainstem nuclei and higher-order regulatory hubs such as the hypothalamus and insular cortex. It can be subdivided into a dorsal and a lateral pathway that overlap in large part.
The dorsal column connects the dorsal medullar seed (DMS) to the locus coeruleus (LC), the periaqueductal grey (PAG) and the hypothalamus (HT). Streamlines can be further tracked to the medial prefrontal cortex (mPFC). An overview of the entire pathway obtained from deterministic bundle-specific tractography is displayed in Fig. 5 A. Group density representation maps underlying the bundle-specific tractography are presented in Supplementary  Figure 8. In addition, a lateral column exists that connects the lateral medullar seed (LMS) with the locus coeruleus, and then unites with the dorsal column . An overview of both pathways and their anatomic relation is shown in Fig. 5 B. Group density representation maps underlying the bundle-specific tractography are presented in Supplementary Figure 9. Whereas the dorsal and lateral columns are characterized by their connection to periaqueductal grey, hypothalamus and medial prefrontal cortex, projections from the antero-and postero-dorsal insula directly target the locus coeruleus and furthermore extend into the lateral or dorsal medullar seed. An overview of theses insular connections obtained from deterministic bundle-specific tractography is displayed in Fig. 6 A. Group density representation maps underlying the bundle-specific tractography are presented in Supplementary Figure 10. Fig. 6 B ultimately displays the topographical relationship between the insular connections the dorsal-and ventral column of the subcortical autonomic system.
1 The organigram of the central autonomic system: Fig. 7 provides an overview of the entire CAS. The cortical part of the autonomic system encompasses all cerebral lobes, whereas the subcortical autonomic system creates the link between prefrontal and insular cortex and the brainstem. The insular cortex turned out to be a central hub in the system, as it served as an interface for all of the sub-systems defined in Figs. 1-6 .

Integrating the cortical autonomic system into current anatomic and functional concepts
The insula turned out to be the central hub within the central autonomic system, especially as it forms the interface between the cortical and subcortical part of the network. Thus, it makes sense to focus on insular connectivity to integrate the newly tracked pathways into existing knowledge of networks. With respect to cortical regions, insular connectivity was extensively studied using a DTI based probabilistic fiber-tracking approach by ( Cloutman et al., 2012 ). There, three higher-order pathways were defined that share similarities to the tracts of the CAS in our study: 1. An "anterior pathway ", that originates in the anterior insula and is connected to pre-and orbitofrontal, mid-and inferior-temporal and middle occipital cortex shares similarities to the occipito-posterior-temporo-frontal pathway (OPTFP) and the rostralanterior part of the limbic pathway. This is also consistent with the insertion zone of these tracts in the ventral anterior insula and insular apex ( Fig. 8 ). With respect to the established nomenclature of association tracts, fibers of the OPTFP, which connect lingual gyrus and posterior middle and inferior temporal gyrus with the insula and prefrontal cortex, encompass parts of the inferior occipito-frontal fasciculus (IFOF; orbitofrontal to occipital; Kier et al., 2004 ;Lawes et al., 2008 ;Fernández-Miranda et al., 2008 ) and of the extreme capsule fiber system (FTecF, Petrides 2014 ), 2. A "posterior pathway " that is connected to the temporal pole and superior temporal gyrus, the supramarginal gyrus and middle occipital cortex resembles part of our parieto-anteriortemporal pathway (PATP). The insertion zones in the posterior rim of the insula and the insular apex ( Fig. 8 ) also fit well into this concept. The PATP may encompass parietal branches of the IFOF ( Martino et al., 2010 ), fibers of the uncinate fascicle (UF) and parts of the arcuate fascicle ( Catani et al., 2012 ). Furthermore, the PATP is characterized by fibers following the anterior commissure ( Catani and Thiebaut de Schotten, 2008 ) and hypothalamic projections (vide infra). 3. A "transitional pathway " that originates within the mid-dorsal insula and connects the cortex of the operculum to the structures of the temporal lobe. Connections from the operculum (insula) to hippocampus, parahippocampal gyrus and superior temporal gyrus are also present within the midcingulate limbic pathway that has its insertion zone in the mid-basal insula ( Fig. 8 ). However, the mid-cingulate limbic pathway is also characterized by additional prominent cingular fibers and fibers of the anterior commissure ( Catani and Thiebaut de Schotten 2008 ).
Apart from tractography, brain networks can be unrevealed by their functional coupling using fMRI for detection of coherent oscillations ( Deco et al., 2011 ). These resulting "resting-state " or "intrinsicconnectivity " networks cover several modalities, including a salience, default mode or an executive control network ( Seeley et al., 2007 ;Smitha et al., 2017 ). Especially the salience network also includes subcortical brainstem centres and shows a fair overlap with the ROIs used for streamline selection in our study (e.g. amygdala, anterior cingulate cortex, insula, hypothalamus and PAG; Seeley et al., 2007 ;Smitha et al., 2017 ). Apart from the detection of "salient " external stimuli, the salience network is thought to contribute to social behavior and self-awareness through the integration of multimodal information ( Menon, 2015 ). The involvement of similar structures for homeostatic regulation and emotional-motivational processing reflects the close coupling between physiological body states, feeling and emotion ( Critchley, 2005 ;Guo et al., 2016 ;Lutz et al., 2009 ;Pasquini et al., 2020b ). Within the salience network, the ventral anterior insula occupies a key position, as it becomes activated in social-emotional contexts ( Kurth et al., 2010 ) and is thought to mediate feeling states in response to environmental salient stimuli ( Phillips et al., 2003 ). Analyses of functional connectivity reveal a coupling to limbic but also frontal and temporal areas ( Uddin, 2015 ). In line with these findings, the ventral anterior insula is involved in a functional connectivity mode that involves the anterior cingulate cortex, temporal and medial orbitofrontal cortex ( Pasquini et al., 2020c ). This pattern fits well to the involvement of the ventral anterior insula in our OPTFP, PATP and the rostral part of the limbic pathway ( Fig. 8 ) and highlights the functional and anatomical overlap between the central autonomic system and the salience network.

Integrating the subcortical autonomic system into current anatomic and functional concepts
The essential autonomic reflex centers are located in the rostral medulla oblongata ( Benarroch, 1993 ;Dampney, 2016 ;Loewy, 1981 ;Silvani et al., 2016 ): a dorsal group of nuclei (i.e. nucleus of the vagus nerve, the nucleus ambiguus and the nucleus of the tractus solitarius) receives afferent input from baro-and chemoreceptors, nasopharyngeal receptors and skeletal muscles. Furthermore it contains parasympathetic preganglionic neurons that mediate the parasympathetic response. The rostral (ventro-) lateral medulla oblongata hence contains efferent sympathoexcitatory premotor neurons that trigger the systemic sympathetic response via activation of the intermediolateral nucleus of the thoracic spinal cord ( Mai and Paxinos, 2012 ). In our study, the dorsal and lateral medulla could be connected over the locus coeruleus and the periaqueductal grey towards the hypothalamus and further to the medial prefrontal cortex. Thus, this fibre system encompasses at least in part a composite of the dorsal longitudinal fasciculus (DLF; hypothalamus to PAG) and the medial longitudinal fasciculus (MLF; PAG to medulla; H. Zhang et al. 2012 ). Connections between mPFC, hypothalamus and PAG have been also described by retrograde tracing experiments in macaques Ongür et al., 1998 ). The PAG, in turn, is tightly coupled to rostral medullary nuclei. Within this system, parallel circuits are thought to process distinct types of stressful stimuli that become coupled to a particular autonomic response . Furthermore, connections between the locus coeruleus and the periaqueductal grey have been described in rats ( Ennis et al., 1991 ). The locus coeruleus exerts an excitatory effect on preganglionic sympathetic neurons within the rostral ventro-lateral medulla and is reciprocally coupled with the hypothalamus ( Samuels and Szabadi, 2008 ). Thus, this system routs informations about internal states (e.g., anxiety, arousal) or pleasant/stressful stimuli that are represented within cortical or subcortical regions to medullar autonomic reflex centres, where it becomes transduced to an autonomic response.

Connecting the cortical and subcortical autonomic network
With respect to anatomy and connectivity, detailed knowledge about the autonomic coupling between cortex and subcortical structures exists for non-human primates ( Benarroch, 2019 ;Evrard, 2019 ). In our study, the insula is connected with the brainstem by two tracts that target the anterior and posterior part of the dorsal insula ( Figs. 5 and 8 ). The dorsal insula bears a primary interoceptive representation receiving various sensory input that is routed from thalamic nuclei. With respect to topology, anterior parts process visceral information derived from the solitary tract, whereas the posterior part receives nociceptive information. Thus, the anterior and posterior connecting tracts may contain ascending viscero-and nociceptive afferents. However, these tracts directly reach the insula without being conducted over posterior thalamic nuclei. The lack of thalamic relay argues against an afferent nature of these tracts (Sherman and Guillery, 2006). Instead, they could be considered as a descending system enabling a "top-down " (i.e., subject-driven) control for brainstem signals as proposed for pain modulation ( Lu et al., 2016 ). However, this hypothesis requires experimental verification. Interestingly, the insular brainstem projections bypass the periaqueductal grey ( Fig. 5 ). Although there is evidence of fibers directly connecting the PAG with the posterior insula in rats , this connection may be too sparse to yield a robust tract using our tracking algorithm.
In non-human primates, efferent connections to the hypothalamus and brainstem nuclei arise from the ventral-anterior insula ( Benarroch, 2019 ;Evrard, 2019 ). As a target region of the parietoanterior-temporal pathway (PATP; Fig. 1 ), the hypothalamus is directly embedded within the cortico-insular network. In line with this view, hypothalamic connections between hypothalamus and amygdala (in cats; Hopkins and Holstege 1978 ), insular cortex (in cats; Clascá et al., 1989 ) and temporal regions (in monkeys; Adey and Meyer 1952 ) have been described. As mentioned above, the hypothalamus is connected to the periaqueductal grey , locus coeruleus ( Samuels and Szabadi, 2008 ) and the medulla oblongata (J. J. Zhang et al., 2012 ). The insular insertion zone of hypothalamic fibers within the ventral-anterior insula ( Fig. 8 ) fits also well into the concept of a "visceromotor " agranular insula ( Evrard, 2019 ). Thus, apart from the insula, the hypothalamus is an important interface region between CAN and SAN: it forms a relay of insular visceromotor efferents that furthermore operates under the control of prefrontal areas .

The role of the temporal lobe and laterality
The reconstruction of the central autonomic network is conditionally defined by initial selection of the ROIs. However, our global tractography algorithm rebuilds the entire connectome as a whole and not only reconstructs connecting streamlines between ROIs . Thus, streamline-bundles also serve to identify regions that do belong to the network, despite not being included in the initial sample of ROIs. Apart from the hypothalamus that is discussed above, the temporal lobe turned out to be an integral part of the CAN ( Fig. 6 ). Alterations of autonomic function (e.g., changes in blood pressure and heart rate) are well known symptoms in temporal lobe epilepsy ( Wannamaker, 1985 ) that also affect the inter-ictal phase ( Druschky, 2001 ). Human fMRI studies indicate an activation of the temporal cortex during Valsalva maneuver (for review see Sklerov et al., 2019 ) and an activation of the hippocampus, medial and superior temporal gyrus has been linked to parasympathetic regulation ( Beissner et al., 2013 ). Electrical stimulation of the vagal nerve influences excitability of pyramidal neurons within the hippocampal CA1 and CA3 subfield ( Radna and MacLean, 1981 ) and cardiopulmonary chemoreceptor activation modulates the expression of hippocampal theta rhythm ( Yu and Blessing, 2000 ). Thus, the inclusion of temporal regions into the central autonomic system is highly plausible.
Regarding the fronto-occipital and the mid-cingulate limbic pathways, the left sided tracts show a stronger expression when compared to the right sided ones. Laterality within the central autonomic system has been repeatedly demonstrated in functional neuroimaging ( Harper et al., 2000 ;Macefield and Henderson, 2016 ), even if a nonlateral stimulus (e.g., Valsalva maneuver) was applied. However, the side-preference differed between brain regions and also between subjects ( Harper et al., 2000 ). With respect to the anterior insular cortex, lateralization of activity is thought to balance the outflow of the central autonomic system: whereas the dominant (left-sided) hemisphere controls the parasympathetic tone, the non-dominant (right-sided) hemisphere is more related to sympathetic responses (Craig, 2009;Critchley and Harrison, 2013;Sturm et al., 2018 ). However, laterality in activation patterns is not necessarily the consequence of a stronger coupling on the level of tract anatomy. Furthermore, laterality in tract strength has also been observed within the "anterior pathway " that was also generated by a DWI based probabilistic fiber-tracking (see above, Cloutman et al., 2012 ). In contrast to our fronto-posterior-temporooccipital (OPTFP) tract, this "anterior pathway " that connects the insula with frontal and occipital cortical areas was more pronounced in the right hemisphere. Such asymmetries may arise from artefacts in the diffusion data, whose axis of susceptibility induced distortions is sagitally oriented.

Pathophysiological considerations and outlook
The regulation of autonomic processes can be disturbed in case of structural or functional perturbation of the central autonomic system. Thus "neuro-cardiogenic " sequela may occur after ischemic stroke, intracerebral hemorrhage or subarachnoid hemorrhage . These sequela comprise an increment in cardiac enzymes, ECG changes, arrhythmias and cardiac wall motion disorders that are caused by a sympathetic over-activation ( Fukuda et al., 2015 ;Samuels, 2007 ).
So far, research focused on grey substance regions to identify loci that are especially vulnerable regarding neurocardiogenic complications. For example, right - ( Ay et al., 2006 ) or left-sided ( LaowattanaZeger et al. ) ischemic lesions within the insula have been associated with cardiac adverse events or markers of myocardial injury. However, using a voxelbased lesion-symptom mapping approach, the occurrence of post-stroke cardiac arrhythmias was associated with a pattern of widespread lesions scattered across the entire right hemisphere ( Seifert et al., 2015 ). Thus, neurocardiogenic complications may not only occur after lesions in single strategic regions (e.g. insula), but also by the disturbance of the network at multiple points. Even though no focal tissue lesion occurs in the early phase of subarachnoid hemorrhage, neurocardiogenic complications occur even more often compared to intracerebral hemorrhage or ischemic stroke ( van der Bilt et al., 2009 ). Here, transient intracranial hypertension after the aneurysmal rupture ( Wybraniec et al., 2014 ) may lead to hypothalamic hypoperfusion ( Foreman, 2016 ) thereby disturbing the autonomic network.
In addition to neurovascular damage, the integrity of central autonomic circuits can be threatened by neurodegeneration ( Cersosimo and Benarroch, 2013 ). This is especially well investigated in patients suffering from the behavioral variant of frontotemporal dementia (bvFTD), a disease associated with an asymmetric degeneration particularly involving fronto-insular and anterior-cingulate cortices ( Seeley et al., 2012 ). Consecutively, patients are characterized by progressive impairment of social-emotional functions and empathy due to a progressive salience network dysfunction ( Guo et al., 2016 ;Pasquini et al., 2020a ;Rankin et al., 2006 ;Seeley et al., 2012 ). Prevalent left-sided (dominant) hemispheric degeneration is associated with a reduced baseline cardiac vagal tone indicating dysfunctional parasympathetic signalling ( Guo et al., 2016 ), whereas right-sided (non-dominant) affection of amygdala and hypothalamus is accompanied with impaired sympathetic signalling ( Sturm et al., 2018 ). Fiber integrity can be monitored MRI-based in vivo by extraction of mesoscopic diffusivity parameters like the intra-axonal volume or the intra-axonal water fraction (H. H. Zhang et al., 2012 ). Intra-axonal volume/water is inversely correlated to fiber damage Margoni et al., 2019 ), and corresponds to efficacious neuronal function ( Genç et al., 2018 ). Using a Bayesian estimator, extraction of mesoscopic diffusivity parameters can be accomplished with scanning/processing protocols suitable for daily clinical routine ( Reisert et al., 2017 ). Thus, impaired fiber integrity could be selectively recorded in autonomic tracts of patients suffering from cerebrovascular or neurodegenerative disease and correlated to markers of neurocardiogenic damage or dysfunction. This approach would shift the focus away from the view of single regions towards an estimation of network integrity.

Limitations of our study
As a consequence of our method, the following considerations should also be applied to the interpretation of our results: 1. Streamlines do not denote the directionality of projections, thus the denomination of "afferent " or "efferent " pathways was chosen with respect to the context. 2. It is not possible to distinguish between whether a synaptic connection occurs within a ROI, or if a streamline simply passes through it. 3. Since a threshold was applied to the probability of streamline occurrence, tracts with high individual variability in tract anatomy, or a low absolute streamline number, were omitted from our analysis. Therefore, the absence of a particular tract cannot exclude the existence of functionally relevant connections between ROIs. Furthermore, it was our aim to visualize a scaffold that harbors the flow of information within the central autonomic network, rather than carrying out an extensive structural connectivity analysis, or an exhaustive review of brain regions related to autonomic regulation. 4. On the other hand, the occurrence of a streamline does not prove the existence of an underlying neuronal connection ( Maier-Hein et al., 2017 ), thus, more evidence from other, independent techniques are needed to underpin the present findings. Also the notion "fiber density " has to be handled carefully as it typically correlates with the true underlying axonal/neurite density ( Yeh et al., 2020 ), but it is not directly measuring it. 5. The selection of ROI-combinations followed an unconstrained, random and iterative approach. As a start ROIs were applied individually before all possible ROI-combinations were subsequently tested for fiber selection. Combinations of ROIs that were not connected were excluded from further investigations. On the one hand, combining multiple ROIs discards and restricts connections that are involved in different functional or anatomical contexts. On the other hand, combination of multiple ROIs warrants specificity for a particular system.

Declarations
Funding: This work was funded by the BrainLinks-BrainTools programme of the German research Foundation and the German Council of Science and Humanities. The funders had no role in study design, data collection, data analysis or interpretation. The article processing charge was funded by the Baden-Wuerttemberg Ministry of Science, Research and Art and the University of Freiburg in the funding programme Open Access Publishing.
Availability of data, material and code: data, material and code will be provided by the authors upon reasonable request.

Ethical statement
All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. For more information, we refer to the Human Connectome Project homepage ( https://www.humanconnectome.org ).

Declaration of Competing Interest
There are no conflicts of interest.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.neuroimage.2021.117852 .