Mathematical Biology Models of Parkinson's Disease.

Parkinsons disease (PD) is a progressive neurodegenerative disease with substantial and growing socio-economic burden. In this multifactorial disease, aging, environmental, and genetic factors contribute to neurodegeneration and dopamine (DA) deficiency in the brain. Treatments aimed at DA restoration provide symptomatic relief, however, no disease modifying treatments are available, and PD remains incurable to date. Mathematical modeling can help understand such complex multifactorial neurological diseases. We review mathematical modeling efforts in PD with a focus on mechanistic models of pathogenic processes. We consider models of α-synuclein (Asyn) aggregation, feedbacks among Asyn, DA, and mitochondria and proteolytic systems, as well as pathology propagation through the brain. We hope that critical understanding of existing literature will pave the way to the development of quantitative systems pharmacology models to aid PD drug discovery and development.

Parkinsons disease (PD) is the second most‐common progressive neurodegenerative disease (ND) after Alzheimers disease (AD) and affects around seven million people worldwide. The incidence of PD is expected to increase with life expectancy. PD is characterized by tremors and movement rigidity. Neuronal degeneration is observed in a brain region called the substantia nigra pars compacta (SNc). This region contains dopaminergic neurons and their loss results in reduced dopamine (DA) in the striatum, which is responsible for the motor symptoms of PD. Current therapeutic interventions focus on restoring DA levels either through direct administration of a DA precursor (such as L‐dopa) or blocking of DA degrading enzymes (e.g., monoamine oxidase blockers). DA receptor agonists are also used to functionally compensate for loss of DA. Although these treatments have been successful in achieving symptomatic relief in PD, they are not disease modifying and, hence, PD remains incurable.

PD is a complex multifactorial disease resulting from aging, genetic predisposition, and exposure to environmental toxins. Figure 1 represents the current understanding of the complex interaction network associated with PD pathogenesis. We have used recent reviews1, 2, 3, 4, 5 to construct this network. This is not an exhaustive network because we have restricted it to broad themes for clarity. An exhaustive map of PD (PDMap) has been published elsewhere.1

Figure 1

Interaction network of various pathways involved in pharmacodynamic pathogenesis. The network is generated by referring to recent reviews in the field.1, 2, 3, 4, 5 Molecular species are shown in green ovals whereas molecular/cellular processes are shown in yellow rectangles. Positive and negative interactions are identified using sharp and blunt arrows, respectively. Double‐negative feedback motif is identified by red arrows, while double‐positive feedback motifs are identified by blue arrows. Gray arrows and processes shown in lighter shade of yellow indicate interactions that have not been modeled quantitatively. Asyn, α‐synuclein; DA, dopamine; GSH, glutathione; RNS, reactive nitrogen species; ROS, reactive oxygen species. See Table 1 for the list of abbreviations.

Table 1

List of abbreviations

ACh	Acetylcholine
AD	Alzheimer's disease
Asyn	Alpha‐synuclein
BG	Basal ganglia
BST	Biochemical systems theory
CMA	Chaperone‐mediated autophagy
CSF	Cerebrospinal fluid
DA	Dopamine
ECF	Extracellular fluid
FBA	Flux‐balance analysis
GSH	Glutathione
LBs	Lewy bodies
LNs	Lewy neurites
MRI	Magnetic resonance imaging
NDs	Neurodegenerative diseases
NCP	Nucleation conversion polymerization
NP	Nucleation polymerization
ODE	Ordinary differential equations
PK/PD	Pharmacokinetic/pharmacodynamic
PD	Parkinson's disease
PET	Positron emission tomography
QSP	Quantitative systems pharmacology
RNS	Reactive nitrogen species
ROS	Reactive oxygen species
SNc	Substantia nigra pars compacta
TNT	Tunneling nanotubes
UPDRS	Unified Parkinson's Disease Rating Scale
UPP	Ubiquitin proteasome pathway

