A continuous-time mathematical model and discrete approximations for the aggregation of β-Amyloid

Alzheimer's disease is a degenerative disorder characterized by the loss of synapses and neurons from the brain, as well as the accumulation of amyloid-based neuritic plaques. While it remains a matter of contention whether β-amyloid causes the neurodegeneration, β-amyloid aggregation is associated with the disease progression. Therefore, gaining a clearer understanding of this aggregation may help to better understand the disease. We develop a continuous-time model for β-amyloid aggregation using concepts from chemical kinetics and population dynamics. We show the model conserves mass and establish conditions for the existence and stability of equilibria. We also develop two discrete-time approximations to the model that are dynamically consistent. We show numerically that the continuous-time model produces sigmoidal growth, while the discrete-time approximations may exhibit oscillatory dynamics. Finally, sensitivity analysis reveals that aggregate concentration is most sensitive to parameters involved in monomer production and nucleation, suggesting the need for good estimates of such parameters.


Introduction
Amyloid-β (Aβ) is a peptide generated by the proteolytic cleavage of the amyloid precursor protein (APP) by the action of βand γ -secretases [10,28]. Aβ deposits in senile and neuritic plaques and hyperphosphorylated tau proteins in neurofibrillary tangles (NFT) are extracellular and intracellular expressions, respectively, of the Alzheimer's disease (AD) neuropathological phenotype, together with selective neuronal loss in hippocampal and neocortical regions. The accumulation of Aβ plaque in the brain is considered as one of the primary (postmortem) diagnostic criterion of AD [29,58,68]. Aβ is a 39-to 43-residue proteolytic product of a parental amyloid precursor protein (APP) that localizes to the plasma membrane, trans-Golgi network, endoplasmic reticulum (ER) and endosomal, lysosomal and mitochondrial membranes [69]. It is produced in the brain throughout life and it accumulates in the cerebral cortex in the elderly and to an excessive degree in AD, but its role in the etiology of Alzheimer's is unproven [60].
When Aβ is initially cleaved, it forms single pieces or monomers. The monomers are then thought to combine into clusters to form Aβ oligomers, and finally the oligomers form insoluble fibrils. This process is collectively known as aggregation. These fibrils further accumulate and become deposited as plaques. It is not clear how long this entire process takes, but some researchers hypothesize that it occurs over many years. Aβ can be produced by numerous types of cells such as neurons, astrocytes, neuroblastoma cells, hepatoma cells, fibroblasts and platelets, suggesting, along with its conserved sequence among different species, that this peptide should have an important function in normal cell development and maintenance [48].
One theory stipulates that the neurodegenerative effects of AD arise from Aβ. This is commonly known as the amyloid hypothesis and is one of the main models of AD pathogenesis [8,30,31,69,71,72]. For the case of Aβ linkage to AD, it has been found that Aβ becomes toxic once aggregated [32,64,70,74,78]. Moreover, numerous studies show a strong correlation between soluble Aβ oligomer levels and the extent of synaptic loss [13,14,24,28,41,46], further suggesting that the soluble oligomers are the causative agents of AD [47,49,76,77].
Despite mounting support from biochemical, genetic and transgenic animal studies in favour of the amyloid hypothesis, the hypothesis itself remains controversial [27,35,36,54,79]. Counter arguments to the amyloid hypothesis have been presented by many researchers. One of the central challenges is the fact that multiple clinical trials of anti-Aβ drugs have failed to reduce the symptoms of AD [40]. In fact, various immunotherapies targeting Aβ in AD model mice were effective in decreasing Aβ deposition in the brain, but it did not lead to the improvement of actual symptoms [57]. A more profound challenge to the amyloid hypothesis was put forth by Hyong-gon Lee et al. [38]. The authors suggest that Aβ, far from being the harbinger of disease, actually occurs secondary to more fundamental pathological changes and may even play a protective role in the diseased brain. Other researchers have put forward the hypothesis that the main factor underlying the progression of AD is the aggregation of tau protein and not Aβ [39].
Whether Aβ is the causal factor of AD or a downstream response to some as yet unidentified causative agent, the association of Aβ aggregation with AD means that understanding this process is of considerable importance. Consequently, in this paper we focus on modelling the process of Aβ aggregation. In its simplest forms, Aβ plaque formation can be described by protein aggregation, involving the misfolding of Aβ into soluble and insoluble assemblies [43,80]. Kinetic studies have suggested that the misfolding of monomeric Aβ has been shown to precede the formation of oligomers, which then serve as seeds for accelerated fibril growth [61], as illustrated in Figure 1. The two phases of Aβ aggregation are: (i) the nucleation phase, in which monomers undergo misfolding and associate to form oligomeric nuclei, and (ii) the elongation phase, in which the oligomeric nuclei rapidly grow by further addition of monomers, resulting in the formation of larger fibrils. The nucleation phase occurs gradually and at a slower rate than the elongation phase which proceeds it since elongation is more thermodynamically favourable [43]. A sigmoidal curve can thus describe this process.
Mathematical models of Aβ kinetics provide a clearer mechanistic understanding of amyloid fibril growth. A number of mathematical models, both purely theoretical and experimentally driven, have been developed to describe the aggregation of Aβ. These include kinetic (ODE) models that consider the concentration of monomers and fibrils Figure 1. The two phases of Aβ aggregation: from monomers to fibrils, adapted from [43]. We can expect to see a lag phase during the nucleation process followed by a phase of rapid growth (green curve). The addition of more oligomers (seeds) speeds up the process and induces faster aggregate formation (blue curve). In contrast, the lack of oligomers introduces lag time and slows down the aggregation process (red curve).
of length j [16,42], more complicated Smoluchowski-type (PDE) systems that incorporate the transportation and diffusion of fibrils [5,11,15,23,25,33,66], models that consider the impact of (hypothetical) treatment strategies [18,19,53,65], systems models to consider the role of Aβ in inflammation and neuron death [67], as well more experimentally driven models for AD including [52,59,63]. Previously, the authors [20] developed a discrete-time model for the aggregation of Aβ from monomers to diamers, . . . up to pre-oligomers. This model was developed under the assumptions that aggregation occurs due to the addition of monomers, the aggregation process is not reversible and oligomers are formed from exactly six monomers. Hence they studied a discrete-time model with five states, i.e. a five-dimensional nonlinear discrete dynamical system.
In this paper, we develop a more general model that extends the model developed in [20] in several ways. First, we incorporate oligomer and fibril stages into the model. These additional stages allow us to model the two phases in the aggregation process: the slow nucleation stage that results in the formation of oligomers and the faster elongation stage that converts oligomers into fibrils. We note that, unlike the traditional kinetic equations (also called nucleation-polymerization models) that only model fibrils of length j, the stage-structured approach used here allows us to distinguish between the oligomer and fibril stage. This approach was motivated by studies that have shown kinetic models may not be best suited to reproduce the aggregation observed in laboratory studies. For instance, it was demonstrated for tau aggregation that the kinetic equations underestimate average fibril length and are inadequate at capturing oligomer concentrations over time. Meanwhile, a model that accounts for the conversion of oligomers to fibrils via a reaction of first order in monomer concentration was shown to provide a better fit [73]. Relatively few models for Aβ aggregation include this distinction between the two stages. Exceptions include [20,39] which considered more complicated Smoluchowski-type (PDE) systems. Second, since various estimates have been given for average oligomer size (such as 6 [26], 8 [34], 2 or more [51]), here we do not assume that oligomers are formed by exactly six monomers. Instead, we leave the model flexible by assuming that oligomers are made of at least n + 1 monomers. This results in a system of differential equations with n + 2 states. All together, these assumptions mean that the model is adaptable and may be modified for different types of experimental data.
The paper is organized as follows. In Section 2, we develop a continuous-time model for Aβ aggregation. This model includes n equations for monomers and intermediate aggregates up to size n (where we use size to mean the number of monomers contained in an aggregate), as well as two equations for the oligomer and fibril stages. In Section 3, we provide a detailed analysis of the existence and stability of the equilibria of this model. Our main tool for locating the eigenvalues of the Jacobian matrix of the system is Gershogrin's Theorem [22,62]. Finally, in Section 4, we provide some numerical simulations to further study the continuous-time model as well as two discrete-time approximations which preserve the local dynamics. A summary of our results is provided in the Conclusion section.

