Braiding Braak and Braak: Staging patterns and model selection in network neurodegeneration

The term ‘neurodegenerative disease’ refers to a family of maladies primarily affecting the brain’s neurons; many neurodegenerative diseases result in cognitive decline and, ultimately, a diagnosis of dementia. A hallmark of such diseases is a large concentration of toxic protein aggregates throughout the brain. These toxic proteins interfere with normal brain functions and are associated with neuronal impairment, neuronal loss, brain atrophy, and overall cognitive decline. Here, we address questions related to model selection for computational studies of Alzheimer’s disease (AD); AD is the most contemporarily prevalent cause of dementia. The practical constraints of clinical AD experiments, most especially in humans, have lead to an interest in the development of mathematical models of pathology evolution, particularly those that can accurately capture large-scale features of AD. Many of these models take advantage of the prion hypothesis which asserts that the progression of proteopathy, in AD, follows from a prion-like mechanism.

Prion-like mathematical models can be broadly characterized as either probabilistic (Vogel et al., 2020, 2021) or as continuous (Fornari, Schäfer, Goriely, & Kuhl, 2019a; Raj, Kuceyeski, & Weiner, 2012; Thompson, Chaggar, Kuhl, & Goriely, 2020). Both model types often invoke nonlinearity and are typically discretized, and solved, on either a simplicial mesh (Weickenmeier, Kuhl, & Goriely, 2018), or on the network defined by the brain’s connectome (Fornari et al., 2019a; Raj et al., 2012; Schäfer, Mormino, & Kuhl, 2020; Thompson et al., 2020). Probabilistic and continuous models have their individual advantages. The former can often reveal unique insights in complex datasets while the mathematical structure of the latter can be rigorously analyzed or, in some simple cases, closed form solutions can be derived. Misfolded protein aggregates, in neurodegenerative diseases, often display structured patterns of progression, and reproducing these staging patterns is a desirable trait for any mathematical model. Here, we will discuss the impact of model selection, for continuous prion-like mathematical models, with respect to an observed structured staging sequence. As a case study, we will consider the network progression of misfolded tau protein (τP) in AD. We compare results to staging patterns motivated by a recent study (DeVos et al., 2018) that uses a six-stage Braak pattern for τP in AD. However, our braid surface approach is generalizable to more complex hierarchical patterns (Cho et al., 2016; Delacourte et al., 1999) and well suited to broader investigations, including applications to the study of staging phenomena in other neurodegenerative diseases.

The primary contribution of this manuscript is a methodical, quantitative study of model selection features, as they pertain to structured staging, in the setting of a continuous mathematical model of a prion-like proteopathy processes defined on structural networks. This manuscript introduces two novel investigative tools, braid diagrams and braid surfaces, facilitating this process. Our approach is of practical interest for neurodegenerative diseases for two main reasons. First, there is often a lack of the longitudinal radiotracer data, in which staging would be evident, that would provide for accurate, automatic model selection via, for instance, a parameter inference approach. Such studies are important for short-term predictions, but they may make model choices based on a limited number of time points and arrive at a model that does not capture a long-term expected staging progression. Second, we show that there are several other factors of a network neurodegeneration model, beyond mathematical parameters, which are not always considered in the model selection process and can play a significant role in staging. In particular, we show that the choice of structural connectome scale, tractography method, thresholding method, graph Laplacian weights and, sometimes slight, variation in model parameters can significantly alter the landscape of observed, even accessible, staging patterns.

Our study highlights the nuanced role played by each choice in the network neurodegeneration modeling pipeline, providing important insight, and quantitative tools, that inform model selection. Even though we choose a particular example of an AD staging as a point of study, the methods and observations are general and can be applied to the study of any hierarchical staging process that can appear in the study of neurodegenerative diseases or other dynamical process on networks. The source code, and connectomes, to create the braid diagrams and braid surfaces discussed in the manuscript are freely available online (Putra, Chaggar, Thompson, & Goriely, 2021).

Let Ω_j, for j = 1, 2, … J be a nonoverlapping collection of nodes; that is, Ω_j and Ω_k are disjoint subsets of V when j ≠ k. Let T ∈ [0, 1] be an arbitrary, but fixed, threshold value. As p evolves, according to Equation 1a, the average concentration is computed in each region Ω_j and the time when Ω_j first reaches the threshold T is recorded. This process produces an ordering of the regions Ω_j and the ordered sequence Ω_j1, Ω_j2, …, Ω_jJ is called an observed staging pattern. The generalized staging problem is to ascertain the scope of observable staging patterns, subject to variations in Θ. The observed staging patterns can then be compared to one or more desirable staging patterns.

We note that other authors (Abdelnour, Voss, & Raj, 2014; Pandya, Mezias, & Raj, 2017; Pandya et al., 2019; Raj et al., 2012, 2015) have used normalized forms of the graph Laplacian for their models of neurodegenerative disease progression and a natural question for model selection is the choice of a suitable Laplacian for the underlying process.