To date, around 15 genes have been identified with links to PD. Plotegher and coworkers4 recently published a list of these genes with their associated functions. PD pathogenesis involves processes such as aggregation of a protein named α‐synuclein (Asyn), oxidative stress, and dysfunction of proteasomes and lysosomes. Three feedback motifs have been identified for a long time; they all involve the misfolding of Asyn. One of these is the double‐negative feedback interaction between misfolded Asyn and proteasomal/lysosomal machinery (highlighted in red in Figure 1). Although proteolytic mechanisms are responsible for clearing misfolded proteins, misfolded Asyn is known to inhibit proteasomes and parts of lysosomal function.6 Two double‐positive feedback interactions are also highlighted in Figure 1 (in blue). Misfolded Asyn can permeabilize DA‐containing vesicles, leading to increased cytoplasmic DA concentration. DA can associate with native Asyn leading to its misfolding. Misfolded Asyn is known to cause increased mitochondrial damage, which, in turn, increases oxidative stress leading to increased production of reactive oxygen species and reactive nitrogen species (ROS/RNS). Increased ROS/RNS leads to further Asyn misfolding. Even though we have highlighted the shortest‐path feedback interactions here, several longer path interactions can also be identified. For example, increased cytoplasmic DA leads to increased ROS, which can lead to Asyn misfolding, or compromized lysosomal function due to misfolded Asyn can cause defects in mitophagy, which, in turn, leads to increased ROS/RNS and consequently to more misfolded Asyn, or increased neuroinflammation in response to misfolded Asyn leads to increased ROS/RNS, which, in turn, leads to increased Asyn misfolding. Apart from these feedback mechanisms, several other factors are associated with PD. For example, increased concentrations of metal ions such as iron (Fe2+) and copper (Cu2+) in PD brains are known to cause Asyn misfolding and increased ROS.6 Age‐related decline in protein clearance mechanisms and mitochondrial function as well as increase in inflammation are known to affect PD pathogenesis.

Understanding PD requires an interdisciplinary approach involving experimental and modeling studies. Mathematical models of PD have developed concomitantly with accumulation of experimental insight and address several of the mechanistic aspects of PD pathogenesis. A systematic review of modeling efforts in various NDs has recently been published.7 In this review, we focus on various approaches in PD modeling. PD models may be broadly categorized into two classes: (i) mechanistic models and (ii) phenotypic models. The latter class includes models that quantitatively describe some aspect of motor symptoms, such as tremors or gait disturbances or electroencephalography characteristics, using signal processing. The predominant aim of these models is to identify quantitative differences between healthy and diseased subjects. In further modification, neural networks in certain brain structures, such as the basal ganglia (BG) and spinal cord, are modeled to simulate observed PD phenotype. Sarbaz and Pourakbari8 recently published a detailed review of phenotypic models.

In this review, we focus on the former class of models (i.e., mechanistic models of molecular pathways underlying PD). The PD models have collectively considered all molecular species and several interactions depicted in Figure 1. Through an exhaustive literature search, we have found 32 models that specifically model PD. We categorize these models based on the aspect of PD pathology they study, summarize the insights gained, and provide a critical discussion of the models. Where necessary, we review modeling literature that served as a precursor to these PD models and that from related NDs. Last, we conclude with identifying further modeling opportunities. Because the focus of this article is on mechanistic PD models, we have not extensively reviewed the developments in the disease progression and drug exposure‐response models in the field. We briefly refer to these modeling efforts toward the end of the article.

Mechanistic Models in PD

Mechanistic models of PD are classified into (i) Asyn aggregation models, (ii) pathogenesis models, and (iii) pathology propagation models (Figure 2). Of these, the Asyn aggregation models are the most parsimonious and have the greatest experimental support. Pathogenesis models range in size from a few variables to several hundred variables. Pathology propagation has received much less attention from the modeling community.

Figure 2

Classification of 32 pharmacodynamic (PD) models based on aspect of PD pathology addressed by the models. Further description of each category of models can be found in text. Figure also shows the number of models in each category. CMA, chaperone‐mediated autophagy; DA, dopamine; NCP, nucleation conversion polymerization; NP, nucleation polymerization; ROS, reactive oxygen species; UPP, ubiquitin proteasome pathway. See Table 1 for the list of abbreviations.

Asyn aggregation models

One of the first genes to be associated with familial PD was SNCA — the gene coding for Asyn. Asyn is a 14 kDa protein, which is relatively unfolded in its native state. However, it can assume a misfolded conformation rich in beta‐sheets. This misfolded protein is capable of aggregation into lower or higher molecular weight oligomers or be precipitated into insoluble fibrils or large inclusions. Postmortem studies of PD brains display characteristic histopathological features called Lewy bodies (LBs) and Lewy neurites (LNs). Asyn aggregates are a major component of LBs/LNs.2

Several other molecules, such as amyloid‐beta and tau proteins as well as microtubules and actin, share this property of Asyn. Aggregation is studied in vitro using fluorescence spectroscopy techniques, in which a fluorescence marker responds to changes in protein secondary structure (for example, from relatively unstructured to beta‐sheet core), and can be quantitatively correlated with appearance of oligomers or aggregates. Such experiments can then be extended to study effects of mutations or other (bio)chemicals on the aggregation kinetics. Several in vitro experiments using purified “aggregation‐prone” proteins show a sigmoidal time course of aggregates formation, as shown in Figure 3. The lag‐phase in the scale of seconds to minutes is associated with slow nucleation process, which generates “seeds,” whereas the growth phase indicates rapid elongation of seeds. The eventual plateau may indicate exhaustion of monomeric species or attainment of equilibrium.

