Personalized brain models link cognitive decline progression to underlying synaptic and connectivity degeneration

Abstract

Cognitive decline is a condition affecting almost one sixth of the elder population and is widely regarded as one of the first manifestations of Alzheimer’s disease. Despite the extensive body of knowledge on the condition, there is no clear consensus on the structural defects and neurodegeneration processes determining cognitive decline evolution. Here, we introduce a Brain Network Model (BNM) simulating the effects of neurodegeneration on neural activity during cognitive processing. The model incorporates two key parameters accounting for distinct pathological mechanisms: synaptic degeneration, primarily leading to hyperexcitation, and brain disconnection. Through parameter optimization, we successfully replicated individual electroencephalography (EEG) responses recorded during task execution from 145 participants spanning different stages of cognitive decline. The cohort included healthy controls, patients with subjective cognitive decline (SCD), and those with mild cognitive impairment (MCI) of the Alzheimer type. Through model inversion, we generated personalized BNMs for each participant based on individual EEG recordings. These models revealed distinct network configurations corresponding to the patient’s cognitive condition, with virtual neurodegeneration levels directly proportional to the severity of cognitive decline. Strikingly, the model uncovered a neurodegeneration-driven phase transition leading to two distinct regimes of neural activity underlying task execution. On either side of this phase transition, increasing synaptic degeneration induced changes in neural activity that closely mirrored experimental observations across cognitive decline stages. This enabled the model to directly link synaptic degeneration and hyperexcitation to cognitive decline severity. Furthermore, the model pinpointed posterior cingulum fiber degeneration as the structural driver of this phase transition. Our findings highlight the potential of BNMs to account for the evolution of neural activity across stages of cognitive decline while elucidating the underlying neurodegenerative mechanisms. This approach provides a novel framework for understanding how structural and functional brain alterations contribute to cognitive deterioration along the Alzheimer’s continuum.

Introduction

Cognitive decline is a clinical condition that affects several aspects of everyday life, impairing memory, attention and social functioning [1]. When cognitive decline severely compromises the individual’s ability in performing everyday tasks, it is referred as dementia, typically associated with serious underlying conditions such as Alzheimer’s disease [2].

Several mechanisms of structural alterations have been proposed to contribute to cognitive decline, including brain network disconnection [3], alterations in network topology [4], and in synaptic transmission [5]. However, the precise connections between the decline in cognitive performance and the structural brain network degeneration remain elusive [6].

Brain-inspired neural networks offer a robust framework for understanding the structural bases of cognitive processes [7,8,9,10,11]. Leveraging modern imaging techniques, these models incorporate biologically realistic features to investigate the effects of brain architecture on cognitive output [7, 11, 12]. By systematically manipulating various parameters, such as neural connectivity or synaptic transmission, researchers investigated how functional and cognitive processes depend on the underlying structural organization [7, 13]. However, the use of brain-inspired neural networks has been limited to either investigate cognitive output [7, 12,13,14,15] or replicating features of neural activity [16,17,18].

Despite biological brain networks are subject to structural defects and neurodegeneration processes, these effects are often overlooked by current network models [19,20,21]. While substantial research has focused on the structural bases of optimal cognition, there is a notable gap in studies investigating the interplay between structural network alterations and decline in cognitive function [22]. Understanding these interactions is crucial for elucidating how cognitive processes are affected by defects in network architecture, providing a biologically realistic framework connecting structural degeneration to alterations in neural activity and in cognitive processes [22].

Recent research highlights the potential of electroencephalography (EEG) recordings as a promising tool for investigating neural activity during cognitive processes [23] and monitoring conditions such as cognitive decline [24,25,26,27]. Cognitive tasks engage attentional processes and elicit characteristic event-related potentials (ERPs) that reflect neural responses to stimulus presentation [28]. Cognitive decline is associated with alterations in neural activity and ERPs, notably reductions in peak amplitudes and prolonged latencies of ERPs components [23]. Another measure obtainable from EEG is functional connectivity (FC), which assesses how activity in one brain region can influence other areas [29]. While brain-inspired neural networks have been widely used to mimic the optimal performance of biological neural systems, no study has yet attempted to investigate and replicate neural activity, such as EEG, under conditions of cognitive decline.

We developed a Brain Network Model (BNM) to replicate neural activity during cognitive processing related to task execution (see Fig. 1 and Methods). The model leverages a connectivity architecture derived from human tractography data, combined with a reservoir computing approach [30] to create a digital version of the task performed by participants. By incorporating virtual neurodegeneration processes associated with pathological cognitive impairment [31], the BNM successfully reproduced human EEG responses across various stages of cognitive decline.

Fig. 1

Reproducing individual neural activity during task execution at different levels of cognitive decline (A): EEG recordings were collected during task execution from participants in several stages of cognitive decline. B The BNM simulates EEG activity during virtual task execution (see Methods, “Virtual version of cognitive task” paragraph). C Cognitive task is modelled by a virtual stimulus propagating along the ventral stream of visual information. D, E Both experimental (D) and simulated (E) EEGs were pre-processed to extract ERP components. F Relevant EEG features such as encoding functional connectivity and N1 waveforms were computed from extracted ERP in both experimental and simulated signals. G Experimental and simulated ERPs are compared to tune the parameters of the BNM for each participant, using a gradient descent algorithm. The algorithm determines the combination of model parameters that minimizes difference between experimental and simulated signals. The procedure is repeated for each participant. H The determination of personalized parameter combinations allowed us to reproduce individual EEG signals. Signals were reproduced across several EEG channel groups, specifically occipital channels, prominent in encoding processes (top) and central channels, prominent in decision processes (bottom). I Personalized parameters also allowed to reconstruct a BNM for each participant, investigating structural configurations associated with each cognitive condition. Abbreviations: $cp$: connectivity parameter (parameter of connectivity degeneration); EEG: electroencephalography; ERP: Event-Related Potential; exp: experimental, IT: Inferior Temporal Cortex; $lp$: local parameter (parameter of synaptic degeneration); sim: simulated; V1: Primary Visual Cortex, V2: Secondary Visual Cortex, V4: Visual Area V4

Through model inversion, personalized BNMs were reconstructed for each participant based on individual EEG recordings. Crucially, personalized models adopted different network configurations for patients in each cognitive condition, with neurodegeneration levels proportional to the severity of cognitive decline. Our model revealed a neurodegeneration-driven phase transition in simulated neural activity, predicting the presence of two distinct patterns in neural activity which we later observed in experimental signals. Our work represents the first example of a model capable of replicating individual neural responses to cognitive task execution across different stages of cognitive decline, identifying potential structural culprits of the condition.

Methods

Participants recruitment

Participants were recruited as part of the PREVIEW project [32] (clinical trial identifier: NCT05569083, first submitted 2022–08–24), a longitudinal study on cognitive decline and prodromal stages of Alzheimer’s disease. Detailed information on the study design can be consulted in Mazzeo et al. [32]. We included participants at two stages of cognitive decline: subjective cognitive decline [33] and mild cognitive impairment [34].

Subjective cognitive decline is characterized by self-reported cognitive deterioration, not detectable through standardized cognitive tests [33]. Mild cognitive impairment involves a clinically detectable decline in cognition that does not significantly interferes with daily activities [34]. A total of 126 self-referring patients with cognitive complaints participated in the study, categorized as 86 SCD and 40 MCI. Classification of SCD followed the terminology proposed by the Subjective Cognitive Decline Initiative (SCD- I) Working Group, while MCI classification adhered to the National Institute on Aging-Alzheimer’s Association (NIA-AA) work-groups criteria [34].

Inclusion criteria were: (i) diagnosis of either SCD or MCI; (ii) Mini Mental State Examination (MMSE) score > 24, (including correction due to age and education levels); (iii) normal scores on the Activities of Daily Living (ADL) and the Instrumental Activities of Daily Living (IADL) scales.

Exclusion criteria included: (i) history of severe head trauma; (ii) concomitance of neurological and/or systemic diseases, psychosis or major depression; (iii) history of substance abuse (including alcohol); (iv) use of medications known to affect EEG oscillations (e.g., benzodiazepines, antiepileptic drugs); (v) overt Alzheimer diagnosis according to the criteria specified by National Institute on Aging-Alzheimer’s Association (NIA-AA) [35].

