The role of clearance in neurodegenerative diseases

Alzheimer’s disease, the most common form of dementia, is a systemic neurological disorder associated with the formation of toxic, pathological aggregates of proteins within the brain that lead to severe cognitive decline, and eventually, death. In normal physiological conditions, the brain rids itself of toxic proteins using various clearance mechanisms. The efficacy of brain clearance can be adversely affected by the presence of toxic proteins and is also known to decline with age. Motivated by recent findings, such as the connection between brain cerebrospinal fluid clearance and sleep, we propose a mathematical model coupling the progression of toxic proteins over the brain’s structural network and protein clearance. The model is used to study the interplay between clearance in the brain, toxic seeding, brain network connectivity, aging, and progression in neurodegenerative diseases such as Alzheimer’s disease. Our findings provide a theoretical framework for the growing body of medical research showing that clearance plays an important role in the etiology, progression and treatment of Alzheimer’s disease.


Introduction
Neurodegenerative diseases such as Alzheimer's disease (AD) are progressive disorders associated with a gradual loss of cognitive faculties and, ultimately, a state of dementia. These diseases are characterized by the presence of proteins, present under healthy conditions, but misfolded in the pathological state. Misfolded protein isoforms can be toxic to cells in the brain and form aggregates that are thought to lead to cell death and cognitive decline. The prion-like hypothesis of neurodegenerative diseases [36,48,35,12,22,43] posits that toxic proteins can transmit their misfolded state to otherwise healthy proteins and induce, autocatalytically, the formation of aggregates. Further, misfolded proteins appear to propagate mainly along neuronal pathways [14,25,46,78]. In Alzheimer's disease, the key evidence for the prion-like hypothesis is the presence of Amyloid-β (Aβ) and tau proteins (τP), whereas Parkinson's disease are associated with α-synuclein, SOD1 and TDP-43 are hallmarks of amyotropic lateral sclerosis, as are τP and TDP-43 for frontotemporal lobar degeneration [65].
Low levels of toxic proteins in a brain region do not lead directly to pathological aggregations as various clearance mechanisms continuously remove proteins and other solutes from the brain tissue [65,24]. These mechanisms can be classified as cellular (degradation) clearance, bloodbrain barrier (BBB) mediated clearance, and/or cerebrospinal fluid (CSF) mediated clearance [66]. Mounting experimental or clinical evidence suggest that lifestyle factors such as sleep [79,40,59,18] or normal aging [38] may also alter clearance systems. Since different neurodegenerative diseases are typically characterized by different misfolded protein pathologies, the interaction between a clearance system and a particular protein is often protein-specific and many basic processes remain poorly understood. For instance, while BBB clearance is well established, a complete understanding of CSF-mediated clearance is still elusive for most proteins relevant to neurodegenerative disorders [3,42,6,2,66,56]. Misfolded proteins may interfere with normal brain clearance, and the failure of one or more of the brain's clearance systems could play a role in the pathogenesis of neurodegenerative diseases [65,66,47,28].
Mathematical modelling is now emerging as a promising technology for creating improved understanding of mechanisms underlying neurodegenerative disease [45,41]. In particular, network diffusion models posed on a brain structural connectome [54,55,49] illustrate the potential of graph-based mathematical methods in studying neurodegeneration. These early models have later been generalized to include mechanisms such as aggregation, fragmentation, and clearance observed in-vivo [20,21,67,68]. However, current models fail to include a feedback mechanism that takes into account how clearance pathways are themselves affected by the toxicity of misfolded proteins [5,10,9].
Here, we propose a network neurodegeneration model that couples a reduced-order toxic propagation equation, derived from a Smoluchowski system of oligomer aggregation [68], with a first-order evolution equation accounting for the adverse effects of toxic proteins on the brain's clearance systems. The model is used to examine the brain's resilience to an invasion of toxic proteins, and explores the relation between clearance, brain topology, neurodegenerative disease etiology and progression. We establish a regional measure of resistance to misfolded proteins and demonstrate that network connectivity may protect the brain from invasion. Further, we show that clearance may delay the onset and progression of pathology and that spatial variations in clearance can alter toxic protein progression. Finally, we demonstrate that canonical subtypes in Alzheimer's disease can be induced by changes in local clearance. Overall, our work suggests that variations in clearance may play a key role in the formation and propagation of neurodegenerative pathology, accordance with physiological observations.