Figure 3

Schematic representation of molecular aggregation pathway and the corresponding kinetic observation of appearance of aggregates over time.

Mathematical models of polymerization of biomolecules were developed in early 1970s by Oosawa and Kasai and Oosawa and Asakura.9, 10 These self‐assembly reactions are believed to begin with a nucleation event, in which native monomers are altered to form “seeds” for elongation or polymerization. This class of models is called nucleation polymerization (NP) models. Several extensions of the NP theme have been proposed, such as inclusion of fragmentation and association of fibrils, role of molecular chaperones, as well as additional kinetic steps. Readers are referred to a comprehensive review of protein aggregation models by Morris et al.11 The aim of these models is to extract kinetic parameters corresponding to various processes involved in aggregation. For this reason, these models are often parsimonious in terms of the molecular events modelled. For example, processes such fragmentation of fibrils and generation of seeds from fragmentation may not be included in the model unless experimental data warrant this additional complexity.

A number of groups have studied Asyn aggregation using a combination of modeling and experimental approaches. In the earliest work, Morris and coworkers12 used previously published Asyn aggregation data and fitted them with a two‐step NP model known as the Finke–Watzky model, and studied the effect of temperature and pH on aggregation kinetics.13 More recent experiments with Asyn have revealed that two kinds of Asyn oligomers exist—ones that are formed initially and are unstable, and ones that form slowly and are stable to proteinase‐K treatment and may form the substrate for fibril formation.14, 15 Thus, the NP model was extended to include an additional “conversion” step between nucleation and polymerization. The resulting model was called nucleation‐conversion‐polymerization (NCP) model and included monomers, two types of oligomers, and fibrils.16, 17, 18 Several authors have used experimental studies and kinetic modeling to investigate the effects of metal ions,19 lipid membranes,17 and mutations20 on Asyn aggregation kinetics.

Typically, analytical solutions of ordinary differential equation (ODE)‐based NP models are used for data fitting. Finding analytical solutions to basic and advanced NP models has been an active area of research.16, 21, 22, 23 Furthermore, one group has developed a free web‐based global data fitting algorithm for fitting molecular aggregation data.24 The algorithm is based primarily on NP models and allows for a certain finite variations on the NP theme. As of writing of this manuscript, their algorithm is available at http://www.amylofit.ch.cam.ac.uk/login.

All aggregation models are based on kinetic data generated from in vitro experiments. Such experiments use high concentrations of Asyn monomers (from 0.5−200 μM), presumably to achieve measurable aggregation in experimentally viable time frames. These Asyn concentrations are substantially greater than in vivo Asyn concentrations. For example, Asyn concentrations in human brain homogenates from PD patients are around 10 nM (and around 6 nM for healthy subjects).25 Asyn concentrations in cerebrospinal fluid (CSF) and plasma are even lower.26, 27 It is not clear if Asyn aggregation kinetics differs qualitatively at low concentrations. Note that under experimental settings monomeric Asyn suffers substantial depletion within an hour (for example, see figure 3 b in ref. 18). This is not expected in vivo. Moreover, DA and ROS are known to cause Asyn misfolding. Their effect on the kinetics of aggregation has not been studied yet. Last, dissociation reactions are largely ignored in Asyn aggregation models. Although this is justified on experimental time scales, dissociation may play a role in vivo at larger time scales associated with disease progression.

Pathogenesis models

PD pathology is believed to result from complex interplay and feedbacks among Asyn, ROS, protein degradation machinery, and DA (Figure 1). PD pathogenesis models focus on some or all of the above interactions in various levels of detail. Figure 4 shows a Venn diagram in which the models are categorized based on which cellular processes they consider.

Figure 4

Venn diagram showing classification of pathogenesis models based on the cellular processes they model. The references for the models that used a description of α‐synuclein aggregation are colored in red. CMA, chaperone‐mediated autophagy; DA, dopamine; ROS, reactive oxygen species; UPP, ubiquitin proteasome pathway.

Along with the diversity in cellular pathways modeled, pathogenesis models also differ in the modeling formalisms used. Figure 5 shows classification of models based on modeling formalisms they use. We find that Asyn aggregation models are typically based on ODEs. However, the pathogenesis models use ODEs, network models combined with ODEs, stochastic simulation algorithms, biochemical systems theory (BST; see Glossary), and flux‐balance analysis (FBA). In this section, we use mechanism‐based categories of pathogenesis models and refer to modeling formalisms where necessary.

