Abstract
The accumulation of amyloid-β (Aβ) and hyperphosphorylated-tau (hp-tau) are two classical histopathological biomarkers in Alzheimer’s disease (AD). However, their detailed interactions with the electrophysiological changes at the meso- and macroscale are not yet fully understood. We developed a mechanistic multiscale model of AD progression, linking proteinopathy to its effects on neural activity and vice-versa. We integrated a heterodimer model of prion-like protein propagation and a brain network model of Jansen–Rit neural masses derived from human neuroimaging data whose parameters varied due to neurotoxicity. Results showed that changes in inhibition guided the electrophysiological alterations found in AD, and these changes were mainly attributed to Aβ effects. Additionally, we found a causal disconnection between cellular hyperactivity and interregional hypersynchrony contrary to previous beliefs. Finally, we demonstrated that early Aβ and hp-tau depositions’ location determine the spatiotemporal profile of the proteinopathy. The presented model combines the molecular effects of both Aβ and hp-tau together with a mechanistic protein propagation model and network effects within a closed-loop model. This holds the potential to enlighten the interplay between AD mechanisms on various scales, aiming to develop and test novel hypotheses on the contribution of different AD-related variables to the disease evolution.
Significance Statement
This research presents a closed-loop model of Alzheimer’s disease (AD) mechanisms, bridging the gap between protein distribution and neural activity. Contrary to prior beliefs, the study reveals that interregional hypersynchrony and cellular hyperactivity are not directly linked. Notably, the model identifies neural inhibition as a potential causal factor in neurophysiological AD alterations and posits early depositions of Aβ as a determinant of the spatiotemporal profile of proteinopathy. The significance of this mechanistic disease framework lies in its potential to produce insights into AD evolution and to guide novel treatment strategies. It underscores the importance of further experiments and modeling efforts to refine our understanding of AD, offering hope for more effective treatments and personalized care in the fight against dementia.
Introduction
Alzheimer’s disease (AD) is a progressive neurodegenerative disorder widely characterized at several levels of analysis from molecules to cells to the whole brain. The interactions between these levels of analysis remain poorly understood, contributing to the limited treatment options (Breijyeh and Karaman, 2020; van Dyck et al., 2023). A mechanistic understanding of AD is essential for developing novel therapeutic interventions.
The main pathophysiological changes in AD involve the misfolding and accumulation of amyloid-β (Aβ) and a hyperphosphorylated version of the tau protein (hp-tau) in the brain (Masters et al., 2015). These two proteins have different toxic effects, while Aβ generates hyperactivity through the disruption of GABAergic inhibitory synapses (Palop and Mucke, 2010) and the reduction of glutamate reuptake (Zott et al., 2019), hp-tau disrupts the synaptic connectivity of neurons by reducing the number of dendritic spines in pyramidal cells (Merino-Serrais et al., 2013). Aβ tends to aggregate into plaques, while hp-tau aggregates into neurofibrillary tangles (NFTs). Following the prion hypothesis of AD, misfolded Aβ and hp-tau propagate through the brain acting as seeds (i.e., prions) that trigger the misfolding and aggregation of their normal counterparts (Frost and Diamond, 2009; Walker, 2018; Gomez-Gutierrez and Morales, 2020). It has been shown that the hyperactivity produced by Aβ orients the prionic propagation of hp-tau in the brain (Rodriguez et al., 2020) following the Braak stages from the entorhinal cortex to the neocortex, producing a network disruption that has been linked to cognitive decline (Braak and Braak, 1991; Braak et al., 2011; Therriault et al., 2022). The higher the activity levels in a region, the faster the propagation of hp-tau to those regions. Additionally, hyperactivation has also been linked to an enhanced extracellular secretion and deposition of Aβ generating a cyclical phenomenon that reinforces bidirectionally neural hyperactivity and Aβ concentration (Kamenetz et al., 2003; Cirrito et al., 2005; Stargardt et al., 2015; Yamamoto et al., 2015; Busche and Hyman, 2020).
At the whole brain level, electrophysiological recordings have shown a slowing of α frequency peak along the disease evolution (Babiloni et al., 2009; Garcés et al., 2013). Additionally, several studies have reported an early increase in α power and hypersynchrony between parietal regions (Nakamura et al., 2017, 2018), followed by a disruption in those network’s functional connectivity (FC) and a rise in anterior hypersynchrony (López-Sanz et al., 2016; Pusil et al., 2019). These observations suggest a temporal dissociation between anterior and posterior in the electrophysiological changes that occur in AD (Maestú et al., 2021). Finally, the whole network gets disrupted due to the effect of hp-tau accumulation and further generation of NFTs (Palop et al., 2006; Ranasinghe et al., 2014; Ahnaou et al., 2017; Pereira et al., 2019). Whether the hypersynchronization is linked to the previously mentioned cellular hyperactivity remains unknown (Maestú et al., 2021). Here, we hypothesize that hyperactivity causes an increase in interregional synchrony (i.e., FC) due to a higher rate of neural communication that enhances the coordination of distant neuronal groups.
In recent years, significant computational modeling efforts have employed brain network models (BNMs) to investigate neural activity and its relationship to several aspects of the disease (Stefanovski et al., 2019, 2021; Ranasinghe et al., 2022; van Nifterick et al., 2022; Alexandersen et al., 2023). These models reproduce neural activity using sets of differential equations known as neural mass models (NMMs) that are interconnected through empirically derived structural connectivity (SC) networks based on tractography. These models are expected to serve as a tool for getting deeper insights into the mechanisms of AD and to be able to predict its appearance and evolution.
In this study, we extended a previously published multiscale model of AD evolution based on neuroimaging data that includes molecular, cellular, and interregional features of the brain (Alexandersen et al., 2023). Specifically, it integrates a heterodimer model for the generation and propagation of AD-related proteins, neuronal population activity, and interregional synaptic coupling. We aimed to explore the mechanisms that give rise to the observed neurophysiological changes in AD and to evaluate the relative impact of the different proteins involved. To pursue these goals, first, we explored an isolated neural population and the effects produced by changing the parameter candidates affected by AD (i.e., excitation and inhibition-related parameters). Second, we simulated the neurotoxicity model and adjusted its parameters to reproduce the empirically observed phenomena in AD (i.e., frequency slowing, relative α decrease, increase/decrease in excitation, and increase/decrease in FC). Finally, we use the resulting network model to analyze the impact of different biological mechanisms giving rise to AD progression. The mechanistic explanations derived from the model are key to understanding the effects of delivered treatments, fostering the development of new ones, and enhancing early detection methods.
Methods
Data acquisition
Magnetic resonance imaging (MRI) scans were acquired for 20 healthy participants (17 females) with mean age 60.5 (sd 4.17) at the Centre for Cognitive and Computational Neuroscience, Complutense University of Madrid. Diffusion-weighted images (dw-MRI) were acquired with a single-shot echo-planar imaging sequence with the parameters: echo time/repetition time = 96.1/12,000 ms, NEX 3 for increasing the signal-to-noise ratio, slice thickness = 2.4 mm, 128 × 128 matrix, and field of view = 30.7 cm yielding an isotropic voxel of 2.4 mm; 1 image with no diffusion sensitization (i.e., T2-weighted b0 images) and 25 dw-MRI (b = 900 s/mm2).
All participants provided informed consent.
SC
In a preprocessing stage, FSL eddy was used to correct for eddy current distortion. The correction was conducted through the integrated interface in DSI Studio (http://dsi-studio.labsolver.org). Then, the eddy-corrected dw-MRI data were rotated to align with the antero-posterior commisure line. The restricted diffusion was quantified using restricted diffusion imaging (Yeh et al., 2016). The diffusion data were reconstructed using generalized q-sampling imaging (Yeh et al., 2010) with a diffusion sampling length ratio of 1.25. The tensor metrics were calculated using dw-MRI with a b-value lower than 1750 s/mm2. A deterministic fiber tracking algorithm (Yeh et al., 2013) was used with augmented tracking strategies (Yeh, 2020) to improve reproducibility. The whole brain was used as seeding region. The anisotropy threshold was randomly selected. The angular threshold was randomly selected from 15° to 90°. The step size was randomly selected from 0.5 voxels to 1.5 voxels. Tracks with lengths shorter than 15 or longer than 300 mm were discarded. A total of 5 million seeds were placed.
A combination of Desikan–Killiany (Desikan et al., 2006) and ASEG (Fischl et al., 2002) atlases were used as a volume parcellation atlas, and two SC matrices were calculated (Fig. 1D): weights using the count of the connecting tracks, and tract length using the average length of the streamlines connecting two regions. The final SC consisted of a subset of these networks including 40 regions (Table 3) representing the cingulum bundle (Bubb et al., 2018), one of the most prominent white matter structures that interconnect frontal, parietal, temporal, and subcortical regions, and whose impairment has been related to AD and cognitive decline (Yu et al., 2017; Bubb et al., 2018; Chen et al., 2023).
Modeling framework overview. A, The closed-loop neurotoxicity model integrates a model of proteinopathy dynamics and a BNM to simulate neural activity, interacting in a closed-loop: proteinopathy affects neural activity and vice-versa. B, Scheme of the Jansen–Rit (JR) NMM with three interacting subpopulations where triangles represent synaptic contacts between subpopulations, boxes represent the sigmoidal transformation of voltage into firing rate, and colored hexagons represent the transformation of firing rate into voltage. C, Implementation timeline for the closed-loop model. First, a baseline simulation is gathered to get a baseline firing rate. Then, the closed-loop model simulation begins by integrating the proteinopathy with dt = 0.25 (years) and simulating the BNM to update hyperactivity damage (
BNM
SC matrices served as the skeleton for the BNMs implemented in the virtual brain (Sanz-Leon et al., 2015) where regional signals were simulated using JR NMMs (Jansen and Rit, 1995). The JR is a biologically inspired model of a cortical column capable of reproducing α oscillations through a system of second-order coupled differential equations (Table 1 for a description of parameters; Fig. 1B):
JR-BNM parameters used in simulations
This model represents the electrophysiological activity (in voltage) of three subpopulations of neurons: pyramidal neurons (y0), excitatory interneurons (y1), and inhibitory interneurons (y2) and their derivatives (y3, y4, and y5, respectively). These subpopulations are interconnected through an average number of synaptic contacts (Cep, Cip, Cpe, and Cpi), and integrate external inputs from other cortical columns. The communication is implemented in terms of firing rate, using a sigmoidal function for the conversion of the voltage inputs (Eq. 7). Parameters regarding postsynaptic potential amplitudes (He, Hi) and time constants (
The input (Ii(t)) represents two main drivers of activity in the NMMs: interregional communication and an intrinsic noisy input. The latter is defined as a local and independent Gaussian noise
A global coupling factor g is implemented to scale linearly tracts’ weights. Both g and s are scaling factors that apply to all nodes and were adjusted to reproduce the qualitative features of the AD evolution.
Closed-loop neurotoxicity model
We have adapted an already published model for AD protein propagation, first shown by Thompson et al. (2020), that reproduces the evolution of tau and Aβ concentrations in the AD brain, and a further extension by Alexandersen et al. (2023) that linked this evolution to changes in neural activity patterns using NMMs. Here, we adapt this latter model by implementing a JR NMM and introducing a feedback mechanism that connects the NMMs’ activity to the production and diffusion of toxic proteins in the AD brain (Fig. 1A).
Proteinopathy dynamics are described by the heterodimer model, one of the most common hypotheses that describe the prion-like spreading of toxic proteins. This hypothesis suggests that a healthy (properly folded) protein misfolds when it interacts with a toxic version of itself (misfolded; the prion/seed) following the latter’s structure as a template (Garzón et al., 2021). Therefore, this model includes healthy and toxic versions of amyloid-β (Aβ,
Parameters used in the closed-loop neurotoxicity model
The impact of protein concentration on the BNM is established through two sets of transfer functions. First, a set of damage functions that are used to translate the concentration of toxic proteins into damage variables (Alexandersen et al., 2023):
Second, a set of equations that translate the damage just mentioned into NMM parameter changes. We adapted these equations to the JR NMM. It was assumed that each protein (
∙
∙ Soluble hp-tau (
The effects of
To this point, we have completed half of the closed-loop model. The second half regards the impact of the simulated neural activity on protein propagation. Neural hyperactivity fosters extracellular secretion and deposition of Aβ generating a cyclical phenomenon that reinforces bidirectionally neural hyperactivity and Aβ concentration (Kamenetz et al., 2003; Cirrito et al., 2005; Stargardt et al., 2015; Yamamoto et al., 2015; Busche and Hyman, 2020). Additionally, hyperactivity orients the prionic propagation of
To capture these effects in our model, we introduced an additional damage variable for hyperactivity (Eq. 21), based on the averaged increase in firing rate per region compared to baseline:
In-silico experiments and simulations
In this work, we performed simulations for several complementary in-silico experiments. First, we explored the parameter spaces of a single JR NMM to evaluate the effects of proteinopathy on neural activity, by varying the AD parameter candidates (He, Cip, Cep, and wij). As this first experiment’s simulations were performed with a single node, we accounted for the reduction of interregional connectivity (wij) produced by hp-tau using the mean of the intrinsic noisy input to the node (p), a JR variable that substitutes the interregional afferences not being considered in the simulation. In total, we performed four sets of simulations in this in-silico experiment: the first one focused on the parameters related to Aβ effects, and the remaining three focused on parameters related to tau effects. The simulations were performed for 20 s of neural activity discarding the initial 12 to avoid transients, and we recorded information regarding the spectral frequency peak, spectral power in different frequency bands, and firing rate.
Second, we built and simulated the neurotoxicity model to deepen the mechanisms of AD progression. The initial conditions of the proteinopathy dynamics were fixed using a toxic protein seeding following literature on AD staging: the entorhinal cortex for tau (Braak and Braak, 1991) with an initial
We used the neurotoxicity model in several ways: first, we adjusted the model by exploring the impact of the limits of change for the parameter candidates
The limits of change for the model parameters determine the extent of the proteinopathy’s impact on JR parameters and, therefore, on the neural activity. We explored those parameter spaces to select a set of values that could reproduce qualitatively the empirical observations reported in the AD continuum: the slowing of α frequency, and the rise and decay of relative α power, cellular firing rate, and FC. The adjusted parameters are used here as the default values for our model and are included in Tables 1 and 2. Second, we described the outputs of the model with the default parameters including an averaged description of parameter trajectories (Eqs. 16–19), spectra, firing rate, and PLV of the simulated brain regions; third, we disentangled the contribution of Aβ and hp-tau to the changes in inhibition by isolating their respective effects on the model evolution. This was achieved by setting to zero alternatively the parameters that define the impact of Aβ and hp-tau on inhibition in our model (i.e.,
Spatiotemporal profile assessment
We evaluated the spatial propagation of hp-tau using the Braak stages as a reference and the impact of Aβ seeding over the temporal antero-posterior differentiation in terms of firing rate and FC. Braak stages were adapted from Therriault et al. (2022) considering the regions included in our cingulum bundle network (Table 3): stage one (rI) included the entorhinal cortex; stage two (rII) included the hippocampus; stage three (rIII) included the parahippocampus and amygdala; stage four (rIV) included insula, posterior cingulate, inferior temporal, and inferior parietal; stage five (rV) included other cortical regions.
Regions of the cingulum bundle included in the BNM (Bubb et al., 2018)
We simulated our closed-loop neurotoxicity model by implementing different seeding strategies. The baseline seeding (i.e., fixed strategy) was the same used in previous studies (Thompson et al., 2020; Alexandersen et al., 2023) as stated above. The remaining strategies consisted in randomizing the seeding of Aβ, randomizing the seeding of hp-tau, and randomizing both. The randomization respected two conditions: the number of seeded regions had to be the same as in the fixed strategy and the seeding had to be symmetrical between hemispheres. For the antero-posterior differentiation, we also limited the randomization of Aβ seeding to posterior areas or anterior areas.
One hundred simulations were performed per seeding strategy. For each timestep, we averaged the concentration of hp-tau in the regions of the Braak stages to determine what Braak regions were presenting higher levels of the protein, and we sorted them into Braak sequences following Putra et al. (2021). Then, we extracted the dominant Braak sequence per simulation as the most repeated sequence in time. For the antero-posterior differentiation, we averaged posterior and anterior neural dynamics and we extracted the time to peak (given the rise and decay in the metrics) of each average.
Code accessibility
The models described above were implemented in Python 3.9 and executed in a Surface Pro 7 (8Gb RAM) running Windows 11 and a high performance computing cluster. The code used in the paper is freely available online at github.com/jescab01/ADprogress_2023 and as Extended Data in this manuscript.
Extended Data 1
Code developed in Python 3.9 to simulate and visualize the closed-loop neurotoxicity model. Download Extended Data 1, ZIP file.
Results
The proposed model for the physiological changes in AD is based on extensive literature that describes the effects of Aβ and hp-tau on neural tissue. These changes include the disruption of glutamate reuptake, the reduction of GABAergic synapses, and the reduction of pyramidal dendritic spines. To effectively model these changes, we proposed four parameter candidates that could account for the observed effects. Specifically, we selected two parameters to model the effects of Aβ on neural activity including the amplitude of the excitatory postsynaptic potential (He) to account for the disruption of glutamate reuptake, and the number of local synaptic contacts from inhibitory interneurons to pyramidal neurons (Cip) to account for the reduction of GABAergic synapses. Regarding the effects of tau, we selected three parameter candidates that could effectively model the disruption of pyramidal dendritic spines: Cip mentioned above, Cep as the local synaptic contacts from excitatory interneurons to pyramidal neurons, and wij as the interregional SC weights.
Single-node experiments favor inhibition over excitation to explain hyperactivity
In order to gain a better understanding of the impact of toxic proteins on neural behavior, we explored the parameter spaces of the proposed candidates (He, Cip, Cep, and wij) through single-node NMM simulations. Note that we used the mean of the intrinsic noisy input to the node (p) to account for the reduction of wij as simulations were performed with a single node.
In general terms, and in line with previous analysis of the JR model (Grimbert and Faugeras, 2006; Spiegler et al., 2010), results showed two main patterns of oscillation: a fast limit cycle in the α frequency band (8–12 Hz) and a slow limit cycle in θ frequency band (4–8 Hz). These two behaviors were categorically separated in the parameter spaces as can be observed in the frequency and relative power charts of Figure 2. The noisy regions surrounding the main limit cycles correspond to the fixed point states of JR that are found at both extremes of the excitation–inhibition continuum. In these states, the model might not have enough activation to give rise to an oscillatory behavior, or it might be saturated by a too-high excitation. In both states, the noise acts as a driver giving rise to low-amplitude and low-frequency noisy oscillations, and the more stable the fixed point (i.e., more extreme excitation–inhibition balance) the lower the frequency and amplitude. Accordingly, note that the noisy regions coincide with the low-power regions in the absolute power plots of Figure 2.
Single-node in-silico experiment. The behavior of a single NMM was simulated varying parameter values. Four different experiments are presented in columns exploring different combinations of parameters. In rows, measures derived from the simulations: oscillatory frequency peak of the node, mean firing rate, absolute power, and the relative power in two frequency bands: α and θ. The red circle points out the default parameter values. A 2D representation of each parameter change is shown in Extended Data Figure 2-1. Sample time series for the different oscillatory regimes are shown in Extended Data Figure 2-2.
Figure 2-1
2D representations of single node experiments varying only one parameter at a time (i.e., He, Cip, Cep, p). Dashed lines showing the default parameter value. Markers' colour scales are the same as in the heatmaps of Figure 2. Download Figure 2-1, TIF file.
Figure 2-2
Timeseries and spectra for three different regimes shown in Figure 2: slow limit cycle, fast limit cycle and fixed point state. Two traces from independent simulations are shown. Download Figure 2-2, TIF file.
Focusing on the limit cycle regions, the parameter candidates related to intranode excitation (i.e., He and Cep; see first and second–third columns of Fig. 2, respectively) showed quadratic relationships in their parameter spaces, producing curved isopleth borders separating the slow and fast oscillatory behaviors (i.e., θ and α, respectively). These curves represented a convex relationship between the parameters and the neural behavior in terms of frequency, power, and firing rate. The position of their original values derived from biology (Jansen and Rit, 1995) differed with respect to the vertex of the quadratic curve: the original He was found at lower values than the vertex, while Cep was found at higher values.
For the other two parameters (i.e., Cip and p), the observed changes inside the limit cycle regions were monotonical, and better described as logistic relationships in which the changes happen quickly at the edge between behaviors (Extended Data Fig. 2-1). These relationships were direct in the case of p with both frequency peak and firing rate, in contrast to the inverse relationship found for Cip (i.e., lower inhibition raised the firing rate and the frequency of oscillation). Note that lowering inhibition out of the limit cycle regions had two different effects on firing rate and spectra. The firing rate kept rising monotonically as inhibition lowered, while the frequency peak and power tended to be reduced getting into the noisy low amplitude and slow region commented above.
In this model, two parameters could explain the hyperactivation produced by Aβ: Cip lowering down due to the reduction of GABAergic synapses and He rising due to the disruption of glutamate reuptake. Assuming that the original values of these parameters, which were based on biological observations (Jansen and Rit, 1995), are a good approximation of the underlying ground truth, our experiment would favor a hypothesis in which the hyperactivity found in AD is explained more likely by the reduction of Cip than by the increase of He given that the latter would induce an initial reduction of activation levels that is not expected in AD (see the effects of lowering Cip and rising He in Fig. 2, second row, first column).
Closed-loop neurotoxicity model: from proteins to neural activity
To explore the link between proteinopathy and neural activity, we extended a previously published multiscale model of AD evolution (Thompson et al., 2020; Alexandersen et al., 2023) by implementing a biologically plausible NMM (Jansen and Rit, 1995) and establishing a bidirectional interaction between neural activity (Eqs. 1–6) and proteinopathy (Eqs. 9–12). The proteinopathy starts from an initial distribution of toxic seeds, and during the temporal evolution of the model, the BNM integrates the effects of those proteins over neural tissue, which in turn influences the production of Aβ and the propagation of hp-tau based on firing rate.
The regional evolution of the proteinopathy follows a sigmoidal shape with a final slight decay in Aβ toxic concentration (see average curves in Fig. 3A). The model captures the time delays between the rise of Aβ and hp-tau reported in the literature (Therriault et al., 2022). The effect of proteinopathy on neural activity was modeled through changes in the NMM parameters mediated by a damage variable (Fig. 3B). These changes were monotonical and followed the sigmoidal shape of the proteinopathy and its damage. The impact of Aβ over the neural tissue was modeled as a reduction of Cip and an increase in He, while the impact of hp-tau was modeled as reductions of the local and interregional connectivity parameters explored before (Cip, Cep, wij) as shown in Figure 3C and Extended Data Figure 3-1. The changes in parameters due to the proteinopathy generated a rise and decay in firing rate over time (Fig. 3D) that affected the production of Aβ and the distribution of hp-tau through a hyperactivity damage variable (Fig. 3E).
Evolution of the parameters, state variables, and neural activity outputs from the closed-loop neurotoxicity model. Results are averaged over regions included in eight categories including frontal, parietal, temporal, and cingulate regions. A, Concentration of toxic proteins over time during model execution. B, Damage variables that determine the impact of the proteinopathy over the neural activity (i.e., JR-BNM parameters). C, JR-BNM parameters change over time. D, Functional measures derived from the simulation of the BNM including spectral frequency peak and power, PLV, and firing rate. The latter is used in the calculation of q(ha). E, Damage variable q(ha) that determines the impact of hyperactivity on the production and propagation of toxic proteins. The average of the parameter changes shown in (C) are plotted over the single-node experiments’ heatmaps in Figure 3-1.
Figure 3-1
Parameter trajectories of the simulated closed-loop model with default parameters on heatmaps from single node experiments. Trajectories are represented as curves with varying colours. The colour represents time: starting in red and ending with blue. Note that heatmaps are extracted from single-node simulations in which just two parameters are varied at a time, while the trajectories imply a 4-dimensional parameter change. Therefore, the underlying heatmaps should be interpreted just as an orientation of what might happen when changing parameters in one direction. Download Figure 3-1, TIF file.
The specific limits of change of the parameters (included in Table 2) were fitted by iterative explorations of the parameter spaces generated with the neurotoxicity model to achieve the reproduction of some of the main pathophysiological changes of AD evolution (Maestú et al., 2021): frequency slowing, reduction in relative α power, hypo-/hyperactivity, and hypo-/hypersynchrony in α (Fig. 4). We noticed that the reduction of α FC does not translate into an enhancement of FC in lower frequency bands (i.e., θ/δ; Extended Data Fig. 4-1). Additionally, we could not reproduce a clear temporal distinction between antero-posterior effects of the proteinopathy on FC (i.e., posterior regions are affected first) as it has been previously proposed (Fig. 4A, PLV). However, we also found slight antero-posterior differences in magnitude between those regions for FC and firing rate.
Functional measures from the closed-loop neurotoxicity model. Metrics averaged from the simulation shown in Figure 3. A, We extracted spectral frequency peak, relative band power, firing rate, and α FC from the simulation over 40 years of the neurotoxicity model. A quadratic relationship is found over time for all these metrics (i.e., measures increase early in time and decrease later). B, Temporal excerpts of the same simulation including 3D representations of the BNM with node size as degree, node color as firing rate (orange/red is a higher rate), and edge colors as PLV. Last row shows averaged normalized spectra for all the regions included in the model. FC evolution in other frequency bands is shown in Figure 4-1. Regional spectra at several time points are shown in Figure 4-2.
Figure 4-1
Averaged FC (PLV) of the simulated closed-loop model with default parameters for different frequency bands: delta (2 - 4 Hz), theta (4 - 8 Hz) and alpha (8 - 12 Hz). Download Figure 4-1, TIF file.
Figure 4-2
Normalized spectra per region over time (yrs.). Note the transition of all spectra to more noisy states with lower frequency peaks. Download Figure 4-2, TIF file.
Hyperactivity is not directly linked to hypersynchrony
The parameter space explorations that allowed us to fit the model are used in the following to establish relationships between different variables and outcomes derived from the neurotoxic model simulation. To obtain these parameter space explorations, we simulated the model with different values of the parameter candidates’ limits of change.
The model behavior could be analyzed by dividing it into two stages: a first stage in which we find an increase in relative α power, firing rate, and FC; and a second stage in which all these measures decay (Fig. 5). These two stages are not synchronized for all measures. We observed a delayed peak in firing rate at the time interval of years 18–25, while all other measures peak earlier in time between years 10 and 18. From year 18 to 25, the firing rate keeps rising while power and FC have already lowered. This observation suggests a dissociation between firing rate and FC that has not been previously established.
Parameter spaces for the limits of change of the JR-BNM parameter candidates. For each parameter candidate (rows) and limit value explored (heatmaps y-axis), one simulation of 40 years (heatmaps x-axis) is performed. For each simulation, four neural activity measures are extracted including averaged frequency peak, relative α power, firing rate, and FC (columns). When exploring a parameter candidate, all other parameters are kept fixed with their default values that are highlighted in red and included in Table 2. These default values were used in the simulations for Figures 3 and 4. Note how taking apart from zero the limits of change impact differently the resulting behavior of the model. Two additional parameter spaces were computed to define values for g and s (Fig. 5-1).
Figure 5-1
Parameter spaces to select working point. Selected parameters were g=25 and s=20 m/s (see red highlights). Note how the rising of g leads to a situation in which no FC rise is observed, similar to the reduction of s, due to high early FC levels. Also, lowering g leads to the prebifurcation regime of the JR NMMs, a situation in which the spectral frequency peak lowers towards the delta band. Download Figure 5-1, TIF file.
The changes in inhibition lead the evolution of AD
Two parameter candidates could be responsible for the initial rising stage: the increase of He and the reduction of Cip. Looking at the parameter spaces, only the changes to the Cip limits showed a disruption in the rising phase when reduced less than a 20% (Fig. 5, second row). On the other hand, changes to the limits of He did not show any significant effect in this stage, even using an upper limit equal to 0% (i.e., keeping the parameter fixed in time) the rise in all measures could still be observed (Fig. 5, first row). This would suggest that the rising stage is mainly controlled in our model by the changes in inhibition (i.e., the reduction of Cip parameter) representing the Aβ dependent reduction of GABAergic synapses.
The decaying stage required a more complex analysis. In the previously exposed single-node experiments, we observed that lowering inhibition showed a transition from the slow limit cycle, to the fast limit cycle, to the noisy low frequency and amplitude oscillation in which excitation–inhibition imbalance saturates the regional response and leads the JR NMMs to operate in a fixed point state (Fig. 2 and further explanations above). In that context, the reduction of Cip itself could explain the decay stage for frequency peak and relative α power by entering the noisy region where lower inhibition (i.e., further imbalance of excitation–inhibition) leads to a more stable JR fixed point state and lower spectral frequencies. Interestingly, those spectral changes were time-locked with the changes in α band FC (Figs. 4 and 5). Despite explaining the decay stage for those three measures, the reduction of Cip could not explain the decay for firing rate. Lowering inhibition produced a monotonic ascent as could be expected from Figure 2 even in the noisy low-amplitude region of the single-node parameter spaces.
The decay of firing rate could then be explained by another two parameter candidates related to the tau-derived reduction of pyramidal dendritic spines: the reduction of Cep and the reduction of SC weights. The parameter spaces showed that only the reduction of local excitation (i.e., Cep) interacted with the firing rate (Fig. 5, third column). Disabling the change for Cep (i.e., Cep change limit of 0%) showed a plateau in firing rate during the decay stage, while activating its reduction showed a progressive effect on the decay of the firing rate (Fig. 5, third row, compare
Amyloid-β controls changes in inhibition
From the analysis above, we can derive that the disruption of inhibition is a relevant driver of AD progress. Given that both Aβ and hp-tau contribute to the disruption of inhibition by reducing GABAergic synapses and reducing dendritic spines, respectively, we wondered whether they contribute differently to the inhibitory-related effects reported above.
To disentangle their contribution, we computed an additional set of parameter spaces in which we evaluated independently the impact of the inhibitory disruption caused by each protein (Fig. 6). We observed a clear similarity between the parameter space of Aβ-only effects and the previously shown full model in Figure 5. The isolated effect of tau could also generate a slight rise and decay in power, frequency, and FC delayed from the Aβ effect and shorter in time, with a more moderate change in firing rate. This would suggest that both proteins could produce to some extent the effects observed in AD due to the disruption of inhibition, although with different timings and hyperactivity levels.
Parameter spaces for Cip isolating the effects of Aβ and hp-tau. First row shows the results of simulating the neurotoxicity model with different limits of change for Cip isolating the effects of Aβ on inhibition, while the second row shows the results of simulating the same parameter space isolating the effects of hp-tau on inhibition.
Toxic seeding determines the spatiotemporal profile of the evolution
The spatiotemporal sequence of hp-tau spreading in AD has been accurately characterized through the Braak stages. We used this characteristic spatiotemporal sequence to evaluate the spatial accuracy of the tauopathy simulated in our model, and the impact of the selected toxic proteins’ seeding regions on its spatial evolution. We implemented four seeding strategies including fixed seeding, Aβ randomized, hp-tau randomized, and both randomized. We performed 100 simulations per strategy and we measured the concentration of hp-tau in the regions of the Braak stages extracting the dominant spatial sequence of propagation per simulation.
Results showed that fixed (i.e., as usual) and Aβ random seeding strategies reproduced effectively the Braak staging in 100% and 70% of the simulations, respectively (Fig. 7A). In these simulations, just a slight temporal shift separated stages II (i.e., hippocampus) and III (i.e., parahippocampus and amygdala), and this way 30% of the simulations with Aβ-random seeding showed a sequence in which hp-tau reached stage III before stage II. Interestingly, randomizing tau seeding produced a variety of Braak sequences that lowered the dominance of the theoretical sequence to less than 5% of the simulations, at the same levels of randomizing both seeding locations. Examples of these sequences in time are shown in Figure 7B. These observations posit the importance of the localization of hp-tau seeds (i.e., early accumulations) to determine the evolution of its propagation.
Dominance of Braak sequences per seeding strategy. A, The percentage of simulations per seeding mode in which different Braak sequences of hp-tau propagation were dominant. Fixed implied seeding a usual, Aβ random implied randomizing Aβ seeding and using the fixed version for hp-tau, vice-versa for tau random, and random implied randomizing both Aβ and hp-tau seeding. B, Samples of the most representative Braak sequences showing the normalized concentration of
In addition, we wondered to what extent the antero-posterior temporal differentiation of the neurophysiological changes observed along the AD continuum could be related to the seeding of Aβ and hp-tau. We performed an additional set of simulations randomizing the seeding of Aβ and hp-tau and measuring the differences in firing rate and FC between anterior and posterior regions.
We observed a stable precession of posterior regions in terms of firing rate for the cases of fixed seeding (i.e., following Palmqvist et al., 2017), randomizing Aβ and hp-tau in posterior regions, and randomization them over all regions (Fig. 8). Only in the case of anterior seeding, the model showed a slightly faster peak in anterior regions. This suggests that hyperactivity tends to appear earlier in posterior regions disregarding where Aβ and hp-tau start to accumulate. Note that Aβ and hp-tau seedings affected differently the levels of hyperactivity shown by the regions: non-Aβ-seeded and tau-seeded regions lowered their maximum firing rates as shown in Extended Data Figure 8-1.
Impact of Aβ and hp-tau seeding on the antero-posterior temporal differentiation for firing rate (top) and FC (PLV; bottom). The neurotoxicity model shows a rise and decay for several variables including firing rate and FC. We measured the time elapsed to that peak for anterior and posterior regions as a proxy to understand where the neurophysiological changes happened before. Samples of the temporal evolution for each seeding strategy and for both firing rate and FC are represented in Figure 8-1.
Figure 8-1
Samples of the evolution for firing rate (left column) and FC (i.e., PLV; right column) in the antero-posterior differentiation experiments. In rows, each of the four seeding implemented strategies. The effect on FC is limited to a temporal shift of the curves, however, the seeding affects the level of hyperactivity reached by the anterior or posterior regions. Download Figure 8-1, TIF file.
In contrast, the antero-posterior temporal differentiation for FC relied more on the seeding of Aβ. When the seeding was located in posterior regions they peaked earlier while, on the contrary, when anterior regions were seeded they tended to peak earlier. Interestingly, these two strategies generated a general delay in the time to peak both in the seeded region and in the complementary. Note that changing hp-tau seeding reduced the time to peak but did not affect the antero-posterior differentiation relative to the fixed strategy.
Discussion
In this study, we developed a multiscale closed-loop neurotoxicity model of AD integrating a BNM that reproduced brain activity with another modeling the production and propagation of toxic proteins in the brain, both influencing each other through biologically plausible links. Our results showed that the disruption of inhibition was favored to explain cellular hyperactivity, and indeed that inhibition was key to explaining the AD-related changes in frequency, power, and FC. These changes in inhibition were primarily attributed to Aβ effects. Additionally, our results suggested a disconnection between cellular hyperactivity and interregional hypersynchrony. Finally, we observed a significant association between the spatiotemporal profile of the evolution of AD to the initial localization of the accumulation of AD proteins.
Single-node experiments favored the lowering of inhibition (related to GABAergic synapses’ disruption) over the rising of excitation (related to impaired glutamate reuptake) to explain the characteristic hyperactivity produced by toxic Aβ on neural tissue. This idea is in line with previous animal studies showing a direct association between the levels of GABAergic synaptic activity and the hyperactivity found in the vicinity of amyloid plaques (Busche et al., 2008, 2015). The conclusions derived from the single-node experiments are based on the assumption that the original definition of the JR parameters approximates well the characteristics of neural populations. Another definition of parameters could lead to different results. However, note that JR parameters were derived from empirical studies, for instance, Cip and Cep were suggested in light of histological research establishing proportions of synaptic connections between different cell types in several neocortical regions (Jansen and Rit, 1995). In any case, these experiments propose testable predictions regarding the extent to which these two Aβ-derived physiopathological phenomena may contribute to neural activation in AD.
Additionally, single-node experiments showed different effects for the changes implemented to the parameter candidates in terms of neural activity measures. Inside the limit cycle space, we observed quadratic relationships with firing rate and power for the local excitation parameters (i.e., He and Cep), while a monotonical relationship was found for the inhibitory parameter (i.e., Cip; Fig. 2). These observations suggest two different roles for excitation and inhibition. While inhibition would predispose the system to order neural activation into synchronized oscillations, excitation is needed to trigger this mechanism. Therefore, both increasing and decreasing excitation levels could lead to a rise in cellular firing rate: while low excitation may become ineffective to trigger inhibition, high excitation could overcome the capacity of inhibition to organize action potentials and generate oscillatory behavior. This type of contradictory effect (i.e., rising excitation parameter could lower firing rate, and vice-versa) has been reported in empirical studies under the name of inhibitory stabilization of network dynamics, in which, for instance, the stimulation of inhibitory neurons produce a counterintuitive suppression of inhibitory firing (Sadeh and Clopath, 2020; Sanzeni et al., 2020). This mechanism is characteristic of networks with strong recurrent connectivity and protects the system from unstable dynamics such as epileptic seizures.
In agreement with single-node experiments, the simulation of the multiscale neurotoxicity model gave a prominent role to inhibition leading the neural activity changes typically observed in AD. Despite both Aβ and hp-tau can be associated with changes in inhibition—due to the disruption of GABAergic synapses, and to the reduction of dendritic spines that reduces the number of inhibitory inputs to pyramidal cells, respectively—in the model, Aβ was linked to the main inhibitory changes that lead AD evolution. However, it was also shown that hp-tau alone could generate inhibitory effects in the same direction but weaker, shorter, and delayed in time.
The relationship between cellular hyperactivity and interregional hypersynchrony in AD has been suggested (Sepulcre et al., 2017; Koelewijn et al., 2019; Maestú et al., 2021; Stam et al., 2023) but not directly tested to our knowledge. The rationale behind this association is based on the studies linking Aβ independently with both hyperactivity and hypersynchrony. On one hand, some studies report higher cellular activation, or lower firing thresholds, in regions affected by Aβ (Busche et al., 2008; Palop and Mucke, 2010). On the other hand, some studies report spatial correlations between Aβ accumulation (and hp-tau later) and the levels of FC, which are mainly observed within the default mode network (Mormino et al., 2011; Schultz et al., 2017; Quevenco et al., 2020; Ranasinghe et al., 2020; Schoonhoven et al., 2023). From these observations, it has been derived that Aβ produces cellular hyperactivity, and that cellular hyperactivity, in turn, produces FC hypersynchrony. In our study, we show that hyperactivity and hypersynchrony might not be directly linked, as we found a reduction of FC while Aβ accumulation was still enhancing cellular hyperactivity, suggesting a temporal mismatch between the two. In contrast, FC levels were strongly associated in time with changes in spectral power. We believe that this dissociation between spectral power and firing rate is essential to understand the hyperactivity/hyperconnectivity question. Whenever hyperactivity overcomes the capacity of inhibition to order neural firing into an oscillatory behavior, the regional activity may become noisier and therefore more difficult to synchronize interregionally. However, when hyperactivity moves the excitation/inhibition ratio toward a more balanced state, larger oscillations appear that may give rise to FC changes in a network.
The simulation of the closed-loop neurotoxicity model showed a sigmoidal shape in the evolution of the concentration of toxic proteins with a final slight decay for
In this work, we used the entorhinal cortex as the seeding region for hp-tau propagation. However, recent studies suggest that the classical Braak staging of tau pathology in AD (Braak and Braak, 1991) may not offer a comprehensive view and may disregard alternative tau epicenters and propagation trajectories (Franzmeier et al., 2020; Lamontagne-Kam et al., 2023). From the epicenters, hp-tau spreads to functionally connected regions (Franzmeier et al., 2022; Schoonhoven et al., 2023) producing different network disruptions and clinical profiles. In a study by Vogel et al. (2021), they identified four distinct tau propagation patterns, each capable of explaining between 18 and 33% of the cases. We believe that our framework poses a good opportunity to advance the understanding and prediction of personalized tau propagation trajectories. Our approach considers not only structural and FC networks but also accounts for the progressive changes in neural activity due to proteinopathy.
An aspect that has not been captured by our model is the effect of chronic neuroinflammation. It may affect the normal functioning of glial cells, disrupting the clearance function of both microglia and astrocytes and, therefore, fostering protein accumulation and neuroinflammation (Heneka et al., 2015; Henstridge et al., 2019). Although it is not the only additional aspect that could be implemented in the model due to its impact on AD progression (e.g., vascular damage leading to reduced protein clearance, infections such as herpes virus, genetic factors, etc.), we expect that by including these glial effects, we could reproduce the slight decay in hp-tau concentration that is found in the literature, but not reproduced in our results. Further modeling studies including additional AD-related factors are needed.
Regarding simulated neural activity, our model qualitatively reproduced some of the most relevant biomarkers along the AD continuum, including the slowing of the α frequency, the rise and decay of relative α power, firing rate, and α FC, and the Braak stages (Maestú et al., 2021). However, we could not directly reproduce two FC-related phenomena reported in the literature: the rise of FC for lower frequency bands (Nakamura et al., 2017; Schoonhoven et al., 2022), and the spatiotemporal dissociation in FC between posterior and anterior brain regions (López-Sanz et al., 2016; Nakamura et al., 2017, 2018; Pusil et al., 2019). The latter was dependent on spatial seeding of Aβ that, by including the orbitofrontal cortex in addition to the precuneus and posterior cingulate cortex (Palmqvist et al., 2017), balanced antero-posterior changes in time. We believe that assuming this spatial seeding to be correct, the result could be expected empirically. However, we wonder whether the technical difficulties in detecting accumulations of Aβ could hinder the detection of earlier accumulations in posterior regions that could explain the electrophysiological evidence mentioned above.
Previous modeling research has explored the mechanisms underlying frequency slowing in AD, including the work by Alexandersen et al. (2023) that served as the basis for the research presented here. They found that the oscillatory changes in the AD continuum were more related to local neurodegeneration than to the disruption of interregional SC. Our results reproduced most of their observations, supporting the local neurodegeneration hypothesis. Additionally, van Nifterick et al. (2022) presented an exploration of six excitatory-inhibitory mechanisms that could lead to the slowing of AD. They found that both hyperexcitation and disinhibition (related to the Aβ neurotoxic effects) could lead to oscillatory slowing. Our model integrated three of these mechanisms: rise of the excitatory postsynaptic potential (He), lowering of the inhibitory coupling (Cip), and lowering of excitatory coupling (wij). Their results fit into the single-node experiments parameter spaces presented here, which offer a comprehensive picture in which non-linear relationships between the levels of excitation–inhibition and the spectral frequency changes can be observed. Finally, other studies have previously modeled the frequency slowing observed in AD by manipulating the time constants of the NMMs. Ranasinghe et al. (2022) fitted the model parameters to match empirical spectra from patients, finding relationships between Aβ accumulation and inhibitory time constants changes, and between hp-tau accumulation and excitatory time constants. Also, Stefanovski et al. (2019) modeled the inhibitory disruption through changes in the inhibitory time constants linked to Aβ accumulation. These models point to inhibitory disruptions due to Aβ accumulation, in line with the results presented here. We did not include effects on time constants in our modeling approach considering that the empirical evidence for AD at the cellular level has reported more often changes in the magnitude of excitation/inhibition instead of changes in neural time constants.
This work provides a comprehensive and self-contained mechanistic explanation of how different factors contribute interactively to the course of AD and, therefore, it becomes relevant for urgent clinical needs. First, the observed hyperexcitation is the treatment target of the approved symptomatic drug Memantine (Reisberg et al., 2003, 2006) and other treatments currently under investigation such as the anti-epileptic drug Levetiracetam (Vossel et al., 2021). Several clinical trials in AD have determined that selecting the correct individuals at the correct disease stage for a particular intervention is crucial, our model aims to help in this process in the future. Furthermore, both increased neural activity and altered synchrony can be targeted with invasive or non-invasive brain stimulation (Chang et al., 2018). Deep brain stimulation in the fornix has positively affected cognition, but mainly if the stimulation site precisely activates memory networks (Ríos et al., 2022). The promising effects of transcranial magnetic stimulation depend on a particular adaptation between the used frequencies and their effects on the excitation–inhibition balance (Weiler et al., 2019). Neuromodulation cannot only improve but also deteriorate the function of a diseased brain, and the prediction of an individual stimulation regime requires a deep understanding of the underlying mechanisms. The same applies to pharmacological treatments: while the recent success of Lecanemab in slowing down cognitive decline is promising and raises hopes for the amyloid hypothesis after decades (van Dyck et al., 2023), the observed effects are minimal, and the long-term outcome is still unclear. One avenue toward more efficient use of potentially disease-modifying treatments is, again, the choice of the correct timing for the therapy. The closed-loop model in this work can even explain several potential reasons that hinder an effective slowing of the neurodegenerative process: it explains how both Aβ and hp-tau can drive the vicious circle independent from each other, and it formalizes how the self-sustaining process evolves even without the presence of the initiating factor. Altogether, only a holistic understanding of AD will open the possibility to once even preventing the development of dementia, and computational neuroscience can powerfully contribute to this.
Footnotes
The authors declare no competing financial interests.
We thank Spanish Ministry of Universities (FPU2019-04251) and European Union eBRAIN-Health Consortium (HORIZON-101058516).
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