Modeling toxic protein and clearance interplay
Network neurodegeneration models are motivated by two primary observations [76]. First, the propagation of toxic proteins is biased by axonal bundles connecting cortical and subcortical regions. Second, the prion-like hypothesis of neurodegeneration postulates a local, autocatalytic reproduction mechanism for misfolded proteins.

Network models of toxic protein propagation
A natural mathematical starting point is to consider a concentration p of toxic proteins that diffuse between brain regions along axonal bundles with local dynamics governed by autocatalytic replication and clearance. We assume that the protein distribution evolves on the connectome. A structural connectome is a network G = (V, E), whose node set V contains N nodes corresponding to brain regions of interest (ROIs) and whose edge set E represents axonal fibers between these regions. Given a structural connectome G, p i = p i (t) denotes the toxic protein concentration associated with the node i at time t and represents the average concentration in the corresponding ROI. A general network neurodegeneration model including toxic protein transport and local clearance may then take the form [20]: find the concentration p i for i ∈ V such thatṗ where the superimposed dot denotes the time-derivative, λ i is the clearance level at node i, L ij are the entries of the graph Laplacian L associated with G, ρ is a transport coefficient, R is a reaction relation governing the local dynamics within the ROIs, and p i,0 is the initial value of the concentration on i.
The graph Laplacian L is defined by L = D − W where W is a weighted adjacency matrix of G and D is the diagonal matrix D = diag(d 1 , . . . , d N ), where d i is the weighted degree associated with node i defined by Other normalized forms of the graph Laplacian may and have been used in the literature, but this standard form simultaneously conserves mass and enforces Fick's constraint, thus guaranteeing that no transport occurs between regions with same concentration [53, Supplementary S1]. Specifically, we select the weighted adjacency matrix with entries W ij = n ij / 2 ij , where n ij is the number of fibers along the axonal tract, connecting node i to node j, and ij is the average fiber length between the same two nodes. The quantities n ij and ij arise from fiber tractography [13] and are provided in the connectome data [37]. This choice of weights is a consistent generalization of conservative finite volume methods applied to diffusion problems [69,27].
Here, we will use idealized graphs G as well as a human brain connectome with N = 1015 nodes ( Fig. 1) generated from data from 426 individual Human Connectome Project patients using the Lausanne atlas, a multi-resolution anatomical parcellation [37,13]. coronal view sagittal view Figure 1: A structural human brain connectome with N = 1015 nodes. Each node is assigned a color based on its classification into 83 disjoint anatomical region of interests.

A coupled model of toxic protein transport and brain clearance
Network models can and have been used to investigate various aspects of Alzheimer's disease, Parkinson's disease, supranuclear palsy, frontotemporal dementia, and other neurodegenerative diseases. Most early studies consider diffusion without reaction, i.e. R = 0 [54,1,55,50,74,51,49], or use network approaches as a means of interpreting voluminous imaging data sets [31,73,72]. Recent work has focused on the autocatalytic nature of protein dynamics leading to a local expansion of the toxic population in agreement with the prion-like hypothesis [21,67,53]. In such models, a nonlinear function for R(p i ) is chosen to model autocatalytic exponential growth at small concentration and saturation at larger concentration as observed in longitudinal studies [33].
Here, we further generalize network neurodegeneration models to take into account dynamic and heterogeneous brain protein clearance. Brain clearance relies on several different modalities such as cellular degradation, transport via the blood, and cerebrospinal fluid-mediated clearance pathways [66,29,47], and may vary in a complex manner from region to region. Moreover, brain clearance has been proposed to play a distinct role in the etiology and progression of neurodegenerative diseases [65,66,47,28], but the precise mechanisms are understudied and poorly understood. In particular, we model the effect of toxic proteins on brain clearance homeostasis and the potential implications for neurodegenerative pathology [28,5,10,9].
At the local level, it has recently been shown, using a Smoluchowski model for the dynamics of proteins, that the effect of reducing clearance is to create an instability. Close to this instability, the dynamics takes the universal form of a transcritical bifurcation [68]. Explicitly, continuing to denote by λ i the level of clearance at node i, the normal form of the bifurcation close to where α describes the population expansion and can be obtained from microscopic models [68]. For λ > λ crit the fixed healthy point p = 0 is stable, but it loses stability as λ decreases below λ crit . For λ < λ crit , the fully toxic state instead becomes a stable fixed point. In light of these findings, we define R in (1) by (3) and the global evolution of the toxic protein concentration is thus governed by, for i = 1, . . . , N : Vice versa, the presence of toxic protein oligomers affects the brain clearance pathways [5,10,9,11,32,15,8,26]. We model a deteriorating clearance due to the presence of toxic proteins by a first-order rate law:λ for i = 1, . . . , N , where β i > 0 is a kinetic constant for each i, λ ∞ i is the minimum regional clearance value. We assume that λ i (0) ≥ λ ∞ i . In summary, the full dynamic system is given by In the case where λ i (0) = λ ∞ for all i, this system reduces to the standard Fisher-Kolmogorov system [20,75,53,61].