A second condition that we impose on the transport is that transport is driven by difference in concentrations. If the concentration is equal at two nodes, there is no driving force to create an imbalance between these nodes. By analogy with diffusion processes based on Fick’s law, we call this condition the Fick’s condition. In terms of the graph Laplacian on an undirected network, it reads simply L_a,b · 1 = 0. We note that there are other possible assumptions when the transport process is viewed as a probabilistic event as shown in Masuda, Porter, and Lambiotte (2017).

We show in Supporting Information S1 that the only possible graph Laplacian, in the class defined by L_a,b, that satisfies both mass conservation and the Fick’s condition is the standard graph Laplacian L = L_0,0. We also show how to generalize this result when regions of different volumes are considered and Fick’s condition suitably generalized.

We want to emphasize that this unavoidable constraint on a physical model of transport does not invalidate previous studies based on the normalized graph Laplacian as these have shown great predictive and explanatory powers when applied to actual patient data. Yet, the choice of using a normalized graph Laplacian is not innocent and should be fully justified as it has some direct implications on the assumptions used in the model, even if these assumptions are not usually given explicitly.

We consider networks G = (V, E) based on two families of multiresolution structural connectomes generated from Human Connectome Project (HCP) data (see Methods, Structural connectomes). Each edge, e_ij, is associated with two values: the number of fibers n_ij constituting the connection between (anatomical region) nodes i and j; and the fiber length ℓ_ij. We consider the choice of weights as part of the model selection and study three possible weights based on the literature: the length-free weighting (Abdelnour et al., 2014; Raj et al., 2012, 2015), the ballistic weighting (Fornari et al., 2019a, 2019b), and the diffusive weighting (Thompson et al., 2020). The formulas for these three weighting schemes are listed in Table 1.

A primary contribution of this manuscript is the introduction of braid diagrams and braid surfaces; these powerful tools present a direct and visual assessment of an otherwise complex, nonlinear process evolving in time. We begin by describing a braid diagram. Let G = (V, E) be a fixed network and suppose that Ω₁, Ω₂, …, Ω_J are a fixed set of nonoverlapping node regions. Suppose further that T₁, T₂, …, T_N are (biomarker) threshold values in the unit interval [0, 1]. A braid diagram is a graph whose abscissa is the index of the regions Ω_j and whose ordinate corresponds to the threshold values T_k. As a dynamical system, such as Eq. 2a, evolves on G, the time t_j,k at which each region, Ω_j, first achieves each threshold, T_k, is recorded. If a given threshold is never achieved in a particular region, the recorded time is prescribed as t_j,k = ∞. The collection of time values, for a fixed threshold index, establishes an ordering of the regions Ω_j. For each threshold, the ordering can then be visualized as a braid diagram.

We are also interested in staging outcomes as the system parameters are varied. This information is contained in a braid surface that generalizes a braid diagram. For each possible staging, we assign a color. The braid surface is a two-dimensional plot that assigns for each value of one parameter and one biomarker the corresponding staging color. Computationally, this surface is generated by a simple algorithm. First, one discretizes the continuous values of a parameter, such as β in Eq. 2a, of interest. Then, for each discretized parameter, the underlying system (e.g., Eq. 2a) is solved, a braid diagram is constructed and the observed staging pattern is determined for each threshold value; if an observed staging pattern has not yet been encountered, it is added to a list. At the end of the process, every pairing of discrete parameter and threshold has been assigned to an observed staging pattern. The set of observed staging patterns are assigned to colors, and these colors are plotted to visualize the braid surface. The τP seed staging braid surface for the example network, and weighting matrix, of Figure 1 is shown in Figure 2. The green corresponds to the β-parameter, and threshold, region where a I → III → II → IV staging pattern is observed while the red corresponds to the parameter region where a I → III → IV → II is produced. The braid surface shows that the staging dynamics, for the four-vertex system with diffusive weights, are simple; we will observe complex staging dynamics on structural connectome graphs of the brain.

The problem of the structured staging of τP NFT in AD has been well studied. The most famous study being the seminal work of Braak and Braak (Braak, Alafuzoff, Arzberger, Kretzschmar, & Del Tredici, 2006; Braak & Braak, 1991) that looked at the progression through six regions. The six-region progression view is still in contemporary use (DeVos et al., 2018; Vogel et al., 2020), though some authors have also proposed refinements to ten regions (Delacourte et al., 1999). More recently, the hierarchical progression of flortaucipir tracer has also been studied for progression patterns; authors have studied standardized uptake value ratio (SUVR) progression using a canonical 6-region (Schöll et al., 2016) and an extensive 25-region (Cho et al., 2016) pattern.