Figure 5

Classification of pharmacodynamic pathogenesis models based on the modeling formalisms they use. BST, biochemical systems theory; FBA, flux‐balance analysis; ODE, ordinary differential equation.

Reactive oxygen species models

Of the 15 different genes so far associated with PD, more than half lead to mitochondrial defects.4 Aging and environmental toxins also lead to mitochondrial dysfunction.1, 28 Early work in PD modeling has focused on mitochondrial energy metabolism and ROS generation as pathogenesis mechanisms.

One of the earliest ODE models of PD pathogenesis was proposed by Raichur and coworkers.29 Generation of ROS and RNS via various pathways was modeled using minimal descriptions. Apart from ROS, the authors also studied the effect of iron concentration and ubiquitin proteasome pathway (UPP) dysfunction on the levels of Asyn aggregates. Their model was further expanded by adding glutathione (GSH) metabolism and mitochondrial dysfunction. This model predicted that loss of mitochondrial efficiency could affect GSH synthesis via adenosine triphosphate (ATP), which can contribute to increased ROS/RNS, further affecting mitochondrial efficiency. Authors hypothesized that this increased ROS burden may be responsible for neurodegeneration.30

In 2009, Cloutier and coworkers31 presented a detailed model of brain energy metabolism in which they modeled glycolysis and mitochondrial energy metabolism in neurons and astrocytes and compared their predictions with in vivo data from rat brains. In a later work, the authors incorporated the effect of aging in this model by time‐dependent reduction in efficiency of mitochondrial complex 1. They argued that reduction in ATP levels in response to transient high demands might lead to reduced cellular housekeeping. Thus, impaired energy metabolism was proposed as the etiological reason behind accumulation of toxic effects of PD.32 Energy‐related metabolites in mouse brain slices (healthy and PARK2 knockout) were studied using metabolomics and the data were used to train the energy metabolism model. This work showed that the cells can maintain robust control of ATP levels despite the genetic mutation,33 and, as such, reduction in ATP levels may not be the etiological reason behind PD pathology.

ODE‐based brain energy metabolism models lack the spatial details, such as effects of diffusion and locus of synaptic activity, with respect to capillaries. Incorporating these effects may result in predictions that differ spatially, and averaging spatially detailed models with ODEs may only be valid under certain parameterizations.34, 35

A different hypothesis for PD pathogenesis, which involves feedback between ROS and misfolded Asyn, was tested through modifications of the energy metabolism model. Authors included Asyn aggregation and minimal description of its proteolytic clearance together with the aforementioned feedback. This model with 33 variable ODEs displayed a bistable‐switch‐like behavior in response to various factors, such as aging, environmental toxins, and mutations. The lower steady state characterized by low levels of ROS and misfolded Asyn was denoted as the healthy state and the higher steady state with higher levels of ROS and misfolded Asyn was called the diseased state.36 A reduced version of the same model with only two variables (ROS and misfolded Asyn) was able to recapitulate the behavior of the larger model and allowed for analytical insights into the bistability and bifurcation behavior.37 In vivo experiments in which ROS stress response to paraquat (an herbicide known to induce Parkinsonism) was observed in rat brains provided support for the bistable switching behavior, albeit over a short timescale.38 Simulations of the model by Raichur and coworkers29 also hint at the presence of bistability, although the authors did not comment on it at the time.

Ubiquitin proteasome pathway, chaperone‐mediated autophagy, and lysosomal clearance

PD pathogenesis involves a negative feedback between misfolded Asyn and proteasomes (Figure 1). Two genes with function in proteasome, namely Parkin and UCHL1, are known to be mutated in familial PD.6 Recent work has shifted the focus from proteasomes to lysosomes and autophagosomes.4

In this section, we discuss models involving feedback between proteolytic pathways and Asyn aggregation. A theoretical model of three ODEs with minimal description of Asyn aggregation and its interaction with proteasomes was proposed.39 The model predicted bifurcation behavior as a function of the ratio between Asyn fibrils and free proteasomes. Homeostasis could be maintained for lower ratios, but at higher ratios, free proteasome levels oscillated with prolonged periods of low concentration. Authors hypothesized that such prolonged proteasome depletion could lead to accumulation of Asyn oligomers leading to PD.39 Thus, PD pathogenesis was predicted, although the proteasome inhibiting function of Asyn aggregates was not explicitly modeled in this work. A common feature of this model shared with a few previous ROS models is that both types of models predict pathogenesis via the system undergoing a bifurcation.