Stability and criticality in the homogeneous coupled model
If the initial conditions and other model parameters are homogeneous i.e. p i,0 = p 0 , λ i,0 = λ 0 , λ ∞ i = λ ∞ , and β i = β for i = 1, . . . , N , the system (6) reduces to a homogeneous system equivalent to that of a single node:ṗ

Fixed points and stability
This dynamical system (7) admits a class of fixed points corresponding to the absence of a toxic protein load combined with any clearance λ ∞ ≤ λ ≤ λ(0): Another fixed point of (7) corresponds to the case in which p * > 0 and the clearance is at its minimal value: The properties of the Jacobian J of (7) determine the stability of these fixed points. Specifically, we have that: and its eigenvalues are given by The eigenvalues of J evaluated at (p * ,2 , λ * ,2 ) both have negative real parts; this fixed point is thus unconditionally stable. For the other fixed points, note that the Jacobian's determinant vanishes at (p * ,1 , λ * ,1 ). However, its local dynamics can be evaluated directly by differentiating (7a) with respect to p and evaluating the result at p * ,1 = 0. From these calculatations, we note that the fixed point(s) (0, λ) are stable if and only if λ > λ crit , and unstable otherwise.

Clinical characterization of fixed points
Combining these fixed points and their stability with their clinical interpretation, we define three regional (nodal) homogeneous disease states: healthy, susceptible and diseased as follows. Each disease state corresponds to fixed points of (7). I) In the healthy state, (p, λ) = (0, λ) with λ > λ crit . By (7b), even a small toxic load p leads to a reduction in clearance at a rate β > 0. However, (7a) dictates that p decreases towards to zero as long as λ > λ crit . This state is thus equipped with a degree of resilience to toxic protein seeding and is as such protected from the onset of pathological neurodegeneration.
II) In the susceptible state, (p, λ) = (0, λ) with λ ≤ λ crit . In this state, clearance is dysfunctional, and a small change in the toxic protein concentration will send the system to the diseased state.
In this case, the clearance is at its minimal value and the toxic protein concentration is saturated.
The regional characterizations provide a direct perspective on a clinical characterization of the full system (6). First, observe that any region whose state is (9) is, effectively, an incubator of toxic protein seeds. Even if the neighbors of a diseased region are otherwise in the healthy state, (7b) dictates that toxic protein seeds, from the diseased region, will erode the otherwise healthy neighboring clearance values. Toxic seeds will continue to originate in the unstable node, and migrate to adjacent neighbors, until the clearance of the diseased region's neighbors satisfy λ ≤ λ crit and, in turn, they transition to diseased regions themselves. Once a node has transitioned to (9), a toxic infection propagates outwards from it and the cascade of deficient clearance, and subsequent invasion, permeates throughout the brain. Thus, we say that the whole brain, as represented by the non-homogeneous system (6), is in the healthy state if all regions satisfy (I), that it is in the susceptible state if at least one region satisfies (II) and that it is in a diseased state if at least one region satisfies (III).

