Linking Molecular Pathways and Large-Scale Computational Modeling to Assess Candidate Disease Mechanisms and Pharmacodynamics in Alzheimer's Disease

Introduction: While the prevalence of neurodegenerative diseases associated with dementia such as Alzheimer's disease (AD) increases, our knowledge on the underlying mechanisms, outcome predictors, or therapeutic targets is limited. In this work, we demonstrate how computational multi-scale brain modeling links phenomena of different scales and therefore identifies potential disease mechanisms leading the way to improved diagnostics and treatment. Methods: The Virtual Brain (TVB; thevirtualbrain.org) neuroinformatics platform allows standardized large-scale structural connectivity-based simulations of whole brain dynamics. We provide proof of concept for a novel approach that quantitatively links the effects of altered molecular pathways onto neuronal population dynamics. As a novelty, we connect chemical compounds measured with positron emission tomography (PET) with neural function in TVB addressing the phenomenon of hyperexcitability in AD related to the protein amyloid beta (Abeta). We construct personalized virtual brains based on an averaged healthy connectome and individual PET derived distributions of Abeta in patients with mild cognitive impairment (MCI, N = 8) and Alzheimer's Disease (AD, N = 10) and in age-matched healthy controls (HC, N = 15) using data from ADNI-3 data base (http://adni.loni.usc.edu). In the personalized virtual brains, individual Abeta burden modulates regional Excitation-Inhibition balance, leading to local hyperexcitation with high Abeta loads. We analyze simulated regional neural activity and electroencephalograms (EEG). Results: Known empirical alterations of EEG in patients with AD compared to HCs were reproduced by simulations. The virtual AD group showed slower frequencies in simulated local field potentials and EEG compared to MCI and HC groups. The heterogeneity of the Abeta load is crucial for the virtual EEG slowing which is absent for control models with homogeneous Abeta distributions. Slowing phenomena primarily affect the network hubs, independent of the spatial distribution of Abeta. Modeling the N-methyl-D-aspartate (NMDA) receptor antagonism of memantine in local population models, reveals potential functional reversibility of the observed large-scale alterations (reflected by EEG slowing) in virtual AD brains. Discussion: We demonstrate how TVB enables the simulation of systems effects caused by pathogenetic molecular candidate mechanisms in human virtual brains.


INTRODUCTION
Neurodegenerative diseases (NDD) gain increasing socioeconomic relevance due to an aging society (WHO, 2011;Wimo et al., 2011Wimo et al., , 2017Xu et al., 2017). The Alzheimer's Association's latest report estimates the yearly cost of Alzheimer's disease (AD) treatment in the U.S. at $277 billion (Alzheimer's Association, 2018). By 2050 this number is expected to rise as high as $1.1 trillion. According to the report, early diagnosis could save up to $7.9 trillion in cumulated medical and care costs by the year 2050. While the prevalence of AD-the most common cause of dementia and the most common NDD in general-increases, its cause is still not understood, nor is there a cure. Our understanding of their pathogenesis and classification remain insufficient. Therefore, we aim to integrate clinical data from molecular biology and neurology, using nonlinear systems theory. Our aim is to build predictive models for health-outcome and cognitive function by individual virtual brain simulations using The Virtual Brain (TVB; thevirtualbrain.org) platform (Ritter et al., 2013;Sanz Leon et al., 2013). TVB integrates various empirical data in computational models of the brain that allow for the identification of neurobiological processes that are more directly linked to the causal disease mechanisms than the measured empirical data. Biomedical sciences are currently lacking a mapping between the degree and facets of cognitive impairments, biomarkers from high-throughput technologies, and the underlying causal origins of NDD like AD. The imperative for the field is to identify the features of brain network function in NDD that predict whether a person will develop dementia. The heterogeneity of NDD makes it difficult to develop robust predictions of cognitive decline. This can be addressed by large prospective studies where there is potential for participants to develop NDD. It is difficult in general to predict individual disease progression and this is a particular challenge in complex nonlinear systems, like the brain, where emergent features at one level of organization (e.g., cognitive function) can come about through the complex interaction of subordinate features (e.g., network dynamics, molecular pathways, gene expression). The Virtual Brain takes into account the principles of complex adaptive systems and hence poses a promising tool for identifying mechanistic predictive biomarkers for NDD. Due to the high dimensionality of brain models and the even greater complexity of the to-be-simulated brain states, selecting the used modeling approach carefully for a specific question of interest is essential.
The candidate biological mechanism under investigation in the present study is related to amyloid beta (Abeta), a protein that is an oligomeric cleavage product of the physiological amyloid precursor protein (APP) (Bloom, 2014;Selkoe and Hardy, 2016). The soluble oligomers have the tendency for polymerization (Sadigh-Eteghad et al., 2015;Selkoe and Hardy, 2016). Due to their non-physiological configuration they aggregate and accumulates in brain tissue-a process that starts already in early preclinical stages of AD, i.e., many years before the onset of symptoms-typically in the fifth decade of life (Braak and Braak, 1997)-as shown in rodent models (Busche et al., 2012) and human studies (Klunk et al., 2007;Jack et al., 2009). Aggregated Abeta and its intermediates, soluble Abeta oligomers, can act directly neurotoxic (Hardy and Selkoe, 2002;Prasansuklab and Tencomnao, 2013;Selkoe and Hardy, 2016) and have been found intra-or extra-cellularly (Hardy and Selkoe, 2002;Walsh and Selkoe, 2007;Selkoe and Hardy, 2016). Those findings led to the hypothesis that the deposition of Abeta poses an initial step in the pathology of AD while Abeta has been suggested as a key feature in the pathogenesis of AD leading to major changes in the functionality and structure of the brain (Klunk et al., 2007;Jack et al., 2009;Villemagne et al., 2009). The goal of the present study is to incorporate the hypothesized qualitative and quantitative effects of Abeta on neuronal population dynamics into our brain network models, i.e., adding mathematical models that describe how molecular changes alter population activityso called cause-and-effect models. We will focus here on the disrupted inhibitory function of interneurons and consecutive hyperexcitability caused by Abeta-while we are aware of various other factors with potential roles for AD etiology, such as vascular changes (Love and Miners, 2016;Storck and Pietrzik Claus, 2018;Bannai et al., 2019), neuroinflammation (Heneka et al., 2015a,b;Wang and Colonna, 2019;Zhou et al., 2019), genetics (Mahley, 2016;Hudry et al., 2019;Takatori et al., 2019), environmental factors (Alonso et al., 2018;McLachlan et al., 2019) and concomitant proteinopathies others than Abeta pathology (Robinson et al., 2018a,b). Beside Abeta there is a second molecular hallmark associated with the pathogenesis of AD: the phosphorylated Tau "tubulin-associated unit" protein (Bloom, 2014;Guo et al., 2017;Tapia-Rojas et al., 2019) which contributes to microtubule stability in the neural cytoskeleton (Guo et al., 2017). One major argument in favor of the more prominent involvement of Abeta in the pathogenesis of AD, in contrast to Tau, is its higher specificity to AD and its appearance in the early familial variants of AD, where the molecular pathway is better understood (Blennow et al., 2006;Klunk et al., 2007;Villemagne et al., 2009). Therefore, most therapeutic strategies in the past targeted Abeta. Yet recently three clinical trials with antibodies against Abeta had to be terminated in phase III: aducanumab (Biogen, 2019;Chiao et al., 2019), crenezumab (Salloway et al., 2018;Roche, 2019), and solanezumab (Doggrell, 2018;Honig et al., 2018) did not meet the expectations to act in a disease-modifying manner slowing down the cognitive decline (Selkoe and Hardy, 2016). Nevertheless, there are still studies ongoing, e.g., with BAN-2401, an antibody against soluble monomeric Abeta protofribrils (as aducanumab) (Logovinsky et al., 2016;Osswald, 2018;Panza et al., 2019). A relevant percentage of clinically diagnosed AD patients show additional brain pathologies beside Abeta and Tau in autopsy (Robinson et al., 2018a). Even in the cases of neuropathological AD diagnosis (i.e., secured Abeta and Tau pathology in histology), 55% of cases also exhibited a pathology of alpha synuclein (which we would expect in synucleinopathies like Parkinson's disease) and 40% showed transactive response DNA binding protein 43kDa (TDP-43), a protein which we would expect in frontotemporal dementia or amyotrophic lateral sclerosis (Robinson et al., 2018b). Brain tissue of people who did not had relevant neurodegenerative brain changes in histological exams after death were showing Abeta in 50% and Tau pathology in 93% of the cases when using sensitive immunohistochemistry methods (Robinson et al., 2018b). Although Abeta and Tau are widely accepted as involved parts in the pathogenesis of AD and also define the disease entity (Jack et al., 2018), it remains unclear if they might be only epiphenomena of other contributing factors. This study hypothesizes a mechanistic role of Abeta in the disease process and builds a link between the molecular pathway alteration that leads to Abeta phenomenon of disinhibition and neural slowing in EEG (Figure 1). Our mechanistic modeling approach can help to understand the complex inter-dependencies between the involved factors in AD and will improve through iterative refinement.
Near Abeta plaques, a shift in neural activity has been observed (Busche et al., 2008). In AD mouse models with overexpression of APP and Presenilin-1, the number of hyperactive neurons was increased near Abeta plaques. This shift in the neuronal activity was associated with decreased performances in memory tests. Neuronal hyperactivity could be reduced by GABA agonists, suggesting pathology due to impaired inhibition. In neocortical and dentate gyri, pyramidal cells have been found to increase network excitability in vivo in an AD mouse model with overexpression of Abeta, that led to membrane depolarization and increased firing rates. A study by Hazra et al. (2013) investigated an AD mouse model by stimulation of the perforant pathway. AD mice showed increased amplitude and larger spatial distribution of response after stimulation. The reason for this increased network excitability was due to impaired inhibitory neuron function, i.e., the inhibitory neurons of the molecular layer of the dentate gyrus in hippocampus were in part unable to produce action potentials, which resulted in a slower postsynaptic firing rate. Ulrich (2015) added Abeta to layer V pyramidal cells of rats. In their experiments they could show a decline in inhibitory postsynaptic currents (IPSCs), attributed to postsynaptic GABA A receptor endocytosis after Abeta application. In a recent study by Ren et al. (2018) Abeta was found to increase excitability of pyramidal cells in the anterior cingulate cortex of mouse brain. The reason for hyperexcitability was again due to disturbed inhibitory input. Abeta seems to interact with the dopaminergic D1 receptor system. The D1 receptor regulates GABA release in fast-spiking (FS) inhibitory interneurons. By adding a D1 receptor antagonist to the cells they could reverse the effect of Abeta, increase IPSCs and decrease pyramidal excitability whereas D1 agonists had similar effects as Abeta. The underlying working model is that Abeta leads to dopamine release in dopaminergic neurons that activates D1 receptors at FS inhibitory interneurons and thus inhibits GABA release. As a consequence, the amplitude, frequency and total number of IPSPs is decreased. The instantaneous decrement of postsynaptic amplitude and frequency is also known as a toxic effect of Abeta in the glutamatergic system (Ripoli et al., 2014). Hence for the present modeling approach we decided to implement this Abeta dependent impaired inhibitory function. From the literature above, potential models for this disinhibiton could be either a lower IPSP amplitude or a lower firing rate or a combination thereof.
One already established drug that assesses the pathology of hyperexcitation is memantine, an N-methyl-D-aspartate (NMDA) antagonist. Memantine is recommended for the symptomatic treatment of severe AD as a mono-and combination therapy with cholinesterase inhibitors and should be also considered as possible treatment in moderate AD in the current version of the UK National Institute for Health and Care Excellence (NICE) guidelines of dementia management (Pink et al., 2018). However, normally it is considered as an alternative or addition to cholinesterase inhibitors (Pink et al., 2018). In contrast, memantine has shown in a current metaanalysis its efficacy to improve cognitive function and reduce behavioral disturbances in AD patients compared to placebo (Kishi et al., 2017). The effect was particularly caused by the moderate-to-severe AD patients (Chen et al., 2017;Kishi et al., 2017) and was also observable in combination therapies with acetyl cholinesterase inhibitors, with a significant superiority for the combination of memantine and donepezil compared to any cholinesterase monotherapies (Kishi et al., 2017). It therefore is also recommended as possible first-line therapy in AD (Kishi et al., 2017). In our study, we will evaluate "virtual memantine" interacting with the Abeta-derived hyperexcitation.
Changes in electroencephalography (EEG) are described in AD as a general and progressive slowing of brain oscillations. In AD, cognitive decline and 18 F-fluorodeoxyglucose (FDG) PET signal decreases are linked with increased left temporal power in the delta and the theta frequency bands, whereas temporo-parieto-occipital alpha band coherence decreases and delta coherence increases (Loewenstein et al., 1989;Rice et al., FIGURE 1 | Biology-infered cause-and-effect model: alteration of the molecular Abeta pathway in AD cause hyperexcitation in the neural mass model. An altered pathway from soluble Abeta monomers to oligomers to insoluble plaques leads to potentially neurotoxic Abeta accumulation (Hardy and Selkoe, 2002;Prasansuklab and Tencomnao, 2013;Selkoe and Hardy, 2016) that can be quantified by PET. Region specific Abeta burden leads to disinhibition in the neural mass model (Busche et al., 2008;Hazra et al., 2013;Ripoli et al., 2014;Ulrich, 2015;Ren et al., 2018)-thus building a bridge between molecular pathways and brain network modeling. For the evaluation of the used mathematical model, see the discussion section and Figures 12, 13. Parts of the figure are modified from Deco et al. (2017Deco et al. ( ). 1990Malek et al., 2017). Moreover, the spatial appearance of slow rhythms and hypometabolism in FDG PET have been linked (Dierks et al., 2000;Babiloni et al., 2016). A recent study produced similar findings in magnetoencephalography (MEG): A global increase of theta and a frontal increase of delta were correlated with entorhinal atrophy and glucose hypometabolism (Nakamura et al., 2018). In summary, a global slowing has been reported for AD, in particular a shift from alpha to theta and delta activity (Loewenstein et al., 1989;Rice et al., 1990;Dierks et al., 2000;Babiloni et al., 2016;Malek et al., 2017;Nakamura et al., 2018).
As a consequence of these findings, we will focus in our modeling approach on three main aspects of AD: 1. Spatial heterogeneous Abeta distribution in the brain 2. Hyperexcitation caused by impaired inhibitory function 3. Slowing of neural frequencies.
For Abeta, we propose a change in local neuronal excitability. Therefore, we construct a model of a healthy "standard brain" with an averaged structural connectivity (SC) with inferred micro-scale characteristics of excitation in those areas where a deposition of Abeta is found. We will infer this information about the local distribution of Abeta from individual AV-45 (florbetapir) positron emission tomography (PET) images. AV-45 is a PET tracer which binds to Abeta (Clark et al., 2011;Ossenkoppele et al., 2015;Morris et al., 2016;Martinez et al., 2017). Ante-mortem Abeta PET imaging can be related to post mortem Abeta pathology in brain tissue (Murray et al., 2015), corresponding to the THAL phases of Amyloid deposition (Thal et al., 2002)-as well as Tau PET (Schöll et al., 2016) can be related to the BRAAK stages of neurofibrillary tangles Braak, 1991, 1997;Braak et al., 2006). This has led to updated diagnostic criteria for AD, wherein Abeta and Tau PET can be used equivalently to neuropathology for AD diagnosis (Jack et al., 2018).
We investigate three clinical diagnostic groups of age-and gender-matched healthy controls (HC), individuals with mild cognitive impairment (MCI) and AD patients [see method section Alzheimer's Disease Neuroimaging Initiative (ADNI) Database and Table 1]. For the simulated EEG and the underlying local neural activity frequency we expect a slowing in rhythms and particular a shift from alpha to theta activity with disease progression. Finally, we will simulate the effect of an antiexcitotoxic drug, the NMDA antagonist memantine for which we expect a reversal of the observed EEG slowing.
We will in the following provide an overview of the fundamentals of the here employed brain simulation technique. The particular strength of computational connectomics (Ritter et al., 2013;Kringelbach et al., 2015;Deco et al., 2017) or brain network modeling (BNM) is to unite various kinds of information in a single biophysically plausible framework (Breakspear, 2017). BNM are typically structurally informed (or constrained) by (a) geometric information of the brain, e.g., via T1 magnetic resonance imaging (MRI), and (b) the structural connectivity (SC) derived from the tractography of diffusion MRI that is supposed to represent the white matter fiber tracts (Jirsa et al., 2002;Schirner et al., 2015). The static three-dimensional scaffold of the brain is brought to life through It is a subset of the suitable ADNI-3 participants, that had 3T imaging and all necessary image modalities. Only data from Siemens scanners was used (because this was the biggest subset of scanners).
the implementation of mathematical models, which generate activity at each brain region or node of the network, the so-called neural masses or population models (Spiegler and Jirsa, 2013;Sanz-Leon et al., 2015;Cabral et al., 2017). Population models are reduced descriptions of microscopically detailed neuronal networks (Wilson and Cowan, 1972;Zetterberg et al., 1978;Hindmarsh and Rose, 1984;Jansen and Rit, 1995;Wong and Wang, 2006;Stefanescu and Jirsa, 2008;Sanz-Leon et al., 2015)inferred for example with methods of mean field theory (Deco et al., 2008;Jirsa, 2009;Bojak et al., 2010). They describe the so called meso-scale of the brain (Deco et al., 2008;Wright and Liley, 2010), i.e., population activity as captured with imaging methods like EEG, MEG and fMRI. Some neural mass models (NMM) are linked to (and still reflect to a certain degree) neurophysiological processes at the microscopic scale while others mathematically describe the observed lumped biological behavior not differentiating between underlying neurophysiological processes (phenomenological models). Time delays in the interaction between nodes (Jirsa and Kelso, 2000;Jirsa et al., 2002;Spiegler and Jirsa, 2013;Sanz-Leon et al., 2015) are critical for the spatiotemporal organization of the evolving activity patterns in the brain (Petkoski et al., 2016(Petkoski et al., , 2018. Measured functional brain data such as EEG, MEG or functional MRI (fMRI) are used to tune the mathematical models-i.e., to fit selected free parameters of the model-to faithfully reproduce selected empirical features (Honey et al., 2007;Ghosh et al., 2008;Sotero and Trujillo-Barreto, 2008;Bojak et al., 2010;Jirsa et al., 2010;Ritter et al., 2013;Sanz-Leon et al., 2015;Kunze et al., 2016). By performing a systematic model parameter exploration, using e.g., brute force exhaustive parameter space searches, Monte-Carlo methods or weighted optimization algorithms, we can identify the optimal parameter configuration to portray the empirical functional phenomena. Thereby, we obtain indices of the brains individual function in relation to the explored parameters. This approach opens various possibilities to not only describe dependencies (i.e., correlations), but to make statements about potential underlying causal processes, i.e., mechanisms.
In this study we used TVB, an open source neuroinformatics platform (Ritter et al., 2013;Sanz Leon et al., 2013;Sanz-Leon et al., 2015;Stefanovski et al., 2016) (www.thevirtualbrain.org) for large-scale BNM simulations. We have already established the software TVB, and applied it to normative datasets, stroke, epilepsy, brain tumors, and neurodegenerative disease. For example, in stroke recovery, TVB models of patients were built using the patient's structural neuroimaging data, and the dynamics of local populations were tuned to fit the patient's functional neuroimaging data (Falcon et al., 2015(Falcon et al., , 2016. The obtained parameters for excitatory/inhibitory (EI) balance of local neuronal populations predicted the patient's response to rehabilitation up to 1 year after therapy. Our work on epilepsy was able to infer seizure propagation with a model based on the patient's own diffusion weighted MRI and stereotaxic EEG Proix et al., 2017). Moreover, positive surgical outcome was strongly associated with the epileptogenic zone that was excised as predicted by the patient's TVB model. Previous work with AD patients (n = 16), controls (n = 73), and persons with amnestic MCI (n = 35), all from the Sydney Memory and Aging Study, confirms the benefit of using the model parameters to characterize cognitive status (Zimmermann et al., 2018).
TVB provides several types of NMMs. In the present study, we selected a NMM that can simulate EEG and enables us to implement disinhibition. The wiring pattern of cortical circuitry is characterized by recurrent excitatory and inhibitory loops, and by bidirectional sparse excitatory connections at the largescale (Schüz and Braitenberg, 2002). Several NMMs therefore feature projection neurons aka pyramidal cells with long axons projecting to distant cortical regions and local excitatory and inhibitory feedbacks (Lopes da Silva et al., 1974;Freeman, 1975;Jansen and Rit, 1995). The NMM by Jansen-Rit comprises an elementary circuit of three interconnected NMMs (Figure 2) describing a cortical area (or column). It has been used to explain both epilepsy-like brain activity (Wendling et al., 2000(Wendling et al., , 2002 and various narrow band oscillations ranging from the delta to the gamma frequency bands (David and Friston, 2003) including intracranial EEG (Spiegler and Jirsa, 2013). The Jansen-Rit model has been explored extensively on a single population level (Wendling et al., 2002;David and Friston, 2003;Spiegler et al., 2011) and in BNMs (Sotero et al., 2007;Merlet et al., 2013;Kunze et al., 2016). The Jansen-Rit model has a rich dynamic repertoire, which was extensively described before (Spiegler et al., 2010).
Specifically we chose the Jansen-Rit model for the present study due to the following considerations: (1) The Jansen-Rit model comprises three interacting neural masses (representing different cellular populations) in each local circuitry: pyramidal cells, inhibitory, and excitatory interneurons ( Figure 2B). This is unique and opens the possibility to simultaneously model disinhibition, i.e., an impairment of the inhibitory neuronal subpopulation in one neural mass, and an anti-NMDAergic effect, i.e., a downscaled . This is intended to lead to an increased activity and higher output of the pyramidal cell population. The excitatory impulse response function (IRF) is specified as h e (t) = tH e exp(-t/τ e )/τ e , the inhibitory IRF is specified as h i (t, β) = tH i exp(-t/τ i (β))/τ i (β) (Equations 1, 2). These IRFs can be translated into second-order ordinary differential equations, see Equations 3-5. For explanation of the used variables, see Table 2 [figure modified from Spiegler et al. (2010)]. (D) Virtual EEG as the simulation output (projection of oscillating membrane potentials to the scalp surface) reveals a shift from alpha to theta activity in AD participants. Shown is a 5 second period of exemplary EEG channel at location T7 in participant 21 (HC, above) and 4 (AD, below). The ordinate is showing the dimensionless correlate for electric potential Φ.
The exemplary timeseries shows a typical simulation result in the study: in the alpha mode, which was the starting point of the Jansen-Rit model without the effect of Abeta, it produces monomorphic alpha activity with amplitude modulations (above). Mainly exclusively in the AD virtual brains a much more irregular theta rhythm appears (below). 2 | Used parameters for each Jansen-Rit element in the large-scale brain network (Jansen and Rit, 1995 transmission from excitatory interneurons to pyramidal cells, at the same time. (2) The ratio of excitatory and inhibitory time constants τ e /τ i in the Jansen-Rit model is suitable to model the effect of Abeta on the inhibitory interneurons (by affecting the transmission from inhibitory interneurons to pyramidal cells, Figures 2B,C) and is also known to have an effect on the simulated neural frequency (Wendling et al., 2002;Spiegler et al., 2010). Because oscillations emerge in the Jansen-Rit model of a brain region due to the interplay of positive and negative feedback loops (excitatory and inhibitory interneurons), a change in one of the time constants does not necessarily slow down or speed up rhythms. However, if both excitatory and inhibitory time constants, τ e and τ i are scaled simultaneously and uniformly, the local equilibrium of interaction between the neural masses remains the same but the time signature such as frequency changes [see Figure 9 in Spiegler et al. (2010), and see Chapter 3.1.3 "Model Equivalence" and Chapter 3.1.4 "Normalization" in Spiegler (2011)]. To conclude, higher τ i does not necessarily lead to slower rhythms and vice versa. (3) Jansen-Rit can simulate physiological rhythms observable in local field potentials (intracranially), stereo-EEG (sEEG), scalp EEG, and MEG (Jansen and Rit, 1995;Spiegler et al., 2010;Sanz-Leon et al., 2015).
Our hypothesized effect of local Abeta deposition as inferred from subject-specific AV-45 PET is a decrease of local inhibition (Busche et al., 2008;Grienberger et al., 2012;Limon et al., 2012;Verret et al., 2012;Hazra et al., 2013;Ripoli et al., 2014;Ren et al., 2018), which leads to a relatively stronger local excitation. This theory allows us translation of the Abeta distribution into the altered dynamics of a population model (Equation 14 and Figure 3). We use an averaged healthy SC to control the effect of individual differences in the connectome. I.e., in our simulations the distribution of Abeta is the only individual factor and can therefore be seen as the cause of any differences between the participants. The hypothesized microscale (synaptic), spatially distributed effect is assumed to develop an effect at the population (mesoscale) level and to eventually propagate to the large-scale of the whole brain. A schematic illustration of this concept is provided in Figures 1, 2.
The ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner. The primary goal of ADNI has been to test whether serial MRI, PET, other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of MCI and early AD. For up-to-date information, see http://www.adni-info.org. In the presently ongoing trial, ADNI-3, the measurements contain T1, T2, DTI, fMRI, Tau PET, Abeta PET, and FDG PET for the participants. The total population of ADNI-3 will contain data of about 2,000 participants (comprising AD, MCI, and HC, see http://adni.loni.usc.edu/adni-3/). As inclusion criterion for AD patients the diagnosis criteria of NINCDS-ADRDA from 1984 were used, which contains only clinical features (McKhann et al., 1984). Inclusion criteria for both HC and MCI were a Mini Mental State Examination (MMSE) score between 24 and 30 as well as age between 55 and 90 years. For MCI in addition, the participant must have a subjective memory complaint and abnormal results in another neuropsychological memory test. To fulfill the criteria for AD, the MMSE score had to be below 24 and the NINCDS-ADRDA criteria for probable AD had to be fulfilled (McKhann et al., 1984). Imaging and biomarkers were not used for the diagnosis. For the full inclusion criteria of ADNI-3 see the study protocol (page 11f in http://adni.loni.usc.edu/wp-content/themes/freshnewsdev-v2/documents/clinical/ADNI3_Protocol.pdf). An overview of the epidemiological characteristics of the participants included in this study can be found in Table 1.

Structural MRI
We calculated an individual brain parcellation for each included participant of ADNI-3. We followed the minimal preprocessing pipeline (Glasser et al., 2013) of the Human Connectome Project (HCP) for our structural data using Freesurfer (Reuter et al., 2012) (https://surfer.nmr.mgh.harvard. edu/fswiki/FreeSurferMethodsCitation), FSL (Smith et al., 2004;Woolrich et al., 2009;Jenkinson et al., 2012) and connectome workbench (https://www.humanconnectome.org/software/ connectome-workbench). Therefore, we used T1 MPRAGE, FLAIR and fieldmaps for the anatomical parcellation and DWI for tractography. This consists of a Prefreesurfer, Freesurfer, and Postfreesurfer part. We skipped the step of gradient non-linearity correction, since images provided by ADNI already are corrected for this artifact. Also, the MNI templates were used at 1 mm resolution instead of 0.7 mm. In the Freesurfer pipeline we skipped the step of downsampling our data from 0.7 to 1 mm 3 , and all recon-all and intermediate steps were performed with the original image resolution. We then registered the subject cortical surfaces (32 000 vertices) to the cortical parcellation of Glasser et al. (2016) using the multimodal surface matching (MSM, see Robinson et al., 2014) tool. For the registration we used cortical thickness, MyelinMaps, cortical curvature and sulc from the subject and template surface. We mapped the parcellation on the surface back into the gray matter volume with connectome workbench. This volume parcellation surfed as the mask for the connectome and PET intensity extraction.

PET Images
We used the preprocessed version of AV-45 PET. These images had following preprocessing already performed by ADNI: Images acquired 30-50 min post tracer injections: four 5-min frames (i.e., 30-35 min, 35-40 min...). These frames are co-registered to the first and then averaged. The averaged image was linearly aligned such that the anterior-posterior axis of the subject is parallel to the AC-PC line. This standard image has a resolution of 1.5 mm cubic voxels and matrix size of 160 · 160 · 96. Voxel intensities were normalized so that the average voxel intensity was 1. Finally, the images were smoothed using a scannerspecific filter function. The filter functions were determined in the certification process of ADNI from a PET phantom. We used the resulting image and applied the following steps: Rigid aligning the PET image to participants T1 image (after being processed in the HCP structural pipeline). The linear registration was done with FLIRT (FSL). The PET image was than masked with the subject specific brainmask derived from the structural preprocessing pipeline (HCP). To obtain the local burden of Abeta, we calculated the relative intensity to the cerebellum as a common method in the interpretation of AV-45-PET, because it is known that the cerebellum does not show relevant AV-45 PET signals and can therefore act as a reference region for interindividual comparability between patients (Clark et al., 2011;Meyer et al., 2018). The intensity of gamma radiation, which is caused by a neutralization reaction between local electrons and the emitted positrons of the nuclear tracer is measured for each voxel in the PET image and divided to the cerebellar reference volume: the standardized uptake value ratio (SUVR). We therefore receive in each voxel a relative Abeta burden β which is aggregated according to the parcellation used for our present modeling approach (see below). Thus, we obtain a value β a for the Abeta burden in each brain region a. The cerebellar white matter mask was taken from the Freesurfer segmentation (HCP structural preprocessing). The image was then partial volume corrected using the Müller-Gärtner method from the PETPVC toolbox (Thomas et al., 2016). For this step the gray (GM) and white matter segmentation from Freesurfer (HCP structural preprocessing) was used. Subcortical region PET loads were defined as the average SUVR in subcortical GM. Cortical GM PET intensities were mapped onto the individual cortical surfaces using connectome workbench tool with the pial and white matter surfaces as ribbon constraints. Using the multimodal parcellation from Glasser et al. (2016) we derived average regional PET loads.

DWI
We calculated individual tractography only for included HC participants of ADNI-3 to average them to a standard brain template (see section Virtual Human Standard Brain Template Out of Averaged Healthy Brains below). Preprocessing of the diffusion weighted images was mainly done with the programs and scripts provided by the MRtrix3 software package (http:// www.mrtrix.org).
The following steps were performed: Dwidenoise. Denoising the DWI data using the method described in Veraart et al. (2016).
Diw2mask. brainmask estimation from the DWI images. Dwiintensitynorm. DWI intensity normalization for the group of participants.
Dwi2response. The normalized DWI image was used to generate a WM response function. We used the algorithm described by Tournier et al. (2013) in this step.
Average_response. An average response function was created from all participants. Dwi2fod. Using the spherical deconvolution method described by Tournier et al. (2007) we estimated the fiber orientation distribution using the subject normalized DWI image and the average response function. From the DWI data a mean-b0 image was extracted and linear registered to the T1 image. The inverse of the transform was used to bring the T1 brain masked and aparc+aseg image (from HCP structural preprocessing) into DWI space. The transformed aparc+aseg image was used to generate a five tissue type image.
Tckgen. Anatomical constrained tractography  was performed using the iFOD2 algorithm (Tournier et al., 2010). Tracks in the resulting image were weighted using SIFT2 algorithm (Smith et al., 2015). We mapped the registered parcellation from Glasser back into the volume. The cortical and subcortical regions than were used to merge the tracks into a connectome.

EEG Forward Solution in TVB
After structural preprocessing with the HCP pipeline we used the individual cortical surfaces and T1 images to compute the person specific Boundary Element Models in Brainstorm (Tadel et al., 2011). Scalp, outer, and inner skull were modeled with 1922 vertices per layer. Using the default "BrainProducts EasyCap 65" EEG cap as locations for the signal space and the cortical surface vertices as source space. The leadfield matrix was estimated using the adjoint method in OpenMEEG with the default conductivities 1, 0.0125 and 1 for scalp, skull and brain, respectively. Because we are performing region-based simulations only (i.e., no vertexwise modeling), the leadfield matrix was simplified by summing the coefficients of vertices that belong to the same region. EEG signal was generated by matrix multiplication of the neural time series with the lead field matrix.

Virtual Human Standard Brain Template Out of Averaged Healthy Brains
We use the SCs of all ADNI-3 participants of the group HC, derived from the diffusion-weighted and structural MRI, to average them to one connectome matrix. Two of the HC participants included in the average template were excluded for simulations because it was impossible to compute their leadfield matrices for EEG calculation. Therefore, we use an arithmetic mean C µ = ( n i =1 C i )/n = (C 1 + C 2 + . . . + C n )/n, wherein C µ is the averaged SC matrix, n is the number of HC participants and C i is the individual SC matrix. The SC matrix and the organization of the corresponding graph can be found in Figure 4. As it can be seen in Figure 4B, general characteristics of physiological SCs as symmetry, laterality, homology, and subcortical hubs are maintained in the averaged connectome. By choosing an averaged SC instead of individual SCs, it was possible to control all factors except of the individual Abeta distribution supporting our intention to compare the simulated activity that resulted from a "pathogenic" modification by Abeta.

Cause-and-Effect Model of Abeta in the Jansen-Rit Model
The dynamics of the Jansen-Rit model show a rich parameter dependent behavior (Spiegler et al., 2010). A bifurcation analysis of the single population Jansen-Rit model (in contrast to network embedded interacting populations) catalogs and summarizes the repertoire of the model. Bifurcation here refers to a qualitative change in the system behavior with respect to parameter changes. Qualitative changes can be for instance the shift from waxing and waning alpha rhythm as observed in resting human brains to spike wave discharges as observed during epileptic seizures. Bifurcation diagrams explore the qualitatively different states [divided by bifurcations, see Supplementary Figure 1, from Spiegler et al. (2010)]. The bifurcation analysis revealed an important feature of the Jansen-Rit model, which is bistability, that is, the coexistence of two stable states for a certain parameter range (i.e., regime). The bistable regime allows the coexistence of two self-sustained oscillatory states for the standard parameter configuration (Jansen and Rit, 1995) and Table 2 of which one state generates rhythmic activity in the alpha band and the other one produces slower big spike-wave complexes in theta rhythm. Changes in the kinetics of excitatory and inhibitory PSPs (i.e., changes of time constants) change the model behavior in a way which makes it suitable to scale, that is, to speed up or to slow down local dynamics (Spiegler et al., 2010)-and therefore to scale the global frequency, too. The results of the systematic parameter exploration of the excitatory and inhibitory time constants are summarized in Supplementary Figure 2. For our study, we constructed the model to portrait a wide range of physiological neural rhythms by using a fast limit cycle with alpha and beta frequencies and a slow limit cycle with theta and delta frequencies. To achieve this dynamic behavior of two limit cycles, we used first a very low input on the pyramidal cells (firing rate 0.1085/ms) and no input on the inhibitory interneurons to not overlay the Abeta effects we introduce here. Therefore in the "healthy" condition without the effect of Abeta, the system operates near the subcritical Andronov-Hopf and the saddlesaddle bifurcations (leftmost region in Supplementary Figure 1). For the time constants, we used the area of alpha rhythm (blue area in Supplementary Figure 2) as control condition without any effect of Abeta. The detailed parameter settings can be found in Table 2.
The information about the local Abeta burden is derived from the individual AV-45 PET. As there exists no established clinical standard for SUVR cut-off thresholds differentiating normal form pathological Abeta loads. To scale the possible neurotoxic effect in a realistic way, we need to approximate at what point Abeta toxicity occurs. Following the literature, a 96% correlation to autopsy after Abeta PET was achieved via visual assessment of PET images. The corresponding SUVR cutoff was 1.2 (Clark et al., 2011). Another study showed a higher cut-off point at SUVR ≥ 1.4 for a 90% sensitivity of clinically diagnosed AD patients with an abnormal Abeta PET scan (Jack et al., 2014). We use here the higher cut-off threshold of SUVR ≥1.4. Consequently, we propose a cause-and-effect model for Abeta that is mapping molecular changes to computational brain network models: The inhibitory time constant τ i in each point is a function of β a . The higher Abeta SUVR, the higher τ i and therefore the filter action for the synaptic transmission is slower. We decided for this implementation via a synaptic filter slowing because of several reasons: 1. We are focusing on disease linked alterations of EEG frequencies. Hence, we intended to assess a model feature that is already known to be frequency-effective, i.e., it can vary resulting simulated EEG frequencies. From former explorations of the Jansen-Rit-model we know that the neural frequencies are influenced by the ratio of excitatory and inhibitory time constants (Spiegler et al., 2010). 2. Cellular studies are supporting the hypothesis of altered inhibition as a cause for hyperexcitation (Hazra et al., 2013;Ripoli et al., 2014;Ren et al., 2018)-hence we decide for mapping Abeta on the inhibitory time constant leading to a disturbed Excitatory-Inhibitory (E/I) balance. 3. By using a time-effective feature, we intended to differentiate the micro-scale neurotoxic effect of Abeta on synaptic level (Ripoli et al., 2014;Ulrich, 2015;Ren et al., 2018) from connectivity-effective phenomena on a larger scale, which could e.g., be modeled by an alteration of connection strength.
A detailed exploration of the effects that we introduce by this model can be found in the discussion section. We develop a transform function to implement the PET SUVR in parameters of the brain network model. Specifically, we postulate a sigmoidal decrease function that modifies the  (Glasser et al., 2016), 9 left subcortical regions, 9 right subcortical regions, 1 brainstem region. It gets obvious the difference between interhemispherical commissural fibers (lower weights, with a slightly pronounced diagonal between homologous regions) and intrahemispherical association fibers (higher weights). Moreover, we can observe the strong connection pattern of the subcortical areas (above region 360). (B) Graph of the underlying SC. As a threshold, only the strongest 5% of connections were kept for binary transformation to the adjacency matrix for the graph. Node positions are derived from the inner structure of the graph by a "force" method (Fruchterman and Reingold, 1991), assuming stronger forces and therefore smaller distances between tightly connected nodes. It can be seen that the laterality is kept in the graph structure (also for subcortical regions) and the whole graph is highly symmetric. Node size linearly represents the graph theoretical measure of structural degree for each node. Most important hubs are subcortical regions. The shown features of symmetry, laterality, homology, and subcortical hubs indicate that the averaged SC still kept its physiological characteristics.
default value for inhibitory time constant τ i (Equation 14 and Figure 3). We assume the healthy brain without super-threshold Abeta burden operates in a region of the parameter space, which is close to a network criticality. A criticality describes an area in the parameter space, where subtle changes of one variable can have a critical impact on others (Strogatz, 2015) (in this case bifurcations, see Supplementary Figure 1. The thresholding "cut-off " value β off -differentiating normal form pathological Abeta burden-was chosen according to the literature, stating that only after a certain level of tracer uptake a region is considered pathological (β off = 1.4, see above). The maximum possible Abeta burden value β max was chosen to be the 95% percentile of the Abeta regional SUVR distribution across all participants. The midpoint of the sigmoid was chosen such that it was half the way between β off and β max . The steepness was chosen such that the function converges to a linear function between β off and β max .

Brain Network Model Construction and Simulation
For the reasons stated in the above introduction, for our simulation approach we selected the Jansen-Rit model (Zetterberg et al., 1978;Jansen and Rit, 1995;Wendling et al., 2000;David and Friston, 2003;David et al., 2006;Spiegler et al., 2010Spiegler et al., , 2011Sanz-Leon et al., 2015;Kunze et al., 2016). The differential equations are presented in Equations 3-5 (Jansen and Rit, 1995). The employed parameter values can be found in Table 2.
As a general approach, the impulse response function (IRF) of a neural mass allows to transform an incoming action potential into a PSP by using a linear time-invariant system. The IRF is the transfer function of the system, which is convoluted with the incoming input (action potentials) to calculate the output (PSPs). The general form of the IRF is the systems output to a (infinitesimal short and high) Dirac impulse and can be estimated experimentally by using short impulses or step functions (Lopes da Silva et al., 1974;Freeman, 1975).
The excitatory IRF h e (t) is therefore specified as h e (t) = tH e exp(−t/τ e )/τ e , where τ e is the excitatory time constant (the time until the PSP reaches its maximum), H e is a coefficient of the PSP amplitude and t is time.
Similarly, the inhibitory IRF h i (t, β) is specified as with the same variables as above. As we will describe below in detail, the inhibitory IR is a function of the spatially distributed Abeta burden β, which affects the time characteristics τ i (β) and therefore h i (t, β) is a function of time and space. These IRFs can be translated into second-order ordinary differential equations by interpreting them as Green's functions. See Spiegler et al. (2010) for a detailed explanation of the dimensional reduction used here.
The differential equations that describe the network of three neural masses are now presented in Equations 3-5. The variables used for the simulations are listed in Table 2: Excitatory projections υ 1 onto pyramidal cells at location a in discretized space (a = 1, 2, . . . , N: N = 379 regions): Inhibitory projections υ 2 onto pyramidal cells at location a in space: Projections υ 3 of pyramidal cells a in space: dυ 3,a (t)/dt = υ 6,a (t) dυ 6,a (t)/dt = S(υ 30,a )(t)/τ e −2υ 6,a (t)/τ e − υ 3,a (t)/τ 2 e , wherein c 31 , c 13 , c 23 are the local connectivity weights between the three neural masses. Equation (4) shows the spatial dependency of the activity of inhibitory interneurons projected onto the pyramidal cells by τ i (β a ). Taking into account the biologically plausible configuration of the Jansen-Rit model shown in Figure 2, we use several mathematical simplifications to reduce the model's dimensionality without loss of generality. Taking into account that Equation (1) is valid for all excitatory input at dendrites irrespective of the source allows for using one single IRF at the pyramidal cells and, thanks to linearity, translates the summation of excitatory postsynaptic potentials into a sum of incoming firing rate, that is, m 3T,a (t) + c 31 S(c 13 υ 3,a (t)) in Equation (3), describing the excitatory projections onto pyramidal cells υ 1 . This simplification is without restrictions, simply exploits the linearity of the operators and reduces the dimensionality by 2. Furthermore, to adjust notation, the postsynaptic potentials caused by the inhibitory neural mass at pyramidal cells are denoted as and its kernel is as The projecting variable of one brain region at location a to other regions in the network is the sum of excitatory and inhibitory postsynaptic potentials at the local neural mass of pyramidal cells transferred into a firing rate using a sigmoid. The general form of this transfer function is with, λ = υ, S υ, max = 2e 0 and S υ , min = 0 for the potential-tofiring-rate transfer. Incoming mean firing rates m 3T,a (t) at the pyramidal cells at location a from other brain regions b = 1, 2, . . . , N, where N is the number of 379 regions are given by where m 3T,0 is baseline input m 3T,0 = const. for ∀t and all locations ∀a. The global coupling factor G is a coefficient of the incoming activity and therefore scales the connections w a,b incoming at location a from all b provided by the SC. Because of this, G is the crucial factor that moderates the influence of the network to each neural mass and therefore mediates the difference between an uncoupled systems (G = 0) and a strong connected system. In all populations, the state variable [υ 1 , υ 2 , υ 3 ] a (t) are the mean membrane potentials and the derivatives thereof with respect to time t, namely [υ 4 , υ 5 , υ 6 ] a (t) represent the mean currents.
To model how the local Abeta load β a , measured by the Abeta PET SUVR is affecting the inhibitory time constant we introduce a transfer function (Figure 3). The primary assumption of this transfer function is a dependency of the E/I balance on the local Abeta concentration as described above. With higher concentration of Abeta, we assume dynamic changes in the inhibitory population of the neural mass that lead to local hyperexcitation. To model this inside the existing Jansen-Rit equations, potential candidate parameters would be H i and τ i as well as c 32 . The coefficient H i is not a suitable candidate because it has no direct physiological correlate. The coupling coefficient c 32 corresponds best to synaptic transmission from inhibitory to pyramidal cells and therefore can be mainly seen as a receptor surrogate. The time constant τ i acts as a filter for IPSPs and correlates best with the evidence of decreased IPSP firing rate (Busche et al., 2008;Grienberger et al., 2012;Limon et al., 2012;Verret et al., 2012;Hazra et al., 2013;Ripoli et al., 2014;Ren et al., 2018) and is moreover well explored for the Jansen-Rit model (Wendling et al., 2002;Spiegler et al., 2010). As described above, we expect this transfer function to behave in an asymptotic way for Abeta concentrations below and above a specific range. We determined the lower border at β a,off = 1.4 and the upper border at the 95th percentile in our data at β a,max = 2.65. By exploring the effects of τ i in a single region model, we determined the effective range 14.29 ms ≤ τ i < 50 ms. Based on this range we defined the following three-conditional linear function wherein τ i,min and τ i,max are the maximum and minimum values for τ i , m = (τ i,max -τ i,min )/(β a,max -β a,off ) = 28.6 and c = m · β a,off -τ i,min = 25.7. Since this function is not differentiable in β a,off and β a,max , we used the sigmoid function Equation (12) instead, which is continuous and differentiable. Moreover, a sigmoid can be interpreted as the cumulative (of a logistic distributed) activity acquired by the PET of a small brain volume (voxel) with a low spatial resolution of about 2.5 mm and above (Moses, 2011). Therefore, the Abeta transfer function is defined as wherein r βa is the slope of the sigmoid, β 0 is the midpoint of the sigmoid and the coefficients are chosen to fit the conditions explained before. In this function, τ i appears as its reciprocal value τ i −1 as it is implemented in the code of TVB. Because τ i is a time in ms, the inverse of τ i is a rate of potential change, and does not directly correspond to a firing rate. The Abeta load affects the inhibitory rate following a sigmoid curve. The rate ranges between S min and S max and the time constant ranges consequently between 1/S max and 1/S min . To simulate the model using TVB, physical space and time are discretized. The system of difference equations is then solved using deterministic Heun's method with a time step of 5 ms. We used a deterministic method to avoid stochastic influences since the simulation was performed in the absence of noise.
The system was integrated for 2 min and the last minute was analyzed in order to diminish transient components in the time series due to the initialization and settle the system into a steady state.
We explore a range of 0 ≤ G ≤ 600 which provides an overview about the possible population level behaviors at different states of network coupling. Because the coupling factor G has a crucial influence on the external input on the neuronal populations, this allows different regions to operate in different dynamical regimes, as it can be seen in the bifurcation diagrams of Supplementary Figure 1. Global coupling factor G that was sampled between G = 0 (i.e., isolated regions) and G = 600 with a step size of G = 3. The initial values were taken from 4,000 random time points for each state variable in each region. The length of 2 min for the simulations was chosen with the aim to diminish possible transient components due to the initialization of state variables at t = 0. For analysis we used only the second minute of the simulated signals. No time delays are implemented in the large-scale network interactions since they are not required for the emergence of the here evaluated features and setting them to zero increases reduces required computation resources.

Spectral Properties of the Simulated EEG
In TVB, we simulate EEG as a projection of the oscillating membrane potentials inside the brain via its electromagnetic fields to the skin surface of the head (Sanz-Leon et al., 2015) using the individual lead field matrices which take into account the different impedances of white matter, gray matter, external liquor space, pia and dura mater, the skull and the skin. Our leadfield matrices considered the impedances of three compartment borders: brain-skull, skull-scalp and scalp-air (Jirsa et al., 2002;Bojak et al., 2010;Litvak et al., 2011;Ritter et al., 2013). The postsynaptic membrane potential (PSP) considered for the projection is the one of the pyramidal cells, as they contribute the mayor part to potential changes in EEG (Kirschstein and Köhling, 2009). The PSP is calculated by summing the synaptic input from excitatory and inhibitory subpopulations to the pyramidal cells. The baseline PSP was derived as the mean PSP across time for every region. For the LFP or EEG peak frequency, we computed the power spectrum using the "periodogram" function of the Scipy python toolbox (Jones et al., 2001). From the spectrogram the "dominant rhythm" was identified as the frequency with the highest power.

Abeta-Inferred Dynamics Lead to Individual Spectral Patterns
We analyzed the dominant frequency in the simulated EEG and regional neural signal (referred to as local field potential (LFP) (Figures 5G-J).
We observed a physiologically looking irregular behavior with two frequency clusters in the alpha and in the theta spectrum ( Figure 5G). This behavior is expressed in the area of lower global coupling G for all 10 AD participants and in 3 out of 8 MCI and 4 out of 15 HC participants. The irregular time series and the broad continuous frequency spectra (Figure 5B) of network regime in 0 < G < 150 are indicative for deterministic chaos. Such chaotic network regimes in a BNM have already been reported using the same local dynamic model [ Figure 2 in Kunze et al. (2016)]. Beside this emerging chaotic behavior in our simulations other phenomena occurred in the parameter space exploration: a state of hypersynchronization between regions (Figures 5H,J) and a state of a "zero-line" with no oscillations that clearly does not reflect a physiological brain state (Figures 5I,J). Abscissa is frequency, ordinate is an estimate of the spectral power (dimensionless equivalent of amplitude per 1 Hz). Colors are representing the EEG frequency bands from delta to beta, indicated with Greek letters (note that this is regional neural activity, not EEG). Corresponding time series on the right: neural activity at single regions, each showing 5 s. Abscissa is time, ordinate is a dimension-less equivalent of the electric potential. (A) Shows an irregular, amplitude modulated alpha to beta rhythm, (B) an irregular theta with some delta and alpha inside. In (C) we can observe a monomorphic spike signal with a theta/delta frequency of 3 Hz and higher order harmonies. In order to locate the individual simulations in the spectrum of possible dynamics, meaning in the range of possible Abeta load, we examined extreme values of Abeta distribution. The virtual brains with a mean Abeta load of zero (Supplementary Figure 3A) and with the maximum Abeta load at all regions (Supplementary Figure 3B), we see as expected for the Abeta-free system a behavior similar to the low-Abetacontaining HC participants. This is not surprising, because when the HC subjects do not have a high Abeta signal, the dynamics will converge to those with zero Abeta, which is in fact then only determined by the underlying standard SC and therefore remains the same for all participants. However, the homogeneous application of maximum Abeta burden does not lead to an AD-like pattern but shows a zero-line at the whole spectrum.
To give a mathematical explanation of those phenomena, we related each participants Abeta-burden to the corresponding inhibitory time constant τ i and used former analyses of the uncoupled local Jansen-Rit model (Spiegler et al., 2010) to estimate the bifurcation diagrams for the coupled system in this study (Figure 6). Shown diagrams allow to predict and explain the occurrence of alpha and theta rhythms or zero-lines depending on the underlying Abeta burdens. The variation of τ i by local Abeta burden fundamentally influences the systems bifurcations by shifting the bifurcation point along the range of external input to the pyramidal cells. As a consequence, different values of Abeta lead to a variable occurrence of two limit cycles and a stable focus. Therefore, for a single region with constant external input on pyramidal cells, depending on Abeta the region might be in an alpha limit cycle, in a theta limit cycle, in a bistable condition where both cycles are possible or in a stable focus. Figure 7 displays how the mean dominant rhythms differ between the groups. In the range below G = 100 we find a slowing in the AD group. Since in the range of lower G all three groups exhibit realistic frequency spectra and no zero-lines we consider this range of G as "physiological." Significant differences appear between AD and non-AD for ranges of high and low G and also for high alpha and low theta rhythms (Figure 7). The heterogeneous distribution of Abeta (in contrast to an averaged homogeneous distribution) plays a crucial role in the development of this AD-specific slowing. This is indicated by simulations with the mean averaged Abeta of each participant mapped on all regions. The simulations revealed a regionally more homogenous behavior in all groups (Supplementary Figure 4). Moreover,  Spiegler et al. (2010). Contrary to the implementation we used for our present TVB modeling approach, here the bifurcation diagrams explore the behavior in an uncoupled system, and in accordance with (Spiegler et al., 2010) the IPSP amplitude coefficient Hi changes inversely to τ i to keep the product of synaptic gains and dendritic time constants constant. The default input m 3T,0 on pyramidal cells starts at a firing rate of 108.5/s. Because of the potential-to-firing-rate transfer function (Equation 11), global scaling factor G is affecting both the input currents and the firing rates. For higher values of G, the input on pyramidal cells is expected to increase. First Columns, panels (A-D): Bifurcation diagram with the default time constant of 14 ms. This appears in the simulation if the Abeta SUVR is below the clinical cut-off 1.4, because then the time constant is unaffected according the transfer function in Equation 14. An Abeta burden below 1.4 SUVR corresponds to THAL phase 0 (Murray et al., 2015), i.e., that we expect no Abeta pathology in such a brain region. In this situation, there is only one limit cycle existing, which produces a frequency in alpha range (A). After increasing the input on the pyramidal cells, the alpha cycle collapses and transforms to a stable focus, where no oscillations appear in the absence of noise. This is the "zero-line" in our results. ]. This homogeneity explains the low variance of rhythms shown in the lower G ranges of (B), because all regions are in the same limit cycle and in the absence of artificial noise there is no possibility for an amplitude modulating factor. The "alpha regime" appears for all regions with an Abeta PET SUVR below 1.95. This corresponds to neuropathologic THAL phases 0, 1, 2, and 3, i.e., the regions will have no severe Abeta pathology (Murray et al., 2015). Second column, panels (E-H): Bifurcation diagram with a time constant of 22 ms, which corresponds to an intermediate Abeta load and a bistable dynamical regime which occurs for time constants between 20 and 28 ms. This corresponds to Abeta PET SUVRs between 1.95 and 2.15 (THAL phase 4 or 5), i.e., moderate-to-severe Abeta pathology (Murray et al., 2015). (E) Starting at the blue line (initial condition in alpha cycle), with an increased input on the pyramidal cells (e.g., by the network) it gets possible to reach the second limit cycle, which produces a theta rhythm and coexists with the alpha cycle while the pyramidal input is in a specific range (120/s−170/s). When the input is increased too much (e.g., by many connections of the network or by increased coupling factor G), the theta cycle disappears and the system jumps back to the alpha cycle and later on to the stable focus, which shows no oscillations in the absence of noise. This can explain some of the spectral behaviors we observed typically in the AD group (F,H): It starts with chaotic rhythms in alpha (blue line) and theta (red line) and in the shown AD participant 1 then gets synchronized to either alpha or theta. With higher couplings, the frequency gets more probably synchronized to alpha (green line), because higher G indicates a higher pyramidal input and therefore a higher attraction of the alpha cycle. (G) Remarkably for the shown participant is the fact that the bistable behavior is caused by a very small amount of regions in bistable regime, which propagate the theta rhythm to most other regions in the area 200 < G < 300. Third column, panels (I-L): Bifurcation diagram with a time constant of 50 ms, which correlates to a 95th percentile Abeta load and above. Those high Abeta burdens lead to a theta dynamical regime, which occurs for time constants between 28 and 50 ms. In clinicopathology, this corresponds to Abeta PET SUVRs above 2.15-about 90% of those patients are can be classified as THAL phase 5, i.e., they have severe Abeta pathology (Murray et al., 2015). In comparison to panel (E), the alpha limit cycle disappeared in panel (I). Therefore, we expect only theta rhythms or an activity at the stable focus. The theta cycle now begins shortly above the initial condition of pyramidal input without the alpha cycle in between. For an initial input of 108.5/s the system is in a stable focus. This may explain why in the simulation with maximum Abeta load at all regions (so each with a time constant of 50 ms) we see a zero-line without alpha at lower G values (Supplementary Figure 1). (J,L) A state of theta-only rhythm appeared in few AD participants at higher Gs (blue line). In the spectral behavior of AD participant 7, we can moreover observe a strong bistable pattern with chaotic frequency distributions for G < 300. This is likely caused by the high amount of bistable regions (K), while the synchronization to theta in higher G is an effect of the high proportion of regions in theta regime.

Simulated EEG Slowing in AD Is Caused by Heterogeneous Abeta Distribution
with homogeneous distribution of Abeta the slowing in AD participants does not appear: we don't see a significant change in the theta band ( Figure 7B). This is a strong indicator for the importance of the individual Abeta distribution and a proof for the necessity of heterogeneous excitotoxic effects for the creation of neural slowing. However, the absence of slowing in the simulations with homogeneous Abeta distribution does not proof the importance of a specific spatial Abeta pattern. In contrast, it only shows that there must be few regions with very high Abeta pathology to slow down the system (see Figures 6F-G).
FIGURE 7 | AD-specific slowing in EEG and LFP and influence of the heterogeneous pattern of Abeta distribution to the spectral behavior. (A,B) The panels show the "spectrograms," more precise the amount of regions with a dominating frequency averaged for all G values and the subjects of each group. Below, black bars are indicating significant differences for all 90 examined frequencies by a Kruskal-Wallis test (compared were the means of the amount of regions in each group having this particular frequency). In (A), for the empirical Abeta distribution pattern, the red dotted line (AD) diverges from the non-AD participants with a strong presence of dominating theta (peak at 4 Hz) and the absence of zero-line rhythm (except of very few regions, see arrow). Significant differences only appear between AD and each HC and MCI, namely for high alpha/low beta and for theta/delta (black bars). At f = 1.2 Hz (red bar), the significance level is also achieved when using a strict Bonferroni correction (p < 0.05/90). In contrast, (B) shows the same plot if the spatial distribution was "blurred": There is no visual difference between the behavior of the three groups, and also no theta rhythm is existing in the simulations. All groups have a dominating zero-line behavior averaged across the full G range (see arrow). However, there are some frequencies that significantly differ between AD and each HC and MCI in alpha / beta range, which could be also visually related to small peaks at the plots beside. In theta and delta, where we would expect to see the slowing, there is no significant difference at all. Due to readability, for (A,B) the y-axis was limited to the amount of 100 regions. In (A), the zero-line peak of HC and MCI ends at 211, in panel (B) all zero-line peaks end at 323. The different spectra lead to different G-dependent mean frequencies for the groups, which significantly differ in areas of high and low G: (C,D)-comparison of EEG and LFP between groups. Mean dominant rhythms across all simulated EEG channels (C) and region-wise simulated LFPs (D) for all analyzed global coupling values. The frequencies of AD patients are significantly different in EEG as well as in the regional neuronal population signal. Filled shapes and thin lines represent the quantiles at 0.95 and 0.05 for each group. (C) For EEG one can see that the 95%-quantile of AD and HC as well as MCI is not overlapping in the physiological area of lower G, where AD tends to slower frequencies. In a Kruskal-Wallis test, the difference between the means of all channel frequencies per subject in the three groups is significant for AD and non-AD at 0 < G < 60 (each AD to HC and AD to MCI: p < 0.0001). They are also significantly different in the area of higher G, where AD is faster-at 450 < G < 470 (each AD to HC and AD to MCI: p < 0.0001). (D) For simulated regional neural signal the slowing effect is less prominent. The broader range of frequencies for AD is represented by the high and low limit of the 95%-quantile. This can be related to the two frequency clusters in AD at alpha and theta, which are not frequently apparent in non-AD (as in Figure 5). In a Kruskal-Wallis test, the difference between the means of all regional frequencies per subject in the three groups is only continuously significant for AD against HC at 400 < G < 450 (AD compared to HC p < 0.0001). For the other comparisons, only isolated G values deliver significant differences in the area of low G (HC and MCI) and intermediate G (AD and MCI). Because of the big amount of tests necessary to test all global coupling values, none of the tested G values achieved Bonferroni corrected significance. However, because we assume that neither the frequencies at (A,B) nor the G values at (C,D) are independent variables (which is also the reason for the "grouped" clusters of significance at alpha and theta and G = 50 and G = 450), a Bonferroni correction is not necessary.

Intra-individual Ratio of High vs. Low Abeta Burden Across All Regions Determines Simulated EEG Frequency Spectrum-Distinct Spatial Configurations of Abeta Do Not Matter for Slowing
We next examined how LFP/EEG slowing is related to the underlying Abeta burden (Figure 8). We revealed significant linear dependencies for all groups between Abeta burden and frequency. We found a strong inverse dependency for AD (R 2 = 0.625), i.e., an Abeta-dependent EEG slowing. In contrast, for non-AD participants the relation was revers, i.e., higher values of Abeta caused EEG acceleration.
To test if specific regions are more important for the observed phenomena, we had to overcome the bias that only specific regions were strongly affected by Abeta. I.e., for the empirical Abeta distribution we cannot say e.g., for a region with high Abeta if it shows EEG/LFP slowing only because of its high Abeta value or because of its specific spatial and graph theoretical position in the network. Therefore, we next performed simulation with 10 random spatial distributions of the individual Abeta PET SUVRs for the 10 AD participants. In these simulations, the neural slowing appeared similarly to the empirical spatial distributions of Abeta (Supplementary Figure 5), which indicates a minor role of the distinct spatial patterns of Abeta. Instead, the ratio of regions corresponding to the three different dynamical regimes (alpha, theta, and bistable) determined the simulated frequency spectrum (Supplementary Figure 6). For an optimal value of G with 100 < G < 150, the ratio of regions with an Abeta value in theta regime best corresponded to the ratio of regions with theta frequency in LFP. Moreover, the number of regions in different regimes enables to predict the individual spectral behavior across G. This analysis shows the crucial role of G for the simulation dynamics. There might exist different optima of Gdependent on what phenomenon in the simulation is of interest. But for a specific phenomenon, in this case the correspondence of underlying Abeta PET to frequencies, we can find a narrow optimum of G wherein a specific behavior occurs.
The results of random spatial distribution of Abeta PET SUVRs were also used for a parameter space exploration (Figure 9). The analysis reveals that (1) alpha rhythms are only apparent for low time constants with τ i < 30 ms, but for the full spectrum of G, more probable for lower G values; (2) relevant amounts of bistable rhythms are only apparent for 17 ms < τ i < 39 ms and G > 120; (3) theta rhythms are present across almost the full spectra of G and τ i , with an equal appearance across G, but with a local minimum at τ i ≈ 18 ms, where the system is dominated by alpha and bistable rhythms. This exploration demonstrates two major insights. First, it confirms the crucial role of τ i for the appearance of alpha or theta rhythms as we expect it out of the (non-coupled) bifurcation diagrams of Figure 6. Network effects are present (e.g., there are theta rhythms for low values of τ i ), but play a minor role here. Second, the value of G does not significantly affect the probability of theta rhythm, except of an alpha-theta shift for low τ i < 20 ms and higher G > 160. This is caused by the coexistence of stable focus in alpha regime and theta limit cycle in theta regime for high pyramidal input (Figures 6A,I).

Neural Slowing Propagates to Central Parts of the Network Independently of the Spatial Abeta Distribution
In the analysis of spatial distribution in relation to the organization of the underlying SC network (Figure 10), it can be seen that unless Abeta is distributed more peripherally, the FIGURE 8 | Abeta-dependent slowing of LFPs is specific for AD participants. Meanwhile there is a significant linear dependency between Abeta and LFP frequency for all groups, only for AD a higher burden of Abeta leads to a decrease of frequency. HC and MCI show inverse correlations. Plotted are density plots showing the dependency between the local Abeta loads and LFPs. (A) HC group, (B) MCI group, and (C) AD group. The matrices are containing the resulting regional peak frequencies for all examined coupling values G for all participants. Linear regressions (black lines) revealed highly significant regression coefficients (p < 0.0001). A strong linear dependency between mean Abeta and LFP, that explains the greater part of the variance, is only apparent in the AD group (C). 37.5% of the variance yet cannot be explained by this linear dependency. Moreover, only for AD the dependency leads to slower frequencies for higher Abeta SUVRs, meanwhile HC and MCI have slightly faster frequencies for higher Abeta SUVRs. Visually one can see at least four contributing patterns in the AD group (C): (1) the linear decrement of frequency for higher Abeta, shown by the regression line, (2) the two frequency clusters (orange spots) at alpha and theta, (3) some regions with the zero-line behavior, particular those with low Abeta (thin line at the left, with SUVR of about 1.5), and (4) a broad variability of frequencies for regions of the same Abeta SUVR (horizontal distribution). These phenomena cannot be explained completely by a linear dependency and moreover not by a linear system at all. The criticality that divides the dynamics into three different frequency modes (zero, alpha and theta) is a phenomenon of the Jansen-Rit model as a non-linear system ( Figure 6 and Supplementary Figure 6) and the broad frequency distribution is (probably) a network effect.
FIGURE 9 | Alpha and bistable rhythms only appear in a specific part of the parameter space between G and τ i . This parameter space exploration was done by coupled simulations and therefore includes network effects. Frequency (by color) is presented dependent on global coupling G (x) and inhibitory time constant τ i (y). Projections to G and τ i are shown beside the matrix plot, here the frequencies are classified into alpha rhythm (f > 8 Hz), theta rhythm (f < 5 Hz) and bistable rhythms (5 < f < 8). No relevant proportion of zero-lines appeared in the simulations. The difference to empirical EEG classes (with slightly lower borders for theta, meaning more exactly a theta/delta rhythm) are reasonable here because of the knowledge of only two different limit cycles in the examined configuration of the Jansen-Rit model ( Figure 6). This is also the reason for the classification of frequencies between 5 and 8 Hz as bistable. The exploration was non-systematically performed by using all regions of random distributed Abeta SUVR values of the 10 AD participants, with 10 iterations of randomization per participant. However, except single values of τ i , the full spectrum of τ i could be explored. Single empty columns are filled with neighbor columns for better readability. In principle wee the an "isle" of alpha for low coupling and low time constant, while the rest of the dynamics is dominated by theta and delta. A full frequency spectrum (also green and yellow colors) is only apparent near the borders of the alpha isle in higher coupling.
Abeta-dependent effect of neural slowing is focused to central parts of the network. Even a random distribution of Abeta SUVRs leads to this effect ( Figures 10E,F), indicating that this is a network effect. Probably this phenomenon is caused because the slowing effects are not only affecting the region itself, but also its local circuitry and neighbored regions. Hubs with a high degree and many close neighbors are therefore more probable of being affected by slow rhythms propagated by other regions. To relate this to empirical facts: We know from our data ( Figure 10A) that Abeta is not deposited in hubs, but more in peripheral regions of the networks. This shows, however, how the consecutive pathologic slowing effect is afterwards focused to central and important parts of the networks. A weak peripheral affection of the inhibitory system therefore disturbs the full system seriously.

Virtual Therapy With the NMDA Antagonist Memantine
The former analyses have shown that Abeta-mediated simulated hyperexcitation can lead to realistic changes of simulated brain imaging signals in AD such as EEG slowing (Figures 5, 6). We therefore wanted to know if an established way to protect the brain of the hyperexcitation, which is the NMDA antagonist memantine, can lead to functional reversibility.
The idea in our model is now that in theory memantine acts anti-excitotoxic via its NMDA antagonism and should therefore be able to weaken the hyperexcitation we introduced to the system by Abeta (Figure 11). As mentioned above, the local coupling parameter c 31 represents the main part of the glutamatergic transmission and can therefore also be seen as a surrogate of NMDAergic transmission ( Figure 11A). We homogeneously increased the default value of c 31 stepwise to observe the effects on the system. The analysis of the Jansen-Rit equilibria supports the concept of lower excitation introduced by lower c 31 (Figures 11B,C). The equilibrium manifold is the manifold of fixed points (the equilibrium) that is projected onto the PSP at the pyramidal cells as a function of two parameters, that is, the local excitatory-to-pyramidal coupling coefficient c 31 and the input on pyramidal cells m 3T,0 . The manifold is the FIGURE 10 | Theta rhythms affect central parts of the network independently of the spatial distribution of Abeta. (A) Abeta PET SUVR for AD participants: the distribution is diffuse along the cortex with no strong affection of subcortical hubs. This well corresponds to the neocortical stage C of Abeta distribution Braak, 1991, 1997;Braak et al., 2006). (B) There is no linear dependency between the Abeta SUVR and the structural degree, as the graph above already indicates. In contrast to that, (C) shows the distribution of theta rhythm, computed as the proportion of each regions simulations (201 for different values of G for 10 subjects) with dominant theta rhythm (here simplified as a frequency that is below 8 Hz and not zero, so more precise the theta-delta-band). The patterns are not consistent with those of panel (A). This indicates that not the distinct region affected by Abeta is crucial, but more its local circuitry. Moreover, one can observe that regions with a higher degree often have a high appearance of theta rhythm (D) and show a linear dependency with R 2 = 0.183, in contrast to the distribution of Abeta (B), which hasn't shown such a dependency. This phenomenon is stable also for the random spatial distribution of Abeta SUVRs (E): Here we see even a stronger dependency (R 2 = 0.29) between structural degree and theta rhythm (F). This is remarkable because (unless the spatial distribution is random) the "pathologic" theta is focused on the hubs. This indicates that there must be network effects which concentrate the appearing theta to those regions with higher degree. object onto which the system is moving or collapsing dependent on the parameters-in a way the equilibrium that underlies the dynamics of the system. The virtual memantine leads to a partial reversibility of the altered dominant frequencies in AD compared to HC/MCI ( Figure 11D). Virtual memantine increases the mean dominant EEG frequency. These simulated functional effects do not imply reversibility of neurodegeneration, but they illustrate how pharmacological intervention can theoretically counteract those processes. This observation provides first a potential mechanistic explanation of the pharmacodynamics of memantine. Second, it shows that TVB in general and the Abeta-hyperexcitation model of this study in particular are able to test the efficacy of treatment strategies such as drugs and have therefore the potential to be used for the discovery of new treatment options. Finally, it supports the concept of this study, where the impaired inhibitory function is modeled by an increased synaptic delay and it indirectly indicates that higher Abeta (by increasing τ i ) has led to a local hyperexcitation. It is to mention, that in an uncoupled network both the decrease of c 31 (memantine) and the increase of τ i (by Abeta) would have the same effect (Figures 11B,C, 13C,D). In a coupled simulation, the effects are in contrast antagonistic. One reason for this seems to be, that the effect of virtual memantine is focused to central parts of the network (Figures 11E,F)-the same parts, where slowing ( Figure 10) and Abeta-derived hyperexcitation (Figure 12) are occurring. The homogeneously applied memantine evolves its action, guided by the topology of the SC network, along the same path as the hyperexcitation is distributed. The effects of altered delay of GABA transmission can be reversed by adjusting NMDA transmission at another subset of the local population model. This illustrates that theoretically an alteration of the inhibitory transmission dynamics may lead to disinhibition causing hyperexcitation in downstream populations, which is reversible by reduction of excitatory input.

DISCUSSION
Local Abeta-mediated disinhibition and hyperexcitation are considered candidate mechanisms of AD pathogenesis. In TVB simulations, the molecular candidate mechanism has led to macro-scale slowing in EEG and neural signal with a particular shift form alpha to theta previously observed in AD patients (Loewenstein et al., 1989;Rice et al., 1990;Dierks et al., 2000;Babiloni et al., 2016;Malek et al., 2017;Nakamura FIGURE 11 | Modeling NMDA antagonism by virtual memantine. We modified the local dynamics for the AD group by homogeneously decreasing the coefficient c 31 , which represents the coupling from excitatory population to the pyramidal cell and therefore is a potential surrogate for NMDA receptor activity (A). The coefficient c 31 was decreased by 25% to model the effect of memantine and was applied homogenously to all regions for the 10 AD participants. (B) Equilibrium manifold of the Jansen-Rit model, wherein υ 30 is a function of the model input m 3T,0 and c 31 for the median Abeta load of AD participants β = 2.1447. Note the decrease of PSP at pyramidal cells υ 30 with increasing c 31 for a constant input level-this can also be seen in the top view of the same three-dimensional plot in panel (C), where our view is parallel to the z axis. The slope of the manifold with the input decreases with increasing c 31 . Blue areas indicate the lower branch of the equilibrium manifold and red areas the upper branch (white areas are unstable). This demonstrates, as we suggested, that when maintaining the same input level, a lower c 31 leads to a lower PSP at the pyramidal cells. (D) Mean EEG frequency for the three groups HC (blue), MCI (green), AD (red, with shadowed area for the range between 5th and 95th percentile) and AD with memantine (red dotted line). The virtual application of memantine shifts the AD group to the level of HC and MCI (arrow) and out of the variance of AD without memantine. (E) Change in the relative PSP (local PSP relative to mean of all regions) due to virtual memantine in AD patients. The decrement of activity seems to be focused to the hubs and to central parts of the network. Interestingly, this is the same topological area wherein relative Abeta-derived hyperactivity takes place ( Figure 12). (F) Shows a significant and strong linear dependency (R 2 = 0.783) between the effect of memantine and the relative hyperactivity of regions. The homogenously applied virtual memantine therefore acts selectively in those regions, where hyperexcitation is already there. These regions in central network parts are also those, where slowing effects appears ( Figure 10). The implementation of virtual memantine provides therefore a link between its supposed anti-excitotoxic effect (because it acts selectively in relatively hyperactive regions) and the reversed functional phenomenon of slowing (because slowing, hyperactivity and the memantine effect are all focused to the same parts of the network). et al., 2018). These observations cannot be directly inferred by the hyperexcitation implemented in our model. Because we standardized all other factors and used a common SC for all simulations, this approach enables to examine the effects of locally altered E/I balance on an individual whole-brain level but without any other confounding factors.
We showed that the slowing in simulated EEG and LFP is specific for the AD group (Figures 7, 8). This offers an explanation, how the shift from alpha to theta, that is observable in EEG of AD patients (Loewenstein et al., 1989;Rice et al., 1990;Dierks et al., 2000;Babiloni et al., 2016;Malek et al., 2017;Nakamura et al., 2018), could be explained on a synaptic level-namely by an impaired inhibition. This computational modeling result supports the findings of specific toxicity of Abeta to inhibitory neurons (Ripoli et al., 2014;Ulrich, 2015;Ren et al., 2018).
We demonstrate the computational principles underlying this Abeta dependent slowing of EEG/LFP (Figure 6,  Supplementary Figure 6). Dependent on the Abeta burden alpha, theta or bistable regime develop caused by an alteration of the inhibitory time constant that leads to changes of the systems bifurcation behavior (Figure 6, Supplementary Figures 1, 2, 6).
The simulated LFP/EEG slowing in AD patients crucially depends on the spatially heterogenous Abeta distribution as measured by PET-the slowing disappears when using a homogenously distributed mean Abeta burden instead for simulation (Figure 7). To exhibit the slowing effect few regions with high Abeta burden are required-while the FIGURE 12 | Local hyperexcitation is introduced by Abeta and spatially linked to LFP slowing. (A-D) Shows the relative firing rate of pyramidal cells, i.e., the difference between the local firing rate and the global firing rate (averaged across all regions for one simulation / one value of G). Because the maximum firing rate is 5/s, the relative firing rate ranges from −5/s to 5/s. The firing rate can be calculated by using the potential-to-firing-rate transfer function from Equation 11. Color indicates the natural logarithm of the number of regions at each point of the histogram. (A) In the absence of Abeta, where all regions have an inhibitory time constant τ i = 14 ms, the firing rate shows low heterogeneity. There are neither hypoactive, nor hyperactive regions-the whole systems activity is near the "baseline of the brain" (mean firing rate of all regions). (B) For HC participants, the Abeta burden of some regions already introduces heterogeneity. Although, most regions are still near the baseline, as also in (C) for the MCI participants. In panel (D), we can see that the AD participants have enough Abeta to introduce a strong local heterogeneity of firing rate.
There are now about 4 clusters, which are stable across G: one near the baseline [but in contrast to (A-C) not directly at the baseline, but below], another hypoactive and two hyperactive clusters of regions. The spatial distribution of the relative activity is demonstrated in (E): Here we can observe that the hyperexcitation of pyramidal cells is focused to the hubs and to central parts of the network, which is also shown in (G). This pattern is similar to what we could observe for the theta rhythm, which is also focused to central network parts ( Figure 10) and to the effect of memantine ( Figure 11). There is moreover a significant linear dependency between the relative PSP and the probability of theta rhythms (F) as well as a significant and strong linear dependency (R 2 = 0.594) between the relative PSP and the natural logarithm of the structural degree (G). specific location of these regions seems not to be relevant (Supplementary Figure 5). The crucial factor for AD-specific slowing behavior in our simulations is the presence of very few regions that are strongly affected by Abeta (Figure 6F and middle column of Supplementary Figure 6).
Independently of the location of high Abeta burdens in the simulated brain, slowing emerges at the core, i.e., hubs of the structural connectome (Figure 10). This indicates that the central parts of the system are impacted functionally by the Abeta burden. Moreover, it shows that while Abeta is often distributed in peripheral parts of the structural connectome, its functional consequences affect the important hubs. This could provide a possible explanation why a peripheral distribution of Abeta leads to severe disturbances of cognitive function.
Abeta leads by the disturbance of E/I balance to more local hyperexcitation (Figure 12). Because the range of activity is broader, we have more hypoactive and also more hyperactive regions (Figure 13). The local hyperexcitation is strongly corrleated with local LFP slowing (Figure 10) and also focused to the hubs of the network (Figure 12).
We also showed that the drug memantine that is known for improving brain function in severe AD can be modeled by a decreased transmission between the excitatory interneurons and the pyramidal cells and is able to achieve a "normalized" brain function in silico, too (Figure 11). Its effect is evolves selectively in hyperactive regions and in those parts of the network, where slow rhythms appear. This moreover demonstrates the potential of TVB to test and develop new treatment strategies.
One major limitation of this study is the lack of direct validation of the simulated electrophysiological phenomena. Neither EEG nor LFP data was available in the ADNI-3 cohort. Although EEG slowing in AD is an established concept (Loewenstein et al., 1989;Rice et al., 1990;Dierks et al., 2000;Babiloni et al., 2016;Malek et al., 2017;Nakamura et al., 2018), FIGURE 13 | Evaluation of the used cause-and-effect model. (A) Inhibitory convolution kernel h i (t) as a function of Abeta. The kernel is flattened with increasing Abeta and the area under the curve increases as follows: AUC = H i • τ i (β). Therefore, longer τ i by higher Abeta leads to a slowed down filter action (because of the delayed maximum of the IPSP). As a side effect, because of the higher AUC in the inhibitory transmission, the overall inhibition in the system increases-this can also be seen in (C) and (D): Equilibrium manifold (the manifold of fixed points onto which the system is collapsing) of the Jansen-Rit model in three-dimensional perspective (C) and in top view in parallel to the z axis (D), wherein υ 30 is a function of the model input m 3T,0 and β. The reaction of the system to input is changed by Abeta. Note the decrease of PSP at pyramidal cells υ 30 with increasing β for a constant input level-this can also be seen in (C). The slope of the manifold with the input decreases with increasing β. Blue areas indicate the lower branch of the equilibrium manifold and red areas the upper branch (white areas are unstable). This demonstrates, as we suggested because of (A), that when maintaining the same input level, a higher Abeta burden leads to a globally lower PSP baseline at the pyramidal cells. Because of the larger AUC at inhibitory transmission, we need more input to achieve the same global activity. As we would not expect the brain to increase or decrease its global activity level and are more interested in the local activity, this can be seen as a mathematical artifact as well as a limitation of the model. To evaluate local hyperexcitation, we therefore used the relative PSP to the baseline in Figure 12 instead of the absolute PSP. (B) The bifurcation diagram of local bifurcations shows the PSP of pyramidal cells as an explicit function of Abeta. We can see that the richest dynamic repertoire appears in the range around β = 2. Here appears multistability, which can be seen because of the folded Hopf curve. The local bifurcations support our hypothesis of diversified regional activity introduced by Abeta: above β = 2, two stable foci coexist, while one of them is positioned at a low PSP and the other one at a high PSP. This enables us to expect a diversified activity to nodes with lower and nodes with higher activity-as well, because the diagram shows only the local bifucrations, it does not proof it. (E) This leads together with our other results (in particular Figures 10-12) to the following mechanistic theory: longer inhibitory time constants, introduced by Abeta, lead to relatively more hypoactive and more hyperactive regions because of a disturbance of the local E/I balance in the vicinity of Abeta (Figures 12A-D). The hyperactivity is focused to the hubs and central network parts (Figures 12E-G). Because we know that this pattern does not evolve from the spatial Abeta distribution ( Figures 10A,B), we can say that the hyperactivity is centralized (to central network parts) and the hypoactivity is vice versa peripherialized (to the periphery of the graph). This phenomenon goes along with the distribution of slow rhythms in our results ( Figure 10). Therefore, our model linked local disturbed E/I balance with a slowed filter action of the inhibitory signal to hyperactivity along with slowing in central parts of the brain network. Parts of the figure are modified from Deco et al. (2017). future studies will have to validate the simulated data directly with individual EEG.
The second important limitation is the implementation of disturbed E/I balance by the inhibitory time constant. Although the longer time constants lead to slowed filter action ( Figure 13A) and local hyperexcitation at important network structures (Figure 12), the global activity is decreased for the same input level. To overcome this model limitation, it would be necessary to correct the input level, e.g., by increasing the default input m 3T,0 with higher mean Abeta burdens or by increasing a coefficient inside the IRF (Equation 2) to keep the AUC and therefore the transmitted energy at inhibitory transmission constant. This should be examined in future studies to evaluate the effect of such a correction.
However, this would only be necessary if the global activity level would be a target of interest for another research question. Because of the feedback loops in a coupled brain network, each introduction of over-or dis-inhibition will lead both to hypo-and hyper-active regions. An analysis of hyperactivity needs therefore always a control activitybecause hyperactivity can be meant spatially, temporally, or dependent on other factors. In our model, we could introduce spatially distributed hyperactivity (Figure 12) that was linked to local slowing (Figure 10), network topology (Figures 12E-G) and could be antagonized selectively by virtual memantine (Figure 11).
Another limitation is the small sample size of 33 participants. Future studies will have to consider much more participants, which will then help to formulate even more general conclusions. However, because of emergent effects in the brain simulation, differences between the groups were often very relevant and significant. An overview of all used statistical test in this study can be found in Supplementary Table 8.
However, we present a first proof of concept for linking molecular changes as detected by PET to large-scale brain modeling using the simulation framework TVB. This study therefore can work as a blueprint for future approaches in computational brain modeling bridging scales of neural function. For the research on AD pathogenesis, this study provides a possible mechanistic explanation that links Abeta-related synaptic disinhibition at the micro-scale to AD-specific EEG slowing. In general, our study can be seen as proof of concept that TVB enables research on disease mechanisms at a multi-scale level and has potential to lead to improved diagnostics and to the discovery of new treatments.

DATA AVAILABILITY
The raw data for this study is available in ADNI. The codes used in this study are available on request to the corresponding author.

ETHICS STATEMENT
This study has been approved from the Ethics Board of the Charité -Universitätsmedizin Berlin under the approval number EA2/100/19.

AUTHOR CONTRIBUTIONS
All authors have made substantial intellectual contributions to this work and approved it for publication. LS and PR had the idea to this study. LS, PT, ASp, and PR developed the concept and study design. LS wrote the manuscript, conducted the analysis and interpretation of results, and developed the figures. PT performed the MRI and PET image processing and supercomputer simulations. PT, ASp, M-AD-C, ASo, VJ, AM, and PR contributed to the interpretation of the results, figure development, and writing of the manuscript.

FUNDING
Data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on