Altered network stability in progressive supranuclear palsy

Highlights • We investigated network dynamics in the tauopathy progressive supranuclear palsy• Abnormal temporal properties of large-scale networks are related to phenotype• Progressive supranuclear palsy paradoxically increases frontoparietal state time• Reductions in neural signal complexity relate to altered network dynamics• Dynamic network and topological changes occur distally to primary sites of atrophy


Introduction
The human brain optimises efficiency by balancing integration and segregation of information transfer among neural assemblies. The activity and connectivity of regional specialisation is dynamic ( Deco et al., 2017 ;Friston et al., 2012 ;Honey et al., 2007 ;Shine et al., 2019 ;Tognoli and Kelso, 2014 ), even on the suprasec-work dynamics we use task-free functional MRI (fMRI) from two independent cohorts.
The clinical features of PSP, together with its established imaging and pathological findings, qualify it as a model disease to investigate network dynamics. PSP has prominent cognitive and behavioral features, including a dysexecutive frontal syndrome, apathy, impulsivity and language impairment, in addition to the movement disorder of axial akinetic-rigidity and impaired postural reflexes ( Burrell et al., 2014 ;Steele et al., 1964 ). The disruption to static functional connectivity in PSP ( Brown et al., 2017 ;Gardner et al., 2013 ;Rosskopf et al., 2017 ;Whitwell et al., 2011 ) affects frontal cortical regions associated with cognitive control and behavior, alongside striatal degeneration and loss of dopaminergic and noradrenergic projections from the brainstem to forebrain ( Murley and Rowe, 2018 ). The latter are critical to balance network integration and segregation ( Shine, 2019 ), with catecholaminergic deficits related to dynamic connectivity, cognitive performance and disease severity ( Eldar et al., 2013 ;Kaalund et al., 2020 ;Shine et al., 2018 ).
Network dynamics need to be interpreted in the context of neural complexity ( Honey et al., 2007 ;McDonough and Nashiro, 2014 ). Complexity varies with the timescale analyzed, with the potential for scale dependent relationships between integrative or synchronous activity, state switching and complexity ( McDonough and Nashiro, 2014 ;McIntosh et al., 2014 ;Wang et al., 2018 ). This relationship may be due to interference between neural complexity and regional phase relationships, decreasing the likelihood of synchrony between brain regions ( Ghanbari et al., 2015 ). Alternatively, sufficient signal complexity may be required to establish long range dependencies, leading to a positive relationship between connectivity and complexity conditional on timescale ( McDonough and Nashiro, 2014 ;Wang et al., 2018 ).
Entropy measures can assess complexity in the relatively short, non-linear and noisy time series typical of fMRI ( Grandy et al., 2016 ;Pincus and Goldberger, 1994 ;Turkheimer et al., 2015 ). Sample entropy measures the likelihood that repeated patterns are present in data: signals with a repetitive structure have lower entropy ( Richman and Moorman, 20 0 0 ). Multiscale entropy (MSE) extends sample entropy by assessment at multiple timescales, with the advantage that random noise can be differentiated from complex signal: random fluctuations increase entropy at fine time scales, but with increasing the timescale entropy decreases ( Costa et al., 2005 ).
Large-scale network dynamics can be quantified by hidden Markov modelling (HMM), in terms of a finite number of mutually exclusive states between which the brain switches over time ( Vidaurre et al., 2017 ). While states are inferred at a group level, information is obtained about the order and timing of individuals' state transitions.
We used these methods to investigate the impact of PSP on network dynamics, as a function of changes in signal complexity, brain structure and functional reorganisation. We examined two contemporary but independent cohorts of PSP, and controls, from the Cambridge Centre for Parkinson-plus (CCPP) and the UK national PSP Research Network (PROSPECT-MR). In each cohort, we analyzed HMM and MSE of task-free functional magnetic resonance. We then tested whether the network properties in PSP varied as a function of disease severity (PSP-rating-scale) and PSP phenotype (Richardson's syndrome, cortical-and sub-cortical variants).

Participants
Forty-five participants with PSP (possible or probable, according to MDS-PSP criteria ( Höglinger et al., 2017 )) and 27 controls were recruited at the Cambridge University Centre for Parkinson-plus (CCPP). 49 study participants with PSP and 37 controls were recruited to Progressive Supranuclear Palsy-Corticobasal Syndrome-Multiple System Atrophy-UK (PROSPECT-MR) study ( Jabbari et al., 2019 ). 34 participants (28 PSP, 6 controls) were removed following assessment of motion and image quality ( section 2.2 ). For both cohorts clinical assessment included the PSP-ratingscale (PSPRS) ( Golbe and Ohman-Strickland, 2007 ) and Addenbrooke's Cognitive Examination (ACE: Addenbrooke's Cognitive Examination-Revised for CCPP ( Mioshi et al., 2006 ), Addenbrooke's Cognitive Examination-III for PROSPECT-MR ( Hsieh et al., 2013 )). Summary scores and demographic details are outlined in Table 1 . PSP is a heterogeneous syndrome with variant presentations other than the classical Richardson's syndrome ( Jabbari et al., 2019 ). The clinical phenotype is related to the distribution of tau and focal grey matter loss, allowing us to test whether variation in the topographical distribution of disease burden ( Ling et al., 2014 ;Sakae et al., 2019 ;Tsuboi et al., 2005 ) or atrophy ( Jabbari et al., 2019 ) results in distinct changes in network dynamics. In keeping with Jabbari et al, clinical phenotype was categorized as PSP Richardson's syndrome (PSP-RS), PSP-subcortical (i.e., PSP-P with predominant parkinsonism or PSP-PGF with progressive gait freezing) or PSP-cortical (PSP-F with frontal presentations, PSP-CBS with corticobasal features or other focal cortical syndromes) ( Jabbari et al., 2019 ). Clinical phenotypes for both cohorts are included in Table 1 .
Image preprocessing used FSL's FEAT for registration to the structural image, motion correction, 100Hz high-pass temporal filtering, 5mm FWHM spatial smoothing. Denoising was performed using FSL's FIX with a training set of 10 subjects from each disease group per cohort. Additional removal of motion artefact used wavelet despiking (brainwavelet.org).
Given the sensitivity of estimates of network dynamics to participant motion ( Laumann et al., 2017 ;Nickerson et al., 2017 ;Power et al., 2012 ), we excluded thirty-three participants (27 PSP, 6 Control) with greater than 1 standard deviation from the whole sample mean for any of four data quality indices obtained from an fMRI dataset of 408 controls or participants with neurodegenerative diseases (maximum spike percentage ( Patel et al., 2014 ), median spike percentage, maximum framewise displacement ( Power et al., 2012 ), maximum spatial standard deviation of successive volume difference ( Smyser et al., 2010 )), and one participant with PSP for incomplete data. Summary measures by group are in supplementary table 1. We took the average of the four metrics post standardization as a covariate of no interest in further analyses.

Structural MRI
We extracted cortical thickness and subcortical grey matter volumes for 246 nodes of the Brainnetome Atlas ( Fan et al., 2016 )  and volumes for four brainstem substructures (midbrain, pons, medulla and superior cerebellar peduncle) ( Iglesias et al., 2015 ) using Freesurfer 6.0 ( Dale et al., 1999 ). We compared differences between participants with PSP and controls in thicknesses and volumes with permutation testing, family-wise error correction for multiple comparisons, and a statistical threshold of p < 0.05. Age and total intracranial volume were included as nuisance variables. We compared network dynamic metrics with atrophy measures, focusing on subcortical volumes and frontal cortical thicknesses given that neuropathological changes occur earlier and sequentially in these regions ( Kovacs et al., 2020 ).

Hidden Markov modelling
To investigate changes in network dynamics in PSP we first performed a group independent component analysis using FSL's MELODIC tool on preprocessed fMRI from controls and participants with PSP. We chose a model order of 30 as a balance between excessive network fragmentation ( Ray et al., 2013 ) and predetermining HMM outputs. The same component maps were used for both datasets. Participant specific timecourses for each component were generated by the first stage of a dual regression. From standardized per participant component timecourses a multivariate Gaussian HMM with 8 brain states was inferred using the HMM-MAR toolbox ( Vidaurre et al., 2017 ). A common covariance matrix was modelled for all states, so that the dynamic structure was driven primarily by activation changes. This approach was adopted to ensure that the model would effectively capture network dynamics in a cohort where there may be significant between-subject differences in static connectivity (see ( Vidaurre et al., 2021 ) for a summary of different possible approaches to modelling with HMM). The order selection in an HMM is predetermined; it has previously been shown that 8 states capturing large scale networks can provide robust behavioral inferences ( Vidaurre et al., 2018 ).
Different metrics can be used to quantify the temporal characteristics of the HMM states, including switching rate (the frequency with which states transitions occur) and fractional occupancy (the proportion of time a state is active). We assessed between-group differences in these metrics. Given the interdependence of fractional occupancy rates we performed a principal component analysis (PCA) to compare with severity measures.

Multiscale entropy
To investigate changes in complexity we calculated MSE for the same component timeseries used to infer the HMM, adapting LOFT's Complexity toolkit . We averaged over a fixed number of timescales (3 scales PROSPECT-MR, 4 for CCPP due to the longer time series), and calculated sample entropy on the time series constructed for each scale ( Costa et al., 2005 ). MSE is then sum of sample entropy across all timescales. We selected pattern length of 1 and pattern matching threshold of 0.35 given evidence that these parameters provide robust results . We took the average MSE calculated across the 30 component timeseries for further analyses. We assessed betweengroup differences and the correlation across-subjects between complexity and fractional occupancy.

Graph measures
We performed graph theoretical analysis using Maybrain software ( https://github.com/RittmanResearch/maybrain ) and Net-workX ( Hagberg et al., 2008 ), with the Brainnetome parcellation. Association matrices were constructed by taking the wavelet crosscorrelation between each region using a maximal overlap discrete wavelet transform and Daubechies filter performed using the waveslim package in R. The second band of 4 was used corresponding to a frequency range of 0.0675-0.125Hz ( Achard and Bullmore, 2007 ).
To test our hypothesis that changes in network dynamics were related cortical network topological remodelling in PSP in response to subcortical tau burden, we focused on the following metrics: weighted degree , measuring the number and strength of nodal functional connections; clustering coefficient , the proportion of triangular connections formed by each node over the proportion of all possible such connections; and path length , the average shortest topological distance between nodes of the graph. The combination of path length and clustering coefficient defines randomness or regularity of the graph, with random graphs exhibiting short path length and small clustering coefficient ( Watts and Strogatz, 1998 ). Path length and clustering coefficient were assessed across the brain, and weighted degree in cortical and subcortical regions and between groups. Metrics except for weighted degree were binarised after thresholding and normalised against 10 0 0 random graphs with identical degree distribution and random connections. A network density threshold of 5% was used. We also report results at density thresholds of 1%-10% for significant results to ensure robustness.

Statistical approach
We conducted initial analysis in the two cohorts separately. This was to allow analysis of MSE at higher scales in the CCPP cohort, thereby increasing our ability to differentiate complex signal from randomness and to contrast with HMM metrics, given that with fewer than 50 timepoints error of sample entropy es- timates may increase ( Yang et al., 2013 ;Yang et al., 2018 ). Replication can ensure that results are robust, an important factor given concerns that apparent changes in network dynamics from resting state fMRI are attributable to analysis techniques, head motion and sleep ( Laumann et al., 2017 ).
Statistical tests used a general linear model with permutation testing (10 0 0 0 permutations), with family-wise error correction for multiple comparisons and contrasts using FSL's PALM ( Winkler et al., 2014 ) and a statistical threshold of p < 0.05. The exceptions were moderation analysis, comparisons with graph metrics and direct tests of slope, which were performed in R (R Core Team, 2018 ). Participant motion, age and sex (estimated total intracranial volume for contrasts involving measures of volume) were included as nuisance variables.
Given demographic differences in both cohorts we repeated our analysis on a subset of participants, using the matchit function in R to select participants. We initially removed participants outside the common support region, followed by nearest neighbour matching with a caliper of 0.5 if groups remained unbalanced.

Demographics
We present analyses from 24 participants with PSP and 22 controls from CCPP, and 42 participants with PSP and 36 controls from PROSPECT-MR. Demographic details are outlined in Table 1 . There were significant differences in age in both cohorts and in gender in PROSPECT-MR.

Network dynamics
CCPP: We used temporally concatenated participant timeseries from ICA components to fit an HMM with 8 brain states. Mean activation maps for these states are shown in Fig. 1 A. There was no difference in switching rate between controls and participants with PSP (t = 0.37, p = 0.59). We performed a PCA of fractional occupancy rates, which are collinear and compositional. Three components with eigenvalues greater than 1 explained 75% of the variance and were taken forward for further analysis. The first component was significantly more negative in PSP than controls (t = 4.0, p = 0.0 0 08) ( Fig. 1 B). States with the highest positive loadings (states 5 and 7) had prominent subcortical and posterior activation, while states with the most negative loadings (states 1, 4 and 6) largely constituted the executive control and salience networks.
PROSPECT-MR: Mean activation maps for the 8 HMM inferred PROSPECT-MR brain states are shown in Fig. 1 D. States 1-4 were the closest Dice coefficient matches (supplementary figure 1) for both positive and negative maps. There were anatomical differences in the mean activation maps between the two cohorts, particularly in identified anti-correlations, resulting in divergence between designations for positive and negative maps. Therefore, we adjudicated matching for the remaining states by visual inspection.
There was no difference in switching rates between the two groups, but participants with PSP had altered fractional occupancy (see Fig. 1 E). Two components with eigenvalues greater than 1 explained 74% of the variance and were taken forward for further analysis. Component scores for the second component were significantly decreased in PSP (t = 3.1, p = 0.006). States with the most positive loadings (states 5 and 7) had prominent subcortical, posterior and motor region activations, while states with negative loadings overlapped with executive networks (states 1 and 6).
Given the demographic differences in both cohorts we repeated our analysis using a subset of participants matched for age and sex (see supplementary material). Our key findings were unchanged in both cohorts. In summary, in both groups participants with PSP spent a greater proportion of time in states of activity in executive networks, and away from states of activity in networks with posterior and subcortical activations.

Complexity
We investigated whether complexity differed between PSP and Controls, to provide further insight into temporal dynamics in the disease.
CCPP: We calculated MSE for each participant at 4 timescales using the same component timeseries used to infer the HMM. MSE was significantly reduced in PSP (t = 2.3, p = 0.022, Fig. 2 A).
PROSPECT-MR: MSE was calculated over 3 rather than 4 timescales, due to shorter time series. In contrast to our locally collected data, we did not find the reduction in MSE in PSP to be significant (t = 1.1, p = 0.27, Fig. 2 D)

Network dynamics and clinical severity
We tested whether the distinct changes in network dynamics in PSP were related to clinical severity, as measured by the PSPRS and ACE.
CCPP: The first fractional occupancy component negatively correlated with the PSPRS (r = -0.60, t = -3.1, p = 0.022, see Fig. 1 C). There was also a relationship between principal component 3 and ACE which was not significant after correction for multiple comparisons (r = 0.47, t = 2.4, uncorrected p = 0.03).
PROSPECT-MR: Component scores for the second fractional occupancy component correlated with PSPRS (r = -0.52, t = -3.7, p = 0.002, Fig. 1 F). There was also a relationship between the second fractional occupancy component and ACE which was not significant after correction for multiple comparisons (r = -0.34, t = -2.1, uncorrected p = 0.046).

Network dynamics and complexity
We investigated the relationship between signal complexity and network dynamics. We asked whether a) the distinct changes in fractional occupancy were related to complexity and b) switching rate related to complexity, and whether these relationships interacted with diagnosis.
3.5.1 Fractional occupancy and MSE . In both groups there was a significant relationship between MSE and the fractional occupancy component that differed between people with PSP and Controls (CCPP r = 0.44, t = 3.2, p = 0.004, Fig. 2 C; PROSPECT-MR r = 0.44, t = 4.2, p = 0.0 0 03, Fig. 2 F). Moderation analysis did not show any significant group differences in these relationships. In PROSPECT-MR the slope in controls was driven by a single outlier; post outlier removal the relationship between component and complexity differed by diagnosis (PROSPECT-MR r 2 = 0.08, F = 9.8, p = 0.003).