In one of the more recent works, the Asyn mediated inhibition of proteasomes was explicitly modeled.40 The model included misfolding of DA‐bound Asyn, although minimal description for DA metabolism was used. Protein degradation machinery was included as three pathways—lysosomes, chaperone‐mediated autophagy (CMA), and UPP. The pathways differed in the kinds of Asyn species they degrade, and the kinds of species that irreversibly block their function. No significant detail of individual pathways (such as ubiquitination, etc.) was included. Authors used stochastic simulations of their model to qualitatively test hypotheses/experimental observations regarding effects of Asyn synthesis, DA levels, and increased lysosomal degradation on levels of toxic oligomeric species. This model built on the previous models of UPP, in which aggregation of a generic (nonspecific) protein was included in a model of UPP and lysosomal pathway.41, 42 The parameters for these models were trained on in vitro data collected on a time scale of 7 days and, as such, it is not clear how the pathogenesis progresses over the timescales relevant for the clinical progression of PD. Furthermore, by nature, these studies only address genetic aspects of PD pathogenesis. Given the stochastic nature of the simulations, it would have been interesting to explore whether stochastic fluctuations alone (or in combination with a genetic component) could produce qualitatively different model outputs, thus switching from healthy to diseased states. This kind of analysis is worth exploring with respect to sporadic disease onset.

DA metabolism

Loss of dopaminergic neurons and resulting reduction in DA in brain are hallmarks of PD. Dopaminergic neurons are involved in synthesis, storage, release, and reuptake of DA. Increased cytoplasmic levels of DA can lead to generation of ROS as well as misfolding of Asyn (Figure 1). Qi and coworkers43, 44 built a model of DA metabolism using the BST formalism. Simultaneously, another group of authors used an ODE‐based approach to model DA metabolism.45 Both groups investigated effects of key enzymes or transporters on DA homeostasis. In further work, Qi and coworkers46 also investigated the effect of rotenone and paraquat (two herbicides known to cause Parkinsonism) on DA metabolism.

A model of dopaminergic neuron using FBA considered steady‐state fluxes of several variables including DA, Asyn, ROS, and proteasomal machinery.47 This model showed in a qualitative manner that stressors, such as neurotoxins and increased production of Asyn and DA, can lead to PD pathology. The BST approach with minimal kinetic information was used by Sass and coworkers48 to study the interplay of DA metabolism, Asyn, and proteasomal and lysosomal pathway. Like Büchel et al.,47 the authors showed qualitatively that disruption of modeled cellular pathways results in PD. A curious connection among insulin resistance, inflammation, DA, and Asyn aggregation was modeled in a theoretical work by Braatz and Coleman.49

All the models utilizing the BST or FBA approach contain several variables and processes. Lack of kinetic information regarding these processes has led the authors to use relative species concentrations and parameter values in these models. The model outputs are also typically fold‐changes in variables at steady states (diseased or healthy) rather than absolute concentrations. Despite this limitation, qualitative agreement with experimental data has been demonstrated.43, 48 Smaller ODE‐based models, on the other hand, have used levels of Asyn aggregation or ROS or proteasomes/lysosomes to identify a “PD phenotype.” This is predominantly informed by in vitro experiments in which levels of such markers are studied. Thus, PD phenotype is understood in terms of markers of the pathways that are modeled. However, even in the case of ODE models, the comparison with experimental/clinical observations remains largely qualitative.

We would like to note that only a fraction of pathogenesis models include a submodel with kinetic description of Asyn aggregation (Figure 4, highlighted in red). None of these models use the kinetic parameters derived from Asyn aggregation models discussed previously. A subset of those models, which do not include Asyn aggregation kinetics, do include some kind of the Asyn aggregates (such as oligomers or fibrils) as a model species. In these cases, the feedback between misfolded Asyn and ROS or proteasomes is modeled.37, 39, 49 A few other models focus on non‐Asyn mediated pathogenesis mechanism (in particular the role of DA homeostasis). We feel there is a clear possibility to enrich the pathogenesis models with realistic Asyn aggregation kinetics and parameters derived from aggregation models. It is possible that the individual steps of monomer to oligomer to fibril formation are less impactful for the pathogenesis as compared with the feedback of misfolded Asyn species on relevant systems. Once adequately parameterized, modeling can be used to ascertain the relative contributions of each process and then perform informed model reduction.