Additionally, we included a control group of 19 age-matched healthy subjects (CTR, age: 64.06 ± 4.93) who volunteered for the study. Patients in the SCD and MCI groups were predominantly females (Supplementary Table 1), following the standard gender distribution of the pathology. There was no significant difference in gender composition across groups (χ² test = 4.73, p = 0.094).

Participant recruitment and EEG recordings were conducted in accordance with the Declaration of Helsinki guidelines and the standards set by the Committee on Human Experimentation of the Careggi University Hospital in Florence, Italy. The study received approval by the local Institutional Review Board (reference 15691oss). Relevant demographic information about participants is provided in the Supplementary Materials (Supplementary Table 1).

Three-choice vigilant task structure

Participants performed a Three-Choice Vigilance Task (3CVT) while undergoing EEG recording. The test was designed to assess two fundamental aspects of cognition: sustained attention and associative memory. In the 3CVT, participants were required to discriminate between a primary geometric shape consisting in an upward triangular target (referred to as TARGET) and two secondary shapes: a diamond and an inverted triangular target (referred to as NON-TARGET). The task consisted in two phases: an encoding phase, where the stimulus was presented to the participant, and a decision phase, where the participant indicated whether the presented shape was a TARGET or a NON-TARGET.

During task, each shape was presented individually for 0.2-s, with random-duration breaks between presentations. The mean number of shape presentations was 733 ± 45. Participants responded using a keypad, pressing the left button for TARGETs and the right button for NON-TARGETs. TARGET stimuli were presented with 70% probability, while NON-TARGET stimuli were presented with 30% probability. A brief training session preceded the test to familiarize participants with the test and minimize errors. The total task duration was of 20 min, and a black screen signalled the end of the task. Test performance was evaluated based on accuracy, reaction time and F-measure, which is an aggregated measure of both accuracy and reaction time.

Experimental EEG device and signal acquisition

EEG data were collected from all participants meeting the study’s inclusion criteria [32]. Recordings were conducted while participants were sitting comfortably during the 3CTV task, using the 64-channels Galileo-NT system (E.B. Neuro S.p.a.). The sensor topography followed the extended 10/20 system [36]. Unipolar signals were recorded at a sampling rate of 512 Hz. Electrode impedances, constantly monitored during EEG acquisition, were maintained in the 7–10 kOhm range. Portions of the signal where impedance fell outside this range were discarded.

EEG pre-processing and ERP computation

The EEG pre-processing pipeline was executed on MATLAB using the commercially-available EEGLAB Toolbox [37]. Pre-processing steps included (i) band-pass filtering in the 1—45 Hz using a Butterworth filter, (ii) identification and removal of noisy channels, (iii) average re-referencing and (iv) Independent Component Analysis (ICA)-based detection and elimination of artefactual components. Each trial of the 3CVT task experiment lasted 1000 ms, divided in two parts: stimulus presentation in the first 200 ms (encoding phase), and participant response in the following 800 ms (decision phase).

Event-related potentials (ERPs) were epoched in the [0 ms,750 ms] window, with 0 ms aligned to stimulus presentation (this choice was made to the align experimental stimulus with the virtual one injected in the model). We analyzed the N1 ERP component [23, 38] in occipital channels (PO7, PO8, O1, O2, Oz). The N1 component, the first negative deflection after stimulus onset, occurs in the encoding timeframe ([0 ms, 200 ms], corresponding to the timeframe in which the stimulus is presented). The N1 is the prominent ERP component related to encoding processes, and it is usually centred in the [100 ms,150 ms] window. To investigate decision processes, we also analyzed the [300 ms, 500 ms] window in central channels (FC1, FCz, FC2, C1, Cz, C2), associated with the P2 component [23], the second positive deflection after stimulus onset (note that in some experimental settings the component associated with decision processes can be referred as P3, due to the presence of an additional positive deflection occurring after the N1, that we did not observe in analyzed signals). While this positive deflection usually occurs at around 300 ms after stimulus presentation, we anticipated the P2 to peak at around 400 ms, due to age effects [39]. N1 and P2 ERP amplitudes were determined by computing the integral in their respective windows. We also analyzed the N1 negative peak amplitude (N1 peak), while we neglected analysis of latencies (timings of N1 peak value) since they did not present significant differences between groups (H = 3.76, p = 0.91).

Functional connectivity analysis

Extra-occipital Functional connectivity (FC) is a measure of seeded FC [29] computed as the percentage of channels recruited beyond the occipital channel cluster during the encoding timeframe. In the main text, extra-occipital FC is simply referred as FC for the sake of brevity. An FC value close to 0% indicates no high positive correlations (we set the threshold for high positive correlation equal to r = 0.75) between the occipital channels and other channels. A value of 100% indicates that the whole scalp is correlated (with r > 0.75) with the occipital channels. Only significant correlations (Bonferroni-Holm p < 0.05) were included, with computation performed using the Spearman non-parametric measure. Spearman was chosen due to its low computational impact and its suitability for computing FC between macro-areas, such as the occipital region and the rest of the brain, thereby mitigating volume conduction effects [29]. For simulated signals, we computed FC using a similar methodology between local activity of brain regions [29].

We tested interaction between the N1 component values and categorical FC values with two-ways ANOVA, finding significant interaction (see Results, “Two different regimes of functional connectivity determine neural activity evolution along the cognitive decline continuum” paragraph). This allowed us to define Low and High FC subgroups, that we identified by using the median FC value as boundary.

Brain network model

The brain network model (BNM) utilized in this study has been extensively discussed in a previous publication from our group [40]. In summary, the brain was modelled as a network of 76 interacting region using The Virtual Brain platform or TVB [41]. Local neural activity was depicted by a modified version of the Jansen-Rit neural mass model [42], and the network connectivity was derived from the TVB connectome, which combines high precision and reliability [43]. The network connectivity matrix can be found in Supplementary materials (Supplementary Fig. 1). Local activity was projected onto a virtual grid of electrodes on a virtual scalp to simulate EEG signals. The mathematical framework for EEG simulations [44] consists of solving Maxwell equation in the dipole approximation across several layers of different conductance (cerebral tissue, cranium, scalp). The progression of neurodegeneration is regulated by two parameters, $cp$ and $lp$, that respectively describe connectivity and synaptic degeneration in the brain. Increasing either $cp$ or $lp$ results in a network organization characterized by more severe structural alterations. EEG signals were derived from local activity following the methodology introduced in Sarvas et al. [44], by computing electric field propagation in the dipole approximation [43] through a boundary element method model of the interfaces between three layers of different conductance, representing brain, skull and scalp.

Model equations

Local activity of network nodes was modeled with the Jansen-Rit neural mass model [42]. The Jansen-Rit model portrays activity in neural structures using a mean-field approach. Network nodes are depicted as ensembles of three sub-populations, inspired by the architecture of the brain cortex. These sub-populations include pyramidal cells (excitatory), stellate cells (excitatory), and inhibitory interneurons (inhibitory). Mathematically, the activity of these sub-populations for each node is described by solutions of a set of stochastic differential equations. The neural activity of a node is transmitted as input to the equations of other nodes, with coupling specified by the network connectivity matrix. The equations of the model are the same of the standard Jansen-Rit model [42]:

$$\left\{\begin{array}{c}\frac{d}{dt}{y}_{0}(t)={y}_{3}(t);\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\frac{d}{dt}{y}_{3}(t)=\frac{A}{{\tau }_{e}}Sigm\left[{y}_{1}(t)-{y}_{2}(t)\right]-\frac{2}{{\tau }_{e}}{y}_{3}(t)-{\left(\frac{1}{{\tau }_{e}}\right)}^{2}{y}_{0}(t)\\ \frac{d}{dt}{y}_{1}(t)={y}_{4}(t);\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\frac{d}{dt}{y}_{4}(t)=\frac{A}{{\tau }_{e}}\left\{p\left(t\right)+{C}_{2}\text{Sigm}\left[{C}_{1}{y}_{0}\left(t\right)\right]\right\}-\frac{2}{{\tau }_{e}}{y}_{4}(t)-{\left(\frac{1}{{\tau }_{e}}\right)}^{2}{y}_{1}\left(t\right) \\ \frac{d}{dt}{y}_{2}(t)={y}_{5}(t);\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\frac{d}{dt}{y}_{5}(t)=\frac{B}{{\tau }_{i}}{C}_{4}Sigm\left[{C}_{3}{y}_{0}(t)\right]-\frac{2}{{\tau }_{i}}{y}_{5}(t)-{\left(\frac{1}{{\tau }_{i}}\right)}^{2}{y}_{2}(t)\end{array}\right.$$

(1)

Here, ${y}_{0}$, ${y}_{1}$, and ${y}_{2}$ represent the post-synaptic activities (measured in mV) of the pyramidal, stellate, and interneuron sub-populations, respectively, while ${y}_{3}$, ${y}_{4}$, and ${y}_{5}$ represent their derivatives. Each post-synaptic activity is described by an equation of the type $\ddot{y}\left(t\right)=\alpha ax\left(t\right)-2a\dot{y}\left(t\right)-{a}^{2}y(t)$, corresponding to a damped oscillator with external driving force, where the $a$ term has the dimension of ${s}^{-1}$ and $\alpha$ is the amplitude of external input. The $p\left(t\right)$ term represents the input, written as a firing rate (measured in s⁻¹). This term combines the activity coming from other network nodes $\mu \left(t\right)$ and a stochastic term $\sigma (t)$ modeling physiological neural noise: $p\left(t\right)= \mu \left(t\right) + \sigma (t)$. The $\sigma (t)$ term was modelled as an additive Gaussian white noise. The sigmoid term is a nonlinear gain term, used to convert the average post-synaptic activity into a mean firing rate.

$$\text{Sigm}\left(v\right)=\frac{{\upsilon }_{max}}{1+{e}^{r\left(v-{v}_{0}\right)}}$$

(2)

In this equation, ${\upsilon }_{max}$ is the maximum firing rate of the given sub-population, ${v}_{0}$ is the value of the potential corresponding to a firing rate of 50% the max value, and $r$ is the sigmoid slope at ${v}_{0}$ (measured in mV⁻¹). The local connectivity constants ${C}_{1}$, ${C}_{2}$, ${C}_{3}$, ${C}_{4}$ are proportional to the number of synapses linking the sub-populations together. The constants $A$ and $B$, expressed in mV, are the maximal activity amplitudes for excitatory and inhibitory neurons, respectively. The ${\tau }_{e}$ and ${\tau }_{i}$ parameters, expressed in ms, are the time constants of excitatory and inhibitory synapses, combining dynamic features shaping the time-course of synaptic transmission.

These last two parameters were altered to model neurodegeneration effects on synaptic transmission (see next paragraph). Specifically, tuning the $\frac{{\tau }_{e}}{{\tau }_{i}}$ ratio is a straightforward method for implementing the excitation/inhibition imbalance associated with pathological cognitive decline [45]. Model parameters are summarized in Supplementary Table 2.

The model captures the local dynamics of each single node of the network. The nodes are then linked using the network connectivity matrix reported in Supplementary Fig. 1. In the network connectivity matrix, each pair of nodes is characterized by a connective weight ${C}_{weight}$ that mathematically transcribes the coupling between brain regions, derived from high-fidelity imaging data [43]. The weights ${C}_{weight}$ of the TVB connectome are assigned integer values ranging from 0 to 3, respectively transcribing null, weak, medium and strong anatomical connections.

The ${C}_{weight}$ connective values of the native network connectivity matrix were altered to describe neurodegeneration and neuroplasticity effects.

Virtual neurodegeneration and simulated EEG signals

In our computational Brain Network Model, we represented the brain as a network of nodes, with each node corresponding to a brain region, and edges representing the structural connections between these regions. Edges values were derived from a high-fidelity network connectivity matrix [43] whose entries represent connective weights ${C}_{\text{weight }}$ between regions. Initially BNM parameters were determined matching empirical data obtained from healthy individuals. Subsequently, we introduced virtual network degeneration to simulate the most common form of neurodegeneration underlying cognitive impairment, particularly Alzheimer-type dementia. We introduced the alterations in the nodes corresponding to brain regions involved in initial phases of the disease (reported in Supplementary Fig. 2), derived from well-known Braak stages of disease progression [46].

The model parameters used to describe network degeneration are the synaptic degeneration parameter (referred to as local parameter or $lp$, since this parameter affects only the node to which it is applied) and the parameter of connectivity degeneration ($cp$). The $lp$ parameter affected local neural dynamics by accelerating excitatory synapses while slowing inhibitory synapses, following the equation:

$$\begin{array}{cc} & {\tau }_{i}\to {\tau }_{i}^{healthy}+lp\times \left({\tau }_{i}^{\text{max}}-{\tau }_{i}^{healthy}\right)\\ & {\tau }_{e}\to {\tau }_{e}^{healthy}+lp\times \left({\tau }_{e}^{\text{min }}-{\tau }_{e}^{healthy}\right) \end{array}, 0<lp<1$$

(3)

where ${\tau }_{i}$ and ${\tau }_{e}$ are the timescales of inhibitory and excitatory synapses, respectively. Small values of both variables indicate fast synapses and vice versa. The ${\tau }_{i}^{\text{max}}$ and ${\tau }_{e}^{\text{min}}$ represent the maximum value of ${\tau }_{i}$ and the minimum value of ${\tau }_{e}$, respectively and ${\tau }_{e}^{healthy}$- ${\tau }_{i}^{healthy}$ represent standard healthy values (see Supplementary Table 2).

Modeling synaptic neurodegeneration by altering the characteristic times of excitatory and inhibitory synapses was meant to capture early alterations in synaptic transmission inducing excitation/inhibition imbalance [45, 47]. In later stages, local disconnection can also cause alterations in local connectivity parameters of the model, namely ${C}_{1}$, ${C}_{2}$, ${C}_{3}$, ${C}_{4}$, that we did not alter to focus on the early stages of neurodegeneration.

The $cp$ parameter affected the network connectivity matrix by reducing long range connections, while simultaneously increasing short range connections, according to equations:

$$C_{\text{weight}}\rightarrow \left\{ \frac{C_{\text{weight}}^{HC}-cp\times C_{\text{weight}}^{max}, \text{if}\;C_{\text{length}}>C_{\text{length}}^{\text{th}}}{C_{\text{weight}}^{HC}+np\times cp\times C_{\text{weight}}^{\text{max}},\text{if}\;C_{\text{length}}< C_{\text{length}}^{th}}, \right. 0< cp < 1$$

(4)

Biophysically, $lp$ portrays hypo-inhibition and hyperexcitation of neural ensembles [45], while $cp$ portrays white matter atrophy and neuroplastic effects [48, 49], which reportedly increase connections between neighbouring areas after failure of high-degree nodes, or “hubs” of the brain network. Both parameters range between 0 and 1, with 1 representing the most severe condition. The magnitude of neuroplastic mechanisms is represented by $np$ and its default value is 1. The values of ${\tau }_{i}^{\text{max}}$ and ${\tau }_{e}^{\text{min}}$ in Eqs. 3 and 4 (${\tau }_{i}^{\text{max}}$=40 ms and ${\tau }_{e}^{\text{min}}$=8.9 ms) are deduced from biophysical constraints described in our previous work [40]. The value ${C}_{\text{weight }}^{\text{max}}=2$ was chosen so that the maximum connectivity damage (at $cp$ =1) would severe connections up to medium ones, without eliminating strong connections, codified in the structural connectivity matrix by the value ${C}_{\text{weight }}=$ 3. The only strong connections disrupted by Eq. 4 are those between affected nodes. This happens since $cp$ acts on the structural connectivity matrix by uniformly decreasing the afferent connections from cortical regions (rows of the structural connectivity matrix) and then the efferent connections (columns of the structural connectivity matrix) of affected nodes. This causes the submatrix corresponding to connective weights between affected nodes to be reduced twice by $cp$, causing strong connections to become eventually severed by neurodegeneration when $cp$ reaches the 0.75 value, since at this value the connections in the submatrix are reduced by $cp\times {C}_{\text{weight }}^{\text{max}}\times 2,$ equal to 3 when cp=0.75.

The nodes affected by Eqs. 3 and 4 are fourteen, derived from known stages of Alzheimer’s type dementia [46], which is the leading cause of pathological cognitive decline. The names and locations of such nodes can be found in Supplementary Fig. 2. Since the connective weights ${C}_{\text{weight}}$ of the network connectivity matrix are graded in integer steps from 0 to ${C}_{\text{weight }}^{\text{max}} = 3,$ a $cp$ value of 0.5 corresponds to rescission of weak connective fibers (${C}_{\text{weight }}=1,$ mainly from the entorhinal and temporal cortex), while a $cp$ value of 1,0 corresponds to rescission of medium connective fibers ${(C}_{\text{weight }}=2)$.

Virtual version of cognitive task

Our approach leverages a reservoir computing architecture, where the nodes of the network function as a dynamic reservoir that operates a non-linear transformation on the stimuli injected in the input layer (Supplementary Fig. 3). Connections between nodes are fixed and determined by human tractography data [43, 50]. The output, which corresponds to the simulated N1 components and FC, depends on the values of the $cp$ and $lp$ parameters, that were trained to minimize a loss function (reported in the Personalized Brain Network Models and parameter estimation paragraph). For each patient, $cp$ and $lp$ parameters were updated to minimize the distance between the virtually generated and experimentally observed ERPs, thus forming the readout layer of the model.

The experimental cognitive task involved recognizing a shape (a triangle) presented to the participant. To model the onset of visual stimulation and the encoding processes required for recognizing the presented shape, we injected a virtual stimulus into the node of the network corresponding to the primary visual cortex. Subsequently the stimulus propagated along regions of the ventral pathway of visual information [51]. Note that in this way we modelled the correct recognition of the TARGET shape, and therefore we considered only neural signals coming from trials in which the patients correctly recognized the shape.

Virtual stimulation was simulated through a feed-forward propagating square wave input, spreading across the regions of the ventral pathway of visual information, namely: Primary visual cortex (V1) > secondary visual cortex (V2) > Visual Area V4 (anatomically close to the Caudal portion of Inferior Temporal Cortex) > Rostral portion of Inferior Temporal Cortex. Both the regions and the time delays in the feed-forward propagation of the stimulus were deduced from previous findings by Riesenhuber and Poggio [51]. Nodes corresponding to these regions represent the input layer of our BNM. The activity elicited by the stimulus is then propagated to the whole network, with an intensity determined by the entries of the structural connectivity matrix. The virtual stimulus was injected into the corresponding nodes of the network as an additive term in $p\left(t\right)$ in Eq. 1.

Other stimulus parameters, namely duration time and amplitude of the square wave, were carefully selected through parameter exploration. The amplitude range was set from 0.01 to 5.0 μV, while the duration range went from 10 to 200 ms (the latter matching the experimental duration of stimulus presentation). The chosen values were selected to best match the average experimental N1 peak value, which emerged as the most robust neural activity correlate of cognitive decline progression (see Results, “Neural activity presents non-monotonic evolution with cognitive decline severity” paragraph). The stimulus with selected parameters consisted in a square wave of fixed duration and amplitude injected in V1. Stimuli injected in other nodes of the input layer consisted in square waves of lower amplitude and increased duration when compared to the stimulus injected in V1. Precise values of stimuli parameters can be consulted from the GitHub repository reported in the Code and Data availability paragraph. Gaussian white noise, generated with a random seed, was superimposed onto this stimulus to account for biological neural noise. To account for biological variability of EEG signals, we conducted 20 simulations with different noise seeds and then computed the mean virtual ERP by averaging the results. This procedure was repeated for each virtual ERP generation.

Graph-theoretical analysis

We investigated the topology of the BNM using standard graph-theoretical quantities [52, 53] such as Efficiency and Assortativity index. Efficiency quantifies how effectively information is transmitted across a network by measuring the average inverse shortest path length between nodes; while Assortativity describes the tendency of nodes in a network to connect with other nodes that have similar properties, such as degree. The general equation for efficiency is [52, 54]:

$$Efficiency=\frac{1}{n(n-1)}\sum\limits_{i\neq j}1/d(i,j)$$

(5)

where $n$ represents the number of edges involved in the sum, and $d(i,j)$ is the minimum distance between nodes $i$ and$j$. In the case of a cortical structural connectivity matrix, $1/d\left(i,j\right)$ can be represented by the connective weight ${C}_{\text{weight}}$ between nodes $i$ and $j,$ also considering indirect connections (e.g. suppose that $i$ and $j$ are disconnected but are both connected to $k$ with ${C}_{\text{weight }}(i,k)$ = 3 and ${C}_{\text{weight }}\left(j,k\right)=2$: we have that the Efficiency of transmission between $i$ and $j$ is equal to $\frac{1}{d\left(i,j\right)}= {}^{6}\!\left/ \!{}_{5}\right.$ with $d\left(i,j\right)={1/C}_{\text{weight }}\left(i,k\right)+ 1/{C}_{\text{weight }}\left(j,k\right)= {}^{5}\!\left/ \!{}_{6}\right.).$

Assortativity was computed using the Randic index, a generalized measure of assortativity in network theory [52]. We chose this index as it has previously been utilized to investigate structural network degeneration in prodromal phases of dementia and Alzheimer’s disease [4]. The Assortativity equation is as follows:

$$Assortativity=\sum\limits_{i\neq j} \left(\sum\limits_{k}\left(C_{weight}^{i,k}C_{weight}^{j,k} \right)^a \right)$$

(6)

This corresponds to the sum of weighted degrees product over each pair of nodes of the network connectivity matrix, with the exponent $\alpha$=0.25.

Comparison of experimental and simulated event related potentials

Experimental and simulated ERP components were normalized within their respective windows [100 ms, 150 ms] for N1 and [300 ms, 500 ms] for P2. We then calculated the average difference between the experimental and simulated ERPs, comparing this difference with the standard deviation $({\sigma }_{sd}^{exp})$ of experimental signals in the respective window. This procedure was conducted separately for each group. Since the model could not replicate the peak values of the P2 component, an offset was added to simulated P2 components to match simulated signals to experimental recordings.

Personalized brain network models and parameter estimation

We employed a stochastic gradient-descent-based algorithm to find the values of degeneration parameters that best replicated individual EEG recordings from each participant. The loss function of the algorithm was determined by the difference between N1 peak value and FC between simulated and experimental signals, according to the equation:

$$Loss=\sqrt{Weight\times(FC^{exp}-FC^{sim})^2+(N1^{exp}-N1^{sim})^2}$$

(7)

Were $N1$ represents the N1 peak value, $FC$ represents the FC value, and $Weight$ is a quantity determined to minimize the loss function, with optimal value found to be $Weight$=0.1. The exponents $exp$ and $sim$ represent experimental and simulated values, respectively. Both $cp$ and $lp$ values were updated with steps of 0.02 in the gradient descent algorithm. The gradient descent algorithm also included 5 stochastic thermalization steps at each update of the parameter values, to avoid getting stuck in local minima. Thermalization steps were uniformly distributed in the (−0.05, 0.05) interval.

Differently from our previous work [40], we tested several possible values for the magnitude of neuroplastic effects ($np)$ and the number of affected nodes. Particularly, we considered $np$ values between 0.5 and 1.5, and all configurations in which disconnected nodes were less than the original fourteen (disconnected nodes thus ranged between one and fourteen). Participants from High and Low FC subgroups were fitted to parameters from different sides of the $cp$=0.5 line ($cp$<0.5 for the High FC subgroup, $cp$>0.5 for the Low FC subgroup). The optimal values for both $np$ and number of disconnected nodes were determined by minimizing the cumulative loss function of Eq. 7 computed on all participants. Briefly, we determined the personalized BNMs from all possible configurations, constructed with each possible combination of $np$ and number of disconnected nodes. We then computed the cumulative loss function for each configuration as the sum of loss functions for all participants. The optimal values of $np$ and number of disconnected nodes were those that minimized the cumulative loss function. Optimal value for $np$ was found to be 1, while the optimal affected nodes are reported in Supplementary Fig. 2.

Statistical tests

Statistical comparisons were conducted using the Kruskal–Wallis H-test statistics, followed by post-hoc analysis with the Mann–Whitney U-rank test. Statistical results were corrected with Bonferroni-Holm technique, setting significance threshold for corrected p-value at p < 0.05.

For each EEG signal recorded on occipital and central channels, we measured differences across groups for each timestamp using KW H-test. The interaction between FC subgroups and N1 integral values was evaluated using a two-ways ANOVA.

The bimodality of FC and $cp$ distributions was assessed with the Dip-Hartigan’s test. Significance for Dip-Hartigan’s tests was determined by bootstrapping 250 times to create a null distribution, comparing observations with the 95-percentile value of the null distribution (p = 0.05). To determine the distribution of personalized brain network parameters across groups, we used a kernel density estimation routine. Comparisons of the peak population histograms in group-wise distributions of brain network parameters were conducted by bootstrapping 100 times for each group and computing the population inside peak area. A significance level of p < 0.05 was determined if the mean value of areas obtained from bootstrapping group A distribution fell outside of the 2.5th to 97.5th percentile range of group B. If the 2.5th to 97.5th percentiles of group A and B did not overlap, p < 0.0025 was determined. We conducted a power analysis based on Mann Whitney U test to assess the robustness of statistical comparisons between CTR group and other groups, since the CTR group presented a limited number of subjects. The effect size was estimated as the standardized rank difference between groups, and power was computed using a two-sided test with a significance level of 0.05.

All analyses, except for pre-processing discussed in the dedicated section, were performed in Python, using commercially available packages such as Numpy, Pandas, Seaborn and Scikit-Learn. Post-hoc analysis and two-ways ANOVA test were conducted with the Statsmodel package, while all other statistical analyses were performed using the Scipy package.

Outlier management

Extreme outliers, defined as datapoints distant more than 3 ${\sigma }_{sd}$ from the closest datapoint, were excluded. Here ${\sigma }_{sd}$ represents the standard deviation of the sample. With this approach, we encountered only one outlier in experimental N1 component for the MCI group, reducing the group size to n = 39 patients (Fig. 2), and two outliers in the CTR group in the pipeline for personalized BNM determination, reducing the group to n = 17 subjects (Fig. 5).

Fig. 2

Neural activity presents two distinct regimes of functional connectivity. A Stimulus-onset aligned ERP in occipital channels, for CTR (blue), SCD (green) and MCI (red) participants. Dashed lines indicate the encoding timeframe over which N1 is computed. Shaded areas indicate standard error of the mean. B N1 integral values across groups, Dots indicate single participant values, markers indicate significant difference (p < 0.05), Boxes represent lower and upper quartile values in both plots, while whiskers are set at a value corresponding to 1.5 times the interquartile range. C Boxplot of total FC values across groups, notation is the same as for (B). D Distribution of FC values for the CTR group, black vertical line represents the median FC value. E Distribution of FC values for the SCD group, with clear bimodal distribution, as confirmed by DH test. Notation is the same as (D). F Distribution of FC values for the MCI group. Notation is the same as (D). SCD is the only group displaying bimodality in FC values. G Boxplot of N1 integral for the High FC subgroup. Notation is the same as in (B). H Average N1 components in occipital channels for the High FC subset, divided by groups. N1 displays monotonic decrease with cognitive decline severity. Same notation as in (A). I Boxplot of N1 integral for the Low FC subgroup. Dots indicate single participant’s values; markers indicate significant difference (p < 0.05). J Average N1 components in occipital channels for the Low FC subset, divided by groups. The N1 ordering is the same as the non-monotonic one observed in (A). Abbreviations: CTR: Control; FC: Functional Connectivity; MCI: Mild Cognitive Impairment; SCD: Subjective Cognitive Decline

Participants age

Patients in the MCI group were older (on average) than the other groups (Supplementary Table 1). We ensured that age did not present a confounding factor for our analysis. This was done by measuring correlation between relevant metrics in neural activity and participants age, revealing no relevant relations between the two quantities.

Results

Neural activity presents non-monotonic evolution with cognitive decline severity

Our goal was to implement a brain network model reproducing neural activity during cognitive task execution across various stages of cognitive decline. We recorded and analyzed EEG signals (see Methods, “EEG pre-processing and ERP computation” paragraph) from participants in different stages of cognitive decline during the execution of a cognitive task. We included healthy control subjects (CTR) Subjective Cognitive Decline patients (SCD) and Mild Cognitive Impairment patients (MCI). Task performances decreased monotonically with condition severity (see Vergani et al. [55] for more details), group averaged performances are reported in Supplementary Table 3.

We analyzed EEG recordings, extracting the P2 and N1 components (see Methods, “EEG pre-processing and ERP computation” paragraph), related to proficiency in decision and encoding processes, respectively [23].

We computed P2 amplitudes as the integral in the [300 ms, 500 ms] window from central channels. Since the P2 amplitude values presented no relevant differences across groups, they were discarded from further analyses and from the fitting of model parameters.

We computed the N1 amplitude as the integral in the [100 ms, 150 ms] window from occipital channels ((see Methods, “EEG pre-processing and ERP computation” paragraph and Fig. 2A), comparing values across groups (Fig. 2B). The SCD group displayed significantly smaller N1 integral values compared to the HC group, but not compared to MCI (H = 10.13, p = 0.006; post-hoc comparison: CTR vs SCD, U = 1165.0, p = 0.011; SCD vs MCI U = 1410.0, p = 0.10). Since the CTR group presented significantly fewer subjects than the SCD one, we assessed the statistical power of these results, performing a post-hoc power analysis based on the Mann–Whitney U test (see Methods, “Statistical Tests”, paragraph). The analysis yielded an estimated power of 0.74, corroborating the results of the test. Interestingly, the N1 integral did not differ significantly between CTR and MCI (U = 505.0, p = 0.087), indicating a non-monotonic evolution of this feature with cognitive decline severity.

We then compared the average encoding functional connectivity between occipital and extra-occipital channels (from now on FC, see Methods, “Functional Connectivity analysis” paragraph) across the three groups. Overall differences between groups were significant (Fig. 2C, H = 10.23, p = 0.006), primarily due to disparities between SCD and MCI (post-hoc analysis U = 1224.0, p = 0.016), while differences between other groups were non-significant (SCD vs CTR: U = 1078.0, p = 0.057; CTR vs MCI: U = 404.50, p = 0.817).

Notably, FC also exhibited a non-monotonic evolution with condition severity, as it was observed in N1 integral values. We also analyzed the shape of group-wise FC distributions (Fig. 2D-F). Interestingly, the SCD group exhibited a clear bimodality in FC values (DH test = 0.78, p = 0.0001, evaluated by a Dip-Hartigan’s test, see Methods, “Statistical tests” paragraph and Fig. 2E) in contrast to both CTR (DH test = 0.70, p = 0.56, Fig. 2D) and MCI (DH test = 0.75, p = 0.089, Fig. 2F) groups. We assessed the correlation between N1 integral values and age, which revealed no relevant relation (r = 0.09, p = 0.30). We also checked the correlation between extra-occipital FC and age with similar results (r = 0.14, p = 0.09).

Two different regimes of functional connectivity determine neural activity evolution along the cognitive decline continuum

We found a significant interaction between cognitive decline severity and FC in determining the N1 integral values (Two-ways ANOVA, F = 16.29, p < 0.0001, see Methods, “Statistical tests” paragraph).

We divided participants into High and Low FC subgroups, using the median value as boundary. In the High FC subgroup, no significant differences in terms of N1 integral values were found across groups (Fig. 2G-H). In contrast, significant differences were observed in the Low FC subgroup (Fig. 2I-J, H = 8.64, p = 0.013). Post-hoc analysis revealed significant differences only between MCI and SCD (U = 146.0, p = 0.039) and not between other group pairs.

In the High FC subgroup, the ranking of N1 integral values was CTR > SCD > MCI (Fig. 2H), consistent with the expected monotonic decrease in ERP components associated with increasing cognitive decline severity (38). However, in the Low FC subgroup, the ranking was CTR > MCI > SCD (Fig. 2J), exhibiting the same non-monotonic evolution along the cognitive decline continuum observed in Fig. 2A.

Simulated network connectivity degeneration leads to a phase transition between two distinct regimes of network topology and functional connectivity

We adapted a BNM previously implemented to simulate EEG signals at rest [40] to reproduce EEG activity during cognitive task execution. Task condition was translated into the model by injecting a cascading virtual stimulus in nodes of the network corresponding to the ventral pathway of visual information (see Methods, “Virtual version of cognitive task” paragraph). The BNM represents neurodegeneration related to cognitive decline through two parameters, $lp$ and $cp$, which capture synaptic and connectivity degeneration, respectively (see Methods, “Virtual neurodegeneration and simulated EEG signals” paragraph).

Experimental results showed that FC plays a key role in determining neural activity changes during cognitive decline. We utilized the BNM to study how neurodegeneration affects FC in response to varying stimulus intensities. For low connectivity degeneration (cp = 0.3, lp = 0.6, Fig. 3A), FC rapidly increased, doubling at 0.05 μV and saturating at 1.00 μV (mean FC = 0.88). With higher connectivity degeneration (cp = 0.7, lp = 0.6, Fig. 3B), saturation was lower (mean FC = 0.80), occurring at 1.75 μV. These findings suggest that greater connectivity degeneration reduces FC's response to increasing stimulus intensity. Synaptic degeneration did not affect FC response to stimulus intensity.

Fig. 3

Simulated network connectivity degeneration determines a phase-transition in network topology with cascading effects on functional connectivity. A Average FC for various level of stimulus intensity at fixed $cp$ and $lp$, $cp$<0.5. FC steeply increase with stimulus intensity (see inset). B Average FC for various level of stimulus intensity at fixed $cp$ and $lp$, $cp$>0.5. C Simulated FC according to $cp$ and $lp$ values. FC is computed as percentage of relevant functional connections from occipital regions to extra-occipital regions. Each point represents the value of the FC computed from a simulation with a given combination of $cp$ and $lp$. Values are normalized for better visualization. D Efficiency in the brain network. Blue edges represent high Efficiency links, while red edges represent low Efficiency links. E Assortativity in the brain network. Notation is the same as (D). F Efficiency values according to $cp$ for different values of number of disconnected regions. Red represents the value used for simulations. At $cp$=0.5 there is a first-order phase transition dependent on the number of disconnected regions. Red curve represents the value used for simulations. G Efficiency values according to $cp$ for different values of neuroplasticity. Same notations as in (F). H Linear regression between magnitude of phase transition in Efficiency and number of disconnected nodes. Red box represents the value selected for simulations. I Logarithmic fit capturing the relationship between the magnitude of phase transition in Efficiency and neuroplasticity. Notation is the same as in (H). J Assortativity according to $cp$ for different values of number of disconnected regions. Same notation as in (F). K Assortativity according to $cp$ for different values of neuroplasticity. Same notations as in (F). L Linear regression between average Assortativity and number of disconnected nodes. Error bars represent standard deviation values. Red box represents the value selected for simulations. M Logarithmic fit between neuroplasticity and Assortativity values. Notation is the same as in (L). Abbreviations: $cp$: connectivity parameter; FC: Functional Connectivity; $lp$: local parameter

Next, we fixed the stimulus intensity at 5 μV and investigated how simulated FC depended on $lp$ and $cp$ parameters (Fig. 3C). The effect of $cp$ on simulated FC values was pronounced. Specifically, $cp$=0.5 determined a first-order phase transition between two distinct FC regimes for all $lp$ values. While stimulation resulted in a rapid increase to a saturation value at both sides of the $cp$=0.5 line, this value was significantly lower in the $cp$>0.5 zone (H = 14,079.75, p < 0.0001).

To quantify how the phase transition depended on underlying network structure, we investigated how the $cp$ parameter affected network topology, in terms of both Efficiency and Assortativity (see Methods, “Grap-theoretical analysis” paragraph). Remarkably, we observed a first-order phase transition in network Efficiency (Eq. 5, Fig. 3F) at $cp$=0.5, which coincided with the phase transition in simulated FC. Conversely, Assortativity (Eq. 6, Fig. 3J) exhibited an inverted U-shaped evolution with a peak at around $cp$=0.2.

The $cp$ parameter condenses two connectivity alterations (see Methods, “Virtual neurodegeneration and simulated EEG signals” paragraph): (i) a decrease in long-range connections in selected nodes [56] and (ii) connectivity rewiring due to neuroplastic mechanisms [57]. The last one causes a reinforcement of short-range connections (from now on called “neuroplasticity” for the sake of brevity). We investigated the role played by these two alterations in shaping the phase transition, assessing their separate effects on both Efficiency and Assortativity.

The phase transition in Efficiency occurred for number of disconnected nodes > 1 (Fig. 3F) and for all values of the neuroplasticity term ($np$ in Eq. 4, Fig. 3G). To quantify the relation between number of disconnected nodes and the magnitude of phase transition (i.e. Efficiency discontinuity, simply computed as the difference in Efficiency at different sides of the transition value), we conducted a linear regression analysis, finding a highly significant positive correlation (r = 0.98, p < 0.00001, Fig. 3H).

Interestingly, a similar result was found for FC, whose phase transition was also strongly affected by number of disconnected nodes, with an exponential dependence (Supplementary Fig. 4).

Neuroplasticity, on the other hand, had minimal impact on shaping the phase transition, with Efficacy discontinuity values ranging from 1.09 to 1.18 for all tested neuroplasticity values (range for number of disconnected nodes was 0.04–1.18). The relationship between Efficiency discontinuity and neuroplasticity was accurately described by a discontinuous logarithmic function (r² = 0.98, Fig. 3I) with a jump for neuroplasticity = 1.0.

The inverted U-shaped evolution of Assortativity while increasing $cp$ was present for all tested values of number of disconnected nodes (Fig. 3J) or neuroplasticity (Fig. 3K). However, the number of disconnected nodes significantly reduced mean Assortativity levels (linear regression analysis, r = −1.00, p < 0.00001, Fig. 3L). Conversely, neuroplasticity increased Assortativity levels, albeit only logarithmically (r² = 0.95, Fig. 3M). These results underscore the interplay between network topology (as measured by network Efficiency) and functional connectivity, as they both exhibited a phase transition for the same value of connectivity degeneration.

The connectivity phase transition leads to distinct evolutions of neural activity across lp spectrum

The experimental N1 ERP component strongly interacted with FC values, with participants in the Low and High FC subgroups exhibiting distinct N1 evolutions along the cognitive decline continuum (Fig. 2G-J). To determine if the model could capture such dynamics, we conducted a grid search on simulated N1 integral values (similar to what we did for FC, see Fig. 3C) to assess how this quantity depended on $lp$ and $cp$ parameters.

Simulated N1 component displayed a transition at $cp$=0.5 (Fig. 4A, for the same value that caused a transition in FC). For $cp$<0.5, N1 integral values decreased monotonically in absolute value while increasing $lp$ (linear regression, r = 0.86, p < 0.00001, Fig. 4B). For $cp$>0.5, N1 integral values displayed an inverted U-shaped evolution with $lp$, which can be fit to a Gaussian function (r² = 0.33, Fig. 4C).

Fig. 4

Simulated signals with high and low connectivity degeneration accurately reproduce observed N1 evolutions in Low FC and High FC subgroups. A Absolute values of N1 integral from simulated signals according to $cp$ and $lp$ values. Each point represents the value of the N1 integral computed from occipital channels from a simulation with a given combination of $cp$ and $lp$. Values are normalized for better visualization. B Simulated N1 integral values according to $lp$, for each $cp$<0.5. Values were normalized between 0 and 1 for each $cp$ separately. Red line represents the best fitting line of the linear regression analysis (r = 0.86, p < 0.00001). Darker blue represents higher $cp$ values. C Simulated N1 integral values according to$lp$, for $cp$>0.5. Notation is the same as in (B). Red curve represents the best fitting function, a gaussian (r² = 0.33). While for $cp$<0.5, N1 integral monotonically decreased with $lp$, for $cp$>0.5 the evolution resembles an inverted-U shape. D Simulated N1 components in occipital channels for $lp$ values of escalating intensity, $cp$<0.5. E Simulated N1 components in occipital channels for $lp$ values of escalating intensity, 0.5 < $cp$<0.75. F Peak values of simulated N1 components for $lp$ values of escalating intensity and $cp$<0.5 (with the same parameters used in Panel D) compared with experimental peak values of N1 components in the Low FC subgroup. As the $lp$ value increases, simulated N1 component peaks exhibited the same monotonic decrease as experimental N1 component peaks with increasing cognitive decline severity. G Peak values of simulated N1 components for $lp$ values of escalating intensity and 0.5 < $cp$<0.75 (with the same parameters used in Panel E) compared with experimental peak values of N1 components in the High FC subgroup. Simulated and experimental signals present the same inverted U-shaped evolution. Abbreviations: $cp$: connectivity parameter; CTR: Control; FC: Functional Connectivity; $lp$: local parameter; MCI: Mild Cognitive Impairment; SCD: Subjective Cognitive Decline

The monotonic decrease of N1 integral values in the $cp$<0.5 zone mirrored the observed experimental decline of N1 along the cognitive decline continuum in the High FC subgroup. Conversely, within the $cp$>0.5 zone, simulated signals generated from increasing $lp$ values paralleled the inverted U-shaped evolution observed in low FC participants as their cognitive decline severity increased (MCI > CTR > SCD).

These findings suggest that the phase transition at $cp$=0.5 might be the structural basis for the distinct neural activity patterns observed in the two FC subgroups. To explore this hypothesis, we selected parameter combinations to replicate the experimentally observed N1 patterns in both Low FC and High FC subgroups.

Strikingly, BNMs with $lp$ values of escalating intensity in the $cp$<0.5 replicated the N1 evolution observed in the High FC subgroup (Fig. 4D, compare with Fig. 2H). Conversely, $lp$ values of escalating intensity at $cp$>0.5 replicated the N1 evolution observed in the Low FC subgroup (Fig. 4E, compare with Fig. 2J). These results demonstrate distinct N1 evolutions determined by synaptic degeneration levels at the two sides of the $cp$ phase transition. The similarity between simulated signals on different sides of the $cp$=0.5 line and experimental signals from the two FC subgroups can also be appreciated in terms of N1 peak values (Fig. 4F-G). For values of $cp$>0.75 the increase in $lp$ does not affect neural activity, effectively creating a third possible regime (see Supplementary Fig. 5). Model parameters used in the simulations of all three regimes are reported in Supplementary Table 4.

The model replicates neural activity during task execution across cognitive decline stages

The model captured the evolution of both FC and N1 components across cognitive decline stages, suggesting possible neurodegeneration mechanisms. We then tested the model’s ability to reproduce quantitatively the experimental EEG signals during cognitive task execution, by comparing ERP components extracted from simulated signals with those extracted from experimental recordings. Despite the experimental P2 components presented no significant differences across groups, we included this component (along with N1) in the analysis.

We selected model parameters that reproduced the group-averaged values of ERP components (comparison between simulated and experimental ERP components is reported in Supplementary Figs. 6–7, while parameter values used in the simulation in Supplementary Table 5). The fidelity of simulated ERP components was determined by computing differences with respect to their experimental counterparts point by point (using as a metric the Euclidean distance). We then compared the differences with the standard deviation ${\sigma }_{sd}$ of experimental signals.

For N1, we found that the differences were 0.57 ${\sigma }_{sd}^{CTR}$ for the CTR group, 0.26 ${\sigma }_{sd}^{SCD}$ for the SCD and 0.81 ${\sigma }_{sd}^{MCI}$ for the MCI group (Supplementary Fig. 6). For P2, the differences were 0.77 ${\sigma }_{sd}^{CTR}$ for the CTR group, 2.33 ${\sigma }_{sd}^{SCD}$ for the SCD and 2.11 ${\sigma }_{sd}^{MCI}$ for the MCI group (see Supplementary Fig. 7). This indicates high accuracy in replicating the N1 component, while larger discrepancies are found for the P2, reflecting more complex interaction with high-order brain areas associated with this component (23) not included in the model.

Synaptic degeneration determines cognitive decline severity in networks close to the connectivity phase transition

We proceeded to reconstruct the personalized BNM for each participant, by identifying the combination of model parameters that best replicated individual EEG recordings at each side of the phase transition (see Methods, “Personalized Brain Network Models and parameter estimation” paragraph).

BNMs clustered at the connectivity phase transition (55/145 participants, ${\chi }^{2}$=8.45, p = 0.004). BNMs constructed from CTR subjects predominantly exhibited either extremely low $cp$ values or $cp$ values close to the transition boundary (Fig. 5A). In contrast, SCD BNMs exhibited an additional peak at the extremal value $cp$=1, (Fig. 5B), while MCI BNMs displayed scattered values with a peak near the phase transition (Fig. 5C). Notably, group-wise $cp$ distributions mirrored the bimodality of FC values, with the SCD group being the only bimodal group observed (CTR: DH test = 0.099, p = 0.096; SCD: DH test = 0.078, p = 0.0003; MCI: DH test = 0.05, p = 0.64).

Fig. 5

Personalized Brain Network Models present configurations with levels of synaptic degeneration proportional to cognitive decline severity. A Distribution of $cp$ parameters for BNMs constructed from CTR subjects. B Distribution of $cp$ parameters for SCD BNMs. This is the only bimodal distribution among all groups, consistently with experimental FC results. C Distribution of $cp$ parameters for MCI BNMs. D Assortativity of CTR BNMs. Blue dots represent Assortativity values of CTR BNMs, with dot width representing numerosity. Red dots are Assortativity values of the network at different $cp$ values. E Assortativity of SCD and MCI patients. Brown dots represent Assortativity values of SCD and MCI BNMs combined, notation is the same as in (D). F Distribution of Assortativity values for CTR BNMs, the distribution is extremely narrow. G Distribution of Assortativity values for SCD BNMs, displaying high dispersion. H Distribution of Assortativity values for MCI BNMs. I Distribution of personalized BNM parameters for CTR subjects. The distribution is unimodal, centred at low values of network degeneration. J Distribution of personalized BNM parameters for SCD patients. The distribution is bimodal, with a peak at low values and a peak at intermediate values of network degeneration. K Distribution of personalized BNM parameters for MCI patients. The distribution is unimodal, centred at intermediate values of network degeneration. L Population histogram of BNMs near the peak at low $lp$ levels ($lp$=0.12) divided per group. Histogram was computed from the personalized values of $cp$ and $lp$ parameters (see Methods, “Personalized Brain Network Models and parameter estimation” paragraph). Peak central $cp$ and $lp$ values are reported. M Population histogram of BNMs near the peak at intermediate $lp$ levels ($lp$=0.56) divided per group, same notation as in (L). N Population histogram of BNMs near the peak at high $lp$ levels ($lp$=0.96) divided per group, same notation as in (L). The peak analysis showed how participants from different groups clustered at levels of synaptic degeneration proportional to cognitive decline severity. Significance was computed by bootstrapping original values (see Methods, “Statistical tests” paragraph). Significance notation: * stands for p < 0.05, ** stands for p < 0.0025. Abbreviations: $cp$: connectivity parameter; CTR: Control; cv: coefficient of variation; $lp$: local parameter; MCI: Mild Cognitive Impairment; SCD: Subjective Cognitive Decline

Personalized $cp$ values also allowed for the reconstruction of personalized levels of Assortativity (Fig. 5D-E). Despite the dispersion in individual $cp$ parameters, BNMs constructed from CTR subjects exhibited a narrow distribution of Assortativity (coefficient of variation, cv = 0.06, Fig. 5F). Conversely, pathological groups displayed higher dispersion (SCD: cv = 2.34, Fig. 5G; MCI: cv = 0.55, Fig. 5H).

We investigated the group-wise distributions of both connectivity degeneration and synaptic degeneration parameters ($cp$ and $lp$). Utilizing a kernel density estimation routine, we determined the principal peak of the distribution for each group (Fig. 5I-K, see Methods, “Statistical tests” paragraph).

Crucially, personalized BNMs presented different levels of synaptic degeneration ($lp$ parameter) proportional to cognitive decline severity. While all group distributions presented a peak at $cp$ values near the transition value $cp$=0.5, groups of increasing cognitive decline severity peaked at progressively higher values of $lp$ parameters. Specifically, CTR BNMs presented its peak at $lp$=0.12 (Fig. 5I), while the SCD presented peaked at $lp$=0.56 (Fig. 5J), and MCI at $lp$=0.96 (Fig. 5K).

To corroborate this result, we observed that under the peak at $lp$=0.12 the percentage of BNMs reconstructed from CTR subjects was significantly higher than those reconstructed from the other two groups (CTR: 36.8 ± 6.4%; SCD: 20.9 ± 2.5%; MCI: 25.0 ± 3.9%, H = 619.1, p < 0.00001, Fig. 5L). Similarly, the percentage of SCD BNMs was significantly predominant in the peak at $lp$=0.56 (CTR: 26.3 ± 5.9%; SCD: 43.0 ± 3.1%; MCI: 25.0 ± 4.0%, H = 809.6, p < 0.00001, Fig. 5M), and the percentage of MCI BNMs was significantly predominant in the peak at $lp$=0.96 (CTR: 10.5 ± 4.1%; SCD: 15.1 ± 2.2%; MCI: 22.5 ± 4.1%, H = 578.5, p < 0.00001, Fig. 5N). Statistical tests and standard deviation for such percentages were determined by bootstrap technique (see Methods, “Statistical tests” paragraph).

These results corroborated our hypothesis about synaptic degeneration as the primary structural culprit determining cognitive decline severity at different sides of the phase transition. This establishes crucial connections between structure and function, bridging network neurodegeneration mechanisms to cognitive decline.

Discussion

Our BNM successfully replicated individual EEG responses to cognitive task execution across as a function of biologically-inspired neurodegeneration parameters. This marks the first example of a model reproducing individual neural activity, including relevant features such as N1 and P2 components as well as FC, especially across several stages of cognitive decline. The model used in this work has been previously validated on resting state EEG data to model alterations of power spectral density (particularly the alpha band) and resting-state functional connectivity [40]. The ability of capturing electrophysiological alterations both in the resting state and in task condition highlights the potential of the model in bridging between neural activity alterations and underlying neurodegeneration mechanisms.

Our model revealed a neurodegeneration-induced phase transition between two distinct neural activity regimes during task performance, successfully replicating experimental findings of distinct FC subgroups with unique N1 component trajectories. By incorporating structural parameters describing several neurodegenerative mechanisms [45, 46, 48], it uncovered the structural mechanisms driving these dynamics, offering insights into how varying neural activity patterns emerge from the underlying network architecture.

Specifically, connectivity degeneration $cp$ led to a phase transition in network Efficiency [52], a quantity representing the efficacy of information transmission between different brain regions [53]. This quantity is known to be affected by cognitive decline [4]. The phase transition in Efficiency triggered a cascading phase transition in functional connectivity, resulting in FC bimodality in simulated signals.

The N1 component evolutions on either side of the transition line closely mirrored the distinct patterns observed in participants from the Low and High FC subgroups along the cognitive decline continuum. This highlights the interplay between connectivity and synaptic degeneration as the key structural factors driving the different neural activity patterns among FC subgroups.

Interestingly, the $cp$=0.5 value that determined the phase transition in the BNM coincides with rescission of fibers in the posterior cingulum, a white matter bundle connecting the entorhinal cortex with cingular cortex (regions of the network relative to this bundle are the parahippocampal cortex, ventral and polar temporal cortex, see Supplementary Fig. 2). The posterior cingulum is a hub of the brain network [49], known to be heavily affected in cognitive decline [58, 59]. This suggests that the non-monotonic evolution of N1 component in the Low FC group (and thus in the $cp$>0.5 zone, see Fig. 4) may derive from neuroplastic mechanisms adopted after posterior cingulum damage. The model also identified a third regime of neural activity evolution at $cp$=0.75 This is due to the fact that at this $cp$ value, even strong connections between affected regions become severed (see Eq. 4, at $cp$=0.75 the reduction of ${C}_{\text{weight }}$ is equal to 1.5, thus strong connections between affected regions, that are affected twice by cp, are eliminated), which could cause a cascading effect into simulated neural activity, introducing a new regime of N1 evolution. Alterations in structural connectivity for increasing $cp$ values can be visualized in Supplementary Fig. 8.

A FC phase transition depending on a parameter describing connection probability between brain regions has been similarly identified by Stam et al. [60]. Their model neglected neuroplastic rewiring of structural connections, which did not influence the phase transition in our findings.

The determination of personalized, participant-specific BNMs allowed us to identify network configuration characteristic of each cognitive condition. Notably, BNMs of CTR subjects adopted $cp$ values that did not alter the Assortativity of the healthy brain network (at $cp$= 0).

Crucially, BNMs near the phase transition presented synaptic degeneration values proportional to cognitive decline severity (Fig. 5I-K). This underscores synaptic degeneration as one of the main culprits of cognitive decline, consistent with existing literature [5]. Lusch et al. [21] and Tuladhar et al. [19] have investigated the relationship between alterations of synaptic transmission and cognitive decline. In both studies, neurodegeneration was implemented as a rescission of synaptic weights between units of a deep neural network, leading to impairment in cognitive performance.

Taken together, our results identified several crucial connections between structural network alterations, neural activity and cognitive decline, notably: (i) FC during cognitive task execution is closely related to network Efficiency; (ii) the presence of distinct FC subgroups is determined by different levels of connectivity degeneration (see Fig. 3C); (iii) the presence of distinct evolutions of neural activity along the cognitive decline continuum is determined by the interplay of connectivity and synaptic degeneration (see Figs. 4 and 5); (iv) FC bimodality in the SCD group is caused by bimodality of connectivity degeneration values; (v) Cognitive decline severity is determined by magnitude of synaptic degeneration.

It is worth noting that this study has several limitations. First, the three groups included in the ERP analysis are unevenly composed, with the CTR group being notably smaller than the others. This imbalance limits the generalizability of our results, and needs to be addresses by testing our approach on large cohorts. Also, the MCI group is older than the other groups (See Supplementary Table 1). The age difference between the SCD and MCI groups was expected given the natural progression of the disease, which is more prominent in older individuals. In contrast, the age difference between the CTR and MCI groups resulted from the constraints of the recruitment process of healthy subjects, which led to a younger CTR group. However, we accounted for this by ensuring that relevant metrics of neural activity used in the experimental analysis and simulations presented no correlations with age of the participant (see Methods, ‘Participants age’ paragraph).

Additionally, the study primarily focuses on ERP neural markers and FC. Given that cognitive decline is a multifaced condition characterized by complex mechanisms it is crucial to investigate other relevant aspects of its progression. Furthermore, while the choice of the task is motivated by the fact that 3-CVT elicits multiple brain functions such as memory retrieval and attention and has recently been proposed as a potential MCI biomarker [28], it represents only one of many possible tasks for assessing cognitive decline [61]. Future iterations of the model should incorporate additional cognitive tasks to capture a broader range of cognitive processes.

Lastly, although the BNM provides valuable insights into the structural alterations underlying cognitive decline, it remains a simplified representation of the complex biological mechanisms involved in cognitive task execution. Future refinements should aim to improve biological realism by integrating additional neurophysiological and pathological factors. Moreover, future studies should validate these findings with larger and more balanced datasets while expanding the investigation into different cognitive states and network dynamics.

Overall, the model accurately reproduced human EEG responses during task execution in different stages of cognitive decline, linking neurodegeneration to evoked neural activity and decline in cognitive performance.