Critical toxic seeding
The healthy and susceptible states differ by whether λ > λ crit or not. Consider the phase plane in Fig. 2, where orbits for different initial toxic loads p 0 are shown. For λ 0 > λ crit , each toxic seeding event, yielding p > 0, will degrade the clearance capacity until λ reaches λ crit . For sufficiently small seeds, the toxic proteins are cleared and p tends to 0. However, for larger seeds, p tends to the diseased p * ,2 > 0. Hence, there is a value p(0) = p crit , the critical toxic seeding, that sends the system to the susceptible state (0, λ crit ). Conversely, this value is an indicator of the system's resilience to toxic seeding events.
To derive an analytical expression for p crit in terms of the model parameters, we examine orbits in the phase plane on the form p = p(λ) and integrate. More precisely, eliminating time derivatives between (7a) and (7b) yield After integrating with respect to λ and inserting (p(0), λ(0)) = (p 0 , λ 0 ), we obtain where r = α/β. The critical toxic seeding is given by the orbit intersecting the λ-axis (p(λ) = 0) at exactly λ = λ crit with p 0 = p crit , and is thus given by The critical seeding provides important insight into the dynamics of the homogeneous coupled system. Initial seeding values below the critical threshold p 0 < p crit will result in a healthy steady state (0, λ * ,1 ) where λ 1, * is the largest, strictly positive root of p(λ) = 0. Conversely, an initial seed of p 0 > p crit result in the diseased state.

Network connectivity increases brain resilience via diffusion and clearance
We now turn from the homogeneous, single-node case to the network case. The analysis of Section 3 demonstrates that clearance effectively contributes to reduce toxic protein load. Next, we will show that a node's local connectivity can increase its resilience against toxic proteins by relying on the clearance of neighboring regions.

Perturbation analysis of critical network toxic seeding
Consider a node i and its neighborhood G i defined as the set consisting of i and the indices of all connected nodes: j ∈ G i iff L ij = 0. For small ρ 1 (i.e. slow diffusion), the graph Laplacian term in (6) adds a regular perturbation to the homogeneous system (7). Thus, to investigate the effects of connectivity on the critical toxic seeding, we expand the concentration p k and clearance λ k for each node k in G i with respect to ρ, as Substituting these into (6), equating powers of ρ and dropping the O(ρ 2 ) terms yield two sets of equations for each k ∈ G i :ṗ As initial conditions for (13), first consider a toxic seed p S at node i only and a uniform initial clearance λ 0 > λ crit : Integrating (13) with (15) allows us to express (14) as a simple inhomogeneous system given bẏ for j = i. The solution of this equation is, for all nodes j in G i with i = j, The critical seeding value p crit , subject to the ρ perturbation, is yet to be determined. By definition, if the seeding node is seeded at a level below critical (p i (0) = p S < p crit ), then p i,0 decreases, monotonically by (6a), to the asymptotic state p 1, * i = 0. As a result, (17a) implies that and we conclude that p j (t) < p S for small perturbations ρ. Therefore, when p S < p crit , the toxic protein concentration p j and clearance λ j in node j reach the steady state (0, λ j, * ) with λ j, * > λ crit . The task that remains, then, is to ascertain how p crit changes as a function of ρ. Due to the initial condition (15), the evolution equation (13) gives p j,0 (t) = 0 for j = i so that (14a), for k = i, is given byṗ The set of equations (13) and (19) alongside (14b) can now be solved numerically using standard ordinary differential equation solution algorithms, and the perturbed solutions be reconstructed using (12), to quantify the critical network toxic seeding for different connectome configurations. We investigate the impact of alterations in diffusivity and connectivity on the critical network toxic seeding. First, to quantify the effect of diffusion, we estimate the initial condition p S needed to reach the asymptotic state (0, λ crit ) via the afore numerical procedure, for different ρ. For each experiment, we consider a fixed network consisting of one node with five neighbours, let w ij = 1, λ 0 = 2.0, α = 2.1, and β = 1, and consider a range of ρ. The results demonstrate that the critical network toxic seeding increases with increased diffusivity (Figure 3).
Second, we are interested in the effect of brain connectivity on the critical toxic network seeding. Letting the node degree measure regional brain connectivity, we consider a series of numerical experiments with a node i and increasing node degree d i i.e. an increasing number of connected neighbours. Again, we observe that the critical toxic network seeding p crit increases  Figure 4: The impact of connectivity on the critical network toxic seeding -for a network consisting of one node with N neighbors. Plot of clearance λ versus toxic seeding p S . ρ = 6.4 × 10 −2 , while other parameters as in Figure 3. with the degree (Figure 4). This analysis suggests that the brain's connectivity may protect regions by allowing them to share the burden of clearing toxic loads with their neighbors. These conceptual observations may be interpreted in the context of neurofibrillary tangle (NFT) staging [53,72]. Indeed, (6) dictates that as toxic protein proliferates through neighboring regions, clearance is reduced below λ crit and the formation of NFTs can result. Toxic infection will therefore generally take hold most rapidly in neighbors with lower clearance. Thus, the evolving distribution of clearance, and local toxic population growth, may affect the specific regional sequence of NFT staging. In addition to differences in NFT staging, patient-specific regional variations in clearance may also offer an explanation as to how extra-entorhinal seeding locations might emerge, as hypothesized in [72].