The construction of the continuous-time model
To model the aggregation process of Aβ, we choose to model the slower nucleation phase as size structured, that is we consider the concentration of monomers M 1 , diamers (2 monomers) M 2 , triamers (3 monomers) M 3 , . . . , etc. Meanwhile, we model the oligomer and fibril states as average stages O and P, respectively. We assume that oligomers are made of at least n + 1 monomers, resulting in a system containing n + 2 stages: the monomer stage M 1 , the intermediate aggregate stages M 2 , . . . M n , the oligomer stage O, and the fibril stage P. We assume that aggregation (both nucleation and elongation) occurs via monomer addition. For simplicity, we do not consider disassociation (i.e., monomer removal from aggregates) or fibril fragmentation. We note that disassociation, which occurs at a slower rate than aggregation mechanisms, is typically neglected in analysis of kinetic models [9,16]. Meanwhile, fragmentation was shown not to be a dominant mechanism for aggregate formation [17]. Given these assumptions, the possible types of reactions may be illustrated as where O a and P a give the average oligomer and fibril size, respectively. The first equation describes the primary nucleation process where we have assumed that the number of monomers making up the smallest stable aggregate, typically denoted by n c , is 2. The last two equations describe multi-step processes where we have assumed that there is a rate determining step that is slower than the others. This is the same formulation that is used in [73] to describe these reactions.
Since monomers are produced in the brain, we assume that its production process is represented by a saturating function f ( , which is the logistic differential equation first introduced by Verhulst [75], where δ is the growth rate of monomers and γ is the carrying capacity. This choice of function was motivated by the assumption that monomer concentration cannot grow unbounded. However, we note that the following results would also hold for a general nonlinearity f (M 1 ) satisfying the conditions f (0) = δ > 0 and f (M 1 ) < 0. We also assume that there is degradation/clearance of each stage, denoted by μ i . While we do not distinguish between degradation and clearance here, in general degradation occurs at a faster rate than clearance [29]. All together, these assumptions result in the following model: The factors (O a − n) and (P a − O a ) in the monomer equation represent the average number of monomers that need to be added to an M n aggregate to form an oligomer and the average number of monomers that need to be added to an oligomer to form a fibril, respectively. These factors ensure that the law of conservation of mass of the system is satisfied. The law states that the mass of the chemical components before the reaction is equal to the mass of the components after the reaction, that is, no monomers are lost or created through the aggregation process, as demonstrated by Lemma 2.1.

Lemma 2.1:
The following equality holds: Note that the above lemma shows that if the source and sink terms are zero, that is δ = 0, μ i = 0, i = 1, . . . , n, μ O = 0 and μ P = 0, then Equation (2) reduces to: which implies that n i=1 iM i + O a O + P a P = C, i.e. the total number of monomers in all the aggregates is constant in time.

Existence and stability of equilibria
We now study the existence and stability of the equilibria of model (1). This model contains an extinction equilibrium and a unique positive equilibrium, as shown in Theorem 3.1. Proof: The fact that E * = (0, 0, . . . , 0) is always an equilibrium of the system is clear. To prove (ii), we note first that any equilibrium must satisfy the following: From (4) we find that the coordinates M 2 , . . . , M n , O and P for any equilibrium point must satisfy: From these equations it is easy to verify that dM i dM 1 > 0, for i = 2, . . . , n, dO dM 1 > 0 and dP dM 1 > 0. Now, to show M 1 has a positive fixed point coordinate, we divide the first equation of (4) by M 1 to obtain Next, we study the local stability of the two equilibria of the model, E * and M * . The Jacobian matrix of (1) evaluated at E * is given by Clearly, the eigenvalues of J(E * ) are given by Next we compute the Jacobian matrix at the interior equilibrium M * , J(M * ), to obtain ⎛ To determine the stability of M * , we use the following result due to Gershgorin [62].
To this end, we assume that and n−1 Using (5)-(7) is easy to verify that From this and assumption (9) if follows that The following theorem summarizes the above stability analysis.
Proof: Applying Theorem 3.2 we will show that all the Gerschgorin disks lie on the lefthand side of the complex plane, i.e., all eigenvalues have negative real part. To this end, consider the first Gershgorin disk S 1 centred at a 11 = − δ γ M * 1 − 2K 1 M * 1 . Any eigenvalue that resides in S 1 must satisfy From this we have the following two-sided inequality for the real part of λ, Simplifying we get: by assumption (10). Thus such an eigenvalue must have a strictly negative real part. Next, note that for i = 2, . . . , n Thus, Applying assumption (9) and inequality (13), we obtain From (9) and (10) it follows that Thus any eigenvalue in the disk S i must have a strictly negative real part. Similar arguments can be applied to S n+1 and to S n+2 centred at a n+1,n+1 = −K O M * 1 − μ O and at a n+2,n+2 = −μ P , respectively, to show that any eigenvalue in these two disks must have a negative real part. From the above it follows that all eigenvalues of the matrix which are in S = ∪ n+2 i=1 S i have negative real parts. Hence, M * is locally asymptotically stable.

Remark 3.4:
One can also observe that J(M * ) is a block matrix with two blocks: the first (n + 1) × (n + 1) elements and a 1 × 1 block formed by the element J n+2,n+2 . Thus one of the eigenvalues is given by λ n+1 = −μ P < 0 and hence it is sufficient to apply the above Gershgorin theorem to the block matrix obtained from the first (n + 1) × (n + 1) elements of J(M * ) and show that the (remaining) n + 1 eigenvalues have negative real parts to obtain the conclusion of Theorem 3.3. The argument for showing that these eigenvalues have negative real parts is identical to that in the proof of Theorem 3.3.

Remark 3.5:
The relation between the chemical reaction and the activation energy is governed by the Arrhenius equation [45]. The activation energy is the minimum energy needed for the reaction to occur, expressed in joules per mole. The Arrhenius equation states that the rate constant K of a chemical reaction depends primarily on the activation energy E a and temperature T: where the pre-exponential term A is related to the geometric and other secondary factors that may affect the frequency of collisions, and k B is the Boltzmann constant [45]. For near isothermal reactions, the primary factor is given by the activation energy. The accumulated electron cloud gets stronger as more monomers are added up. Therefore, the energy level of M i+1 is lower than that of M i , for i = 1, . . . , n, that is, M i+1 is more stable than M i . For this reason, the reverse reaction is ignored. We may conclude that E a 1 > E a 2 > · · · > E a n . Given that the activation energy plays the primary role in determining the reaction rate constant, we obtain from the Arrhenius equation the relation:

Numerical investigation
In this section, we use numerical simulations to further explore the dynamics of model (1) for parameter values lying outside of the region of stability predicted by Theorem 3.3. In addition, we also consider the dynamics of two discrete-time approximations of this model. Deriving discrete-time approximations to continuous-time aggregation models that preserve non-negativity, conservation of mass and the local dynamics of the continuous-time model is a challenging problem. In the appendix, we derive the two discrete-time models considered here and show that they preserve non-negativity and the local stability conditions for model (1). However, these models do not preserve the conservation of mass property of the continuous-time model and, as a result, we observe that they may produce cyclic dynamics when the continuous-time model produces stable equilibria. Importantly, we note that the model obtained using non-standard finite difference discretization (NFSD) is equivalent to the model developed in [20] with the additional modelling assumptions introduced in this paper. Therefore, the numerical simulations demonstrate that the continuous-time model is better able to describe the dynamics of Aβ aggregation than the previously developed model in [20]. Here we use simulations to address two questions. First, since the stability conditions obtained in Theorem 3.3(ii) are not sharp, what type of dynamics occur when these conditions are not met? Do we still have stable equilibria or are more complicated dynamics possible? Second, how do the solutions to the two discrete-time approximations compare to the solution to the continuous-time model? We showed in Theorems A.3 and A.4 that the discrete-time models preserve the local dynamics. This indicates that these solutions should agree well when conditions (9)-(10) are met. However, these solutions are not guaranteed to agree for parameter values falling outside of the region defined by these conditions.
Biologically, the parameters involved in modelling the aggregation process of Aβ may have very different order of magnitude ranges. For instance, the nucleation rate may be as slow as 10 −6 M −1 s −1 while the elongation rate may be as fast as 10 4 M −1 s −1 [17]. Given this variability, where available we obtained parameter estimates, given in Table 1, from which we based our numerical simulations. However, we caution that the values used in the numerical simulations are solely for illustrative purposes and may not accurately reflect the rate constants for Aβ associated with AD as some estimates, such as those in [17,42], apply to forms of Aβ with enhanced aggregation.
The estimates for the elongation rate given in Table 1 are based on the addition of a single monomer. Since the elongation rate K O in this model involves the addition of P a − O a monomers, we modified this estimate based on the relative sizes of the fibril and oligomer stages. In [29], the authors estimated the degradation rate of oligomers as 1/10μ 1 . Therefore, to estimate the degradation rate of fibrils, we assume, based on relative sizes, that μ P = 1/10μ O = 1/100μ 1 . Since we may expect the degradation rates of aggregates to decrease with aggregate size, while the nucleation rates increase with aggregate size (see Remark 3.5), to pick values of μ i and K i , for 2 ≤ i ≤ n, we choose relationships that satisfy these. Specifically, we use the relations Notice that with this choice of K i , condition (9) is always satisfied. We were unable to obtain estimates for the production rate of monomers δ or the carrying capacity of monomers γ . Therefore, in our simulations, we choose an arbitrary value of δ that makes the extinction equilibrium unstable, δ − μ 1 > 0, and vary the value of γ to change the value of which gives a locally asymptotically stable positive equilibrium when this quantity is negative (see condition (10)).
In the simulations presented in Figures 2-5, we use the parameter values δ = 50, K 1 = 10 −4 , K O = 0.1, n = 6, O a = 10, P a = 100, together with the relationships given in (16) and the initial condition M 1 = 10, M 2 = · · · M 6 = O = P = 0. In Figure 2(a), we show the time series dynamics for models (1), (A1) and (A2) when γ = 75 resulting in D = 8.3639 > 0. This first figure demonstrates that the region of stability for the positive equilibrium extends beyond the stability region established in Theorem 3.3, as determined by condition (10). The transience for all three models follows a sigmoidal curve and, as is to be expected, the larger aggregates take longer to reach equilibrium, with the fibril stage taking much longer than the others. In Figure 2(b), we show the time-series dynamics when γ is increased to γ = 5000 resulting in D = 9.0206. In this figure, we observe that all models still converge to a stable positive equilibrium. However, both discrete-time models exhibit transient oscillations. Increasing γ even further to γ = 8000, resulting in D = 9.0244, we observe in Figure 3 that the discrete-time approximations no longer agree with the continuous-time model. Instead, model (1) equilibrates while the two discrete-time models converge to stable cycles that oscillate around this equilibrium value. Though we were unable to find stable oscillations in the solution to the continuoustime model for the parameter ranges given in Table 1, it is possible for this model to exhibit stable cycles for parameters outside of these ranges. This is demonstrated in Figure  3(b) where γ = 30, 000 and all other parameters are the same as Figure 3(a) except we have taken K O = 0.01 and μ 1 = · · · = μ n = μ O = μ P = 1. In Figure 4, we provide the transient solution to M 1 for the continuous-time model for each of the four scenarios considered in Figures 2 and 3. We note that the continuous-time model exhibits some transient oscillations for the parameter values in Figures 2(b) and 3(a), but these oscillations are less significant when compared to the oscillations observed in the discrete-time models. In Figure 4(c), we provide a clearer view of the cycles shown in Figure 3(b) for the continuous-time model.
As discussed in the appendix, the two discrete-time approximations (A1) and (A2) do not preserve the mass conservation property of the continuous model (see Lemma 2.1). Specifically, for both discrete-time models, mass is lost in the aggregation process. We suspect that this loss of mass is the cause of the instability observed in the numerical simulations. In Figure 5, we demonstrate how much mass is lost for the scenarios presented in Figures 2 and 3. In Figure 5(a), we plot the absolute error between the total mass concentration over time, as given by the solution to (2), and the mass concentration over time as predicted by models (A1) and (A2). As we can see from this figure, model (A2) does a slightly better job than model (A1) at preserving mass. In Figure 5(b), we give the relative error for both models. Here we observe that, as we move further away from stability condition (10) by increasing γ , the relative error in mass concentration increases.

Sensitivity analysis
Sensitivity analysis plays an important role in parameter estimation and optimization [1][2][3] and uncertainty [4,21]. Sensitivity analysis involves studying the effects of varying model parameters on the model output. For our model, this involves understanding the (n + 2) × (2n + 5) sensitivity coefficient matrix given by where for convenience we denote  Because our model parameters are measured in different units and differ by orders of magnitude, this presents a challenge in interpreting sensitivity results and makes it difficult to compare sensitivities between different parameters. Thus we focus on elasticity analysis as it is a more appropriate measure in such a case [12]. Elasticity analysis looks at the proportional response to a proportional (rather than additive) change in a model parameter. The elasticity coefficient matrix is given by In particular, elasticities are proportional sensitivities that are scaled so that they are dimensionless. Thus one can directly compare elasticities among all model parameters.
Here, we use a simple finite difference to approximate the derivative and numerically solve for the elasticity coefficient matrix. Specifically, we use the following scheme to numerically compute these elasticity coefficients for i = 1, . . . , n + 2, j = 1, . . . , 2n + 5 and small p j > 0. In Figure 6, we present the elasticity of the concentration of aggregates over time with respect to three parameters: K 1 , K 2 and γ . This graph was generated using the same parameter values as Figure 2(a) and gives the absolute value of the elasticities in order to more easily compare their relative magnitudes. From Figure 6, we can observe both the transient effects, which contribute to the time it takes different stages to equilibrate, as well as the asymptotic effects, which describe how model parameters impact final equilibrium concentrations. Similar qualitative behaviour is observed for other parameters but the magnitude (elasticity level) changes from one parameter to another. Thus in Table 2 we present the transpose of the full elasticity coefficient matrix at the final time of the simulation T = 800, i.e. the transpose of C E (T; p). Since n = 6 this matrix is of dimension 17 × 8. This table clearly demonstrates that aggregate equilibrium concentrations are most sensitive to γ , followed by K 1 . In general, we also observe that, when comparing the effects of  different elongation or degradation rates, the model outputs are more sensitive to the rate with a lower index. Exceptions occur for the stage that rate directly impacts, for instance changes in K 4 have a greater impact on M 4 than changes in K 2 or K 3 . For comparison purposes, we also computed the elasticities when γ is increased to γ = 5000 as in Figure 2(b) (not shown). These elasticities show the same general pattern as is given in Table 2 except that the impact of δ is increased, with the elasticities with respect to δ and γ having the same order of magnitude.

Conclusion
In this paper, we developed a continuous-time model for the aggregation of Aβ. This model uses size (in term of the number of monomers) to model the slow nucleation phase and stage to model the faster elongation phase. Motivated by previous studies [73], the model also explicitly distinguishes between oligomer and fibril stages, and their elongation rates.
The main limitation of this new model is that we assume oligomers and fibrils can be well characterized by average sizes for each of these stages. While describing a group according average characteristics is a common approach in stage-structured population modelling, such an approach may not be appropriate when there is a large amount of variability within a group. For instance, model (1) may be more appropriate for describing the early stages of AD progression, while, for later stages of AD, an extension of this model in which we include multiple oligomer and fibril stages may work better. We showed that the continuous-time model (1) has an extinction equilibrium and a unique positive equilibrium, and applied Gershgorin's Theorem to establish local stability. Numerical simulations show that, for parameter values falling within biologically reasonable ranges, this model exhibits the sigmoidal growth, Figure 2(a), that is predicted by theory [37]. These simulations also show that the stability region for the positive equilibrium extends beyond the stability region, determined by condition (10), established in Theorem 3.3. We conjecture that the positive equilibrium may, in fact, be globally asymptotically stable and we aim to establish conditions on the model parameters for global stability in future work. However, we point out that the stability of the positive equilibrium is not always guaranteed, as is demonstrated by the numerical results in Figure  3(b) and Figure 4(c) which show that the model may exhibit more complicated dynamics including periodicity of solutions. It is worth noting, however, that these oscillations were obtained using parameter values outside of the ranges given in Table 1. Meanwhile, numerical simulations suggest that the continuous-time model does not appear to exhibit oscillations when parameter values are chosen from the biological ranges provided in Table 1.
Other approaches, including discrete-time models [20], have been used to model Aβ aggregation. In order to compare how model predictions may vary based on the approach used, here we also developed two discrete-time approximations to the model (1) using non-standard schemes. We showed that these models preserve the local dynamics of the continuous-time model. Numerical simulations verify that, even outside of the stability region determined by (10), these difference schemes may provide good approximations to model (1). However, as we move away from this condition, the discrete-time models start to exhibit oscillatory transience and eventually stable cycles. This loss of stability appears to coincide with a loss of conservation of mass, as shown in Figure 5. This suggest that one should use caution when using discrete-time approximations to model aggregation processes.
The onset and progression of Alzheimer's disease can vary greatly between individuals due to factors that are not well understood [6]. In Section 4.1, we performed sensitivity analysis of the model outputs to understand how variability in the parameters describing the aggregation process of Aβ may influence the time series dynamics. When comparing the various elongation rates, we observe that, even after accounting for differences in order of magnitude of parameter values by using elasticity measurements, the impact of these rates on the final equilibrium values may vary by orders of magnitude with the elongation rates of smaller aggregate sizes having a larger impact on equilibrium values. In particular, this suggests that it may be important to have distinct elongation rates for the smaller versus larger aggregates. This was observed experimentally for tau aggregation [73], where the traditional kinetic equations were found to underestimated average fibril length and be inadequate at capturing oligomer concentrations over time. Meanwhile, a model that accounts for the conversion of oligomers to fibrils via a reaction of first order in monomer concentration, as we have in model (1), was shown to provide a better fit despite categorizing aggregates into only one of three stages: monomer, oligomer or fibril. We observe a similar pattern was comparing the degradation rates of the different aggregate sizes. All together, we observe that the parameters related to the production δ and prevalence γ of (abnormal) monomers as well as the nucleation rate K 1 , which leads to the creation of new aggregates, have the greatest effect on both the transient concentrations of aggregates and the final equilibrium concentrations. The sensitivity analysis presented here, together with knowledge on which model parameters are more likely to vary between individuals, may further help identify causes of individual variation in AD progression.
Finally, we note that the model presented in this paper is a first step toward providing a more detailed description of Aβ aggregation that will serve as a basis for more realistic models that account for additional aggregation mechanisms such as dissociation, fragmentation and secondary nucleation. Unlike the traditional kinetic equations that are often used to model Aβ, here we do not assume that the elongation rates for aggregates of different sizes are the same [42]. Without this assumption, the model in this paper has the disadvantage of being less readily analysable. However, the distinction between these rates may prove important for reproducing the results of kinetic studies [73]. In addition, the model is formulated to be flexible enough to allow very detailed information on smaller aggregates which can be experimentally measured. The results of the sensitivity analysis presented in this paper demonstrate the need for having good estimates for some of the model parameters that are highly elastic including parameters that affect the monomer production and monomer nucleation phase. In future work, we plan to estimate model parameters using dynamic light scattering (DLS) data which allows the approximation of aggregate size and molecular weight distribution as a function of time.

Disclosure statement
No potential conflict of interest was reported by the authors.

Appendix. Discrete-Time approximations to the continuous-time model
Various types of mathematical models have been used to describe the aggregation of Aβ or of other proteins. Here we contrast two modelling strategies, specifically continuous versus discrete time, by constructing two discrete-time approximations to the continuous-time model (1) using non-standard discretization methods.
The first method is based on a nonstandard finite difference discretization (NSFD) (see, e.g. [55]) which preserves the non-negativity of all variables, while the second approximation follows from an exact discretization scheme (NEDS) [44]. Though this latter method is not guaranteed to preserve non-negativity, solutions will remain non-negative provided that the removal terms μ i are not too large, as is the case here. The main advantage of the discrete-time approximations presented below in comparison with those available in the literature (including Runge-Kutta methods with and without adaptive time step size) is that these approximations are guaranteed to preserve the local dynamics of the continuous system even when the time-step size is large. Specifically, both discrete-time models have the same equilibria and local stability results as the continuous-time model (1). This is important from a practical point of view as often the step-size represents a generation time or the period of some empirical measurements and hence is fixed.
In our first approach, we approximate the derivatives by with h = 1. We then apply a non-standard discretization method on the continuous model (1) as follows. We approximate the quadratic term − δ γ M(t) 2 in the first equation of (1) and all the removal Proof: (i) The Jacobian matrix of system (A1) evaluated at E * is given by Thus E * is asymptotically stable if 1+δ 1+μ 1 < 1, which is equivalent to δ − μ 1 < 0, and unstable if 1+δ 1+μ 1 > 1, which is equivalent to δ − μ 1 > 0. Meanwhile, the Jacobian matrix of system (A2) evaluated at E * is given by Thus E * is asymptotically stable if δ − μ 1 < 0 and unstable if δ − μ 1 > 0.
(ii) Next we compute the Jacobian matrix at the interior fixed point M * . For models (A1) and (A2), this Jacobian has the form J(M * ) = u ψ ρ T .