Much of the staging literature refers to the evolution of τP NFT or of SUVR quantities, though a recent study (DeVos et al., 2018) has advanced the notion that τP seeds precede NFT pathology in a similar structured manner. The model that we study, that is, Eq. 2a, tracks the progression of both τP seeds (via Eq. 2a) and τP NFT (via Eq. 2b). Therefore, we adopt the staging regions given by DeVos et al. (2018); we will observe how the many, sometimes surprising, aspects of model selection can strongly influence the set of observed regional progressions for a computational model of proteopathic τP staging in AD. The same method can easily be generalized to the study of other staging regions such as those advanced in Cho et al. (2016), Schöll et al. (2016), Vogel et al. (2020), or Delacourte et al. (1999).

We studied the staging problem for the model (Equation 2a) of τP proteopathy in AD, using the braid surface approach. To study this problem, a five-region adaptation, to match the structural connectome ROI, of the six regions used in DeVos et al. (2018) was adopted and is shown in Table 2.

We refer to staging progression, determined by a concentration threshold 0% < T ≤ 100%, through these five regions with Roman numerals and use the I → II notation. To present the braid surface results, we first enumerate all observed regional staging patterns and assign each one a color as shown in Table 3. From the original point of view (Braak & Braak, 1991) of τP Braak staging, one would expect that the I → II → III → IV → V progression is a likely candidate of interest, at least for NFT. However, based on Alzheimer’s Disease Neuroimaging Initiative (ADNI) data, we also identified three other observed progression sequences that may be of interest for both τP seeds and NFT. The set of identified staging patterns of potential interest are called Suggested Patterns and summarized in Table 3 (top).

We compared two different types of connectome tractography methods: both deterministic tractography, prepared by The PIT Bioinformatics Group (2019) using the MRtrix software, and probabilistic tractography, prepared for this study using the FSL software. These generating methods used parcellation ROIs defined by the Lausanne multiresolution atlas and constructed from HCP data. The braid surfaces for the deterministic connectomes and the probabilistic connectomes are, in general, quite different. However, there are some notable similarities. For instance, both connectome types show great sensitivity as the parameter β increases to large values. In this diffusion-dominated regime, the observed staging pattern becomes nondistinct with multiple overlap for small changes of parameters. Another common feature is that there are large areas, of T and β parameter space, where the τP seed staging and τP NFT staging are consistent, meaning that both display the same staging pattern. A final notable similarity is that both connectome types are able to produce the progressive Braak () pattern for NFT staging.

Despite these similarities, there are also clear differences between the staging behavior observed on the deterministic and probabilistic connectomes. In particular, we examined ADNI SUVR data subject to an ROI selection from a recently published study (DeVos et al., 2018) on τP seeding along the Braak pathway. The SUVR staging pattern suggested () by our analysis appears prominently on the deterministic connectomes (Figure 3 and Figure 4) but is nearly absent on the probabilistic connectome counterparts (Figure 5, Figure 6, Figure 7, Figure 8, and Figure 9). In fact, the deterministic connectomes can express all four suggested staging patterns (, , , and ) at high resolution for both τP seeds (Figure 3) and NFT (Figure 4) while the probabilistic connectomes almost always express fewer patterns.

We also observed (Figures 5 through 9 and Supporting Information S2) significant qualitative differences in the braid surface landscapes when various thresholding methods were used in preparing the probabilistic connectomes (Table 4). Therefore, we conclude that connectome preparation plays an important role in all downstream results and is an integral part of model selection.

Several graph Laplacian weighting schemes have been suggested to model disease dynamics. But, it remains unclear as to which of these weighting choices is most biologically sound. Our results suggest that this choice also affects the observed computational staging patterns in two ways. The first of these concerns τP seed staging within the parameter regime where diffusion takes a more prominent role (i.e., large values of β). The move from length-free (LW) to ballistic (BW) weights tends to decrease the surface area where the observed computational staging is sensitive to small changes in the biomarker abnormality threshold; this surface area is again reduced when moving from BW to diffusive (DW) weights. This can be seen in both the deterministic (Figure 3) and probabilistic (Figures 5 through 9, top two rows) braid surface results. In this sense, the length dependence in the weights ‘stabilizes’ the τP seed staging, appearing with an increase in β, in Eq. 2a.

The second general observation concerns the prevalence of the progressive Braak () and SUVR-suggested () staging patterns. In the former case, the general trend is that the choice of a DW graph Laplacian tends to increase the prevalence of the progressive Braak () NFT staging (Figure 4 and Figures 6 through 9, bottom two rows). However, we remark that this observation can depend on the choice of thresholding technique as we observed in the case of the density filtered probabilistic connectome where the LW had this effect (Figure 5, bottom two rows and Supporting Information S2 Figure S3, bottom two rows). However, if one were comparing computational results directly to SUVR staging (e.g., Schöll et al., 2016; Vogel et al., 2020), using the ROI of DeVos et al. (2018), the observation is more complex. In the case of deterministic streamlining, the prevalence of the observed SUVR-suggested staging () can decrease, within the ranges of β we considered, as one moves from LW to BW and again to DW (Figure 4, right-most column). However, we do note that a number of undesirable staging patterns are observed, for LW, when ln(β) > 0 on the highest resolution connectome (Figure 4C); thus, extra consideration should be given to the mathematical model parameters in this case. For probabilistic connectomes, the appearance of the SUVR-suggested staging pattern () depends more, in general, on thresholding technique and the parcellation resolution; the choice of weights had only a weak impact in this specific regard. Our results suggest that, in the absence of a clear biological impetus for selecting particular graph Laplacian weights, a braid surface analysis can play an important, practical role in selecting amenable weights for a particular study.