Parameter Value Unit
For instance, examining the relative values of the weighted graph Laplacian degree d ii /d max (Figure 5), we note that the (right) entorhinal cortex (EC) is among the set of poorly supported regions. Thus, toxic τP seeds originating in the EC may tend to linger there. This observation is particularly interesting in the context of AD, consistent with observations from studies of τP staging [7,16,11] in AD, and motivates further research into modeling the region-specific balance and toxic protein load and of τP-related clearance mechanisms.

Dynamic brain clearance alters toxic protein progression
Finally, we investigate the full model (6) at the organ level by direct simulation. The simulations use a common set of model parameters (Table 1) where ρ and α were selected to produce maximal rates of toxic protein increase approximately on par with recent modeling studies employing AD imaging data [61,60,62]; λ crit was chosen in line with experimental studies of aggregation kinetics [68]; λ ∞ i reflects the assumption that the regional clearance can reach a low  but non-zero value; and β i was chosen to be one, which is consistent with a typical time-scale for disease progression of about 30 years. The computational results suggest that the distribution of clearance, throughout the various regions of the brain, may play a significant role in delaying disease onset, in producing the varied patterns of disease progression and can also serve as a mechanism that may explain some canonically studied AD subtypes.

Clearance delays disease onset and progression
We begin by mimicking a progression of τP in AD by placing an initial average toxic seeding p 0 = 0.1 in each of the bilateral entorhinal cortices, alongside an initial clearance there of λ 0 = λ ∞ . All other regions of the brain were initialized with p 0 = 0.0 and λ 0 = γ sim λ crit (see Table 1). A series of simulations for seven different values of γ sim ranging from 5% to 65% were performed. Note that the toxic protein concentration at each connectome node will saturate to the asymptotic value of p * = 0.23 with this model set-up. The computational results show that an increase in the level of healthy, homeostatic clearance delays the onset and progression of neurodegeneration ( Figure 6). The lowest and highest initial clearance rates tested (λ 0 = 5% λ crit , λ 0 = 65% λ crit ) yield onset times of t = 37.6 and t = 93.2 years, respectively, corresponding to a relative increase of 148%. A nonlinearly increasing relationship was noted between the initial clearance and onset time across the simulations ( Figure 6, top right). Improved brain clearance, especially before neurodegenerative onset, may 5% 15% 25% 35% 45% 55% 65% thus have significant benefits to brain health.
Regional toxic burden was also seen to vary with initial clearance (Figure 6, bottom). In particular, and in line with the whole-brain average concentrations, the toxic protein load in any fixed region, at the median arrival time of t = 53.1, decreases with increasing initial clearance. In addition, a higher initial clearance is associated with a limbic and temporal predominance. Decreasing values of initial clearance are associated with increased temporal, parietal and frontal burdens. Furthermore, the observed progression of τP burden, as a function of initial clearance from right to left in Figure 6 (bottom row), is similar to experimentally observed τP NFT progression [35, Fig. 1f].

