An afferent white matter pathway from the pulvinar to the amygdala facilitates fear recognition

Our ability to rapidly detect threats is thought to be subserved by a subcortical pathway that quickly conveys visual information to the amygdala. This neural shortcut has been demonstrated in animals but has rarely been shown in the human brain. Importantly, it remains unclear whether such a pathway might influence neural activity and behavior. We conducted a multimodal neuroimaging study of 622 participants from the Human Connectome Project. We applied probabilistic tractography to diffusion-weighted images, reconstructing a subcortical pathway to the amygdala from the superior colliculus via the pulvinar. We then computationally modeled the flow of haemodynamic activity during a face-viewing task and found evidence for a functionally afferent pulvinar-amygdala pathway. Critically, individuals with greater fibre density in this pathway also had stronger dynamic coupling and enhanced fearful face recognition. Our findings provide converging evidence for the recruitment of an afferent subcortical pulvinar connection to the amygdala that facilitates fear recognition. Editorial note: This article has been through an editorial process in which the authors decide how to respond to the issues raised during peer review. The Reviewing Editor's assessment is that minor issues remain unresolved (see decision letter).


Introduction
Decades ago, rodent research uncovered a subcortical pathway to the amygdala that rapidly transmits auditory signals of threat even when the auditory cortex is destroyed (Ledoux, 1998). Since this discovery, researchers have sought an equivalent visual pathway that might explain how it is that people with a lesioned primary visual cortex can still respond to affective visual stimuli that they cannot consciously see (Tamietto et al., 2010). The superior colliculus, pulvinar, and amygdala have been identified as nodes of a human subcortical route to the amygdala that bypasses the cortex (Morris et al., 1999). These subcortical areas consistently coactivate in cortically blind patients (Pegna et al., 2005) -as well as in healthy adults (Vuilleumier et al., 2003;Morris et al., 1999)when they view emotional stimuli, such as angry or fearful faces. Magnetoencephalography studies using computational modelling have investigated whether the activation of these subcortical nodes is causally related. These studies have consistently found evidence for a forward connection between the pulvinar and amygdala (McFadyen et al., 2017;Garvert et al., 2014;Rudrauf et al., 2008). The dynamic causal relationship between the superior colliculus and the pulvinar, however, remains unexplored in the human brain . The pulvinar also has several functional and cytoarchitectural subregions (Barron et al., 2015) and it is unclear how these connect to the superior colliculus and the amygdala and what roles these subregions may play in mediating transmission along the subcortical route (Koller et al., 2018;Pessoa and Adolphs, 2010). As such, the hypothesis that the subcortical route rapidly transfers information from the retina to the amygdala without interference has been heavily criticised (Pessoa and Adolphs, 2010;Pessoa and Adolphs, 2011). Furthermore, the pulvinar is highly connected with a widespread network of cortical regions that may contribute to transmission along the subcortical route (Bridge et al., 2016;Zhou et al., 2016). Hence, it remains unknown whether the functional activation of the human superior colliculus, pulvinar and amygdala during affective processing bears any relation to an underlying structural pathway (Pessoa and Adolphs, 2010).
Recent animal research has revealed several potential direct subcortical pathways that have a causal relationship with fearful behaviour in response to visual threats (Zhou et al., 2017;Wei et al., 2015;Shang et al., 2015). In the absence of relevant postmortem human research, however, our anatomical knowledge of the human subcortical route to the amygdala can only be derived from tractography of diffusion-weighted images (DWI). Tamietto et al., 2012 examined DWIs from a blindsight patient whose left primary visual cortex was destroyed. The white matter structure of the subcortical route was estimated for the patient and for ten healthy, age-matched controls. Critically, the subcortical route had greater fractional anisotropy in the patient's damaged hemisphere, suggesting a neuroplastic increase in structural connectivity to compensate for the disrupted cortical pathways (Tamietto et al., 2012). In a similar study, Rafal et al. (2015) used tractography to investigate the subcortical route in 20 healthy humans and eight macaques. The subcortical route was reconstructed in both hemispheres for 19 of the 20 human participants and 7 of the eight macaques (Rafal et al., 2015). Notably, this sample of human participants was recently expanded and re-examined, further demonstrating that individuals with greater fractional anisotropy along the subcortical route also had a stronger bias toward threat when making saccades to scenes (Koller et al., 2018).
Diffusion tractography may grant insight into the strength of anatomical connectivity between regions, but it cannot reveal the direction of information transfer nor can it be used as direct eLife digest : Being able to quickly detect and respond to potential threats is essential for survival. Fear and threat trigger a range of responses in the body, which are controlled by different regions in the brain. For example, a structure located deep within the brain called the amygdala is connected to other parts of the brain that regulate hormones, senses and muscles. The amygdala is highly responsive to signs of threat, and research in rodents has shown that it plays a role in transmitting sounds that indicate danger. However, so far it has remained unclear if this was also the case for visual information.
This is particularly challenging to study in humans because it has been difficult to image the deeper regions in the human brain. Now, McFadyen et al. reconstructed the pathways between the deeper brain regions important for processing vision and the amygdala using the brain scans of 622 participants. Then, they tested whether there was any connection between these pathways and the ability to recognise emotional expressions. To do so, fMRI brain scanning was used to measure the blood flow in the brain of volunteers looking at 40 faces that were either happy, sad, angry, fearful or neutral.
The results showed that when people were looking at pictures of fearful and angry faces, the blood flow between visual areas and the amygdala increased, especially in individuals with stronger connections, such denser nerve fibres, between the involved regions. The denser those fibres were, the better the people were at recognising when a face was fearful.
These discoveries suggest that the amygdala also plays a role in transmitting signals from deepbrain visual areas indicating danger and is likely to be one of the first areas to trigger a fear response in the brain. People with autism respond less to fearful faces, while people with anxiety respond more. Future research could investigate if the pathways to the amygdala differ in these people. evidence alone for the anatomical existence of a neural pathway. The anatomical presence and the direction-specific neural flow of emotional visual information along the subcortical route has never been concurrently investigated in humans to definitively show that the subcortical route is a direct, afferent pathway specifically associated with fear (Pessoa and Adolphs, 2010;Pessoa and Adolphs, 2011). Such a finding would have important implications for the very foundation of visual threat perception, given this pathway's potential for rapid information transfer (McFadyen et al., 2017;Silverstein and Ingvar, 2015) and unconscious processing (Tamietto et al., 2010). Here, we aimed to comprehensively investigate this putative amygdala pathway in a large sample of over 600 healthy human adults from the Human Connectome Project (HCP) dataset using a multimodal imaging approach to encompass structure, function and behaviour. First, we used DWI to reconstruct the subcortical route from the superior colliculus to the amygdala, via the pulvinar, and estimated its fibre density in a large sample. Next, we modelled the direction-specific flow of haemodynamic responses to faces, testing whether a functional subcortical route is recruited to transmit information toward the amygdala. Finally, we asked whether the fibre density of the subcortical route predicts both fearful face recognition as well as the strength of dynamic coupling between the superior colliculus, pulvinar, and amygdala.

