Axial T1-weighted structural fast spoiled gradient-echo scans (TE = 1.5 ms, TR = 6.3 ms, TI = 400 ms, 15° flip angle, 230 × 230 × 156 isotropic 1 mm voxels) and high angular resolution diffusion imaging data (55 directions, b = 1,000 s/mm², 72 1.8-mm-thick interleaved slices, 0.8594 mm × 0.8594 mm planar resolution) were acquired on a 3T GE Signa EXCITE scanner from 73 fully consented young healthy volunteers under a previous study approved by Weill Cornell’s institutional review board; for details of study protocols see Ref. (25). The exclusion criteria were pregnancy, a history of neurological or psychiatric diagnosis, seizure, or drug or alcohol abuse. Demographic characteristics of these young healthy subjects are shown in Table 1. Tractograms were extracted from these 73 young healthy subjects to create the normative connectome for the study. The measure of connectivity used in this paper is the anatomical connection strength (ACS) as proposed in Ref. (26). The overall differences in regional connectivity were normalized by a scaling factor equal to the sum of the connections. Thus, the connection strength, c_i,j, between ROIs i and j, was defined as the ACS score of streamlines connecting the two regions i and j. We refer to this network as graph G = {V,E} whose nodes ν_I ∈ V represent the ith GM region, and edges e_i,j ∈ E represent white matter fiber pathways whose connection strength is c_i,j. Connections are assumed to be bidirectional since directionality is not deducible from diffusion tensor imaging (DTI) tractography data. Gray matter regions from the T1-weighted images were parcellated using FreeSurfer volumetric pipeline (27) and a Desikan-Killiany atlas (28) into 68 cortical regions, 34 from each hemisphere and 18 subcortical regions. Six subjects were eliminated from the 73-subject dataset due to FreeSurfer failure or insufficient MR contrast.

Data used in preparation of this article were obtained from the 4RTNI database (see text footnote 1). Neurological and physical examinations, cognitive testing, behavioral testing, and other clinical data from PSP subjects between the ages of 45 and 90 were collected over three longitudinal visits (baseline, 6 months, and 12 months). The conditions for exclusion were any significant neurological disease other than PSP, including PD, multi-infarct dementia, HD, brain tumor, multiple sclerosis, or history of significant head trauma followed by persistent neurological deficits or known structural brain abnormalities. Demographic characteristics of these subjects are shown in Table 1. For our study, we used 60 PSP subjects after elimination due to unknown age, gender, and failed FreeSurfer quality control. Age-matched healthy controls (NC) data used in this study were obtained from the public Alzheimer’s disease neuroimaging initiative (ADNI) database² consisting of multimodal imaging data on AD and healthy subjects. T1-MPRAGE, FLAIR, DTI, ASL-perfusion (ASL), and resting state functional MRI (rs-fMRI) MRI sessions were acquired at 3T. The parameters for the PSP scans were chosen to match those used by the ADNI, including those being developed for ADNI-2 for ASL, rs-fMRI, and DTI at the time. Two-sample t-test was used to generate a random sample of 150 age-matched controls from ADNI database. Demographic characteristics of these subjects are shown in Table 1. This randomly generated sample of 150 age-matched healthy controls was used to generate a group atrophy vector for the study.

68 cortical and 18 subcortical volumes from 3T T1-weighted baseline MRI images were extracted using FreeSurfer software for both the cohorts. Estimated Total Intracranial Volume generated by FreeSurfer was used as an estimate for intracranial Volume (ICV) as a normalization measure to correct for head size in this study. ICV has previously been used in several studies for normalization particularly for neurodegenerative disorders for better estimation of regional atrophy that is caused by pathology (31, 32). Image processing steps were visually inspected for white–gray matter boundary and skull-stripping errors to ensure they had been carried out correctly. Subjects were eliminated from the dataset due to FreeSurfer failure or insufficient MR contrast. A vector of regional atrophy was created by using a two tailed t-test between PSP and NC mean ICV corrected regional volumes such that t_PSP = {t_PSP(i)|i ∈ [1,N]} (N = 86). The t-statistic was converted to the natural range [0,1] using the logistic transform, following (17). These atrophy measures were then used to test the propagation modeling analyses.