Network neurodegeneration studies have evolved systems, such as Eq. 2a, on both low-resolution (Fornari et al., 2019a; Goriely et al., 2020; Raj et al., 2012, 2015; Schäfer et al., 2020) and high-resolution (Thompson et al., 2020) connectomes. We studied how connectome resolution alters observed computational staging patterns. In the case of connectomes generated with deterministic streamlining, increasing the connectome resolution improved the overall observed staging results. For instance, we see a stabilizing effect on τP seed staging, when ln(β) > 0, (columns of Figure 3) with increased connectome resolution. Moreover, the combined area of the progressive and SUVR-suggested staging patterns of interest ( and , respectively) tend to increase in prominence with increased resolution; conversely, the transposed prefix stagings ( and ) are more pronounced at lower resolutions. This observation also holds for NFT staging (Figure 4).

In the case of probabilistic streamlined connectomes, the thresholding method and the effect of resolution are intertwined (Figures 5 through 9). It is invariably clear that resolution has a pronounced effect on observed staging, but those effects varied depending on the thresholding method. In general, the low-resolution connectomes consisted primarily of a single staging pattern (typically either or ) when growth was sufficiently strong (circa ln(β) < −5) and increasing the connectome resolution increased heterogeneity in addition to introducing the progressive Braak staging pattern ().

Together, our results show that connectome resolution plays a deciding role in the landscape of observable staging patterns. In particular, the seemingly ubiquitous practice of the use of low-resolution connectomes may partially explain discrepancies between modeling results and comparisons to imaging data in the latter Braak regions (e.g., IV and V).

Deterministic connectomes were constructed using deterministic streamlining (see Methods). We report the observed staging for both τP seeds (Figure 3) and for NFT (Figure 4). We present these results for the lowest (Scale-33), median (Scale-125), and highest (Scale-500) resolution connectomes in addition to all three graph Laplacian weighting schemes. Both figures summarize results for the model parameter value β in the range −10 ≤ ln(β) ≤ 2. Further results are reported in the Supporting Information (S2).

Probabilistic connectomes were constructed from HCP data. We focus on the effects of the parameter β and of several connectivity thresholding methods on observed staging results. For brevity, we present only one thresholding level per thresholding method, but other thresholding levels are presented in the Supporting Information (S2).

In the context of studying prion-like spreading of misfolded proteins in AD, different studies have made use of different model choices. For instance, within the class of reaction-diffusion models such as Eq. 2a, several types of graph Laplacian weights have been used; including those that are free of the influence of length (Pandya, Kuceyeski, & Raj, 2017; Raj et al., 2012, 2015), weights that model prion-like transport as a velocity (Fornari et al., 2019a, 2019b), or as diffusive (Thompson et al., 2020) in nature. Various connectome resolutions, for the Lausanne parcellation, have also been used, both low resolution (Fornari et al., 2019a, 2019b; Goriely et al., 2020; Pandya, Kuceyeski, & Raj, 2017; Raj et al., 2012, 2015; Schäfer et al., 2020) and high resolution (Thompson et al., 2020). Several works have compared, or fitted, spreading models to atrophy (Raj et al., 2015) and SUVR data (Pandya, Kuceyeski, & Raj, 2017; Schäfer et al., 2020; Vogel et al., 2020, 2021). However, to our knowledge, a systematic methodical investigation into the general staging problem for AD has not been advanced, even for the simpler class of prion-like progression models such as Eq. 2a.

In this manuscript we have undertaken a methodical investigation of the generalized staging problem for a simple network model of AD; this important problem concerns the ordered progression of a marker of interest propagating over a network. Using an adaptation of a six-stage sequence of ROIs (DeVos et al., 2018) along the Braak pathway, we have systematically investigated how various aspects of model selection, including those mentioned above, may alter the progression of τP seeds and NFT across the connectome. To do so, we have introduced, and applied, the novel tools of braid diagrams and braid surfaces. Though we have focused on a continuous network neurodegeneration model, our method also applies to probabilistic spreading models as well (Iturria-Medina et al., 2014; Vogel et al., 2020).