Spatial variations in clearance alter toxic protein progression
Neuroimaging studies suggest that the clearance capacity within the brain varies regionally and may be altered by age-related factors [66]. Assessments of perfusion [80], ubiquitination [63,64,39] and perivascular CSF circulation [58,17,18] point to specific regional differences as well as temporal variation in the major modes of brain clearance [66]. Here, we will demonstrate that regional variations in initial clearance can cause striking differences in the propagation of protein pathology. Since the connectivity of brain structural connectomes is complex, we first use an idealized geometry in order to demonstrate that clearance perturbs the flow of protein pathology by producing a toxic front that moves orthogonal to the gradient of the clearance field. We next extend these ideas to the connectome in Section 5.3.
We first consider four test cases defined over a uniform lattice of the unit square comprised of 100 equally spaced grid points in each direction. Each test case is initialized with the same toxic seeding concentration, p 0 = 0.01, at the origin node (0, 0) with initial clearance set to uniform Gaussian diagonal linear λ ∞ . At all other nodes p i (0) = 0. We define four different distributions for λ i (0): labelled as uniform, Gaussian, diagonal and linear, shown in Figure 7 and described further below. For the uniform case, we set λ i (0) = λ ∞ , ∀i. For the Gaussian case, the initial clearance is set according to a Gaussian distribution with λ i (0) = λ ∞ + |γ|, where γ is a real value selected from a normal distribution with with mean µ = λ ∞ and a standard deviation of σ = λ ∞ /2. For the diagonal case, we set λ i (0) = λ ∞ in a central diagonal band and λ i (0) = 3λ ∞ otherwise. The linear case sets an initial clearance field by λ i (0) = λ ∞ at all nodes along the y = 0 line, and increasing linearly towards a maximum of λ i (0) = 3λ ∞ along the line y = 1. Other parameters are taken from Table 1 and we set the total simulation time to be t = 100. The corresponding simulated toxic protein progressions are shown in Figure 8. To increase the visibility of the flow front, toxic protein concentrations near zero are transparent while those coinciding with the onset value, determined by half of the maximal saturation value of p i = 0.23 (purple), are opaque. The toxic propagation first develops orthogonal to the clearance gradient, if such exists, subject to the underlying graph topology. With the uniform initial distribution of clearance, there is no gradient and the toxic front is constrained only by the topological connectivity of the graph. Similarly, the (discrete) gradient of the initial Gaussian clearance field is also Gaussian with µ = 0, and σ ≈ λ ∞ / √ 2. Hence, a clear sense of orthogonality is lacking also in this second test case, and we observe that the toxic front spreads in all connected directions. The diagonal test has a sharp gradient in the initial clearance field which is clearly reflected in the toxic protein propagation pathway. Finally, the fourth test case exhibits a constant clearance gradient oriented along the y-axis and the resulting toxic front advances first along the x-axis before propagating upwards.
These results both echo and extend the observations of Section 5.1. The pattern with uniform or Gaussian initial distributions are similar, but with a smaller time scale in the Gaussian case due to a lower mean clearance. The case with diagonal or linear initial distributions extend this perspective by demonstrating that the direction of the initial clearance gradient can significantly alter the evolution of pathology and that pathology spreads most rapidly in the direction orthogonal to the gradient of clearance. Overall, these results strongly suggest that variations in brain clearance may significantly alter the patterns of toxic protein deposition in neurodegenerative diseases and may have further implications for the various trajectories [72] of τP deposition related to Alzheimer's disease.

