Impaired long-range excitatory time scale predicts abnormal neural oscillations and cognitive deficits in Alzheimer’s disease

Abstract

Background

Alzheimer’s disease (AD) is the most common form of dementia, progressively impairing cognitive abilities. While neuroimaging studies have revealed functional abnormalities in AD, how these relate to aberrant neuronal circuit mechanisms remains unclear. Using magnetoencephalography imaging we documented abnormal local neural synchrony patterns in patients with AD. To identify global abnormal biophysical mechanisms underlying the spatial and spectral electrophysiological patterns in AD, we estimated the parameters of a biophysical spectral graph model (SGM).

Methods

SGM is an analytic neural mass model that describes how long-range fiber projections in the brain mediate the excitatory and inhibitory activity of local neuronal subpopulations. Unlike other coupled neuronal mass models, the SGM is linear, available in closed-form, and parameterized by a small set of biophysical interpretable global parameters. This facilitates their rapid and unambiguous inference which we performed here on a well-characterized clinical population of patients with AD (N = 88, age = 62.73 +/- 8.64 years) and a cohort of age-matched controls (N = 88, age = 65.07 +/- 9.92 years).

Results

Patients with AD showed significantly elevated long-range excitatory neuronal time scales, local excitatory neuronal time scales and local inhibitory neural synaptic strength. The long-range excitatory time scale had a larger effect size, compared to local excitatory time scale and inhibitory synaptic strength and contributed highest for the accurate classification of patients with AD from controls. Furthermore, increased long-range time scale was associated with greater deficits in global cognition.

Conclusions

These results demonstrate that long-range excitatory time scale of neuronal activity, despite being a global measure, is a key determinant in the local spectral signatures and cognition in the human brain, and how it might be a parsimonious factor underlying altered neuronal activity in AD. Our findings provide new insights into mechanistic links between abnormal local spectral signatures and global connectivity measures in AD.

Introduction

Alzheimer’s disease (AD) is the most common form of dementia, progressively impairing the cognition and behavior of the affected individual. It has been proposed that the effect of AD neurodegeneration on cortical neuronal networks is partially reflected by the abnormal mechanisms of cortical neural synchronization and coupling [1]. Neural synchronization refers to the simultaneous activity of neuronal groups in the brain. Repetitive spiking activities of neural populations manifests as oscillations ranging from slow delta to fast gamma frequencies. Synchronization of neural oscillations may represent both local synchrony , typically estimated from regional level power spectral density (PSD) of the electrophysiological signal, and long-range synchrony, estimated from pair-wise coherences between signals originating at different locations. Neurodegeneration disrupts both local and long-range synchrony [2, 3]. Because functional deficits precede the structural deficits in AD, it is likely that local and long-range synchrony deficits occur before the onset of clinical symptoms, may worsen as the disease progresses, and may even play a role in disease manifestation [4, 5]. It is therefore important to understand how synchrony within and between brain regions is disrupted in AD and is associated with cognitive impairment.

Electro- or magneto-encephalography (E/MEG), that capture temporal functional activity scales with millisecond precision, studies have shown thats both local and long-range synchrony are abnormal in AD [6, 7].The relative PSD of patients with AD is significantly increased in the delta and the theta, while reduced in the alpha frequency band [2, 8,9,10,11,12,13], often referred to as an oscillatory slowing. This oscillatory slowing has also been observed in MEG studies by our group and others [6, 14,15,16,17,18,19,20,21]. Interestingly, these aberrant local synchrony patterns are associated with the pathological processes in AD [22]. These data have naturally led to a search for a common underlying neural mechanism whose impairment might account for observed spatial and spectral shifts in E/MEG imaging data in AD patients.

In this paper, we hypothesize that global alterations in neural mechanisms for synchrony can account for the abnormal spatiospectral profile of neural oscillations across the brain as observed in MEG, and predict cognitive decline in patients with AD. We test this hypothesis using a biophysical model of whole-network level brain activity, called the Spectral Graph Model (SGM), which can capture these phenomena with a parsimonious set of biophysically interpretable parameters. SGM posits that the anatomical network of fiber projections is a key substrate that underlies the emergence of spatially- and spectrally-patterned alterations in functional activity [23]. At its core, this model incorporates the effect of short- and long-range communication between cortical neuronal populations supported by the anatomical network. To test the proposed hypothesis, we performed a thorough parameter inference of the SGM on individual subjects’ source reconstructed MEG data and assessed parameters from a well-characterized clinical population of AD patients and a cohort of age-matched healthy controls. Consistent with our hypothesis, we demonstrate that a global slowing of long-range excitatory time scale is predictive of AD spatiospectral patterns and cognitive decline.

Methods

Data description

Eighty-eight patients with AD (diagnostic criteria for probable AD or mild cognitive impairment due to AD) [24,25,26] and 88 age-matched controls were included in this study. Each participant underwent a complete clinical history, physical examination, neuropsychological evaluation, brain magnetic resonance imaging (MRI), and a 5-10-minute session of resting MEG. All participants with AD were recruited from research cohorts at the University of California San Francisco-Alzheimer’s Disease Research Center(UCSF-ADRC). Healthy control participants were recruited at UCSF-ADRC as well as from several ongoing studies at the Biomagnetic Imaging Laboratory at UCSF. Informed consent was obtained from all participants and the study was approved by the Institutional Review Board (IRB) at UCSF (UCSF-IRB 10-02245). The mean (std) age of controls and patients with AD was 65.07 (9.92) and 62.73 (8.64) years, respectively. 51 (58%) of 88 controls, and 53 (60.2%) of patients with AD were females. The mean (std) MMSE score of patients with AD was 22.14 (5.55), while the mean Clinical Dementia Rating-Sum of Boxes (CDR) score of patients with AD was 4.90 (2.75).

Clinical assessments and MEG, and MRI acquisition and analyses