Switching rates and MSE:
In the CCPP group MSE did not correlate significantly ( Figs. 2 B) with switching rate (r = 0.25, t = 1.6, p = 0.11). In the PROSPECT-MR group a significant relationship was found (r = 0.57, t = 6, p = 0.0 0 01, Fig. 2 E). Moderation analysis did not show any significant group differences in these relationships.

Network dynamics in PSP versus structure, topology, and clinical presentation
We tested the hypothesis that people with PSP spend a greater proportion of time in inefficient states due to cortical remodelling in response to focal disease. Given that distribution of atrophy and pathology differ by phenotype, we tested whether network dynamics vary by clinical phenotype. Since the two datasets showed overlapping changes in brain state occupancy in PSP, we performed a combined analysis to investigate these hypotheses, focusing on fractional occupancy. We assessed whether fractional occupancy components: (i) relate to frontal cortical thickness and subcortical volumes; (ii) relate to regional topological remodeling; and iii) differed depending on PSP clinical phenotype.
We first sought to investigate the distribution of atrophy across the two cohorts. For participants with PSP we found significant areas of grey matter atrophy compared to controls in the midbrain and subcortical regions defined by the Brainnetome Atlas ( Fig. 3 ). We also found cortical atrophy primarily in the frontal lobe and peri-Rolandic regions, but also in the temporal and parietal cortex. There were no regional differences when comparing participants with PSP from the two cohorts. Fig. 3. Significant areas of grey matter volume reduction in PSP v controls in a combined analysis across the two cohorts, where regions are nodes of the Brainnetome Parcellation. p < 0.05 after family-wise error correction for multiple comparisons. There were no regional differences in direct comparison of participants with PSP from the two cohorts. In a combined analysis of the two cohorts the first two principal components differed between PSP and controls. Mean activation states with PCA loadings > |0.3| are shown. (C) and (D) Component 1 correlated with subcortical volume but not frontal cortical thickness, although with no significant difference in slope. Component 2 did not correlate significantly with either frontal cortical thickness or subcortical volume. E and F) Loadings in component 2 were associated with reduced clustering coefficient and reduced weighted degree in PSP but not controls. No relationships were found with component 1.
The two principal components derived from HMM analysis ( Figs. 4 A and 4 D) were taken forward for analysis differed between PSP and controls (component 1 t = 2.6, p = 0.023; component 2 t = -2.8, p = 0.014).
Despite the relationships between network dynamics, focal atrophy and topological changes, we found no significant difference by clinical phenotype (component 1 F = 3.4 p = 0.083; component 2 F = 0.57, p = 0.82).

