Twenty-three patients affected by neurodegenerative diseases were recruited at the IRCCS Mondino Foundation. The study was approved by the Local Ethical Committee and carried out in accordance with the Declaration of Helsinki. Written informed consent was obtained from all subjects. The protocol was approved by the Local Ethical Committee of the IRCCS Mondino Foundation. Patients underwent a complete diagnostic workup including clinical and neuropsychological assessment (see section below) MRI, and, when available, cerebrospinal fluid (CSF) biomarkers (amyloid-β and τ protein) assessment following the harmonized protocol of the RIN network [Italian Network of the Institutes (IRCCS) of Neuroscience and Neurorehabilitation] (Nigri et al., 2022). Subjects were classified into two main groups: 16 AD patients (13 females, 70 ± 8 years) and 7 FTD patients (1 female, 69 ± 5 years), further classified into distinct phenotypes. In particular, AD patients were additionally classified into: typical AD (10 subjects; Dubois et al., 2014); AD logopenic variant (2 subjects; Dubois et al., 2014); AD frontal variant (ADfv, 1 subject; Dubois et al., 2014); AD posterior cortical atrophy (ADpca, 1 subject; Dubois et al., 2014). One patient was classified as having corticobasal syndrome (CBS, 1 subject; Hassan et al., 2011), and one with dementia with Lewy bodies (DLB, 1 subject; McKeith et al., 2017). On the other hand, FTD patients were classified into: behavioral FTD (FTDbv, 5 subjects; Rascovsky et al., 2011); Primary Progressive Aphasia non-fluent variant (PPAnf, 1 subject; Gorno-Tempini et al., 2011), and Primary Progressive Aphasia semantic variant (PPAsv, 1 subject; Gorno-Tempini et al., 2011). Pharmacological therapy was also recorded.

Ten healthy controls (HC, 6 females, 67 ± 3 years) were enrolled on a voluntary basis as reference group. All HC underwent clinical assessment to exclude any cognitive impairment. For all subjects, exclusion criteria were: age >80 years, a diagnosis of significant medical, neurological and psychiatric disorder, pharmacologically treated delirium or hallucinations and secondary causes of cognitive decline (e.g., vascular metabolic, endocrine, toxic, and iatrogenic). Supplementary Table 1 shows demographic, clinical, and neuropsychological data.

All subjects underwent a neuropsychological examination based on a standardized battery of tests to assess their global cognitive status (Mini-Mental State Examination, MMSE) and different cognitive domains: memory (verbal: Rey’s Auditory Verbal Learning Test, RAVLT; visuo-spatial: Rey–Osterrieth complex figure recall), phonemic and semantic fluency, visuo-constructional abilities (Rey–Osterrieth complex figure copy), attention (Trial Making Test part A, TMT-A) and executive functions (Frontal Assessment Battery, FAB; Trial Making Test part B and B-A; Stroop color-word test interference, time and errors; Raven’s Colored Progressive Matrices, CPM47).

Raw scores were corrected for the effect of age, education, and sex according to the reference norms for the Italian population. Accordingly, corrected scores were classified into five Equivalent Scores (ES), from 0 to 4, with an ES of 0 reflecting a pathological performance, based on percentiles (Capitani and Laiacona, 1997). Domain scores, calculated by averaging the ES of the single tests, were obtained for memory, language-fluency, visuo-constructional abilities, attention, and executive functions, respectively.

All subjects underwent MRI examination using a 3T Siemens Skyra scanner with a 32-channel head coil. The MRI protocol was harmonized within the RIN network including both diffusion weighted imaging (DWI) and resting-state fMRI (rs-fMRI) (Nigri et al., 2022). For DWI data a two-shell standard single-shot echo-planar imaging sequence (EPI) [voxel size = 2.5 mm × 2.5 mm × 2.5 mm, TR/TE = 8,400/93 ms, two shells with 30 isotropically distributed diffusion-weighted directions, diffusion weightings of 1,000 and 2,000 s/mm², 7 non-diffusion weighted b = 0 s/mm² images (b₀ images) interleaved with diffusion-weighted volumes] was implemented, and 3 non-diffusion weighted images with the reversed phase-encoding acquisition were additionally acquired for distortion correction. For the rs-fMRI data, GE-EPI sequence (voxel size = 3 mm × 3 mm × 3 mm, TR/TE = 2,400/30 ms, 200 volumes) was set. For anatomical reference, the protocol included a whole brain high-resolution 3D sagittal T1-weighted (3DT1) scan (TR/TE = 2,300/2.96 ms, TI = 900 ms, flip angle = 9°, voxel size = 1 mm × 1 mm × 1 mm).