Clearance variation may promote AD subtypes
Recent studies have used τP progression to define notions of AD subtypes and have assumed that these different pathologies stem from different seeding regions [19,72]. Here, we test the alternative hypothesis that subtypes can arise from the same seeding region but with regional differences in clearance.
To define AD subtypes from histopathology, post-mortem NFT distributions of nearly two thousand patients of Braak stage V or later were collected in previous studies [44,77,34]. For each brain, these studies counted NFTs in the hippocampus (HP) and in the association cortex (ASC), where the latter was defined by superior temporal, middle frontal and inferior parietal regions (Figure 9). The ratio of HP to ASC NFTs (scores) was then computed, and the overall cohort distribution of values was determined [44]. The AD subtype classification is as follows: Hippocampal sparing: In this subtype, NFTs invade the association cortex more than the hippocampal region. It is defined by scores of less than the 25th percentile in the cohort distribution.
Typical AD: Here, NFTs invade both the association cortex and the hippocampal region and is defined by scores between the 25th and 75th percentiles in the cohort distribution.
Limbic predominant: NFTs invide the hippocampal region more than the association cortex and this subtype is defined by scores larger than the 75th percentile in the cohort distribution.
We will here demonstrate that simple variations in the initial distribution of clearance can elicit variations in the observed patterns of toxic protein progression and explain these AD subtypes.   To compare the effect of clearance on the distribution of NFTs, we follow previous studies [55,67,53,60] and augment (6) with a measure of (nodal) NFT production, denoted by q i (t), reflecting damage accumulation following the arrival of toxic proteins. Given a toxic protein concentration p i (t), the (post-processed) NFT aggregation marker is defined as the solution of the damage equation:q The variable q i is a local damage variable that increases from 0 to 1 as the disease progresses.
To measure the influence of clearance on the HP and ASC regions, we use the open-source NetworkX software package [23] to define influential regions for each. We consider the prevalence of the connection strengths between the nodes of a given composite ROI (either the HP or the ASC) and the nodes of its immediate neighbors; and by assessing how frequently a region appeared in shortest paths that originated in the EC and terminated in the HP or ASC composite ROIs (see Table 2).
We investigate the progression of AD tauopathy by solving (6) and subsequently (20) for thirteen simulation scenarios. Each case uses the default model parameters (Table 1) along with initial clearance values for the HP influential region nodes set to λ i (0) = λ crit , and the initial conditions (p i , λ i (0)) = (0.01, λ ∞ ) in the bilateral entorhinal cortices. The other initial clearance values λ i (0) are defined as follows. We consider 13 equispaced values of M ∈ [0.7, 1]. For each M , the initial clearance in the ASC influential regions is set to λ(0) = M λ crit . The initial clearance in all other regions is set to λ i (0) = 1.8λ crit . Next, from the evolution of the damage q i in each region, we obtain q HP as the average of (20) over the nodes of the hippocampal region and analogously for q ASC . For each simulation, we record the onset time at which either q HP or q ASC first reach 50%. Finally, each simulation result is classified according to the post-mortem methodology of determining the ratio of q HP /q ASC and assigning the subtype category based on its quartile range [44,34,77].
Our computations reveal that higher initial clearances in the ASC influential regions (M ∈  Figure 10: Simulated AD subtypes with average NFT aggregation marker (top row) in the HP (q HP in green) and ASC (q ASC in blue). For each M , we compute the time of onset which is defined as the first time that either q HP or q ASC reaches 1/2. The value of the NFT aggregation marker (20) in the ASC (middle row) and HP (bottom row) are then plotted, at time of onset, for each of the corresponding NFT aggregation plots (top row).  Figure 10 shows four representative examples along the simulated type spectrum. Moreover, we find that the onset times cluster into similar groupings ( Figure 11). These results reproduce the clinical observation that the hippocampal sparing variant reaches onset before typical AD, which itself reaches onset before the limbic predominant variant [44,34,77]. We conclude that regional variations in brain clearance may explain AD subtypes -an observation that motivates further studies in this direction.

Conclusion
Network neurodegeneration models have been widely used in both the study of patient data [31,73,61,72,60] and in the examination of potential mechanisms underlying neurodegenerative disease pathology [54,55,50,20,21,67,68]. Experimental evidence suggest that clearance systems in the brain may play a fundamental role in neurodegenerative disorders [65,66,79,56,3,42,47,28]. It is therefore natural to include such effects in the framework of network models and study the possible role of clearance theoretically.
Our model starts with a healthy brain that has sufficient clearance to eliminate small amount of toxic proteins. However, as toxic proteins increase, the clearance system becomes increasingly damaged and the brain is subject to full invasion. In the absence of evolving clearance, the region that have sufficient clearance will always be protected. Our model of coupled clearance and neurodegeneration provides insights for future research. Our analysis suggests that the brain may exhibit clearance-dependent regional homeostasis and that the topology of the brain may provide resilience against a toxic protein infection taking hold. Our simulations, motivated by the progression of τP in AD, further suggest that increasing clearance in healthy brain may be instrumental in significantly delaying AD onset, that the progression of toxic pathology depends on regional clearance levels and that patient-specific distribution of brain clearance may play a key role in the manifestation of AD subtypes. Our study has a number of limitations that provide a direct opportunity for continued research. First, neurodegenerative diseases often involve fundamental interactions between brain proteins [30,4,57,71] and clearance mechanisms [70,66,52]. However, our model only includes one toxic protein species and combines distinct brain clearance systems into a single (regional) term. Extensions of our model to multi-species, as proposed in [67], or multi-clearance models, all represent avenues for further theoretical development. In addition, our model has yet to be assessed using regional clearance or toxic concentration values garnered from neuroimaging data sets, or derived directly from an experimental setting. Despite these limitations, our network model coupling clearance to toxic protein pathology represents a clear and promising avenue for investigating the many consequences of evolving brain clearance in the etiology and progression of neurodegenerative diseases.