Moreover, almost all pathogenesis models use a substantially high native Asyn concentration (100 μM in ref. 36, 37; 100–200 μM in ref. 29; or 2 μM in ref. 40). As noted before, experimentally observed Asyn concentration in the brain is lower by several orders of magnitude. This has implications for models with bifurcations or bistability (e.g., ref. 37), in which the bistable behavior is sensitive to concentrations of native Asyn. Recalibration of relevant parameters in pathogenesis models is necessary to reflect physiological Asyn concentrations. Furthermore, in case of models with bistability, other critical parameters may also need to be revised if bistable behavior is to be preserved at lower Asyn concentrations.

Additional models

Hodgkin and Huxley50 proposed an ODE‐based model of action potentials in squid giant axon. This model has been used as a basis for neuronal action potentials in several neural systems, including the SNc neurons. In 1986, Gerold Porenta51 published a model of neurotransmission using the Hodgkin–Huxley model of membrane potential, with an addition of neurotransmitters such as DA, acetylcholine (ACh), and gamma‐aminobutyric acid. Through analysis of this model, he showed that disease parameters result in a shift in the eigenvalues of the model. Measurements of local field potentials in the brains of patients with PD show typical oscillatory activity. This correlates with bradykinesia. Several groups have modeled these so‐called beta oscillations and proposed different hypotheses for the differences between healthy and diseased states. For example, Francis and coworkers52 modeled neurotransmission using sodium, potassium, and calcium currents. Their model recapitulated the tonic pacemaking oscillations of SNc neurons. They argued that increased energy costs are associated with calcium channels in these neurons, thus making them susceptible to energy stress. Using a neural network of 100 neurons, McCarthy and coworkers53 showed that oscillations arise due to inhibitory interaction of a group of striatal neurons. They are enhanced due to increase in ACh caused by DA depletion in PD. Moran and coworkers54 have modeled membrane potentials and connectivity in a diverse set of neurons, including BG, thalamus, cortex, striatum, and subthalamic nucleus. They argued that chronic DA depletion alters the connectivity in BG and thalamo‐cortical circuits, thereby generating the enhanced beta oscillations. In a more recent article, Pavlides and coworkers55 hypothesize that resonance and feedback between BG neurons and motor cortex is responsible for the oscillations. All the models of beta oscillations are based on the assumption that the neural networks are in a DA‐depleted state. In this sense, they mechanistically model the emergence of symptoms downstream of DA depletion but not the mechanisms of DA depletion itself. In a recent study, Roberts and coworkers56 developed a computer‐based mechanistic model of BG and combined it with clinical data for various pharmacological interventions. Their model provides a mechanistic link among neural anatomy, physiology, and pharmacology, and, therefore, is an example of a quantitative systems pharmacology (QSP) model.

Pathology propagation models

Braak and coworkers57 showed over a decade ago that PD pathology spreads along anatomically connected brain regions, much like prion diseases. Even though longitudinal data have not been gathered, it is expected that this pathological spreading occur over a time‐scale of years. Pathology propagation has been shown in several animal models of PD.58 Recently, exosomes, axonal spread, tunneling nanotubes (TNTs), gliosis (see Glossary), leakage from degenerating neurons and glymphatic flow have all been implicated in the spread of pathology.59 This diversity of mechanisms combined with brain geometry and anatomic connectivity considerations makes it difficult to develop mechanistic propagation models.

Consequently, this aspect of PD has received much less attention from the modeling community in comparison to Asyn aggregation and pathogenesis. To our knowledge, a single group of modelers have studied intra‐axonal transport of Asyn.60, 61 Kuznetsov and Kuznetsov60 modeled Asyn transport as either a purely advective process (Asyn is transported by molecular motors as part of multiprotein complexes) or as a combination of advective and passive diffusive processes. They adjusted boundary conditions to construct healthy and diseased neurons and predicted that accumulation of Asyn in axon terminals could lead to formation of aggregates. In a further article,61 they divided a neuron into two subcellular compartments, namely soma and synapse. They collapsed the diffusive/advective processes within an axon into a single “transport” term for native Asyn, added minimal descriptions of Asyn aggregation (Finke–Watzky model) and degradation, and tested hypotheses regarding likely causes of Asyn aggregation in PD.

The aim of Kuznetsov and coworkers62 was to study PD pathogenesis. Thus, despite being the only work to consider Asyn transport, the authors have not explored the potential of their models to simulate pathology propagation according to the Braak hypothesis. Arguably, one‐dimensional axonal Asyn transport (which may be an oversimplification of propagation in the complex brain geometry) may be adequate to study pathology propagation from the enteroneural system to the brain and may be useful to explore gut‐to‐brain propagation of PD.