We modeled taupathy progression as a diffusion process on graph G. From Ref. (16) the transmission of pathology to all brain regions via the whole brain disease vector x(t) = {x(ν,t), ν ∈ V} and: (1)dx(t)dt=−βHx(t), where β is a global diffusivity constant and H is the well-known graph Laplacian. H=I−D−1C, where D is a diagonal matrix whose diagonal entries contain the degree of each node, degree being defined as the sum of weighted connections emanating from the node. Note, in order to accommodate regions having widely different out-degrees, we have used above the degree-normalized version of the Laplacian matrix (17).

Two cerebellar regions were removed, leaving regional PSP atrophy statistics on 84 cerebral regions. The NDM is described by x(t) and Φ(t), two 84 × 100 tables that represent the pathology (unobserved) and atrophy (measured empirically) in all 84 ROIs over 100 points in time. Pearson correlation strength (R statistic) and associated p-values were calculated comparing the empirical atrophy vector, t_PSP, with the ND atrophy pattern at all 100 points in time.

Figure 1 illustrates “glass brain” renderings of the spatial distribution of PSP atrophy. The spheres are located at the centroid of each of 84 brain regions, their size is proportional to the t-statistic of PSP atrophy after logistic transform and color coded by lobe per: frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green. The relative order of regional atrophy in this cohort roughly mirrors the spatial pattern of atrophy. HTH, pallidum, entorhinal, inferiortemporal, and superior frontal are the most atrophied regions as seen by the largest spheres in Figure 1.

The properties of the directional connectome built using both human and mouse connectome are described in Figure 2. First, we show that the connectivity patterns (Figures 2A,B) as well as the statistics of both the directional and the original non-directional connectomes are very similar (Figures 2C,D), as they should be. The directionality, defined such that −1 represents a purely retrograde connection, and +1 represents a purely anterograde connection, is shown as a histogram in Figure 2E. Clearly, most connections are bidirectional (value 0) and a small minority is directional, with equal numbers in both directions. Mainly subcortical structures display strong directionality (bottom and right portions of the matrices displayed in Figures 2A,B) while most corticocortical connections are bidirectional.

Next each region was computationally “seeded” in turn and NDM was played out over time on the canonical healthy connectome. Figure 3A shows spread of maximum Pearson correlation strength Rmaxi corresponding to the best fit between empirical data and the NDM seeded at region i. The distribution of Rmaxi is being shown as a histogram. Since Rmaxi serves as a measure of the likelihood of each region being a seed, this information is spatially depicted in the “glass-brain” insets. Table 2 shows top 20 regions with maximum correlation strength for each region seeded in turn. Clearly, the HTH and other limbic structures serve as the best seed regions for PSP. Figure 3B shows the R–t curves between the model evolution (Eq. 2), seeded at each region in turn and empirical PSP atrophy. For seed regions that are plausible, this would yield an intermediate peak in R where the NDM best matches empirical data, then diffusing uniformly and resembling actual data less and less. The highest R was achieved by the HTH, a region known to suffer early atrophy and could most likely be a seed to PSP for the connectome-based network diffusion.

Next in a model-free analysis, we establish the role played by proximate anatomic features in governing the regional patterns of PSP atrophy. We considered HTH as anatomic predictors based on highest R achieved from repeated seeding analysis above. Each covariate is a 84-long vector, covering the entire brain. Since the group atrophy t-statistic is generally bilateral, we removed lateralization effects by averaging the left- and right-hemispheric values of these vectors, giving predictor vectors of size 42 × 1. Linear bivariate correlation analyses of regional PSP atrophy with these proximate predictors are shown in Figure 3C. Dots are color coded as per lobes (frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green). We see a non-significant negative correlation between PSP atrophy and connectivity from bilateral HTH (R = −0.11, p < 0.01).

