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

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 PrP^C) into an infectious agent protein (usually denoted by PrP^Sc). 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 PrP^Sc by spontaneous conversion of PrP^C. The PrP^Sc protein then catalyzes the conformational change of PrP^C by forming a heterodimer, which in turn, dissociates at faster rate into two molecules of PrP^Sc. 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 PrP^Sc would determine whether the system maintain the healthy PrP^C state (i.e. stable concentration of PrP^C) or turn (by an autocatalysis mechanism) into a state of exponential production of PrP^Sc. 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 PrP^Sc (26). This behavior clearly contrasts with the observations of fibrils of pathogenic PrP^C 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 PrP^Sc aggregates can either elongate by the incorporation of a PrP^C monomer or break into two new aggregates. These mechanisms of elongation and fragmentation of PrP^Sc 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).

A detailed analysis [see Ref. (28)] of the Eq. A1 revealed a mechanism whereby aggregation could be initiated from a PrP^Sc 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 PrP^Sc aggregates below a critical threshold size would become unstable and quickly dissociate into pieces that return to a PrP^C 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 PrP^Sc 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 R₀ as the number of secondary infections produced on average by one infectious prion (i.e., PrP^Sc). Thus, if R₀ ≤ 1, the prion replication stops and the disease-fee equilibrium state turns out globally asymptotically stable. On the contrary, if R₀ > 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 PrP^Sc in the peripherical nervous system and its appearance in the central nervous system. This latency cannot be explained by the rate at which PrP^C is converted into PrP^Sc 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 PrP^Sc 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 template-assisted 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 PrP^Sc 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.

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 multi-step 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 nucleation-dependent 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 early-times 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 R₀ 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 steady-state level if the polymerization ratio R₀ 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).