Our findings have several implications for model selection. The most surprising is that the particulars of the connectome itself play an important role in the staging problem. That is, when studying τP progression, the construction of the structural connectome should be considered as an inextricable part of the model. To illustrate this, consider two particular staging patterns: the progressive Braak pathway NFT staging pattern () and the staging pattern suggested () by our ADNI flortaucipir data analysis. It is worth noting that the difference between the progressive Braak () versus the SUVR-suggested () staging is whether a biomarker first appears in the hippocampus (II) or the posterior parahippocampal gyrus (III). In DeVos et al. (2018), the histopathologically identified Braak stage II brains expressed τP seeding in the parahippocampal gyrus (III) at a slightly higher average level than in the hippocampus, potentially explaining the difference between the SUVR-suggested () and progressive () staging patterns. That is, () and () may essentially reflect the same effective staging but from two different points of view. Only the () pattern was reliably expressed on the often-used, low-resolution (Scale-33) structural connectome constructed with deterministic streamlining. However, both patterns appeared in the high-resolution (Scale-500) case (Figures 3 and 4, left column vs. right column). Conversely, the high-resolution connectomes generated with probabilistic tractography, and subsequently thresholded for sparsity, almost never produced the () pattern suggested by our ADNI data analysis (Figures 5 through 9, second and fourth rows), though they did produce the progressive staging (). The SUVR pattern () did, however, appear robustly for the Scale-33 connectome when naive thresholding was used (Figure 4). For the Scale-500 naively thresholded connectome, the () pattern does appear but is only consistent with seed staging in a very small region of parameter and threshold space (see Figure 4, second and last rows). Figure 4 suggests that naive thresholding may be a promising strategy for probabilistic connectomes; however, we found that the () and () patterns were sensitive to the threshold used (see Supporting Information S2) and were almost completely absent, at both scales, for the nearby thresholds we examined.

The consequences of these observations are nontrivial. For instance, a study comparing the NFT model Eq. 2b to ADNI flortaucipir data on the low-resolution probabilistic tractography connectomes, thresholded with a density filter, would not be capable of achieving the staging pattern () obtained by our SUVR analysis (see Figure 5, third row), invariably leading to a misfit between model and data. Overall, our observations imply that the methods used to construct the structural connectome, for example, the choice of tractography method, parcellation resolution, and any subsequent thresholding, can affect all downstream results and limit the landscape of realizable staging patterns.

The graph Laplacian plays a distinctive role in reaction-diffusion models of prion-like propagation such as Eq. 2a; it provides the mechanism for the transport of misfolded proteins over the brain’s structural network. As such, the choice of graph Laplacian weights has varied in the literature; particular weight selections are often motivated by a biological argument or an appeal to data fit. Overall, our results suggest that, while the choice of weights often has a minimal impact on the number of observable staging patterns of interest (, , , ), they did affect the prevalence of those patterns, in particular, the diffusive weights (DW) scheme tended to favor the expression of the progressive Braak () staging pattern over the options of length-free weigt (LW) and ballsitic weights (BW).

In the case of τP seed staging, LW and BW showed significantly more sensitivity, in the observed staging, than the DW scheme as ln(β) drew nearer to zero. The staging sensitivity to the weighting scheme can be illustrated by computing the standard deviation of the staging time series. Plots of these standard deviations, for τP seeding on the Scale-500 deterministic connectomes (see Figure 3) are shown in Figure 10. Regions, in Figure 3, where the staging is sensitive to small perturbations in β or T are precisely those where the standard deviation of the staging time series is low. This implies that the observed sensitivity is a race condition where the differences in threshold-attainment time are nearly identical.

As length plays a more dominant role in the weighting scheme, the standard deviation of the staging time series increases, staging patterns are no longer sensitive to small perturbations and τP seed staging becomes pronounced. However, it is not currently known if τP seeds should follow a clear hierarchical progression, though the appearance of τP seeds have been shown to proceed NFT pathology (DeVos et al., 2018). In the case of NFT staging, the DW scheme generally promoted the progressive Braak pathway NFT staging pattern () across all connectomes. This pattern, however, was not the SUVR-suggested () pattern according to our analysis. It is worth noting that the patient data necessary to achieve more statistically significant results, for latter stages, was not present (see Methods, Data preparation and Methods, Identifying additional computational staging patterns of interest). Moreover, as we mentioned previously, the II → III versus III → II paradigm, differentiating the progressive staging () from the SUVR-suggested staging (), may be related to the prevalence of seeds (DeVos et al., 2018) in the posterior parahippocampus at histopathological Braak classification stage II. Nevertheless, our results suggest that direct comparisons to flortaucipir data, at least for the ROI we considered, may be best conducted (see Figures 3F and 4F) with the BW scheme and on high-resolution connectomes created with deterministic streamlining and with a choice of mathematical parameters given by −2 ≲ ln(β) ≲ 1.

Conversely, if a progressive Braak pattern () is desirable, this can most reliably be achieved, in general, on high-resolution connectomes with the DW scheme, sufficiently dominant growth (i.e., β sufficiently small) and a typical biomarker threshold of at least 50%. As a point of practice, this should be checked with a braid surface analysis, as we did observe that thresholding methods can perturb this trend on the probabilistic connectomes that we considered. For example, Figure 5 shows that the density filter threshold favors the LW and BW schemes over the DW scheme for the progressive staging (); likewise, Figure 7 shows that the high-salience scheme more readily expresses the progressive staging () on lower resolution connectomes.