The spread of PD pathology due to misfolded Asyn has led to the characterization of PD as a prion disease.63 A few alternative approaches have been used to model propagation of prion‐like proteins in other NDs and prion diseases. Most notably, these involve brain magnetic resonance imaging (MRI) combined with network approaches. Raj and coworkers64 used MRI data to construct neuronal connectivity tracts. They then allowed a generalized “disease factor” (which could be any prion‐like protein) to diffuse through this connected neural network. This approach allowed the authors to correctly predict brain regions most affected by atrophy in AD.64, 65 In parallel, Zhou and coworkers used functional MRI to study connectivity patterns of brain regions most affected by selected NDs and showed that these patterns can give rise to network “epicenters,” which are most susceptible to atrophy associated with those NDs.66 Nowak and coworkers67 used a two‐dimensional array with random distribution of cells in silico and allowed prion pathology propagation using simple epidemiologic sensitive and infected cell models. Their work showed that the propagation probability increases with concentration of prion proteins per cell. In a recent theoretical article, authors used a model of a neural network together with processes of protein synthesis, misfolding, and transport. Transport of misfolded proteins was modeled as axonal active transport (anterograde and retrograde), trans‐synaptic transport, and passive diffusion. Authors showed that origination of seed and diffusive spread were crucial in determining pathogenic propagation patterns.68 Yau and coworkers69 studied cortical thinning in PD patients and showed that the atrophy pattern correlated with neuronal connectivity to a “disease reservoir,” adding support to the Braak hypothesis. Almost all of these models are agnostic with respect to the pathologic protein they describe. The characteristic pathology propagation strongly depends upon the network connectivity. In this respect, it is worth noting that although the prion‐like behavior of Asyn is widely accepted, its culmination into PD pathogenesis is believed to require further cell‐specific factors, such as high bioenergetic demands of highly branched neurons.70 Thus, combination of network‐diffusion approaches with cell‐specific pathophysiological models may be necessary to recapitulate spread of PD pathology.

Conclusion and Future Directions

In this article, we have presented a historic review of disease modeling literature in PD. A wealth of literature has provided robust kinetic understanding of the Asyn aggregation process. Models of PD pathogenesis have been developed with diverse modeling formalisms and several etiological mechanisms from experimental literature. These interdisciplinary approaches have proved useful in generating testable hypotheses and providing insights into PD pathogenesis. Models with diffusion in neural networks have been developed for AD and other NDs, and have implications for PD research.

The exact neurodegenerative mechanisms in PD are not yet clearly understood. Several possibilities exist, such as ROS‐toxicity, neuro‐inflammatory toxicity, or toxicity from misfolded protein burden. However, these cellular insults have not been quantitatively linked to cytotoxicity. Consequently, quantitative mechanistic models do not link “diseased” states of these variables or processes with cytotoxicity. Instead, a diseased state is identified by increased burden of misfolded Asyn. Quantitative comparisons of measurable variables (such as misfolded Asyn or ROS) in models with experimental models of PD are necessary to calibrate these models.

Several feedback loops have been identified in PD pathogenesis (Figure 1). The pathogenesis models typically only include one such feedback system (e.g., ROS‐Asyn feedback or Asyn‐UPP feedback). It would be of interest to build a dynamic model to integrate the multiple feedback loops to rightly capture the multifactorial nature of PD, provided it can be adequately parameterized. Such integrated models would be extremely useful for understanding the relative importance of various components involved in the pathogenesis for comparing effectiveness of different potential interventions, and for exploring the synergistic value of combination therapies.

It has been proposed several times that Asyn oligomers are the toxic species, whereas insoluble aggregates, such as fibrils and LBs/LNs, are nontoxic forms. Formation of LBs/LNs may, therefore, be neuroprotective as a mechanism for oligomer clearance, although it is hard to imagine that the neuronal existence of LBs/LNs is not detrimental. Cells may have a dichotomous response to Asyn misfolding, with certain cells forming LBs/LNs, whereas other cells are involved in vigorous clearance of misfolded Asyn, thereby preventing formation of large aggregates. No models so far have explicitly studied the effect of proteolytic parameters on propensity to form LBs/LNs. Given the small numbers of LBs/LNs per cell, stochastic modeling and noise‐induced bistability may be explored in this context.71

Indeed, several mechanistic models of PD have shown existence of a bifurcation. In combination with stochastic switching, this is an attractive mechanism to explain sporadic onset of multifactorial diseases. It is imaginable that multiple factors, including environmental toxin exposure, age, and genetics, may predispose individuals to PD. This may occur by tuning the system parameters so as to bring the system close to a bifurcation point. Noise in the system could then induce spontaneous stochastic switching from a healthy to diseased state, leading to apparently idiopathic onset.