All the processing pipelines are the same as that for a previous study [20]. Patients with AD were assessed via MMSE and a standard battery of neuropsychological tests. Patients with AD were assessed via a structured caregiver interview to determine the Clinical Dementia Rating.

MEG scans were acquired on a whole-head biomagnetometer system (275 axial gradiometers; MISL, Coquitlam, British Columbia, Canada) for 5-10 min, following the same protocols described previously [20, 27]. Tomographic reconstructions of source-space data were done using a continuous 60-second data epoch, an individualized head model based on structural MRI, and a frequency optimized adaptive spatial filtering technique implemented in the Neurodynamic Utility Toolbox for MEG (NUTMEG; http://nutmeg.berkeley.edu). We derived the regional power spectra based on Desikan-Killiany atlas parcellations for the 68 cortical regions depicting neocortex and allocortex, the latter including the entorhinal cortex. Regional power spectra were derived from FFT and then converted to dB scale.

Resting state MEG data acquisition

Each subject underwent MEG recording on a whole-head biomagnetometer system consisting of 275 axial gradiometers (MISL, Coquitlam, British Columbia, Canada), for 5-10 min. Three fiducial coils including nasion, left and right preauricular points were placed to localize the position of head relative to sensor array, and later coregistered to each individual’s respective MRI to generate an individualized head shape. Data collection was optimized to minimize within-session head movements and to keep it below 0.5 cm. 5-10 min of continuous recording was collected from each subject while lying supine and awake with eyes closed (sampling rate: 600 Hz). We selected a 60s (1 min) continuous segment with minimal artifacts (minimal excessive scatter at signal amplitude <10 pT), for each subject, for analysis. The study protocol required the participant to be interactive with the investigator and be awake at the beginning of the data collection. Spectral analysis of each MEG recording and whenever available, and the simultaneously collected scalp EEG recording were examined to confirm that the 60-s data epoch represented awake, eyes closed resting state for each participant. Artifact detection was confirmed by visual inspection of sensor data and channels with excessive noise within individual subjects were removed prior to analysis.

Source space reconstruction of MEG data and spectral power estimation

Tomographic reconstructions of the MEG data were generated using a head model based on each participant’s structural MRI. Spatiotemporal estimates of neural sources were generated using a time-frequency optimized adaptive spatial filtering technique implemented in the Neurodynamic Utility Toolbox for MEG (NUTMEG; https://nutmeg.berkeley.edu/). Tomographic volume of source locations (voxels) was computed through an adaptive spatial filter (8-mm lead field) that weights each location relative to the signal of the MEG sensors [28, 29]. The source space reconstruction approach provided amplitude estimations at each voxel derived through the linear combination of spatial weighting matrix with the sensor data matrix [28]. A high-resolution anatomical MRI was obtained for each subject (see below) and was spatially normalized to the Montreal Neurological Institute (MNI) template brain using the SPM software (http://www.fil.ion.ucl.ac.uk/spm), with the resulting parameters being applied to each individual subject’s source space reconstruction within the NUTMEG pipeline [29].

To prepare for source localization, all MEG sensor locations were coregistered to each subject’s anatomical MRI scans. The lead field matrix (forward model) for each subject was calculated in NUTMEG using a multiple local-spheres head model (three-orientation lead field) and an 8-mm voxel grid which generated more than 5000 dipole sources, all sources were normalized to have a norm of 1. The lead field matrix describes transformation from a unit dipolar source to observed sensor patterns. Columns of the lead field matrix (also called gain vectors or forward kernels) describe the spatial sensitivity pattern of a single dipolar source across all sensors. Here we use 3 orientations for each dipolar source and that is what is meant by a three-orientation lead-field matrix for each source. The MEG recordings were projected into source space using a beamformer spatial filter. Source estimates tend to have a bias towards superficial currents and the estimates are more error-prone when we approach subcortical regions, therefore, only the sources belonging to the 68 cortical regions were selected for further analyses. We used a spatial resolution of 8mm between dipolar sources. Specifically, all dipole sources were labeled based on the Desikan-Killiany parcellations, then source time-courses for each orientation of dipoles within each region were averaged together to obtain regional time-courses. In this study, we examined the broad-band (1-35 Hz). Power spectra were derived by applying FFT on the time-course data and then converted to the dB scale.

The dual signal subspace projection (DSSP) preprocessing denoising was run on 9 controls and 13 patients with AD where large artifacts were seen in MEG data [30]. No other data noise reduction procedures were done.

Magnetic resonance image acquisition and analysis

Structural brain images were acquired from all participants using a unified MRI protocol on a 3 Tesla Siemens MRI scanner at the Neuroscience Imaging Center (NIC) at UCSF. Structural MRIs were used to generate individualized head models for source space reconstruction of MEG sensor data. Structural MRI scans were also used in the clinical evaluations of patients with AD to identify the pattern of gray matter volume loss to support the diagnosis of AD.

Extraction of spectral peaks

Spectral features of the MEG PSD were extracted using the FOOOF toolbox [31]. We used the FOOOF.fit() function to extract the first two peaks of the periodic component of the PSD. If only one peak was present, the second peak was assigned the same value as the first one. We set the lower limit of the peak_width_limits to be 2 Hz. All the other settings were the same as default.

Model

SGM provides a closed-form solution of the steady-state frequency response of different brain regions. Here, we use a Desikan-Killiany parcellation scheme [32] to estimate the brain regions. The SGM is characterized by 7 parameters, which are either global or local but spatially-invariant. These parameters include the spatially-invariant local synchrony-related time constants $\tau _{\textrm{e}}, \tau _{\textrm{i}}$ and spatially-invariant but local neural gains $g_{\textrm{ei}}, g_{\textrm{ii}}$ as a measure of overall synaptic strength at the local scale for both excitatory and inhibitory neuronal subpopulations; as well as a global excitatory time constant $\tau _{\textrm{G}}$ at the long-range scale representing the long-range network connections, global coupling constant $\alpha$, and speed of transmission of signals among regions v. SGM is characterized by an additional local excitatory gain parameter $g_{\textrm{ee}}$ but it is set to 1 to ensure parameter identifiability. Each region is assumed to consist of local excitatory and inhibitory neuronal subpopulations that interact with each other and regulate the long-range excitatory neuronal populations. The long-range populations are assumed to be connected to each other via the structural connectome. Here, we use a template structural connectome from the Human Connectome Project (HCP). Hence the model entails no features that may change from region to region, except of course features from the heterogeneously connected anatomical network. The structural connectivity matrix is shown in Supplementary Fig. S5 and the distance matrix is shown in Supplementary Fig. S6. To infer the SGM parameters, we fit SGM output to the frequency spectra obtained from MEG for healthy controls and AD subjects. The model used here is similar to the SGM developed previously [33,34,35,36], and is described in detail in the Supplementary document.

The model solution can be obtained in a closed form in the frequency domain as a function of angular frequency $\omega$ as:

$$\begin{aligned} \varvec{X}(\omega ) = \sum \limits _{k=1}^N \frac{\varvec{u}_k(\omega )\varvec{u}_k(\omega )^\mathrm{H}}{\mathrm{j} \omega + \tau _{\mathrm{G}}^{-1}\lambda _k(\omega ) F_{\mathrm{G}}(\omega )}H_{\mathrm{local}}(\omega )\varvec{P}(\omega ) \, , \end{aligned}$$

(1)

where, $\varvec{X}(\omega )$ is a vector of the long-range pyramidal neuronal population signal of every brain region, $\varvec{u}_k(\omega )$ are the eigenmodes and $\lambda _k(\omega )$ are the eigenvalues obtained by the eigen-decomposition of a complex Laplacian matrix, N is the total number of brain regions (86 in this case), $H_{\textrm{local}}(\omega )$ is the transfer function of the local excitatory and inhibitory interactions, $\varvec{P}(\omega )$ is the Fourier transform of white Gaussian noise which is an input to the local excitatory and inhibitory signals, $F_{\textrm{G}}(\omega )$ is the Fourier transform of a Gamma-shaped neural ensemble response function (derived in the Supplementary material). Equation (1) is the closed-form steady-state solution of the long-range signals at a specific angular frequency $\omega$. We use this modeled spectra to compare against empirical MEG spectra and subsequently estimate model parameters. In practice, only a few eigenmodes $k\in [1,K], K\ll N$ are needed to obtain sufficiently strong fits to empirical data, including especially the lowest eigenmodes [34].

Model parameter estimation

The model parameter estimation procedure is same as described previously [36]. Modeled spectra was converted into PSD by calculating the norm of the frequency response and converting it to dB scale by taking $20\text {log}_{10}()$ of the norm. Pearson’s r between modeled PSD and the MEG PSD was used a goodness of fit metric for estimating model parameters. Pearson’s r between modeled and MEG PSD was computed for all 68 brain regions. Its average r across all regions is referred to as the spectral correlation. Next we calculated the spatial correlation by obtaining the regional distribution of alpha band (8-12 Hz) raw power of both model $\varvec{x}$ and MEG $\varvec{y}$. Then, the spatial correlation was defined as $\varvec{x}^T \, \Vert (\varvec{C}+w\varvec{I})\Vert \, \varvec{y}$, where $\varvec{C}$ is the row degree normalized structural connectivity matrix, $\varvec{I}$ is the identity matrix, w is an empirical weight, and $\Vert (\varvec{C}+10\varvec{I})\Vert$ is the row normalized version of $\varvec{C}+ 10\varvec{I}$. The objective function for optimization and estimation of model parameters was the sum of spectral and spatial correlations. We used a dual annealing optimization procedure in Python for performing parameter optimization [37].

Parameter initial guesses and bounds for estimating the static spectra are specified in Table 1. We defined three different bounds on the neural gain terms to ensure that the model is stable, based on prior work on model stability [35]. First, we supplied a larger bound on the neural gains for optimization. If the optimal model parameter was outside the stability boundary, we repeated optimization with a smaller bound. We repeated this procedure 3 times to ensure that the final optimal model parameters correspond to the stable model solutions. We used a dual annealing optimization procedure in Python for parameter optimization [37]. The dual annealing optimization was performed for three different initial guesses, and the parameter set leading to maximum sum of spectral and spatial correlations was chosen for each subject. The dual annealing settings were: maxiter = 500. All the other settings were the same as default.

Table 1 SGM parameter values, initial guesses, and bounds for parameter estimation for static spectra fitting

Statistical analyses

Statistical tests were performed using SAS software (SAS9.4; SAS Institute, Cary, NC) and the statsmodels package in Python. To compare the neuronal parameters between the controls and patients, we used a linear mixed-effects model (PROC MIXED), to compare model parameters ($\tau _{\textrm{G}}, \tau _{\textrm{e}}, \tau _{\textrm{i}}, g_{\textrm{ii}}, g_{\textrm{ei}}, \alpha , v$), including age as a covariate into the models. We reported the estimated least-squares means and the statistical differences of least-squares means based on unpaired t-tests. We corrected for multiple testing using a Bonferroni correction on the 7 comparisons. We also developed univariate linear regression models to examine the associations between model parameters and MMSE and CDR scores in AD. In these models, the dependent variables included MMSE and CDR (in separate models), and the predictor variables included the model parameters we found significant between AD and controls ($\tau _{\textrm{G}}$, $\tau _{\textrm{e}}$, and $g_{\textrm{ii}}$). We also corrected for multiple testing using a Bonferroni correction on the 3 comparisons. Next, we developed multivariate linear regression models with dependent variables as MMSE and CDR (separately), and the predictor variables included all the significant parameters $\tau _{\textrm{G}}$, $\tau _{\textrm{e}}$, $g_{\textrm{ii}}$, and age as covariates.

Classification between AD and controls

We trained a random forest for classifying AD and controls. Here, we used the SGM parameters and age as features of the model. For training and testing, we employed a 5-fold stratified cross validation method. We divided the dataset into 5 folds and used the 4 folds for training, and the 5^th fold for testing the model. We repeated this procedure 100 times. While training the model, no information of the testing fold was provided. With the 4 folds of the training dataset, further 5-fold cross validation was performed to estimate the tuning parameter of the random forest. Here, we only tuned for the max depth with the following options for max depth: None, 2, 3, 4. All other hyperparameters were kept as default in the sklearn package in Python. After estimating the tuning parameter, the model was trained using the entire training dataset and then tested on the 5^th fold. The mean AUROC of the test dataset was finally reported. The feature importance was estimated as the average of the feature importance from the random forest classifier that was trained 100 times.

For training the classifier with hippocampal volume or MMSE as features, we only used a subset of subjects since this data was not available for all the subjects (56 controls and 51 patients with AD with hippocampal volume data; 66 controls and 87 patients with AD with MMSE data). Since the controls and patients with AD were not age-matched when using hippocampal volume or MMSE as features, we did not include age as an additional feature for classification. For training the classifier with MEG PSD as features, we extracted the exponent of the aperiodic component and the center frequency and power of the first two peaks of every region’s PSD using the FOOOF toolbox [31]. This resulted in a total of $5\times 68 = 340 + 1$ (for age) features. We set the lower limit of the peak_width_limits in FOOOF to be 2 Hz. All the other settings were the same as default.

Results

Global timescales of empirical MEG electrophysiological recordings are longer in patients with AD

First, we show that the global timescales of empirical MEG recordings are longer in patients with AD. To this end, we collected MEG recordings for 88 patients with AD and 88 age-matched healthy controls. To evaluate the global timescales of the MEG recordings, we first calculated the autocorrelation function of the bandpass filtered MEG time-series for every brain region and then took a mean across all the brain regions (Fig. 1A). From this averaged autocorrelation function, we obtained an exponential decay time constant ($\tau$ in $e^{-t/\tau }$) that captures the timescale of a time series. Precisely, this is the value of the lag time for which the autocorrelation function is $e^{-1}$ (see [38] and the references therein for more details). As seen in Fig. 1B, the time constant is significantly longer for AD based on a Kolmogorov-Smirnov test, implying that the MEG signals of patients with AD decay slower. To supplement this result, we obtained the mean PSD over all regions (Fig. 1C) and compared the center frequencies corresponding to the first and the second peaks of the PSD in C (Fig. 1D). As seen in Fig. 1D, the center frequencies of both the first and the second peaks of the PSD are lower in AD, based on a Kolmogorov-Smirnov test. This decrease in the center frequencies of the peaks also indicates a slowing of the timescales in AD. To investigate the underlying mechanisms of this slowing, we inferred SGM parameters for the two cohorts as explained in the methods section and Supplementary Fig. S1.

Fig. 1

A: Autocorrelation function from MEG timeseries. Each autocorrelation function is a mean of autocorrelation functions of all the brain regions, and therefore each line corresponds to a single subject. The timescale is the value of lag that corresponds to an autocorrelation function value of $e^{-1}$. B: Distribution of timescales obtained from autocorrelation function. Based on a Kolmogorov-Smirnov test, time constant of MEG recordings for patients with AD is significantly larger than the timescale of MEG recordings for healthy controls ($p<0.001$, Cohen’s D effect size = 0.42). C: Power spectral density (PSD, in dB scale) extracted from the MEG recordings. Each PSD is a mean of all regions and then centered to the mean and scaled to unit variance for every subject separately. D: Center frequencies of the first and the second peaks of the PSD in C for every subject. Center frequencies of both the first and the second peaks are lower in AD ($p<0.001$, Cohen’s D effect size = 0.21 for the first peak, $p<0.001$, Cohen’s D effect size = 0.32 for the second peak). The peaks were extracted using the FOOOF toolbox [31]

SGM reliably reproduces the spectral and spatial patterns of power spectral density

The predicted spectra from SGM reliably captured the empirical MEG spectra from patients with AD and age-matched controls (Fig. 2A; The mean (std) spectral correlations were 0.72 (0.08) and 0.78 (0.09) for controls and AD, respectively, shown in Fig. 2C). Compared to age-matched controls, patients with AD showed a lower alpha peak and a higher spectral power within the low-frequency delta-theta range (2-7 Hz), in their empirical spectral recording from MEG. This characteristic spectral change is clearly replicated in the predicted spectra derived from SGM. The spatial distribution of spectral power density of the alpha band, as expected, showed a postero-anterior distribution in both controls and patients. The spatial patterns of the predicted spectra from SGM reproduced this postero-anterior distribution with high fidelity (Fig. 2B; the mean (std) spatial correlations were 0.60 (0.09) and 0.66 (0.09) for controls and AD, respectively, shown in Fig. 2D). The region-wise spectral correlations and the frequency band-specific spatial correlations are also shown in Supplementary Fig. S2. These correlations were greater than 0.5 for more than 90% of subjects, regions, and frequency bands overall. Sample empirical versus modeled PSD for a single subject is shown in Supplementary Fig. S3. Sample spatial patterns in the beta frequency band are shown in Supplementary Fig. S4.

Fig. 2

A: Comparison of empirical (left) and SGM (right) frequency spectra for controls (top) and patients with AD (bottom). The darker lines correspond to the PSD averaged over all regions and subjects. The shaded region corresponds to the 90% confidence interval for the mean PSD (over all subjects) for different regions. B: Spatial distribution of the empirical (left) and SGM (right) alpha frequency band, for subjects with mean spatial correlations in controls (top) and patients with AD (bottom). The color scale of each spatial distribution was chosen based on their dynamic range. C: Spectral correlations of model fitting for controls and patients with AD. D: Spatial correlations of model fitting for controls and patients with AD

Patients with AD have altered network time constants and neural gains

Next, we compared the network parameters derived from the SGM between patients with AD and age-matched controls. Recall that these parameters are either global or local but assumed spatially-invariant. These parameters include the spatially-invariant local synchrony-related time constants $\tau _{\textrm{e}}, \tau _{\textrm{i}}$ and spatially-invariant but local neural gains $g_{ei}, g_{ii}$ as a measure of overall synaptic strength at the local scale for both excitatory and inhibitory neuronal subpopulations; as well as a global excitatory time constant $\tau _{\textrm{G}}$ at the long-range scale representing the long-range network connections, global coupling constant $\alpha$, and speed of transmission of signals among regions v. Using a general linear model with age included as a covariate, we found that patients with AD have significantly elevated long-range excitatory time constant ($\tau _{\textrm{G}}$; controls mean = 7.50 ms, confidence interval = (6.32 ms, 8.68 ms), AD mean = 13.90 ms, confidence interval = (12.72 ms, 15.09 ms), Cohen’s D effect size = 1.16), local excitatory time constant, ($\tau _{\textrm{e}}$; controls mean = 11.88 ms, confidence interval = (10.36 ms, 13.41 ms), AD mean = 15.01 ms, confidence interval = (13.48 ms, 16.53 ms), Cohen’s D effect size = 0.41) and local inhibitory neural gain ($g_{\textrm{ii}}$; controls mean = 0.26, confidence interval = (0.16, 0.36), AD mean = 0.46, confidence interval = (0.36, 0.56), Cohen’s D effect size = 0.42; Fig. 3A, B, and C). The highest effect size among the parameter comparisons was found in $\tau _{\textrm{G}}$ between AD and controls. Collectively these results indicate that while global network parameters are altered in AD, long-range excitatory connections may reflect such changes with greater sensitivity than other parameters.

Fig. 3

A, B, C: Statistical significance testing of difference in model parameters between AD and controls, with age as a covariate. Distribution of A: long-range parameters $\tau _{\textrm{G}}$ (long-range excitatory time constant), v (speed), and $\alpha$ (coupling constant); B: local time constants $\tau _{\textrm{e}}$ (excitatory) and $\tau _{\textrm{i}}$ (inhibitory); and C: local neural gains $g_{\textrm{ii}}$ (inhibitory gain) and $g_{\textrm{ei}}$ (gain of signals from the coupling between excitatory and inhibitory neurons). P-values are reported after correcting for multiple testing using a Bonferroni correction. $*$: $p<0.05$, $***$: $p<0.001$. D, E, F: Univariate associations of D: $\tau _{\textrm{G}}$, E: $\tau _{\textrm{e}}$, and F: $g_{\textrm{ii}}$ with MMSE in patients with AD. G, H: Classification of AD vs controls with a random forest classifier with SGM parameters and age as features of the classifier. G: ROC curve for classification of AD versus controls. H: Feature importance plot of SGM parameters. I: Spectral correlations when optimizing for only certain model parameters while keeping the others fixed at the average of the optimized model parameters of both AD and controls. “None” implies that all the model parameters were fixed at the average. The second point on the x-axis with the label $\tau _{\textrm{G}}$ implies that only $\tau _{\textrm{G}}$ was allowed to be optimized while the other model parameters were fixed at the average values. The third point on the x-axis with the label $\tau _{\textrm{e}}$ implies that both $\tau _{\textrm{G}}$ and $\tau _{\textrm{e}}$ were allowed to be optimized while keeping the other model parameters fixed at the average values. All the subsequent points on the x-axis correspond to similarly including more model parameters in optimization, based on their importance in the classification of AD vs controls

Altered long-range excitatory connections are correlated with global cognitive deficits in patients with AD

To investigate the association between altered global network parameters and cognitive deficits in patients with AD, we examined the correlations between $\tau _{\textrm{G}}$, $\tau _{\textrm{e}}$, and $g_{\textrm{ii}}$ with global cognitive decline measured by Mini Mental State Exam (MMSE), and overall disease severity measured by clinical dementia rating sum of boxes (CDR), in patients with AD. We first tested for univariate associations between the model parameters and MMSE and CDR separately, using linear regression. After adjusting for multiple testing (Bonferroni), $\tau _{\textrm{G}}$ showed significant negative associations with MMSE (Fig. 3D) where higher $\tau _{\textrm{G}}$ predicted greater cognitive deficits in MMSE. Next, we tested for the association between $\tau _{\textrm{G}}$ and MMSE including $\tau _{\textrm{e}}$, $g_{ii}$, and age as covariates in a multivariate linear regression model. This multivariate analysis also showed a significant negative association between $\tau _{\textrm{G}}$ and MMSE only ($p = 0.007$ for the association between $\tau _{\textrm{G}}$ and MMSE, model $r = 0.402$, model adjusted $r^2 = 0.121$, $F = 3.961$). Similar to the univariate results, none of the parameters were significantly associated with CDR in a multivariate regression model after adjusting for multiple testing. Details of statistics are mentioned in the Supplementary Tables S1 and S2.

Altered global network parameters can distinguish between AD and controls with high accuracy

Next, we examined the sensitivity and specificity of altered global network parameters to distinguish between patients with AD and controls. To this end, we trained and tested a random forest classifier including the model parameters and age as the classifier features. The average AUC of the ROC curves from the testing folds is 0.85, with a standard deviation of 0.02 (Fig. 3G). The other classification metrics included: accuracy $=$ 0.78, precision $=$ 0.79, recall score $=$ 0.75, and f1 score $=$ 0.77, on average. The confusion matrix is shown in Supplementary Fig. S7. We also obtained the feature importance score of the features used in training the model, shown in Fig. 3H. The time constant $\tau _G$ was the most important feature in classifying AD versus controls. Collectively, these results indicate that altered global network parameters are reliable indices to identify patients with AD from their age-matched counterparts and that long-range excitatory connections are the most sensitive indicators of AD-related global network deficits.

We also compared the performance of the SGM parameters against other measures, shown in the Supplementary Fig. S8. First, we used MEG PSD features directly to train a classifier. The mean AUC of the ROC curves from the testing folds of this classifier was 0.9. This is expected because SGM is fitting to this data and therefore cannot perform better than the data it is fitting to. Second, we used MMSE as a feature to train a classifier separately. The mean AUC of the ROC curves of this classifier was 0.95. This is also expected because MMSE scores closely reflect the Clinical Dementia Rating scores which were used to classify AD from controls in the first place. Third, we used mean hippocampal volume as a feature to train a classifier separately. The mean AUC of the ROC curves of this classifier was 0.73. Note that the aim of this study is not to obtain the best classifiers of AD, but to identify “biophysically interpretable” markers of AD. Therefore, even though the classification mean AUC drops when using SGM parameters instead of MEG PSD or MMSE as features, using SGM parameters have an added advantage of interpretability. Also note that the classifier needs to be trained on a higher dimensional data when using MEG PSD (340 features), as compared to using SGM parameters (7 features).

Minimal set of altered model parameters capture the empirical power spectral density patterns in AD

In order to assess the importance of model parameters in capturing the empirical PSD, we evaluated the spectral and spatial correlations after optimizing for certain model parameters, based on their importance from Fig. 3H, while keeping the remaining model parameters as the average of all the optimized model parameters for AD and controls together. First, we evaluated the correlations when none of the model parameters are optimized for and are all the average of the optimal parameters obtained previously. Next, we optimized only for $\tau _{\textrm{G}}$ while keeping all the other model parameters fixed since $\tau _{\textrm{G}}$ was the most important parameter in the classification of AD vs controls. Subsequently, we optimized for both $\tau _{\textrm{G}}$ and $\tau _{\textrm{e}}$ while keeping the remaining model parameters fixed since $\tau _{\textrm{e}}$ was the second most important feature in classification. We repeated this procedure till we included all the model parameters for optimization. The spectral correlations from this evaluation are reported in Fig. 3I. As seen in the figure, we see a sharp increase when $\tau _{\textrm{G}}$ is allowed to vary while keeping the other model parameters fixed. Upon including the subsequent model parameters, we do not see a substantial increase in the spectral correlation. This result strengthens our prior observation on the importance of $\tau _{\textrm{G}}$ in differentiating AD from controls. Note that we did not see any substantial difference in the spatial correlations.

Discussion

In this study, we determined local and long-range neuronal parameters of a computational model of brain activity that can account for abnormal neurophysiological activity in AD observed in high spatio-temporal resolution MEG imaging. We used SGM, which is ideally suited for this exploration providing a computational link between structure and function in the brain. The neuronal time constant associated with long-range excitatory connections was the most sensitive biophysical property that mediates abnormal global network dynamics in patients with AD. The long-range excitatory time constant not only predicted the global cognitive deficits in AD patients but also classified AD versus controls with high accuracy. To our knowledge, this is the first report of a single global parameter change that can reliably reproduce spatial and spectral activity patterns in patients with AD and is also correlated with cognitive deficits. These findings provide critical insights about potential mechanistic links between abnormal neural oscillations and cellular correlates of impaired neuronal activity in AD.

Biophysical significance of significantly altered neuronal parameters

The parameters that were differentially distributed in AD were the excitatory time constants $\tau _{\textrm{G}}$ and $\tau _{\textrm{e}}$, and inhibitory neural gain $g_{\textrm{ii}}$. Each parameter has a distinct biophysical meaning, and clear implications in AD pathophysiology, as discussed below.

Long-range time constant. The most important differential parameter was the long-range excitatory time constant $\tau _{\textrm{G}}$, which was (1) increased in AD; (2) capable of recapitulating the spectral shift seen in AD patients; (3) the most important feature in classifying AD from controls. Higher $\tau _{\textrm{G}}$ in AD indicates the slowing of long-range brain-wide communication of neural activity, implicating primarily the large layer-specific pyramidal glutamatergic neurons [39]. These pyramidal neurons are known to be selectively vulnerable in AD [40]. This result is also in concordance with a recent study demonstrating long-range axonal connectivity disruption in AD transgenic mice [41]. Increased pyramidal neuron time constants have also been reported in AD mouse models [42].

We also found that $\tau _{\textrm{G}}$ is also associated with global cognitive deficits in AD patients. This implicates impairment in the synaptic processing of long-range excitatory neurons, a potential factor contributing to increased $\tau _{\textrm{G}}$ in AD. Our recent study showed that alpha hyposynchrony is correlated with the degree of global cognitive dysfunction in patients with AD [27]. Another MEG-based study also demonstrated that oscillatory slowing predicts general as well as domain-specific cognitive function in patients with AD [21]. While such associations can be obtained using neuroimaging data directly, here we were able to identify a specific biophysically grounded parameter, $\tau _{\textrm{G}}$, that can potentially explain these biological relationships. Linking biophysical processes to clinical scales has historically been extremely challenging for conventional machine learning approaches due to the mismatch in dimensionality between input features (thousands) and output features (a handful of clinical measures). Increased long-range time constant in AD capable of recapitulating the spectral shifts in AD and correlated with MMSE, therefore, may be the first report of a single biophysical correlate accounting for clinical deficits in patients with AD.

Local excitatory time constant. Interestingly, the local excitatory neural time constant was also significantly different in patients, whereas local inhibitory time constant was not. Local excitatory-inhibitory imbalances in AD have been demonstrated in numerous basic science studies [20, 43]. Our results are therefore broadly consistent with these studies, although with smaller effect sizes compared to the long-range time constant parameter. Among the local parameters, the strongest relationship was found with the local excitatory time constant $\tau _{\textrm{e}}$, which was higher in AD subjects than in controls, consistent with our previous findings. While higher $\tau _{\textrm{e}}$ implies the slowing of excitatory signals at the local level, we previously demonstrated that increased $\tau _{\textrm{e}}$ is distinctly associated with tau accumulations in AD [20]. The relationship between spatially invariant long-range excitatory time constant and regional tau accumulation in AD remains to be elucidated.

Local inhibitory neural gain. We found that the inhibitory neural gain $g_{\textrm{ii}}$ is higher in patients with AD than in controls. Higher $g_{\textrm{ii}}$ implies a higher synaptic strength of the inhibitory neuronal connections within a region. Alteration in the neural gain term indicates a neuronal excitatory-inhibitory imbalance; such an imbalance has been reported in various preclinical AD models [44,45,46]. However, the current finding is at odds with some prior findings in AD. Previously, we tested the hypothesis that local changes in model parameters could recapitulate regional spectral shifts in AD patients. Using a local circuit model fitting procedure (as compared to the current global network model) our group had previously reported that the local subpopulation in AD showed reductions in inhibitory gains [20]. It may be considered counterintuitive that patients with AD would have stronger inhibitory neuronal connections when interpreting $g_{\textrm{ii}}$ as a measure of hyperexcitability. This discrepancy warrants further study in the future, but we may offer some possibilities. The inhibitory gain $g_{\textrm{ii}}$ interacts with other model parameters of SGM in a complex manner perhaps even more than what is seen in the local decoupled model. While hyper-excitability is consistently seen in AD, its relationship to $g_{\textrm{ii}}$ remains an open question. SGM can achieve hyperexcitability even if $g_{\textrm{ii}}$ is high, depending on the other model parameters. It is further possible that a few specific regions experience reduced inhibitory gain while other, more widespread areas experience a net increase thereof. Nevertheless, it should be noted that the effect size of $g_{\textrm{ii}}$ was only moderate in comparison to $\tau _{\textrm{G}}$’s effect size in the current study.

While the SGM parameters were derived only from neuroimaging data, the range of these parameters are similar to the estimates from other means. The bounds for the time constants, speed, and long-range connectivity coupling constant are explained in detail in [35] and below. The time constants in SGM correspond to the average neural ensemble response functions that capture delays not just due to membrane capacitance or dendritic arborization, but also due to all the local circuit delays. Thus, we expect the time constants to be of the order of tens of milliseconds. Neural responses in the ferret V1 were reported 20 ms after a short virtual stimulus [47]. Moreover, cortical depolarization evoked by a brief deflection of a single barrel whisker in the mouse was reported to spread to parts of sensorimotor cortex within tens of milliseconds [48, 49]. Excitatory/inhibitory synaptic time constants derived from cortico-cortical evoked potentials and mathematical modeling were approximately 5 and 7 ms, respectively [50]. The bounds for the neural gain terms and alpha were set to ensure that the model is stable.

Global network effects versus local circuit effects

Two competing hypotheses can potentially account for the spatial and spectral abnormality patterns in AD: spatially variant effects of local circuits, and spatially invariant global network effects which is the focus of this study. Due to the highly specific spatial topography of AD pathology [51], prior literature has broadly focused on the neural correlates of local circuits as the primary means of describing observed electrophysiological data [43]. In a recent study, we tested the hypothesis that local changes in model parameters could recapitulate regional spectral shifts in AD patients. We also reported that alterations in local excitatory and inhibitory parameters are distinctly associated with tau and amyloid-$\beta$ accumulations in AD patients [20]. The current study addresses a very different hypothesis: that observed spatiospectral changes in AD patients may be explained by global changes in the network, as compared to spatially-varying changes in local neural masses. While the two hypotheses on global network versus local circuit effects in AD are not mutually exclusive, our key contribution here is to show that global changes are sufficient to recapitulate the observed spatial and spectral abnormalities in AD. A previous modeling study also found differences in both coupling and local circuits [52]. Even though AD may induce both local and global changes, it is possible that the latter may dominate, as previously noted from a modeling perspective [53,54,55,56,57] and from our results indicating long-range $\tau _{\textrm{G}}$ as the most important parameter. Heterogeneity of the Amyloid-$\beta$ load was previously found to be essential to simulate the slowing of rhythm [58], but we have demonstrated that spatial variations of any kind in our model are not needed to capture the spectral and alpha-band spatial patterns. This certainly leaves room for the possibility that the SGM will be enriched by including the spatial patterns of Amyloid-$\beta$ and tau. Future explorations of the respective contributions of local versus global network changes in AD will be critical.

This local versus global distinction also means that our current results are not directly comparable to prior spatially-variable modeling results. Both the long-range ($\tau _{\textrm{G}}$) and the local ($\tau _{\textrm{e}}$ and $g_{\textrm{ii}}$) parameters in our SGM model showed significant group differences, but we did not reproduce other local changes reported, e.g., in [20]. Nevertheless, $\tau _{\textrm{e}}$ being higher in AD in both the local as well as global study indicate a common underlying mechanism involving excitatory neuronal subpopulations at both local and global level in AD.

Relationship to previous modeling works

Even though no mathematical model can capture the complex brain structure-function relationship completely, many can aid in identifying mechanisms that cannot be inferred with neuroimaging data alone. Indeed, various model-based markers of AD have also been shown in the past. For example, the Virtual Brain Modeling platform has been used to estimate local and global parameters of a neural mass model for fMRI and to subsequently differentiate between AD and controls [52]. While the literature on fMRI studies in AD is vast, comparable depth is lacking in the use of higher frequency data like MEG. In our work, we focus on MEG because it provides us with a high temporal resolution and can give insights into oscillatory signatures, especially the spectral and spatial patterns thereof, that are directly linked to cellular mechanisms. A neural mass modeling approach attributed slowing of alpha in AD using MEG to neuronal hyperactivity, though without directly fitting to the empirical data [59]. Another modeling approach examined different stimulation strategies to preserve functional network integrity in AD and found that stimulating excitatory neurons were the most successful [60]. Another virtual brain simulation approach integrated local field potential simulations with regional amyloid-$\beta$ and tau uptake as empirical features to classify healthy controls, MCI, and AD and obtained an average F1 score of 0.743 [61] – our study reports a higher F1 score of 0.77 for classification of AD from controls with just a few parameters as features of a random forest classifier.

A key difference from prior modeling approaches is that our SGM is a linear model with a small set of biophysically interpretable global parameters. Therefore, it can be obtained in a closed-form solution in the frequency domain, and model parameter inference is more tractable. We employed SGM because prior studies indicate that the emergent macroscopic activity is independent of the microscopic activity of individual neurons [53,54,55,56,57, 62], and is primarily governed by the long-range connections [63,64,65,66]. Indeed, it was already demonstrated that SGM outperforms a Wilson-Cowan neural mass model in fitting the empirical MEG spectra [33]. A recent comparison showed that linear models outperformed non-linear models in predicting resting-state fMRI time series. This was attributed to the linearizing effects of macroscopic neurodynamics and neuroimaging due to spatial and temporal averaging, observation noise, and high dimensionality [67]. Given that the vast majority of computational models involving neural masses involve highly non-linear concepts like multistability, metastability, and other complex dynamics [68,69,70,71,72], it may be questioned whether AD-induced changes in brain macroscopic dynamics can even be reliably measured and robustly inferred. Instead, we expect that while neural activity in AD and health might be highly dynamic and non-linear, its macroscopic spatial and frequency patterns are known to be far more stable across individuals [20, 27, 54, 73, 74]. This is a key motivation for our use of the linear and deterministic SGM, which has demonstrable tractability and only a few free parameters capable of predicting spectral and regional profiles of MEG activity. To our knowledge, this is the first study identifying a parsimonious biophysically interpretable marker of AD and cognitive decline in AD.

Structural network harmonics are responsible for pathology transmission

It was previously shown by our group that the eigendecomposition of the graph Laplacian can be used to describe the spread of pathology as it ramifies within the brain’s anatomic connectivity network. It was demonstrated that only the eigenmodes corresponding to the lowest eigenvalues - named “persistent modes” are involved in AD pathology progression [75]. Since any aberration in local synchrony explored here must arise from the underlying progression of pathology in the AD brain, it is expected that the same or similar eigenmodes responsible for pathology progression may also be involved in aberrant synchrony. The SGM too can be decomposed into a small set of eigenmodes (see Eq. 1). Remarkably, it was recently shown by our group that the lowest few eigenmodes of the SGM capture a large portion of the spatial distribution of alpha-band power [34], and are also important in explaining low-frequency long-range synchrony from fMRI [76]. This striking resemblance of eigenmodes of both pathological and electrophysiological processes establishes a conceptual bridge that has been hitherto unknown.

Alternative neuroimaging modalities for capturing synchrony

Resting-state functional MRI (fMRI), which measures the slow fluctuations of blood oxygenation signal in the brain as a proxy for neural activity [77], is also widely used to identify abnormal synchrony in AD. Leveraging graph-theoretic analyses, many studies now routinely describe the alterations in fMRI measures during the course of AD pathophysiology [51]. However, graph theoretic statistics of resting-state fMRI have shown inconsistent differences between patients with AD and healthy controls [78]. In addition, fMRI is limited in its ability to capture fast temporal scales of neuronal activity [79]. Electrophysiological techniques such as E/MEG address this limitation by capturing temporal activity scales with millisecond precision, though with a lower spatial resolution. Leveraging their temporal resolution, E/MEG can be used to infer the dynamic neural activity directly [80].

Limitations

To determine the shape of the power spectra we used Pearson’s R as the cost function. Future studies should aim at capturing the magnitude as well as selected spectral features. Here we employed the same template structural connectome from HCP for both cohorts, as it allowed us to pinpoint the biophysical alterations solely due to functional alterations. A prior study has also demonstrated that white matter network organization is preserved in AD [81]. However, it will be ideal to obtain individual structural connectomes in all individuals. In addition, we observed that the SGM fits better to spectral and spatial patterns from AD rather than from controls. This may be attributed to the spectral shape of AD – it has a clearer exponential fall-off that is easier to fit to. In comparison, the spectral shape of controls has an additional peak in the beta band superimposed on the exponential fall-off. Lastly, we note that even though no mathematical model can capture the complex brain structure-function relationship completely, many can aid in identifying mechanisms that cannot be inferred with neuroimaging data alone.

Since the most important indicator of AD was a long-range global parameter, a key question arises whether this result is because of potential inherent bias in the model for long-range effects. This bias towards long-range effects will be more prominent in cases where differences in spatial patterns are the most important features to capture with the model. In the case of this paper, the most prominent pattern changes that we focused on were the spectral pattern changes, especially the spectral slowing in AD. Such spectral changes can be impacted by both global and local parameters. Indeed, the local parameters were also significantly different between AD and controls although with a smaller effect size. Also, while we fit our model to match the spatial pattern of the alpha frequency band as well, there were no prominent differences in the spatial patterns between AD and controls. Therefore, it is unlikely that the inherent bias towards global effects would have impacted the model parameter estimation in this case.

Conclusion

This work suggests that a global impairment in the excitatory long-range pyramidal neuronal population is an important indicator of AD, and is also associated with global cognitive decline in patients with AD. Intriguingly, our work is able to recapitulate the spatial and spectral patterns of AD-related functional activity without introducing any spatial heterogeneity; indeed, the SGM model is entirely global and spatially-invariant. This raises the possibility that a global increase in the long-range excitatory time constant might be a sufficient factor underlying observed spatiotemporal alterations of neuronal activity in AD. This modeling approach highlights a parsimonious framework for identifying cellular biomarkers of abnormal electrophysiological oscillations and cognitive deficits in AD, that can aid in guiding future clinical trials.