Preprocessing of diffusion and fMRI data was performed according to Monteverdi et al. (2022). Briefly, DWI data were denoised, and corrected for motion and eddy currents distortions (FMRIB Software Library and FSL)¹ (Andersson and Sotiropoulos, 2016), then white matter, gray matter (GM), subcortical GM and CSF were segmented from the co-registered 3DT1 volume (MRtrix3)² (Patenaude et al., 2011). 30 million streamlines whole-brain anatomically constrained tractography (Smith et al., 2012) was performed within MRtrix3, estimating fibers orientation distribution with multi-shell multi-tissue constrained spherical deconvolution (CSD) and using probabilistic streamline tractography (Tournier et al., 2012). fMRI preprocessing was carried out combining SPM12³, FSL and MRtrix3 commands in a custom MATLABR2019b script. Marchenko–Pastur principal component analysis (MP-PCA) denoising (Ades-Aron et al., 2020) was firstly performed, followed by slice-timing correction, realignment, co-registration to the 3DT1 volume, polynomial detrending, nuisance regression of 24 motion parameters (Friston et al., 1996) and CSF temporal signal (Muschelli et al., 2014), and temporal band-pass filtering (0.008–0.09 Hz).

An ad hoc anatomical atlas in MNI (Montreal Neurological Institute) space was created combining 93 cerebral (AAL) (including cortical/subcortical structures) and 33 cerebellar (SUIT) labels (Diedrichsen et al., 2009). We then performed a mapping between our ad hoc atlas and the Buckner and Yeo (Buckner et al., 2011; Thomas Yeo et al., 2011) cerebral and cerebellar functional atlases to select the gray matter anatomical nodes of six networks known to support specific functions: (i) integrative networks: DMN, frontoparietal network (FPN), limbic network (LN), and attention network (AN); (ii) motor and sensory networks: visual network (VN) and somatomotor network (SMN) (Figure 1). For each subject, the gray matter parcellation of our combined anatomical atlas was applied to the whole-brain tractography to extract a whole-brain structural connectivity (SC) matrix, with the normalized number of streamlines as edges and cortical/subcortical/cerebellar areas as nodes. The subset of nodes defining each network and their connections were extracted from whole-brain SC obtaining specific network SC matrices, used as input to TVB (as detailed below). In addition, both static and dynamic experimental FC (expFC and expFCD, respectively) were reconstructed from rs-fMRI data for each of the six brain networks, to capture not only synchronous fluctuations of BOLD signals but also their spatiotemporal-dynamics during resting-state (Hansen et al., 2015). The expFC matrix was created by extracting the time-course of BOLD signals for each node and computing the Pearson’s correlation coefficient (PCC) of the time-course of pairs of atlas-defined brain regions. Matrix elements were converted with a Fisher’s z transformation and thresholded at 0.1206 (Palesi et al., 2020). FCD is the dynamic representation of FC over the time and reflects time-variant changes of resting state recordings. To obtain expFCD, expFC was computed over a sliding window of 40 s (expFCsw), shifted incrementally by 1 repetition time, which for our data it means to have 178 expFCsw (Battaglia et al., 2020). Then, each expFCsw was vectorized by considering the upper triangular entries and the vectorized expFCsw were correlated with each other generating the expFCD. Thus, expFCD was calculated as a time-versus-time matrix, containing the Pearson correlation between each expFCsw and all expFCsw, centered at all other time points along the total acquisition window, quantifying, therefore, time-evolving dynamics.