Results
Reconstructing the subcortical route using tractography The first step in our investigation was to evaluate the evidence for an anatomical subcortical route to the amygdala in the healthy human brain. We exploited high-quality neuroimaging data from a large sample of 622 participants made available by the HCP . We then reconstructed the white matter structure of the subcortical route using two complementary tractography methods for cross-validation. We began with global tractography, a Bayesian approach to reconstructing whole-brain fibre configurations that best explain DWI data (see Materials and Methods for details). We discovered that the superior colliculus (SC) was connected to the pulvinar (PUL; fibre counts for PUL; left: M = 13.23, SD = 5.56, right: M = 13.00, SD = 5.59, minimum of 2 fibres per participant). The pulvinar and the amygdala were also connected (fibre counts for left: M = 5.33, SD = 2.79, and right: M = 6.75, SD = 2.90), with most participants having at least one connecting fibre (zero fibres for left PUL-AMG for eight participants -only 1.28% of total sample). Thus, this relatively conservative method of fibre reconstruction (as it takes into account the entire brain) can reliably detect evidence for a subcortical route across a large sample of participants.
After covarying out head motion and removing four participants with outlying standardised residuals (z-score threshold ±3), we established that there were significantly more fibres connecting the SC and PUL (M = 13.119, 95% CI = [12.738, 13.500]) than the PUL and AMG (M = 6.040, 95% CI = Figure 1. Global and local probabilistic tractography reconstructions of the subcortical route to the amygdala. Fibres reconstructed between the superior colliculus and the pulvinar are shown in (A) and between the pulvinar and amygdala in (B). Both (A) and the top row of (B) show 3D-renders of the major ROIs (amygdala in pink, pulvinar in orange, superior colliculus in green), as well as the left anterior thalamic radiation (orange) and the right corticospinal tract (pink). The reconstructed fibres for all participants are rendered by semi-transparent white streamlines. Streamlines of a single example participant are shown below in boxes. The bottom half of (B) shows a different, top-down perspective of the pulvinar to amygdala pathway, illustrating the right inferior frontal occipital fasciculus (light blue) and the left inferior longitudinal fasciculus (dark blue). (C) shows the global fibre counts (left graph) and average fibre density (right graph) for global and local tractography, respectively (SC-PUL in yellow, PUL-AMG in purple). Group average is indicated by solid white line, with dotted lines indicating the 95% confidence interval. (D) shows the average number of fibres terminating in each subregion of the pulvinar and amygdala, for both global and local tractography. 3D renders of the subregions (colour-coded to match the graphs) Figure 1 continued on next page [5.867, 6.214 To uncover more anatomical features of the reconstructed fibres, we used subregion-specific masks of the amygdala (basolateral, centromedial, and superficial) and the pulvinar (anterior, medial, superior, inferior, and lateral; see Materials and methods for ROI specification details) to determine where the reconstructed fibres terminated. This masking approach revealed that the global tractography fibres present between the SC and PUL connected predominantly to the inferior and anterior pulvinar (see Appendix 1-tables 3 to 5 for detailed statistics). Between the PUL and AMG, fibres terminated almost exclusively in the inferior PUL and then predominantly in the basolateral AMG. Hence, the inferior pulvinar served as the connecting node between the SC and the basolateral AMG for the globally reconstructed subcortical pathway.
To assess the validity of our findings we used a second tractography method, namely 'local' probabilistic streamline tractography, as used by Rafal et al. (2015) to reconstruct the subcortical route to the amygdala (Rafal et al., 2015). We generated streamlines between our regions of interest (ROIs) and found that the superior colliculus connected to the pulvinar (streamline counts for left: M = 1403.32, SD = 417.16, right: M = 111.59, SD = 358.60, minimum of six streamlines per participant) and the pulvinar connected to the amygdala (left: M = 575.42, SD = 203.03, right: M = 575.42, SD = 248.85, minimum 66 streamlines per participant). To evaluate whether these streamlines counts were reconstructed significantly above chance, we compared the numbers with those produced by a null distribution algorithm (Morris et al., 2008). We found that the number of streamlines was significantly different from chance for each connection, as determined by a series of paired two-sided ttests (see Appendix 1-table 9), suggesting that the DWI data produced meaningful streamlines between our ROIs.
We employed a recently developed method, SIFT2, which estimates the apparent fibre density of the streamlines connecting two regions of interest. This method more accurately represents the true underlying white matter structure (Smith et al., 2015). The apparent fibre density of the streamlines generated using local tractography followed the same pattern as the global tractography fibre counts, such that there was greater fibre density for the SC-PUL connection (M = 5.793,95% CI = [5.663,5.923 provide strong convergent evidence for a subcortical white matter pathway to the amygdala in the human brain. Like in the global tractography, we investigated the overlap between the locally generated tracks and known white matter fasciculi. The pattern of results was the same, with up to 60% of fibres traversing the anterior thalamic radiation, corticospinal tract, and inferior longitudinal and fronto-occipital fasciculi. For the SC-PUL pathway, the majority of track density was found in the anterior thalamic radiation (left: M = 52.85%, SD = 11.70%, range = 16.55% to 96.96%; right: M = 58.21%, SD = 16.24%, range = 5.49% to 89.30%) and the corticospinal tract (left: M = 12.89%, SD = 8.78%, range = 0.18% to 37.84%; right: M = 32.35%, SD = 12.49%, range = 0.56% to 63.64%; see Figure 1). For the PUL-AMG pathway, the majority was found in the corticospinal tract (left: M = 32.29%, Figure 1 continued are shown on the left. Graphs of the number of fibres terminating in pulvinar subregions are shown for tracts connecting the superior colliculus and pulvinar (first and second graphs) and the pulvinar and amygdala (third and fourth graphs), while fibres terminating in amygdala subregions are shown for tracts connecting the pulvinar and amygdala (fifth and sixth graphs). Error bars represent 95% confidence intervals. DOI: https://doi.org/10.7554/eLife.40766.003 SD = 15.58%, range = 0.59% to 65.25%; right: M = 37.38%, SD = 16.21%, range = 0.32% to 65.99%), followed by the anterior thalamic radiation (left: M = 16.47%, SD = 9.20%, range = 5.24% to 59.52%; right: M = 7.00%, SD = 3.49%, range = 1.46% to 50.09%), and then the inferior longitudinal (left: M = 5.69%, SD = 7.92%, range = 0% to 50.75%; right: M = 3.96%, SD = 6.33%, range = 0% to 47.51%) and fronto-occipital (left: M = 1.54%, SD = 4.21%, range = 0% to 72.79%; right: M = 7.96%, SD = 10.99%, range = 0% to 79.60%) fasciculi. Mean track densities were lower than 0.20% and 0.01% in other fasciculi for SC-PUL and PUL-AMG, respectively. We also examined which subregions of the pulvinar and amygdala the seeded probabilistic tracks terminated in. For the SC- The relationships between recognition of angry (purple), fearful (red), and sad (blue) faces (x-axes) and the predicted local fibre densities for the left (top row) and right (bottom row) PUL-AMG connection (y-axes), as produced by a multivariate ANOVA. Fearful face recognition accuracy was significantly related to local fibre density for the left and right PUL-AMG connections (b and p values shown for parameter estimates from multivariate ANOVA, *p < .05). DOI: https://doi.org/10.7554/eLife.40766.004 PUL pathway, the greatest number of streamlines terminated in the anterior PUL, followed by the inferior pulvinar (see Appendix 1-tables 6 to 8 for detailed statistics), consistent with the global tractography. Also like the global tractography, the local tractography fibres between the PUL and AMG terminated almost exclusively in the inferior PUL. For the AMG, however, fibres terminated predominantly in the basolateral subregion in the left hemisphere (consistent with the global tractography) but in the centromedial amygdala on the right.

Greater fibre density predicts better fearful face perception
We wanted to translate our work in humans to animal research that has demonstrated clear relationships between the anatomical presence of a subcortical route and fearful behaviour (Zhou et al., 2017;Wei et al., 2015;Shang et al., 2015). To this end, we examined behavioural data from an out-of-scanner task, the Penn Emotion Recognition task, that assessed a different component of face processing than the in-scanner task (analysed below). In the Penn Emotion Recognition task, participants were serially presented with 40 faces that were either happy, sad, angry, fearful, or neutral (8 faces presented in each category). Participants were most accurate with identifying the emotional expression of happy faces (M = 7.96, SD = 0.21), followed by neutral (M = 7.22, SD = 1.18), and then fearful faces (M = 7.02, SD = 1.03). Recognition was poorest for angry (M = 6.86, SD = 0.98) and sad faces (M = 6.82, SD = 1.12).
We then investigated the association between these scores (see Figure 2A) with the fibre density of the subcortical route. We chose not to include happy or neutral expressions in our analysis because the data were substantially negatively skewed (skewness for: happy = À5.821; neutral = À2.053; angry = À0.719; fearful = À1.188; sad = À1.090). Thus, we entered fibre density measures for the SC-PUL and PUL-AMG pathways into two separate multivariate regressions (one per tractography method, to reduce collinearity) with recognition accuracy scores for fearful, angry, and sad faces as covariates, plus head motion as a control covariate. We removed outliers (four for global tractography, 15 for local tractography) with z-scored residuals ± 3. While there were no significant multivariate relationships between global tractography and emotion recognition (see Appendix 1-tables 10 and 11 for detailed statistics), there was a significant relationship between local tractography and recognition accuracy for fearful faces (F(4,598) = 2.501, p=0.042, Wilk's L = 0.984, n p 2 = 0.016; see Figure 2B). This was driven predominantly by fibre density of the left (b = 0.140, p=0.004) and right (b = 0.143, p=0.012) PUL-AMG connections' local fibre density. The local fibre density of the left and right SC-PUL did not contribute significantly to the model. These results suggest that the fibre density of the PUL-AMG half of the subcortical route is associated with fearful face recognition more so than with other negative (sad) or threatening (angry) emotional expressions.

Subcortical and cortical BOLD signal to emotional faces
We used dynamic causal modelling to infer the dynamic (or effective) connectivity between each node of the subcortical route and determine the directionality of the functional interactions occurring along the anatomical pathway mapping described above. First, it was necessary to establish any differences in functional activation within these nodes. To do this, we used the 'Emotion' task from the HCP battery of fMRI experiments, in which participants performed a matching task using images of faces or shapes. We contrasted activation in face versus shape blocks and reviewed the results at the whole-brain group level across all 622 participants, p<0.05 FWE (see Figure 3). This revealed a network of significant BOLD clusters spread across occipital, temporal, frontal, parietal, and subcortical areas, replicating previous work with this dataset . Critically, the most significant cluster included the left and right amygdala as well as the left and right fusiform gyri (FG) and inferior occipital gyri (IOG). These latter two regions are key nodes in the cortical visual processing stream for faces, which may feed information forward to the amygdala (Tamietto et al., 2010). We used the SPM Anatomy Toolbox to confirm the anatomical positions of our functional activations. In the absence of an anatomical template for the superior colliculus and the pulvinar, we masked the map of statistically significant voxels (p < 0.05, FWE) with our a priori manual anatomically defined superior colliculi mask and functionally-defined pulvinar masks from Barron et al. (2015); see Materials and Methods for ROI generation). This revealed significant voxels in each area (proportion of significant voxels within each mask: left SC = 65.08%, right SC = 73.55%, left PUL = 36.13%, right PUL = 51.49%). Therefore, the faces-vs-shapes fMRI HCP task established functional activation in the three subcortical nodes of interest, as well as two nodes of a potential cortical pathway to the amygdala for conveying information about emotional faces.
A forwards-only subcortical route is engaged during face processing After observing significant BOLD signal in regions along the subcortical route as well as in other visual cortical areas, we next asked whether these regions were dynamically connected. We designed a space of testable models that mapped onto the functional hypothesis of a subcortical route to the amygdala that operates alongside a cortical visual pathway and is modulated by faces. Due to the presence of the IOG and FG in the whole-brain corrected fMRI activation and their known roles in face processing (Johnson, 2005), we defined several plausible functional cortical connections to the amygdala. These consisted of reciprocal pathways between IOG and FG, FG and amygdala, IOG and amygdala, as well as pulvinar and IOG (see Figure 4). Note that, while previous research has defined motion-related area V5/MT as a significant component of the pulvinar's subcortical visual network (Zhou et al., 2017), we did not observe strong involvement of this area in the faces-vs-shapes fMRI task (12% probability in cluster 37 with two voxels). Hence, we omitted area V5/MT from our model space. We named models containing both a cortical and subcortical route to the amygdala as 'Dual' models, whereas models in which the subcortical route was absent were named 'Cortical'. Our model space also included different sources of visual input, namely to the superior colliculus, pulvinar, or both, given that the pulvinar also receives direct retinal input (Cowey et al., 1994) as well as input via the superior colliculus (Berman and Wurtz, 2011). This gave us six families of models: 1) Cortical with SC input, 2) Cortical with PUL input, 3) Cortical with SC and PUL input, 4) Dual with SC input, 5) Dual with PUL input, and 6) Dual with SC and PUL input. We considered all possible combinations of forwards and reciprocal (forwards and backwards) cortical and subcortical connections, giving us a comprehensive model space of 102 models.
Of the available 622 participants, we conducted dynamic causal modelling on a subset of 237 participants who had sufficient above-threshold activation in all ROIs (these were defined by the subcortical masks used thus far and by spheres surrounding the coordinate of peak group BOLD signal in the IOG and FG; see Materials and Methods for more details). We conducted Bayesian model selection on the model space (grouped by families) to estimate how well the models explained the data, taking into account model complexity. We used the random effects implementation to account for potential individual differences in the recruitment of a subcortical pathway for viewing faces (Stephan et al., 2009). The winning family was the 'Dual with SC and PUL input' family (expected probability = 67.34%, exceedance probability = 100%) and the winning model across the entire model space was also within this family (expected probability = 21.24%, exceedance probability = 98.01%, protected exceedance probability = 98.18%; see Figure 4). This model included reciprocal cortical connections between IOG and FG and between FG and the amygdala. It also included a forwards-only subcortical route from the superior colliculus to the pulvinar to the amygdala, with input to both the superior colliculus and the pulvinar. The Bayesian Omnibus Risk score was p=1.78Â10 À124 , indicating a very small chance that the winning model was indistinguishable from all models tested (Rigoux et al., 2014). We replicated this finding (same winning family and winning model) on a subsample consisting of only the unrelated (i.e. non-sibling) participants within this group (49 participants; see Appendix 1).
The winning model revealed that the functional network that best explained the BOLD responses in our sample of 237 participants included visual inputs to the superior colliculus and pulvinar, forward connections from superior colliculus to the amygdala via the pulvinar, and recurrent interactions between IOG and FG, as well as between FG and amygdala. To extrapolate this finding to the general population and assess the consistency of dynamic coupling at each individual connection, we performed inferential statistics (t-tests) on the parameter estimates of each connection within the winning model (i.e. connectivity strength, in their natural space). We looked at the connectivity modulation parameters that represent the change in connection strength caused by the effect of faces. We removed extreme outliers (>3 SDs from mean) participants from each connection (M = 5.25, range = 3 to 8 participants excluded from sample of 237) and found that all connectivity modulations were significant (one sample t-tests against a test value of 1; see Appendix 1-table 13 for detailed statistics) except for the backward connection from left and right FG to left IOG (see Figure 5). These results suggest that the modulation of these connections by faces was consistently strong and so we can infer that a subcortical route for processing faces is likely present within the general population.
Greater fibre density relates to stronger effective connectivity Our findings from tractography, fMRI, and dynamic causal modelling provide convergent evidence for a subcortical route to the amygdala in humans. The final question we set out to answer was whether this converging evidence was correlated, such that participants with stronger structural connectivity also had stronger effective connectivity. In other words, we asked whether the structural connectivity along the subcortical amygdala route enables functional interactions amongst the nodes that lie within it. We computed eight partial correlations (with head motion as a control covariate) to examine the relationship between each parameter estimate and the corresponding global fibre count and local summed weights per connection (left and right SC-PUL and PUL-AMG). After removing multivariate outliers (leaving N = 213; see Appendix 1-table 12 for details), we discovered that participants with more global fibres also had greater modulatory activity for the right (r = 0.180, p = 0.004, p bonf = 0.032; see Figure 6) but not the left (r = 0.095, p = 0.084, p bonf = 0.672) PUL-AMG connection. The SC-PUL connection was not significantly related to its corresponding DCM parameters for global (left: r = À0.022, p = 0.627, p bonf = 1.000; right: r = 0.002, p = 0.488, p bonf = 1.000) or local (left: r = À0.101, p = 0.928, p bonf = 1.000; right: r = À0.028, p = 0.659, p bonf = 1.000) tractography. Note that we successfully replicated this finding within a subsample of unrelated participants (49 participants; see Appendix 1). Thus, our study is the first to successfully harmonise functional and structural information about the subcortical pulvinar connection to the amygdala.

Discussion
The elusive subcortical route to the amygdala has posed a unique challenge in studies of the human brain, due to its depth and its fast activation. Evidence has accumulated over recent years across many studies using various neuroimaging modalities showing that this pathway may underlie primitive threat-related behaviour. These studies, however, often take a unimodal approach on typically small samples, making it difficult to relate the specific structural connections between superior colliculus, pulvinar, and amygdala to observed functional brain activity and behavioural output. Our study, which used a large sample of participants from the HCP, supports the existence of a subcortical pulvinar connection to the amygdala in the healthy human adult brain that facilitates dynamic coupling between these regions and also enhances fear recognition. We reconstructed the subcortical route to the amygdala using sophisticated tractography methods and found that the white matter fibre density of the pulvinar-amygdala connection significantly predicted individuals' ability to recognise fearful faces. We then computationally modelled the functional neural networks along this structurally connected network that were engaged while people viewed emotional faces. We found that it was more likely for the network to include a subcortical visual route to the amygdala than a cortical route alone. Finally, we revealed converging evidence from structural and effectivity connectivity, such that the fibre density of the right pulvinar to amygdala pathway was positively correlated with the strength of the dynamic coupling (i.e. effective connectivity) between these regions.
This study marks the first time that structural and effective connectivity have been concurrently investigated in the one large sample to address the controversy on the existence and functional role of the putative subcortical route to the amygdala. Up to 60% of its fibre density overlapped with major fasciculi, including the corticospinal tract, anterior thalamic radiation, inferior longitudinal fasiculus, and inferior fronto-occipital fasciculus. Tractography of diffusion images is susceptible to both false positives and false negatives and thus is seldom used in isolation to determine the existence of particular neuroanatomical pathways (Jbabdi and Johansen-Berg, 2011). We established the validity of our tractographically reconstructed subcortical route by directly relating our measures of fibre density to both behaviour and effective connectivity, as well as by using two different tractography methods. Had the fibre density measures been simply due to noise, we would not have expected these theoretically relevant relationships with fearful face processing to emerge within this large sample of individuals. Notably, these intermodal relationships were only found for the pulvinaramygdala connection, despite there being greater fibre density between the superior colliculus and the pulvinar and this connection being present in the winning dynamic causal model. One explanation for this is that we had relatively less BOLD signal-to-noise ratio in the superior colliculus due to its small size and proximity to major blood vessels in the brain stem (Wall et al., 2009), thus weakening the likelihood of finding consistent covariance of its functional coupling with fibre density. Another explanation, particularly regarding the behaviour-tractography relationship, is that the pulvinar plays a significant functional role in the subcortical route to the amygdala. Research on macaques has demonstrated the pulvinar's response to emotional faces Maior et al., 2010) and its role in modulating attention  and so we would indeed expect the strength of the pulvinar-amygdala connection to be more predictive of fearful face recognition. Future research could more deeply investigate the relative contribution of each half of the subcortical route to emotional face processing by using an optimised fMRI approach (Wall et al., 2009) and contrasting different types of stimuli -for example, low vs. high spatial frequency (Gomes et al., 2017) or moving stimuli (Berman and Wurtz, 2011).
Our decision to reconstruct the two halves of the subcortical route separately was motivated by our interest in the relative contribution of each connection to face-related processing (as described above) but was also a limitation imposed by anatomically-constrained tractography, where reconstructed fibres are terminated at boundaries between grey and white matter (Smith et al., 2012). Given that the pulvinar is made up of thalamic cell bodies (grey matter), the likelihood of reconstructing a continuous streamline of axon bundles traversing the pulvinar's grey matter may have been restricted by these boundary constraints. Previous studies that have not imposed these constraints have successfully traced a continuous pathway from the superior colliculus to the amygdala via the pulvinar (Rafal et al., 2015;Tamietto et al., 2012), supporting animal research showing that inferior-lateral pulvinar neurons receiving superior colliculus afferents also have efferent connections to the lateral amygdala (Day-Brown et al., 2010). Our investigation into pulvinar and amygdala subregions support these findings, such that we found the superior colliculus to project predominantly onto the inferior (and anterior) pulvinar, which was the same subregion to receive the vast majority of fibres from the amygdala (see Figure 1). Furthermore, pulvinar fibres terminated predominantly within the basolateral amygdala, which is known to process visual information about threat and faces (Hortensius et al., 2016). Further studies could use both anatomically constrained tractography and this subregion-specific approach with ultra-high-resolution imaging to better differentiate grey-white matter boundaries and more accurately determine if and where a continuous, subcortical route might traverse the pulvinar.
While our results suggest that the inferior pulvinar may serve as a disynaptic connection point between the superior colliculus and amygdala, the continuity of information flow along the subcortical route is still a disputed feature due to the strong cortical influences on the pulvinar (Bridge et al., 2016;Pessoa and Adolphs, 2011). This dispute has also arisen from prior work investigating the spatial frequency content of information conveyed along the subcortical route. Research on blindsight patients has found evidence only for low spatial frequencies which suggests that such information originated from magnocellular cells in the superior colliculus (Burra et al., 2017;Méndez-Bértolo et al., 2016). On the other hand, work in healthy participants has found no such spatial frequency preference, which suggests that rapid pulvinar-amygdala transmission might include input from other parvocellular pathways (McFadyen et al., 2017). We did not exhaustively explore the extent to which the cortex contributes information to the pulvinar-amygdala connection. The winning effective connectivity model, however, did not include cortical connections between the pulvinar and the inferior occipital gyrus. Hence, it is unlikely that the primary visual cortex contributed (either via direct anatomical connections or functional coupling along the ventral visual stream; Pessoa and Adolphs, 2010) to the information transmitted along the subcortical route. The winning model did, however, include input to the superior colliculus as well as directly to the pulvinar, which could reflect direct retinal input or input from areas not explicitly included in the model, such as the parietal cortex, temporal cortex, or the LGN (Bridge et al., 2016), that may transmit both low and high spatial frequency information. Furthermore, it remains to be shown how interactions between the pulvinar and other cortical areas, such as the inferotemporal cortex (Zhou et al., 2016), may directly influence activity along the pulvinar-amygdala connection.
Our findings open avenues for future studies on how this subcortical pathway might influence threat-related behaviour. While our findings demonstrated that greater pulvinar-amygdala fibre density related to better fearful face recognition, it remains to be seen how this might compare with structural connectivity of other cortical networks. In other words, would the fibre density of this subcortical connection explain fearful face recognition above and beyond, say, structural connections between the inferior temporal or orbitofrontal cortex and the amygdala (Pessoa and Adolphs, 2011) or between the thalamus and the superior temporal sulcus (Leppänen and Nelson, 2009)? Evidence from blindsight patients suggests that this subcortical connection ensures redundancy and compensation, such that it strengthens when cortical connections are destroyed (Tamietto et al., 2012). Taking this in conjunction with our findings, we might consider that the pulvinar-amygdala connection contributes to fear recognition in faces (and effective connectivity underlying face perception) in healthy participants but can increase or decrease its influence depending on the functioning of other networks. Such increases and decreases are already evident in certain clinical populations. For example, structural connectivity between the superior colliculus, pulvinar, and amygdala is weakened in individuals with autism compared to healthy controls (Hu et al., 2017), and BOLD signal to fearful faces is reduced in these areas (Kleinhans et al., 2011;Green et al., 2017), unless participants are explicitly instructed to fixate on the eyes (Hadjikhani et al., 2017). On the other hand, people who suffer from anxiety show hyperactive activity along the subcortical route compared to non-anxious individuals (Hakamata et al., 2016;Tadayonnejad et al., 2016;Nakataki et al., 2017). How and why this subcortical visual pathway to the amygdala is altered in these clinical populations remains a significant and relatively unexplored avenue of research.
We observed hemispheric lateralisation of the pulvinar-amygdala connection, such that both the local and global tractography showed greater fibre density along the right than the left, and there were stronger tractography-behaviour and tractography-connectivity relationships for the right than the left. Early studies on the subcortical route observed specifically right-sided BOLD responses during non-conscious fearful face viewing (Morris et al., 1999;Morris et al., 1998), and a previous tractography study has also found that only the fractional anisotropy of the right subcortical route was significantly related to threat-biased saccades (Koller et al., 2018). There is mounting evidence for right-sided specialisation for ordered (Wyczesany et al., 2018) and disordered (McDonald, 2017) emotion processing, particularly for non-conscious signals transmitted along the subcortical route (Gainotti, 2012). Thus, our results lend support to this theory by demonstrating evidence for the right pulvinar-amygdala connection's stronger fibre density and its relationship to emotional face viewing and fearful face recognition. Our understanding of this lateralisation may be deepened by future exploration of left-vs. right-sided structural connectivity and function along the subcortical route during conscious vs. non-conscious emotion processing in healthy participants.
One limitation of the present study is the discrepancy between how local and global measures of fibre density related to other measures; namely, that local tractography covaried with fearful face recognition scores while global tractography covaried with effective connectivity. While the reconstructed fibres shared many similarities (e.g. the pattern of findings for each connection across hemispheres and subregions, as well as the overlap with major fasciculi; see Figure 1) even after accounting for head motion, it is possible that the local tractography's relatively greater susceptibility to noise may have decreased its relationship to corresponding effective connectivity parameters. Indeed, global tractography has been shown to better reflect local connection architecture (Jbabdi and Johansen-Berg, 2011), such as the subcortical connections we have investigated. Such discrepancies between global and local tractography have been reported in other work (Anastasopoulos et al., 2014) and so further research (particularly those that only recruit a single tractograpy method) will benefit from specific investigations into why these discrepancies might arise.
In conclusion, our study has made substantial progress towards settling the long-held debate over the existence and function of a subcortical route to the amygdala in the human brain. Our multimodal neuroimaging approach, leveraged by computational modelling, provides convergent evidence for a fundamental and conserved pulvinar-amygdala pathway that is specifically involved in fear. We demonstrate that the white matter tracts that form the subcortical structural pathway from the pulvinar to the amygdala enables functional, dynamic interactions involved in emotional face perception. Critically, we show that structural connectivity between the pulvinar and the amygdala leads to better recognition of fearful expressions.

Participants
We used the data from the publicly available Human Connectome Project (HCP) S900 release, collected between 2012 and 2016, containing data from 897 consenting adults . Ethical permission to use this data and the associated restricted access data (including variables such as specific age information) was obtained from the University of Queensland Human Research Ethics Committee. Out of these participants, 730 young adults had complete MRI and dMRI data, as well as fMRI data for the faces-vs-shapes task (Van Essen et al., 2012). Of these, we excluded 95 people due to positive drug/alcohol tests and an additional 13 for abnormal colour vision. This resulted in a final sample of 622 participants aged between 22 and 36 years (M = 28.81, SD = 3.68 years), 259 of whom were male and 363 female, with 569 right-handed and 53 lefthanded. Within our sample, 495 participants were related to one or more other participants (328 families in total). This included 53 pairs of monozygotic twins, 50 pairs of dizygotic twins, and 289 participants with one or more non-twin siblings in the sample. The remaining 127 participants were unrelated. We acknowledged that the many siblings in the HCP sample might spuriously decrease the variance in our neural measures (due to the structural and functional similarity between siblings, for example) and thus influence our statistics. Because of this, we replicated some of the analyses from the full sample on the subsample of unrelated participants (see Appendix 1).

dMRI preprocessing
We used the minimally processed images provided by the HCP. For the T1 images, this included gradient distortion correction, bias field correction (using FSL: Jenkinson, Beckmann, Behrens, Woolrich, and Smith, 2012; http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/), and cortical segmentation (using FreeSurfer: Dale, Fischl, and Sereno, 1999; http://surfer.nmr.mgh.harvard.edu/). For the diffusion images, this included intensity normalisation across runs, echo planar imaging (EPI) distortion correction and eddy current/motion correction (for full details, see Glasser et al., 2013). The latter stage produced motion parameters (three translations and three rotations), across which we computed the mean relative displacement between dMRI volumes. These values per participant were used in all further analyses to account for any confounding effect of motion (Baum et al., 2018). We conducted all further processing using MRTrix 3.0.15  and FSL 5.0.7.

Global intensity normalisation
First, we corrected for low-frequency B1 field inhomogeneities in the dMRI volumes. We then conducted global intensity normalisation across each participant's corrected diffusion-weighted image so that we could later perform quantitative analyses of fibre density (i.e. apparent fibre density; Raffelt et al., 2012). This step normalises the median white matter b = 0 intensity (i.e. non-diffusionweighted image) across participants so that the proportion of one tissue type within a voxel does not influence the diffusion-weighted signal in another. Given our large sample size, we selected a subset of 62 participants (approximately 10% of the sample) to create a representative fractional anisotropy (FA) population template and white matter mask. We then used the population template and white matter mask to normalise the white matter intensity of all 622 participants' dMRI volumes.

Response function estimation
We segmented each participant's T1 image into five tissue types (cortical grey matter, subcortical grey matter, white matter, CSF, and pathological tissue) using the Freesurfer parcellation image provided by the HCP. We then estimated response functions (i.e. the signal expected for a voxel containing a single, coherently-oriented fibre bundle) for grey matter, white matter, and CSF using the Multi-Shell Multi-Tissue (MSMT) Constrained Spherical Deconvolution (CSD) algorithm (Jeurissen et al., 2014). After completing this step for all participants, we averaged their response functions to produce representative response functions per tissue type. We then conducted MSMT CSD on each participant again using the group averaged response functions, producing individual multi-tissue fibre orientation distributions (FODs).

fMRI processing fMRI acquisition
As with the dMRI data, the HCP acquired whole-brain gradient-echo echo planar imaging (EPI) data using the Connectome Skyra (TR = 720 ms, TE = 33.1 ms, flip angle = 52˚, bandwidth = 2,290 Hz/ Px, in-plane field of view = 208Â180 mm, 72 slices at 2 mm thick, voxel size = 2 mm isotropic, echo spacing = 0.58 ms) with a multiband factor of eight. They collected data in a one-hour session (either on the same day as the dMRI or one day before/after) along with two or three other functional tasks in the HCP battery. For the faces-vs-shapes task, there were two runs, one with right-to-left phase encoding and the other with left-to-right phase encoding, each with 176 frames at a duration of 2 min and 16 s.

fMRI task
The HCP developed the 'emotion task' (i.e. faces vs. shapes) from the paradigm presented by Hariri et al. (2002). Participants were presented with three visual stimuli at a time (one image at the top and two at the bottom) using E-Prime (Schneider et al., 2002). Participants were then instructed to make a button press indicating which of the two images at the bottom (left or right) matched the image at the top. Images were either a face (angry or fearful) or a shape (circle, horizontal oval, or vertical oval). The stimuli were presented on screen for 2000 ms separated by a 1000 ms inter-stimulus interval. A cue was presented at the beginning of each block to indicate the block type (i.e. 'shape' or 'face'), where each block contained either six faces trials or six shapes trials. Finally, a fixation cross was presented for eight seconds at the end of each of each run. The last block of each run only contained the first three trials due to a technical error that occurred early in HCP data collection. As the first block was always a shape block, our analysis was conducted on three shape blocks and 2.5 face blocks.

fMRI preprocessing
We used the minimally preprocessed fMRI data provided by the HCP corrected for gradient distortion, motion, and field map-based EPI distortion. The HCP intensity normalised the data and spatially transformed it to MNI152 space using FSL (see Glasser et al., 2013) for full details on preprocessing pipeline). We further increased the signal-to-noise ratio of the fMRI data in SPM12 (SPM12, www.fil.ion.ucl.ac.uk/spm) by applying spatial smoothing using a 4 mm Gaussian kernel (Hillebrandt et al., 2014).

Regions of interest
We chose the superior colliculus, pulvinar, and amygdala as our ROIs. We created masks of these ROIs in standard MNI space using FSL. For the amygdala (AMG) binary mask, we used the probabilistic Harvard-Oxford Subcortical atlas at a threshold of at least 50% probability. For amygdala subregions, we used the basolateral, centromedial, and superficial amygdala regions in the Juelich Histological Atlas (Amunts et al., 2005) at a threshold of at least 40% probability. For the pulvinar (PUL), we were interested in the structure as a whole, as well as its subregions (results for the latter are detailed in Appendix 1). To do this, we used the parcellated pulvinar mask generated by Barron et al. (2015), who isolated five distinct pulvinar clusters based on functional co-activation profiles in fMRI data from 29,597 participants across 7772 experiments (Barron et al., 2015). For the pulvinar as a whole ROI, we merged the five clusters together and used FSL to manually fill any holes in the resultant binary mask. Finally, we manually created binary masks for the left and right superior colliculi (SC) in the absence of an atlas-based mask by drawing the boundaries of the superior colliculus over the MNI152 single participant T1 template with reference to an anatomical atlas (Tamraz and Comair, 2004) and filling the centre. We then used FSL to warp these masks into native diffusion space for each participant's tractography analysis. All our ROIs in MNI space are freely available online from the Open Science Framework: doi:10.17605/OSF.IO/KBPWM.

dMRI analysis
In this study, we implemented two tractography methods that use different approaches to white matter reconstruction for cross-method validation. We first used the multi-tissue model of global tractography. This method takes a Bayesian approach to reconstructing a full-brain fibre configuration using a generative signal model to best explain the underlying data. It is less sensitive to noise that may accumulate for longer distance tracts in other 'local' tractography methods throughout their stepwise approach (Christiaens et al., 2015;Reisert et al., 2011). Hence, for comparison, we computed probabilistic ('local') tractography between our regions of interest (Tournier et al., 2010). This method also uses a Bayesian approach to account for one or more distributions of fibre orientations within each voxel, thus incorporating uncertainty into the model (Zhou et al., 2017). To acquire a biologically accurate measure of apparent fibre density (Raffelt et al., 2012) along the resultant streamlines, we used the Spherical-Deconvolution Informed Filtering of Tractograms version 2 (SIFT2) method to weight each streamline by a cross-sectional area multiplier directly related to the underlying data (Smith et al., 2015;Raffelt et al., 2012). For both the global (producing 'fibre count' as a variable) and local tractography with SIFT2 (producing 'summed weights' as a variable), we computed 2 (hemisphere: left, right) by 2 (connection: SC-PUL, PUL-AMG) repeated-measures ANOVAs to quantitatively examine the properties of these pathways.

Global tractography
Global tractography is a data-driven Bayesian approach to estimating the whole-brain fibre configuration that best explains the underlying diffusion-weighted images. As opposed to local streamline tracking, global tractography accounts for the spatial continuity of fibres and thus is better able to discriminate crossing and fanning fibre geometries (Christiaens et al., 2015). Furthermore, because the simultaneously-reconstructed fibre configurations are optimised with respect to the data at hand, the density of the final tractogram quantitatively represents the apparent fibre density (AFD; i.e. the proportion of space occupied by white matter fibres (Raffelt et al., 2012).
We conducted global tractography on the global-intensity-normalised DWI volumes for each participant using the group-averaged multi-tissue response functions. After 250 million iterations to optimise a full brain reconstruction, we filtered the tractogram using the ROI masks described above to isolate fibres that terminated in 1) both the superior colliculus and pulvinar masks, and 2) both the pulvinar and amygdala masks. We also used the masks for the five individual functionally-defined pulvinar subregions to isolate the subregion-specific fibres connecting to the superior colliculus and to the amygdala.

Local tractography
As a less conservative approach than global tractography, we also conducted local probabilistic tractography using our ROIs as seeding and terminating regions. We used the iFOD2 algorithm, iteratively planting a seed point 25,000 times (or until at least 10,000 tracks had been selected and written) in each voxel of the seeding ROI . We applied the anatomically-constrained variation of this technique, whereby each participant's five-tissue-type segmented T1 image provided biologically realistic priors for streamline generation, reducing the likelihood of false positives (Smith et al., 2012). We edited the final streamlines so that only those that terminated at white-grey matter boundaries in our ROIs remained.
Using these methods, we traced streamlines between the two halves of the subcortical route (i.e. SC to PUL, PUL to AMG). We then reversed the seeding location (i.e. PUL to SC, AMG to PUL) and based all statistics on the average between the forwards and backwards seeding directions to reduce any influence of possible asymmetries in seed ROI volume. We applied SIFT2 to these streamlines to enable us to quantitatively assess the connectivity. SIFT2 makes this possible by weighting the streamlines by a cross-sectional multiplier such that the sum of these weighting factors better represents the underlying white matter fibre density (Smith et al., 2012).

fMRI analysis General linear modelling
Using the spatially smoothed fMRI data, we convolved the onset of each Face and Shape block with a canonical hemodynamic response function (HRF) using SPM12. We closely modelled this first-level general linear model (GLM) analysis on the work by Hillebrandt et al. (2014), such that we did not slice time correct the multiband data due to the fast TR. We partitioned the GLM into sessions (leftto-right and right-to-left encoding) and we included 12 head motion parameters as multiple regressors (six estimates from rigid-body transformation, and their temporal derivatives). We generated statistical parametric maps (SPMs) of the expected BOLD signal for faces minus shapes and shapes minus faces.
We then entered the faces minus shapes contrast into a second-level analysis (a one-sample ttest) across all participants. After examining the estimated BOLD signal to Faces at the whole-brain level (p<0.05, family-wise error corrected), we applied the superior colliculus, pulvinar, and amygdala a priori defined masks to more specifically estimate functional activation in these anatomicallydefined areas.

Dynamic causal modelling
We implemented Dynamic Causal Modelling (DCM) to infer the causal direction of information flow between neural regions using a biophysically informed generative model (Friston et al., 2003). First, we examined the map of significant activation produced by the fMRI analysis of the 'Emotion' (i.e. faces vs. shapes) HCP task. Based on this and our a priori hypotheses, we defined the left and right superior colliculus, pulvinar, amygdala, inferior occipital gyrus (IOG), and fusiform gyrus (FG) as our ROIs. For the two gyri, we used MNI coordinates of the most significant peak from the group level analysis (left IOG: À22-92 À10, right IOG: 28-90 À8, left FG: À38-50 À20, right FG: 40-52 À18). We then placed spheres with a radius of 12 mm around these four coordinates to search for the participant-specific local maxima within each participant's session-specific SPM for the faces minus shapes contrast (adjusted for the t effects of interest, p < 0.05 uncorrected). Note that for the purposes of extracting the fMRI data for the DCM nodes, one does not need corrected p-values (Hillebrandt et al., 2014). Next, we defined the ROIs by a 6 mm radius sphere around the participant-and session-specific local maxima. For the subcortical areas of interest, we defined the initial search radius by the anatomically defined ROI masks (as described above) instead of significant peaks from the group analysis to confine our search within subcortical grey matter.
We used a 'two-state' DCM model, which accounts for both excitatory and inhibitory neural populations (Hillebrandt et al., 2014;Marreiros et al., 2008). Our model space was dictated by our specific, theory-driven hypotheses about subcortical and cortical visual pathways to the amygdala, as well as by the significant regions of the BOLD signal observed at the group level in our GLM analysis. Both face and shape blocks contributed to input parameters within each model. All endogenous and intrinsic connections in each model were modulated by the effect of faces over shapes.
To specify a DCM, each participant needed to have above-threshold activation (at p < 0.05, uncorrected) within each ROI across both scanning sessions. This was the case for 237 out of the 622 participants (see Appendix 1-figure 1). The ROIs with the highest numbers of below-threshold participants were the left and right superior colliculi (261 and 246 participants, respectively), followed by the left and right pulvinar (46 and 32 participants, respectively), and finally the left and right amygdala (40 and 25 participants, respectively). This may be due to the bilateral superior colliculi's relatively smaller volume as well as lower statistical power (its mean t-statistic was approximately 10.57 compared with 33.35 for IOG and 30.14 for FG). Critically, the group of 237 participants with above-threshold BOLD responses in all ROIs did not differ significantly from the other group of 385 participants in the global or local tractography results (main effect of 'group' and interactions with 'group' were all p > 0. 162  .717) during the fMRI task, age (p = 0.782), or gender (p = 0.359). Therefore, using the information available to us, we had no evidence to assume that our DCM sample was biased by any confounding variable. The final model space consisted of 102 models (see Figure 5), where the first Cortical family contained six models, the second and third Cortical families contained 12 models each, and the Dual families (families 4, 5, and 6) contained 24 models each. The different families correspond to different input types (superior colliculus only, pulvinar only, or superior colliculus and pulvinar) and the different models within these families arise from different combinations of forward and backward connections. Each of the final 102 DCMs were modelled separately for both fMRI sessions. Both hemispheres were included in each model with no cross-hemispheric connections. To determine which model best explained the data, we conducted family-wise Bayesian Model Selection (Stephan et al., 2009;Rigoux et al., 2014), which penalises models for complexity according to the free energy principle (Friston et al., 2006). We used the random effects implementation to account for potential individual differences in the recruitment of a subcortical pathway for viewing faces (Stephan et al., 2009).

Code availability
All computer codes that were used to produce the results (from raw HCP data to track counts, fibre density, BOLD signal and DCM files) is freely available online via GitHub (McFadyen, 2018; copy archived at https://github.com/elifesciences-publications/hcp-diffusion-dcm) and the Open Science Framework (doi:10.17605/OSF.IO/KBPWM).

Data availability
The data analysed in this study came from the publicly-available Human Connectome Project S900 release: https://www.humanconnectome.org/study/hcp-young-adult/document/900-subjects-datarelease. Restricted access was obtained through the HCP to acquire specific participant ages (in years) and drug/alcohol information. Ethical permission was granted by the University of Queensland Human Research Ethics Committee. No figures display raw data. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability
The data analysed in this study came from the publicly-available Human Connectome Project S900 release: https://www.humanconnectome.org/study/hcp-young-adult/document/900-subjects-datarelease. Restricted access was obtained through the HCP to acquire specific participant ages (in years) and drug/alcohol information. Computer code used to produce the results is available on GitHub (https://github.com/jjmcfadyen/hcp-diffusion-dcm; copy archived at https://github.com/elifesciences-publications/hcp-diffusion-dcm).
The following datasets were generated: Appendix implemented in MRtrix 3. Outliers with a difference more than three standard deviations from the mean were removed from the 622-participant sample for each paired t-test.
Bootstrapping was set to 1000 iterations.

Multivariate diffusion-behaviour relationships
We conducted two separate multivariate analyses of covariance, one for global and one for local tractography measures of fibres, to examine the relationship between emotional expression recognition (fearful, sad, and angry expressions) and fibre density along the subcortical route.
Partial correlations between fibre density and effective connectivity parameters We conducted a series of right-sided Pearson's partial correlations between the fibre density of each connection (left and right SC-PUL and PUL-AMG connections, global and local tractography -giving eight in total) and its corresponding DCM parameter estimate (modulatory effect of Faces > Shapes over region coupling). Head motion was entered as a control variable. Multivariate outliers were detected according to Mahalanobis distance (df = 8, c 2 criterion = 15.507, p=0.05), resulting in 24 participants being excluded from the analysis (N = 213). Note that outliers were substantially influencing the results (see table below for comparison).

DCM parameter estimates
We conducted one-sample t-tests on the exponentiated DCM parameter estimates (B matrix) against a test value of 1 to examine whether modulatory connection strength was consistently greater than the prior within our sample of participants. Outliers were removed that were more than 3 SDs from the mean of each variable (note that this did not change the pattern of results).