Investigations of τP staging have been carried out histopathologically (Braak et al., 2006; Braak & Braak, 1991; Delacourte et al., 1999), using SUVR (Cho et al., 2016; Schöll et al., 2016), and in fitting models of spreading to SUVR data (Schäfer et al., 2020; Vogel et al., 2020, 2021). As computational power is readily available, the prion-like hypothesis has enabled a new generation of mathematical models facilitating the latter types of studies. However, it is currently not known how various model selection choices may promote, or limit, the expression of particular staging patterns. We have found that every aspect of model selection can alter the overall progression of τP staging. This includes parcellation resolution, tractography method, choice of thresholding method, choice of threshold, graph Laplacian weights, and the mathematical model parameters. We have shown that the the generalized staging problem couples the global topology of the structural connectome to the mathematical model governing the local dynamics of a biomarker evolving on that network and that comprehensive model selection cannot, in general, be ensured from a decoupled perspective. We have introduced braid diagrams and braid surfaces as a means to investigate the complex staging landscape, and changes thereto, for this coupled problem and to observe the implications of model selection choices.

The generalized staging problem in AD reveals the need for more experimental data and poses new theoretical questions. We propose that further study into the progression patterns of τP seed staging, extending the results of DeVos et al. (2018), could promote a more reliable selection of parameter ranges for the balance of growth and diffusion (β) for τP seeds and, due to the consistency we observed, for τP NFT. Moreover, an increased number of SUVR data for patients in the latter Braak stages would help to belay uncertainty in determining the progression patterns suggested by data studies. Mathematically, the study of the coupling between the connectome topology, the weighted distances between ROIs, and the model parameters is a promising candidate for theoretical research.

Braid surfaces have been introduced to facilitate model selection, in the context of staging patterns developing on two sets of undirected, multiresolution structural connectomes; both sets of connectomes were originally constructed from the data of participants in the HCP. All of the structural connectomes considered in this manuscript were constructed using the Lausanne multiresolution atlas parcellation (Daducci et al., 2012) with five levels of potential resolution: the coarsest scale (Scale-33), three intermediate scales (Scale-60, Scale-125, Scale-250), and a fine scale (Scale-500). The edges, at all scales, include information regarding the number of fibers (n_ij) and fiber length (ℓ_ij) associated to each edge (e_ij) connecting region i to region j. The first set of connectomes (see Figure 11) were constructed using MRtrix (Kerepesi, Szalkai, Varga, & Grolmusz, 2017; Tournier, Calamante, & Connelly, 2012) and a deterministic streamlining with 20,000 streamlines and randomized seeding. These connectomes span five scales; the lowest and highest scales are shown in Figure 11. These connectomes are publicly available, as the dataset named ‘Full set, 426 brains, 20,000 streamlines’, from The PIT Bioinformatics Group (2019).

The second set of connectomes were constructed for this study in order to assess degrees of certitude regarding observations for the deterministic connectomes mentioned above. The Connectome Mapping Toolkit (Daducci et al., 2012) was used to parcellate a high-resolution MNI reference template; the FSL (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012; Smith et al., 2004; Woolrich et al., 2009) PROBTRACKX algorithm was employed, for the probablistic tractography, with 10,000 streamlines per voxel. The sparsity of the resulting connectivity matrices was low (approximately 7–12%); these matrices are commonly thresholded before use in computational models. To study the effect of thresholding, we considered five different thresholding techniques in our comparative analysis. The thresholding techniques used are summarized in Table 4, and a description of the method can be found in the corresponding citation. The naive thresholding method removes edge e_ij if the corresponding connectivity coefficient (n_ij) is below a prescribed threshold value.

For comparison, we have constructed connectomes for the lowest (Scale-33, 50 patients) and highest scale (Scale-500, 25 patients) resolutions of the Lausanne multiresolution atlas. The braid surface source code, and the full set of thresholded connectomes, used in this study is available online (Putra et al., 2021).

We have discussed several hierarchical staging patterns for τP progression in AD. As our mathematical model accounts for both τP seeds and τP NFT, we selected a six-stage model that has been related to both quantities (DeVos et al., 2018). These six regions, along the Braak pathway, are enumerated in Table 5.

The regions in Table 5 can be mapped to anatomical ROI of the Lausanne atlas used in the construction of the structural connectomes for this study. We therefore adapt the stages of Table 5 to the Lausanne atlas by considering the staging process of Table 2. The primary difference is that stages V and VI, of Table 5, are combined into a single terminal stage and approximated by the presence of deposition in the cuneus, pericalcarine cortex, lateral occipital cortex, and lingual gyrus.