Traditionally, clinical studies have used a scoring system called the Unified Parkinson's Disease Rating Scale (UPDRS) to monitor the progression of PD in the presence or absence of a particular therapeutic intervention. The UPDRS is used to give a score based on self‐reported as well as clinically observed aspects of cognitive and motor functions. Simplistic linear or logistic regression models, together with a saturating drug effect, have been proposed to simulate the time course of UPDRS scores (see review ref. 72). However, it may be difficult to establish a meaningful link between mechanistic PD models and clinical UPDRS scores. One example of such a link is the work by Roberts and coworkers.56 Alternatively, one may need to use biomarkers or clinical end points, which can be modeled mechanistically to build a QSP model. To this end, blood‐based biomarkers73 or brain imaging using positron emission tomography (PET) to study dopaminergic integrity74, 75 are suitable candidates.

It has long been recognized that different aspects of disease symptoms captured by UPDRS progress at different rates.76 Recent application of item response theory to modeling UPDRS77, 78 allows item‐specific analysis in an integrated fashion. This approach has the potential for developing the link between biomarkers and clinical symptoms, as well as for identifying “pharmacological fingerprints” for drug‐class‐specific effects on the whole spectrum of the symptoms and their progression rates. Through system perturbation, such knowledge may eventually help elucidate the link between disease biology and various aspects of the clinical manifestation. A multiscale model of PD, which incorporates molecular level and small timescale changes to neuronal connectivity to evolution of clinical UPDRS scores, would require an approach that integrates diverse modeling formalisms. It is conceivable that a combination of differential equation‐based molecular level models and an agent‐based cell‐level model and network‐based neural connectivity models may be necessary. Building such multiscale models with hybrid modeling formalisms is challenging and presents exciting methodological opportunities. Furthermore, the modeling effort would need to be complemented with dedicated gathering of relevant data in clinical or preclinical settings. A recent analysis of PD markers from the Parkinson's Progression Markers Initiative study has identified several genetic and nongenetic markers that predict the rate of motor function decline.79 Such analyses can enrich multiscale modeling efforts and vice‐versa.

Last, newer therapeutic approaches in PD are directed toward inhibiting Asyn aggregation directly (small or large molecules against Asyn) or indirectly (iron chelation, or targeting products of other PD‐risk genes, such as glucocerebrosidase).80, 81 Indeed, a few ongoing clinical trials (e.g., for monoclonal antibodies against Asyn — MEDI1341 and BIIB054) are using plasma/CSF Asyn levels or brain PET scans as outcome measures. Few other disease‐modifying clinical trials with Asyn antibodies or inhibitors of LRRK2 (a lysosomal function gene that is mutated in PD) are ongoing. Mechanistic modeling approaches presented in this article may be useful to assist such clinical trials. QSP models are useful to compare the effects of these different strategies on desired biomarkers of the disease.82 PD pathogenesis models could be expanded to test the likelihood of success of various potential treatments.

Glossary

Tunneling nanotube

TNTs are membrane‐covered conduits between neighboring cells. They can range in diameter from 50−700 nm and can transport small molecules to large organelles like mitochondria depending upon their size.

Gliosis

Gliosis involves proliferation and migration of various glial cells, such as astrocytes, microglia, and oligodendrocytes. Gliosis occurs in response to central nervous system injury and has a role in neuroprotection and inflammation.

Biochemical systems theory

Biochemical systems theory is an ODE‐based modeling formalism in which reactions are represented as power‐law expansions of reactant species. Under the BST framework, the dynamics of the ith species X is given as

where j represents one of the biochemical processes affecting the dynamics, μ denotes the stoichiometric coefficient, γ the rate constant, and f represents the kinetic order. The BST approach allows greater flexibility in capturing nonlinearity of biochemical network because it allows use of noninteger and negative kinetic orders. The solutions are accurate at an “operating point” (for example, a steady state) and are approximated by linearization in the vicinity of the operating point.83 The method requires estimating a large number of parameters and has limited capability to simulate dynamic response.

Flux‐balance analysis

FBA is a type of metabolic control analysis used to calculate steady‐state metabolic fluxes through large‐scale metabolic networks. Steady‐state assumption, together with optimality (e.g., conservation of resources) leads to a linear system, which requires very few parameters and has low computational cost.

Funding

This work was funded by GSK Early talent postdoctoral scheme.

Conflict of Interest

As Editor‐in‐Chief for CPT: Pharmacometrics & Systems Pharmacology, Piet van der Graaf was not involved in the review or decision process for this article. The authors declared no competing interests for this work.