The TVB workflow [reported in Monteverdi et al. (2022) for the whole brain] was applied to each one of the six selected brain networks (Figure 2). The Wong-Wang neural mass model (Deco et al., 2014; Supplementary Figure 1), implemented with an optimized C code (Schirner et al., 2022), was chosen to simulate local microcircuits activity, resulting from two populations of interconnected excitatory and inhibitory neurons coupled through NMDA and GABA receptor types. In our TVB simulations, this neural mass model was associated to each node of the network, while the SC matrix was used for the nodes interconnection. A set of parameters had to be tuned globally for each network: the global coupling (G), which is a scaling factor that represents the connections strength, and three synaptic parameters, i.e., the excitatory (NMDA) synapses (J_NMDA), the inhibitory (GABA) synapses (J_i), and the recurrent excitation (w₊). The neural activity simulated with TVB was fed into the Balloon-Windkessel hemodynamic model (Stephan et al., 2007) to reconstruct resting-state BOLD fMRI time-courses over 8 min length and compute simulated FC (simFC) and FCD (simFCD). Parameters were adjusted iteratively using expFC and expFCD of each network as targets to optimize model fitness and the validity of the result was assessed by iterating the optimization using different initial conditions (Supplementary Figure 2; Good et al., 2022). For the simFC vs. expFC comparison, model parameters were tuned until the PCC between experimental and simulated data reached the highest value. For the simFCD vs. expFCD comparison, differences between experimental and simulated FCD were assessed using the Kolmogorov–Smirnov (KS) distance: lower KS values corresponded to a lower distance of frame-by-frame FCD properties, meaning that model and experimental matrices were closest to each other. Thus, to achieve the optimal TVB simulation it was necessary to find both the highest PCC and the lowest KS values. To this aim, an overall cost function was defined as (1 − PCC) + KS and lowest cost function values implied the best fit both to static and dynamic functional data (Kong et al., 2021).

Statistical tests were performed using SPSS software version 21. Optimal TVB parameters derived for each subject and for each network were tested for normality (Shapiro–Wilk) and then two control tests were performed to assess: (i) whether different networks presented a different E/I balance within the same clinical group (i.e., evaluation of the inter-network E/I balance); and (ii) whether inter-networks E/I balance changed in healthy vs. pathological subjects. Two statistical tests were used: (i) univariate general linear model followed by bias-corrected accelerated Bootstrap (Pek et al., 2018) to correct for age and gender differences in the groups and take into account non-Gaussian data distributions; and (ii) multivariate general linear model between the mean difference (i.e., the difference between the mean value) of TVB parameters in each network compared to the other networks in different clinical groups. Then, a multiple regression analysis was performed to investigate the relationship between individual scores of the 5 cognitive domains (memory, language-fluency, visuo-constructional abilities, attention, and executive functions) and the optimal TVB parameters. Neuropsychological scores in each cognitive domain were considered as dependent variables while model parameters derived for each network were used as predictors in a backward approach. The regression algorithm automatically removed one or more predictors to identify which of them significantly (p < 0.05) explained neuropsychological scores variance.

Meaningful TVB parameters were given as an input to clustering analysis. To avoid overfitting in the study design, the clustering algorithm first performed a feature selection reducing the number of TVB parameters (i) through a semi-supervised approach using LASSO regression model with TVB parameters as independent variables and the diagnostic class as dependent variable; (ii) via PCC between the survived TVB parameters and the diagnostic class; (iii) through Variant Inflation Factors to find out just three meaningful but not correlated features. Then the number of clusters was derived using Gap statistics and the K-means algorithm was applied to label each subject into one cluster defining a personalized fingerprint (Redolfi et al., 2020).

Model optimization was performed in each of the six brain networks considered in this work. Global coupling (G) and mesoscopic network parameters (J_i, J_NMDA, and w₊) were adjusted iteratively to fit the experimental data. The reliability of the procedure was assessed by an extensive exploration of the parameter space and by iterating the optimization using different initial conditions (Supplementary Figure 2; Good et al., 2022). Model optimization yielded subject-specific sets of model parameters describing connectivity and E/I balance in each network. TVB parameters revealed differences between networks of healthy and pathological subjects (Supplementary Figure 3 and Supplementary Table 2) that will be further analyzed and explained below.