Discussion
The principal results of this study are that (i) the exemplar tauopathy of PSP changes network dynamics, with a higher proportion of time spent in frontoparietal activation states; (ii) these changes in network dynamics are related to complexity as measured by multi-scale entropy, (iii) the changes in network dynamics correlated with clinical severity and regional atrophy; and (iv) altered network dynamics occur in the context of widespread changes to network topology. The effect of PSP phenotypic variance is expressed in terms of the relationship between network dynamics, clinical severity and focal atrophy, in cortical versus subcortical regions.
In two independent datasets, people with PSP spent more time in states whose spatial distributions mirror executive control networks. Given that in health, occupancy of networks associated with higher order cognition correlates positively with cognitive function ( Vidaurre et al., 2017 ), these results may seem surprising. Time in these networks did not correlate with frontal atrophy but did show a negative relationship with clustering coefficient and weighted degree only in PSP. This suggests a loss of small-world properties towards greater network randomness ( Bassett and Bullmore, 2016 ;Watts and Strogatz, 1998 ), with occupancy of a remodeled and more random network no longer relating to its effective functioning.
PSP causes severe disruption to connectivity between the subcortex/brainstem and cortical regions ( Brown et al., 2017 ;Whitwell et al., 2011 ), changes that may account for our finding of reduced time in states with activity in subcortical regions, which correlated with subcortical atrophy. So, why do participants with PSP spend more time in executive control networks and less time in states representing the default mode network with negative activations in regions of the task-positive network? Regional structural network properties influence brain state transitions, with access to frontoparietal cognitive control networks depending on nodes within weakly connected regions ( Gu et al., 2015 ). Key nodes in a dysfunctional random network may no longer effectively determine state transitions, causing altered network dynamics in regions remote from the primary pathology. Specifically, the results shed light on the interplay in disease between atrophy, remodeling of network topology and network dynamics.
We found that network dynamics in PSP varied by disease severity, both in terms of relationship to atrophy and to the PSPRS: the latter is sensitive to disease progression and predicts survival ( Bang et al., 2016 ;Golbe and Ohman-Strickland, 2007 ). We did not find that this dependency translated into differences between PSP phenotypes, although our study was powered to only detect large subgroup differences. Nonetheless, given that the intrinsic network architecture of the brain present at rest shapes that seen during tasks ( Cole et al., 2014 ;Smith et al., 2009 ), and that network structure predicts behavioral traits ( Arbabshirani et al., 2017 ;Beaty et al., 2018 ;Meer et al., 2020 ;Rosenberg et al., 2016 ), we hypothesize that the observed changes in brain state transitions in PSP underpin cognitive symptoms of the disease.
We have shown that measuring complexity provides a complementary method to assess temporal dynamics in disease. In the CCPP dataset we found reduced complexity in PSP, with the greatest differences seen at the highest scale, indicating a true reduction in complexity and not randomness. This result was not seen in PROSPECT-MR, perhaps because these images were acquired with fewer time points. We found that complexity correlated with both switching rate and fractional occupancy component, suggesting that changes in global signal influence network dynamics in PSP. This raises the possibility of additional aetiological factors to those outlined above, including the profound neurotransmitter deficits which occur early in PSP ( Murley and Rowe, 2018 ) and alter global fMRI signal ( Turchi et al., 2018 ) and dynamic connectivity ( Shine et al., 2018 ) in health.
There has been controversy as to whether variability in dynamic functional connectivity in task-free fMRI is more than motion ( Laumann et al., 2017 ) . We used stringent exclusion criteria to limit the impact of movement artefact. A key advantage of this study is that results are replicated in two datasets, an important part of the solution to non-generalizable results in neuroimaging due to low statistical power and analytic flexibility ( Poldrack et al., 2017 ). Our approach did however result in significant numbers of exclusions, particularly from our CCPP dataset (due to the longer time series). It may be that these excluded participants have distinct clinical phenotypes, and therefore our conclusions would not apply to all individuals with PSP. Modelling network dynamics via an HMM requires forced choices in analysis, notably the number of states inferred. If brain states consist of a hierarchy of structures ( Vidaurre et al., 2018 ) it is likely that HMMs can provide multiple related solutions. Indeed, this may account for some of the anatomical differences observed in our two datasets. We believe that this challenge is best tackled by independent replication and relating findings to clinical scores.

