# NETWORK SPREAD MODELS OF NEURODEGENERATIVE DISEASES

EDITED BY : Ashish Raj and Yasser Iturria Medina PUBLISHED IN : Frontiers in Neurology and Frontiers in Neuroscience

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-768-7 DOI 10.3389/978-2-88945-768-7

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# NETWORK SPREAD MODELS OF NEURODEGENERATIVE DISEASES

Topic Editors:

Ashish Raj, University of California, San Francisco and Weill Cornell Medicine, United States

Yasser Iturria Medina, McGill University, Canada

Image: GiroScience/Shutterstock.com

It is a real pleasure to introduce the excellent papers contributed by leading experts to our Research Topic proposing various network models of dementia spread. The focus is strongly on disease-specific mathematical modeling rather than general graph theory. The emerging field of network neuroscience visualizes the brain as a graph consisting of nodes representing regions and edges as connections between them. This complex network supports efficient communication along neural projections, but also, unfortunately, the transmission and progression of Alzheimer's and other neurodegenerative disorders. If we could know the brain's network organization, could we then predict how degenerative processes might develop on this network? As these papers demonstrate, the answer is, increasingly, yes.

Citation: Raj, A., Medina, Y. I., eds. (2019). Network Spread Models of Neurodegenerative Diseases. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-768-7

# Table of Contents


Yu-Chi Huang, Tun-Wei Hsu, Chau-Peng Leong, Han-Chin Hsieh and Wei-Che Lin

*18 Spatiotemporal Propagation of the Cortical Atrophy: Population and Individual Patterns*

Igor Koval, Jean-Baptiste Schiratti, Alexandre Routier, Michael Bacci, Olivier Colliot, Stéphanie Allassonnière and Stanley Durrleman for the Alzheimer's Disease Neuroimaging Initiative


Jordi Manuello, Andrea Nani, Enrico Premi, Barbara Borroni, Tommaso Costa, Karina Tatu, Donato Liloia, Sergio Duca and Franco Cauda

*70 Predictive Model of Spread of Progressive Supranuclear Palsy Using Directional Network Diffusion*

Sneha Pandya, Chris Mezias and Ashish Raj


Neil P. Oxtoby, Sara Garbarino, Nicholas C. Firth, Jason D. Warren, Jonathan M. Schott, Daniel C. Alexander and For the Alzheimer's Disease Neuroimaging Initiative

# Editorial: Network Spread Models of Neurodegenerative Diseases

Ashish Raj 1,2 \* and Yasser Iturria-Medina<sup>3</sup>

*<sup>1</sup> Department of Radiology, Weill Cornell Medicine, New York, NY, United States, <sup>2</sup> Department of Radiology, University of California at San Francisco, San Francisco, CA, United States, <sup>3</sup> Montreal Neurological Institute, McGill University, Montreal, QC, Canada*

Keywords: connectome, neurodegeneration, graph theory, Alzheimer's disease, protein aggregation, trans-neuronal transmission

**Editorial on the Research Topic**

#### **Network Spread Models of Neurodegenerative Diseases**

The emerging field of network neuroscience visualizes the brain as a graph consisting of nodes representing regions and edges as connections between them. This complex network supports efficient communication along neural projections, but also, unfortunately, the transmission and progression of neurodegenerative disorders like Alzheimer's disease (AD). If we could know the brain's network organization, could we then predict how degenerative processes might develop on this network? The answer is, increasingly, yes. This Research Topic collects a series of papers on various network models of dementia spread, focusing on disease-specific mathematical modeling rather than general graph theory.

#### Edited by:

*Wendy Noble, King's College London, United Kingdom*

#### Reviewed by:

*Delphine Boche, University of Southampton, United Kingdom*

\*Correspondence: *Ashish Raj asr2004@med.cornell.edu*

#### Specialty section:

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

Received: *24 October 2018* Accepted: *14 December 2018* Published: *08 January 2019*

#### Citation:

*Raj A and Iturria-Medina Y (2019) Editorial: Network Spread Models of Neurodegenerative Diseases. Front. Neurol. 9:1159. doi: 10.3389/fneur.2018.01159*

#### MODELING NEURODEGENERATIVE DYNAMICS ON BRAIN NETWORKS

Disturbances in global and local network organization are well documented in neurodegenerative diseases including AD (1), frontotemporal dementia (FTD) (2) and amyotrophic lateral sclerosis (ALS) (3). However, far more interesting than static networks is the potential to understand their dynamics—ongoing brain changes throughout disease (4). Broadly, this can happen in two ways: First, aberrant connectivity and network degeneration (5–7) via demyelination and axonal injury, secondary Wallerian degeneration, loss of signaling, axonal, and dendrite retraction. Second, disease factors can directly propagate along (possibly unchanging) neural connections, underpinned by "prion-like" protein aggregation followed by their trans-synaptic transmission (8–14). Thus, instead of being primarily impaired in degeneration, the network serves mainly as a conduit for disease transmission. Which type of dynamics predominates in neurodegeneration is a matter of lively debate.

In this Issue, Carbonell et al. provide a condensed historical review summarizing mathematical modeling of complex misfolded proteins mechanisms, both at the local/regional level and at the whole brain network level. They describe many recent network spread models, including models of cooperative spread (15), and communication cascades (16). In particular, an epidemic spreading model of network spread (17) was successfully validated on PET Amyloid-ß patterns in AD patients. An analytical Network Diffusion Model (NDM) mathematically derived the behavior of protein transmission as a graph heat equation under a connectivity-driven mechanism (18). This model is the basis of 2 papers in this Topic–Mezias and Raj and Pandya et al.

The accompanying paper by Mezias and Rajset out to apply the NDM, originally used on human data, to the problem of predicting the progression of Aβ pathology in transgenic mouse models. They show that while NDM is capable of predicting Aβ spread, it is not necessarily better than a model of spatial spread—one that does not involve network transmission. This is in contrast to the same authors' earlier publication showing a strong network mediation effect of tau spread (19), revealing a potentially important distinction between their respective mechanisms.

Pathological commonality and overlap is observed in various dementias: misfolded tau, A-beta and alpha-synuclein are present to varying degree in most degenerative diseases: AD, semantic dementia, FTD (20), ALS, dementia with Lewy bodies and posterior cortical atrophy (21). Pandya et al. explore this issue in a rare tauopathy called Progressive Supranuclear Palsy (PSP), which principally affects brainstem and striatal areas. Using the NDM Pandya et al. were able to recapitulate empirical PSP atrophy patterns. Although human imaging does not enable measurement of projection polarity, they presented a cross-species approach that transfers some directionality information from mouse to homologs regions in humans. Using this directional connectome, they found that anterograde and retrograde transmission give somewhat different spatiotemporal patterns of spread.

This raises important questions that can only be resolved by future advances in connectomics. Interestingly, the contributing paper by Neitzel et al. addresses this point directly in their perspective article. Using multimodal imaging data, including PET and functional MRI, they give a new perspective on how these techniques can be used to infer directionality of network connections in the human brain. As and when these techniques become increasingly refined and widely adopted, we predict that the issue of pathology spread on directional networks might assume critical ramifications.

The paper by Manuello et al. takes a meta-analysis approach for analyzing voxel-based morphometry data to understand the phenomenon of network spread. They report that in AD, gray matter alterations do not occur randomly across the brain but, on the contrary, follow identifiable patterns of distribution. This alteration pattern exhibits a network-like structure composed of co-altered areas that can be defined as co-atrophy network. The benefit of this type of data-driven approach is that it does not rely on specific hypotheses about network spread, unlike the papers by Mezias and Raj and Pandya et al.

A data-driven methodology was also adopted by Koval et al. which used the temporal alignment and combination of several short-term observation data to reconstruct the long-term atrophy

#### REFERENCES


history of AD. This model provided a description of both the spatiotemporal patterns of cortical atrophy at the group level and the variability of these patterns at the individual level in terms of propagation pathways, speed of propagation, and age at propagation onset. Oxtoby et al. explored how the pathology propagates through the connectivity network via a data-driven event-based model. By analyzing the changes in the elderly brain's anatomical connectivity over the course of AD, the authors clarified both the location and the sequence of changes to white matter connections. The results supported that degeneration of anatomical connectivity in the human brain may be an early and valuable biomarker when studying neurodegenerative diseases. Two other data-driven studies, Weber et al. and Huang et al., illustrate the methodological importance of considering focal axonal swelling and functional connectivity analysis for quantifiying memory deterioration rates in traumatic brain injury and evaluating clinical effects in infarction patients with dysphagia, respectively.

#### SUMMARY AND OUTLOOK

The present contributions encompass several cutting-edge techniques for investigating the phenomenon of networked spread in neurodegeneration, both model-based and data-driven. Together, they provide a framework for understanding the way the disease moves around within the brain network, that might adequately explain archetypal patterns of regional specificity in various dementias. Moving away from current statistical or descriptive graph theory, these papers trace underlying network dynamical processes. The emerging frontier of network spread modeling has the potential for wide applicability in diagnostics and therapeutic interventions.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### ACKNOWLEDGMENTS

AR acknowledges partial support by the following grants: R01 NS092802 and R01 EB022717 from the National Institutes of Health.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Raj and Iturria-Medina. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Clinical Effects and Differences in Neural Function Connectivity Revealed by MRI in Subacute Hemispheric and Brainstem Infarction Patients With Dysphagia After Swallowing Therapy

Yu-Chi Huang1†, Tun-Wei Hsu2†, Chau-Peng Leong<sup>1</sup> , Han-Chin Hsieh<sup>1</sup> and Wei-Che Lin<sup>3</sup> \*

<sup>1</sup> Department of Physical Medicine and Rehabilitation, Kaohsiung Chang Gung Memorial Hospital, Chang Gung University College of Medicine, Kaohsiung, Taiwan, <sup>2</sup> Department of Diagnostic Radiology, Taipei Veterans General Hospital, Taipei, Taiwan, <sup>3</sup> Department of Diagnostic Radiology, Kaohsiung Chang Gung Memorial Hospital, Chang Gung University College of Medicine, Kaohsiung, Taiwan

#### Edited by:

Yasser Iturria Medina, Montreal Neurological Institute and Hospital, McGill University, Canada

#### Reviewed by:

Karmele López-de-Ipiña, University of the Basque Country (UPV/EHU), Spain Kazutaka Takahashi, University of Chicago, United States

> \*Correspondence: Wei-Che Lin u64lin@yahoo.com.tw

†These authors have contributed equally to this work and first author.

#### Specialty section:

This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neuroscience

Received: 20 July 2017 Accepted: 29 June 2018 Published: 20 July 2018

#### Citation:

Huang Y-C, Hsu T-W, Leong C-P, Hsieh H-C and Lin W-C (2018) Clinical Effects and Differences in Neural Function Connectivity Revealed by MRI in Subacute Hemispheric and Brainstem Infarction Patients With Dysphagia After Swallowing Therapy. Front. Neurosci. 12:488. doi: 10.3389/fnins.2018.00488 Background: Early detection and intervention for post-stroke dysphagia could reduce the incidence of pulmonary complications and mortality. The aims of this study were to investigate the benefits of swallowing therapy in swallowing function and brain neuro-plasticity and to explore the relationship between swallowing function recovery and neuroplasticity after swallowing therapy in cerebral and brainstem stroke patients with dysphagia.

Methods: We collected 17 subacute stroke patients with dysphagia (11 cerebral stroke patients with a median age of 76 years and 6 brainstem stroke patients with a median age of 70 years). Each patient received swallowing therapies during hospitalization. For each patient, functional oral intake scale (FOIS), functional dysphagia scale (FDS) and 8-point penetration-aspiration scale (PAS) in videofluoroscopy swallowing study (VFSS), and brain functional magnetic resonance imaging (fMRI) were evaluated before and after treatment.

Results: FOIS (p = 0.003 in hemispheric group and p = 0.039 in brainstem group) and FDS (p = 0.006 in hemispheric group and p = 0.028 in brainstem group) were both significantly improved after treatment in hemispheric and brainstem stroke patients. In hemispheric stroke patients, changes in FOIS were related to changes of functional brain connectivity in the ventral default mode network (vDMN) of the precuneus in brain functional MRI (fMRI). In brainstem stroke patients, changes in FOIS were related to changes of functional brain connectivity in the left sensorimotor network (LSMN) of the left postcentral region characterized by brain fMRI.

Conclusion: Both hemispheric and brainstem stroke patients with different swallowing difficulties showed improvements after swallowing training. For these two dysphagic stroke groups with corresponding etiologies, swallowing therapy could contribute to different functional neuroplasticity.

Keywords: dysphagia, stroke, swallowing therapy, videofluoroscopy, magnetic resonance imaging

# INTRODUCTION

Dysphagia is a common disorder after stroke (Paciaroni et al., 2004). Approximately 25∼45% of patients show difficulties in swallowing after acute stroke while dysphagia is associated with a high risk of aspiration pneumonia, malnutrition, and mortality after acute stroke (Barer, 1989; Johnson et al., 1993; Martino et al., 2005; Walter et al., 2007; Falsetti et al., 2009). Previous studies have reported that 15% of stroke patients depend on tube feeding in order to maintain their nutritional requirements 6 months after stroke (Croghan et al., 1994; Dennis et al., 2005; Smithard et al., 2007). Consequently, the early detection and management of post-stroke dysphagia could reduce the incidence of subsequent co-morbidities and help to reduce the length of duration of hospitalization (Dziewas et al., 2004; Brotherton and Judd, 2007; Maeshima et al., 2014). Therefore, early swallowing therapy for dysphagia is very important in preventing further pulmonary complications and mortality after stroke (Doggett et al., 2001; Foley et al., 2008).

Brainstem or hemispheric lesions after stroke can lead to motor or sensory deficits in swallowing function. Brainstem stroke may also impair the sensory input or motor function of the mouth, tongue, cheek, pharynx, larynx, vocal fold, or cricopharyngeal muscles during the pharyngeal phase of swallowing (Horner et al., 1991; Moon et al., 2012). Hemispheric stroke may also impede the motor control and coordination of bolus preparation, mastication, or pharyngeal peristalsis during the oral and pharyngeal phases during swallowing. Therefore, different types of stroke may cause swallowing impairments in different ways in stroke patients with dysphagia (Robbins et al., 1993; Moon et al., 2012; Kim et al., 2014).

Several swallowing therapies have been utilized in the management of post-stroke dysphagia, including thermal stimulation, postural compensation, food consistency modifications, oropharyngeal exercises, swallowing maneuver and electrical stimulation (Carnaby et al., 2006). Recently, some investigators have used neuromuscular electrical stimulation for dysphagia patients after stroke and found this mode of treatment to be beneficial in the treatment of dysphagia (Lee et al., 2014; Poorjavad et al., 2014; Byeon and Koh, 2016). In earlier studies, some researchers reported that these therapeutic techniques might provide positive effects on neuroplasticity for stroke patients with dysphagia. Voluntary swallowing is regulated by complex sensorimotor cortical and subcortical network while the swallowing reflex is controlled by swallowing centers in the brainstem (Soros et al., 2009). However, the actual mechanism underlying swallowing function in the brain has yet to be elucidated. In stroke patients, damage to both of the associated hemispheres or brainstem result in swallowing disorders (Hamdy et al., 1996, 1997, 1998a,b; Li et al., 2009, 2014).

Previous studies have used clinical observations, EEG, magnetoencephalography (MEG) or transcranial magnetic stimulation (TMS) to explore the cerebral areas associated with swallowing dysfunction (Robbins et al., 1993; Hamdy et al., 1997; Luan et al., 2013; Jestrovic et al., 2016). However, over recent years, functional magnetic resonance imaging (fMRI) has become a more popular tool with which to evaluate cerebral cortical function during volitional and reflexive swallowing in humans (Hamdy et al., 1999; Mosier K. et al., 1999; Mosier K. M. et al., 1999; Kern et al., 2001; Martin et al., 2001; Malandraki et al., 2011; Dehaghani et al., 2016). Li et al. used fMRI to investigate cerebral cortical activation during swallowing with tasks in acute dysphagic stroke patients involving of unilateral hemisphere (Li et al., 2009). On the basis of these results, they suggested that fMRI is a useful method to investigate the spatial localization of changes in the neuroactivity of the bilateral hemispheres during swallowing tasks which may be related to the functional recovery of dysphagia (Li et al., 2009). The strongest activations were noted in the sensorimotor cortices, insula and cingulated gyrus of the intact hemisphere. Furthermore, they reported decreased connectivity in bilateral swallowing associated brain network while applying rs-fMRI study to explore the alternations of functional and structural connectivity in stroke patients with dysphagia (Li et al., 2014). Another fMRI study investigated recovered swallowing in dysphagic stroke patients and also showed overall decreased fMRI-activation in the associated swallowing network, but increased activities in the contralesional primary somatosensory cortex (S1) (Mihai et al., 2016).

Diffusion Tensor Imaging (DTI) is commonly used to study normal white matter anatomy and structural connectivity. This technique might have the potential to also assess plastic changes in the white matter in response to intense therapy in stroke patients recovering from a motor deficit. However, some of the DTI-derived absolute measures have shown pre-vs. post-therapy changes that were in the order of variability seen in normal subjects. Furthermore, significant remodeling of the ipsilesional and contralesional corticospinal tract was observed in chronic stroke patients undergoing an intense therapy program (Betzler et al., 2009). Consequently, in order to assess the sensitivity of this technique to detect possible structural changes as a function of therapy, it is necessary to carry out an assessment of the observed variance in DTI-derived measurements (Betzler et al., 2009).

To the best of our knowledge, there are limited published literature pertaining to the relationship between the changes in swallowing function and functional connectivity maps of fMRI, and in hemispheric and brainstem stroke patients following swallowing therapies. The aims of the present study were to investigate the benefits of swallowing therapy for swallowing function and brain neuro-plasticity and to explore the relationship between swallowing function and fMRI findings following swallowing therapy in hemispheric and brainstem stroke patients with dysphagia.

### MATERIALS AND METHODS

#### Participants

Between January 2011 and July 2013, we collected 31 acute stroke patients with dysphagia who met the inclusion criteria from the inpatient rehabilitation unit at one medical center. (clinical trial No.: NCT03048916) Stroke was diagnosed by the attending neurologist according to the patient's neurological insults and findings from brain computed tomography (CT) or magnetic resonance imaging (MRI) findings. The inclusion criteria for the participants were recent hemispheric or brainstem stroke (with a duration of less than 3 months since stroke), dysphagia reported by a physician who had assessed choking or coughing during swallowing, or patients were failed to complete the 100 ml water test during bedside swallowing assessment (Chen et al., 2016). We used functional oral intake scale (FOIS) to present patient's swallowing performance, which may reflex the severity of swallowing; patients were included in the study if their FOIS was equal to or less than 4. The exclusion criteria in the study were as follows: impaired communication ability due to aphasia; cognition impairment; a history of other neurological deficits leading to dysphagia; use of an electrically sensitive biomedical device (such as a cardiac pacemaker or a metal clip in the brain). Only seventeen stroke patients with dysphagia underwent detailed clinical assessments, videofluoroscopic swallowing studies (VFSS), and fMRI examinations before and after swallowing therapy. The study protocol was reviewed and approved by the Institutional Review Broad in our hospital. Informed written consent was obtained from each participant.

## Clinical Assessments and Swallowing Therapy

For each patient, the following clinical characteristics were recorded during admission into the rehabilitation unit: age; gender; swallowing therapy type; duration since stroke onset; the National Institute of Health Stroke Scale (NIHHS); affected hemisphere and lesion location. Finally, we recruited 11 patients with hemispheric stroke into the hemispheric group and 6 patients with brainstem stroke into the brainstem group; all of these patients underwent the entire evaluation and treatment process.

During hospitalization, these subacute stroke patients with dysphagia received traditional swallowing therapy or a combined therapy including both traditional swallowing therapies and neuromuscular electrical stimulation (NMES) therapy. Every patient was treated 3 times per week (60 min per session). Total ten treatment sessions with traditional or combined therapy were performed for each patient. The traditional swallowing therapy was performed by one experienced speechlanguage therapist, including oropharyngeal exercises, thermal stimulation, food consistency modifications, compensatory techniques, and swallowing maneuver based on VFSS findings and clinical presentations. NMES therapy was administered by one licensed physiatrist who used the VitalStim device (Chattanooga Group, Hixson, TN, USA) with a dual channel and two bipolar electrodes for each channel (700 µs pulse width, 80 Hz, and 0-25 mA wave-amplitude). The two sets of electrodes were placed on the patient's anterior neck area (Freed et al., 2001). Wave amplitude setting was dependent upon the patient's level of tolerance. The physiatrist gradually increased the amplitude until the patient felt a tingling sensation on the anterior neck and a muscle contraction. The current intensity was determined and maintained during one NMES treatment session based on each patient's tolerance level. The patients all underwent traditional swallowing training and NMES by the same physiatrist while receiving the combined therapies.

# Clinical Swallowing Function and Videofluoroscopy Evaluations

The functional oral intake scale (FOIS) for dysphagia was evaluated before and after intervention by a speech-language therapist who was blinded to all study procedures. The FOIS (Crary et al., 2005) is widely used for clinically assessing oral intake in stroke patients. Seven swallowing functional levels were defined according to their oral intake conditions ranging from nothing by mouth (level 1) to total oral diet with no restriction (level 7). The VFSS was performed by another speech-language therapist who was also blinded to the study interventions. The 8-point penetration-aspiration scale (PAS) (Rosenbek et al., 1996) and functional dysphagia scale (FDS) (Han et al., 2001) were assessed according to VFSS findings. The 8-point PAS indicated the severity of aspiration while swallowing. Aspiration was defined as any materials entering the larynx below the vocal folds. The 8-point PAS scale ranged from a score of 1 (normal swallowing without material entering the airway) to a score of 8 (severe airway compromise with material entering the airway and passing below the vocal folds). The FDS includes 11 items relating to oral and pharyngeal function and was investigated during the VFSS according to lip closure, bolus formation, residual matter in the oral cavity, oral transit time, triggering of pharyngeal swallowing, laryngeal elevation and epiglottis closure, nasal penetration, triggering, residue in valleculae and pyriform sinus, pharyngeal coating, and pharyngeal transit time. For each patient, the achievable score is 100 points depended upon the results of the VFSS. A higher FDS score indicated a worse swallowing condition in the VFSS. Another experienced and blinded speech-language therapist interpreted and scored the 8-point PAS and FDS before and after interventions.

## Acquisition of MRI Data

Functional imaging data were acquired using a 3.0 T GE Signa MRI scanner (Milwaukee, WI, USA). Resting-state images from 300 contiguous echo planar imaging whole brain functional scans were acquired (TR: 2 s; TE: 30 ms; FOV: 240 mm; flip angle: 80◦ ; matrix size: 64 × 64; thickness: 4 mm). During the resting experiment, the scanner room was darkened and the participants were instructed to relax, with their eyes closed, without falling asleep. A three-dimensional (3D) high-resolution T1-weighted anatomical image was also acquired using an inversion recovery fast spoiled gradient-recalled echo pulse sequence (TR: 9.5 ms; TE: 3.9 ms; TI: 450 ms; flip angle: 20◦ ; field of view: 256 mm; matrix size: 512 × 512).

# Resting-State Functional MRI (rs-fMRI) Preprocessing

Rs-fMRI data were preprocessed using Statistical Parametric Mapping (SPM8, Wellcome Department of Cognitive Neurology, London, UK; http://www.fil.ion.ucl.ac.uk/spm/) and Data Processing Assistant for rs-fMRI (DPARSF) and in-house software for Matlab (The MathWorks Inc. Natick MA, USA).

In the first step, the first 5 volumes were discarded to reach a steady-state magnetization and allow the participants to adapt to the scanning noise. We excluded any data involving a head motion of more than 2.0 mm maximum displacement in any of the x, y, or z directions, or 2.0◦ of any angular motion throughout the course of the scan. Data were also visually inspected for movement-related artifacts. The standard Montreal

#### TABLE 1 | Patients characteristics.

Neurological Institute template provided by SPM was further used for normalization with re-sampling to 3 mm cubic voxels and a Gaussian kernel of 6 mm (full width at half maximum) for spatial smoothing. The waveform of each voxel was finally used to remove linear trends of the time courses and for temporal


Patients characteristics including age at fMRI study date, gender, swallowing therapy, duration after stroke onset, NIHSS, affected hemisphere, lesion location. Gender: M, male; F, female; Affective hemisphere (Aff. Hem.): L, left; R, right; NIHSS, National Institute of Health Stroke Scale.

TABLE 2 | A comparison of findings in the VFSS and clinical oral intake conditions within and between hemispheric and brainstem groups.


Within-group and between-group comparisons of FDS, 8-point PAS, and FOIS were performed by the Wilcoxon signed-rank test and Mann-Whitney test. VFSS, videofluoroscopic swallow study; FDS, functional dysphagia scale; PAS, penetration-aspiration scale; FOIS, functional oral intake scale; IQR, interquartile range. p: comparison between before and after treatment in each group; p1, comparison between two groups before treatment; p2, comparison between two groups after treatment. \*p < 0.05.


TABLE 3 | A comparison of the changes in clinical outcome measures between the hemispheric and brainstem groups.

Between-group comparisons of the changes in FDS, 8-point PAS, and FOIS were performed by the Mann-Whitney test. FDS, functional dysphagia scale; PAS, penetrationaspiration scale; FOIS, functional oral intake scale; IQR, interquartile range. \*p < 0.05.

band-pass filtering (0.01–0.08 Hz) to reduce high-frequency physiological noise.

#### Independent Component Analysis of rs-fMRI Networks

To investigate rs-fMRI networks, we performed independent component analyses (ICA) using group ICA for fMRI toolbox (GIFT; http://mialab.mrn.org/software/gift). This toolbox supports a group ICA approach, which first concatenates the individual data across time, followed by computation of the subject-specific components and time courses. Data dimensionality (the number of components) was estimated using the minimum description length criteria tool in GIFT (Calhoun et al., 2001), which suggested that 30 was the optimal number of independent components (ICs). The dimensions of the functional data were then reduced using principal component analysis and the ICs were then estimated using the Infomax algorithm (Bell and Sejnowski, 1995). Stable estimation was achieved by rerunning the ICA analyses 20 times using the ICASSO (Himberg et al., 2004; Chenji et al., 2016) toolbox implemented in GIFT. Twenty-four ICs with an ICASSO stability index of <0.9 were then estimated in our patient sample and the spatial maps of the components were converted into z-value maps (Stevens et al., 2009). Further analysis of ICs were selected based on the largest spatial correlation (Calhoun et al., 2001) with specific rs-fMRI network templates reported in previous studies (Shirer et al., 2012), which correspond to known anatomical and functional segmentation. Group statistical maps of subject IC patterns representing rs-fMRI networks were entered into oneand two-sample random effects analyses in SPM8.

#### Statistical Analysis

SPSS package (version 17.0, Chicago, IL, USA) was used to run statistical analysis for group differences in demography, clinical characteristics, and VFSS findings. The Mann-Whitney test was used to analyse data relating to age, the duration since stroke, and NIHSS, and while Fisher's exact test was used to analyse gender, swallowing therapy types, and the affected hemisphere between the brain and brainstem groups. Within-group and between-group comparisons of FDS, 8-point PAS, and FOIS were performed with the Wilcoxon signed-rank test and Mann-Whitney test. The significance threshold of group difference was p < 0.05, FDR-corrected for multiple comparisons. Functional connectivity (FC) between-group two sample t tests were masked with a within-group mask threshold at p < 0.05, and corrected for multiple comparisons using a combination of an uncorrected height threshold of p < 0.05 with a minimum cluster size. Cluster size was determined over 1000 Monte Carlo simulations using the AlphaSim program distributed with the REST software tool (http://resting-fmri.sourceforge. net/). Anatomical labeling was defined by the Anatomical Automatic Labeling atlas (AAL) (Tzourio-Mazoyer et al., 2002). In voxel based statistical analysis, some significant clusters between groups were found. To evaluate their interactions with the clinical symptoms, the connectivity value was extracted from those significant clusters for further correlation analysis. The correlation between changes in clinical parameters and the alteration of functional connectivity was investigated using Spearman's rank correlation.

## RESULTS

#### Demographic and Clinical Characteristics

We initially collected 31 acute stroke patients with dysphagia but excluded 14 patients who showed poor cooperation during the VFSS/ fMRI study or swallowing training. Twenty-three out of 31 stroke patients received full-course swallowing training during hospital stay. Five of those 23 participants were not suitable to receive fMRI assessments and one of 23 patients could not completely take all 3 kinds of food in VFSS. In this study, we did not find any exacerbations of the medical conditions including aspiration or pneumonia during swallowing training and no harm related to all procedures for each participant. Consequently, a total of 17 patients completed all of the procedures and interventions required by this study.

Of the 17 patients included in the study, there were 11 patients with hemispheric stroke assigned to the hemispheric group and 6 patients with brainstem stroke assigned to the brainstem group. The demographic and clinical characteristics of these participants are shown in **Table 1**. In the hemispheric group (2 women and 9 men; median age: 76 years), 7 patients received combined therapy and 4 patients received traditional swallowing therapy. In the brainstem group (1 woman and 5 men; median age: 70 years), 4 patients received combined therapy and 2 patients received traditional swallowing therapy. There were no significant differences in terms of age, swallowing therapy type, duration since stroke onset, NIHSS, and affected hemisphere. Within-group comparison (**Table 2**) showed significant differences in terms of thin-liquid FDS (p = 0.012), total FDS scores (p = 0.006), FOIS (p = 0.003) in the hemispheric group after treatment, and significant differences in soft diet FDS (p = 0.043), thin-liquid FDS (p = 0.046), total FDS score (p = 0.028), PAS (p = 0.041), and FOIS (p = 0.039) in

default mode network; HVN, high visual network; LGN, language network; AN, auditory network; LSMN, left sensorimotor network; PCCN, precuneus network; CSMN, cerebellar sensorimotor network; aSN, anterior salience network; RECN, right executive control network.

the brainstem group after treatment. Between-group comparison showed significant differences in thick-liquid FDS (p = 0.049) and FOIS (p = 0.039) after treatment. There was also a significant difference between the 2 groups of patients after intervention in terms of their relative changes on the 8-point PAS (p = 0.028) (**Table 3**).

# Components of the Resting-State Functional Networks

Rs-fMRI networks were identified by using 14 templates in all ICA components for each of the 17 subjects. Correlation coefficients between the spatial templates and ICs of ICA analysis were as follows: visuospatial network, VSN (IC01), 0.292; TABLE 4 | Functional connectivity differences in rs-fMRI networks: a group comparison.


AAL, Automated Anatomical Labeling; MNI, Montreal Neurological Institute stereotaxic coordinate; IC, independent components; vDMN, ventral default mode network; SMN, sensorimotor network; VSN, visuospatial network; dDMN, dorsal default mode network; HVN, high visual network; AN, auditory network; PCCN, precuneus network; PVN, primary visual network.

posterior salience network, pSN (IC03), 0.474; left executive control network, LECN (IC04), 0.469; primary visual network, PVN (IC07), 0.398; right sensorimotor network, RSMN (IC09), 0.284; dorsal default mode network, dDMN (IC10), 0.471; ventral default mode network, vDMN (IC15), 0.263; high visual network, HVN (IC16), 0.330; language network, LGN (IC17), 0.330; auditory network, AN (IC21), 0.414; left sensorimotor network, LSMN (IC22), 0.264; precuneus network, PCCN (IC23), 0.368; cerebellar sensorimotor network, CSMN (IC24), 0.272; anterior salience network, aSN (IC26), 0.313; and right executive control network, RECN (IC28), 0.492. The results of the one-sample ttest (P < 0.05, FDR corrected) for patients with hemispheric and brainstem stroke were determined using SPM toolbox and are shown in **Figure 1**.

#### Functional Connectivity in rs-fMRI Networks: Group Comparison

We compared FC between patients with hemispheric and brainstem stroke in a voxel-wise manner. Significance differences arising from two-sample t-tests (between the two groups) and paired t-tests (between the two states) of spatial map regions are presented in **Table 4**.

The hemispheric group with post-intervention exhibited significantly increased intra-network functional connectivity in the left superior parietal lobule of the VSN, the right middle frontal cortex of dDMN, and the right cerebellum of SMN.

The brainstem group with post-intervention exhibited significantly increased intra-network functional connectivity in the left cerebellum and right calcarine of the PVN. In contrast, the brainstem group with post-intervention exhibited significantly reduced intra-network functional connectivity in the right postcentral cortex of the SMN, left superior parietal lobule of the vDMN, left precuneus of the vDMN, right fusiform of the HVN, left middle temporal lobe of the AN, left postcentral cortex of the SMN and left inferior parietal of the PCCN.

#### Relation Between rs-fMRI Networks and Clinical Parameters

Brain rs-fMRI showed that in hemispheric stroke patients, changes in the FOIS were related to changes in the functional connectivity of the ventral default mode network of the left precuneus (**Table 5** and **Figure 2A**). However, in brainstem stroke patients, changes in the FOIS were related changes in the functional connectivity of the left sensorimotor network of the left postcentral region. (**Table 5** and **Figure 2B**).

#### DISCUSSION

To the best of our knowledge, this represents the first research study to explore the relationship between brain connectivity network changes in rs-fMRI and clinical swallowing functional recovery following swallowing therapy in subacute hemispheric and brainstem stroke patients with oropharyngeal dysphagia.


TABLE 5 | Correlation between the changes of clinical parameters and altered functional connectivity maps after treatment in hemispheric and brainstem groups.

IC, independent components; vDMN, ventral default mode network; RSMN, right sensorimotor network; CSMN, cerebellar sensorimotor network; VSN, visuospatial network; dDMN, dorsal default mode network; HVN, high visual network; AN, auditory network; PCCN, precuneus network; PVN, primary visual network; LSMN, left sensorimotor network; PAS, penetration-aspiration scale; FOIS, functional oral intake scale. \*p < 0.05.

Our findings revealed that swallowing therapy led to a significant improvement in total FDS, FOIS scores in patients with hemispheric and brainstem stroke, and in 8-point PAS score in patients with brainstem stroke. Our research also showed that in patients with hemispheric stroke, the reduced brain connectivity of the ventral default mode network of the left precuneus was related to improvements in FOIS after swallowing therapy. In patients with brainstem stroke, the reduced left sensorimotor connectivity of the left postcentral network was associated with improvements in FOIS after swallowing therapy.

Both hemispheric and brainstem stroke patients with oropharyngeal dysphagia received swallowing therapy and achieved clinical benefits in terms of their swallowing function. In addition, brainstem stroke impaired the movements of the larynx, pharynx, vocal fold, or cricopharyngeal muscles during the pharyngeal phase of swallowing (Martin and Sessle, 1993; Carnaby et al., 2006; Chenji et al., 2016). In the present study, we observed a worse 8-point PAS score (based on VFSS findings) before intervention in patients with brainstem stroke. Therefore, this particular cohort of patients showed significant improvement after training; this was due to the patients developing better muscle coordination during the pharyngeal phase and reducing the severity of their aspiration. These mechanisms also reduced further pulmonary complications.

In earlier studies, associated cortical and sub-cortical regions involving in swallowing had been demonstrated (Teismann et al., 2009, 2011; Luan et al., 2013; Kober et al., 2015) in stroke patients with dysphagia. Li et al. (2009) reported that disrupted functional brain networks related to swallowing motor control in acute stroke patients with dysphagia according to the findings of rs-fMRI study. Some researchers investigated brain imaging using fMRI for acute and chronic stroke patients with dysphagia and they proposed that recovered swallowing were associated with increase of cerebral activation in the contralateral cortex of the intact hemisphere and ipsilateral anterior cerebellum (Li et al., 2009; Mihai et al., 2016). Another research study (Chenji et al., 2016) showed that increased levels of DMN connectivity in patients with greater disability and that this connectivity decreased when performing goal-oriented tasks; consequently, DMN connectivity could be viewed as a baseline functional network. We considered that better swallowing function could lead to reduced levels of DMN connectivity while swallowing. In our present study, we found similar results (Li et al., 2009; Chenji et al., 2016) in that our patients with hemispheric stroke showed improved swallowing function following therapy and that this was associated with reduced brain function network connectivity.

In normal subjects, Babaei et al. (2013) and Malandraki et al. (2009) used brain fMRI to demonstrate increased activation in the cingulate, insula, sensorimotor cortex, prefrontal and parietal cortices, cerebellum, and thalamus during swallowing. Furthermore, Vahdat et al. (2011) observed that alterations in the functional connectivity of the sensorimotor network may be related to motor learning and the importance of the role of somatosensory region in swallowing had been gradually discovered (Babaei et al., 2013; Dehaghani et al., 2016; Mihai et al., 2016). In the present study, stroke patients received a variety of treatments to improve their swallowing, including motor functional training, thermal stimulation, postural compensation, oropharyngeal exercise, compensatory techniques, and swallowing maneuver. We also observed some correlations between clinical swallowing improvements and the SMN network after swallowing therapy in patients suffering from brainstem stroke.

There were several limitations to our study that should be considered when interpreting our conclusions. Firstly, only a small number of patients participated in the study. Secondly, the detailed techniques involved with swallowing interventions were not recorded and compared between hemispheric and brain stem stroke patients. Thirdly, we did not apply longer swallowing interventions for moderate to severe oropharyngeal dysphagia in stroke patients or follow-up the long-term effect of swallowing therapy on swallowing function and neuroplasticity in stroke patients with oropharyngeal dysphagia. In addition, the number of males was 4–5 times as females were enrolled in this study, and the results might be biased by the gender ratio. Finally, we did not recruit stroke patients with oropharyngeal dysphagia in the control group to account for their spontaneous neural recovery, which could have contributed to the alternations in brain fMRI.

#### CONCLUSION

In summary, there was both significant recovery of swallowing function after swallowing therapy for subacute stroke patients with different oropharyngeal dysphagia. The clinical

improvements in FOIS were associated with default mode network in hemispheric stroke and with sensorimotor network in brainstem stroke characterized by fMRI, respectively. For these two dysphagic stroke groups with corresponding etiologies, swallowing therapy could contribute to different functional neuroplasticity, which might serve as an indicator for further prognostic evaluation.

#### ETHICS STATEMENT

Institutional Review Board in Chang-Gung Memorial Hospital approved the study protocol, and all of the participants or their guardians provided written informed consent. All of the participants or their guardians provided written informed consent.

#### AUTHOR CONTRIBUTIONS

Y-CH and W-CL designed the study, participated in data collection, data analysis and interpretation, writing and revising manuscript and final approval of manuscript. T-WH participated in patient enrollment, data collection, data analysis, and data interpretation. H-CH conducted the data statistical analysis, data interpretation, and manuscript drafting. C-PL participated in patient assessment, data analysis and data interpretation. All

#### REFERENCES


authors reviewed the manuscript, contributed to its revision, and approved the final version submitted.

#### FUNDING

This work was supported by funds from the National Science Council (NMRPG896021-3:NSC 99-2314-B-182A-021- 1090 MY3). This study was supported by grants from the Taiwan National Science Council.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Huang, Hsu, Leong, Hsieh and Lin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# **Spatiotemporal Propagation of the Cortical Atrophy: Population and Individual Patterns**

*Igor Koval 1,2 \*, Jean-Baptiste Schiratti 1,2, Alexandre Routier <sup>1</sup> , Michael Bacci <sup>1</sup> , Olivier Colliot 1,3,4, Stéphanie Allassonnière<sup>2</sup> and Stanley Durrleman<sup>1</sup> for the Alzheimer's Disease Neuroimaging Initiative†*

*1 Inria Paris-Rocquencourt, INSERM U1127, CNRS UMR 7225, Sorbonne Universités, UPMC Univ Paris 06 UMRS 1127, Institut du Cerveau et de la Moelle épinière, ICM, Paris, France, <sup>2</sup> INSERM UMRS 1138, Centre de Recherche des Cordeliers, Université Paris Descartes, Paris, France, <sup>3</sup>AP-HP, Pitié-Salpétriere Hospital, Department of Neurology, Paris, France, <sup>4</sup>AP-HP, Pitié-Salpétriere Hospital, Department of Neuroradiology, Paris, France*

Repeated failures in clinical trials for Alzheimer's disease (AD) have raised a strong interest for the prodromal phase of the disease. A better understanding of the brain alterations during this early phase is crucial to diagnose patients sooner, to estimate an accurate disease stage, and to give a reliable prognosis. According to recent evidence, structural alterations in the brain are likely to be sensitive markers of the disease progression. Neuronal loss translates in specific spatiotemporal patterns of cortical atrophy, starting in the enthorinal cortex and spreading over other cortical regions according to specific propagation pathways. We developed a digital model of the cortical atrophy in the left hemisphere from prodromal to diseased phases, which is built on the temporal alignment and combination of several short-term observation data to reconstruct the long-term history of the disease. The model not only provides a description of the spatiotemporal patterns of cortical atrophy at the group level but also shows the variability of these patterns at the individual level in terms of difference in propagation pathways, speed of propagation, and age at propagation onset. Longitudinal MRI datasets of patients with mild cognitive impairments who converted to AD are used to reconstruct the cortical atrophy propagation across all disease stages. Each observation is considered as a signal spatially distributed on a network, such as the cortical mesh, each cortex location being associated to a node. We consider how the temporal profile of the signal varies across the network nodes. We introduce a statistical mixed-effect model to describe the evolution of the cortex alterations. To ensure a spatiotemporal smooth propagation of the alterations, we introduce a constrain on the propagation signal in the model such that neighboring nodes have similar profiles of the signal changes. Our generative model enables the reconstruction of personalized patterns of the neurodegenerative spread, providing a way to estimate disease progression stages and predict the age at which the disease will be diagnosed. The model shows that, for instance, APOE carriers have a significantly higher pace of cortical atrophy but not earlier atrophy onset.

**Keywords: Alzheimer's disease, cortical atrophy, brain networks, spatiotemporal propagation patterns, individual variability**

#### *Edited by:*

*Yasser Iturria Medina, McGill University, Canada*

#### *Reviewed by:*

*Alexandra Young, University College London, United Kingdom Bruno Michel Jedynak, Portland State University, United States*

#### *\*Correspondence:*

*Igor Koval igor.koval@icm-institute.org*

*† Data used in preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at: http://adni.loni.usc.edu/wp-content/ uploads/how\_to\_apply/ADNI\_ Acknowledgement\_List.pdf*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 28 July 2017 Accepted: 26 March 2018 Published: 04 May 2018*

#### *Citation:*

*Koval I, Schiratti J-B, Routier A, Bacci M, Colliot O, Allassonnière S and Durrleman S (2018) Spatiotemporal Propagation of the Cortical Atrophy: Population and Individual Patterns. Front. Neurol. 9:235. doi: 10.3389/fneur.2018.00235*

# **1. INTRODUCTION**

Neuroimaging studies have shown an alteration of the brain structure during the course of Alzheimer's disease (AD) (1, 2). These lesions appear during the prodromal phase of the disease (3– 5) whose observation have been limited due to the absence of clinical symptoms and diagnosis. The importance of the structural changes before the clinical symptoms led to hypothetical models (6), which have been later refined thanks to the gathering of multiple scientific evidences. These modifications took the form of a structural change of the brain in particular an important neuronal loss and an atrophy of the brain cortex (7, 8). The study of the temporal evolution of the cerebral cortex reveals an atrophy of the gray matter (9). This cortical atrophy presumably relates the traces of the progression of the lesions over the brain surface. A fine-scale modeling of the atrophy propagation is likely to give a wider understanding of the disease evolution, as the structural markers seems reliable to assess the conversion to the AD stage, potentially carrying subtle indicators of the disease progression in early phases.

The spatiotemporal propagation of these alterations encloses two entangled components. On the one hand, the spatial characterization of the lesions over the brain surface at each time, and, on the other hand, a temporal dynamic of these alterations that may differ from one region to another. Characterizing the proper dynamics of these lesions relies on the possibility to reconstruct the whole time-line of AD, at both a spatial and temporal level, out of short-term observations that are not temporally aligned. Another challenging aspect consists in the variability inherent to the individual patterns of atrophy that requires to consistently compare the subject-specific spreads of alterations. Accounting for the interindividual variability in term of lesion propagation should allow to reconstruct individual patterns of propagation, paving the way to possible personalized model of atrophy, that potentially informs on subject-specific age of conversion or disease stage.

Recently, large datasets have opened the opportunity to investigate data-driven models that have refined and validated these hypotheses to some extend, in particular event-based models (10– 12) that considers the propagation as a series of events, allowing to define a sequence of disease stages. They characterize the overall variability of the events ordering at a population level. However, these models are not well suited to relate for the temporal delays of the alterations at a population level, neither to determine individual cortical atrophy. Multimodal observations, including positron emission tomography (PET) scans, magnetic resonance imaging (MRI) and biomarkers, have been gathered within longitudinal databases, i.e., repeated observations of patients during significant periods of time. The underlying intention is to provide multiple individual snapshots of the disease—patients examined during short-term periods—in order to reconstruct the long-term history of the pathology (13, 14) at a group and individual level. Moreover, it offers the possibility to describe and interpret the observed data contrary to quantiles or percentiles that require arbitrary reference distributions. A challenging aspect of AD patient comparison is the fact that, even though AD is related to age, the latter is not a good proxy of the disease stage (15–17) leaving us without any easy way to align all the individual on the same time-line. In Ref. (18), the authors introduced mixed-effect model that consider each individual trajectory as a variation of a mean scenario of evolution, with a time-warp function that is able to realign the subjects on the same time-line (19). It allows to characterize a spatial and temporal variability of propagation in the sense that it defines a group-average trajectory of propagation with the possibility to reconstruct individual observations thanks to personalized parameters. Nevertheless (18), constrain the model to parallel profiles of progression which does not hold when looking at signals that have various dynamics. Moreover, the model does not take into account the spatial correlations between the data whereas (17), which focus on spatiotemporal patterns of progression for images, exhibited that this led, in the case of a non-linear mixed-effects model, to poor estimations of the subject-specific parameters and individual trajectories.

To account for the spatial structure of the signal, networks have been introduced (20, 21), representing the brain areas as the graph nodes. In this paper, the networks correspond to a graph representation of a signal spatially distributed, namely, the cortical thickness mapped on a mesh representation of the cortex. The node values are the cortical thickness values over time on the related brain area. Extracting and projecting patients cortical thickness on the common mesh allows to compare their atrophy on the same atlas to exhibit similar patterns. As we expect the signal propagation to be spatially smooth with a similar temporal profile of change for neighbor nodes, we consider that a subset of the graph nodes act as control nodes. They define an evaluation function such that the signal at each node is an interpolation of the signal at the control nodes, enabling to smooth the high frequencies (22). The proximity between nodes is defined by the distance matrix which informs on the distance between any pair. Moreover, the number of nodes of this vertex-based graph can be tuned based on the desired application, potentially the same as the resolution of the input data, e.g., a voxel for MRI or PET data.

The aim of this paper is to introduce a model of the cortical atrophy propagation during the long-term course of AD thanks to a graph representation of the neuroimaging data. This model is able to personalize the reconstruction of the propagation to individual longitudinal measurements, allowing to describe the stages of the disease, potentially in the future. The model is described as a general framework for any longitudinal data spatially distributed on a common graph and it is instantiated to exhibit the propagation of the cortical atrophy on the left hemisphere of the brain, across nearly 2000 regions, thanks to longitudinal observations of 154 mild cognitive impaired (MCI) patients that were later diagnosed with AD. While exhibiting an average pattern of propagation, this mixed-effects model allows to reconstruct individual observations through time.

#### **2. MATERIALS AND METHODS**

### **2.1. Sketch of the Method**

Prior to detail our method, we would like to sketch the key ideas and notations of our work to ease and guide the reading. First, we consider *I* patients; each patient *i* is observed *J<sup>i</sup>* times, his *j*th visit Koval et al. Spatiotemporal Propagation of the Cortical Atrophy in AD

being at age *tij*, and each observation led to an MRI scan as shown on the left hand side of **Figure 1**. Segmentation of the cortical thickness, out of the neuroimaging observations, are mapped onto a mesh, as presented on the middle part of **Figure 1**. The last step corresponds to a subsampling process that led to a graph *G* of K nodes, characterized by a distance matrix *D*. At each node, the individual observations define a time-series describing the evolution of the signal through time.

In a second time, we assume that, at each node *k* of the graph *G*, there exists a function *t 7→γk*(*t*) that describes a characteristic evolution of the signal at this node, as shown on **Figure 2**. The timeseries of individual *i* at node *k* derives from a continuous function *ηik*(*t*), which is assumed to be a spatial and temporal variations of the representative trajectory *γk*(*t*), illustrated on **Figure 3**. The temporal variation corresponds to the time realignment of individual *i* on the common time-line. It adjusts the individual dynamics to a mean pace of evolution, thanks to personalized parameters *τ <sup>i</sup>* and *αi*. *τ <sup>i</sup>* stands for the individual time-shift to the mean disease onset, allowing an early (*τ <sup>i</sup> <* 0) or delayed (*τ <sup>i</sup> <* 0) age at diagnosis. The parameter *α<sup>i</sup>* integrates the patient-specific possibility to have a faster (*α<sup>i</sup> >* 1) or slower (*α<sup>I</sup> <* 1) pace of atrophy compared to the mean scenario of changes. On the other side, the spatial variation corresponds to the adjustment from the mean cortical thickness to individual data. It accounts, for instance, for the difference in size or in spatial thickness distribution at the same disease stage.

We consider that the characteristic signal *γk*(*t*) at node *k* belongs to a family of curve, here the straight line curves, parametrized by the cortical thickness *p<sup>k</sup>* and the rate of atrophy *vk*. To account for the spatial structure of the signal and the large number of nodes, a subset of nodes, referred to as control nodes, is selected to control the interpolation of the cortical and atrophy values over all the nodes. The distribution of the control nodes depends on the size of the kernel bandwidth such that the kernels densities map almost uniformly the feature space.

**FIGURE 1** | Data preprocessing that projects the cortical thickness of the raw MRI observation (left) on a mesh, namely, the FSAverage atlas constituted of 163.842 nodes per hemisphere (middle) before subsampling it and averaging the signal onto a 1,827-node graph (right).

The model introduces population parameters, that allow to define a characteristic spatiotemporal trajectory of the atrophy, and individual parameters, that not only enable to reconstruct individual trajectories but also permit the statistical study of the distribution of spatiotemporal atrophy patterns. These parameters are estimated thanks to the Monte-Carlo Markov-Chain Stochastic Approximation Expectation-Maximization (MCMC-SAEM) algorithm, which handles non-linear mixed-effects models, with theoretical guarantees and consistent results in practice.

# **2.2. Subjects and Data Preprocessing**

Data used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). The ADNI was launched in 2003 as a public–private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD). For up-to-date information, see www. adni-info.org.

The Alzheimer's Disease Neuroimaging Initiative (ADNI) dataset contains longitudinal MRI data for patients that are, at each visit, either cognitively normal (CN), mild cognitive impaired (MCI) patients, or, AD subjects. We selected all the subjects that presented a monotonous decline from MCI to AD, called the MCI converters, removing those that may convert from AD back to MCI or CN. Although AD patients get through an MCI phase, we could not keep CN to MCI patients as they might just as well convert to another dementia. Also, the patients that underwent from CN to MCI and then to AD are not numerous enough to give robust estimation of early stages (CN to MCI). Thus, we kept only the MCI to AD visits of such patients. Altogether, the paper focuses on 154 MCI patients that represents 787 visits, each individual being examined 5 times on average, from 2 to 7 times.

Each visit led to a T1-weighted MRI acquisition, as shown on the left side of **Figure 1**. The longitudinal pipeline of FreeSurfer (23) was used to extract the cortical thickness of the left hemisphere of the brain which was then projected on a common atlas, namely, FSAverage (24), which is a three-dimensional mesh composed of 163,842 nodes for each hemisphere represented on the central part of **Figure 1**. This common fixed-graph allows to compare the cortical thickness between visits or patients, node to node.

The data acquisition and interindividual alignment led to a considerable noise, especially in terms of variability in the measures for close nodes. To smooth this noise and to reduce the computational time, we subsampled the initial graph into a new graph of 1,827 nodes. To do so, we selected 1,827 nodes uniformly distributed over the whole FSAverage graph; the other nodes were then associated to one of the 1,827 nodes thanks to a geodesic distance *d* on the graph (i.e., the length of the shortest path on the surface mesh between the nodes) using the Fast Marching Algorithm on the mesh (25). Therefore, it constitutes collection of nodes referred to as patches. The value of each node of the subsampled graph is the average value over the corresponding patch, each being constituted of approximately 89 initial nodes of the FSAverage graph. The resolution of this vertex-based approach is lower than the initial one, shown on the right hand side of **Figure 1**, but still holds the brain topology while smoothing part of the acquisition noise. In our case, each observation can be considered as a vector of size 1,827 where the *k*th coordinate is related to the *k*th node of the common fixed-graph *G*. The latter is also described by the distance matrix *D* between the 1,827 nodes. It was obtained using the geodesic distance *d* between the 1,827 nodes on the initial graph FSAverage, whose edges are weighted by a physical length. Finally, for all *i*, *j ∈* {1, *. . .* , 1,827}, we set *Dij* = d(**x***i*, **x***j*) where **x***<sup>i</sup>* and **x***<sup>j</sup>* are two nodes of the graph.

characterizes the individual spatiotemporal trajectory.

In the following, we will present a data-driven model which allow to track the propagation of any signal spatially distributed, supposedly the cortical thickness. We consider a longitudinal dataset **y** = (**y***i,j*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, 1 *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji* of *I* individuals, each patient *i* being observed *J<sup>i</sup>* times during the study at ages (*tij*) <sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji*. We suppose that there exist a common fixed-graph *G* defined by a set *V* = (**x**1, *. . .* , **x***K*) of K nodes and a distance matrix *D* which accounts for the distance between the nodes. Any node **x***<sup>k</sup> ∈* R 3 corresponds to a coordinate of a point in space. Each observation **y***ij* = = (**y***ij*1, *. . .* , **y***ijK*) *∈* R *k* corresponds to the measured signal spatially distributed over the N nodes of *G*, represented by a point in the multivariate space R *k* , schematically represented on **Figure 3A** for *K* = 3, as if there were only 3 vertices in the mesh. Therefore, it defines a network whose nodes are valued with the signal of interest. It follows that the collection (**y***ij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji* of the

observations of a particular subject defines a network that embeds a time-series on each node of *G*, indexed by the patient age at each observation (*tij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji*.

#### **2.3. Model**

#### 2.3.1. From Short-Term Data to Long-Term History

We assume there that the repeated observations of a subject are sampled from a continuous function *t 7→η<sup>i</sup>* (*t*) = (*ηi*1(*t*), *. . .* , *ηiN* (*t*)), where *ηik*(*t*) describes the decrease of cortical thickness of this *i*th individual at vertex *k*, such that

$$\begin{aligned} \forall i \in \{1, \ldots, I\} & \quad \forall j \in \{1, \ldots, J\_l\} & \quad \forall k \in \{1, \ldots, K\} \\ \eta\_{ijk} &= \eta\_{ik}(t\_{\vec{\eta}}) + \epsilon\_{i\vec{\eta}k}, \end{aligned} \tag{1}$$

where *ϵijk ∼ N*(0*, σ*<sup>2</sup> ) corresponds to the model noise, whose variance is *σ* 2 .

The function *t 7→ηik*(*t*) describes the evolution of the timeseries at node *k* for the individual *i*. Thus, the vector function *t 7→η<sup>i</sup>* (*t*) = (*ηi*1(*t*), *. . .* , *ηiN* (*t*)) describes the continuous evolution on the graph for a particular individual, i.e*.*, the spatiotemporal propagation of the signal over the whole brain. It corresponds to a spatiotemporal trajectory in the space of measurements. The trajectory *t 7→η<sup>i</sup>* (*t*) is therefore able to reconstruct the existing observations (**y***ij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji*, defined at the related time-points (*tij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji*, as shown on **Figure 3B**, but also generate an observation at any time t, potentially in the future.

The repeated data of each individual is a particular window in the long-term course of the disease that potentially overlaps with other patients. We aim to realign along a common time-line these short-term sequences by carefully analyzing the spatiotemporal patterns within each short-term snapshot. Nevertheless, to do so, we also need to account for the interindividual variability in cortical thickness measurements and trajectories of propagation across the network. The interindividual variability prevents us from considering any individual propagation as a good representation of the disease evolution.

Consequently, we assume that there exists a mean scenario of propagation, defined by a group-average spatiotemporal trajectory *t 7→γ* <sup>0</sup>(*t*), represented on **Figure 3C**, such that each individual trajectory *t 7→ηi*(*t*) is a temporal and spatial variation of this mean scenario of changes, detailed in section 2.3.2. This typical scenario of change describes the mean pattern of spatiotemporal propagation of the signal and writes *γ*0(*t*) = (*γ*1(*t*), *. . .*,*γK*(*t*)) where for all *k ∈* {1,*. . .*, *K*}, *t 7→γk*(*t*) characterize the typical temporal evolution of the cortical thickness on the brain region related to the node *k*. As represented in **Figure 2**, each node has a different temporal profile of atrophy, accounting for the variation of the cortical thickness over time.

#### 2.3.2. Individual Estimation

Translating the generic framework introduced by Schiratti et al. (18) into this case requires to exhibit individual parameters that characterize the individual spatial and temporal variations to the mean, namely the *space shifting* and the *time reparametrization*.

#### *2.3.2.1. Time Reparametrization*

We introduce a time-warp function *ψ<sup>i</sup>* (t) that corresponds to a time reparametrization that adjust the individual dynamics on a common time-line, which here is the average spatiotemporal trajectory *γ*0. For any patient *i* with observations (*yij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji* at time-points (*tij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ji*, *ψi*(*tij*) = *α*i(*tij–t*0*–τ <sup>i</sup>*) + *t*<sup>0</sup> where *t*<sup>0</sup> is a common reference time of the reparametrization, *α<sup>i</sup>* encodes for the individual pace of propagation and *t*<sup>0</sup> + *τ <sup>i</sup>* describes subjectspecific time-shift to the mean disease onset. As such, if the acceleration factor *α*<sup>i</sup> is greater than 1, it corresponds to a faster pace of cortical atrophy whereas *α*<sup>i</sup> *<* 1 indicates a slower pace of atrophy. In the same way, the larger the value of the timeshift *τ <sup>i</sup>* is, the later the disease occurs. Therefore, it leads to write *ηi*(*t*) = *γ*0(*ψi*(*t*)) + *εij*. It adjusts the pace at which the trajectory is followed for the *i*th individual.

#### *2.3.2.2. Space Shifting*

In the space of measurements R *K* , we consider individual observations and the mean trajectory *γ*0(*t*) as shown in **Figure 3A**. In order to account for the spatial variability of the individual trajectories, we assume that there exists, for any individual *i*, a vector *w<sup>i</sup> ∈* R *K* called the space shift, that characterizes the spatial variations from *γ*0(*t*) to the observations as shown in **Figure 3B**. For any point on *γ*0(*t*), *γ*0(*t*) + *w*<sup>i</sup> is assumed to be on the individual trajectory. Therefore, it is possible to translate all the points (*γ*0(*t*))*<sup>t</sup> <sup>ϵ</sup>* <sup>R</sup> to (*γ*0+ *wi*)*<sup>t</sup> <sup>ϵ</sup>* <sup>R</sup> as shown in **Figure 3C**. This collection defines the individual trajectory *ηi*(*t*). This space shift must be orthogonal to the trajectory as it ensures the identifiability of the model. In fact, if the direction *w<sup>i</sup>* was not orthogonal to the trajectory, then the projection of *w<sup>i</sup>* on the geodesic *γ*<sup>0</sup> would interfere with the individual time realignment induced by the dynamic parameters (*αi*,*τ <sup>i</sup>*).

Using mathematical tools from the Riemannian geometry beyond the scope of this study (26) shows that the *k*th coordinate of the individual spatiotemporal trajectory writes *ηik*(*t*) = *γk*( *wik <sup>γ</sup>*˙ *<sup>k</sup>*(*t*0) + *ψi*(*t*)). As the space shift must be estimated in R *K* , *w<sup>i</sup>* is supposed to be a linear combination of few independent components, in the spirit of independent component analysis (ICA) (27). It leads to consider *A* a *K × N<sup>s</sup>* matrix of *N<sup>s</sup>* independent directions, and (*sij*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, 1 *<sup>≤</sup> <sup>j</sup> <sup>≤</sup> Ns* parameters to estimate. *s<sup>i</sup>* = (*si*1*, . . . , siN<sup>s</sup>* ) *∈* R *<sup>N</sup><sup>s</sup>* correspond to parameters of individual *i* that characterize his spatial variations from the mean spatiotemporal trajectory. The orthogonality condition, mentionned in the previous paragraph, leads to consider a basis (*B*1*, . . . , B*(*K−*1)*N<sup>s</sup>* ) of matrices, whose columns are orthogonal to the direction of *γ*0(*t*), and parameters (*βl*)<sup>1</sup> *<sup>≤</sup> <sup>l</sup> <sup>≤</sup>* (*K*–1)*Ns* such that *A* = ∑(*K−*1)*N<sup>s</sup> <sup>j</sup>*=<sup>1</sup> *βjBj*. This procedure allows to reduce the dimension of the parameters to estimate for each *w*i, from *K* to the chosen number of sources.

It leads to write:

$$\gamma\_{i\bar{j}k} = \gamma\_k \left( \frac{w\_{i\bar{k}}}{\dot{\gamma}\_k(t\_0)} + \alpha\_i (t\_{\bar{i}\bar{j}} - \tau\_{\bar{i}} - t\_0) + t\_0 \right) + \mathfrak{e}\_{i\bar{j}k}.\tag{2}$$

#### 2.3.3. Curve Parametrization

In this paper, we consider a *straight line model* such that *γk*(*t*) = *vk*(*t–tk*) + *pk*, *v<sup>k</sup>* accounting for the ratio of atrophy and *p<sup>k</sup>* for the thickness value at time *tk*. A linear decay in cortical atrophy is then represented by a straight line trajectory, parametrized by time, in the *K*-dimensional space as shown on **Figure 3**. Note that as shown on **Figure 2**, it is possible to parametrize the same curve with two distinct sets (*p*1, *t*1) and (*p*2, *t*2) preventing from having an identifiable model. We decided to fix the parameter *t<sup>k</sup>* among all the nodes such that for all *k ∈* {1,*. . .*,*K*} *t<sup>k</sup>* = *t*0, the time reference used in section 2.3.2.1, without any loss of generality as *t 7→γk*(*t*) is defined on R. Despite the linear form of each coordinate *t 7→γk*(*t*), the resulting model is non-linear as it includes among others, multiplication of individual and population parameters.

Finally, equation (2) becomes

$$
\rho\_{ijk} = \rho\_k + \omega\_{ik} + \nu\_k \alpha\_i (t\_{ij} - \tau\_i - t\_0) + \mathfrak{e}\_{ijk}.\tag{3}
$$

This model therefore defines a distribution of multivariate straight line trajectories that accounts for the distribution of the individual trajectories.

#### 2.3.4. Spatial Smoothness

The model proposed in this paper deals with data that are spatially distributed on a graph *G* defined by a set of nodes *V* = (**x**1,*. . .*, **x***K*), where each node embeds a spatial coordinate in R 3 . We expect a smoothly varying profile of atrophy across nodes. The proximity between edges is given by the distance matrix *D*.

In order to ensure small variations of the signal, we introduce a subset *V<sup>c</sup>* = (**x***<sup>d</sup>*<sup>1</sup> *, . . . ,* **x***<sup>d</sup>Nc* ) *⊂ V* whose vertices are called control nodes. Instead of estimating (*pk*)<sup>1</sup> *<sup>≤</sup> <sup>k</sup> <sup>≤</sup> <sup>K</sup>* (resp. (*vk*)<sup>1</sup> *<sup>≤</sup> <sup>k</sup> <sup>≤</sup> <sup>K</sup>*) at all the nodes, we consider only the parameters at the control nodes (*p<sup>d</sup><sup>k</sup>* ) 1*≤k≤N<sup>c</sup>* (resp. (*v<sup>d</sup><sup>k</sup>* ) 1*≤k≤N<sup>c</sup>* ). We introduce a estimation function **x** *7→p*(**x**) (resp. **x** *7→v*(**x**)) for all **x** *∈ V* such that, at the control nodes, the function is equal to the parameters: *∀k ∈* {1,*. . .*,*Nc*}, *p*(**x***<sup>d</sup><sup>k</sup>* ) = *p dk* (resp. *v*(**x***<sup>d</sup><sup>k</sup>* ) = *v dk* ). At the other nodes, the function is an interpolation of the parameter value at the control nodes weighted by the distance to each of them. Therefore, the control vertices control the evaluation of the parameters among all the nodes.

We choose a Gaussian kernel *K<sup>b</sup>* as interpolation splines: *∀ x, y ∈ V, Kb*(*x, y*) = exp(*− d*(*x,y*) 2 *b* <sup>2</sup> ) where *d* is the geodesic distance on the mesh and *b* is the kernel bandwidth. This interpolation allows to remove the possible high frequencies, smoothing the signal spatially. Therefore, it leads to write:

$$\forall \mathbf{x} \in \mathcal{V}, \rho(\mathbf{x}) = \sum\_{i=1}^{N\_c} \beta\_p^i K\_b(\mathbf{x}, \mathbf{x}\_{d\_i}) \quad \text{and}$$

$$\forall \mathbf{x} \in \mathcal{V}, \nu(\mathbf{x}) = \sum\_{i=1}^{N\_c} \beta\_\nu^i K\_b(\mathbf{x}, \mathbf{x}\_{d\_i}). \tag{4}$$

The parameters (*β i p*) 1*≤i≤N<sup>c</sup>* (resp.(*β i v*) 1*≤i≤N<sup>c</sup>* ) are the solution of the linear system *p <sup>d</sup><sup>k</sup>* = ∑*<sup>N</sup><sup>c</sup> i*=1 *β i <sup>p</sup>Kb*(**x***<sup>d</sup><sup>k</sup> ,* **x***<sup>d</sup><sup>i</sup>* ) (resp. *v <sup>d</sup><sup>k</sup>* = ∑*<sup>N</sup><sup>c</sup> i*=1 *β i <sup>v</sup>Kb*(**x***<sup>d</sup><sup>k</sup> ,* **x***<sup>d</sup><sup>i</sup>* )).

Given these interpolations, equation (3) writes

$$\mathbf{y}\_{\text{ijk}} = \mathbf{p}(\mathbf{x}\_k) + (A\mathbf{s}\_i)\_k + \nu(\mathbf{x}\_k)\alpha\_i(\mathbf{t}\_{\text{ij}} - \boldsymbol{\tau}\_i - \mathbf{t}\_0) + \mathbf{e}\_{\text{ijk}}.\tag{5}$$

Even though the distance computed for the cortical thickness corresponds to a distance on the brain cortex, it is possible to compute a connectivity distance based on the connectome, or even an appropriate combination of some of these distances. The challenging part is to put into correspondence areas defined by the connectivity matrices and other networks such as the FSAverage Atlas.

The choice of the set *V<sup>c</sup>* of control nodes among the whole set of nodes *V* is mostly determined by the choice of the bandwidth *b*: their uniform distribution is such that there is an approximate distance *b* between them. In the case of the cortical thickness, we have chosen a bandwidth equal to 16 mm which is representative of the spatial variability of the signal.

#### **2.4. Algorithm**

Equation (5) describes a mixed-effects model, introducing population and individual parameter in this high-dimensional non-linear model. We consider that ((*αi*) 1*≤i≤I ,* (*τ* ) 1*≤i≤i ,* (*sij*) 1*≤i≤pI,*1*≤j≤N<sup>s</sup>* ) are random-effects of the model, leading to write *∀ i ∈* {1,*. . .*, *I*} *∀j ∈* {1,*. . .*, *Ns*}:

$$\alpha\_{i} = \exp(\xi\_{i}) \; , \xi\_{i} \sim \mathcal{N}(\mathbf{0}, \sigma\_{\xi}^{2}) , \tau\_{i} \sim \mathcal{N}(\mathbf{0}, \sigma\_{\tau}^{2}) , \text{ and},$$

$$s\_{ij} \sim \text{Laplace}\left(\mathbf{0}, \frac{1}{2}\right) . \tag{6}$$

*α<sup>i</sup>* corresponds to the realization of a log-normal distribution so that it is always positive, preventing the individuals to present an increasing cortical thickness over time. Moreover, the Laplacian distribution of *sij* arises from theoretical considerations as we need the model to be identifiable, i.e., the solution of the problem to be unique. Finally, these random-effects account for the statistical distribution of the individual trajectories. In the following, we consider **z** = ((*αi*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, (*τ* )<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, (*sij*) <sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> I,* <sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> Ns*) as hidden variables.

Given equation (5) and the observations **y**, we would like to estimate the parameters *θ* = (*t*0*,*(*p dk* ) 1*≤k≤N<sup>c</sup> ,*(*v dk* ) 1*≤k≤N<sup>c</sup> ,* (*βk*) 1*≤k≤Ns*(*K−*1) *, σ<sup>τ</sup> , σξ, σ*) as a maximum likelihood estimate (MLE) *θ* \* = argmax *p*(**y**|*θ*). The natural way to perform such estimation in mixed-effects models is the Expectation-Maximization algorithm (28). Unfortunately, the E-step is intractable and it is not possible to sample according to the conditional distribution **ALGORITHM 1** | Estimation of the general and individual cortical thickness decrease with the MCMC-SAEM algorithm.

```
Input: Longitudinal dataset y = (yi,j)i, j of measurement maps, with the
          corresponding ages (ti,j)i, j.
        Initial parameters θ
                             0
                               and latent variables z
                                                        0
                                                         .
        Geometrically decreasing sequence of step-sizes ρ
                                                                  k
                                                                   .
        Sufficient statistics S
                               k
Initialization: set k = 0 and S
                                0 = S(z
                                        0
                                         ).
repeat
      Simulation: foreach block of latent variables zb do
            Draw a candidate z
                                  c
                                  b ∼ pb
                                          (.|z
                                              k
                                              b
                                                ).
            Set z
                  c =
                       (
                         z
                          k+1
                          1
                             , . . . , z
                                     k+1
                                     b-1 , z
                                           c
                                           b
                                             , z
                                               k
                                               b+1, . . . , z
                                                           k
                                                           nb
                                                             )
                                                               .
            Compute the acceptation ratio ω = min [
                                                            1,
                                                                q
                                                                 (
                                                                  z
                                                                   c
                                                                    |y,θk
                                                                         )
                                                                q(z
                                                                   k|y,θk)
                                                                           ]
                                                                             .
      end
      Stochastic approx.: S
                              k+1 ← S
                                        k + ρ
                                              k
                                                [
                                                 S
                                                   (
                                                     z
                                                      k+1)
                                                            − S
                                                                 k
                                                                  ]
                                                                   .
      Maximization: θ
                        k+1 ← θ
                                 ∗
                                   (
                                     S
                                       k+1)
                                            .
      Increment: set k ← k + 1.
until convergence;
output: Estimation of θ*.
          Samples (z
                      s
                       )
                        s
                          approximately distributed following q(z|y,θ
                                                                          ∗
                                                                           ).
```
*p*(**z**|**y**, *θ*). Therefore, we use a stochastic version of the EM algorithm coupled with a Monte-Carlo Markov-Chain method, namely, the Monte-Carlo Markov-Chain Stochastic Approximation Expectation-Maximization (MCMC-SAEM) algorithm that is able to deal with non-linear equations in a high-dimensional setting. The algorithm is proven to convergence (29) if the model belongs to the exponential family. In our case, it corresponds to consider that *p <sup>d</sup><sup>k</sup> ∼ N* (*p, σ*<sup>2</sup> *<sup>p</sup>*)*, v <sup>d</sup><sup>k</sup> ∼ N* (*v, σ*<sup>2</sup> *<sup>v</sup>* ) and *β<sup>k</sup> ∼ N* (*β<sup>k</sup> , σ*<sup>2</sup> *<sup>β</sup>*).

This leads to consider **z** = ((*ξi*) 1*≤i≤I ,* (*τi*) 1*≤i≤I ,* (*si*) 1*≤i≤I ,* (*p dk* ) 1*≤k≤N<sup>c</sup> ,* (*v dk* ) 1*≤k≤N<sup>c</sup> ,* (*βk*) 1*≤k≤Ns*(*K−*1) ) as the extended hidden variables and *θ* = (*t*0*, p, v,*(*βk*) 1*≤k≤Ns*(*K−*1) *, σξ, σ<sup>τ</sup> , σp , σv, σ*) as the parameters of the model. The latter introduces sufficient statistics *S* of the model that are functions of the observations **y** and latent variables **z**. The aim of such functions is to disentangle the maximization of the parameters *θ* and the simulation of the latent variables **z**.

The pseudo-code of the algorithm, reproduced in **Algorithm 1**, shows the different steps of the optimization until convergence. For further information about the steps of the algorithm, the reader is referred to Ref. (29–31) and references therein.

#### **2.5. Simulation Study**

Since we introduce a new approach to deal with longitudinal data spatially distributed, we performed a simulation procedure to show both the legitimacy of the model used, and the effectiveness of the estimation procedure. To this end, we define a graph represented on the top left of **Figure 4**, representing a square mesh of 7 nodes per edge, thus 49 nodes in total. Among them, 9 equally distributed nodes represent the control nodes, in red on the figure. As we simulate data according to equation (3), we choose position and velocities across the node of the graph, as shown on the top right part of **Figure 4**. Then we simulated realizations (*ξi*,*τ <sup>i</sup>*, (*sij*)<sup>1</sup> *<sup>≤</sup> <sup>j</sup><sup>≤</sup> <sup>N</sup>*s) <sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>N</sup>* for 350 patients, from 4 to 12 visits each (2,980 visits in total) such that it represents 350 longitudinal trajectories of biomarkers spatially distributed. These

**TABLE 1** | The table shows the ability of the algorithm to estimate the real value of the model parameters.


datawere used to find the parameters used to simulate them. Thus, we have performed 10 runs of the estimation procedure. In order to account for the stochasticity of the algorithm and the motion of the Markov Chains, the results in **Table 1** are given with their standard deviation over 10 runs. As we need an initial value for the parameters, we initialized the algorithm without specific knowledge about the positions and velocities, contrary to the experience on the cortical atrophy, so it might reflect a worst-case scenario. **Table 1** shows how well the algorithm performs on either control nodes or random nodes, as well as for the individual parameters.

The bottom part of **Figure 4** shows some results of the estimation procedure. On the left hand side, we provide an example of the stochastic estimations of a parameter over the iterations of the algorithm - the figure shows 10 independent runs. The right hand side presents the final estimation of the velocities across the node of the graph, showing that the model is likely to reproduce the real signal. Overall, these results confirm that such procedure seems reasonable to assess the validity of the model and of the estimation procedure in order to estimate the temporal profile of longitudinal data spatially distributed, such as the cortical atrophy.

# **3. RESULTS**

# **3.1. Initialization**

We evaluated the propagation of the cortical atrophy thanks to cortical thickness values of 154 MCI converters (787 observations) distributed on a graph with 1,827 nodes.

The initialization of the MCMC-SAEM algorithm requires initial values of the parameters *θ* and realizations **z**. We would like to draw attention on the realizations ((*αi*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, (*τ <sup>i</sup>*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, (*si*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*) and ((*p dk* ) 1*≤k≤N<sup>c</sup> ,*(*v dk* ) 1*≤k≤N<sup>c</sup> , t*0). The former are chosen equal to 0, leading to initial individual trajectories that are equal to the mean spatiotemporal trajectories. The pattern of atrophy is the same for everyone at the beginning. The latter variables, ((*p dk* ) 1*≤k≤N<sup>c</sup> ,*(*v dk* ) 1*≤k≤N<sup>c</sup> , t*0), are initialized based on the raw data. Besides *t*<sup>0</sup> that is chosen as the mean age of the input observations, for each control node *k*, we computed linear regressions on the longitudinal thickness values of each patient. Then we average the regression coefficients, each corresponding to a given subject, such that we end up with one rate of atrophy *v<sup>k</sup>* per patch. Also, *p<sup>k</sup>* was chosen as the average thickness on a given patch. **Figure 5** shows the map of the initial *v<sup>k</sup>* distributed over the cortical surface which looks reasonable.

**FIGURE 5** | Annual rate of atrophy mapped over the brain surface used as initialization of our algorithm. Given one area, the corresponding rate of atrophy is obtained as the average regression coefficient of the linear regressions applied to each patient independently.

The initializations of **Figure 5** present areas with important cortical decrease over time, such as the temporal lobe and the hippocampus area. On the other hand, the primary visual cortex is less subject to a cortical atrophy. This initialization looks reasonable; however, these linear regressions are not able to reconstruct the individual observations, preventing from a characterization of personalized patterns of atrophy. It avoids describing the temporal and spatial variability of the individual propagations. Moreover, the linear regressions do not take into account the spatial coherence of the propagation as shown by the color-bar on **Figure 5** where some areas present an important increase of the cortical thickness. It may be associated to the important noise within the data which is produced by the data acquisition, the extraction of the cortical thickness, and, the alignment on the same atlas.

Thanks to the model we introduced, we were able to reconstruct a mean (resp. individual) spatiotemporal trajectory, detailed in section 3.2 (resp. 3.3), that takes the form of the input measurements, preventing from working with percentiles or clusters that cannot be compared directly to the real observations. Due to the numerous number of hyperparameters and the stochastic behavior of the MCMC-SAEM, the algorithm was computed several times, each run of 100,000 iterations taking approximately 15 h. The runs led to similar results. In the following, the results are presented for the run that provided the best individual reconstruction, i.e., the smaller standard deviation *σ* of the noise. Its last estimation is of 0.29 mm, where 90% of the input data are between 1.5 and 4 mm.

## **3.2. Population Level**

The model exhibits a long-term characteristic pattern of atrophy propagation from early MCI stage to post AD diagnosis. It corresponds to the group-average trajectory described in section 2.3.1 whose spatial (**w***i*) and temporal (*α<sup>i</sup>* and *τ <sup>i</sup>*) variations corresponds to individual spatiotemporal trajectories. It is important to mention that this trajectory is a mean trajectory in a statistical sense, as its parameters are the mean values of the individual parameters.

**Figure 6** shows the temporal and spatial evolution of the cortical atrophy, from 66 to 78 years old. The brain medial and lateral views shows an important atrophy on the temporal lobe and the medial temporal lobe, especially the fusiform and the parahippocampical gyrus. An important cortical decrease is also discernible on the superior frontal gyrus and at the wider region defined by the inferior parietal lobe and the angular gyrus. On

the other side, the prefrontal cortex, the primary visual cortex, the calcaris sulcus, and the post central gyrus are less subject to atrophy.

These results are supported by **Figure 7** that shows the map of the annual atrophy *v<sup>k</sup>* for the mean spatiotemporal trajectory, distributed over the corresponding brain areas. The areas affected by the cortical atrophy correspond to previous knowledge (32–34) even tough the different measurements and methodologies lack in consensus. The patterns are still debated in order to find the best characterization of AD compared to normal aging or other neurodegenerative diseases. The proposed model may provide results for different populations on the same atlas, facilitating the comparison between diseases or with normal aging.

## **3.3. Individual Reconstruction**

The model is able to characterize personalized patterns of atrophy propagation thanks to a reconstruction of the individual observations. The validation is assessed thanks to the relative error of reconstruction. As mentioned previously, the input data are noisy, at both a temporal and spatial level. As for the temporal part, the 154 patients represent 281.358 temporal profiles (timeseries) over the 1,827 patches, from which only 6.4% present a monotonous profile of decrease. Given all the linear regression computed for the algorithm initialization, the mean (resp. the variance) of the corresponding R-square values is equal to 0.348 (resp. 0.307). On the other side, the spatial noise corresponds to high variation of the signal for neighbor nodes. Given this important noise, the goal of the reconstruction is not to reconstruct

distributed on the graph.

perfectly the data but rather to smooth the propagation over the brain and to capture individual tendencies of atrophy propagation. Thus, the 787 observations involve 1,437,849 reconstruction *yijk*, whose relative error of reconstruction is represented on **Figure 8A** which confirms the hypothesis that the noise is a Gaussian distribution with a zero mean (*p*-value = 4.24.10–109 for a *t*-test comparison with a theoretical distribution of mean equal to zero). As highlighted by **Figure 8B**, that represents the relative error of reconstruction over the 1,827 patches, the error is mostly randomly distributed over the brain surface. It confirms that the reconstruction error does not have a spatial component as it is uniformly distributed over the brain surface. The color-bar was chosen according to the extreme values: it is important to mention that the larger error of reconstruction corresponds to areas that are close to the corpus callosum where the interpolation relies on

the individual data with the characteristic pattern of atrophy. Moreover, the spatiotemporal trajectory *η<sup>i</sup>* of individual *i* is not estimated only at the observed time-points but it is a continuous function of the time, as shown on **Figure 3C**. Therefore, it is possible to reconstruct the observation at any point, potentially in the future.

One of the properties of the model is to exhibit individual temporal parameters, namely the acceleration factor *α<sup>i</sup>* and the time-shit *τ <sup>i</sup>*, which allow to reparametrize the individual

a fewer number of control points.

**FIGURE 8** | The model is able to reconstruct the data at the individual level, while smoothing the signal over the brain surface, with a relative error randomly distributed. **(A)** Histogram of the relative error of reconstruction of all individuals across all nodes. **(B)** Average relative error of reconstruction over each patch, dynamics on a common time-line. As the data used here correspond to the cortical thickness, the realignment is estimated thanks to structural biomarker dynamics. On the other side, the MCI converters have an age at disease onset, *tdiag*, which corresponds to a clinical status. The latter is not straightforwardly related to the structural dynamics of the individual. In that sense, we decided to realign the age at onset *tdiag*, a clinical biomarker, on the same time-line, assessed with the

**FIGURE 10** | In red, the histogram of the observed age at diagnosis *tdiag,i* for the 154 MCI converters. In blue, the histogram of the repatametrized age at diagnosis *ψ<sup>i</sup>* (*tdiag,i*) once aligned on the common time-line. This shows that the age at diagnosis is mapped to a smaller range of time-points, in the model of cortical atrophy, suggesting that conversion to AD occurs at a specific stage of cortical atrophy.

structural biomarkers. The observed age at diagnosis *tdiag,i* are represented by the red histogram on **Figure 10**, which is not unimodal and present an important variance. The realignment of the clinical status is represented thanks to the distribution of (*ψi*(*tdiag,i*))<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, which is centered with a reduced variance. It suggests that the clinical conversion to AD, determined with *tdiag* corresponds to a specific stage of the cortical atrophy.

As the model estimates individual spatiotemporal trajectories, it allows to describe the variability within the population. The distributions of (*αi*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>*, (*τ* )<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>* and (**w***i*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>* account for the distribution of the individual patterns of atrophy. Furthermore, the ADNI dataset provides, for each patient, multiple features, such as the number of alleles of the APOE-*ϵ*4 gene, the gender, the marital status, and the educational level. In the case of the APOE-*ϵ*4 gene, which is known as a genetic risk factor regarding AD (35, 36), we exhibited the distribution of (*αi*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>* and (*τ <sup>i</sup>*)<sup>1</sup> *<sup>≤</sup> <sup>i</sup> <sup>≤</sup> <sup>I</sup>* for the subpopulations defined by the number of alleles of the gene as shown on **Figure 11**. The more alleles, the more likely to have AD (35, 37).

As shown on the left hand side of **Figure 11**, the patients with two alleles (resp. one allele) present a mean time-shift of *−*2.98 years (resp. *−*0.20 years) after the mean scenario, contrary to patient without APOE-*ϵ*4 alleles that present an average timeshift of 1.89 years, meaning that the more alleles, the earlier the atrophy onset occurs. However, we applied Mann–Whitney two-sided statistical tests that lead to insignificantly differences between the subpopulation. On the other side, same tests were conducted for the same subpopulation with the mean acceleration factor whose distributions are presented on the right hand side of **Figure 11**. In this case, the group of individual with no alleles presented an average acceleration factor of 0.780, statistically different from the group of individual with one alleles (resp. two alleles) that presented an average acceleration factor of 1.415 (resp. 1.236) with a *p*-value equal to 0.00104 (resp.0.00511). However,

this acceleration factor is not statistically different between the population with one or two alleles (*p*-value = 0.51518), meaning that these subpopulation have similar rate of atrophy. Additional investigation on the gender, the marital status and the education level did not led to significant differences. It is important to mention that the Mann–Whitney test is sensitive to the number of samples whereas this study focuses on only 154 MCI patients that might lead to insignificant results in some cases, particularly in the case of the educational level (20 categories) or the marital status (unbalanced classes). Finally, it should mentioned that the tests conducted on the individual space shifts *w*<sup>i</sup> and the related sources *s<sup>i</sup>* did not lead to significant results, mainly because these parameters account for the difference in brain size, and thus thickness, between people.

#### **4. DISCUSSION**

The paper presents a mixed-effects model of the atrophy propagation that is able to characterize a typical pattern of propagation, and, that reconstructs individual observations and scenarios of atrophy. The model exhibits brain areas that are the most affected by the cortical atrophy, such as the parahippocampical gyrus, the temporal lobe, and the superior frontal gyrus. The lesions are less important in the primary visual cortex, the prefrontal cortex, and the primary sensomotory cortex. The model allows to account for the different temporal dynamics of the alterations that can be then compared and ordered.

The proposed model offers a wide versatility of instantiation in terms of profile of temporal variations (exponential decay, sigmoid decay) and spatial variations (resolution, number of control nodes, kernel bandwidth) as it defines a generic framework for the estimation of longitudinal signals spatially distributed. It should be compared to other types of graph-related approaches, such as supervoxels (38) or a vertex-cluster method (39). The latter has exhibited clusters of regression that show profiles of atrophy similar to our results. However, such models do not deal with individual characteristics neither directly with imaging data but rather with normalized values or percentiles, which restrict the interpretation. Further efforts should be concentrated on the validation and improvement of our model, possibly with more complex data and signal propagation.

The individual reconstructions also inform about subjectspecific patterns of atrophy propagation, with potential personalized estimation of the cortical atrophy at future time-points. Further investigations have to be conducted to ensure the quality of the new observations the model is able to generate, so that one can exploit the outcome that the model can predict for an individual some years after his of her last visit. This should be done with a proper validation set to determine the population parameters, and a test-set to predict the individual parameters and thus the future observations. Consistent results might provide information about the structural biomarkers related to the progression of AD, such as in Ref. (40).

Another improvement of the model relies in the distance matrix computation. In this paper, the distance between the nodes is related to the distance on the brain surface, hiding potential effects of the neuronal connections. New distances might be computed based on functional connectivity or combination of different distances, in order to associate the functional and structural components of the brain that are supposed to be complementary in the disease process (41–43).

The model has the potential to exhibit the spatiotemporal propagation of any signal spatially distributed over a graph. It can be used in order to compare the patterns of propagation in distinct population, e.g., normal aging or any other neurodegenerative diseases. It is also a first step to define personalized patterns that would help for a future prognosis of the patient stages.

# **AUTHOR CONTRIBUTIONS**

All authors listed have made substantial, direct and intellectual contribution to the work, and approved it for publication.

# **ACKNOWLEDGMENTS**

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 Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer's Association; Alzheimer's Drug DiscoveryFoundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.;Cogstate;Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann- La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.;Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies;Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www. fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

# **FUNDING**

This work has been partly funded by the European Research Council (ERC) under grant agreement No. 678304, European Union's Horizon 2020 research and innovation program under grant agreement No. 666992, and the program "Investissements d'avenir" ANR-10-IAIHU-06.

# **REFERENCES**


progression. *Neurobiol Aging* (2015) 36:S23–31. doi:10.1016/j.neurobiolaging. 2014.04.034


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Koval, Schiratti, Routier, Bacci, Colliot, Allassonnière and Durrleman. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Mathematical Modeling of Protein Misfolding Mechanisms in Neurological Diseases: A Historical Overview

#### *Felix Carbonell1 \*, Yasser Iturria-Medina2,3 and Alan C. Evans2,3*

*1Biospective Inc., Montreal, QC, Canada, 2Department of Neurology & Neurosurgery, McConnell Brain Imaging Centre, Montreal Neurological Institute, Montreal, QC, Canada, 3 Ludmer Centre for NeuroInformatics and Mental Health, Montreal, QC, Canada*

Protein misfolding refers to a process where proteins become structurally abnormal and lose their specific 3-dimensional spatial configuration. The histopathological presence of misfolded protein (MP) aggregates has been associated as the primary evidence of multiple neurological diseases, including Prion diseases, Alzheimer's disease, Parkinson's disease, and Creutzfeldt-Jacob disease. However, the exact mechanisms of MP aggregation and propagation, as well as their impact in the long-term patient's clinical condition are still not well understood. With this aim, a variety of mathematical models has been proposed for a better insight into the kinetic rate laws that govern the microscopic processes of protein aggregation. Complementary, another class of large-scale models rely on modern molecular imaging techniques for describing the phenomenological effects of MP propagation over the whole brain. Unfortunately, those neuroimaging-based studies do not take full advantage of the tremendous capabilities offered by the chemical kinetics modeling approach. Actually, it has been barely acknowledged that the vast majority of large-scale models have foundations on previous mathematical approaches that describe the chemical kinetics of protein replication and propagation. The purpose of the current manuscript is to present a historical review about the development of mathematical models for describing both microscopic processes that occur during the MP aggregation and large-scale events that characterize the progression of neurodegenerative MP-mediated diseases.

Keywords: misfolded protein, prion-like hypothesis, mathematical modeling, neurodegeneration, therapeutic interventions

# INTRODUCTION

Proteins, large molecules composed by amino acids, play a central role in biological processes and constitute the basis of all the living organisms. During different conformational phases of the proteins, their folding into compact three-dimensional structures is a remarkable example of biological self-assembly and complexity (1). Only correctly folded proteins have long-term stability in crowded biological environments, while a folding failure is traditionally associated with a variety of pathological conditions (1). Proteins that fail to configure properly are called misfolded proteins (MP). In particular, these are thought to disrupt the function of cells, tissues, and organs (2, 3) and

#### *Edited by:*

*Justin John Yerbury, University of Wollongong, Australia*

#### *Reviewed by:*

*Laura Ferraiuolo, University of Sheffield, United Kingdom Michele Vendruscolo, University of Cambridge, United Kingdom*

> *\*Correspondence: Felix Carbonell felix@biospective.com*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 20 July 2017 Accepted: 16 January 2018 Published: 02 February 2018*

#### *Citation:*

*Carbonell F, Iturria-Medina Y and Evans AC (2018) Mathematical Modeling of Protein Misfolding Mechanisms in Neurological Diseases: A Historical Overview. Front. Neurol. 9:37. doi: 10.3389/fneur.2018.00037*

have been causally related to multiple conformational disorders, such as Prion diseases, Alzheimer's disease (AD), Parkinson's disease (PD), Creutzfeldt–Jakob disease, amyotrophic lateral sclerosis (ALS), and several other human degenerative disorders (see **Figure 1**) (1, 4–7).

Despite the biological importance of its negative effects, the mechanisms underlying MP seeding, aggregation, propagation, and/or effective toxicity spreading are not totally understood (15–17) and have been the subject of scientific controversy for decades (14, 18, 19). The high complexity of the underlying processes, as well as the difficulties to extrapolate their effects from microscopic (e.g., molecular) to macroscopic (e.g., organs) scales, have become central obstacles toward the identification of conformational disease-specific triggering events. As a consequence, discrete advances have occurred in the development of effective therapeutic interventions.

For several decades, we have seen the emergence of different mathematical approaches aimed to complement our understanding of the biological mechanisms that lead to MP-related diseases. In a broad manner, mathematical models can be methodologically categorized into two different classes. Namely, a large class of models designed to reproduce molecular-level processes (e.g. seeding, aggregation, short-range spatial spreading in any biological tissue), and a relatively small class of models that account for inter-regional macroscopic interactions (e.g., long-range MP propagation in the human brain). The former class of models have been traditionally developed within the field of chemical kinetics, while the latter one is almost restrictive to neuroimaging studies. Unfortunately, both research fields tend to follow divergent paths. It is not hard to note that neuroimaging-based models invariably follow a single (large-scale) perspective, barely referring to pertinent models coming from the chemical kinetics field.

Motivated by the current lack of an integrative methodological perspective, in this article, we provide a comprehensive historical overview of mathematical models aimed to characterize MP seeding, aggregation, and propagation processes. The manuscript is organized in three sections. The first section offers an overview of models describing prion dynamics, particularly seeding, aggregation, and short-range spatial spreading processes. The second section reviews recent advances on the modeling of prion-like dynamics associated with neurological disorders. The third section highlights important aspects on modeling strategies for the development of drugs and therapeutic interventions in neurological diseases.

#### PRION DYNAMICS

#### Prion Aggregation: One-Dimensional Models

Pioneering mathematical models on prion dynamics (20) mainly focused on simulating the biological mechanisms of prion replication/aggregation previously described in Ref. (21, 22). By using ordinary differential equations (ODEs), Eigen (20) proposed a prototype model for the autocatalysis mechanism of Prusiner (21, 23), which intended to explain the conversion of a typical prion protein (denoted by PrPC) into an infectious agent protein (usually denoted by PrPSc). The autocatalysis mechanism, usually called either the template-assisted model or the heterodimer model (24) (see steps 3 and 4 in **Figure 2**), features a low rate production of PrPSc by spontaneous conversion of PrPC. The PrPSc protein then catalyzes the conformational change of PrPC by forming a heterodimer, which in turn, dissociates at faster rate into two molecules of PrPSc. However, the steady-state analysis (20) of the system describing this mechanism revealed unrealistic requirements for the rate constants. Consequently, the heterodimer model would be, for instance, unable to replicate the long incubation period typically found in prion diseases.

In the same study, Eigen (20) extended the heterodimer model to a framework of a cooperative catalysis mechanism. In this new scenario, a threshold effect on the concentration of PrPSc would determine whether the system maintain the healthy PrPC state (i.e. stable concentration of PrPC) or turn (by an autocatalysis mechanism) into a state of exponential production of PrPSc. Despite that extended model allows for a more meaningful values range for the rate constants, it is still unable to simulate realistic scenarios of prion replication. In order to overcome this limitation, Eigen (20) showed the need to suppress the linear autocatalysis component from the model. In fact, another disadvantage of the cooperative catalysis model is its tendency to generate dynamics of globular aggregates of PrPSc (26). This behavior clearly contrasts with the observations of fibrils of pathogenic PrPC actually conforming geometrically linear on a macroscopic level (22). Nevertheless, with the aim of keeping the argument of a threshold effect, Eigen (20) also simulated the Lansbury's mechanism (22) of plaque formation. In this mechanism, prion replication does not explicitly require a catalysis process but relies on nucleated polymerization. In this way, the steady-state analyzes in Ref. (20), in conjunction with the *in vitro* models developed by Harper and Landsbury in Ref. (27), yield a mathematical formalization of the prion replication mechanism that is currently known as the nucleated polymerization model (NPM) (26, 28).

In summary, the NPM is characterized by four key aspects: (i) no replication occurs below a critical threshold of protein concentration (i.e., polymerization is very slow below a critical size); (ii) a lag time before polymerization for protein concentrations just above the critical threshold; (iii) a relatively rapid polymerization for protein concentrations well above the critical size; and (iv) the slow nucleation process can be bypassed by the introduction of exogenous nucleus or seed (27, 29). A formalization of a mathematical model describing the dynamics of the NPM is due to Nowak et al. (28). In particular, that model describes how PrPSc aggregates can either elongate by the incorporation of a PrPC monomer or break into two new aggregates. These mechanisms of elongation and fragmentation of PrPSc are typically expressed by an infinite set of ODEs (see Eq. A1 in Appendix in Supplementary Material). As noted by Nowak et al. (28) and Masel et al. (26), that system of infinite ODEs can be closed by summation to yield an ODEs system of only three equations (Eq. A2 in Appendix in Supplementary Material). Remarkably, Nowak et al. (28) realized, for the very first time, the analogy between this model and popular epidemiological models for describing virus dynamics (30).

Figure 2 | Schematic representation of prion aggregation mechanisms. Native (sphere) prion molecules undergo conformational changes that lead to an abnormal (cube) configuration (Step 1). This event is unfavorable because the abnormal configuration is either unstable (Step 2) or sensitive to clearance. According to the template-assisted model (24), prions in their abnormal configuration interact with native prions (Step 3) and convert them into the abnormal configuration (Step 4). The NPM proposes that abnormal prions can interact with molecules in a similar state (Step 5), the oligomeric species formed are unstable and grow by the incorporation of abnormal prion molecules (Step 6) and dissociate (Step 7) until a stable nucleus is formed. Such a stable prion aggregate can then grow indefinitely from one or both ends and can also break into smaller fragments (Step 8) that act as new nuclei. Figure reproduced from Brundin et al. (25) with the permission of the journal.

A detailed analysis [see Ref. (28)] of the Eq. A1 revealed a mechanism whereby aggregation could be initiated from a PrPSc monomer at a uniform rate, thus neglecting the (highly probable) possibility of aggregation from a nucleated seed. Based on this observation, Nowak et al. (28) introduced a modified model that allows for a more realistic assumption: all PrPSc aggregates below a critical threshold size would become unstable and quickly dissociate into pieces that return to a PrPC state. Therefore, a realistic aggregation mechanism would require a nucleation-based seed of at least a critical threshold size. Unfortunately, that modified model [see Appendix B in Supplementary Material in Ref. (28)] was not closed by summation as in the case of Eq. A1 and required further approximations. Nevertheless, Masel et al. (26) extended Nowak's NPM for covering the possibility of a nucleated seed [see Eq. 8 in Ref. (26)] while also reducing the number of kinetic rates and making the system closed by summation. As in the case of Eq. A2, the Masel's NPM can be also described by an epidemiological-like system of three ODEs (see Eq. A3 in Appendix in Supplementary Material). In comparison with the original formulation of Nowak et al. (28), the Masel's model not only provided important mathematical simplifications but also facilitated the validation of the NPM dynamics with *in vitro* data.

Several modifications and extensions to the NPM has been proposed during the last two decades (31–36). A remarkable one was given by Greer et al. (33), which, rather than modeling the PrPSc dynamics through an infinite set of differential equations, considered a continuum of possible fibril lengths, written down by a transport partial differential equation (PDE). The PDE formalism turned out to be more accessible from the mathematical point of view and yielded three-dimensional epidemiological-like models. In addition, the PDE formalism enables the study of different aspects of prion dynamics, such as the distribution of prion fibrils and fragmentation processes (33). By relying on the Greer's approach, Prüss et al. (34) provided a detailed characterization of the epidemiological-like behavior of the NMP. By making an analogy with the basic reproductive ratio (30) used in epidemiological theory, Prüss et al. (34) defined a constant *R0* as the number of secondary infections produced on average by one infectious prion (i.e., PrPSc). Thus, if *R0* ≤ 1, the prion replication stops and the disease-fee equilibrium state turns out globally asymptotically stable. On the contrary, if *R0* > 1*,* the prion replication persists as an asymptotically stable disease state.

As in the case of the discrete formulation of the NPM, some studies (37, 38) have focused on extending the Greer's continuous fibril length modeling (33) to more general frameworks (e.g., size-dependent kinetic parameters). Despite the mentioned advantages of the continuous formulation, some efforts have been put on investigating its mathematical connection with the discrete model, as well as the biological implications of such generalization (39). A more recent study (40) has demonstrated that the discrete formulation of the NPM can still be exploited in order to gain further insights about the dynamics of prion replication/ aggregation.

#### Prion Propagation: Spatial Spreading

The main motivation for considering neuronal transport of prion proteins seems to be given by early observations in animal models (41) about the latency between the appearance of PrPSc in the peripherical nervous system and its appearance in the central nervous system. This latency cannot be explained by the rate at which PrPC is converted into PrPSc at specific spatial location, but it is more likely related to the rate of spread between neighboring localities (42). The first attempt to simultaneously model processes of prion aggregation and neuronal transport inside the brain was given by Payne and Krakauer (42). Indeed, Payne and Krakauer (42) extended the template-assisted model (20, 23) of prion replication by adding a spatial spreading component consisting on the statistical assemblage of PrPSc molecules with a classical diffusion process. Note that, with the introduction of this spatial spreading component, Payne and Krakauer (42) overcame the main limitation (i.e. no disease latency) of the templateassisted model. However, due to the acceptance of nucleated polymerization as a plausible mechanism for prion aggregation, incorporation of spatial spreading into the NPM seems to be a promising choice. With this aim, Matthäus (43) used the diffusion approach within the discrete NPM in order to model prion spatial spreading over simple (i.e., one-dimensional) domains like the spine or the mouse visual system (see Eq. A4 in Appendix in Supplementary Material). However, classical diffusion is better suited for modeling free movement of particles in homogenous medium, while the whole brain is a highly heterogeneous one. By acknowledging this limitation, Matthäus (43) proved that the solutions of the NPM model with one-dimensional diffusion follow a traveling wave behavior. Despite such characterization provided a clearer interpretation about the prion propagation speed along pre-defined spatial domains, the proposed model cannot be reduced to a simpler epidemiological-like system. Besides, as pointed out by Matthäus in Ref. (43), it is not realistic to extend the isotropic diffusion approach to large multi-dimensional domains. Indeed, modeling diffusion in a homogeneous medium would force prions to spread with equal speed in all spatial directions. This is unlikely to happen in practice since *in vitro* models have proved that prions spread along the neuronal pathways (44), where infection may reach distant cells at the same time or even faster than neighboring cells (43).

By using a simplified approach, Stumpf and Krakauer (45) modified the epidemiological-like configuration of the NPM in order to incorporate spatial "connectivity" features into the temporal evolution of prion kinetics. Specifically, Stumpf and Krakauer (45) used a lattice domain to model the influence of cell connectivity and cell density in several prion diseases (see Eq. A5 in Appendix in Supplementary Material). The main assumption in Stumpf–Krakauer's model (45) is that PrPSc components spatially spreads along axons and dendrites by slow axonal transport, where the rate of spread from cell to cell depend on the connectivity strengths*.* To our knowledge, Stumpf and Krakauer model constitutes the first successful attempt of incorporating brain connectivity features into prion diseases models described by epidemiological-like equations.

Summarizing, Matthäus (43) and Stumpf and Krakauer (45) were the first studies that attempted to fill the gap between two different modeling scales: reaction kinetics at molecular level and spatial spreading along a large (e.g., whole brain) and complex domain. Unfortunately, these pioneering studies have received little acknowledgment in subsequent studies about proteins propagation over large-scale brain networks. Due to its close relationship with current large-scale protein propagation models, we will provide more details about the network approach over the following sections.

# PRION-LIKE DYNAMICS IN NEUROLOGICAL DISEASES

Similarly to prion diseases, several neurodegenerative diseases (e.g., AD, PD, and FTD) are pathologically associated with the presence of MP (e.g., tau, Aβ, α-synuclein; see **Figure 1**) (15, 17, 25, 46, 47). By using *in vitro* models, it has been demonstrated (48–50) that fibril aggregates of α-synuclein, tau and Aβ proteins self-propagate under biochemical mechanisms analogous to those described for prion aggregation/propagation. These observations in conjunction with several *in vivo* animal models (50, 51) established the founding of the so-called prion-like hypothesis of neurodegenerative progression. Under the prion-like hypothesis, the MP "infectivity" propagates from initial seed regions with a relative high concentration of pathogenic proteins to other "non-infected" brain regions. It should not be surprising that the development of mathematical models for aggregation/propagation of Aβ have followed concurrent paths with those of prion evolution.

# Nucleated Polymerization of A**β**

As in the case of prions, the NPM has been accepted as a plausible preliminary mechanism for Aβ aggregation/propagation (22, 27). However, early *in vitro* studies (27, 52, 53) suggested that the actual mechanism of Aβ aggregation might be more complex than the classical NPM. Indeed, Aβ aggregation is a mechanism likely involving the formation of intermediates soluble micelles [also called protofibrils in Ref. (53)], which are in rapid equilibrium with APP monomers. Such interaction yields domains of high local protein concentration that facilitate the formation of fibril nucleus (27, 52, 54–56). As pointed out in Ref. (54), the formation of intermediates micelles had not been previously detected because the methods for quantification of Aβ aggregation at that time (e.g. Turbidity, Thioflavin T fluorescence) were only able to detect large polymeric structures such as Aβ aggregated fibrils. Using a Quasielastic Light Scattering technique, Lomakin et al. (52) proposed an *in vitro* model that facilitates a quantitative monitoring of the kinetics of Aβ fibrils formation, and consequently, the detection of smaller polymeric structures. Based on that study, Lomakin et al. (54) formalized a mathematical model for describing the simultaneous temporal evolution of APP monomers, Aβ micelles and Aβ fibrils. While in the classical NPM the long latency phase was interpreted as the time required for nucleation, the findings of Ref. (52, 54) suggest that it is rather the time required for formation of larger Aβ protofibrils (55–57). Based on this re-interpretation of the NPM, fibrillization of Aβ would be still a nucleation-dependent process that occur under two concurrent nucleation pathways: exogenous seeds and intermediates Aβ micelles.

As in the case of NPM, Lomakin's model involved an infinite set of ODEs, one per each fibril length. That system can be closed by summation to yield a set of four differential-algebraic equations. There, two ODEs correspond to the temporal evolution of the total length of Aβ fibrils and the total concentration of Aβ aggregates, while two algebraic equations relate (by a conservation of mass condition) the number of APP monomers and Aβ micelles with the Aβ fibrils. Although Lomakin's model was able to replicate the temporal evolution of the mass concentration of fibrils and the fibrils length, it did not consider any fragmentation process. Despite this limitation, Lomakin's model constituted a successful attempt to mathematically model detailed mechanisms of Aβ fibrils and intermediates assemblies of different sizes.

Realizing that experiments in Ref. (52) had been conducted under nonrealistic physiological conditions, Murphy and Pallitto (58) carried out a thoughtful *in vitro* study to better characterize the properties of intermediates micelles during the process of fibrils formation. Using light scattering techniques, Murphy and Pallitto (58) were able to monitor the temporal evolution of the average length of fibrils, the average number of monomers in a fibril, as well as to compute the time to appearance of macroscopic Aβ aggregates. Based on this characterization of the Aβ assemblies [e.g., monomers, micelles, filaments (i.e. thin fibrils), (thick) fibrils, and aggregates], Pallitto and Murphy (59) proposed a detailed multi-steps model for Aβ aggregation kinetics. In that model, a nucleation mechanism was not assumed for the initial step of conversion of unfolded monomers into micelles but for the further self-association of micelles into a nucleus (59). This multistep mechanism also includes: (i) a cooperative (i.e., reversible) self-association of micelles into polymeric nucleus, (ii) elongation of nucleus into filaments by aggregation of micelles, (iii) lateral aggregation of filaments into fibrils, and (iv) fibril elongation *via* end-to-end aggregation (57, 59). As is usual in kinetic modeling, this multi-step mechanism translates into an infinite set of ODEs, which in turn, can be closed by summation to yield a set of eight ODEs. Note that the more types of Aβ assemblies included in the model, the greater the number of equations and kinetic constants rates. Fortunately, the master equations approach provides a unified framework for formulating kinetic equations for Aβ assemblies of any size (60, 61).

#### Master Equations for Aggregation: A Modern Approach

The main idea underlying the master equations approach is to use principles of chemical kinetics in order to derive equations that explicitly account for the different microscopic processes involved in the proteins aggregation mechanisms (61). The ultimate goal is, from these (master) equations, to derive integrated rate laws that characterize the kinetics of protein aggregation. Although not explicitly developed in the specific context of Aβ, the master equations approach is currently established as the most general mathematical formulation for describing the kinetics of MP formation (60, 61).

Mathematical modeling of protein aggregation by master equations appeared in pioneering studies of filamentous growth phenomena (62). Under the basic principles of homogeneous nucleation, growth, and dissociation processes, Oosawa and Kasai (62) formulated a master equation for describing the time evolution of a population of filaments with different polymerization numbers. Similar to the NPM, the master equation is expressed by an infinite set of ODEs. It also relates directly to experimental measurements through the number and mass concentrations of the aggregates, which temporally evolve by a closed system of two ODEs, the so-called moment equations. Note that, in contrast to Masel et al. (26), the moment equations do not evolve as an epidemiological-like system. Instead, this system explicitly remarks its dependency on the concentration of free monomers. This discrepancy is given by the fact that, while Masel et al. (26) considered a free production of monomers at a constant rate, the master equations approach imposes the principle of conservation of mass. Thus, in the master equation approach, the time evolution corresponding to the monomer concentration only accounts for monomers consumed by incorporation into aggregates and monomers dissociated from aggregates, while keeping the total amount of monomers at a constant level (i.e., by conservation of mass).

Solving the moment equations yields the desired integrated rate laws that govern the reaction time course [see details in Ref. (61)]. Interestingly, a detailed analysis of the integrated rate law in Ref. (62) revealed that the actual role of the nucleation step was to generate new elongations seeds rather than facilitate the incorporation of monomers into aggregates. Even more, nucleation and growth processes seem to occur simultaneously [see a more detailed discussion in Ref. (61)], thus equally contributing to the length of the well-known lag phase observed on nucleationdependent processes. In addition, for the very particular case of all proteins in initial monomeric configuration, Oosawa and Kasai (62) showed that the early stages of the reaction time course follow a quadratic rate law, which is a feature of growth governed by a primary nucleation pathway. The integrated rate law also revealed the dependency of the reaction time course on a single parameter (in terms of the microscopic rate constants), which ultimately scales with many of the phenomenological observable measurements (e.g. half polymerization time, maximal growth rate) (61). In other words, the integrated rate law of protein aggregation kinetics shows a scaling behavior. This feature has become a very important tool for understanding the protein aggregation process across different temporal and spatial scales [see an excellent discussion about this property in Ref. (63)].

During the last few decades, the Oosawa theory has been extended (60, 64–67) to include mechanisms of fragmentation and heterogeneous (secondary) nucleation (collectively called secondary pathways), where new aggregates could be also created at a rate depending on the concentration of existing aggregates. Extension of the Oosawa's original formulation was mainly motivated by a discovery showing that the mass concentration of actin (65) and hemoglobin S (67) polymers tends to increase more rapidly than a quadratic rate law. As a consequence, the master equation should include extra terms for describing fragmentation processes and a secondary nucleation at a rate proportional to the surface area of existing aggregates (61, 64) (see Eq. A6 in Appendix in Supplementary Material). Similar to Oosawa, Ferrone et al. (66) derived closed equations for the time evolution of the number and mass concentration of the aggregates (see Eq. A7 in Appendix in Supplementary Material). Unfortunately, those equations are not analytically integrable, which poses additional difficulties for the derivation of the corresponding integrated rate law. As an alternative to the classical global analysis of the integrated rate law, Ferrone et al. (66) carried out a perturbation analysis for characterizing the early stages of the reaction time course [see a detailed explanation in Ref. (67)]. Such perturbation analysis revealed that the quadratic early-time growth rate predicted by Oosawa's original theory is still valid when the reaction is dominated by the primary nucleation process. By contrast, when the reaction is dominated by the secondary nucleation process, the early-time growth follows an exponential law.

At that point, the applicability of an integrated rate law including fragmentation and secondary nucleation seemed to be limited to the early stages of the reaction time course. Evidently, an alternative solution would be employing numerical integrators for solving the closed system of equations corresponding to the number and mass concentration of the aggregates. However, it is not recommendable due to the highly non-linear structure of this system. Besides, understanding the role played by the kinetic rates is extremely difficult when only numerical solutions are available. In order to overcome these limitations, Knowles et al. (60) introduced a new technique that has been one of the major contributions to the kinetic theory of amyloid formation. Specifically, Knowles et al. (60) derived analytical solutions for the integrated rate law that extends its validity to the entire reaction time course. The main idea underlying this new technique is to use the earlytimes solution as a starting point in order to solve the non-linear moment equations with an iterative strategy [see details in Ref. (60)]. Thus, closed-form integrated rate laws were presented for the case of fragmentation (60) and monomer-dependent secondary nucleation (64). Similar to the Oosawa theory, these integrated rate laws revealed the dependency of the reaction time course on two parameters that can be easily related to experimentally observed phenomenological variables (60, 61, 64).

Having an analytical formula for the integrated rate law becomes extremely important since it allows global fitting [see an excellent discussion in Ref. (68)] of experimental data under different conditions (e.g., changing monomer concentrations). For instance, *in vitro* data corresponding to the peptides Aβ40 and Aβ42 fit very well to the theoretical model of Knowles et al. (60) which has clearly improved our understanding about the formation of Aβ aggregates (69–71). In fact, the analysis of such experimental data points to the secondary nucleation process as the mechanism responsible for the toxicity related to Aβ42 aggregation (69). Moreover, a similar analysis (70) showed clear differences (in the relative importance of primary nucleation versus the secondary nucleation) between the molecular mechanism of aggregation of Aβ40 and Aβ42.

During the last few years, the master equations approach and the corresponding integrated rate laws have been subject of multiple investigations (72–78). Among them, it is worth remarking the contribution of Cohen et al. (76) which included mechanisms of spatial propagation within the master equations framework. By fitting the model to experimental data, Cohen et al. (76) found that the secondary pathways govern the speed of spatial propagation by diffusion.

To conclude this section, note that, solely from the mathematical point of view, modeling the kinetics of protein aggregation through master equations is still a very active research field (79). Not to mention the profound impact that this novel approach has produced in the quest of therapeutic techniques for reducing the toxicity associated to low molecular weight Aβ aggregates (79).

### Coagulation Theory for Aggregation: A Road to Brain Imaging Modeling

As in the case of the master equations approach, the coagulation theory described by Smoluchowski's equations (80) also covers the general case of self-association among particles assemblies of different sizes. The first references to Smoluchowski's equations in the context of Aβ aggregation appeared in Ref. (58, 59) for describing the axial elongation of fibrils by end-to-end aggregation of shorter fibrils. By using Smoluchowski's equations, Craft et al. (81) proposed a polymerization model where the nucleation process appears implicitly incorporated within the mechanism of association of small size polymers (e.g., monomers, micelles, and filaments). This model includes processes of production and dissociation of monomers as well as elongation and fragmentation processes for Aβ polymers. Similar to Ref. (34) for the case of prion diseases, Craft et al. (81) defined a polymerization ratio *R0* as the product of the production and elongation rates divided by the product of the degradation and fragmentation rates. The Aβ burden (i.e., total number of Aβ molecules) falls into a steadystate level if the polymerization ratio *R0* is less than one, and shows an increasing Aβ accumulation if this ratio is greater than one (81). A more formal presentation of the coagulation theory and Smoluchowski's equations in the context of Aβ aggregation *in vitro* can be found in Ref. (82). (see Eq. A8 in Appendix in Supplementary Material).

An important turning point in the field of Aβ aggregation/ propagation mechanism modeling is due to Achdou et al. (83). That study settled the grounds for linking molecular mechanism of early aggregation/propagation of Aβ oligomers with modern imaging techniques for measurements of amyloid deposition *in vivo*. The mathematical approach followed by Achdou et al. (83) was also based on Smoluchowski's equations. However, rather than writing down closed ODEs for the moments of polymer length distribution, Achdou et al. (83) truncated the infinite set of differential equations (see Eq. A9 in Appendix in Supplementary Material) at a large enough number *N*. Under this approximation, large aggregates composed of more than *N* monomers are not supposed to coagulate each other. Thus, the time evolution equation corresponding to the truncation number *N* should be able to describe the summary of all Aβ assemblies composed of more than *N* monomers. In addition, Achdou et al. (83) realized that the Smoluchowski's equations also provide a straightforward framework for incorporating processes of spatial propagation. Note, however, that Achdou's model is only valid on small spatial domains (e.g., domain size comparable to a multiple of the size of a neuron), for which isotropic diffusion is a valid assumption. Under essentially the same assumptions, a straightforward generalization of Achdou's model was given by Franchi and Tesi (84) which added fragmentation terms to equation (A9). As pointed out in Ref. (83) a major limitation is that this model is only able to simulate temporal trajectories up to the scale of microscopic processes occurring at the single neuron level.

The development of modern imaging techniques demands alternative models and the possibility to probe them at greater macroscopic scales. With this aim, a large-scale model was recently proposed by Bertsch et al. (85). Specifically, Bertsch et al. (85) coupled a set of truncated Smoluchowski's equations to a kinetic-type transport equation that models the spreading of neuronal damage by neuron-to-neuron prion-like transmission. A major advantage of such modeling is the ability to incorporate two different temporal scales evolving over the same spatial domain: a rapid temporal scale (e.g. hours) for modeling microscopic processes of Aβ polymer agglomeration (by Smoluchowski's equations); and a slow (e.g., months, years) scale to account for the longitudinal progression of AD (by the transport equation) (85). In that model, Bertsch et al. (85) assumed that oligomers of length greater than *N* can be measured by neuroimaging techniques (e.g., PIB-PET), as well as the parameter that controls the neuronal damage (e.g., FDG-PET). Note that Bertsch's modeling approach is the first study that attempts to build a bridge between the microscopic and macroscopic processes that characterize the impact of Aβ aggregation on the onset and clinical progression of AD. Although still insufficient, some effort has been already put on checking the mathematical correctness and internal consistency of that model (86, 87).

# The Network Approach: Modeling Inter-regional Propagation

As it was already mentioned, modeling spatial spreading of prion proteins and MP by a homogeneous diffusion process is not a realistic choice in large spatial domains like the whole brain. Indeed, prion proteins and MP can spread long nerves and "infect" distant regions (44). Thus, for the very first time, Matthäus (43) used the so-called network approach for covering this scenario. There, a mathematical representation of a network consisted of a set of nodes and edges, where the nodes represented neuronal cells and the edges characterize whether two cells are connected (e.g., in the form of a synapse) or not (43). Using this simple mathematical framework, Matthäus (43) described the spread of prion protein infection along a network of inter-connected neurons by a discrete Susceptible-Infected epidemiological model (30). In this kind of models, the network nodes are classified into susceptible and infected nodes, where the susceptible ones become infected if at least one of their connected neighbors is already infected. Thus, in contrast to a homogeneous diffusional spread, the network approach models a rapid infection spread within clusters of highly connected neurons, and propagation to other clusters *via* long-distance connections (43).

The network approach (43) does not only model the influence of the network topology on the speed of the MPs spread but it is also flexible enough to handle different spatial scales. Correspondingly, Matthäus (88) formulated a system of reaction–diffusion equations by coupling kinetic equations of the heterodimer model with discrete diffusion terms to account for transport on networks. There, the network nodes can represent distant regions covering large spatial domains, thereby overcoming the limitations of Ref. (43). In this approach, the diffusion term at each node is spatially approximated by a sum of flows among neighboring nodes, such that the prion protein concentration is transported to the neighbors of the node and *vice versa* (88). Unfortunately, the reaction–diffusion equations formulated by Matthäus (88) have not been extended to nucleated polymerization mechanisms of prion replication.

While Matthäus (43, 88) modeled microscopic processes at the neuronal level, more recent macroscopic approaches have focused on the large-scale connectivity of the whole brain. In this line, Raj et al. (89) proposed a macroscopic network diffusion model (NDM). In the NDM, the number of MP afferents from a given brain region to any other region uniquely depends upon the disease concentration factors in both regions and upon the anatomical connection strength between them. This model was initially explored with structural atrophy data (89) and posteriorly with FDG PET metabolism (90), reproducing in both cases characteristic spatial distributions of MP effects on a relatively small sample of late-onset AD subjects. In addition, this diffusion model has been recently extended to account for impulsive sources of brain atrophy patterns over the brain connectivity network (91). As a main limitation, the NDM does not consider mechanisms of clearance and production of MP. Instead, the disease-related factors have no causal interpretation and accumulate gradually in the absence of any source and/or system resistance. By considering a more realistic scenario, Iturria-Medina et al. (92) introduced an epidemic spreading model (ESM) that simultaneously accounts for the regional capacity to produce/clean MPs and the topological information of the brain's anatomical network (see Eq. 10). When applied to the study of late-onset AD using Aβ PET data, that model was able to reproduce Aβ deposition patterns at the individual level. In line with recent experimental results (93, 94), the ESM also identifies that reduced Aβ clearance, and not Aβ overproduction as the primary cause of Aβ deposition. Importantly, as highlighted in the ESM study (92), the cognitive and clinical states of the AD patients can only be partially explained by the mechanisms of Aβ production, clearance, and spatial propagation.

# Beyond MPs: An Integrative Modeling Approach

The existence of detailed pathological mechanisms and hypotheses for AD progression (49, 95–98) has motivated the consideration of a more integrative multifactorial modeling approach for MP formation and propagation (98–102). Early remarkable papers published by Edelstein-Keshet and Spiros (98) and Luca et al. (102) looked into detail at the underlying mechanism of deposition, uptake, removal, and degradation of Aβ. In particular, Edelstein-Keshet and Spiros (98) focused on modeling the interaction among Aβ deposits, glial cell, inflammation, and secreted cytokines, as well as the corresponding stress, recovery, and death of neuronal tissue. The numerical simulations carried out in Ref. (98) have helped to fill the gaps between hypothesized interactions and downstream consequences among different processes occurring during the AD progression. For instance, it was shown that an amyloid clearance deficiency could saturate the system and push it to toxic Aβ levels, yielding a state of competition between protective and toxic effects. Importantly, Edelstein-Keshet and Spiros (98) showed, for the very first time, a mathematical model where severity of AD does not need to correlate with sensitivity of neurons to amyloid fibers. In this way, Edelstein-Keshet and Spiros (98) highlighted the necessity of more integrative mathematical formulations that be able to consider AD and other neurodegenerative diseases as multifactorial, inter-dependent processes.

In the same line, Puri and Li (99) presented a (microscopic scale) mathematical approach for describing the dynamics of critical components in AD pathogenesis. Formulated in terms of ODEs, that model describes well-known interactions among microglia, astroglia, neurons, and Aβ. The main feature there is to use neuronal death as a surrogate for senile amyloid plaque formation. Numerical simulations using that model indicates that inflammation might be used as an early biomarker for AD since microglia-related alterations can occurs long before apparent senile plaques formation (99). Similarly, Hao and Friedman (101) recently formulated a set of ODEs for describing microscopic processes of AD that included neurons, astrocytes, microglias, peripheral macrophages, as well as Aβ and hyperphosphorylated tau proteins. Based on numerical simulations, Hao and Friedman (101) suggested that a combination of inflammatory processes by cytokines and accumulation of Aβ plaques are key elements in accelerating the progression of AD. By following an integrative macroscopic approach, Iturria-Medina et al. (100) proposed a multifactorial causal model (MCM) in order to simultaneously account for macroscopic MP effects, regional multifactorial causal interactions, and pathological propagations through physical networks (e.g., axonal and vascular connectomes). The MCM (100) considers that, once a factor-specific event (e.g., Aβ deposition, vascular dysregulation, and structural alterations) occurs in a given brain region or in a set of regions, it can directly interact at the macroscopic level with other biological factors and alter their states. These alterations can also propagate through anatomical and vascular connections to other brain areas, where similar factor-factor and spreading mechanisms may occur in a positive feedback mechanism (see **Figure 3** and Eq. 11).

#### THERAPEUTIC INTERVENTION MODELING

Usually, chemical kinetic models of proteins aggregation are used as a surrogate to therapeutics interventions *in vitro*. The general idea is to simulate how a therapeutic intervention (e.g., drugs, antibodies, and molecular chaperones) might inhibit certain microscopic aggregation processes. Then kinetic rates of protein aggregation can be estimated and compared under both natural and inhibition conditions. For instance, monitoring kinetic rates as function of a hypothetical drug dose might help to extrapolate small drug dosages inherent of *in vitro* environments to dosages more closely resembling *in vivo* conditions (103).

Within the formalism of the simplest NPM (e.g., Eq. A3), two main strategies have simulated a therapeutic intervention on the kinetics of MPs (103). In this pioneering study, Masel and Jansen (103) used mathematical models to simulate the inhibition of amyloid propagation with three main approaches: (i) by lowering the effective monomer concentration of unfolded proteins; (ii) by blocking growing polymer ends; and (iii) by increasing the polymer breakage rate (see **Figure 4**). They found that therapeutics following the second strategy would provide promising results, while the remaining ones may be ineffective or even accelerate the amyloid formation process at low drug dosages (103). Indeed, any attempt of breaking up protein polymers into smaller pieces might yield undesired effects since small oligomers are more prone to propagate and spread neurotoxicity.

As discussed in a previous section, Craft et al. (81) established a non-linear relationship between the polymerization ratio *R0* and the total Aβ burden, where *R0* determines two different regimes: a steady-state of Aβ burden or a super-critical regime of continuous Aβ burden increase. By using an empirical relationship between Aβ burden and clinical dementia scores (CDR), Craft et al. (81) explored a potential therapeutic treatment based on the reduction of the polymerization ratio. Since the polymerization ratio depends on four different kinetic rates, several treatment strategies could be easily simulated in this context: (i) the reduction of Aβ monomers production rate, (ii) the enhancement of fragmentation, (iii) the enhancement of the clearance (i.e., degradation) rate, and (iv) the reduction of the elongation rate. In fact, around the time that paper was originally submitted [although originally submitted in 2001 (81), only appeared published in 2005], several treatment approaches following some of these strategies appeared promising in preclinical studies (104–106).

Note that the assessment of different treatments scenarios is usually carried out in two main steps. First, one assumes that a hypothetical change (e.g., by treatment) on the appropriate kinetic rates and then substitutes those modified rates into the original model in order to evaluate the post-treatment dynamical states. The therapeutic treatments simulations carried out in Ref. (81) suggested three different possible outcomes: (1) a reduction of Aβ burden from a pre-treatment steady state to a post-treatment steady state; (2) a transition from a pre-treatment Aβ burden increasing state to a post-treatment steady state; and (3) a reduction in Aβ burden increasing from a pre-treatment super-critical state to a post-treatment super-critical state. Importantly, as pointed out by Craft et al. (81), the failure of a potential drug to reduce the total Aβ burden may not necessarily be associated to drug inactivity but rather to a late intervention during the super-critical state. Besides, Craft et al. (81) showed that any drug treatment based on clearance rate enhancers might be more effective in reducing the total Aβ burden than those based on polymer fragmentation enhancers.

The previous therapeutic intervention modeling was generalized in Ref. (107) by the simulation of the accumulation and

Figure 3 | Multifactorial causal model. (A) Brain multimodal images and cognitive evaluations. (B) State space vectors (S), characterizing the brain's multifactorial alteration levels with regard to a baseline. (C) Causal multifactorial propagation network capturing the direct interactions among regions (for each biological factor/ imaging modality) or among factors (for each brain region). Diagonal blocks in the matrix correspond to a unique biological factor, with diagonal elements accounting for intra-regional effects and off-diagonal elements accounting for inter-regional alterations spreading across physical connections. Off-diagonal blocks correspond to the direct interactions between two different factors (e.g., glucose metabolism impact on tissue properties, or *vice versa*). (D) System has an output vector (β), representing the influence of the brain's multifactorial state space on the cognitive state. Figure adapted from Iturria-Medina et al. (100) with permission of the journal.

spreading of Aβ among the brain, CSF, and plasma. Note that the proposed compartmental model does not only consider production and degradation of Aβ polymers within the brain but also accounts for sources and losses due to transport to and from the plasma and CSF (107). The numerical simulations carried out in Ref. (107) suggest that potential drugs based on the production of Aβ inhibitors (e.g., by enhancing clearance rate) are likely to reduce Aβ burden in the brain, CSF and plasma. By contrast, any drug based on favoring polymers fragmentation and blocking polymers elongation may also reduce Aβ burden in the brain but may not reduce (or even cause a subtle transient increase) Aβ levels in CSF and plasma. By following similar ideas, Das et al. (108) proposed a two-compartment model for the distribution of the γ-secretase inhibitor between the plasma and the CSF, and its effect on the Aβ concentrations in the two compartments. The steady-state analysis of this model reproduced a primary γ-secretase inhibitor effect that caused a decrease in Aβ concentration in both CSF and plasma. However, the model also captured an overshoot of Aβ in the plasma compartment, which was explained by an off-target effect (attenuation of the Aβ clearance rate) of the γ-secretase inhibitor. Das et al. (108) concluded that any effective Aβ-reducing drug would have to necessarily account for more detailed kinetic mechanisms of Aβ production and clearance.

More recent studies such as Ref. (109, 110) used a stochastic modeling approach for simulating discrete versions of simple ODEs describing MP aggregation processes. Similar to Ref. (81), Proctor et al. (109) showed that a small decrease in the dissociation rate of Aβ monomers is enough to increase the chance of appearance of intermediate Aβ toxic species (e.g., dimers, oligomers). In addition, numerical simulations showed that any potential antibodies therapy against the formation of Aβ dimers would have large benefits as an early intervention strategy. Similarly, Proctor et al. (110) studied a model that accounts for a simultaneous intervention on Aβ and tau pathology. Thus, Proctor et al. (110) was able to show that therapies based on Aβ immunization would not only be able to reduce the amount of senile plaques but also produce small reductions in levels of soluble Aβ species, phosphorylated tau proteins, and neurofibrillary tangles.

Undoubtedly, a renovated interest on therapeutic intervention modeling has been motivated by the recent advances in the

field of chemical kinetics and protein aggregation (111–117). In Ref. (116), Arosio et al. used the chemical kinetic approach to model the interaction between molecular chaperones and different protein species. The main idea there was to identify which microscopic reaction steps (i.e., primary nucleation, elongation, fragmentation, and secondary nucleation) where perturbed by the binding of the molecular chaperones to certain protein species (116). By using the master equations approach and the corresponding reaction profile for the total fibril mass concentration, Arosio et al. (116) compared the estimated kinetic rate constants in the absence and presence of different concentrations of molecular chaperones. Note that this kinetics profile analysis required a new mathematical development [See the Supplementary in the external Ref. (116)] for extending the master/moments equations formalism for modeling the binding between molecular chaperones and different protein species (e.g., monomers, fibril end, and fibril surface). For the particular case of the Aβ42 protein, such analysis revealed that the action of particular molecular chaperone (termed DNAJB6) inhibits the primary nucleation process. By contrast, the presence of another molecular chaperone (termed proSP-C Brichos) produces a reduction in the secondary nucleation rate (116). Note that Arosio et al. (116) also analyzed more complex scenarios where specific molecular chaperones might simultaneously affect different microscopic processes (e.g., elongation and secondary nucleation) that characterize the aggregation of the Aβ42 protein. In summary, Arosio et al. (116) provided a detailed modeling of different combination of mechanisms through which molecular chaperones might suppress amyloid aggregation. These results

open up a research avenue where molecular chaperones and other classes of compound might be used as potential therapeutic agents in MP-related diseases.

Similarly, Habchi et al. (113) used the chemical kinetics approach for developing a rational drug discovery strategy that takes into account the specific microscopic steps in the aggregation of the Aβ42 protein. Analogous to the case of molecular chaperones, a potential drug compound could bind to different species of Aβ42 and selectively affect specific microscopic steps during the aggregation process (113). The strategy then proceeds by monitoring the kinetic profiles of Aβ42 fibril formation in the absence and presence of particular potential drugs. Remarkably, by following the master equations modeling approach, Habchi et al. (113) reported that an anticancer drug (termed bexarotene) disturbs the primary nucleation step in the Aβ42 aggregation, delays the formation of toxic oligomers and completely suppresses Aβ42 deposition. This is a general framework that yields a systematic drug discovery strategy aimed to identify a variety of small molecules (112) and antibodies (115) that not only target the onset of aggregation but also the secondary nucleation step responsible for the proliferation of toxic Aβ42 oligomers.

Unfortunately, large-scale models based on phenomenological imaging-based features of protein aggregation (e.g., models based on the brain network approach) have not yet taken advantage of the outstanding advances on therapeutic interventions using the chemical kinetics approach. In this direction, Iturria-Medina et al. (100) used a theoretical control analysis to predict multifactorial intervention effects required to revert brain biomarkers from an



(*Continued*)

TABLE 1 | Continued


advanced disease stage to a clinical normal state. In particular, Iturria-Medina et al. (100) used a multifactorial causal model as an *in silico* evaluator for comparing the macroscopic effects of multiple possible interventional treatments. Their results provided an efficient ranking of multiple AD interventions, which might explain why recent single-target Aβ-based therapies failed to improve clinical outcomes in AD (118, 119).

#### CONCLUSION

In this paper, we intended to provide an historical overview of the development of mathematical models for aggregation and propagation of MP in neurological diseases. Our main goal was to put in contact the neuroimaging community with important studies of MP chemical kinetics modeling, which have not been traditionally acknowledged but constitute a solid framework for a better understanding of neurological diseases evolution. We have mostly followed a chronological presentation of only those mathematical models that, in our opinion, established turning points from either the mathematical modeling point of view or the ability to describe truly biological processes. As a summary of our presentation, we selected the most important of those contributions and present them in **Table 1**. Note that in the main text, we have barely presented an overview of the mathematical formulation and the corresponding biological interpretation involved in those models. Besides, in order to facilitate our exposition, we have used a unified mathematical notation that, in some cases, might differ from the original formulation (see Appendix in Supplementary Material).

We mainly focused our presentation on those models that simulate microscopic processes of nucleation-dependent mechanisms of MP formation. These kinds of (microscopic scale) models provide a unique theoretical framework for relating microscopic processes to macroscopic kinetic profiles of protein aggregation. Thus, the most accepted model for the formation of protein aggregates relies on a variety of microscopic processes, including primary nucleation, fibril elongation, fibril fragmentation, and secondary nucleation, which are collectively summarized by a macroscopic kinetic profile that follows a characteristic sigmoidal shape. The procedure by which highly reproducible kinetic measurements are fitted to this sigmoidal profile allows for (i) detailed characterization of protein aggregation mechanisms in terms of underlying molecular events and (ii) the development of drugs and early therapeutic interventions that might inhibit some of those molecular events. Among other lessons, we have learned that an increase in the monomer concentration of MPs as well a reduction of the monomers clearance rate yield an increase in the growth rate of amyloid formation. This simple lesson highlighted the importance of systematically incorporating chemical kinetics models into strategies of drugs discovery. In fact, it has been suggested (112) that the failure of current therapeutic strategies against AD can be attributed to a limited understanding of the molecular mechanisms by which the tested compounds interact with different species of protein aggregates.

On the other hand, we also presented macroscopic scale modeling approaches that mainly account for the large-scale connectivity of the brain and the indirect phenomenological mapping of the underlying molecular mechanisms of protein aggregation. In line with the network degeneration hypothesis (25), the macroscopic NDM of Ref. (89) supported that MP propagation follows disease-specific anatomical patterns. Similarly, the ESM of Ref. (92) highlighted the importance of considering intra-regional MP generation/clearance and the inter-regional spreading through the anatomical connections. In addition, by using a multifactorial causal model, Iturria-Medina et al. (100) concluded that late-onset AD it is not caused by a unique dominant biological factor (e.g., vascular or Aβ deposition) but by the complex interplay among multiple relevant biological interactions. Taken together, those large-scale mathematical models point to a lack of an integrative perspective as the main cause for the failure of therapeutic strategies against AD.

Undoubtedly, there is still a gap to fill for properly modeling underlying microscopic processes of protein aggregation and their effects on the progression of the associated neurological diseases, as measured by *in vivo* imaging techniques and the assessment of the patient's cognitive condition. Indeed, despite recent efforts (85), most of the large-scale models for protein aggregation still need to incorporate additional components in order to deal with two different temporal scales and spatial domains. Namely, the small-scale where microscopic aggregation processes occur relatively fast and the large-scale where protein deposits accumulate over a long-time period. We hope that future studies about pharmacokinetic modeling of *in vivo* protein binding using PET imaging might shed light on those unresolved issues and yield systematic drug discovery strategies under a broader, integrative perspective.

## AUTHOR CONTRIBUTIONS

FC and YI-M conceived the review and wrote the manuscript. YI-M prepared the figures. AE revised the manuscript. All authors contributed to constructive discussions during the manuscript preparation.

#### ACKNOWLEDGMENTS

The authors would like to thank Dr. Simone Zehntner and B.Sc. Daniel Marchand from Biospective Inc. for their useful

#### REFERENCES


comments and suggestions during the preparation and revision of this manuscript. Likewise, we would like to thank the suggestions made by the anonymous referees during the review process, which undoubtedly improved our final presentation.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fneur.2018.00037/ full#supplementary-material.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Carbonell, Iturria-Medina and Evans. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Perspectives on How Human simultaneous Multi-Modal imaging Adds Directionality to spread Models of Alzheimer's Disease

*Julia Neitzel1,2, Rachel Nuttall2,3 and Christian Sorg2,3,4\**

*1Department of General and Experimental Psychology, Ludwig-Maximilians-Universität (LMU), München, Germany, <sup>2</sup> TUM-Neuroimaging Center (TUM-NIC), Klinikum rechts der Isar, Technische Universität München (TUM), München, Germany, 3Department of Neuroradiology, Klinikum rechts der Isar, Technische Universität München (TUM), München, Germany, 4Department of Psychiatry and Psychotherapy, Klinikum rechts der Isar, Technische Universität München (TUM), München, Germany*

Previous animal research suggests that the spread of pathological agents in Alzheimer's disease (AD) follows the direction of signaling pathways. Specifically, tau pathology has been suggested to propagate in an infection-like mode along axons, from transentorhinal cortices to medial temporal lobe cortices and consequently to other cortical regions, while amyloid-beta (Aβ) pathology seems to spread in an activity-dependent manner among and from isocortical regions into limbic and then subcortical regions. These directed connectivity-based spread models, however, have not been tested directly in AD patients due to the lack of an *in vivo* method to identify directed connectivity in humans. Recently, a new method—metabolic connectivity mapping (MCM)—has been developed and validated in healthy participants that uses simultaneous FDG-PET and resting-state fMRI data acquisition to identify directed intrinsic effective connectivity (EC). To this end, postsynaptic energy consumption (FDG-PET) is used to identify regions with afferent input from other functionally connected brain regions (resting-state fMRI). Here, we discuss how this multi-modal imaging approach allows quantitative, whole-brain mapping of signaling direction in AD patients, thereby pointing out some of the advantages it offers compared to other EC methods (i.e., Granger causality, dynamic causal modeling, Bayesian networks). Most importantly, MCM provides the basis on which models of pathology spread, derived from animal studies, can be tested in AD patients. In particular, future work should investigate whether tau and Aβ in humans propagate along the trajectories of directed connectivity in order to advance our understanding of the neuropathological mechanisms causing disease progression.

Keywords: Alzheimer's disease, spread of pathology, effective connectivity, metabolic connectivity mapping, simultaneous MR-PET imaging

# INTRODUCTION

Alzheimer's disease (AD) is characterized by the extracellular accumulation of misfolded amyloid-β peptides (Aβ), i.e., Aβ plaques, and intracellular aggregates of hyperphosphorylated tau proteins, i.e., neurofibrillary tangles (NFTs) (1). With disease progression, Aβ plaques and NFTs increase in number, yet following distinct spatio-temporal trajectories as revealed by postmortem

#### *Edited by:*

*Ashish Raj, Cornell University, United States*

#### *Reviewed by:*

*Victor Tapias, Cornell University, United States Jacobus F. A. Jansen, Maastricht University Medical Centre, Netherlands*

> *\*Correspondence: Christian Sorg christian.sorg@tum.de*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 13 October 2017 Accepted: 12 January 2018 Published: 26 January 2018*

#### *Citation:*

*Neitzel J, Nuttall R and Sorg C (2018) Perspectives on How Human Simultaneous Multi-Modal Imaging Adds Directionality to Spread Models of Alzheimer's Disease. Front. Neurol. 9:26. doi: 10.3389/fneur.2018.00026*

neuropathological investigations and molecular imaging (2–4). NFTs first emerge in the locus coeruleus and transentorhinal cortex—around the time when first symptoms arise—which subsequently spread to the hippocampus and other limbic regions before finally emerging in isocortical areas (2). Recent tau positron emission tomography (PET) imaging has largely confirmed this pattern of spread (4, 5). Conversely, several years before first symptom onset, Aβ plaques are initially found in the neocortex and subsequently spread to subcortical brain areas at advanced disease stages (6). This pattern has been essentially replicated by amyloid-PET imaging (3).

Animal models suggest that the spread of these pathologies depends critically on the directionality of neural connections (7). **Figure 1** provides a schematic illustration of these spreading processes. Tau pathology appears to disseminate in an infectionlike or prion-like fashion, whereby a self-propagating "infectious" tau protein emerges in intracellular compartments, spreads along the axon, and trans-synaptically induces pathological changes in nearby normal counterparts (8, 9). Aβ has been suggested to spread in an activity-dependent manner: Aβ aggregates trigger aberrant synaptic activity, resulting in hyperactivity (10, 11), which in turn induces increasing rates of Aβ pathology in remote but directly connected regions *via* axons, likely *via* induced hyperactivity (12).

Neuroimaging has further delineated these spreading pathways in humans (5, 13–16) showing that a regions vulnerability to pathological changes depends on connectivity strength, rather than proximity, to the initially affected areas. Myers et al. (13) found that areas with high functional connectivity (FC) during rest, especially the posterior default mode network (DMN), were associated with higher Aβ burden using a within-subject spatial correlation approach. However, a more direct link between

pathology spread and directed connectivity—as suggested by animal models—have not been established in AD patients due to the lack of methods to identify directed connectivity pathways in humans.

Developments in simultaneous multi-modal imaging now offer new approaches to investigate directional signaling, or effective connectivity (EC), *in vivo*. Specifically, a new measure of intrinsic EC (iEC) has recently been established that exploits the simultaneous acquisition of energy metabolism and FC measures on a hybrid MR-PET scanner (17). In this paper, we discuss how this new approach adds a novel quantitative measure of spreading directionality in AD patients.

#### Metabolic Connectivity Mapping (MCM) Provides a Measure of Signaling Directionality in AD Patients

Numerous studies have used undirected FC, defined as statistical dependencies between the activity signals of two brain regions, to investigate pathways of pathology propagation [e.g., Ref. (16)]. However, correlation analyses do not provide information on the influence that one region exerts over another. To understand the signaling hierarchy across distributed networks of regions, measures of EC, i.e., the directed, causal, activity-dependent relationship between regions, are usually more insightful (18).

A novel approach to identify EC in humans integrates undirected FC with local energy consumption based on simultaneously acquired 18F-fludeoxyglucose (FDG) PET and resting-state functional magnetic resonance imaging (17). This method, called "metabolic connectivity mapping," reveals ongoing or iEC (**Figure 2**). The underlying principle of this method is that most energy is spent on signaling processes, 75% of which is consumed postsynaptically (19). At the macroscopic level, it can be assumed that an increase in local metabolism reflects an increase in afferent EC from source regions. In more detail, the directionality of a single functional pathway linking two regions is investigated by taking the cluster FC time series for one region (the potential seed region), which is correlated with the time series of each voxel in another region (the potential target region), reaping one score of FC for each voxel in the target region. On a voxel-wise

level, these scores of FC are then correlated with FDG activity, which, if correlated, infer that this region is the target of this functional pathway. If the same analysis with the seed and target regions switched also shows a significant correlation between FC and FDG values, this is a bidirectional pathway. This analysis is repeated for all region pairs, resulting in a voxel-wise, wholebrain mapping of EC.

Riedl and colleagues (17) have already applied this method to infer the healthy signaling hierarchy in states of externally directed attention (eyes open condition) versus internally directed attention (eyes closed condition). The authors observed bidirectional communication between early and higher visual areas of occipito-parietal lobes plus top-down signaling from a frontoparietal "dorsal attention" network, independent of condition. As soon as participants opened their eyes, parts of the salience network (including insular and cingulate cortices) exert additional top-down influences on the calcarine sulcus. These data support the idea that MCM reveals dynamically changing signaling pathways and, critically, captures the direction of communication among neural networks.

#### Looking at Other Methods to Infer EC in Humans, Their Application to AD, and Methodological Issues

Other researchers have used statistical approaches to infer EC from undirected fMRI data, including Granger causality mapping (GCM) (20, 21), dynamic causal modeling (DCM) (22), and Bayesian network (BN) learning (23). These methods have been used to investigate changed network dynamics in AD patients, reporting disrupted EC in the DMN, though with certain caveats.

Granger causality mapping is based on the assumption that causes precede, and help to predict, consequences. Vector autoregressive models are used to analyze causal interactions between two brain regions, in which the blood-oxygen-level-dependent (BOLD) signal of one region Y at a particular time is modeled as a linear weighted sum of its own past BOLD fluctuations and that of another region X. Activity in area X is said to "Granger" cause activity in area Y if the past of X contains information that helps to predict the future of Y, over and above the information already in the past of Y itself [for review, see Ref (24)]. Applied to AD, GCM revealed altered EC among DMN regions. While the connection strength to the posterior cingulate cortex was markedly reduced in AD patients compared to healthy controls, the medial prefrontal cortex showed stronger coupling with bilateral inferior parietal regions (25). Contrasting results were found by another GCM study reporting relatively preserved posterior cingulate cortex connectivity in AD patients (26). Disease-related changes in GCM have also been found in other networks besides the DMN, e.g., in the executive control network (27). Notably, several assumptions underlie the application of linear autoregressive models to fMRI. While a detailed account can be found in Ref. (28), the strongest criticism that has been raised concerns spurious "causality" that is in fact the result of naturally occurring time-lags among different brain regions. For example, GCM applied to simulated fMRI time-series data was shown to perform relatively poorly, which "suggests that the directionality results may not be trustworthy" (29).

In contrast to GCM, DCM does not estimate EC directly from the observed activity among different brain regions, but instead infers causality from hidden (unobserved) neuronal states that cause those observations. These hidden states are described in terms of bilinear differential equations, which define how the present state of a particular region influences the dynamics of another under experimental manipulation. In order to infer causal interactions at the neural level, DCM integrates a hemodynamic forward model that describes the transformation from neural activity to the measured BOLD signal. Finally, a Bayesian model selection is used to identify the most likely among a set of competing DCMs by comparing the probability of observing the data under a particular model [for technical details, see Ref (30)]. Up to now, only one research group implemented DCM in AD patients (31). In this work, strength of EC was computed during a simple motor task. Compared to healthy control participants, AD patients had significantly reduced EC between the left and right primary somato-motor cortices. The relative lack of DCM studies in the AD literature might be attributed to some restrictions inherent to this approach. The most fundamental issue is that the assumptions held by the hemodynamic forward model, i.e., the mapping between neuronal activity (hidden states) and measured BOLD response, are most probably violated in AD patients due to the damaged vasculature. In brief, neuronal activity drives vasodilation and thereby increases blood flow, which inflates blood volume and reduces the concentration of deoxyhemoglobin. The latter enters the hemodynamic response equation [for more details, see Ref. (32)]. A growing body of evidence indicates that Aβ not only effects neurons but also cerebral blood vessels (33, 34). Decreased arterial blood flow has been found in healthy old carriers of the APOE ε4 genotype, individuals with mild cognitive impairment and AD patients [reviewed by Zhang et al. (35)]. Consequently, the interpretation of DCM results obtained within the AD spectrum is less straightforward; reduced EC could point to altered neuronal interactions and/or AD-related changes in the neurovascular coupling. Furthermore, parameter estimates are wholly dependent on which set of brain regions are included in the DCM, since it is neither mathematically nor computationally feasible to efficiently search over the full range of all possible regions. Therefore, the resulting patterns of EC are only a parsimonious model of the "real" causal architecture. The problem of missing or novel nodes not considered in the predefined model could be quite serious in AD, where atrophy might profoundly alter inter-regional connectivity (36).

Unlike the aforementioned EC methods, BN approaches aim to train a suitable EC model from the data alone, without the need for prior knowledge and considering the entire brain (23). A BN modal is a directed acyclic (no loops that start or end at the same node) graph that consists of nodes representing neuronal regions and edges that symbolize inter-region connectivity. Conditional probability densities are used to determine the functional network structure. BN-inferred EC patterns of AD patients show a global disruption of connectivity from the hippocampus to other main hubs of the DMN, e.g., to the posterior cingulate and medial prefrontal cortex, while coupling between left and right hippocampi were abnormally increased in patients compared to controls (37). Despite the advantages that BN methods encompass compared to other network modeling techniques, the test results obtained from Smith et al.'s simulation study is not all positive (29). Although BN methods were found to excellently detect network connections, estimated directionality was close to chance performance. One restriction in this respect is that BN cannot model reciprocal connections.

# MCM Can Capture the Spread of Pathology in Whole-Brain and Quantitative Terms

We propose that MCM is a promising new tool that, based on the benefits of multi-modal MR-PET imaging, allows one not only to map iEC changes in AD patients but also to link such changes with pathology spread. Especially in the context of AD, this approach offers several advantages compared to other EC methods. First, MCM is a data-driven or model-free approach which requires comparably little pre-assumptions. In fact, unidirectional as well as reciprocal connections can be captured between any regions spanning the whole brain. This is a favorable property considering that little is known about how AD targets the EC structure of the human brain. Second, MCM is less error-prone to naturally occurring as well as AD-related, inter-regional variations in the neurovascular coupling that cause inhomogeneity in the measured BOLD signal. The reason for that is, EC is not directly estimated from the BOLD response, but from correlating the BOLD time series between distinct regions making it invariant to the signal amplitude. In terms of between-subject variations, Riedl and colleagues (17) showed that MCM can reveal robust and condition-specific changes of EC in a group of healthy participants. Thus, the authors concluded "that the assessment of changes in EC may be more robust to vascular heterogeneity." Third, capturing signal directionality from two imaging modalities with similar voxel size also has distinct advantages regarding sensitivity. Since the data are collected simultaneously and independently from the same patient, preprocessing steps commonly applied before statistical analyses, which spatially distort the data, e.g., spatial normalization and smoothing, can be omitted. Instead, the new approach allows EC mapping in individual subject space and may be even sensitive for single-subject analyses. A final, practical advantage of MCM is that EC can be assessed during the resting state, free of any cognitive demand. Mapping iEC opens up novel opportunities for linking the brain's endogenous signaling hierarchy in AD patients with molecular theories about pathology propagation for which experimental evidence has as of yet been restricted to animal models.

Despite being a highly promising method, it is important to highlight the limitations of MCM. First, there is a large difference between FDG-PET imaging and fMRI in terms of temporal resolution: the former can only acquire one saturated image after a period of 30-min scanning, which can cause problems when analyzed in conjunction with a relatively temporally precise and dynamic measure such as fMRI-based FC (38). It is important to adopt a study design that measures stable FC across extended periods when using MCM, so as to ensure similar time scales across both imaging modalities (17).

Second, vascular heterogeneity in terms of vascular density and cerebral blood flow has been shown to influence BOLD-FC (39–41), which can lead to spatial inhomogeneities in the measured BOLD signals and hence may induce false-positives/negatives in the spatial FDG-FC voxel-wise correlations. However, as mentioned previously, since MCM utilizes FC rather than the BOLD signal directly, concern over this potential limitation is somewhat reduced (17).

Finally, one must keep in mind that MCM can only obtain a proxy of EC, since it uses energy consumption as an indication of signaling direction. Recent studies have shown strong support for the underlying assumption that energy consumption is mostly conducted directly at neurons (42), but the findings for a possible role of astrocytes in glucose uptake suggest that the underlying mechanisms of neuroenergetics may not be so clear cut (43). The BOLD signal is also a proxy measure of neuronal activity, but the neuronal basis of the BOLD signal has been widely supported (44–46). The established drawbacks of PET in terms of resolution and sensitivity and its utility in the study of AD pathology should also be taken into consideration when applying MCM to investigate EC and spread models of AD pathology, which have been extensively discussed in other articles (47–52). Additionally, other multi-modal imaging techniques such as fMRI with MR spectroscopy or flumazenil-PET may also offer interesting insight into AD pathology and FC [for reviews, see Ref (53–55)] but, unlike FDG-PET/fMRI, they do not yet offer the key aspect of directionality of functional pathways, along which animal models have shown amyloid-β and tau pathology to spread (7).

#### Application of MCM in AD Patients and Other Neurodegenerative Conditions

Specific approaches to testing spread models of pathology are outlined below. The general logic of these approaches is to compare maps of pathology characteristics, derived from imaging AD patients, and maps of iEC characteristics and changes in these maps in pre-stage AD patients, such as mild cognitive impairment or subjective cognitive impairment. On the one hand, PETbased pathology imaging has demonstrated significant amounts of pathology in these pre-stage AD patients, on the other hand, FC, which forms the basis of EC is largely preserved, facilitating reliable EC. The ultimate question, then, would be to what extent the pathology patterns can be explained by EC pathways. As a simple example, we suggest that, for a pair of regions sharing intact unidrectional EC and a significant gradient of pathology, some variance in this pathology gradient across patients can be explained by variance in the strength of EC beyond underlying functional or structural connectivity. A further example might be a longitudinal approach, in which the increase of a region's pathology is explained by iEC into this region at the time of first measurement.

Furthermore, the application of MCM to the investigation of other neurodegenerative conditions seems promising. Despite their clinical heterogeneity, many neurodegenerative diseases share a common neurological signature—the misfolding and accumulation of specific proteins. Besides AD, this is the case in Parkinson's disease characterized by α-synuclein; sporadic amyotrophic lateral sclerosis and rare fronto-temporal dementia showing aggregates of TAR DNA-binding protein 43 or in Huntington's disease with huntingtin aggregates. Cell culture and/or animal studies more and more firmly demonstrate that these misfolded proteins share the ability of self-perpetuating neuron-to-neuron spreading, implying that neuronal connections probably play a critical role in disease propagation [see Ref (7, 56, 57), for recent reviews]. First evidence for a directiondependent spreading mechanism have been particularly shown for α-synuclein. Pathological changes in Parkinson's disease appear in a prototypical sequence starting in the lower brainstem and olfactory bulb, from where they proceed to the midbrain and the substantia nigra, before being found in the basal forebrain and ultimately in the neocortex (58). Moreover, α-synuclein's ability to propagate transneuronally along defined neuronal pathways has been confirmed in transgenic mice. After intracerebral injection of brain-derived, pathological α-synuclein, the asymptomatic recipient animals developed Parkinson's diseaselike lesions which were also observed in interconnected regions far beyond the injection sites (59). Estimating direction of neuronal communication in humans by MCM may hence allow testing such an infection-like spreading model in Parkinson's disease patients.

#### REFERENCES


# CONCLUSION

Better knowledge of the mechanisms that cause propagation of Aβ and tau pathology from an initially isolated target site to remote regions of wider brain networks will pave the way for more precise diagnostics and novel treatment strategies. Given the clear predictions of animal models that AD pathology spreads in the direction of neuronal pathways, future research should aim to explicitly test this idea in AD patients. MCM has been demonstrated to be a capable tool for detecting iEC, a proxy for directed connectivity, in healthy participants. Applied to AD patients, this multi-modal imaging approach allows future studies to test whether the spread of tau and Aβ pathology in humans follows the hypothesized trajectories of iEC.

## AUTHOR CONTRIBUTIONS

All authors contributed to the conceptualization and the writing of the article.

#### FUNDING

This work was supported by German Federal Ministry of Education and Science Grant No. BMBF 01ER0803 (to CS) and by the German Research Foundation (DFG) and the Technical University of Munich within the funding programme Open Access Publishing.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer VT and handling editor declared their shared affiliation.

*Copyright © 2018 Neitzel, Nuttall and Sorg. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original*  *author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# The Pathoconnectivity Profile of alzheimer's Disease: a Morphometric coalteration network analysis

*Jordi Manuello1,2†, Andrea Nani1,2,3†, Enrico Premi <sup>4</sup> , Barbara Borroni <sup>4</sup> , Tommaso Costa1,2\*, Karina Tatu1,2, Donato Liloia2 , Sergio Duca1 and Franco Cauda1,2*

*1GCS-fMRI, Department of Psychology, Koelliker Hospital, University of Turin, Turin, Italy, 2 FOCUS Laboratory, Department of Psychology, University of Turin, Turin, Italy, 3Michael Trimble Neuropsychiatry Research Group, Birmingham and Solihull Mental Health NHS Foundation Trust, Birmingham, United Kingdom, 4Neurology Unit, Department of Clinical and Experimental Sciences, University of Brescia, Brescia, Italy*

#### *Edited by:*

*Ashish Raj, Cornell University, United States*

#### *Reviewed by:*

*Gianfranco Spalletta, Fondazione Santa Lucia (IRCCS), Italy Amy Kuceyeski, Cornell University, United States*

> *\*Correspondence: Tommaso Costa tommaso.costa@unito.it*

*† These authors have contributed equally to this work.*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 28 June 2017 Accepted: 21 December 2017 Published: 25 January 2018*

#### *Citation:*

*Manuello J, Nani A, Premi E, Borroni B, Costa T, Tatu K, Liloia D, Duca S and Cauda F (2018) The Pathoconnectivity Profile of Alzheimer's Disease: A Morphometric Coalteration Network Analysis. Front. Neurol. 8:739. doi: 10.3389/fneur.2017.00739*

Gray matter alterations are typical features of brain disorders. However, they do not impact on the brain randomly. Indeed, it has been suggested that neuropathological processes can selectively affect certain assemblies of neurons, which typically are at the center of crucial functional networks. Because of their topological centrality, these areas form a *core set* that is more likely to be affected by neuropathological processes. In order to identify and study the pattern formed by brain alterations in patients' with Alzheimer's disease (AD), we devised an innovative meta-analytic method for analyzing voxel-based morphometry data. This methodology enabled us to discover that in AD gray matter alterations do not occur randomly across the brain but, on the contrary, follow identifiable patterns of distribution. This alteration pattern exhibits a network-like structure composed of coaltered areas that can be defined as *coatrophy network*. Within the *coatrophy network* of AD, we were able to further identify a core subnetwork of coaltered areas that includes the left hippocampus, left and right amygdalae, right parahippocampal gyrus, and right temporal inferior gyrus. In virtue of their network centrality, these brain areas can be thought of as *pathoconnectivity hubs*.

Keywords: brain alterations, coatrophy network, pathoconnectivity hubs, Alzheimer's disease, tauopathy, gray matter atrophy, voxel-based morphometry

#### INTRODUCTION

Widespread alterations of gray matter commonly characterize brain disorders. It has been suggested that neuropathological processes can selectively affect certain assemblies of neurons (1), which typically are at the center of crucial functional networks (1–7). Because of their topological centrality, these areas or network hubs form a *core set* that is more likely to be affected by neuropathological processes (1, 8–16). In particular, neurodegenerative diseases exhibit structural alterations that seem to distribute across the brain in network-like patterns (17, 18). These patterns, which we propose to call *morphometric coalteration networks* or, in the case of gray matter decreases, *coatrophy networks*, can be thought of as a form of pathological anatomical covariance (19, 20) and appear to develop according to the organization of brain connectivity (3, 4, 7). Studies aiming to investigate the networks formed by coaltered cerebral areas in the pathological brain are providing new insight for a better transdiagnostic and neurobiological understanding of the mechanisms at the root of brain disorders (21–23).

This is particularly true in the case of Alzheimer's disease (AD). So far great efforts have been made in order to identify a prototypical pattern of gray matter atrophy due to AD, and to put it into correlation with clinical symptoms (24). It is now known that cortical thinning of specific brain sites can be already detected even before the appearance of the symptomatology and that the atrophy tends to increase when the condition worsens (25). Although the cortical reduction is commonly found in normal aging (26, 27), the pathological fingerprints of AD are mainly observed in a temporoparietal set of brain areas, including hippocampus, entorhinal cortex, precuneus, and posterior cingulate cortex (28, 29). The involvement of these regions has been repeatedly confirmed by meta-analytical studies, which have additionally found the alteration of the right superior frontal gyrus (30). According to Ferreira et al. (31) the left medial temporal lobe is the most impaired area in AD, even in the preclinical phases of the disease, so much so that the impairment of this area can be a good predictor of the clinical worsening of AD. A study of the relationship between the cortical thinning in AD and large-scale structural organization of the brain has revealed that AD reduces both the nodal centrality of temporal and parietal heteromodal association cortices and the positive correlation of thickness values normally found bilaterally between the parietal regions. In contrast, authors reported an increase of positive correlation among brain areas that are part of the default mode network (DMN) (32).

Recently, investigations into the cognitive deficits caused by AD have taken advantage of the methodology of network analysis (33, 34). According to this approach, altered brain areas can be represented by means of a set of nodes, linked together by means of edges representing different statistical values. Studies in this line of research have found that AD increases the correlation between the values of cortical thickness of the fusiform gyrus, temporal pole, parahippocampal gyrus, and cingulum, which are all in proximity to each other. Conversely, a decrease of the correlation has been observed between distant areas (35). Of note, it has been suggested that, by combining different sources of information: (i) large-scale structural networks data, (ii) values of cortical thickness, and (iii) the pace of cortical thinning along time, it could be possible to distinguish patients with AD from healthy controls with an accuracy of 96.1%, as well as predict the conversion of mild cognitive impairment (MCI) into AD 6 months before its clinical onset (36). These studies raise the issue of moving from group analysis to single-subject results, which is an essential aspect when dealing with potential biomarkers for diagnostic purposes and surrogate endpoints for disease-modifying clinical trials. Recent methods of single-subject graph measurements have allowed to link network alterations and cognitive decline. For instance, it has been showed that the more the network becomes disorganized, the worse the clinical condition is (37). Moreover, even in healthy subjects, it has been found an association between Aβ42 CSF low levels and alteration of network properties, which might be interpreted as a very early indication of an underlying pathological process (38). All these results provide evidence that the approach based on network analysis can bring valuable insight to clinical practice (33).

So far, at least four important mechanisms have been proposed to account for the distribution of brain alterations: transneuronal spread, nodal stress, shared vulnerability, and trophic failure (4, 5).

The first mechanism suggests that misfolded proteins (native peptides with an incomplete or incorrect folding, as well as *de novo* polypeptides that become prone to self-aggregation) can diffuse along neuronal pathways (18, 39–41). Increasing evidence indicate that the spread of misfolded proteins presents several similarities to the plasma membrane prion protein intercellular transfer, along axonal fibers, potentially contributing to disease progression (42). This mechanism has been demonstrated in neurodegenerative diseases, such as Alzheimer's, Parkinson's, Huntington's, amyotrophic lateral sclerosis, and tauopathies (43, 44); more recently it has been also generalized to other brain disorders (45).

The second mechanism hypothesizes that the functional stress of the network hubs may result in a greater vulnerability of these areas (1, 4, 14, 46). This susceptibility has been supported in human beings with *in vivo* neuroimaging techniques and voxelbased meta-analyses (14).

The third mechanism suggests that certain brain regions sharing gene or protein expressions may be more vulnerable to neuropathology (4, 47–51), with a potential relationship between gene expressions and connectivity patterns (51, 52).

Finally, the fourth mechanism hypothesizes a disruption in the production of trophic factors, which could bring about the deterioration of neural wiring (4, 5, 53, 54).

If we consider the case of AD, neuropathological signatures, namely amyloid-β (Aβ) plaques and neurofibrillary tangles, are already present in the preclinical phase of the disease, with further spreading during progression. In fact several years before the clinical onset of AD, Aβ, and tau progressively accumulate in the brain with a certain degree of spatial specificity as well as a partial overlap among the two deposits (55). The relationship between tau and amyloid deposits in the cerebral cortex seems to have a hierarchical organization, with tau and Aβ clusters exhibiting distinctive intramodal and intermodal characterizations (56). These findings would support the view of AD as an amyloid-facilitated tauopathy (57). Furthermore, Aβ and tau propagation and the subsequent deposition and cytotoxicity effects appear to occur mainly between anatomically interconnected areas, thus affecting the functional communication among them (58).

The concept of a gradual spread of pathological signs is a crucial aspect put forward by recent theoretical models. Raj et al. (3) have proposed a network diffusion model of disease progression in dementia, according to which the propagation of pathogenic proteins follows the regional concentration gradients under the spatial constraints defined by brain connectivity. Other authors have proposed a stochastic epidemic spreading model to describe intra-brain Aβ propagation and deposition processes, according to which regions with a higher connectivity degree are the main target of Aβ, thus suggesting that brain hubs are the more exposed to the negative effects of these aberrant proteins (40). Finally, in addition to focusing on misfolded proteins and propagation pathways, a further interesting approach suggests the need to investigate the relationship between these two factors (18). This model considers molecular nexopathies as conjunctions of pathogenic protein and brain networks. Key factors are therefore supposed to be structural/functional developmental factors and differential vulnerability of neural connections. Accordingly, long-range axonal connections may be more vulnerable to Aβ, so that functional and structural alterations could occur within the large-scale distributed frontotemporoparietal network, such as the one that supports the DMN processing.

In order to identify and study the coatrophy network of AD, we devised an innovative meta-analytic method for analyzing voxel-based morphometry (VBM) data. This methodology enabled us to address the following issues:


#### MATERIALS AND METHODS

#### Selection of Studies

On March 2017, we performed with the software Sleuth an extensive meta-analytic search in the BrainMap VBM database (www. brainmap.org) (59–61). All the studies that fulfilled the following criteria were retrieved: "Contrast is Gray Matter"; "Context is Disease Effect"; "Observed Changes is Controls > Patients" and "Diagnosis is Alzheimer's Disease." Results were controlled so as to keep only experiments comparing subjects diagnosed with AD against healthy controls. Our search focused on gray matter decreased values only, as the development of AD is strongly characterized by axonal deterioration and neuronal loss that result in brain atrophy (62). Furthermore, thus far just a few VBM studies have investigated gray matter increase in AD, so that these data were not sufficient for obtaining reliable results with our metaanalytical methods.

To ensure a transparent description of the selection process, we followed the "PRISMA Statement" international guidelines (63, 64) (Figure S1 in Supplementary Material). The characteristics of the sample can be viewed in **Table 1**.

## Anatomical Likelihood Estimation (ALE) and the Creation of Modeled Activation (MA) Maps

Voxel-based morphometry data were statistically elaborated with the procedure of the ALE. ALE is a voxel-based meta-analytical technique that models the spatial coherence of different results (101–103). A three-dimensional Gaussian probability distribution is then centered on each focus of every experiment, using the following formula:

$$p(d) = \frac{1}{\sigma^3 \sqrt{(2\pi)^3}} \mathbf{e}^{-\frac{d^2}{2\sigma^2}},$$

in which *d* refers to the Euclidean distance between voxels and the considered focus, while e refers to the spatial uncertainty. The SD can be obtained by means of the full-width half-maximum, such as:

$$\boldsymbol{\sigma} = \frac{\text{FWHM}}{\sqrt{8 \ln 2}}.$$

The combination of these Gaussian distributions constructs a MA map for each experiment. The definite ALE map is finally generated by uniting the MA maps. ALE maps were thresholded at a voxel-level FWD *p* < 0.05, in line with Eickhoff et al. (102, 104, 105). Given a specific threshold for cluster forming, a null distribution of cluster sizes was derived by simulating a long series of experiments with the same characteristics of real data and then by generating an ALE map. The score histogram so obtained was eventually employed to assign a threshold *p*-value.

#### Construction of the Morphometric Coatrophy Network

To identify the distribution of gray matter alterations, we have developed a novel methodology capable of constructing the morphometric coalteration networks associated with brain disorders. Our analysis can in fact discover whether an altered brain area, say A, is statistically related to the alteration of one or more other brain areas (B, C, etc.). Thus, our analysis can construct the morphometric coatrophy network composed of the areas occurring to be altered together and, subsequently, investigate within the coatrophy network (i) how an altered region is statistically associated with other altered regions and (ii) which regions are likely to be involved in a more widespread net of alterations.

#### Node Creation and Labeling

We superimposed the ALE map on the Talairach atlas so as to distinguish automatically the anatomical regions identified by the ALE algorithm. If (at least) 20 voxels of the ALE map were found to be inside a certain area of the atlas, then this area was considered to be altered. We chose this cluster threshold so that less relevant regions could be excluded. We employed a peak detection algorithm to identify the local maxima of the ALE map, and we subsequently selected only those peaks that were greater than the 90 percentile of the value distribution. This set was further reduced by applying a minimum interpeak distance of 10 mm. Finally, we positioned a node, labeled on the basis of the Talairach atlas, in correspondence of every survived peak.

#### Thresholding Values Applied during Nodes Creation and Their Rationale

As described in the previous paragraph, three thresholds were applied during the nodes creation procedure.

The first threshold regulates the minimum number of voxel (i.e., 20 voxels) necessary to consider a brain area as altered. The rationale behind this threshold is to exclude from the coatrophy network nodes representing minimally (or, from a meta-analytical point of view, rarely) altered brain areas, thus improving and simplifying the interpretability of the results without losing highly

#### Table 1 | Selected studies for the meta-analysis.


The Coatrophy Pattern of AD

*Where no information about slice thickness was provided, the voxel-size was expressed. The items are the result of the entire selection process as shown in PRISMA (Figure S1 in Supplementary Material) flow chart.*

**58**

relevant information. However, even considering brain areas in which only one voxel is altered, the results would have not been spurious, since ALE maps were voxel-level thresholded, which implies that each single voxel contains statistically significant information (104) (see Figure S2 in Supplementary Material for the visualization of the network obtained with different threshold values). This choice, however, would have unnecessarily increased the complexity of the coatrophy network.

The second threshold, applied to the peaks-value distribution, allowed us to include in the network only nodes representing those areas for which there is a very high consensus between different experiments (i.e., high ALE value) (104). Even in this case, this threshold could have been removed; all the nodes that can be created with the present methodology represent statistically significant effects, since they can only lie inside the anatomical regions identified by the ALE algorithm, which already has its own statistical thresholding step (see Figure S3 in Supplementary Material for the visualization of the network obtained with different threshold values).

Finally, the interpeaks distance was chosen considering the mean value (10.2 mm; SD = 0.4 mm) of uncertainty in spatial location associated with the reported coordinate discussed in Eickhoff et al. (101).

Therefore, the only effect of those thresholds on our data is to decrease the redundancy of the network, so as to obtain clearer results to be visualized and further analyzed, minimizing the information loss.

#### Coatrophy Distribution

From the set of the nodes as defined in the previous paragraph, we constructed a *N* × *M* matrix or a coalteration matrix, in which each row referred to an experiment, whereas each column referred to a network node. On the basis of a Bernoulli generation data model, we constructed a probability distribution of joint alteration values for each pair of nodes. In other words, for any couple of nodes (*a* and *b*), we were able to describe their four conjoint states of alteration by means of two binary variables: (1) *a* and *b* both altered; (2) *a* and *b* both unaltered; (3) *a* altered and *b* unaltered; and (4) *a* unaltered and *b* altered. Consequently, the following four probabilities were obtained by the frequencies of the different combinations of all experiments:

$$\Theta\_1 = P\left(a = 1, b = 1\right),$$

$$\Theta\_2 = P\left(a = 1, b = 0\right),$$

$$\Theta\_3 = P\left(a = 0, b = 1\right),$$

$$\Theta\_4 = P\left(a = 0, b = 0\right).$$

These formulas refer to the conjoint frequencies of a couple of nodes (*a* and *b*) in all their four possible combinations. **Table 2** shows the marginal probabilities for each couple of nodes.

On the grounds of these four probabilities, we have applied the Patel's *k* index (106)—which has been validated with simulated data by Smith et al. (107)—in order to calculate the degree of coalteration between nodes. This index can measure the probability that two nodes (*a* and *b*) are actually coaltered against the Table 2 | Marginal probabilities between altered and unaltered nodes.


probability that node *a* and node *b* are altered independently of each other. Patel's *k* is calculated as follows:

$$\kappa = \left(\mathfrak{H}\_1 - E\right) / \left( D\left(\max\left(\mathfrak{H}\_1\right) - E\right) + \left(1 - D\right) \left(E - \min\left(\mathfrak{H}\_1\right)\right) \right),$$

where

$$E = \left(\mathfrak{G}\_1 + \mathfrak{G}\_2\right)\left(\mathfrak{G}\_1 + \mathfrak{G}\_3\right),$$

$$\max\left(\mathfrak{G}\_1\right) = \min\left(\mathfrak{G}\_1 + \mathfrak{G}\_2, \mathfrak{G}\_1 + \mathfrak{G}\_3\right),$$

$$\min\left(\mathfrak{G}\_1\right) = \max\left(0, 2\mathfrak{G}\_1 + \mathfrak{G}\_2 + \mathfrak{G}\_3 - 1\right).$$

The numerator refers to the difference between the probability that *a* and *b* are altered together and the expected probability that *a* and *b* are altered independently of each other. The denominator refers to a weighted normalizing constant. Min (ϑ1) refers to the maximum value of the conjoint probability *P*(*a*,*b*), given *P*(*a*) and *P*(*b*), whereas max (ϑ1) refers to the minimum value of *P*(*a*,*b*), given *P*(*a*) and *P*(*b*). Patel's *k* index has values that range from −1 to 1. A value of |*k*| that is close to 1 indicates a high degree of connectivity between nodes. The statistical significance of this index was assessed with a Monte Carlo algorithm that simulated a multinomial, generative model, which took into consideration the alteration of all nodes. This statistical procedure obtained an estimation of *p*(*k*|*z*) by sampling a Dirichlet distribution and by calculating the samples' amount for which *k* > *e*, where *e* was the threshold of statistical significance set at *p* < 0.01.

#### Topological Analysis

We defined our system of interconnected nodes as a network of coatrophy areas and examined it with the network analyzer included in Cytoscape 3.5.1 (108, 109). We were therefore able to achieve a good and reliable description of the net formed by the coatrophy areas under both the aspects of brain structure and functional organization.

#### Node Degree and Edge Betweenness

The node degree was defined as the number of edges linked to a node. We employed this parameter in order to detect the nodes that were more connected within the network, which are commonly considered as brain hubs. In turn, the parameter of edge betweenness was defined as the number of the shortest routes that go through an edge in a graph or a network (110). Thus, edges exhibiting high values of betweenness are supposed to be involved in a large number of shortest routes, so that their elimination is likely to have an impact on communication between many couples of nodes.

#### Network Clustering

Given the great number of nodes as well as the high density of edges within the coatrophy network, we used the *k*-core decomposition algorithm (111, 112)–as it is implemented in the clusterMaker plugin for Cytoscape—to detect a central subnetwork of highly interconnected nodes. This algorithm eliminates all the nodes showing a degree that is lesser than a user-defined

Figure 1 | Gray matter anatomical likelihood estimation (ALE) results. The image summarizes the results of all the experiments considered in this meta-analysis. Colors from red to green show gray matter decreases [ALE maps were thresholded using voxel-level FWD *p* < 0.05 (104) and visualized using Brainvoyager QX].

*k*, thus deriving from the original network the highest connected subgraph.

# RESULTS

# Common Patterns of Morphometric Alterations

The ALE performed on all the data retrieved by our search (57 experiments, 883 subjects, and 691 foci) showed that gray matter alterations caused by AD are mainly located in the right medial frontal gyrus, the right inferior frontal gyrus, the left inferior parietal lobule, the right midcingulate gyrus, the left supramarginal gyrus, the right angular gyrus, the bilateral fusiform gyrus, the right precuneus, the bilateral insula, the right thalamus, the bilateral superior temporal gyrus, the bilateral superior temporal pole, the bilateral hippocampus, the bilateral parahippocampal gyrus, the bilateral amygdala, and the left caudate nucleus (**Figure 1**).

#### Morphometric Coatrophy Network

The left panel of **Figure 2** illustrates the 40 nodes used to build the coatrophy network, while the heat map in **Figure 2** shows the relationship between the elements of each possible couple of nodes measured by Patel's *k* index. **Figure 3** illustrates the whole coatrophy network: the colors' scale ranges from blue to red for the 146 edges and indicates an increase in *k* values. Edges are to be assumed as undirected.

Many nodes densely interconnected characterize the temporal lobe, especially the hippocampus and the parahippocampal gyrus. In contrast, only one node characterizes other brain areas, such as the cingulate cortex and precuneus. Although all the edges that are shown are statistically significant, the ones with the highest *k* value are those involving the left hippocampus, bilateral amygdala, right parahippocampal gyrus, and right inferior temporal lobe (Tables S1 and S2 in Supplementary Material).

Figure 2 | The left panels shows the nodes that entered the coatrophy calculation. The right panel shows the coatrophy matrix. Colors from blue to red indicate increasing Patel's *k* values (i.e., increasing coalteration probabilities).

Figure 3 | Morphometric coatrophy network results. Colors from blue to red indicate increasing Patel's *k* values (i.e., increasing coalteration probabilities).

**Figure 4** reports the organic option of the yFiles Layouts available in Cytoscape 3.5.1 (based on a spring-embedded algorithm) attributed to the coatrophy network. Thick links connect the nodes located in the temporal cortex, parahippocampal gyrus, amygdala, and thalamus. The right precuneus is connected to the rest of the network just through one edge projecting to the left hippocampus, whereas the right cingulate cortex is connected to the network core through the right hippocampus and the right parahippocampal gyrus. In **Figure 4**, colors and dimensions of nodes are proportional to their network degree values. In particular, Amyg\_L\_1 shows the highest degree value (17), followed by Temp\_Inf\_R (16). In turn, Fusiform\_L, Amyg\_L, Temp\_Pole\_Sup\_R, SupraMarginal\_L, and Cingulum\_Mid\_R exhibit the lowest degree value (1). The edges' thickness is proportional to their degree of edge betweenness. The edge linking the nodes Hipp\_R\_2 and ParaHipp\_R\_2 shows the highest value, while the edge between Amyg\_R and ParaHipp\_L\_1 shows the lowest one.

**Figure 5** shows the nodes according to their anatomical position. In order to simplify the visual interpretation, we have merged two or more nodes referring to the same brain area; however, we have kept the edges unchanged. It is worth noting that the coatrophy network of AD is composed of more interhemispheric (75) than intrahemispheric edges (71). Apart from the hippocampus, most of the inter-hemispheric connections link structures in the medial temporal lobes. Furthermore, unilateral nodes in the right inferior temporal gyrus and right precuneus are linked to areas of both hemispheres.

As many nodes populate the hippocampi, we projected them on a 2D template in order to better clarify their spatial localization (**Figure 6**). Five out of the six nodes in the left hippocampus were found to be located in the anterior part, while the remaining one was found to be located in the posterior section. In contrast, the right hippocampus exhibits a more uniform pattern, with two anterior nodes and one posterior.

We also analyzed the connectivity profile of the hippocampi within the coatrophy network so as to better understand their relationship with the other nodes of the network (**Figure 7**). Even though hippocampi have a lot of connections, they are scarcely interconnected (red edges) and, in particular, between the nodes of the right hippocampus there are no direct paths linking them to each other. What is more, the left hippocampus presents a greater number of edges (45) than the right hippocampus (15); however, these edges are generally characterized by a low degree of edge betweenness. In contrast, the 15 edges linking the right hippocampus to the other nodes of the coatrophy network are characterized by a high degree of edge betweenness. Overall, considering the anatomical topology of nodes (**Figure 6**), the left anterior hippocampus appears to be the most densely connected.

Figure 4 | Topological analysis of the coatrophy network of Alzheimer's disease (organic yFiles Layout). Colors and dimensions of nodes indicate their topological degree (smaller node = lower degree; from green to red = from lower to higher values). Thickness of edges indicate the degree of edge betweenness (smaller edge = lower degree).

Figure 5 | Topological analysis of the coatrophy network of Alzheimer's disease. Nodes referring to the same brain areas or strictly close one to the other have been collapsed in a single node.

Given the great number of nodes and the high density of edges of the coatrophy network, we used the *k*-core algorithm to identify the most connected components of the network. The analysis reported a core subnetwork formed by eight interhemispheric nodes (**Figure 8**), including the left and right amygdalae, left hippocampus, right parahippocampal gyrus, and right temporal inferior gyrus. The bilateral presence of nodes within this core subnetwork is consistent with the finding that the coatrophy network is characterized by a large number of interhemispheric edges.

#### DISCUSSION

With an innovative voxel-based meta-analytic method, this study aimed to find out whether gray matter decreases caused by AD distribute throughout specific and identifiable areas rather than affect randomly the whole brain. After constructing a morphometric coatrophy network, we intended to identify which brain areas are more likely to be altered in conjunction with other ones rather than alone. Finally, we examined the potential existence of relevant subcomponents within the coatrophy network.

The gray matter decreases evaluated by ALE involve limbic and temporal areas, in particular the hippocampus and parahippocampal gyrus. This finding is in accordance with most of

Figure 6 | Anatomical localization of the nodes in the hippocampi. Coordinates refers to Talairach space (right sagittal slice *x* = 25, left *x* = 30). Nodes are numerically labeled according to a rostrocaudal criterion.

previous research (30, 113). Nine out of 40 nodes of the coatrophy network are localized within the hippocampus. Specifically, six nodes are in the left hippocampus (five in its anterior part, one in its posterior part) and three in the right one (two anterior, one posterior). This is consistent with the neuropathological studies suggesting that AD is characterized by an earlier and greater involvement of anatomical structures (including hippocampus) in the left hemisphere (114–116). Although there is still debate about the exact functional organization of the hippocampus (117), the neuroscientific community has achieved a substantial consensus on its role in learning and memory (118), which are

Figure 8 | Network clustering with *k*-core decomposition algorithm. Colors and dimensions of nodes indicate their topological degree (smaller node = lower degree; from green to red = from lower to higher values). Thickness of edges indicate the degree of edge betweenness (smaller edge = lower degree). Both node degree and edge betweenness values refer to the original coatrophy network.

both deteriorated cognitive functions in AD. According to Thal et al. (119) the hippocampus (in particular the subfields CA1 and subiculum), along with the amygdala, are pretty soon affected by Aβ plaques during AD evolution (120). In line with AD diagnostic criteria (121) hippocampal and mesial temporal lobe atrophy have been considered as biomarkers of neuronal degeneration, potentially increasing the probability of an underlying AD pathophysiological process. Currently, however, the routinely utilization of hippocampal atrophy in clinical practice is not fully standardized, but preferentially applied in investigational studies and clinical trials. Furthermore, hippocampal atrophy rate could be better accounted for as a sensitive marker of disease progression (122, 123), being able to trace AD natural development and potentially representing an interesting surrogate marker for disease-modifying clinical trials (124, 125). Interestingly, an increased hippocampus and an asymmetry in the shape of the amygdala during the development of AD have been recently demonstrated, with significant correlation to cognitive impairment (126).

According to our analysis, the gray matter coatrophy network of AD appears to be densely interconnected, as suggested by the presence of 146 edges and 40 nodes, 39 of which have at least one connection. The existence of a set of nodes (altered areas) is not a proof *per se* that the disease is spreading. In fact, generally speaking, Patel's *k* is not always able to identify edges between nodes, which means that, even though some areas are altered, there is no apparent temporal coherence in their capitulation to the disease. The fact, though, that our analysis was able to discover a significant number of edges between nodes is proof of the good reliability of our results pointing out that the alteration cooccurrence really happens, as well as of the consistency of our sample.

Our analysis suggests that AD tends to target a somewhat limited set of brain regions, rather than randomly affecting distinct sites. Furthermore, the left hippocampus, bilateral amygdala, right parahippocampal gyrus, and right inferior temporal lobe seem to follow a very similar pace of degeneration (Figure S4 in Supplementary Material).

In order to evaluate the likelihood of each node of the coatrophy network to be coaltered with other ones rather than as an individual spot we calculated their node degree. The highest value pertains to the node of the left amygdala, which is reached by 17 edges, but we found other 13 nodes with at least 10 edges. These nodes are localized in the temporal lobes, right amygdala, parahippocampal gyrus, left hippocampus, and right thalamus. The high degree of pathoconnectivity of these nodes suggests that, when gray matter alteration affects one of them, it is highly probable that many other regions are also found to be altered. It is also true the other way round, that is, when nodes characterized by low degree show atrophy, it is very likely that this process cooccurs in one of the high-degree nodes, rather than in another low-degree node. These results, as well as the *k*-core decomposition, provide evidence that in the coatrophy network of AD certain nodes have the characteristic of being *pathoconnectivity hubs*. Furthermore, the values of the edge betweenness distribution indicate the existence of a dense subnetwork, which is composed of the nodes with the higher degree of pathoconnectivity.

The paucity of connections linking the two hippocampi suggests a limited cooccurrence of alterations between them. The hippocampus is known to be greatly affected by AD, and the MRI volume estimation of this structure is currently considered one of the most reliable *in vivo* biomarker of this disease (62). Our results suggest that both the hippocampi are substantially altered, albeit somewhat independently. According to previous studies, certain molecular alterations typical of AD are more evident in the left hippocampus compared to the right one (127, 128). This discovery might explain the abundance of edges connecting the nodes in the left hippocampus, as well as support the transneuronal spread mechanism in AD. The nodal stress hypothesis could also play a role in virtue of the intense functional activity of this region. Finally, our finding that the anterior part of the hippocampus exhibits a greater number of edges than the posterior part seems consistent with the suggestion that the deterioration of CA1 and subiculum appears to be more correlated with the development of AD than the deterioration of CA3, which appears to be more correlated with healthy aging (11, 120). Recently, the presubicular–subicular complex has been described as one of the earliest site of atrophy in AD, with a significant correlation with memory performances (even in MCI phase), potentially reflecting the ongoing degenerative process through the subiculum passing from entorhinal cortex to dentate gyrus (129, 130).

In addition to the interpretation of the coatrophy network as a whole, some specific aspects deserve a detailed consideration. The first is the relationship between hippocampus and precuneus. In the coatrophy network of AD these regions are linked through an edge exhibiting a very high degree of edge betweenness, which reveals a direct interaction. According to the "hippocampus disconnection hypothesis" proposed by Tahmasian et al. (131), the disruption of functional connectivity between hippocampus and precuneus could induce the characteristic alterations in the hippocampus that we find in AD. Tahmasian et al. (131) have in fact demonstrated that in AD the hippocampus is much less inhibited, and this disinhibition may result in its hypermetabolism. A similar situation could induce neurotoxicity, which might be one of the causes behind gray matter decrease measured with VBM, thus explaining the identification of a significant number of nodes in the hippocampus.

A second interesting aspect is the relationship between the left hippocampus and right inferior temporal gyrus, which was highlighted by *k*-core decomposition. This result is in agreement with the study of Wang et al. (132), which found that the interaction between these two areas is typical of AD. Of note, Wang et al. (132) examined 80 pathological subjects using Bayesian network analysis and prior-defined regions of interest, while the present study has applied a meta-analytical approach on a substantially bigger VBM database of 883 patients diagnosed with AD. This agreement supports the sensitivity of our novel methodology. Furthermore, the slight prevalence of inter-hemispheric connections in the coatrophy network of AD (see Figure S5 in Supplementary Material) is consistent with the deterioration of white matter bundles in AD, in particular concerning the corpus callosum (133–137). Callosal atrophy has been associated with cognitive decline rate as well as to disease progression (138, 139).

Gray matter alterations found in the hippocampus, precuneus, and inferior parietal cortex can be ascribed to the general disruption of the DMN in patients with AD (58, 140). Recently, a study has showed that the DMN dysfunction, as well as the disruption of the interaction between different resting state functional networks, can be attributed to amyloid burden (58). What is more, Chang et al. (141) have found that amyloid burden in the cingulate cortex might promote gray matter atrophy in the other areas constituting the DMN.

Overall, the crucial role played by pathological proteins in AD supports the transneuronal spread hypothesis at the basis of gray matter alterations' distribution (4, 5, 39, 40, 42, 45). However, the complex relationship among different factors (such as amyloid burden, Tau deposition, gray matter atrophy, and disrupted functional connectivity) and the presence of several hub nodes within the coatrophy network of AD suggest that the nodal stress mechanism could as well be involved in the development of the disease (142). Therefore, it is extremely likely that different spreading mechanisms, which are not mutually exclusive, may be involved in the etiology of AD.

#### Limitations and Future Directions

The present investigation and the methodology on which it is based aim to better understand the nature of AD by examining its pathological fingerprints over the brain. To do so, we were able to get access to a very large sample size of patients. If this is an advantage on the one hand, it can also be a limitation on the other, as within this sample it was not possible to determine the average duration of disease, due to unavailability of information in the original studies. This aspect makes it difficult to associate the coatrophy network with a specific stage of AD progression. However, the methodological procedure for defining the areas to be included in the coatrophy network considers primarily the frequency of every single area to be found altered. In case of a neurodegenerative condition such as Alzheimer's we could imagine, generally speaking, a group of patients with a recent diagnosis exhibiting alterations in area A, another group with an intermediate development of the disease exhibiting alterations in areas A–B, and another group with an advanced development of the disease exhibiting alterations in areas A–B–C. Since our methodology privileges the frequency of each area to be found altered, in the final network area A will be more likely to be represented, while area C may be even excluded. Moreover, even if the group of patients exhibiting alterations in A–B–C were greater than the other groups, the pattern A–B–C would be less likely to be represented than the sole area A. For this reason, even if our input data could contain an overrepresented sample of patients in a specific stage of the disease, the resulting coalteration network would not represent the pattern of altered areas which is typical of that stage.

Future studies on longitudinal data analyzed by different methods are needed in order to investigate the sequential formation of the coatrophy network identified in this study, so as to achieve a more detailed picture of the temporal evolution of AD.

# CONCLUSION

This meta-analysis was able to address the following important issues.


The innovative methodological analysis developed in this study for constructing the morphometric coatrophy network of an important neurodegenerative disease such as AD opens a new window into the comprehension of the pathological brain. Increasing evidence is supporting the idea that brain alterations distribute according to a network-like structure. The analysis carried out in this study not only provides support for this hypothesis but also puts forward the significant finding that certain nodes of the coatrophy network may play the role of pathoconnectivity hubs. What is more, our methodology can be equally applicable to study the morphometric coalteration network of any other neuropathological condition. Future investigations into this line of research on databases of different diseases promise to provide valuable insight to the study of the dynamics of brain disorders, so as to achieve a better predictive diagnostic power as well as to improve medical care and treatment.

# AUTHOR CONTRIBUTIONS

JM and AN implemented data collection, analyzed the data, drafted, and revised the article. EP, BB, and KT drafted and revised the article. TC designed the analysis tool, supervised data analysis, drafted, and revised the article. DL retrieved information of the sampled population and implemented bibliographic research. SD revised the article. FC conceived the experiment, supervised data collection, supervised data analysis, drafted, and revised the article.

# FUNDING

This study was supported by the Fondazione Carlo Molo (FC, PI), Turin and by CSP UNITO Excellent Young PI grant CSTO162182 (FC, PI), Turin.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fneur.2017.00739/ full#supplementary-material.

#### REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer AK and handling editor declared their shared affiliation.

*Copyright © 2018 Manuello, Nani, Premi, Borroni, Costa, Tatu, Liloia, Duca and Cauda. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Predictive Model of Spread of Progressive Supranuclear Palsy Using Directional Network Diffusion

*Sneha Pandya1 , Chris Mezias2 and Ashish Raj1,2\**

*1Department of Radiology, Weill Cornell Medicine, New York, NY, United States, 2Department of Neuroscience, Weill Cornell Medicine, New York, NY, United States*

Several neurodegenerative disorders including Alzheimer's disease (AD), frontotemporal dementia (FTD), Parkinson's disease (PD), amyotrophic lateral sclerosis, and Huntington's disease report aggregation and transmission of pathogenic proteins between cells. The topography of these diseases in the human brain also, therefore, displays a wellcharacterized and stereotyped regional pattern, and a stereotyped progression over time. This is most commonly true for AD and other dementias characterized by hallmark misfolded tau or alpha-synuclein pathology. Both tau and synuclein appear to propagate within brain circuits using a shared mechanism. The most canonical synucleopathy is PD; however, much less studied is a rare disorder called progressive supranuclear palsy (PSP). The hallmark pathology and atrophy in PSP are, therefore, also highly stereotyped: initially appearing in the striatum, followed by its neighbors and connected cortical areas. In this study, we explore two mechanistic aspects hitherto unknown about the canonical network diffusion model (NDM) of spread: (a) whether the NDM can apply to other common degenerative pathologies, specifically PSP, and (b) whether the directionality of spread is important in explaining empirical data. Our results on PSP reveal two important findings: first, that PSP is amenable to the connectome-based ND modeling in the same way as previously applied to AD and FTD and, second, that the NDM fit with empirical data are significantly enhanced by using the directional rather than the non-directional form of the human connectome. Specifically, we show that both the anterograde model of transmission (some to axonal terminal) and retrograde mode explain PSP topography more accurately than non-directional transmission. Collectively, these data show that the intrinsic architecture of the structural network mediates disease spread in PSP, most likely *via* a process of trans-neuronal transmission. These intriguing results have several ramifications for future studies.

Keywords: tauopathies, progressive supranuclear palspy, network diffusion, directionality, 4R

#### INTRODUCTION

Several neurodegenerative disorders, including Alzheimer's disease (AD), frontotemporal dementia (FTD), Parkinson's disease (PD), amyotrophic lateral sclerosis (ALS), and Huntington's disease (HD), report aggregation and transmission of pathogenic proteins between cells (1–6). The topography of these diseases in the human brain also, therefore, displays a well-characterized and stereotyped

#### *Edited by:*

*Justin John Yerbury, University of Wollongong, Australia*

#### *Reviewed by:*

*Melissa Calegaro Nassif, Universidad Mayor, Chile Wenquan Zou, Case Western Reserve University, United States*

> *\*Correspondence: Ashish Raj ashish@med.cornell.edu*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 16 August 2017 Accepted: 04 December 2017 Published: 21 December 2017*

#### *Citation:*

*Pandya S, Mezias C and Raj A (2017) Predictive Model of Spread of Progressive Supranuclear Palsy Using Directional Network Diffusion. Front. Neurol. 8:692. doi: 10.3389/fneur.2017.00692*

regional pattern and a stereotyped progression over time. This is most commonly true of AD, as well as other dementias characterized by hallmark misfolded tau or alpha-synuclein pathology, both of which appear to propagate within brain circuits using a shared mechanism. The most canonical tauopathy is AD, but a diverse group of related taupathies are known: FTD, corticobasal degeneration, semantic dementia, and many more (7). Far less studied is a rare disorder called progressive supranuclear palsy (PSP). The hallmark pathology and atrophy in PSP are similar to those in other tauopathies, and its regional patterning is likewise highly stereotyped: initially appearing in the striatum, followed by its neighbors and connected cortical areas (8, 9).

The mechanisms underlying stereotyped patterning and progression in tauopathies are not fully understood, and both cell–cell communications governed by anatomical and functional connections, and cell autonomous molecular factors characterized by gene expression signatures, could conceivably play a role in vulnerability to disease spread (10–15). Foremost among non-cell autonomous factors is network connectivity, which is increasingly considered a plausible and even key driver of vulnerability (16–18).

Assuming that trans-neuronal transmission must proceed along axonal projections, the spatiotemporal dynamics of pathology spread can be given quantitatively and deterministically from the inter-regional anatomic connectivity patterns (19). Modern diffusion-weighted magnetic resonance imaging (dMRI) (20) and post-processing techniques like fiber tractography (21) and connectivity mapping (22) have enabled the computation of interregional connectivity, yielding what is frequently called the human "connectome." Recently, our group proposed a mathematical model of prion-like trans-neuronal spread of neurodegenerative pathology called the network diffusion model (NDM), evolved on dMRI-based structural networks or connectomes. This model demonstrated that observed spatial patterns of neurodegeneration in common degenerative diseases like AD and FTD can be explained simply as a consequence of network spread (16). This model also gives an explanation of selective regional vulnerability in terms of disease epicenters, called eigenmodes, associated with pathology. Since the model is deterministic, it was successfully employed to predict future atrophy patterns of AD subjects from their baseline patterns and connectomes (17).

In this study, we explore two mechanistic aspects hitherto unknown about the canonical NDM of spread: (a) whether the NDM can apply to other common degenerative pathologies, specifically PSP, and (b) whether the directionality of spread is important in explaining empirical data. We, therefore, apply the NDM to quantitatively assess whether trans-neuronal transmission of PSP pathology can recapitulate observed PSP topography. Our empirical data come from an unprecedented dataset of 60 PSP subjects from the 4-Repeat Tauopathy Neuroimaging Initiative (4RTNI) study, a multinational prospective observation study that examines clinical, radiologic, and biological findings of disease progression in tauopathic individuals, including PSP.1 Our interest in PSP arises from its distinct spatial pattern in comparison to AD; hence the ability of NDM to explain PSP pattern would contribute to the emerging notion that all neurodegenerative pathologies follow shared mechanisms of spread.

To assess the second question, we propose a novel construction of a *directional human connectome*, for the first time. Clearly, directionality of tracts or inter-regional anatomic connectivity, as defined by the polarity of individual axonal projects (soma to axonal terminal or *vice versa*) is impossible to determine from dMRI data, as water diffusion along fiber bundles does not respect cell polarity. Instead, we exploit the well-known fact that homologous structures exist between many species, for some of whom we do happen to have anatomic connectivity data from painstaking tracer studies. Using retrograde tracer studies, a detailed mesoscale mouse connectome was reported by the Allen Institute (23). We, therefore, developed a novel technique whereby human and mouse brain atlas parcellations are used to define homologous brain structures, and the Allen mouse directional connectivity is transferred to non-directional human connectome. The importance of studying directional connectomes is that *in vitro* and *in vivo* mouse studies are increasingly revealing directional preference in the transmission of various misfolded proteins; however, conclusive data on directionality of each protein are not currently established (24).

Our results on PSP reveal two important findings. First, that PSP is amenable to the connectome-based ND modeling in the same way as previously applied to AD and FTD. This establishes that PSP might propagate using the anatomic connectome in the same way that is known for AD and FTD. Second, we found that the NDM fit with empirical data are enhanced by using the directional rather than the non-directional form of the human connectome. Specifically, we show that the anterograde model of transmission (some to axonal terminal) explains certain aspects of PSP topography more accurately than the non-directional model. Similarly, certain aspects of PSP are better explained by the retrograde mode of transmission. Overall, both directional models outperform the non-directional model, and retrograde mode is overall the most accurate. These intriguing results have several ramifications for future studies.

#### MATERIALS AND METHODS

#### Building an Anatomic Connectome from Parcellated Atlas and dMRI of Healthy Cohort

Axial T1-weighted structural fast spoiled gradient-echo scans (TE = 1.5 ms, TR = 6.3 ms, TI = 400 ms, 15° flip angle, 230 × 230 × 156 isotropic 1 mm voxels) and high angular resolution diffusion imaging data (55 directions, *b* = 1,000 s/mm2 , 72 1.8-mm-thick interleaved slices, 0.8594 mm × 0.8594 mm planar resolution) were acquired on a 3T GE Signa EXCITE scanner from 73 fully consented young healthy volunteers under a previous study approved by Weill Cornell's institutional review board; for details of study protocols see Ref. (25). The exclusion criteria were pregnancy, a history of neurological or psychiatric diagnosis, seizure, or drug or alcohol abuse. Demographic characteristics of these young healthy subjects are

<sup>1</sup>http://4rtni-ftldni.ini.usc.edu.

Table 1 | Demographic and clinical details of young healthy controls used to construct canonical structural connectome, PSP, and age-matched healthy cohort.


*MoCA, Montreal Cognitive Assesement; UPDRS, Uniform Parkinson's disease Rating Score.*

shown in **Table 1**. Tractograms were extracted from these 73 young healthy subjects to create the normative connectome for the study. The measure of connectivity used in this paper is the anatomical connection strength (ACS) as proposed in Ref. (26). The overall differences in regional connectivity were normalized by a scaling factor equal to the sum of the connections. Thus, the connection strength, *ci*,*<sup>j</sup>*, between ROIs *i* and *j*, was defined as the ACS score of streamlines connecting the two regions *i* and *j*. We refer to this network as graph *G* = {*V*,*E*} whose nodes ν*<sup>I</sup>* ∈ *V* represent the ith GM region, and edges *ei,j* ∈ *E* represent white matter fiber pathways whose connection strength is *ci*,*<sup>j</sup>*. Connections are assumed to be bidirectional since directionality is not deducible from diffusion tensor imaging (DTI) tractography data. Gray matter regions from the T1-weighted images were parcellated using FreeSurfer volumetric pipeline (27) and a Desikan-Killiany atlas (28) into 68 cortical regions, 34 from each hemisphere and 18 subcortical regions. Six subjects were eliminated from the 73-subject dataset due to FreeSurfer failure or insufficient MR contrast.

## Obtaining Directional Human Connectomes by Utilizing Directional Information from Mesoscale Mouse Connectomic Data

#### The Allen Brain Institute (ABI) Mouse Connectivity Atlas

Connectivity data were taken from the supplementary dataset published along with the mesoscale mouse connectome from the ABI (23). The ABI generated their mouse connectivity data using an anterograde viral tracer engineered to express GFP natively using a promoter with high affinity for the ubiquitous transcription factor synapsin (23). The specific vector they used, rAAV1, travels exclusively along axons, across synapses, and through dendrites in the CNS, and almost exclusively travels in an anterograde direction (over 99% of the time). Viral vectors were injected and then mice were sacked after a period of approximately 3 weeks and their brains assessed for the expansion of GFP expressing viruses into other regions and for expansion of the GFP expressing virus within the targeted region. Brains were then sliced and each slice was analyzed by an automated camera system to sum the number of pixels across all voxels in each cytoarchitectonically defined brain region; these values were then normalized by injection volume, to give a per voxel connectivity density measurement between any two areas out of the total 213 regions defined in the data. The cohort of mice used to generate the dataset was of *N* = 1,231 C57/BL6 males, and regions were taken from the previously created ABI Mouse Brain Ontological Reference Atlas, created using specific staining for regional cell types with a cohort of *N* = 469 C57/BL6 male mice. C57BL/6 male mice are a common inbred strain of laboratory mouse and the most widely used mice for use as models of human disease and neurodegenerative disorders. There are several differences in connectivity between mouse lines. No comparable mesoscale whole brain connectome has ever been created for another mouse strain, hence at this point, any large-scale wiring differences in the connectomes of different mouse strains are unknown. Further methodological details and descriptions of the assumptions used to create the dataset can be found in the supplementary information in Ref. (23).

Total projection density between regions was generated by multiplying element-wise by the rows the connectivity matrix times the number of voxels in each seeding region. Next, the 426-region mouse connectome was used to inform the directionality of specific connections in the 86-region human connectome C. Custom MATLAB code then transformed the ABI connectivity matrix CABI of size 426 × 426, whose nodes reflected a parcellated atlas consisting of 426 GM regions, into a custom connectivity matrix of size 86 × 86, based on the human atlas parcelation. To achieve this, the afferent and effect connections of all samples that fell within the (i,j)-th pair of human homolog regions were summed to give the in- and out-degrees of the human regions: degreein and degreeout. The overall directionality was estimated by the anterograde ratio *D*ant in in out degree degree + degree <sup>=</sup> . This directionality ratio was then applied to the non-directional human connectome

*C* to obtain the (anterograde) directed connectome *C*ant = *D*ant.*C*. Most of the neocortical connections are bidirectional, and only a small percentage shows strong directionality (see Results). These connections emanate from primary sensorimotor cortices and subcortical structures, as expected. In order to ensure that mouse directionality was relevant to human data, each strongly directional connection was qualitatively verified from macaque data available in the literature, especially from the CoCoMac database (29), which provides an ordinal measure of directionality (one of three levels of connectivity in either direction).

We applied directional connectivity to all human regions for which there were clear mouse analogs: primary sensory neocortical regions, and subcortical regions. The mouse atlas was transferred to a human atlas by taking directional connectivity ratios (outgoing/incoming) per region-to-region connection for all human region pairings both having mouse analogs. We then multiplied the directed connectivity ratios from mouse analog region pairings with the region-to-region connectivity strength in the Desikan Atlas. Here, we define mouse analogs to human regions in two ways: first, we used region pairings where both human regions had directly comparable mouse regions. For example, the precentral gyrus—postcentral gyrus connection is anatomically analogous to primary motor area—primary somatosensory area connection in mice. Second, we combined mouse regions from the Allen Institute atlas into the larger human regions from the Desikan Atlas when the larger combined region in mice was directly analogous to a human region. For instance, the Allen Institute atlas has the hippocampus broken into the various CA regions and the subiculum; the connectivity profiles of these were summed into one hippocampal region, as per the Desikan Atlas, before taking directional connectivity ratios of connections with the hippocampus. All regions not fitting one of the two above mentioned criteria were ignored and no directionality was imputed on human region-to-region connections where both regions did not have a mouse analog. Prior work has noted similarities between mouse and functional connectomes and the regions of the default mode network (30), suggesting that for a first pass at applying directed connectivity to a human connectome, mouse directional connectivity ratios can suffice. However, all results should be interpreted through the lens of the clear limitations of imputing gross wiring patterns from one species to another.

#### Age-Matched Normal and PSP Cohorts

Data used in preparation of this article were obtained from the 4RTNI database (see text footnote 1). Neurological and physical examinations, cognitive testing, behavioral testing, and other clinical data from PSP subjects between the ages of 45 and 90 were collected over three longitudinal visits (baseline, 6 months, and 12 months). The conditions for exclusion were any significant neurological disease other than PSP, including PD, multi-infarct dementia, HD, brain tumor, multiple sclerosis, or history of significant head trauma followed by persistent neurological deficits or known structural brain abnormalities. Demographic characteristics of these subjects are shown in **Table 1**. For our study, we used 60 PSP subjects after elimination due to unknown age, gender, and failed FreeSurfer quality control. Age-matched healthy controls (NC) data used in this study were obtained from the public Alzheimer's disease neuroimaging initiative (ADNI) database2 consisting of multimodal imaging data on AD and healthy subjects. T1-MPRAGE, FLAIR, DTI, ASL-perfusion (ASL), and resting state functional MRI (rs-fMRI) MRI sessions were acquired at 3T. The parameters for the PSP scans were chosen to match those used by the ADNI, including those being developed for ADNI-2 for ASL, rs-fMRI, and DTI at the time. Two-sample *t*-test was used to generate a random sample of 150 age-matched controls from ADNI database. Demographic characteristics of these subjects are shown in **Table 1**. This randomly generated sample of 150 age-matched healthy controls was used to generate a group atrophy vector for the study.

68 cortical and 18 subcortical volumes from 3T T1-weighted baseline MRI images were extracted using FreeSurfer software for both the cohorts. Estimated Total Intracranial Volume generated by FreeSurfer was used as an estimate for intracranial Volume (ICV) as a normalization measure to correct for head size in this study. ICV has previously been used in several studies for normalization particularly for neurodegenerative disorders for better estimation of regional atrophy that is caused by pathology (31, 32). Image processing steps were visually inspected for white–gray matter boundary and skull-stripping errors to ensure they had been carried out correctly. Subjects were eliminated from the dataset due to FreeSurfer failure or insufficient MR contrast. A vector of regional atrophy was created by using a two tailed *t*-test between PSP and NC mean ICV corrected regional volumes such that *t*PSP = {*t*PSP(*i*)|*i* ∈ [1,*N*]} (*N* = 86). The *t*-statistic was converted to the natural range [0,1] using the logistic transform, following (17). These atrophy measures were then used to test the propagation modeling analyses.

#### A NDM for Taupathy Spread

We modeled taupathy progression as a diffusion process on graph *G*. From Ref. (16) the transmission of pathology to all brain regions *via* the whole brain disease vector *x*(*t*) = {*x*(ν,*t*), ν ∈ *V*} and:

$$\frac{d\mathbf{x}(t)}{dt} = -\beta H\mathbf{x}(t),\tag{1}$$

*w*here β is a global diffusivity constant and *H* is the well-known graph Laplacian.

*H I* = − *D C*<sup>−</sup><sup>1</sup> ,

where *D* is a diagonal matrix whose diagonal entries contain the degree of each node, degree being defined as the sum of weighted connections emanating from the node. Note, in order to accommodate regions having widely different out-degrees, we have used above the degree-normalized version of the Laplacian matrix (17).

#### Directional Laplacian

Next, we define a directional (anterograde) connectome *C*ant. Since this matrix is non-symmetric, we define for each node an in-degree and an out-degree given by row and column sums of the matrix: **d d** row ant ant , , , , , *<sup>i</sup> j i j j i C Ci j* = = ∑ ∑ col . Define *H*ant as the anterograde graph Laplacian matrix

$$H^{\rm unt} = I - diag(\mathbf{d}\_{\rm row}, \mathbf{d}\_{\rm col})^{-\frac{1}{2}} C^{\rm unt}$$

Equation 1 is applied to both the non-directional and anterograde matrices and in each case admits a closed-form solution *x*(*t*) = *e*−β*Htx*0, where *x*0 is the initial pattern of the disease process, and we call term *e*−β*Ht* the *diffusion kernel* since it acts essentially as a spatial and temporal blurring operator on *x*0. The unit of model's diffusion time *t* is arbitrary. Global diffusivity β is unknown; hence, we chose a value that would roughly span tau progression (10–20 years), giving β = 0.15. In future, both *t* and β would be estimated by fitting to longitudinal data, and would then acquire correct dimensions and units.

Subsequent diffusion of tau pathology out of potential seed at the *k*-th region is given at any time point *t* by

$$\mathfrak{x}\_{\rm psp}(t) = e^{-\beta Ht} \mathfrak{e}\_{\hbar},\tag{2}$$

where *ek* is a unit vector with 1 at the *k*-th entry and 0 elsewhere. Note that although the above model involves pathology, what we have available to us empirically regional MRI-derived atrophy. Hence, an underlying assumption in all analyses henceforth is

Frontiers in Neurology | www.frontiersin.org December 2017 | Volume 8 | Article 692

that the empirical atrophy vector is proportional to the pathology vector, hence both are given by *x*PSP(*t*).

#### Statistical Analysis

Two cerebellar regions were removed, leaving regional PSP atrophy statistics on 84 cerebral regions. The NDM is described by *x*(*t*) and Φ(*t*), two 84 × 100 tables that represent the pathology (unobserved) and atrophy (measured empirically) in all 84 ROIs over 100 points in time. Pearson correlation strength (*R* statistic) and associated *p*-values were calculated comparing the empirical atrophy vector, *t*PSP, with the ND atrophy pattern at all 100 points in time.

#### Repeated Seeding

Next, the NDM was run 84 times, each time starting from a different ROI, in order to deduce the most likely seed regions. For each node *i*, we start the model with *x*<sup>0</sup> = *ei*, where *ei* is a unit vector with 1 at the *i*-th location and 0 elsewhere. In the current case, we chose to seed bilaterally, so that two entries in the "unit" vector were assigned 1. This was repeated for each region in turn, and the NDM-predicted atrophy pattern Φ*<sup>i</sup>* (*t*) was calculated from Eq. 3. This gave 84/2 = 42 different Φ (*t*) matrices. For each predictor matrix, the corresponding Pearson's *R* was calculated at all model times *t*, giving *Ri* (*t*). These *Ri* (*t*) values were plotted on common axis, giving what we denote as "*R*–*t* curves." From each *Ri* (*t*), we recorded the maximum value *R<sup>i</sup>* max , which is used here as an indicator of model evidence reflecting the likelihood of the region i being the true region of pathology onset. This effectively "ran the NDM backwards," allowing us to determine which seeded ROI would serve as the most likely origination site to subsequently yield the regional patterns closest to empirical data.

#### Potential Seeding

Given that hypothalamus (HTH) seeding produced the best *R* against empirical data (see Results), we explored seeding from these regions, given by the initial vector *x*<sup>0</sup> = *eHTH*, which yields the model predictor ΦPSP(*t*) = exp(−β*Ht*)*eHTH*. Snapshots of the evolving ΦPSP(*t*) vector were recorded and plotted in glassbrain renderings at selected model times *t* = 7, 15, 22 years with non-directional connectivity, *t*= 13, 26, 39 years with anterograde directional connectome, and *t* = 14, 27, 41 years with retrograde directional connectome.

#### Random Scrambling

In order to build a null distribution for assigning significance to the NDM, we performed two levels of randomization experiments. (1) We ran the NDM on 2,000 randomly scrambled versions of the connectivity matrix C. C was scrambled using a symmetric transformation of the network's nodes by randomly permuting entire rows and columns, and the NDM was evaluated for each shuffled network after bilateral HTH seeding for PSP. This scrambling procedure maintains the edge and node statistics of the true connectivity C. The NDM evaluated on these 2,000 randomly scrambled networks, therefore, constitute null or reference models which supplied significance values to results of the true model. (2) We ran the NDM on 2,000 randomly scrambled PSP atrophy vector. Atrophy values in *t*PSP vector were randomly assigned among the 84 cerebral regions with 2,000 different permutations. This scrambling method maintained the true connectivity C, but replaced true regional atrophy pattern with a random distribution of atrophy.

# RESULTS

## Cross-sectional Spatial Distribution of PSP Atrophy

**Figure 1** illustrates "glass brain" renderings of the spatial distribution of PSP atrophy. The spheres are located at the centroid of each of 84 brain regions, their size is proportional to the *t*-statistic of PSP atrophy after logistic transform and color coded by lobe per: frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green. The relative order of regional atrophy in this cohort roughly mirrors the spatial pattern of atrophy. HTH, pallidum, entorhinal, inferiortemporal, and superior frontal are the most atrophied regions as seen by the largest spheres in **Figure 1**.

# Characterizing the Directional Connectome

The properties of the directional connectome built using both human and mouse connectome are described in **Figure 2**. First, we show that the connectivity patterns (**Figures 2A,B**) as well as the statistics of both the directional and the original nondirectional connectomes are very similar (**Figures 2C,D**), as they should be. The directionality, defined such that −1 represents a purely retrograde connection, and +1 represents a purely anterograde connection, is shown as a histogram in **Figure 2E**. Clearly, most connections are bidirectional (value 0) and a small minority is directional, with equal numbers in both directions. Mainly subcortical structures display strong directionality (bottom and right portions of the matrices displayed in **Figures 2A,B**) while most corticocortical connections are bidirectional.

#### Repeated Seeding of the NDM with Non-Directional Connectivity

Next each region was computationally "seeded" in turn and NDM was played out over time on the canonical healthy connectome. **Figure 3A** shows spread of maximum Pearson correlation strength *Ri* max corresponding to the best fit between empirical data and the NDM seeded at region *i*. The distribution of *R<sup>i</sup>* max is being shown as a histogram. Since *R<sup>i</sup>* max serves as a measure of the likelihood of each region being a seed, this information is spatially depicted in the "glass-brain" insets. **Table 2** shows top 20 regions with maximum correlation strength for each region seeded in turn. Clearly, the HTH and other limbic structures serve as the best seed regions for PSP. **Figure 3B** shows the *R*–*t* curves between the model evolution (Eq. 2), seeded at each region in turn and empirical PSP atrophy. For seed regions that are plausible, this would yield an intermediate peak in *R* where the NDM best matches empirical data, then diffusing uniformly and resembling actual data less and less. The highest *R* was achieved by the HTH, a region known to suffer early atrophy and could most likely be a seed to PSP for the connectome-based network diffusion.

Next in a model-free analysis, we establish the role played by proximate anatomic features in governing the regional patterns of PSP atrophy. We considered HTH as anatomic predictors based on highest *R* achieved from repeated seeding analysis above. Each covariate is a 84-long vector, covering the entire brain. Since the group atrophy *t*-statistic is generally bilateral, we removed lateralization effects by averaging the left- and right-hemispheric values of these vectors, giving predictor vectors of size 42 × 1. Linear bivariate correlation analyses of regional PSP atrophy with these proximate predictors are shown in **Figure 3C**. Dots are color coded as per lobes (frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green). We see a non-significant negative correlation between PSP atrophy and connectivity from bilateral HTH (*R* = −0.11, *p* < 0.01).

#### Network Diffusion Out of Potential Seed Recapitulates Regional Atrophy in PSP Tauopathy Using a Non-directional Connectome

Having demonstrated that bilateral HTH gives the best seeding, we next captured the spatiotemporal evolution of *x*HTH(*t*) (**Figures 4A,B**). The maximum of *R*HTH(*t*) occurs at *t*max= 22, hence *x***HTH**(*t*max) is shown in **Figure 4A**. The ND progression matches quite closely the stereotypical sequencing in PSP, where the disease in due course extends into the subcortical and temporal structures and finally progresses to increasingly involve the cerebral cortex. The maximum correspondence of NDM to empirical data, given by the peak of the *R*–*t* curve in **Figure 3B** occurs at *t* = 22 for PSP. **Figure 4B** shows the evolution of network diffusion process seeded at the HTH starting at early stage (*t* = 7) through mature stage (*t* = 22), the model increasingly resembling empirical atrophy of **Figure 1**. We have selected three equidistant years at *t* = 7, 15, and 22 years between early stage (*t* = 0) with no diffusion through mature stage (*t*= 22). Here, time is arbitrary, and the use of "years" is meant for illustrative purpose. Initial spread is followed by diffusion especially into the pallidum and caudoputamen, followed by other striatal and limbic structures, then to middle temporalfrontal structures, and finally to the cortex.

#### Testing for Significance against Alternate Randomized Models

We evaluated the NDM against alternate network models to show its specificity to PSP atrophy and to the connectome on which it evolves. We evaluated this in two ways and recorded the best *R* achieved by each model. First, we randomly scrambled the healthy average connectivity matrix C 2,000 times, and ran the NDM on each scrambled network for PSP. Second, we randomly scrambled the group *t*-statistic of regional PSP atrophy vector and ran the NDM on original connectivity matrix C. The distribution of Pearson's *R* over 2,000 scrambled matrices is shown in **Figure 5A** and clearly indicates that our correlation results are unlikely to be due to chance. There is a hard limit on the left of this plot at *R* ~ 0.42, which corresponds to the 0-diffusion time value of *R*PSP(*t*) curve in **Figure 3B**. The second random scrambling experiment, where the atrophy vector was scrambled instead of the connectome, gave an *R* distribution shown in **Figure 5B**. This distribution was approximately Gaussian, with mean between 0.1 and 0.2, and SD around 0.1 and 0.2. Random model's *R* was much lower than the maximum *R* of 0.42 achieved by the true model; statistically outside the 95% confidence interval, or *p* < 0.05. Hence, the reported HTH seeded NDM prediction cannot be explained by chance.

## Repeated Seeding of the NDM with Directional Connectivity

Having demonstrated that bilateral HTH gives the best seeding, we next wanted to show if directionality could predict tauopathic

and the non-directional connectomes are also very similar. As expected the values of connections are not different between the non-directional and the directional connectome. (E) Histogram of the directionality defined by (*Ci*,*<sup>j</sup>* − *Cj*,*i*)/(*Ci*,*<sup>j</sup>* + *Cj*,*i*), such that −1 represents a purely retrograde connection, and +1 represents a purely anterograde connection. Most connections are bidirectional (value 0) in E and a small minority is directional, with equal numbers in both directions. Mainly subcortical structures display strong directionality [bottom and right portions of the matrices displayed in (A,B) while most corticocortical connections are bidirectional].

spread either from same or different set of potential seeds. We tested both retrograde and anterograde mode in our model at *t* = 41 and *t* = 40, respectively. With anterograde connectivity there was an increase in maximum *R* (*R* = 0.48 at *t* = 40) as seen in *R*–*t* curve (**Figure 6A**). Our results showed that retrograde mode showed significant improvement in maximum *R* (*R* = 0.53 at *t* = 41) to both bidirectional and anterograde mode (**Figure 6C**). Both with retrograde and anterograde mode, the highest *R* was achieved by the HTH and served as the best seed region. **Figures 6B,D** show glass-brains with maximum *R* achieved by each region seeded in turn and **Table 2** shows top 20 regions with maximum correlation strength for each region seeded in turn with directional connectivity.

#### Directional Network Diffusion Recapitulates Regional Atrophy in PSP

We next captured the spatiotemporal evolution of *x*HTH(*t*) (**Figures 7A–D**) at *t* = 40 and *t* = 41, which corresponds to time taken by diffusion to reach the highest peak with anterograde and retrograde directional connectivity, respectively. The maximum of *R*HTH(*t*) occurs at *t*max = 40 and 41 for anterograde and retrograde mode, hence *x***HTH**(*t*max) are shown in **Figures 7A,C** for each mode, respectively. The ND progression shows better prediction of tau progression in PSP compared to non-directional seeding as seen in **Figure 4**. Both with anterograde and retrograde mode, the spreading pattern extends into the subcortical and temporal structures and finally progresses to increasingly involve the cerebral cortex.

The maximum correspondence of NDM to empirical data, given by the peak of the *R*–*t* curve in **Figure 6A** occurs at *R* = 0.48, *t* = 40 and in **Figure 6C** occurs at *R* = 0.53, *t* = 41 for PSP. **Figures 7B,D** shows the evolution of network diffusion process seeded at the HTH starting at early stage (*t* = 13) through mature stage (*t* = 39) with anterograde mode and through *t* = 14 and *t*= 41 with retrograde mode. With both the modes, the model

Figure 3 | Results for repeated seeding analysis. Each region was seeded in turn and network diffusion model (NDM) was played out for all time points. Pearson's *R* was recorded at each time point between the model and PSP atrophy vector. As the diffusion time increases, more and more of the pathogenic agent escapes the seed region and enters the rest of the network. (A) The point of maximum correlation with measured atrophy using structural connectome was recorded with glass brains of measured *R* shown inset in PSP. (B) NDM seeded at bilateral regions using structural connectome indicates that the hypothalamus (HTH) (shown by red curve) could be a plausible candidate for PSP disease seeding—it has the highest peak *R*, and the characteristic intermediate peak indicative of true pathology spread. The unit of model's diffusion time *t* is arbitrary. For the purpose of demonstration, we consider this arbitrary timescale parameter in terms of "years" but without a robust fitting approach with longitudinal data these time values should be considered only roughly equivalent to years. (C) Once best seed was determined we looked at correlations of mean connectivity from HTH versus empirical PSP atrophy. We see a non-significant correlation between PSP atrophy and connectivity from bilateral HTH (*R* = −0.11, *p* < 0.01). Dots are color coded by lobe—frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green.


Table 2 | Top 20 regions with maximum correlation strength for bilaterally seeded ROIs without and with directional connectivity.

Figure 4 | (A) ND prediction from bilaterally seeded hypothalamus (HTH) with non-directional connectome. Glass brains of network diffusion model seeded at the bilateral HTH at *t* = 22 years shows spatial pattern from HTH using bidirectional connectome. Subcortical striatal structures especially caudoputamen and globus pallidus, and limbic areas like amygdala, hippocampus, and thalamus which are shown in green are the most affected region. (B) Spatiotemporal evolution of empirical PSP pathology. Evolution of HTH seeded network diffusion exhibits the, limbic, striatal, and temporal areas as early affected regions, followed by somewhat slower diffusion into the frontal and parietal regions. This spatial sequencing predicted by the model is a close match with PSP tau progression. We selected three equidistant years at *t* = 7, 15, and 22 years between early stage (*t* = 0) with no diffusion through mature stage (*t* = 22). Here, time is arbitrary, and the use of "years" is meant for illustrative purpose. For both A and B, dots are color coded by lobe—frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green.

Figure 5 | Scrambled networks and PSP atrophy. (A) Histogram of correlation strength between network diffusion model (NDM) and PSP data over 2,000 shuffled networks. (B) Histogram of correlation strength between NDM and 2,000 shuffled PSP data over using unshuffled structural connectome. The true connectome was shuffled by symmetrically permuting its rows and columns randomly, and the NDM was evaluated for each shuffled network after bilateral hypothalamus seeding. The best *R* achieved by each model was recorded and entered into the histogram. The null models are distributed well below the true model, indicating that the latter is highly unlikely to arise by chance (*p* < 0.05).

(NDM) was played out for all time points with directional connectivity. Pearson's *R* was recorded at each time point between the model and PSP atrophy vector. As the diffusion time increases, more and more of the pathogenic agent escapes the seed region and enters the rest of the network. NDM seeded at bilateral regions using directional connectome indicates that the hypothalamus [shown by red curve in (A)] could be a plausible candidate for PSP disease seeding. The highest peak recorded with anterograde directional connectivity is *R* = 0.48 at *t* = 40 and with retrograde directional connectivity is *R* = 0.53 at *t* = 41. The characteristic intermediate peak is indicative of true pathology spread. The unit of model's diffusion time *t* is arbitrary. For the purpose of demonstration, we consider this arbitrary timescale parameter in terms of "years." (B,D) The point of maximum correlation *R* of NDM with measured atrophy using anterograde and retrograde directional connectome, respectively, were recorded and are shown with glass brains in (B,D). For all subplots, dots are color coded by lobe—frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green.

increasingly resembles empirical atrophy of **Figure 1** much so than with non-directional connectivity. We have selected three equidistant years at *t* = 13, 26, and 39 years between early stage (*t* = 0) with no diffusion through mature stage (*t* = 40) for anterograde mode and *t* = 14, 27, and 41 years between *t* = 0 and *t* = 41. Here, time is arbitrary, and the use of "years" is meant for illustrative purpose. As seen in **Figure 7B** with anterograde mode, initial seeding at HTH is followed by diffusion especially into the pallidum and caudoputamen, followed by limbic structures such as hippocampus and amygdala, then thalamus, then to middle temporal–frontal structures, and finally to the cortex. With retrograde mode (**Figure 7D**), initial seeding at HTH is followed by entorhinal, followed by other limbic structures such as hippocampus and amygdala, and thalamus, then to striatal, middle temporal-frontal structures, and finally to the cortex.

#### DISCUSSION

The current study applies network modeling to a rare disease as the first to systematically test hypotheses of disease spread in human subjects living with PSP. Other neurodegenerative diseases like FTD, AD, and ALS syndromes have been well-studied from a human brain network perspective, and have received more attention than equally important and related disorders (15), due to a lack of available human data on this rare disorder. The present analysis takes advantage of an unprecedented dataset of 60 subjects with PSP from a multinational prospective observation study called 4RTNI (see text footnote 1). Warren et al. proposed the term "molecular nexopathy" to refer to a coherent conjunction of pathogenic protein and intrinsic neural network characteristics (33). They enumerated a diverse set of

Figure 7 | (A) ND prediction from bilaterally seeded hypothalamus (HTH) with anterograde directional connectome. Glass brains of network diffusion model (NDM) seeded at the bilateral HTH at *t* = 40 years (corresponding to highest peak *R*) shows spatial patterns of ND progression from HTH using anterograde connectivity. Similar to Figure 4A, even with anterograde directional connectivity the most affected regions are subcortical striatal structures especially caudoputamen and globus pallidus, and limbic areas like amygdala and hippocampus, which are shown in green. (B) Spatiotemporal evolution of empirical PSP pathology with anterograde directional connectivity. Evolution of HTH seeded network diffusion exhibits the, limbic, striatal, and temporal areas as early affected regions, followed by somewhat slower diffusion into the frontal and parietal regions. This spatial sequencing predicted by the model is a close match with PSP tau progression. We selected three equidistant years at *t* = 13, 26, and 39 years between early stage (*t* = 0) with no diffusion through mature stage (*t* = 39). Here, time is arbitrary, and the use of "years" is only meant for illustrative purpose. (C) ND prediction from bilaterally seeded HTH with retrograde directional connectome. Glass brains of NDM seeded at the bilateral HTH at *t* = 41 years (corresponding to highest peak *R*) shows spatial pattern from HTH using retrograde connectivity. With retrograde directional connectivity the next most affected region is entorhinal (big sphere in cyan) and then limbic areas like amygdala and hippocampus which are shown in green. (D) Spatiotemporal evolution of empirical PSP pathology with retrograde directional connectivity. Evolution of HTH seeded network diffusion exhibits the limbic and temporal areas as early affected regions, followed by somewhat slower diffusion into the striatal, frontal, and parietal regions. We selected three equidistant years at *t* = 14, 27, and 41 years between early stage (*t* = 0) with no diffusion through mature stage (*t* = 41). Here, time is arbitrary, and the use of "years" is only meant for illustrative purpose. For all subplots, dots are color coded by lobe—frontal, purple; parietal, red; occipital, orange; temporal, cyan; and subcortical, green.

potential mechanisms by which molecular dysfunction might interact with the neural architecture to produce observed disease topography in neurodegenerative diseases, including dysfunction of synaptic function or maintenance, axonal transport or repair, or a result of downstream trophic or cell–cell signaling. Although a full encapsulation of these various protein-specific mechanisms into a network model like ours will require much more detailed data and studies, here we have made a beginning by explicitly considering one of the key ways in which misfolded protein species might interact with the extant network—i.e., *via* directional (as compared to non-directional) transmission. We explored two mechanistic aspects hitherto unknown about the canonical NDM of spread: (a) whether the NDM can apply to other common degenerative pathologies, specifically PSP, and (b) whether the directionality of spread is important in explaining empirical data.

Our first key contribution is to show that the mathematical graph theoretic NDM of trans-neuronal transmission can apply to other degenerative pathologies, specifically it can recapitulate observed PSP topography. The interest in PSP arises from its distinct spatial pattern in comparison to AD; hence the ability of NDM to explain PSP pattern contributes to the emerging notion that all neurodegenerative pathologies follow shared mechanisms of spread. Our second key contribution is the novel construction of a *directional human connectome*, exploiting the homologies in brain anatomy that exist between species as diverse as human and mouse. The importance of studying directional connectomes is that *in vitro* and *in vivo* mouse studies are increasingly revealing directional preference in the transmission of various misfolded proteins; however, conclusive data on directionality of each protein are not currently established (24). Our third key contribution is to show that the NDM fit with empirical data are enhanced by using the directional rather than the non-directional form of the human connectome. Certain aspects of PSP topography are better explained by the anterograde model of transmission (some to axonal terminal) than the non-directional model, and certain aspects of PSP are better explained by the retrograde mode of transmission. Overall, both directional models outperform the non-directional model, and retrograde mode gives overall the best fit. These intriguing results are further discussed below.

## Focal Seeding Followed by Network Spread Recapitulates PSP Atrophy

We tested whether network spread from focal seed sites would recapitulate PSP regional atrophy patterns. Although trans-neuronal transmission appears the most likely candidate, the exact mode of transmission is unknown. We first simulated spread between regions in a bidirectional manner, such that pathology has an equal chance of spreading in the anterograde or retrograde direction. To minimize model bias, we performed repeated seeding such that each brain region in turn has a chance to serve as the single onset region from which subsequent pathology transmission was modeled using the NDM. The HTH was found to be the best seed region, with the highest model correspondence of *R* = 0.42. The model performance seeded at HTH is significant in comparison with random network scrambling of the structural connectome (**Figure 5**) (*p* < 0.05), confirming that network organization is integral to disease spread in PSP, rather than an effect that has arisen by chance. Surprisingly, it was a better seed than any region in the striatum, even though striatal atrophy is usually the most prominent feature of the topography of PSP. The next best seed region was Pallidum, which is certainly plausible as converging post mortem and neuroimaging studies showing the striatum is the most affected region of pathology in PSP (8, 34).

#### Directional Structural Connectivity Model

Having identified structural connectivity as the most likely mechanism of transynaptic pathology spread in PSP, we sought to further improve our model by adding directionality. Clearly, a direct measurement of directionality of fiber tracts is impossible from current dMRI techniques, as water diffusion along fiber bundles does not respect cell polarity (soma to axonal terminal or *vice versa*). Instead, we exploit the well-known fact that homologous structures exist between many species, for some of whom we do happen to have anatomic connectivity data from painstaking tracer studies. A detailed mesoscale mouse connectome has recently become available from the Allen Institute (23). These data are fully directional, since it is based on retrograde AAV tracer experiments on a large number of mice. We, therefore, developed a novel technique whereby human and mouse brain atlas parcellations are used to define homologous brain structures, and mouse directional connectivity is transferred to non-directional human connectome.

Trans-neuronal transmission can have a distinct directional bias, such that misfolded protein species might follow anterograde or retrograde transport or signaling pathways. This is especially true in subcortical and striatal connections, which are known to be highly directional in comparison to corticocortical connections. Although little work has been done in PSP, tau pathology in AD differentiates between efferent and afferent connections (24). Hence, we tested the hypothesis that directional (whether anterograde or retrograde) process of spread along structural connectivity network will further improve model performance in PSP. The directional NDM results are shown in **Figures 6** and **7**, and confirm that either mode is a better fit than non-directional (**Table 2**). Interestingly, some aspects of PSP topography are better recapitulated with anterograde mode of transmission, especially striatal areas; whereas other aspects are better recapitulated by retrograde transmission, especially temporal and limbic areas. Quite unexpectedly, the HTH is the most likely source of pathology origin in all three modes of transmission, a result that raises many questions that should be addressed in future investigations.

# CONCLUSION

Collectively, these data show that the intrinsic architecture of the structural network mediates disease spread in PSP, most likely *via* a process of trans-neuronal transmission. The additional success of the directional network models suggests a simple process whereby local production of pathology starts off a process of axon-to-soma or soma-to-axon transmission, which due to the nature of the directional network topology, prominently accumulates in striatal and temporal areas. In this study, our focus was to build and apply directional NDMs in a purely data driven fashion, rather than a detailed exploration of the mechanistic aspects that govern directional transmission. We have shown that directionality of transmission may be an important aspect that network models should incorporate in future and related studies involving a wide range of neurodegenerative disorders. If confirmed by future, larger, studies, the apparent selective vulnerability and early seeding of PSP and other tauopathies in specific areas might plausibly entertain an explanation purely in terms of directional transmission, without requiring cell-type or region-specific vulnerability of brain region to pathology. To our mind, this would be the most parsimonious and economical reading of available data.

#### Limitations

Several methodological considerations should be considered when interpreting our results. The first are limitations of the NDM. The NDM is a first-order, linear, parsimonious model of diffusive spread that assumes that the structural connectivity network remains unchanged over the duration of the longitudinal analysis. Moreover, individual subject genetic repeat length, medication history and age of symptom onset were not available. These variables could have implications in the individual group wise analysis, when identifying each subject's seed region or determining individual rate of disease diffusion. Because this is the first study to empirically test multiple network models of pathology spread in PSP, it will benefit from independent replication. Future work elucidating striatal vulnerability as well as the effect of repeat length on disease spread is necessary.

#### REFERENCES


### AUTHOR CONTRIBUTIONS

SP coded, analyzed data, and wrote portions of the manuscript. CM extracted mouse and human connectome and wrote portions of manuscript. AR conceptualized the study, developed mathematical model, supervised image and statistical analysis, and wrote portions of the manuscript.

#### ACKNOWLEDGMENTS

The authors thank the Frontotemporal Lobar Degeneration Neuroimaging Initiative (FTLDNI) (http://4rtni-ftldni.ini.usc. edu/). The investigators at FTLDNI contributed to the design and implementation of FTLDNI and/or provided data, but did not participate in analysis or writing of this report. Data used in preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (http:// adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at: http://adni.loni.usc.edu/wp-content/uploads/how\_to\_apply/ ADNI\_Acknowledgement\_List.pdf.

# FUNDING

This research was supported in part by the following grants: R01 NS092802 and R01 EB022717 from the National Institutes of Health (AR).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Pandya, Mezias and Raj. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

#### *Chris Mezias1 \* and Ashish Raj1,2*

*1Department of Neuroscience, Weill Cornell Medicine of Cornell University, New York, NY, United States, 2Department of Radiology, Weill Cornell Medicine of Cornell University, New York, NY, United States*

While the spread of some neurodegenerative disease-associated proteinopathies, such as tau and α-synuclein, is well studied and clearly implicates transsynaptic pathology transmission, research into the progressive spread of amyloid-β pathology has been less clear. In fact, prior analyses of transregional amyloid-β pathology spread have implicated both transsynaptic and other intracellular- as well as extracellular-based transmission mechanisms. We therefore conducted the current meta-analytic analysis to help assess whether spatiotemporal amyloid-β pathology development patterns in mouse models, where regional proteinopathy is more directly characterizable than in patients, better fit with transsynaptic- or extracellular-based theories of pathology spread. We find that, consistently across the datasets used in this study, spatiotemporal amyloid-β pathology patterns are more consistent with extracellular-based explanations of pathology spread. Furthermore, we find that regional levels of amyloid precursor protein in a mouse model are also better correlated with expected pathology patterns based on extracellular, rather than intracellular or transsynaptic spread.

Keywords: connectomics, neurodegenerative diseases, neurodegeneration, computational modeling, proteinopathy, amyloid, amyloid spread

# INTRODUCTION

Spreading protein pathology is hypothesized by many scientists to underlie the spatiotemporal pattern of lesions (1), regional neuronal and volume loss (2, 3), as well as the progression of the presentation of symptoms (4) in degenerative diseases. Amyloid-β, the misfolded protein cleaved from amyloid precursor protein (APP), is the key constituent of amyloid plaques seen in Alzheimer's disease (AD) and is thought to contribute to the observed toxicity and cell death along with tau (5). More current research has shown both misfolded tau and α-synuclein exhibit progressive protein pathology with prominent lesions progressing in fairly consistent spatiotemporal patterns (6, 7). However, the spread of amyloid-β pathology has not been neatly characterized into stages of pathology, with little evidence of focal seeding as is seen for tau or α-synuclein. However, studies tracking the spatiotemporal development of amyloid pathology have noted early development generally throughout the neocortex (8) with the deposition perhaps initially focused in posterior cortical areas (9). The staging research for amyloid pathology development generally shows more variance than do tau or α-synuclein staging for particular diseases and generally develops in a less focal manner.

#### *Edited by:*

*Claudio Soto, University of Texas Health Science Center at Houston, United States*

#### *Reviewed by:*

*Corinne Lasmezas, Scripps Florida, United States Rodrigo Morales, University of Texas Health Science Center at Houston, United States*

#### *\*Correspondence:*

*Chris Mezias chm2062@med.cornell.edu*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 16 August 2017 Accepted: 20 November 2017 Published: 18 December 2017*

#### *Citation:*

*Mezias C and Raj A (2017) Analysis of Amyloid-β Pathology Spread in Mouse Models Suggests Spread Is Driven by Spatial Proximity, Not Connectivity. Front. Neurol. 8:653. doi: 10.3389/fneur.2017.00653*

**84**

The more focal and ordered staging characterizations for tau and α-synuclein, and perhaps underlying biological progression of pathology, have led to more successful study of the potential mechanisms of pathology spread for these proteins. For example, a whole litany of patient (3, 10) and mouse modelbased (11–15) studies has looked at whole brain and cellular level resolutions and determined that the evidence, regardless of resolution of study, is most consistent with a transsynaptic and intracellular mechanism of tau spread. In fact, our own recent research demonstrates that tau pathology in mouse models progression mirrors the network of mouse fiber tracts (16). Similarly, regarding α-synuclein, there is strong evidence pointing toward transsynaptic and intracellular mechanisms of spread, including transfer of pathology from Parkinson's and LBD patient neurons into grafted and implanted cells (17), the presence of both misfolded and normal state α-synuclein at the synapse (18), and controlled cell-to-cell transfer experiments in mice (18) demonstrating that α-synuclein spread occurs transsynaptically by virtue of initiation into a set of interconnected cells. The success in characterizing the spread mechanisms and abilities of tau and α-synuclein are now leading to efforts to slow down the progression of proteinopathy as a potential avenue for treatment with conditions based on both of these proteins (19, 20). However, the study of mechanisms of cell-to-cell or region-to-region transfer for amyloid-β pathology has produced mixed results.

The spread mechanisms of amyloid-β pathology, in terms of causing a misfold in healthy amyloid-β cleavage product, are well characterized and relatively similar to those for tau. The most accepted model, supported by evidence from patients (21), mouse model, and cell lines (5, 22), is that misfolded amyloid-β can recruit and induce misfolding in healthy protein product amyloid-β. However, unlike the cases of tau and α-synuclein, this finding has not helped gain more clarity in how amyloid-β pathology spreads transcellulary and transregionally. Some studies have proposed extracellular mechanisms of spread (23, 24), while other research has implicated internalization of amyloid-β as a key element in this protein product becoming pathological (5), with other researchers building on this finding and demonstrating potential transsynaptic and intracellular spread (25, 26). Given the failure of clinical trials targeting already misfolded amyloid-β (27, 28), we undertake the present research to help resolve inconsistency in understanding the pathology spread mechanisms of amyloid-β, with an eye toward beginning to make targeting of amyloid-β pathology spread mechanisms an even more viable therapeutic option.

In the present study, we hope to add clarity to the amyloid-β pathology spread debate by analyzing available data on the patterns of regional amyloid pathology severity and spatiotemporal progression in mouse models. We first demonstrate amyloid-β pathology patterns in mouse models are not mirrored by the mouse connectome, while those of tau pathology, accordant with prior research (16), are. Furthermore, a model of connectivitybased amyloid-β pathology spread fails to recreate spatiotemporal patterns of amyloid-β pathology development in a mouse model. However, we demonstrate that a model of extracellular, spatial proximity-based diffusion (SPD), does recreate the regional severity and spatiotemporal progression patterns of amyloid-β pathology. Moreover, spatial proximity to amyloid-β pathology initiating regions serves as a good proxy for predicting the regional severity of said pathology. We conclude that the spread of amyloid-β pathology is predominantly driven by extracellular means and that the basis for regional vulnerability to amyloid-β is how spatially proximal that region is with those already exhibiting pathology.

# MATERIALS AND METHODS

#### Study Selection for Datasets

Studies were found from a literature search of Web of Science and Google Scholar and had to meet the following criteria to merit inclusion: they had to be from 2005 onward to assure their methodology, and measurements were somewhat current with the standards in the field, they had to include data about the regional distribution of amyloid pathology from at least 5 different regions, and they had to be from mouse models using the APPswe mutation and had to assess endogenous pathology development, rather than pathology development induced from exogenous amyloidogenic seeds. We note that this criterion excludes some important studies, such as Ye et al. (24), but uses it nonetheless for consistency in our analyses. Amyloidopathy data here come from Harris et al. (25) (seven total regions and three timepoints), Bero et al. (26) (five total regions, one timepoint), and Lim et al. (29) (five total regions and one timepoint). Tau data cited for comparison purposes in Figure 2 came from Ahmed et al. (12) (5 total regions, 1 timepoint), Bolmont et al. (30) (5 total regions, 1 timepoint), and Clavaguera et al. (11) (11 total regions, 3 timepoints, only the final one used); please refer to the original citations for specific methodological questions about any of the above studies.

## The Allen Brain Institute Mouse Connectivity Atlas and Generating the Directed, Undirected, and Spatial Distance Graphs of the Mouse Brain

Connectivity and spatial distance data were taken from the supplementary dataset published along with the mesoscale mouse connectome from the ABI (31). Total projection volume between regions was generated by multiplying element-wise by the rows the connectivity matrix times the number of voxels in each seeding region. Regional distances used were the Cartesian distances between the centers of mass of any region pairs. For an anatomical illustration of the brain represented as a connectivity graph, see **Figure 1A**; directional connectivity from to and from the entorhinal cortex (EC), used as an example region, can be seen in **Figure 1B**. For an example of a spatial distance graph as a matrix, see **Figure 1C**.

#### Network Laplacian and Its Eigenvectors

For a given connectivity *C* or distance matrix *S*, we define a Normalized Laplacian Matrix, as in Raj et al. (2), *L*:

*L* = −*I D*r c ⋅*C D*⋅

matrices 1,000× when we use DNT or SPD to provide another metric of how unexpected our observed *r*-squared values would be given chance alone.

where *D*r and *D*c are diagonal degree matrices of the sum totals of the rows and columns, respectively, and *I* signifies the identity matrix. In our notation for the connectivity matrices, *C* represents the retrograde transmission mode, and the corresponding anterograde mode is then simply given by *C*T , and the bidirectional mode by (*C* + *C*T)/2. The distance matrix was obtained from a separate dataset appended to Ref. (31) and is referred to in this manuscript using *S*. Normalized Laplacian matrices were generated for each mode with the above equation.

#### Creating a Dynamic Model of Amyloid Pathology Spread Over Time

A previous graph theoretic model of pathology spread in AD throughout a brain network was shown to be predictive of future patterns of disease progression (2). The principle is to seed a graph node corresponding to a brain region with an arbitrary value and then model the diffusion of the disease factor throughout the network *via* the Network Diffusion equation:

$$X(t) = \exp(-\beta Lt) \cdot X(0). \tag{1}$$

This models the long-range patterns of spread of the protein pathology at any time *t* as a product of the initial seeding pattern *X*(0), and the so-called diffusion kernel exp (−β*tL*), with diffusion constant, β (2).

The major differences with previous network diffusion model are twofold: (1) we are for the first time using this model with a directed brain network and (2) we are interested in total pathology accumulation over time, which we model as a summative process:

$$X(i) = \exp(-\beta L \Delta t) \cdot X(i-1) + X(i-1). \tag{2}$$

We use Eq. 2 to calculate, for any point in time, the deposition of tau, amyloid, and neuronal atrophy across the brain regions represented in the network. The Eq. 1 in the present study will be referred to as the standard undirected network diffusion model, or NT, while our modified directional version in Eq. 2 will be referred to as directed network transmission, or DNT. We discuss this equation modeling spatial pathology diffusion in the Section "Spatial Diffusion Modeling As an Alternative Model to Network Diffusion."

#### Spatial Diffusion Modeling As an Alternative Model to Network Diffusion

We additionally created a spatial diffusion model as a comparison or alternative hypothesis to the network diffusion model. The spatial diffusion model was based on the same fundamental network diffusion Eq. 2 explicitly stated above. The difference between DNT (and NT) and spatial diffusion in the present study is that the network for spatial diffusion is a matrix where each entry in the matrix *Si,j* is the reciprocal of the Cartesian distance between the center of mass of each GM region included in the Allen Institute's mouse connectivity atlas. Using this distance matrix, *S*, rather than the connectivity matrix *C*, we ran the diffusion equation stated above in (2) to get a model approximating diffusion based on spatial proximity, which will be referred to as SPD.

### Comparing the DNT/NT and SPD Models with Previously Published Results and Examining the Question of Seeding

We ran the network and spatial diffusion Eq. 2, through the number of iterations, in months, given in each study. If a study measured pathology at 6 months, we ran the model through 6 iterations, and if a study measured pathology at 9 months, we ran the model through 9 iterations. The implicit assumption here is that amyloidopathy spreads in all datasets at the same rate. While these datasets were obtained with different mouse models potentially exhibiting different amyloid strains, we do not have enough data to assess the kinetics or speed of amyloid pathology spread in individual cases. Hence, we decided to impose a minimal and general set of assumptions across the board. In these iterations, we used Δ*t* = 1 month and modified β to achieve the optimal match with the data, as the empirical diffusion constant for various pathologies in AD is *a priori* unknown. Example, βt curves for a range of β values, where the model both shows behavior that is predictive of proteinopathy spread and non-predictive of proteinopathy spread can be seen in **Figure 1C**. We show anatomical examples of a graph that is the basis of NT/DNT in **Figure 1A**, and a model of the brain demonstrating the locations of regional centers of mass, which form the basis of the network used in SPD, in **Figure 1B**. To compare whether DNT/NT or SPD performed better at recapitulating amyloid-β pathology, we calculate Δ*r*-squared values, which are the difference between the *r*-squared value between the amyloidopathic seed (no diffusion) and empirical amyloidopathy, and the same *r*-squared value calculated between the best fit value of the βt model constants of diffusion time and empirical amyloid-β pathology. An example of this calculation can be seen in **Figure 1C**. Given the relatively small sample size of regions quantified from Harris et al. (25), we performed an additional assessment for the significance of regressions that were significant at an alpha of *p* < 0.05: we randomized either the connectivity or interregional distance matrices from the ABA 1,000 times and ran our models to produce simulated data generated with randomized networks, generating a distribution *r*2 -values to use as a comparison with *r*<sup>2</sup> -values from our models run using empirical networks or distances. An example of the results of this process can be found in **Figure 1D**.

#### RESULTS

The current paper attempts to perform a whole-brain scale, meta-analytic study of amyloid-β pathology spread mechanisms. The motivation behind such research is clear: while research into other proteinopathies with prion like qualities, such as tau pathology, has demonstrated a clear mechanism and clear pattern of transmission, such as transsynaptic spread (1), efforts to characterize misfolded Aβ spread are conflicted (24–26). Indeed, understanding misfolded amyloid spread mechanisms as Aβ proteinopathy develops could prove a fruitful avenue of research for future disease course modifying treatments if further amyloid transmission can be prevented. We first demonstrate that given three studies regionally quantifying tau pathology, we see clear evidence of transsynaptic transregional spread, but that such evidence is lacking with respect to amyloid pathology transmission in another three studies. We then demonstrate, using models of both progressive spread over the brain's connectivity network (NT and DNT) and *via* diffusion based on spatial proximity to already affected regions (SPD), that diffusion to spatially proximate regions is a much better characterization of both the development of Aβ pathology and of the regional pattern of Aβ precursor protein, APP, levels. The present results are discussed in detail in the subsections below.

### Amyloid-**β** Pathology, Unlike Tau Pathology, Is Not Well Characterized in Mouse Models by the Brain's Connectivity Network

Prior use of the brain connectome graph metrics and models (such as DNT and NT) have successfully, characterized the spatiotemporal development of regional volume loss (2) and metabolic deficits (3) in patients, and the spatiotemporal proliferation of tau pathology in mouse models (16). We use the first two eigenvectors of the mouse brain connectome to assess whether amyloid-β pathology can be recapitulated by the properties of mouse brain wiring in a way that resembles the accurate characterization of tau pathology. We find that, while, accordant with prior research (16) relative regional severities of tau pathology across all three studies highlighted here [**Figures 2A–C,G,H**; **Table 1**; (11, 12, 30)] are all strongly recapitulated by at least one the first two eigenvectors of all of the mouse connectomes used (anterograde, retrograde, and bidirectional), only one study regionally quantifying amyloid pathology produced any such significant results [**Figures 2D,G,H**; **Table 1**; (26)]. The results comparing the eigenvectors from our mouse connectomes with our other two amyloid datasets (25, 29) showed no significant match between where, based on transsynaptic spread, one would expect the most severe amyloid pathology and where the empirical dataset indicated the most severe pathology, across all mouse connectomes (**Figures 2E–H**; **Table 1**). The present results confirm our prior analyses that tau pathology patterns in mouse models can be well characterized using mouse connectomes but also reveal that no such strong relationship exists between amyloid-β pathology patterns in mouse models and the mouse connectomes.

### Spatiotemporal Amyloid-**β** Pathology Is Better Characterized by a Model of Spatial, Rather Than Connectome Based, Diffusion

We next assess how well our connectivity based DNT/NT model recapitulates the longitudinal spatiotemporal development of

Figure 2 | *u1* and *u2* consistently predict regional tauopathy, but not amyloidopathy, severity. (A–F) The best regressions, by *r* 2 -values, using the anterograde, retrograde, and bidirectional Allen Institute mouse brain connectivity networks, plotted, of *u1* or *u2* vs. deposition data from several studies; (A) is the regression between anterograde *u1* and Ahmed et al. (12), (B) is between bidirectional *u2* and Bolmont et al. (30), (C) is retrograde *u1* vs. Clavaguera et al. (11), (D) is the regression of retrograde *u1* and Bero et al. (26), (E) is between bidirectional *u2* and Harris et al. (25), while (F) is Lim et al. (29) vs. retrograde *u2*. All *r* 2 -values for all *u1* and *u2* regressions vs. each dataset can be found in Table 1. (G) The scatter plot of *u1* and *u2*, for each direction, *r* 2 -values with actual data indicates that at least one direction's eigenvectors, if not multiple, show high correspondence with tau deposition, but that only the data from Bero et al. (26), is correlated with the eigenvectors from the directed connectivity graphs. (H) This bar graph illustrates the consistent pattern of tau and eigenvector correspondence, and the relative inconsistency of correlations between amyloid and eigenvectors. The first *x*-label in (H) refers to the proportion of eigenvectors, regardless of study and first of second eigenvector, yielding a *p* < 0.05, the second those yielding *p* < 0.10, and the third refers to the % of studies used here having any eigenvector yielding significant results at the threshold *p* < 0.05.


*In the table above, all r-squared values between both of the first two eigenvectors of each connectivity network, anterograde, retrograde, and bidirectional, and empirical tau and A*β *pathology patterns are reported, with the study indicated in the leftmost column. ^p* < *0.10, \*p* < *0.05, \*\*p* < *0.01.*

amyloid-β pathology in our one multi-timepoint mouse amyloid dataset from Harris et al. (25). We find that, at no measured timepoint either early or late, does anterograde or retrograde DNT or bidirectional NT successfully predict amyloid pathology patterns (**Table 2**). In fact, at all timepoints, including the final measured timepoint, DNT and NT fail to add any information beyond the EC seeding reported by Harris et al. (25) for predicting the development of amyloid-β pathology (**Figure 3A**). The regional slopes of pathology severity increase show an identical pattern of results (**Figure 3B**).

However, the amyloid-β pathology patterns at all timepoints, including the final measured timepoint, are significantly recreated using the SPD model of pathology spread, which is based on spatial proximity to the EC, where pathology initiated (**Figures 3A,C**; **Table 2**). Moreover, the regional slopes of amyloid-β pathology severity increase are also significantly predicted by SPD (**Figures 3B,D**; **Table 2**). Due to the relatively small sample size of mice and quantified regions from Harris et al. (25), we tested whether the results from the SPD model were due to chance by assessing whether significant results could be obtained using randomized distances between regions; we found that for the final measured timepoint of amyloid-β pathology, SPD using empirical distances between brain regions outperformed 91% of the 1,000 randomizations, and 93% of the randomization when predicting the regional severities of rate of amyloid-β pathology increase (**Figures 3E,F**).

Table 2 | DNT and NT do not predict either Aβ pathology or amyloid precursor protein (APP) levels, at any timepoint measured in Ref. (25), but spatial proximitybased diffusion (SPD) does across all timepoints.


*This table reports* Δ*R-squared values, which is the change in r-squared between baseline and the best-fit each model produce produced. Only SPD produces values* > *0 across all measurements.*

## Regional Levels of APP Are Better Predicted *via* a Model of Spatial, Rather Than Connectome Based, Diffusion

The data from Harris et al. (25) also includes regional quantification of amyloid-β precursor protein (APP) levels using anti-hAPP antibody. While we find that amyloid-β pathology proliferation is not well explained by connectome-based spread models, but rather by spatial proximity-based models, we want to test whether regional levels of APP can be predicted using the connectome. However, our results again indicate that, while retrograde DNT does give a significant prediction of APP levels (**Figures 4A,B**), anterograde DNT and bidirectional NT do not,

Figure 3 | Amyloidopathy is not predicted by DNT or NT but rather by spatial proximity-based diffusion (SPD). (A) The βt vs. *r* 2 -value curve comparing the three DNT/NT and the SPD models' predictions to regional rate of increase and (B) end-state regional amyloid plaque burden data from Harris et al. (25). Note all three DNT/NT models' curves are non-predictive. (C) The plotted regression of the SPD model against the regional rate of increase and (D) end-state regional amyloid burden data. The *r* 2 -value from the SPD model using the Cartesian distance between region matrix performs better than most, but not all, randomized distance matrices for both (E) regional rate of pathology increase and (F) end-state regional amyloid burden data. Color legends are included. All exact *r*-squared values can be found in Table 2.

and SPD outperforms even retrograde DNT at recreating relative regional levels of APP (**Figures 4A,C**; **Table 2**). Furthermore, the significant match between the predictions of retrograde DNT and regional APP levels appears to be driven by a single outlier region (**Figure 4B**), calling the validity of these results into question. Akin to our results about the spatiotemporal development of amyloid-β pathology in the subsection above, we find that regional levels of APP are again best predicted by a model of diffusion based on spatial proximity with regions exhibiting high APP levels, more accordant with extracellular diffusion mechanisms than connectivity based or transsynaptic based proliferation. To emphasize the above results, we also include an anatomical illustration of the regional rates of suspected pathology increase, from the measurement of regional APP levels until the final quantification of amyloid-β pathology can be seen in **Figure 4D**; note that SPD is a better model for spatiotemporal amyloid-β pathology patterns than is DNT/NT.

#### DISCUSSION

While research into other protein pathology spread mechanisms in mouse models, such as work regarding synucleinopathy (18) and tauopathy (11–16) reveal clear pathology development patterns mirroring anatomical connectivity networks and strongly suggest transsynaptic spread, no such clarity exits regarding spread patterns and mechanisms for misfolded amyloid-β. What is clear however is that amyloid pathology does spread outward in a manner that is highly dependent on the brain region into which amyloidogenic seeds are placed, even when using the same amyloid pathology source and transgenic mouse model (32, 33).

We accordingly undertook the present analysis to determine whether we could find clear evidence of transsynaptic spread or pathology development mirroring connectivity networks for amyloid-β, as some prior research has claimed (25, 26), or whether some other pattern and mechanism of spread better explains amyloid-β pathology development (24, 29). We first demonstrate that, especially compared with tau pathology development, amyloid-β pathology, across studies, is not most pronounced in regions that are most likely to accumulate pathology if spread is *via* the mouse anatomical connectivity network (**Figure 2**; **Table 1**). We next demonstrate that, in assessing our lone longitudinal dataset (25) a model based on spread to spatially proximal (SPD), rather than a model assuming transmission to interconnected (NT/DNT) regions, better recreates the spatiotemporal patterns of amyloid-β pathology development and the relative regional severities of amyloid-β plaques (**Figures 3** and **4**; **Table 2**). These results and their further implications are discussed in detail in the subsections below.

## Pathology Spread into Spatially Proximal, but Not Necessarily Highly Interconnected Regions, Implies Extracellular Spread

The present analysis shows results suggestive of amyloid-β pathology spread being predominantly extracellular, rather than intracellular. Modeling assuming transsynaptic spread failed to significantly recapitulate empirical patterns of amyloid-β pathology (**Figures 2**–**4**; **Tables 1** and **2**), while modeling based on spatial proximity significantly recreated observed spatiotemporal amyloid-β pathology development patterns (**Figures 3** and **4**; **Table 2**). While prior studies implicate transsynaptic amyloid-β pathology transmission (25, 26), these findings are contradicted by other research demonstrating a lack of evidence for transsynaptic or even any form of intracellular transneuronal and transregional spread (24, 29). By contrast, studies of other pathological protein species' mechanisms of spread, such as pathological tau, consistently implicate transsynaptic spread (11–16) with little evidence for other mechanisms.

Extracellular diffusion resulting in amyloid-β pathology transmission to the areas of the brain most spatially proximal to those already exhibiting amyloidopathy fits with the known characteristics of amyloid-β as a protein. First, APP, the precursor protein for amyloid-β, is a transmembrane protein with a large extracellular domain, and it is from this large extracellular domain that amyloid-β is formed as a cleavage product from APP (34). Post-cleavage, amyloid-β generally exists in the extracellular space, even as a healthy protein (35). Finally, amyloid-β in its pathological states also forms plaques almost exclusively extracellularly (36). None of this is to say that pathological amyloid-β might not gain capabilities of becoming internalized into neuron or spreading transneuronally or transsynaptically, as some authors have suggested (25). But given our meta-analytic results suggesting no mathematical basis for transsynaptic spread or any spread based on the brain's connectivity network, and suggesting that pathology is more likely to spread into regions spatially close to rather than heavily interconnected with already affected areas, we posit that our present work, strongly implies some mechanism of pathology spread that is extracellular. We hope that our work serves as a starting point for more heavily quantitative investigations of amyloid-β pathology transmission with an aim toward clearing up the inconsistent results about how and where amyloid-β progresses around the brain.

Given that amyloid patterns are known to depend on exogenous seeding site (32, 33), a purely spatial mode of spread would suggest a radial distribution of pathology spread, which was not observed in our data. However, there is very little evidence that initial amyloid seeding is in fact focal; instead early amyloid pathology is found diffusely in the neocortex and then spreads to subcortical areas (8). Quite likely, amyloid pathology could

be caused or initiated by metabolic deficits that target particular regions, leading to amyloid initiation in broad areas of the brain, for instance in the default mod network (9). Due to broad initiation, a radial spread from a focal source may not be manifested under spatial spread. Nonetheless, it is unlikely that spatial spread *via* extracellular spaces is the only process involved in amyloid pathology ramification. Several additional factors might play a role beyond pure spatial spread, and future research will be needed to identify and assess them.

Also of note is that our SPD model accurately recreated regional APP levels from Harris et al. (25). Given that APP is a transmembrane protein (34), spread based on spatial proximity, such as through the extracellular space, seems unlikely. First, it is possible that the anti-hAPP antibody used in the study (25) bound to a non-Aβ fragment of APP following cleavage, sAPPα, which akin to Aβ is also often released into the extracellular space (34) and could therefore plausibly spread into spatially proximal regions. Aβ also has known reactions with full-length transmembrane APP (34), and so spatially spreading pathological Aβ could influence regional APP levels by upregulating APP expression or by inducing the release of transmembrane APP into the extracellular space. These hypotheses for why a model of diffusion into spatially proximal regions accurately captured empirical APP levels, however, are conjecture at this point. Akin to the discussion of factors beyond spatial proximity contributing to the spatiotemporal development of Aβ pathology, above, elucidating why APP levels mirror Aβ pathology spread patterns in at least a pathological mouse model, will require careful future research.

#### Extracellular Spread of Amyloid-**β** Pathology Fits with Known Synaptic Problems Resulting from Its Pathology

Deficits in cellular function directly attributable to amyloid-β in mouse models are generally synaptic deficits (37, 38). Given the propensity for amyloid-β plaques to accrue at synapses (39), it is not surprising that such an accumulation of misfolded protein species would causes deficits in the ability to transmit and receive action potentials (37, 38). However, our findings implying extracellular spread based on spatial proximity to already affected regions, coupled with the above work demonstrating the propensity for amyloid-β pathology to cause synaptic and cell-signaling deficits raise an important question: why should misfolded amyloid-β, if it is spreading and circulating extracellularly, preferentially cause issues at synapses, rather than at other areas along the very large neuronal membranes?

Work whose data are used in the present analysis regarding spatiotemporal amyloid-β pathology development provides a possible explanation. Some recent studies demonstrate amyloid-β pathology appears more likely to accrue in more active regions (26) close to those already exhibiting pathology. Furthermore, this work suggests electrical signaling between neurons may act as an attractant for pathological amyloid-β (26), inducing more amyloid-β to accumulate and therefore form plaques at cellular locations where electrical current leakage is the greatest: at the synapse (39). If, as our results indicate, amyloid-β pathology follows a pattern of transiting to spatially proximal, rather than interconnected areas, likely *via* some form of extracellular diffusion, then electrical current leakage from the synapse could bias amyloid-β pathology spread toward those regions that are both spatially proximal and more likely to have high signaling (and therefore electrical) activity. This could further explain why, even though SPD consistently outperformed DNT/NT, it was still not as strong a model for predicting spatiotemporal amyloid-β pathology development as DNT/NT are for modeling tau pathology (16). Amyloid-β pathology spread patterns might be best explained by a model incorporating both spatial proximity with already affected regions and the average signaling activity of those spatially proximal regions; our research suggests this as a topic for further study in elucidating amyloid-β transmission.

### Extracellular Spread of Amyloid-**β** Pathology Is Complicated by Hypotheses of Interactions with Known Intracellular Pathological Proteins, Tau, and **α**-Synuclein

While amyloid-β is mainly an extracellular protein, other protein pathologies known to be comorbid with amyloid-β pathology are caused by intracellular proteins tau and synuclein. A potential complication for the model of extracellular spread of amyloid-β pathology into the most spatially proximal regions is complicated by evidence demonstrating potential interactions between amyloid-β and tau (30, 40, 41) and known co-occurrence of misfolded amyloid-β plaques embedded within Lewy Bodies composed primarily of misfolded synuclein (42). Synuclein is almost exclusively present at synaptic terminals and tau is an almost exclusively axonal protein, so given the prior work cited above showcasing known amyloid-β propensity for accumulating around and causing issues at synapses (37–39), a predominantly extracellular mechanism of misfolded amyloid-β spread does not necessarily negate the possibility of the interaction of these proteins. However, our present work does not provide a satisfying explanation as to how pathological proteins that are on the one hand predominantly intracellular and on the other are the likely extracellular transiting amyloid-β are likely to interact. We nonetheless feel it is important to point out the limitation to our current work and to encourage future research to pursue this question.

# AUTHOR CONTRIBUTIONS

No animal or human research subjects were directly used in the present research. All authors read and approved of the manuscript.

# ACKNOWLEDGMENTS

We would like to acknowledge our colleagues in the IDEAL lab, especially Dr. Amy Kuceyski and Ms. Eve LoCastro, for their guidance regarding connectivity models and generating graphics.

#### FUNDING

We would like to acknowledge the NIH for their generous support of the IDEAL lab *via* grants, numbers: R01NS092802 and R01EB022717.

#### REFERENCES


signaling pathway. *J Neurosci* (2007) 27(11):2866–75. doi:10.1523/JNEUROSCI. 4970-06.2007


42. Irwin DJ, Lee VM-Y, Trojanowski JQ. Parkinson's disease dementia: convergence of α-synuclein, tau and amyloid-β pathologies. *Nat Rev Neurosci* (2013) 14(9):626–36. doi:10.1038/nrn3549

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer RM and handling Editor declared their shared affiliation.

*Copyright © 2017 Mezias and Raj. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Estimating Memory Deterioration Rates Following Neurodegeneration and Traumatic Brain Injuries in a Hopfield Network Model

#### Melanie Weber\*, Pedro D. Maia and J. Nathan Kutz

*Department of Applied Mathematics, University of Washington, Seattle, WA, United States*

Neurodegenerative diseases and traumatic brain injuries (TBI) are among the main causes of cognitive dysfunction in humans. At a neuronal network level, they both extensively exhibit focal axonal swellings (FAS), which in turn, compromise the information encoded in spike trains and lead to potentially severe functional deficits. There are currently no satisfactory quantitative predictors of decline in memory-encoding neuronal networks based on the impact and statistics of FAS. Some of the challenges of this translational approach include our inability to access small scale injuries with non-invasive methods, the overall complexity of neuronal pathologies, and our limited knowledge of how networks process biological signals. The purpose of this computational study is three-fold: (i) to extend Hopfield's model for associative memory to account for the effects of FAS, (ii) to calibrate FAS parameters from biophysical observations of their statistical distribution and size, and (iii) to systematically evaluate deterioration rates for different memory-recall tasks as a function of FAS injury. We calculate deterioration rates for a face-recognition task to account for highly correlated memories and also for a discrimination task of random, uncorrelated memories with a size at the capacity limit of the Hopfield network. While it is expected that the performance of any injured network should decrease with injury, our results link, for the first time, the memory recall ability to observed FAS statistics. This allows for plausible estimates of cognitive decline for different stages of brain disorders within neuronal networks, bridging experimental observations following neurodegeneration and TBI with compromised memory recall. The work lends new insights to help close the gap between theory and experiment on how biological signals are processed in damaged, high-dimensional functional networks, and towards positing new diagnostic tools to measure cognitive deficits.

Keywords: neurodegenerative diseases, traumatic brain injuries, Hopfield neuronal network, associative memory encoding, memory impairments, axonal swellings

# 1. INTRODUCTION

Neurodegenerative diseases and traumatic brain injuries (TBI) are responsible for an overwhelming variety of functional deficits, both cognitive and behavioral, in animals and in humans. Memory impairment, which is the focus of this work, is a particularly pernicious consequence for those affected. The pathophysiology induced by these brain disorders is usually complex, with key

#### Edited by:

*Yasser Iturria Medina, Montreal Neurological Institute and Hospital, Canada*

#### Reviewed by:

*Roberto C. Sotero, University of Calgary, Canada Paule-J. Toussaint, McGill University, Canada*

\*Correspondence:

*Melanie Weber melwe@uw.edu; pmaia@uw.edu; kutz@uw.edu*

#### Specialty section:

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neuroscience*

Received: *01 July 2017* Accepted: *26 October 2017* Published: *09 November 2017*

#### Citation:

*Weber M, Maia PD and Kutz JN (2017) Estimating Memory Deterioration Rates Following Neurodegeneration and Traumatic Brain Injuries in a Hopfield Network Model. Front. Neurosci. 11:623. doi: 10.3389/fnins.2017.00623* effects of the injuries occurring at small spatial scales that are currently inaccessible by non-invasive diagnostic techniques. Indeed, a hallmark pathology is the abundance of diffuse, Focal Axonal Swellings (FAS) which compromise spike train encodings in neuronal networks. Thus, there is a broad need to understand how neuronal pathologies which develop at a cellular level compromise the functionality of a network of neurons responsible for cognitive function. In this article, we extend a well-established computational model of associative memory, i.e., the Hopfield network model (Hopfield, 1982, 1984), to incorporate FAS pathologies implicated in many brain disorders, providing novel metrics for quantifying memory impairments and functional deficits. While random neuron and synapse loss have been studied previously in memory models, this work is the first to integrate them with a recently developed theory of impaired neuronal responses due to FAS and their effects to spike train coding (Maia and Kutz, 2014a,b, 2017; Maia et al., 2015). While it is obvious that damaging any network will lead to compromised functionality, our quantification of memory decline progression could lead to new metrics for understanding memory impairments in a variety of leading neurodegenerative diseases.

Our computational models have several limitations and are of course only a simplified proxy for the original neural circuity involved in memory tasks. Still, they can provide an interesting tradeoff as we address questions that would be otherwise impossible in a clinical or experimental setting. We define unambiguously and in a precise mathematical way quantities such as memory overlap, successful memory recognition, memory confusion, significance of the memory classification, and noise-handling ability. These variables are all incorporated in our recognition score, a new metric that offers a more complete description of the network' s performance. We can then systematically increment the injury level, the noise level, randomize targeted neurons, change FAS parameter distributions, and address either highly correlated memory sets (faces) or random and uncorrelated ones. Consequently, we are able to quantify the decline in memory performance as a collective function of these variables. Moreover, we run enough simulations to control for each variable, achieve statistically significant results, and estimate the variability between different realizations. We show that as the injury progress, the first affected network functionality is its noise-handling ability, and that explicit memory-confusion occurs only later on. In this regard, the FAS parameters and distributions are crucial to determine when and how fast confusion in associative-memory should be expected.

The seminal work of John Hopfield is still largely used to model memory association (Hopfield, 1982, 1984; Hopfield and Tank, 1985). There, meaningful stimuli are encoded, stored, and later recalled in response to some cue. And although this process improves the interpretation of subsequent stimuli that share common features with previously stored concepts, what ultimately governs the memory recall abilities are the coordinated exchanges of electrical signals between neurons in network structures. Several pathological effects, most notably FAS, can jeopardize this critical electrical activity in the brain, making memory performance particularly sensitive to common disorders. Given its ubiquity across many neurodegenerative diseases and traumatic brain injury, understanding the role of FAS in altering spike train encodings is of paramount importance, particularly on a network level where cognitive functionality occurs. The purpose of this work is to help close this gap and study the effects of FAS to network models that exhibit plausible memory retrieval capabilities.

# 2. RESULTS

Some of our contributions are theoretical and methodological as we combine, for the first, time two models: (a) Hopfield's network model for associative memory (Hopfield, 1982, 1984) and (b) the theory of impaired neuronal responses due to FAS (Maia and Kutz, 2014a,b, 2017; Maia et al., 2015). Full details are provided in the Materials and Methods section and in the supplemental materials. It is important to mention that the previous FAS results by Maia and Kutz concern individual swellings only. Given the geometrical parameters of a particular swelling, they were able to estimate the corresponding type of impaired neural response (transmission, filtering, reflection, or blockage). In other hand, Hopfield's memory model is a network model. To study injuries at its level, we must provide plausible distributions of swelling effects to the injured population. Thus, the FAS parameters (pie-charts) calibrated from the experiments of Wang et al. and of Dikranian et al. (detailed ahead) are also important methodological contributions of this work. In what follows, we will illustrate these ideas for two different memory-retrieval tasks: (i) a face recognition task for highly correlated memories, and (ii) a memory retrieval task for random, uncorrelated memories. For each task, we tried to justify our simulation parameters as much biologically and statistically as possible, but since there were still some arbitrary selections, we made all our codes available to allow further parameter explorations.

# 2.1. Face Recognition Decline on Impaired Associative-Memory Networks

We calibrate a Hopfield-like neuronal network (Hopfield, 1982, 1984) to perform face recognition tasks. See **Figure 1** for a detailed schematics. Our sample memory space consists of five human facial images (1, 044 × 1, 341 pixels), modified from the MIT faces database (Weyrauch et al., 2004) 1 . In our setting, neurons dynamically alternate between multiple activity states to account for multiple shades of grey in the images. This higherdimensional feature space allows for the encoding and storing of more complex images in the network. The neuronal connections (re-weighted during an initial training phase) encode the desired memories as fixed points of the system. We add Brownian noise fluctuations to the dynamics as a proxy for natural stochastic fluctuations. A single realization of the face recognition task consists in presenting a different, noisy version of a memory

<sup>1</sup>Credit is hereby given to the Massachusetts Institute of Technology and to the Center for Biological and Computational Learning for providing the database of facial images. http://cbcl.mit.edu/software-datasets/heisele/facerecognitiondatabase.html

brain injury. (A) For a noisy input with 80% overlap (i.e., 20% initial noise shown in center), the healthy system (left) achieves an overlap of 90% with the correct face (marked yellow). The chart presents the initial (blue) and final (green) Hamming distances between the current network state and the patterns corresponding to the respective facial images. When the recognition task is performed with an injured system (right), we observe a severe decrease in accuracy (overlap significantly smaller than in the healthy network) and confusion between two facial images. Both are characteristic injurious effects of memory impairments arising from FAS and are quantified with a recognition score in our study. (B) Images of FAS—highlighting their morphometric features—adapted from experimental works used to calibrate our model (Dikranian et al., 2008; Wang et al., 2011).

and assess if the system converges to the corresponding fixed point. Our healthy network shown in **Figure 1** converges to the right fixed point approximately 90% of the time, does not converge for the remaining 10%, and never converges to a wrong fixed point. Note that there are two types of errors (nonconvergence and wrong-convergence). Our recognition score R (see Supplementary Materials) is a metric that encapsulates both types of errors and is better suited to account for deficits in network performance as FAS injuries are introduced in the simulations.

**Figure 2** shows recognition scores (heatmaps) for the most similar triplet of faces as we vary the noise level (from 0 to 30%) and the injury level p (from 0 to 30%). Note that recognition is always strong in the upper-left corner (low noise and low injury) and always weak in the lower-right corner (high noise and high injury). A careful examination of the heatmaps allow us to follow (for each face) how the decline in noise-handling ability evolves into memory confusion.

A key innovation of our model is that a swollen neuron may operate in one of the following regimes: total transmission, filtering, reflection, or blockage. The injury level p denotes the percentage of randomly injured neurons in the network, but their regimes are determined by a FAS distribution (pie-charts). The FAS pie-charts change depending on how long it took for the mice to be sacrificed. **Figure 3A** uses FAS pie-charts inferred experimentally (Wang et al., 2011) for mice sacrificed 12 h after TBI. **Figure 3B** uses pie-charts from mice sacrificed 24 and 48 h after TBI. Finally, **Figure 3C** uses pie-charts from another source (Dikranian et al., 2008), for mice sacrificed 30 min, 5, 16, and 24 h after TBI. This allowed us to quantify the rate at which the system's average recognition score R decrease as a function of injury level p:

$$R(\mathfrak{p}) = A - Be^{\mathfrak{p}} \tag{1}$$

The functional form R(p) is produced from a linear regression over N = 1,080 normalized data points obtained from computational experiments. The values for the constants in Equation (1) are [A = 1, B = 0.24] for pie-charts from Wang et al. and [A = 1.21, B = 0.22] for pie-charts from Dikranian et al. See section Materials and Methods and Supplementary Materials for further details.

#### 2.2. Deterioration Rates for Sparsely-Encoded Uncorrelated Memories

After illustrating our model, injury protocols, and performance metrics for a memory-recall task using highly-correlated faces, we move on to a more abstract setup. Consider now a multi-state Hopfield network with M sparse and uncorrelated memories (see **Figure 4**), where M is now close to the theoretical capacity of the system (M ≈ 0.14N, see Hopfield, 1984). It is obvious that any plausible type of impairment to the individual neurons will ultimately lead to deteriorated recognition scores R as a function of p. In our experiments, for a fixed noise level, the normalized decay rates take the functional form

$$R(\mathfrak{p}) = 1 - Be^{\mathfrak{p}}.\tag{2}$$

Previous studies of injury studies in Hopfield networks where done in a binary setup, i.e., a neuron was either fully functional or fully impaired. See, for instance, Ruppin and Reggia (1995). In our setup, this would correspond to having a pie-chart with 100% of the FAS in the blockage regime. To contrast past approaches with our own, we vary the contributions of the two most prevalent types of FAS mechanisms (filtering and blockage) in our simulations. The distinct deterioration rates are shown in **Figure 4A** and summarized in **Table 1**. The discrepancies are significant and may lead to very different estimates regarding the

FIGURE 2 | Effects of FAS on associative memory for a face recognition task. We measure recognition ability and accuracy for a sample set of three facial images (Weyrauch et al., 2004) over a given range of parameters (injury: 0–30%, initial noise: 0–30%). As the degree of injury increases, the noise handling ability of the system drops severely: The coloring of the heat maps change from yellow (significant recognition) in the upper left (small injury and initial noise) to dark orange and red (confusion) in the lower right corner (high injury and initial noise). We observe confusion and a decline in recall accuracy for every image in the set of samples. Regression over all *N* = 1,080 normalized data points indicates an exponentially declining relation between recognition scores and injury level as shown in the diagram.

FIGURE 3 | (A) Mean recognition scores are displayed in red, standard deviations are shaded in gray. The dashed line indicates the recognition function obtained by linear regression. (B,C) We tested our model using experimental results from TBI studies by Wang et al. (2011) and Dikranian et al. (2008). Our results show a modest, but unstable increase in memory performance for the adult brain during the first 24 h after the injury (B, Wang et al., 2011), that drops again to the initial level during the following day. In the infant's mouse brain, we observe a steadily decreasing performance in the first 24 h (C, Dikranian et al., 2008). The observations match the long recovering times that are commonly observed in patients suffering from TBI.

stage of the disease or its progression. We will discuss this in more depth in the following sections.

#### 3. MATERIALS AND METHODS

#### 3.1. Modeling Dynamical Impacts of FAS

Focal Axonal Swellings (FAS) are an ubiquitous pathological feature to several brain disorders (McArthur et al., 2004; Coleman, 2005; Bayly, 2006b; Tang-Schomer et al., 2010, 2012; Xiong et al., 2013). Their presence is known to distort neuronal spike-train dynamics, but precise electrophysiological recordings of pre- and post-FAS spike train dynamics are still unavailable. The recent theoretical models of Maia and Kutz (Maia and Kutz, 2014a,b; Maia et al., 2015), however, provide important estimates of spike train deterioration due to FAS. Using this model, we can characterize the effect of injury on the firing rate activity through a response function. The state variable of a Hopfield node, S<sup>i</sup> , is equivalent to the firing rate of that node. The collection of all nodes is denoted by the vector **S**. Injuries can then be characterized by the transfer function

$$
\tilde{\mathbf{S}} = F(\mathbf{S}, \boldsymbol{\beta}).\tag{3}
$$

where **<sup>S</sup>**˜ is the effective firing rate (state) after the FAS, **<sup>S</sup>** is the firing rate (state) before the FAS, and β is a parameter vector

indicating one of three injury types applied to individual nodes of the network (Maia and Kutz, 2014a,b). The function F(·) maps the pre-injury to post-injury firing rates using biophysically calibrated statistical distributions of injury in both frequency and size (see Supplementary Materials). If no FAS occurs to a given axon, then its state is unaffected (β0). This occurs with probability 1 − p where p is the probability of injury and what we term injury level, i.e., larger p implies more injury. For those neurons with axonal swellings, the following manifestation of spike train deformations have been observed: unimpaired transmission (β1), filtered firing rates (β2), spike reflection (β3), or blockage (β4). The injury type is dependent upon the geometry of the swelling, with blockage being the most severe injury type. From biophysical data collected on injury statistics (Dikranian et al., 2008; Wang et al., 2011), both in swelling size and frequency, we assign a prescribed percentage of each type of injury (βj) to the network using calibrated simulations of spike propagation dynamics (Maia and Kutz, 2014a,b). For a blockage injury (β4), no signal passes the swelling so the effective firing rate of the neuron goes to zero. Thus, **<sup>S</sup>**˜ <sup>=</sup> 0 which prevents the neuron from adapting to the collective dynamics. Filtering injuries were taken to decrease the firing rate, with higher firing rates having a stronger chance of decreasing due to pile-up effects in the spike train (Maia and Kutz, 2014b). Reflection of spike trains effectively filters the firing rate of an axon by a factor of two so that **<sup>S</sup>**˜ <sup>=</sup> 0.5**S**. This is due to the fact that the reflected spike annihilates an oncoming spike (Maia and Kutz, 2014a). Overall then, the method for producing the filtering function F(·) uses the most advanced experimental findings to date with recent computational studies of spike train propagation through FAS (Maia and Kutz, 2014a,b). See the Materials and Methods section and the Supplementary Materials for details on how

(f: filtering in %, b: blockage in %). The computational experiments use random, uncorrelated memories.

TABLE 1 | Deterioration rates for sparsely-encoded uncorrelated memories.


we compute the statistical distribution of injury types that are parametrized by the parameter β<sup>j</sup> .

#### 3.2. TBI/FAS Data from Adult Rats and Infant Mice

As in Maia and Kutz (2017), we use TBI/FAS data of adult rats from Wang et al. (2011) to calibrate the distributions of FAS following TBI in our simulations. This allowed us to reconstruct distributions of total area per swelling and their functional deficits (blockage, reflection, and filtering). See Supplementary Materials for details. Additionally, we use Dikranian's morphometric analysis of TBI experiments in infant mice (Dikranian et al., 2008). See **Table 2** for the experimental parameters of the two studies.

#### 3.3. Simulation Protocols and Parameters

For the simulations, we considered different calibrations of injury mechanisms (see pie charts in Figures S1,S2): Three in the case of Wang et al. and four in the case of Dikranian et al. The distributions were derived by evaluating swelling parameters at different times of sacrifices after the injury occurred (see Supplementary Materials). For each distribution, we performed multiple rounds of simulations, each running over a wide range of injury and noise levels. Recognition scores resulting from the

TABLE 2 | Experimental parameters in animal studies by Wang et al. (2011) and Dikranian et al. (2008).


TABLE 3 | Experimental parameters used in simulations.


analysis of all simulation results were used to fit decay rates. The analysis was performed for two recognition tasks, using (i) highly-correlated facial images and (ii) sparse, random patterns. The running times of the simulations strongly depend on the type and amount of injury introduced into the network. For larger numbers of neurons and significant levels of injury, the simulations are computationally very expensive. The training itself depends on the number of neurons as O(N 2 ). Further details on the model and its implementation can be found in the Supplemental Material.

The choice of simulation parameters (see **Table 3**) was guided by balancing between computational expenses and covering a reasonably large parameter range. All code is publicly available on GitHub for possible future studies.

### 4. DISCUSSION

Overwhelming experimental evidence suggests that FAS is the hallmark manifestation of injuries on neuronal networks. Moreover, there now exists a wealth of experiments characterizing the statistical distribution of FAS as a function of injury level, including size and frequency of swellings (Coleman, 2005). Statistics can even be collected for specific neurodegenerative diseases such as Alzheimer's (Krstic and Knuesel, 2012), Parkinson (Galvin et al., 1999), or Multiple Sclerosis (Hauser et al., 2006), where the swellings occur as a consequence of complicated biophysical and biochemical deterioration of neurons. These neurodegenerative diseases affect a large portion of adults, with Alzheimer's disease alone estimated to be the third leading cause of death, just behind heart disease and cancer, for older people. Likewise, TBI is responsible for millions of hospitalizations worldwide every year (especially among contact-sport practitioners) and is the leading cause of death among youngsters (Faul et al., 2010; Jorge et al., 2012; Xiong et al., 2013).

In this work, we consider impairments caused by FAS in a Hopfield-like computational model for associate memory. In the face-recognition task, the presence of a significant number of blocked neurons in the network blurs the reconstructed concepts (faces) and decreases the accuracy of the recalled information. Both filtering and reflection regimes lead to confusion of correct states with their neighboring ones (in a Hamming distance sense). This decreases the ability of the network to perform denoising tasks. In most cases, all three impaired regimes (filtering, reflection, and blockage) occur simultaneously in the network (Maia and Kutz, 2014a). An interesting consequence is that injured networks often produce erroneous associations of memories. Specifically, it confuses the concepts. Based on simulations of an uninjured network with noise, it could easily be conjectured that the fixed points associated with a given memory would disappear, with the system's dynamics simply preventing it from converging to the correct pattern. Instead, we observe that the noise fluctuations cause the dynamics to converge to an erroneous fixed point. For modest amounts of initial noise (15–30%), the affected system confuses the stored patterns and looses the ability to separate them properly. Confusion of concepts is especially pronounced when blockage and filtering account for the majority of FAS effects.

We calibrate our FAS distributions (pie-charts) from two experimental sources (Dikranian et al., 2008; Wang et al., 2011). We believe that the distributions from Dikranian et al. (2008) are more relevant since they damage the cingulum—recognized as fundamental to memory association—instead of the optic nerve. In addition to a TBI protocol better linked physiologically to memory impairments, they provided a distribution plot for the diameters of the spheroids, which allowed us to generate FAS in a much more direct way. See the Supplementary Materials for more details.

Our second associative-memory task was in a more abstract setup, with memories being random and uncorrelated. Previous studies investigated injuries in this kind of network, but in a binary way, i.e., when neurons are either fully injured or fully functional. Binary setups, however, cannot capture the types of frequency-dependent errors demonstrated in current FAS research (Maia and Kutz, 2014b). As a consequence, they cannot handle the nuanced types of deficits observed in TBI and neurodegenerative diseases, and clinicians would likely underestimate p from the R scores of their patients using the decay coefficient B derived from these sources. It is unlikely that estimates obtained from binary injury protocols where the error produced by a single cell is always maximal could be applicable to concussions (milder and more common forms of TBI) or early dementia, where successful diagnostics are most needed. In fact, in all seven of our derived FAS pie-charts, injured neurons in blocking regime are <50%, whereas a significant percentage of impaired neurons retain some activity level (reflection, filtering). Distinguishing FAS distributions for different stages of neurodegenerative diseases is especially important in the context of diagnostics, where the focus lies on detecting early signs of a progressing condition.

# 5. CONCLUSION

Focal Axonal Swellings (FAS) are a hallmark feature of TBI and neurodegenerative diseases such as Alzheimer's, Parkinson's and Multiple Sclerosis. Our study characterizes a neural network's ability to handle noise and perform recognition tasks at different levels of injury. Our model reproduces symptoms commonly observed in patients suffering from brain disorders (Knutson, 2012) arising from the above discussed causes. In the facerecognition task, the injured network loses accuracy in recalling facial images and confuses faces with similar features. From a Hopfield network viewpoint, less accuracy means that fixed points encoding memories loose their stability. At higher injury levels, the network's dynamic is driven away from the correct face encoding, settling closer to a fixed point associated with another facial pattern. Instability of the fixed points and emergence of wrong attractors results from a large percentage of neurons manifesting FAS, leading ultimately to confusion between previously stored images. Performance decline was formalized through estimating deterioration rates. for the recognition score. We tested the influence of different distributions of FAS mechanisms on deterioration rates for a set of random, uncorrelated memories. Our results show a significant discrepancy between our proposed multi-mechanism FAS model and compared to the previously studied binary lesion models.

We calibrate FAS parameters in our simulations with experimental TBI data; Dikranian et al. (2008) studied FAS in the cortex of infant mice, whereas Wang et al. (2011) investigated swellings in the optic nerve of adult rats. Although Dikranian et al. (2008) provides detailed morphometric data for injured axons in the infant mice, their functional assessments might be still a distant proxy for memory development in human children (McAllister, 2004). TBI experiments in adult rats (Xiong et al., 2013) report functional impairments more analogous to human patients—like deficits in tasks involving context memory (Schacter, 1987), conditional associative learning (Petrides, 1985), planning (Shallice and Evans, 1978), and other cognitive tasks (McDowell et al, 1997). Adult rats, however, may exhibit memory deficits after mild TBI even without many signs of

#### REFERENCES


axonal injury (Lyeth, 1990). In fact, recent studies demonstrate that catecholamines play a central role in the neurochemical activation and regulation of working memory (McAllister, 2004) and such effects were not incorporated in the model. Thus, axonal structural damage may be sufficient but not necessary for the production of neurological and cognitive symptoms associated with TBI. This will be considered in future studies.

There is much room for improvement in our injury model, especially given the variety of effects of different types of brain disorders. The framework introduced in the present article is limited to associate memory. While other forms of memory impairment are also known to occur as a result of neurodegenerative diseases and traumatic brain injury, incorporating these into the model was beyond the scope of this study. Furthermore, we only consider the three basic injurious FAS mechanisms introduced in the Maia and Kutz model and leave out additional physiological effects that might affect the severity of memory impairments. These include the squeezing of neighboring axons when swellings occur in densely connected regions.

Although FAS are universal pathological features, our results should be regarded as a first step towards integrating them into functional neuronal networks and linking cellular impairments to observable psychophysical deficits.

# AUTHOR CONTRIBUTIONS

MW, PM, and JK designed; MW and PM: conducted the study; and MW, PM, and JK: wrote the manuscript.

#### ACKNOWLEDGMENTS

MW was supported by a scholarship of the Konrad Adenauer foundation.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins. 2017.00623/full#supplementary-material


with mild traumatic brain injury. Amer. J. Psychiatry 169, 1284–1291. doi: 10.1176/appi.ajp.2012.12050600


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer PJT and handling Editor declared their shared affiliation.

Copyright © 2017 Weber, Maia and Kutz. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Data-Driven sequence of changes to anatomical Brain connectivity in sporadic alzheimer's Disease

*Neil P. Oxtoby1 \*, Sara Garbarino1 , Nicholas C. Firth1,2, Jason D. Warren2 , Jonathan M. Schott2 , Daniel C. Alexander1 and For the Alzheimer's Disease Neuroimaging Initiative*

*1Progression of Neurodegenerative Disease Group (POND), Centre for Medical Image Computing, Department of Computer Science, University College London, London, United Kingdom, 2Dementia Research Centre, Department of Neurodegenerative Disease, UCL Institute of Neurology, University College London, London, United Kingdom*

Model-based investigations of transneuronal spreading mechanisms in neurodegenerative diseases relate the pattern of pathology severity to the brain's connectivity matrix, which reveals information about how pathology propagates through the connectivity network. Such network models typically use networks based on functional or structural connectivity in young and healthy individuals, and only end-stage patterns of pathology, thereby ignoring/excluding the effects of normal aging and disease progression. Here, we examine the sequence of changes in the elderly brain's anatomical connectivity over the course of a neurodegenerative disease. We do this in a data-driven manner that is not dependent upon clinical disease stage, by using event-based disease progression modeling. Using data from the Alzheimer's Disease Neuroimaging Initiative dataset, we sequence the progressive decline of anatomical connectivity, as quantified by graph-theory metrics, in the Alzheimer's disease brain. Ours is the first single model to contribute to understanding all three of the nature, the location, and the sequence of changes to anatomical connectivity in the human brain due to Alzheimer's disease. Our experimental results reveal new insights into Alzheimer's disease: that degeneration of anatomical connectivity in the brain may be a viable, even early, biomarker and should be considered when studying such neurodegenerative diseases.

Keywords: brain connectivity analysis, data-driven, Alzheimer's disease, disease progression modeling, graph theory analysis, computational model

# 1. INTRODUCTION

There is a growing body of literature using advanced computational and statistical modeling to understand progressive disorders. The majority of these has been applied to the most common neurodegenerative disorder of the brain, Alzheimer's disease (AD)—see Ref. (1), for a review of the field. Discriminative models range from simple diagnosis using supervised machine-learning classifiers (2–4) to unsupervised clustering of disease subtypes (5–7). Such discriminative models typically do not directly include a notion of disease progression. By contrast, generative data-driven models (8–11) estimate disease progression signatures that can aid disease understanding from start to end. Such an understanding is vital for the early identification of individuals who are likely to be responsive to a particular therapy, and also for identifying the earliest pathological changes for designing such therapies. The protracted duration of AD, and especially the decades-long presymptomatic

#### *Edited by:*

*Yasser Iturria Medina, Montreal Neurological Institute and Hospital, Canada*

#### *Reviewed by:*

*Gianfranco Spalletta, Fondazione Santa Lucia (IRCCS), Italy Elis Eleutherio, Federal University of the State of Rio de Janeiro, Brazil*

#### *\*Correspondence:*

*Neil P. Oxtoby n.oxtoby@ucl.ac.uk*

#### *Specialty section:*

*This article was submitted to Neurodegeneration, a section of the journal Frontiers in Neurology*

*Received: 14 August 2017 Accepted: 17 October 2017 Published: 07 November 2017*

#### *Citation:*

*Oxtoby NP, Garbarino S, Firth NC, Warren JD, Schott JM and Alexander DC (2017) Data-Driven Sequence of Changes to Anatomical Brain Connectivity in Sporadic Alzheimer's Disease. Front. Neurol. 8:580. doi: 10.3389/fneur.2017.00580*

**104**

phase, makes this a challenging disorder to study in the general population. Dominantly inherited variants of AD, e.g., Ref. (12)., can be used to alleviate this problem to some extent, but the open question of whether dominantly inherited AD is a suitable model for sporadic late-onset AD remains unanswered.

Alzheimer's disease has been widely described as a disconnection syndrome (13–16). This has contributed to the motivation behind the relatively recent emergence of network models of neurodegeneration. The motivation is to understand mechanisms of disease propagation by relating models of brain networks to observed pathology. The networks are abstract representations of "connections" between brain regions and can be constructed from (1) functional correlations during rest or activity (usually estimated using blood-oxygen-level-dependent contrast MRI (17)), (2) gray-matter covariance (estimated from structural MRI), or (3) anatomical connections (estimated using neuronal tractography from diffusion-weighted MRI).

Network spreading models attempt to explain neurodegeneration in terms of the transneuronal spread of prion-like (18) pathogens such as abnormal proteins. Anatomical connectivity networks are a natural choice for such models as they estimate physical connections between brain regions, rather than the correlations estimated in functional and gray-matter structural covariance networks. Previous work in network spreading models of neurodegenerative diseases has considered how well these models can predict end-stage disease by: correlating patterns of healthy intrinsic (functional) connectivity and gray-matter volume (14); using healthy functional connectivity to compare network-based mechanistic hypotheses of AD progression (19); and using healthy anatomical connectivity to predict atrophy (20), amyloid load (21), or metabolism (22). All of these network spreading models used *static* connectivity of healthy, *young* individuals to build networks with which to predict end-stage AD pathology. This ignores the effects of aging and disease progression on the network substrate being used to predict pathology.

In contrast to previous network spreading models, we consider *elderly* connectivity networks, and we explicitly model the course of disease progression. This enables us to investigate how the brain's anatomical connectivity *changes* with AD progression. We do this in a novel manner by analyzing regional network (graph-theoretic (23, 24)) measures of brain connectivity in groups of healthy and diseased individuals, within the context of an event-based model (8, 10) of disease progression. The model produces a data-driven, fine-grained signature of the sequence of disease-related changes in anatomical connectivity of the human brain, including uncertainty in the sequence. Our innovations beyond earlier event-based models (8, 10) include analyzing biomarkers of brain connectivity and employing a new nonparametric mixture modeling technique (25) for estimating biomarker abnormality that is built upon kernel density estimation (26, 27).

#### 2. MATERIALS AND METHODS

Our analysis can be summarized as three steps. First, we used probabilistic anatomically constrained tractography to construct individual whole-brain connectomes for imaging data from 168 participants from the Alzheimer's Disease Neuroimaging Initiative (ADNI) study, across four age-matched and educationmatched clinical categories (4 × 42 per category, the maximum available in one of the categories). The image analysis pipeline is visualized in the schematic of **Figure 1**. Second, from each connectome, we computed various local network measures (graph-theory metrics) representing the topology of anatomical connectivity in regions of each participant's brain. Finally, we estimated the ordered sequence in which these measures become abnormal using an event-based model (8, 10).

#### 2.1. Data

Data used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database.1 The ADNI was launched in 2003 as a public–private partnership, led by Principal Investigator Michael W. Weiner, M.D. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD).

In February 2017, we remotely accessed the Laboratory of NeuroImaging's Image Data Archive at the University of Southern California and searched for suitable participants to include in our anatomical connectome cohort: ADNI participants whose brains were imaged with both structural magnetic resonance imaging (MRI) and diffusion-weighted imaging (DWI) at a single study visit. We sought age-matched groups across the four diagnoses of Cognitively Normal (CN), Early Mild Cognitive Impairment (EMCI), Late Mild Cognitive Impairment (LMCI), and probable AD (AD). This resulted in 42 participants per group. We

1http://adni.loni.usc.edu.

preprocessing → DWI normalization → rigid registration of T1 MRI to DWI → T1 parcelation → anatomically constrained tractography → connectomes based on density of WM neuronal connections between GM regions of interest. Abbreviations: DWI, diffusion-weighted image; MRI, magnetic resonance image; ACT, anatomically constrained tractography; WM/GM, white/gray matter.

downloaded the unprocessed structural MRI (3.0T, T1-weighted, non-accelerated IR-SPGR; GE Medical Systems) and DWI (Axial diffusion tensor imaging; GE Medical Systems) for the first available suitable visit of these 168 participants from the ADNI2 phase of ADNI (including: 5 rollovers from the first ADNI phase, ADNI1; and a single rollover from the ADNI Grand Opportunity phase, ADNIGO). Herein, we refer to our anatomical connectome cohort as *The168*. We also downloaded associated demographics data and metadata in CSV format. Demographics for *The168* are summarized in **Table 1**, and for all ADNI2 (at ADNI2 baseline) in **Table 2**.

Our disease progression modeling requires both a control and a patient group, which we defined as amyloid-negative CN participants and amyloid-positive AD participants, respectively (see bold figures in **Table 1**). The threshold for amyloid positivity was chosen as an amyloid PET (Florbetapir 18F-AV-45, hereafter AV45) cut-point from the literature (28): AV45 Standardized Uptake Value Ratio *SUVR* ≥ 1.10, which was based on the upper 95% confidence interval for young healthy subjects. Using this criteria, in *The168*, we identified 26 controls and 38 patients out of a possible maximum of 42 each. For the wider ADNI2 cohort (including rollovers from ADNI1 and ADNIGO), we found 106 controls and 122 patients—see the bold figures in **Table 2**.

#### 2.2. Connectomics

Structural connectomes were generated using tools provided in the MRtrix3 software package,2 customized to work with the Geodesic Information Flows algorithm (29) for segmentation and parcelation. The pipeline included DWI denoising (30),

2http://mrtrix.org.

Table 1 | Demographics by diagnosis for our anatomical connectome cohort, *The168*: 168 ADNI2 participants included in this study.


*One-way ANOVA implies that the groups do not differ significantly by age (p* > *0.97), nor education (p* > *0.45). Amyloid positivity was set at cortical mean AV45 SUVR* ≥ *1.10 (28), prior to adjustment for covariates. Numbers in bold font indicate controls and patients (see text for details) used in our disease progression modeling.*

Table 2 | Demographics by diagnosis for all ADNI2 participants at initial visit.


*As in Table 1, amyloid positivity was set at cortical mean AV45 SUVR* ≥ *1.10 (28), prior to adjustment for covariates. Numbers in bold font indicate controls and patients (see text for details) used in our disease progression modeling.*

preprocessing (31, 32), and bias field correction (33); inter-modal registration (34); T1 tissue segmentation (29); spherical deconvolution (35, 36); probabilistic tractography (37) utilizing anatomically constrained tractography (38) and dynamic seeding (39); spherical deconvolution informed filtering of tractograms (SIFT) (40); T1 parcelation (29); and robust structural connectome construction (41). We used the dwiintensitynorm script in MRtrix3 during DWI preprocessing, but found it necessary to modify the subsequent usage of the population\_template script such that DWI masks were not used to create the template (the DWI were pre-masked). Note that we use SIFT to produce biologically plausible tractograms where the streamline density in each voxel matches the fiber orientation distributions estimated from the DWI, instead of simply thresholding the number of fibers connecting two regions.

Our anatomical connectome for each participant is a weighted adjacency matrix that includes only inter-node connections across 130 regions of interest consisting of cortical and subcortical gray-matter regions (including striatal), plus the cerebellum and brain stem. Weights, or connection strengths, are normalized to [0, 1], and so represent within-participant anatomical connection density. The image analysis pipeline is visualized in **Figure 1**. The 130 regions of interest are a subset of those in the labeling protocol used by Geodesic Information Flows (29, 42), which is a modified version of the Desikan–Killiany–Tourville protocol (43).

#### 2.3. Anatomical Connectivity Metrics

For each participant and region of interest, we calculated 12 brain connectivity metrics from the anatomical connectomes, using the Brain Connectivity Toolbox (24) in MATLAB, after appropriate normalization using the weight\_conversion function. The local network metrics fall into the following broad categories:


For each region of interest, the medians of the controls and patients distributions for each network measure were statistically compared using a Mann–Whitney–Wilcoxon rank-sum test. Only measures where *p* < 0.05/12 (Bonferroni-corrected within region) were retained for further analysis within the event-based model of disease progression—we refer to such biomarkers as having "disease signal."

Note that local efficiency and local clustering coefficient return similar information to each other, as do eigenvector centrality and PageRank centrality, and so one may ask whether we are including redundant information in our models. We argue that we are not, as seen in our results where: (1) local efficiency contained disease signal, whereas local clustering coefficient did not and (2) the centrality measures appear in different positions within the model sequences.

# 2.4. Event-Based Model of Disease Progression

The event-based model (EBM) (8, 10) is a data-driven approach for probabilistically sequencing a cross-sectional set of scalar measurements ("biomarkers") in the order in which they become observably abnormal. In this context, an "event" constitutes a biomarker appearing more abnormal/diseased than normal/healthy. The EBM is able to estimate an average sequence of disease progression events from cross-sectional data because the proportion of abnormal measurements within a cohort will decline in concert with the average ordering. That is, the biomarker that changes earliest (the first disease event) will contain the highest proportion of abnormal measurements (from affected and presymptomatic individuals), and so on. The EBM fuses multiple biomarker measurements across individuals, with the simplest versions assuming a single disease progression sequence for all individuals, as done here. Determining biomarker abnormality in a data-driven manner necessitates mixture modeling within biomarkers, for which we use a new method (25), described below in Section 2.4.1. For convenience, we normalize all biomarkers to "c-scores" (standardized z-scores relative to controls) that increase with abnormality. Otherwise, we used the same fitting procedures as in Ref. (10). Crossvalidation of our EBMs was estimated by refitting the sequence (but not the event measures) to 100 separate bootstrap samples from the data.

We planned to build four EBMs of network connectivity changes, corresponding to the three broad categories of connectivity in Section 2.3, plus all biomarkers together:


EBM0 acts as a reference point and for investigating consistency with previous EBMs of AD (10, 45) and includes only non-network biomarkers: average cortical level of amyloid (from AV45 PET) and hypometabolism (from fludeoxyglucose (FDG) PET), test score on the Mini-Mental State Examination (MMSE) (46), and structural MRI volumes of the hippocampus, entorhinal area, ventricles, and whole brain. Biomarkers were adjusted for healthy age, education, and gender using regression (residuals method, controls only). Brain volumes were also adjusted for intracranial volume. We did not fit EBM3 because there was no disease signal (see Section 2.3) in regional clustering coefficients nor eccentricity.

#### 2.4.1. Biomarker Event Measures

The probability of a biomarker event is fundamental for sequencing the biomarker events into a pathological cascade of disease progression. Since not all patients will have experienced later events, and indeed some controls will have already experienced the earliest events in the cascade, it is necessary to fit a mixture model in order to discover the event probability. Previous EBM analyses (8, 10, 45, 47–49) used mixtures of parametric probability distributions such as Gaussian and uniform. Here, we use a new method (25) that fits a mixture of nonparametric kernel density estimate (KDE) distributions.

# 3. RESULTS

#### 3.1. Global Connectivity

**Figure 2** shows that global network connectivity in health and disease did not differ significantly in our cohort, *The168*. That is, we found no significant group-level difference between the 26 controls (AV45-positive CN) and 38 patients (AV45-negative AD) across four global brain network metrics: density (connectedness), efficiency (segregation), transitivity (segregation), and assortativity (resilience). The null hypothesis in each Mann–Whitney–Wilcoxon rank-sum test was accepted using a Bonferroni-corrected significance level of *p* = 0.05/4 = 0.0125: density *p* = 0.0147; transitivity *p* = 1; efficiency *p* = 0.145; assortativity *p* = 0.363.

# 3.2. Event-Based Models

#### 3.2.1. Biomarkers

**Figure 3** is a visualization of disease signal in the full set of 38 biomarkers included in the EBMs. The vertical axis is the standardized "c-score" for each biomarker along the horizontal axis: c-score is biomarker value standardized to controls (c.f., z-scores). Group-average lines are shown with individual data points as green crosses for controls and red dots for patients.

#### 3.2.2. Disease Progression Sequences

The EBM estimates a data-driven probabilistic sequence of biomarker abnormality, which we visualize as plots of grayscale positional variance (horizontal axis) around the maximumlikelihood ordering (vertical axis). The strongest possible ranking of biomarker abnormality would appear as a black diagonal.

Our first experiment was to build "EBM00"—an EBM including only non-network biomarkers, and using available ADNI2 baseline data (summarized in **Table 2**). The results are shown in **Figure 4**, with positional variance estimated from the MCMC fitting procedure (8) shown in the left of the figure, and from bootstrapping shown in the right of the figure (cross-validation). The EBM00 sequence is consistent with current understanding of AD progression (and previous EBMs (10)): early amyloidosis and hippocampal volume loss, followed by cognitive decline, then hypometabolism and broader neurodegeneration.

In our next experiment, we built "EBM0"—the same biomarkers as EBM00, but using only data from our anatomical connectome cohort, *The168* (**Table 1**). The results are shown in **Figure 5**, with the positional variance diagram from fitting on the left, and from bootstrapping on the right. Like EBM00, the EBM0 sequence is consistent with current understanding of AD progression, with only one exception: the apparent late appearance of hypometabolism (FDG). We attribute this to the relatively small number of probable AD patients in *The168* with abnormally

Figure 3 | Standardized "c-scores" for all biomarkers included in our analyses. Non-network markers are to the left, denoted EBM0. Network metrics are to the right, delineated by EBM1 (hubs, i.e., basic measures of centrality) and EBM2 (importance, i.e., advanced measures of centrality and shortest paths). EBM3 was not generated since no segregation/integrationbased network measures contained disease signal. EBM4 includes all biomarkers.

high metabolism, compared to the full ADNI2 cohort (Figure S1 in Supplementary Material). We also note that the positional variance is larger in EBM0 than in EBM00, probably due to the lower numbers of individuals.

Next, we report results from experiments on brain connectivity changes in AD. **Figure 6** shows "EBM1," which includes network biomarkers that measure hubs in the brain network through node degree or strength. **Figure 7** shows "EBM2," which includes biomarkers of network importance as measured by centrality and shortest path metrics. **Figure 8** is "EBM4" which includes both hubs and centrality. ("EBM3" does not exist because no biomarkers of segregation/integration passed the test for disease signal—see Section 2.3.). The network biomarker event measures are presented in Figure S2 in the Supplementary Material.

From **Figure 6** we can infer the ordering the network biomarker event measures (from mixture modeling) are presented in Figure S2 in the Supplementary Material. In which certain anatomical network hubs of the brain deteriorate. We note early involvement of regions in the temporal lobe: the left transverse temporal gyrus (TTG; a.k.a. Heschl's gyrus; auditory cortex) and the left temporal pole (TMP; anterior of the temporal lobe); and a region in the frontal lobe: the right lateral orbital gyrus (LOrG). The model suggests that the later hubs to deteriorate include the middle occipital gyrus (MOG), subcallosal area (SCA), left occipital fusiform gyrus (OFuG), and right middle cingulate gyrus (MCgG). Many, but not all, of these regions are involved in the default mode network (DMN), see, e.g., Ref. (50), and references within.

From **Figure 7** we infer the sequential deterioration of regional importance in the anatomical network of the brain, as measured by local centrality and efficiency. The model suggests that centrality declines first in memory-related parts of the DMN in the temporal and prefrontal lobes: left temporal pole

Figure 4 | EBM for selected standard (non-network) biomarkers, built on available ADNI2 baseline data. Left panel: left-right positional variance of the maximumlikelihood sequence. Right panel: cross-validation of the sequence from bootstrapping.

(TMP), left triangular part of the inferior frontal gyrus (TrIFG), bilateral hippocampus; and last in regions from the occipital (and anterior limbic) lobe(s): right lingual gyrus (LiG; visual cortex), left middle occipital gyrus (MOG), right middle cingulate gyrus (MCgG), and right angular gyrus (AnG).

Finally, **Figure 8** ("EBM4") contains all biomarkers of anatomical connectivity included in this study. The ordering of abnormality in hubs and centrality biomarkers is consistent with EBM1 and EBM2, respectively. **Figure 8** enables comparison of the relative sequential decline among hubs and regional centrality in the AD brain's anatomical connectivity network. We observe that the clearest and earliest decline in anatomical connectivity in the AD brain is both hubs and centrality in the left temporal lobe, all of which are DMN regions known to have memoryrelated storage and functions. Memory dysfunction is the clinical phenotype of typical AD dementia.

# 4. DISCUSSION

#### 4.1. Findings in Context

Brain network hubs are thought to experience increased susceptibility to AD pathology due to higher activity and

Figure 6 | EBM for hubs in the brain's anatomical network, along with selected standard (non-network) biomarkers, built on data from our anatomical connectome cohort, *The168*. Left panel: left-right positional variance of the maximum-likelihood sequence. Right panel: cross-validation of the sequence from bootstrapping. Abbreviations: DegreeZ, degree z-score; L, left; R, right; TTG, transtemporal gyrus; TMP, temporal pole; LOrG, lateral orbital gyrus; TrIFG, triangular inferior frontal gyrus; MCgG, middle cingulate gyrus; IOG, inferior occipital gyrus; MOG, middle occipital gyrus; SCA, subcallosal area; OFuG, occipital fusiform gyrus.

metabolism, such as regions within the DMN, e.g., Ref. (50), and references within. Our experimental results are consistent with this idea, suggesting that the earliest abnormality in AD occurs in hubs (**Figure 6**) and other centrally important regions (**Figure 7**) in the left temporal lobe, and bilaterally in the hippocampus—see **Figure 8**. Specific left temporal lobe regions include the transverse temporal gyrus (part of the auditory cortex), the temporal pole (which may be involved in social and emotional cognition, e.g., because it has been shown to be affected in frontotemporal dementia (51)), and the triangular part of the inferior frontal gyrus (which may be involved in semantic memory (52)).

Our finding of early connectivity changes involving the auditory cortex is potentially of particular clinical relevance to Alzheimer's disease: hearing impairment has recently emerged as a major risk factor for cognitive decline and a focus of intense epidemiological interest (53), while functional alterations in central auditory processing have been identified in both presymptomatic and established Alzheimer's disease (54, 55). Our findings suggest how such alterations may fit within the pathogenic cascade of Alzheimer's evolution.

Our results also suggest that this deterioration of anatomical connectivity in the brain network may be detectable prior to bulk estimates of amyloidosis, such as the mean SUVR across the cortex in AV45 PET. However, prior to claiming that networkbased biomarkers may form sensitive early biomarkers of AD, we would seek confirmation using larger anatomical connectivity cohorts (see future work in Section 4.3).

Our experiments also identified that the anatomical connectivity of some DMN regions is not affected until later in the pathological cascade of AD. These included hubs in the angular gyrus, and hubs and centrality in the medial cingulate gyrus—see **Figure 8**.

We note that the ordering of brain regions involved in our sequence of anatomical connectivity changes to the Alzheimer's disease brain (**Figure 8**) does not follow the ordering of regions in the neuropathological sequence identified by Braak and Braak (56). This suggests that changes in anatomical connectivity may not occur in sync with deposition of abnormal protein in brain tissue. Indeed, our results suggest that connectivity changes may appear before bulk deposition of amyloid across the cortex is detectable *in vivo*.

The brain's anatomical connectivity network is altered by AD. These disease-related connectivity changes may be due to neurodegeneration in gray matter, alterations in deep white matter, or alterations in superficial white matter (57) located near the gray-matter/white-matter boundary. We designed our pipeline to be sensitive to these disease-related alterations by employing anatomically constrained tractography (38) and the SIFT method (40), which should improve the biological accuracy of the tractography-based connectivity estimates.

We have studied the late onset, typical variant of Alzheimer's disease. It will be of interest to apply our technique to explore genetic factors that might influence the progression of Alzheimer's disease; and to young onset Alzheimer's disease where different phenotypic presentations are more commonly seen.

#### 4.2. Novelty of This Work

Our approach differs from previous network spreading models in two very important ways. We investigated *dynamic* pathological disruption of the *elderly* brain's anatomical connectivity network as a function of AD *progression*. Previous network models of AD (14, 19–22) used *static* connectivity patterns estimated from young and healthy individuals, then employed only end-stage patterns of pathology, such as atrophy, to study network spreading mechanisms. Thus, the previous approaches ignore both the effects of normal aging and the effect of disease progression on the network substrate being used to predict pathology.

Our approach provides utility beyond contributing to our understanding of the AD pathological cascade. These include the ability to estimate a longitudinal disease progression signature that is useful for patient staging (10), and for identification of candidate biomarkers for early diagnosis.

In this respect, we believe that our work may be the first data-driven quantification of the widely held idea that AD progression reflects specific anatomical network disruption, not just neurodegeneration.

# 4.3. Future Work

Future work will include seeking a separate, larger anatomical connectivity cohort upon which to validate our results. While our own cross-validation experiments reported above provide some level of confidence in our conclusions, external validation on a larger cohort would be all the more convincing. Validation on animal models of AD is a related idea that might be worth pursuing.

Our analyses could be applied to other neurodegenerative diseases. We will actively seek suitable cohorts from other diseases upon which to examine pathological changes in the connectivity of the human brain. This would subsequently extend quite naturally to include application of the models to differential diagnosis, and to the related problem of within-disease subtyping (7, 58).

On the highest level, our work involved comparison of graphs (59). There are subtle challenges to comparing graphs and metrics derived from them, even when the number of nodes is constant and the average number of connections is the same (60), as was the case in our study. In the future, we will consider normalizing our connectivity metrics with those derived from random graphs or by looking at distances between graphs (60).

## 4.4. Conclusion

We have sequenced the progressive deterioration of anatomical connectivity in the Alzheimer's disease brain. We believe that this is the first attempt to do so in a data-driven manner that incorporates the effects of both aging and disease progression. Aging is accounted for by using healthy elderly brains as controls. Disease progression is estimated using event-based modeling.

Our experimental results reveal new insights into Alzheimer's disease progression: that degeneration of anatomical connectivity in the brain may be a viable, even early, biomarker and should be considered when studying such neurodegenerative diseases.

# AUTHOR NOTE

Data used in preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (http://adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at: http://adni.loni.usc.edu/wp-content/uploads/ how\_to\_apply/ADNI\_Acknowledgement\_List.pdf.

# AUTHOR CONTRIBUTIONS

NO, SG, and DA conceived the study and designed the experiments, in consultation with JS. NO performed the image analysis that produced the connectomes and regional volumes and wrote the manuscript. NO and SG collated the data and performed the event-based modeling analyses with input from NF (in particular the mixture modeling for event measures), JS, and DA. All authors contributed to interpretation of the results and the final version of the manuscript.

# ACKNOWLEDGMENTS

The authors gratefully acknowledge the High Performance Computing team in UCL's Department of Computer Science for providing and maintaining invaluable computational resources.

# FUNDING

NO, SG, JS, and DA are supported by a project that has received funding from the *European Union's Horizon 2020 research and innovation programme* under grant agreement number 666992. NO acknowledges financial support from The Michael J. Fox Foundation for Parkinson's Research, the Alzheimer's Association, Alzheimer's Research UK, and the Weston Brain Institute through "NetMON: Network Models Of Neurodegeneration" (BAND-15-368107, 11042). JW received grant support from the Alzheimer's Society and the NIHR UCLH Biomedical Research Centre. 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 Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer's Association; Alzheimer's Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai, Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche, Ltd. and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer, Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (http://www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/article/10.3389/fneur.2017.00580/ full#supplementary-material.

# REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Oxtoby, Garbarino, Firth, Warren, Schott and Alexander. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*