The mean difference of each network compared to the others was computed in different clinical groups for all the TVB parameters (i.e., G, J_i, J_NMDA, and w₊). Significant mean difference changes were found both for the TVB parameters in several networks (Figure 3) with network changes summarized in Figure 4. In particular, both in AD and FTD, the connectivity strength (G) decreased in LN and increased in DMN compared to HC; in FTD, G of FPN was lower with respect to other networks. Considering mesoscale synaptic parameters, both FTD and AD showed lower excitatory coupling (J_NMDA) in SMN compared to HC; in FTD, J_NMDA was lower in VN and higher in FPN; in AD, J_NMDA in DMN was higher with respect to other networks. Both in AD and FTD, recurrent excitation (w₊) increased in SMN compared to HC; in FTD, w₊ was lower in FPN; in AD w₊ was lower in DMN with respect to other networks. In FTD, inhibitory coupling (J_i) was lower in FPN and higher in DMN; in AD, AN showed higher J_i and LN lower J_i with respect to other networks.

To assess the significance of the observed mean difference changes in TVB parameters, these were used in backward regression to explain the variation of scores associated to different neuropsychological domains assessed in patients. Network-specific levels of global coupling, excitatory coupling, inhibitory coupling, and recurrent excitation (predictors) significantly (p < 0.05) explained a percentage of variance in the cognitive domains, in which the network is involved (Table 1). The explained variance ranged from ∼20 to ∼45%. Therefore, the mean difference changes in TVB parameters were relevant to explain the neuropsychological performance of patients.

The TVB parameters that significantly explained the neuropsychological performance were considered for patients’ labeling using machine learning strategies. From the nineteen parameters identified with backward regression (Table 1) the LASSO algorithm allowed to reduce them to six. Then, G of FPN was excluded, presenting PCC <0.1, and after Variant Inflation Factors three independent and not correlated variables were considered as the most informative features to perform patient’s labeling: J_i of AN, G of the LN and G of the DMN. Gap statistics identified that seven homogeneous classes would be appropriate and the K-means assigned each subject to one of the seven clusters. Each of the identified clusters was characterized by a specific composition of TVB network features (Figure 5A and Supplementary Figure 4). Considering the biophysical meaning of each parameter, they could be described as follows:

Clusters 0 and 3 were associated with the lowest mean MMSE values (20.39 ± 5.21 and 18.57 ± 8.28, respectively) while clusters 1 and 4 were associated with the highest mean values (29.08 ± 1.14 and 29.33 ± 1.16, respectively) (Figure 5B and Table 2). No HC was classified into clusters 0 or 3. Moreover, different disease phenotypes were distributed amongst the clusters (Figure 5B): typical AD subjects spread through clusters supporting a heterogeneous distribution of connectivity values in the LN and DMN networks and inhibition of the AN, but no AD patient was found in cluster 1 and the single AD patient belonging to cluster 4 presented a high MMSE score; on the other hand, cluster 0 contained the DLB phenotype, cluster 1 both the non-amnesic variants of AD (ADlv and ADpca), cluster 3 the logopenic variant and the CBS characterized by low MMSE values and cluster 5 contained the frontal variant. Considering the FTD group, FTDbv were heterogeneous and distributed amongst different clusters, but no FTDbv were found in cluster 3. On the other hand, cluster 0 contained PPAsv and cluster 5 PPAnf. Finally, pharmacological assessment of subjects belonging to different groups indicated that the majority of subjects following an antidepressant or anxiolytic treatment belonged to cluster 0 or 1 (Table 2). In particular, subjects belonging to cluster 0 were following an antidepressant therapy mainly with selective serotonin reuptake inhibitors (SSRIs), with the exception of one patient, treated with vortioxetine. Patients belonging to cluster 1, instead, were taking antidepressant drugs different from SSRIs, such as tricyclic antidepressants (e.g., amitriptyline) and serotonin-norepinephrine reuptake inhibitors (e.g., duloxetine), apart from one HC belonging to this group who was found to be on a SSRIs treatment.

In a first analysis, we compared AD and FTD for their average network model parameters. Model parameters markedly differentiated the mechanisms underlying brain networks dynamics in AD and FTD, with the most typical changes being concentrated in the DMN and LN of AD and in the FPN of FTD.