Conclusions
The quantification of the dynamics of large-scale resting state networks offers an intermediate phenotype with which to understand clinical syndrome. Through hidden Markov modelling we have shown that changes in network dynamics relate to neural signal complexity and to phenotypic variance in progressive supranuclear palsy. Our methodology demonstrates how abnormalities in regional atrophy and topological changes correlate with brain state transitions, and provides a means to directly test the causes and consequences of altered network dynamics.

Verification
We verify that this work has not been published before, is not under consideration for publication elsewhere, and if accepted will not be published elsewhere in the same form. Submission for publication is approved by all authors.

Disclosure statement
J.B.R. serves as an associate editor to Brain and is a nonremunerated trustee of the Guarantors of Brain and the PSP Association (UK). He provides consultancy to Asceneuron, Biogen, and UCB and has research grants from AZ-Medimmune, Janssen, and Lilly as industry partners in the Dementias Platform UK. M.T.H received payment for Advisory Board attendance/consultancy for Biogen, Roche, CuraSen Therapeutics, Evidera, Manus Neurodynamica and the MJFF Digital Health Assessment Board. M.T.H is a coapplicant on a patent application related to smartphone predictions in Parkinson's (PCT/GB2019/052522) pending. All other authors did not declare any funding sources that directly contributed to this study.