Having demonstrated that bilateral HTH gives the best seeding, we next captured the spatiotemporal evolution of x^HTH(t) (Figures 4A,B). The maximum of R^HTH(t) occurs at t_max = 22, hence x_HTH(t_max) is shown in Figure 4A. The ND progression matches quite closely the stereotypical sequencing in PSP, where the disease in due course extends into the subcortical and temporal structures and finally progresses to increasingly involve the cerebral cortex. The maximum correspondence of NDM to empirical data, given by the peak of the R–t curve in Figure 3B occurs at t = 22 for PSP. Figure 4B shows the evolution of network diffusion process seeded at the HTH starting at early stage (t = 7) through mature stage (t = 22), the model increasingly resembling empirical atrophy of Figure 1. We have selected three equidistant years at t = 7, 15, and 22 years between early stage (t = 0) with no diffusion through mature stage (t = 22). Here, time is arbitrary, and the use of “years” is meant for illustrative purpose. Initial spread is followed by diffusion especially into the pallidum and caudoputamen, followed by other striatal and limbic structures, then to middle temporal-frontal structures, and finally to the cortex.

We evaluated the NDM against alternate network models to show its specificity to PSP atrophy and to the connectome on which it evolves. We evaluated this in two ways and recorded the best R achieved by each model. First, we randomly scrambled the healthy average connectivity matrix C 2,000 times, and ran the NDM on each scrambled network for PSP. Second, we randomly scrambled the group t-statistic of regional PSP atrophy vector and ran the NDM on original connectivity matrix C. The distribution of Pearson’s R over 2,000 scrambled matrices is shown in Figure 5A and clearly indicates that our correlation results are unlikely to be due to chance. There is a hard limit on the left of this plot at R ~ 0.42, which corresponds to the 0-diffusion time value of R^PSP(t) curve in Figure 3B. The second random scrambling experiment, where the atrophy vector was scrambled instead of the connectome, gave an R distribution shown in Figure 5B. This distribution was approximately Gaussian, with mean between 0.1 and 0.2, and SD around 0.1 and 0.2. Random model’s R was much lower than the maximum R of 0.42 achieved by the true model; statistically outside the 95% confidence interval, or p < 0.05. Hence, the reported HTH seeded NDM prediction cannot be explained by chance.

Having demonstrated that bilateral HTH gives the best seeding, we next wanted to show if directionality could predict tauopathic spread either from same or different set of potential seeds. We tested both retrograde and anterograde mode in our model at t = 41 and t = 40, respectively. With anterograde connectivity there was an increase in maximum R (R = 0.48 at t = 40) as seen in R–t curve (Figure 6A). Our results showed that retrograde mode showed significant improvement in maximum R (R = 0.53 at t = 41) to both bidirectional and anterograde mode (Figure 6C). Both with retrograde and anterograde mode, the highest R was achieved by the HTH and served as the best seed region. Figures 6B,D show glass-brains with maximum R achieved by each region seeded in turn and Table 2 shows top 20 regions with maximum correlation strength for each region seeded in turn with directional connectivity.

We next captured the spatiotemporal evolution of x^HTH(t) (Figures 7A–D) at t = 40 and t = 41, which corresponds to time taken by diffusion to reach the highest peak with anterograde and retrograde directional connectivity, respectively. The maximum of R^HTH(t) occurs at t_max = 40 and 41 for anterograde and retrograde mode, hence x_HTH(t_max) are shown in Figures 7A,C for each mode, respectively. The ND progression shows better prediction of tau progression in PSP compared to non-directional seeding as seen in Figure 4. Both with anterograde and retrograde mode, the spreading pattern extends into the subcortical and temporal structures and finally progresses to increasingly involve the cerebral cortex.