Model parameters for individual subjects were correlated with behavioral observations. Global coupling and synaptic parameters of each network significantly contributed to explain neuropsychological scores in specific cognitive domains: LN, AN, and VN with memory; DMN and VN with language-fluency; LN with attention; SMN and FPN with visuo-constructional performance; FPN with executive functions. This evidence is in line with several reports on the importance of motor regions in visuo-constructional performance (Chen et al., 2016), the contribution of AN and limbic areas in memory (Epelbaum et al., 2018), the relevance of frontoparietal areas for executive and visuo-constructional control (Melrose et al., 2013; Dixon et al., 2018), the role of DMN integration for semantic fluency (Jockwitz et al., 2017), and the involvement of visual structures in memory and language-fluency (Kucewicz et al., 2019; Vonk et al., 2019).

Thus, the relationship between neurophysiological parameters in brain networks and neuropsychological scores, which has not been investigated before, provides new cues for understanding the physiopathology of AD and FTD.

The most meaningful model biomarkers for patient’s labeling were G in DMN, G in LN, J_i in AN, consistent with known salient aspects of dementia affecting the ability of daydreaming (DMN), emotional control (LN) and attention (AN). Subjects were found to be distributed between seven different clusters revealing correspondence with their cognitive status (assessed with MMSE) and pharmacological treatment.

Patients with different MMSE scores tended to populate different clusters (see Figure 5), broadly separating patients from HC (MMSE >30), highlighting the importance of DMN, LN, and AN connectivity strength and E/I balance to ensure healthy cognitive function. Interestingly, high G between DMN nodes is associated with a worse performance, being hence disruptive and not compensatory. This analysis suggests that the heterogeneity of subject-specific TVB parameters is able to identify AD “subtypes” (Pini et al., 2021; Rauchmann et al., 2021) and FTD variants. Indeed, subjects belonging to atypical forms of AD and FTD variants were assigned to different clusters, capturing specific aspects of these pathologies and mostly mapping clinical severity assessed with MMSE. A finer grained analysis based on clinical phenotypes is not currently possible, given the limited sample size.

Patients’ labeling based on TVB parameters correlated with pharmacological treatment. Most subjects belonging to clusters 0 and 1 were on antidepressant or anxiolytic treatment (cf. Table 2), which may influence the connectivity strength and the E/I balance of cognitive networks. The effect of SSRIs on LN and DMN FC is increasingly recognized (Van Wingen et al., 2014; Li et al., 2021), while the effect of antidepressant treatment with molecules different from SSRIs, such as vortioxetine, tricyclic molecules or SNRIs (Pérez et al., 2018), as well as the influence of antidepressants on GABA and glutamate levels needs further assessment (Spurny et al., 2021). Considering that patients treated with SSRIs belong to cluster 0 while patients treated with other antidepressant classes belong to cluster 1, our results pose a very intriguing question: is there an opposite impact on cognitive networks exerted by antidepressants with different mechanisms of action or does the cognitive networks profile determine pharmacological treatment response? Future work should study TVB parameters longitudinally pre-post treatment to answer this important question with major potential clinical impact.

It should be noted that, in our cohort, patients were not treated with NMDA receptor antagonists (like memantine) (Robinson and Keating, 2006) or acetylcholinesterase inhibitors (like galantamine, rivastigmine, and donepezil) (Marucci et al., 2021), which are also known to act on AD pathophysiology. NMDA receptors are main triggers of synaptic plasticity, also affected by excitotoxicity and cholinergic receptors that, in turns, act on learning (Waxman and Lynch, 2005; Hasselmo, 2006). Since in the Wong–Wang neural mass model J_NMDA is mostly related to slow synaptic mechanisms driven by NMDA receptors (Deco et al., 2014) and receptor density can be remapped onto TVB through parameterization (Deco et al., 2021), an assessment of these receptor-dependent properties could be an important development in future studies.

The small sample size can be seen as a potential limitation in the present study. However, the main aim of this investigation was to assess the ability of TVB to provide a personalized fingerprint of patients, potentially beyond known diagnosis. TVB modeling provides a set of physiological features at single subject level, otherwise not available from standard signal/image acquisition and analysis. Thus, the small sample size does not impact on the TVB ability of uncovering subject-specific features of FC, and E/I profile. The high correlation of TVB parameters with both cognitive performance and pharmacological treatment reveals indeed its exquisite sensitivity to single-subject profiles and opens a broad range of prospective for clinical applications. On the other hand, the application of TVB to a larger cohort of patients bears the potential of improving disease classification of disease subtypes, critical for treatment stratification and for establishing intervention workflows.