In order to assess the implications of various model parameters on observed computational staging patterns, thus facilitating model selection, a choice of preferable staging patterns must be identified. We select five collections of nodes, Ω_i for i = 1, 2, …, 5, such that Ω_k is the set of nodes of the computational connectomes whose anatomical ROI labels are given by the k^th stage of Table 2; thus Ω₁ contains all nodes labeled as belonging to the left and right entorhinal cortices, etc. The staging regions, Ω_i, for the coarsest (top row) and finest (bottom row) connectome resolutions are shown in Figure 12.

The work of DeVos et al. (2018) assessed τP seed staging along the Braak tau pathway and found that tau seeding precedes the presence of NFT pathology. A clear candidate for both τP seed and NFT staging is therefore the progressive regional staging I → II → III → IV → V. To assess other potential staging sequences we used flortaucipir (AV-1451) data from the ADNI; the data was fully preprocessed at the Helen Wills Neuroscience Institute at the University of California, Berkeley, and the compiled results are freely available through the ADNI website. The preprocessed dataset contains 1,184 records; each record corresponds to a subject structural MRI scan and tau PET scan, and the preprocessing steps are explained in the companion document available through the ADNI website.

The Berkeley preprocessing pipeline is as follows: each patient T1 MRI was segmented using FreeSufer 7.1.1; each subject’s normalized-intensity flortaucipir scan was coregistered to their bias-corrected T1 image; partial volume effects were corrected for using the geometric transfer matrix approach (Baker, Maass, & Jagust, 2017; Baker et al., 2017); and mean flortaucipir uptake is reported, in the dataset, for each region of the FreeSurfer segmentation. We further normalized each subject’s regional SUVR score using the inferior cerebellar gray matter as the reference region, as suggested in the dataset documentation.

We proceeded to create two datasets using this preprocessed dataset. For the first set, we computed the volume-weighted SUVR average for each patient record over the whole brain to compute a global score (GS) and for each region of Table 2. Thus, each ADNI patient visit record was assigned a global volume-averaged SUVR score along with five volume-averaged regional scores corresponding to the regions of Table 2. We refer to this dataset as the ‘base dataset’.

We then created a second dataset we refer to as the ‘partitioned dataset’ by following the methodology of Schöll et al. (2016). An overview of the creation process follows: First, the GS was used as a decision variable, in a conditional inference tree (CIT) algorithm. Second, the CIT algorithm was used to compare the GS to the BraakV regional volume-weighted SUVR, producing a GS threshold. This GS threshold was used to separate patient records into two groups who differed significantly based on their BraakV score. The BraakV patient records were then moved, out of the full dataset, into a new data file. This process was then repeated, once again applying the CIT algorithm and using the GS as a decision variable, but now comparing to the BraakIV regional volume-weight SUVR score, thus creating a BraakIV group of patient records, and so on. This process resulted in five separate files, that is, the partitioned dataset, corresponding to a patient record classification according to the regions of Table 2. The thresholds produced by the CIT algorithm can be seen in Table 6. Our application of CITs follows the identical approach of previous authors (Schöll et al., 2016), in a similar SUVR study, to categorize each patient record into the discrete groups corresponding to a predefined set of regional stages.

Our next goal was to see if ADNI flortaucipir data may suggest alternatives to the progressive I → II → III → IV → V (computational) Braak-like sequence, in particular for NFT due to the binding of flortaucipir to paired helical filaments. Our objective here is not to rigorously study hierarchical SUVR behavior as in Cho et al. (2016). Rather, our goal is to conduct a simple inquiry as to whether the ADNI data may suggest alternative staging patterns of potential interest for the computational staging problem and model selection, especially as it pertains to observed NFT staging.

Our first investigation used the partitioned dataset. First, we checked that the CIT algorithm separated the data in such a way that all of regional mean SUVR values, across all possible group-region pairings, were significantly (p < 0.05) different, using a Welch’s t test; this was the case, as expected by an application of CIT-based separation. Next, we computed the Pearson correlation between each stage group’s primary stage and all other stages within that group; the results are reported in Table 7.

Assuming that overt NFT pathology originates in region I, Table 7 suggests that staging sequences with I → III may also be of interest for NFT. Assuming that III proceeds I, the next suggested progression would be III → II. This observation suggests that sequences starting with the prefix I → III → II could be of interest alongside the more progressive stage prefix I → II → III. Finally, likely due to the sparsity of patient records assigned to Group V, none of the correlations of Group V were significant.

Our next investigation uses the base dataset. The primary motivation was to ensure that patient gender and age were not confounding our partitioning-based strategy. Starting with the base dataset, we cross-referenced ADNI patient data and augmented the base set with the age, at the time of scan, in addition to the patient’s gender. We then used regression models covarying for age and gender. We compared each regional volume-averaged SUVR score to that of the other regions; the results are shown in Table 8. Evidence for influence was evaluated using standardized factor scores when the influence was significant. Region III was most affected by Region I (R² = 0.775, B = 0.613, p < 0.001), further supporting interest in the I → III staging prefix. Region II was also most influenced by Region I (R² = 0.486, B = 0.486, p < 0.001) but at a lower value of both B and R². Removing Region I from the model decreased the variational fidelity (R²) of the Region III model by 24.8% and the Region II model fidelity by 13.6%; in the latter case, the Region II score was most affected by the Region III score (R² = 0.42, B = 0.5). Supposing that pathology begins in the entorhinal cortex (Region I), these observations suggest, again, that a I → III → II prefix is also of interest for evaluating computational staging.