The maximum correspondence of NDM to empirical data, given by the peak of the R–t curve in Figure 6A occurs at R = 0.48, t = 40 and in Figure 6C occurs at R = 0.53, t = 41 for PSP. Figures 7B,D shows the evolution of network diffusion process seeded at the HTH starting at early stage (t = 13) through mature stage (t = 39) with anterograde mode and through t = 14 and t = 41 with retrograde mode. With both the modes, the model increasingly resembles empirical atrophy of Figure 1 much so than with non-directional connectivity. We have selected three equidistant years at t = 13, 26, and 39 years between early stage (t = 0) with no diffusion through mature stage (t = 40) for anterograde mode and t = 14, 27, and 41 years between t = 0 and t = 41. Here, time is arbitrary, and the use of “years” is meant for illustrative purpose. As seen in Figure 7B with anterograde mode, initial seeding at HTH is followed by diffusion especially into the pallidum and caudoputamen, followed by limbic structures such as hippocampus and amygdala, then thalamus, then to middle temporal–frontal structures, and finally to the cortex. With retrograde mode (Figure 7D), initial seeding at HTH is followed by entorhinal, followed by other limbic structures such as hippocampus and amygdala, and thalamus, then to striatal, middle temporal-frontal structures, and finally to the cortex.

We tested whether network spread from focal seed sites would recapitulate PSP regional atrophy patterns. Although trans-neuronal transmission appears the most likely candidate, the exact mode of transmission is unknown. We first simulated spread between regions in a bidirectional manner, such that pathology has an equal chance of spreading in the anterograde or retrograde direction. To minimize model bias, we performed repeated seeding such that each brain region in turn has a chance to serve as the single onset region from which subsequent pathology transmission was modeled using the NDM. The HTH was found to be the best seed region, with the highest model correspondence of R = 0.42. The model performance seeded at HTH is significant in comparison with random network scrambling of the structural connectome (Figure 5) (p < 0.05), confirming that network organization is integral to disease spread in PSP, rather than an effect that has arisen by chance. Surprisingly, it was a better seed than any region in the striatum, even though striatal atrophy is usually the most prominent feature of the topography of PSP. The next best seed region was Pallidum, which is certainly plausible as converging post mortem and neuroimaging studies showing the striatum is the most affected region of pathology in PSP (8, 34).

Having identified structural connectivity as the most likely mechanism of transynaptic pathology spread in PSP, we sought to further improve our model by adding directionality. Clearly, a direct measurement of directionality of fiber tracts is impossible from current dMRI techniques, as water diffusion along fiber bundles does not respect cell polarity (soma to axonal terminal or vice versa). Instead, we exploit the well-known fact that homologous structures exist between many species, for some of whom we do happen to have anatomic connectivity data from painstaking tracer studies. A detailed mesoscale mouse connectome has recently become available from the Allen Institute (23). These data are fully directional, since it is based on retrograde AAV tracer experiments on a large number of mice. We, therefore, developed a novel technique whereby human and mouse brain atlas parcellations are used to define homologous brain structures, and mouse directional connectivity is transferred to non-directional human connectome.

Trans-neuronal transmission can have a distinct directional bias, such that misfolded protein species might follow anterograde or retrograde transport or signaling pathways. This is especially true in subcortical and striatal connections, which are known to be highly directional in comparison to corticocortical connections. Although little work has been done in PSP, tau pathology in AD differentiates between efferent and afferent connections (24). Hence, we tested the hypothesis that directional (whether anterograde or retrograde) process of spread along structural connectivity network will further improve model performance in PSP. The directional NDM results are shown in Figures 6 and 7, and confirm that either mode is a better fit than non-directional (Table 2). Interestingly, some aspects of PSP topography are better recapitulated with anterograde mode of transmission, especially striatal areas; whereas other aspects are better recapitulated by retrograde transmission, especially temporal and limbic areas. Quite unexpectedly, the HTH is the most likely source of pathology origin in all three modes of transmission, a result that raises many questions that should be addressed in future investigations.

Several methodological considerations should be considered when interpreting our results. The first are limitations of the NDM. The NDM is a first-order, linear, parsimonious model of diffusive spread that assumes that the structural connectivity network remains unchanged over the duration of the longitudinal analysis. Moreover, individual subject genetic repeat length, medication history and age of symptom onset were not available. These variables could have implications in the individual group wise analysis, when identifying each subject’s seed region or determining individual rate of disease diffusion. Because this is the first study to empirically test multiple network models of pathology spread in PSP, it will benefit from independent replication. Future work elucidating striatal vulnerability as well as the effect of repeat length on disease spread is necessary.