Next, we take a second look at the staging suffix options IV → V versus V → IV. To do so, we regressed the Region IV and Region V scores in terms of the first three regional scores, once more covarying for age and gender. The Region IV model identified Region III as the most influential (R² = 0.43, B = 0.75, p < 0.001). The Region V model was only marginally explained by the data (R² = 0.22), but the most influential factor was, again, the Region III score (B = 0.38, p < 0.001), but with a substantially lower influence than in the Region IV model. These latter observations suggest that a IV → V may be more desirable than a V → IV suffix. We have therefore corroborated, and extended, the conclusions suggested by the partitioning strategy (Table 7).

Next, we corroborate our previous SUVR findings by following a similar approach to Cho et al. (2016) and considering z-scores, computed from the base dataset, for the regions defined in Table 2. First, all of the regional mean SUVR scores were individually standardized. Similar to Cho et al. (2016), we considered a particular region to be ‘involved’ in a fixed patient record if that region’s Z score satisfied Z ≥ 2. We then selected all patient records with an involved Region I (entorhinal cortex). Continuing within this subset of data, we counted the frequency of all other involved regions and computed their mean z score. The results are reported in Table 9. Once more, we see that progression sequences starting with I → III → II may also be of computational interest and, further, that the staging suffix IV → V appears more preferable, for SUVR, than the V → IV alternative.

Finally, we compared the ADNI data findings, above, to the results for τP seeds reported in DeVos et al. (2018) where Braak II brains (classified postmortem) also showed a slight average increase of τP seeds in physiological region III (i.e., region III, Table 5) over that of physiological region II; further suggesting that a computational I → III → II pattern may be of interest for τP seeds. If one considers post-mortem Braak classification as a pseudo-time then, from DeVos et al. (2018, Fig. 3), it is tempting to suggest a computational I → III → II → IV → V pattern of interest for τP seeds. However, the differential between seeding levels in the latter stages is slight. In summary, ADNI SUVR data, and results from DeVos et al. (2018), suggest a broader set, beyond the choice of a strictly ascending progression, of potentially interesting computational staging patterns for the choice of regions in Table 2. These staging patterns are reported in Table 10. We note that either more late stage data (i.e., Group V, Table 6), a finer SUVR staging sequence, such as that of Cho et al. (2016), or both is needed to ascertain a higher degree of certainty in the latter τP computational staging patterns for AD.

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. For up-to-date information, see www.adni-info.org. A complete listing of ADNI investigators can be found at: https://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00208.

Prama Putra: Conceptualization; Formal analysis; Investigation; Software; Visualization; Writing – review & editing. Travis B. Thompson: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Software; Supervision; Visualization; Writing – original draft; Writing – review & editing. Alain Goriely: Conceptualization; Formal analysis; Funding acquisition; Methodology; Project administration; Resources; Software; Supervision; Writing – original draft; Writing – review & editing. Pavanjit Chaggar: Data curation; Formal analysis; Software; Writing – review & editing.

Alain Goriely, Engineering and Physical Sciences Research Council (https://dx.doi.org/10.13039/501100000266), Award ID: EP/R020205/1. Travis B. Thompson, John Fell Fund, University of Oxford (https://dx.doi.org/10.13039/501100004789), Award ID: 000872. Pavanjit Chaggar, Engineering and Physical Sciences Research Council (https://dx.doi.org/10.13039/501100000266), Award ID: EP/L016044/1.

Prion-like:

A prion-like biological agent is one that replicates an undesirable state through an autocatalytic process.

Nonlinearity:

A nonlinearity is any mathematical expression that does not satisfy f (ax + y) = af (x) + f (y) where a is a real number and x and y are input vectors.

Simplicial mesh:

A, often complex, mesh whose elements are triangles (2D) or tetrahedra (3D); widely used to solve partial differential equations on complex domains.

Model selection:

The process of arriving at the set of features that determine the model used to study a process of interest.

Parameter inference:

Typically refers to mathematical approaches to deduce distributions for the parameters of a model, from observed data, using statistical paradigms such as Bayesian inference.

Structural connectome:

A structural connectome is a graph representation G = (V, E) of the brain where the vertices V are anatomical regions of interest and the edges represent connections between regions of interest via axonal bundles.

Graph Laplacian:

The graph Laplacian is a matrix representation of a transport process on a, often weighted and undirected, network graph.

Diffusion-reaction:

A partial differential equation describing the simultaneous diffusion and reactive production of an agent through a medium.

Fisher-Kolmogorov:

A class of diffusion-reaction equations, that admit traveling wave solutions, where the reaction term satisfies a set of conditions. This class includes Fisher’s equation as a special case.