Stochastic master equation for early protein aggregation in the transthyretin amyloid disease

It is significant to understand the earliest molecular events occurring in the nucleation of the amyloid aggregation cascade for the prevention of amyloid related diseases such as transthyretin amyloid disease. We develop chemical master equation for the aggregation of monomers into oligomers using reaction rate law in chemical kinetics. For this stochastic model, lognormal moment closure method is applied to track the evolution of relevant statistical moments and its high accuracy is confirmed by the results obtained from Gillespie’s stochastic simulation algorithm. Our results show that the formation of oligomers is highly dependent on the number of monomers. Furthermore, the misfolding rate also has an important impact on the process of oligomers formation. The quantitative investigation should be helpful for shedding more light on the mechanism of amyloid fibril nucleation.

Scientific RepoRtS | (2020) 10:12437 | https://doi.org/10.1038/s41598-020-69319-x www.nature.com/scientificreports/ condition, to describe the aggregation of monomers into oligomers. And then, we use lognormal closure moment method rather than direct simulation to acquire the time dependent evolution of low-order statistical moments, and we emphasize that the semi-analytic moment method can greatly reduce the computational cost [33][34][35][36] . The paper is organized as following. In "Chemical master equation model" we describe the stochastic method of modeling the aggregation of monomers into oligomers and present moment closure method for computing the time evolution of stochastic models. In "Results", we apply the moment closure method to the stochastic model and examine the numerical results by comparing with the simulated results. Moreover, we investigate the impact of the number of monomers and the misfolding rate on the aggregation dynamics. The conclusions are drawn in last section.

chemical master equation model
Build model of oligomers formation. Amyloid formation is considered to be a complex protein aggregation process, which consists of a range of molecular processes. In the process of protein filaments formation, primary nucleation is usually followed by growth [10][11][12][13][14][15][16][17][18][19] and self-replication through secondary pathways in some cases 11,17,18 . The nucleation phase corresponds to the period where monomers undergo conformational changes and self-associate to form the oligomeric nuclei. This phase is determined by the critical concentration of nuclei and generally considered to be thermodynamically unfavorable 21 . The growth/elongation phase represents the period in which the oligomeric nuclei acting as seeds rapidly grow and form mature fibrils. Considering the nucleation phase is an essential part of the overall aggregation process, as a large variety of oligomeric species are gradually formed 21 , so we ignore elongation but focus on the process of nucleation below.
First, we briefly present the stochastic formulation inspired by Ref. 35 . If a model has L biochemical reactions and N molecular species, it can be written as where s li and r li are the coefficients of change of the molecular species X i involved in the lth reaction. We define the stoichiometric matrix as S = (S lj ) where S lj = r lj − s lj denoting the change of the number of molecular species X j by the lth reaction.
Under well-mixed and dilute conditions, if P(x, t) denotes the probability for the state vector x = (x 1 , x 2 , . . . , x N ) at time t , then the probability evolution of the system (1) can be described by chemical master equation (1) s l1 X 1 + · · · + s lN X N → r l1 X 1 + · · · + r lN X N , 1 ≤ l ≤ L   www.nature.com/scientificreports/ where S l = (S l1 , . . . , S lN ) is the row vector of the stoichiometric matrix and a l (x) is the propensity function for the lth reaction determined by mass action kinetics. For example, in chemical kinetics 38 , for the reaction s 1 X 1 + s 2 X 2 → r 1 X 3 + r 2 X 4 , the reaction rate law expression is where K is the reaction constant. Then, we establish the stochastic mathematical modeling of making oligomers. Denote the number of In modeling, the source term of monomers written as x 1 (t) → x 1 (t) + 1 is also considered and the production rate is assumed to be k p . For the system described in Fig. 1, the stoichiometric matrix is ,then the chemical master equation of the model is The above chemical master equation model describes the initial formation of oligomers in a stochastic sense. We remark that it is essential stochastic generalization of a time discrete deterministic model 32 and a continuous time deterministic counterpart 39 for the aggregation process. Actually, elegant analytical solutions have been developed to the time continuous model without an upper limit to oligomer size 39 , and the time discrete model can be regarded as its simplification with an upper limit to the oligomer size.
Moment closure method. Although the chemical master Eq. (4) is a powerful tool for describing the stochastic dynamics of the system (1), it is difficult to find an exact analytic solution. Instead, some simulation techniques such as Gillespie's stochastic simulation algorithm (SSA) have been presented 40 . But the SSA is relatively expensive in computational cost especially when the system size is large. To overcome this shortcoming, moment closure techniques for approximating the low-order moments have become more and more popular 33,34,[41][42][43] . This is acceptable because the first two order statistical moments (mean and variance or covariance) are sufficient for a decent description of the ensemble dynamics such as averaging behavior and evolution of the noise in the system 41 .
The moment equations corresponding to Eq. (4) is not self-closed since the evolution of the Mth order moments depend on the (M + 1)th order moments. That is, the evolution of the resultant first two order moments involves the third order moments. The so called moment closure is to approximate the involving higher order moments as nonlinear functions of lower-order moments. The frequently adopted closure schemes include Gaussian moment closure, lognormal moment closure, gamma moment closure and binominal moment closure. Here, we apply a kind of lognormal moment closure to Eq. (4) for two reasons. The first reason is that the important features of asymmetry and nonnegativity of the molecular reaction system make that population of species in bio-statistics tends to be lognormal distributed 34 . The second reason is that the lognormal closure scheme, also named the derivative matching closure as presented in Refs. 42,43 , does not necessitate assumptions about the priori distribution. In this method, the nonlinear closure functions are given by matching time derivatives of the unclosed exact moment equations with those of the approximate closed moment equations at some initial point.
Given m = (m 1 , m 2 , . . . , m 7 ) ∈ N 7 ≥0 , we define moment generating function to be ×P(x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , t) The detailed derivation can be found in supplementary method. With the moment generating function (5) in mind, the general form of moment equations is given by rewritten the Eq. (6) and extraction of the coefficients of θ 1 , θ 2 , . . . , θ 7 where a l,i is the multinomial coefficient of a l (x).
Since we mainly focus on the important stochastic quantities namely mean and variance, the system of 2C 1 7 + C 2 7 = 35 equations for the first two order moments are derived shown in supplementary equations, from which we can see these equations depending on third order moments. For calculating, we truncate the infinite hierarchy (7) to a finite-dimensional system by the lognormal closure scheme. Following the lognormal closure scheme 43 , the third order moments can be expressed by the nonlinear function of first two order moments. Assume µ m (m = (m 1 , . . . , m 7 ) ) be one of the third order moments, then we can find a suitable closure function φ (m) (·) with the separable form where µ = [µ m 1 , . . . , µ m s ] is the vector of the first two order moments and γ 1 , . . . , γ s are chosen as the unique solution of the following set of linear equations With the closure functions in Eq. (8) available, all the involving third-order moments can be expressed by the first two order moments, thus a self-closed moment system consisting of thirty-five ordinary differential equations can be deduced from Eq. (7). The detailed deduction based on the lognormal closure scheme (8) can be found from the supplementary. We emphasize that all the deductions are implemented by hand. We then apply the fourth order Ronge-Kutta method to the closed moment system and the first-order and the second-order moments are solved.

Results
It is generally hard to determine the unknown rate constants of oligomer growth. For simplicity, we assume all the polymerization rates are identical, namely K 1 = K 2 = K 3 = K 4 = K 5 = K a and all the polymerization rates are the same, i.e. K 7 = K 8 = K 9 = K 10 = K 11 = K b . Although the identical assumption as an approximation to the nonlinear system is not that realistic, our numerical experience shows that the accuracy of the lognormal closure scheme does not rely on the relative variation of these rate parameters. As oligomeric intermediates in fibril formation are thermodynamically unstable 24,44,45 , we specify equilibrium constant k e = K a /K b around 0.02, which ensures oligomers are suitably unstable compared to monomers. Since monomers are usually not produced in vitro experiments, we set the rate of monomer production k p = 0 . We assume the refolding rate ( K 6 ) is larger than misfolding rate ( K 0 ) in all the simulation. The initial conditions for the seven species are taken as x 2 (0) = . . . = x 7 (0) = 1 and x 1 (0) = 2000 not including Fig. 3. Figure 2 shows the mean and the standard deviation we obtain for monomers, misfolded monomers, diamers, … and hexamers. It can be seen that the results derived from the lognormal closure are in good agreement with the results obtained from the SSA realization. This coincidence implies that the moment method is efficient and accurate in capturing the time evolution of mean and standard deviation. Thus, we only show the results obtained by the moment closure method in the other figures. As seen from the figure, all the reaction species can evolve into a steady state for suitable parameters. Figure 3 shows that the quantitative evolution of the numbers of diamers, …, hexamers under different initial numbers of the monomers. It is shown from the picture that the initial number of monomers has same effect on different tapes of oligomers. That is, as the initial number of the monomers increases, the number of oligomers increases and the difference among initial values also undergoes a transient increase. This result implies that increasing the initial number of the monomers will speed up the process of nucleation phase. This observation supports the statement in Ref. 21 that the aggregation process is highly dependent on the number of monomers.
(e θS l − 1) x e θx P(x, t)a l (x) Scientific RepoRtS | (2020) 10:12437 | https://doi.org/10.1038/s41598-020-69319-x www.nature.com/scientificreports/ Note that the number of the monomers of transthyretin is typically at steady state but may vary from person to person, thus the observation from Fig. 3 may be helpful in explaining why some people are more apt to suffer from transthyretin amyloid disease under same condition to some extent. Figure 4 shows the influence of misfolding rate on the number of oligomers. This qualitative observation in Fig. 4 demonstrates that the misfolding rate has a major impact on the aggregation process, that is, as the misfolding rate increases, all the five types of oligomers grow in number, and the time for these oligomers to reach a given level also becomes shortened. Since the toxic oligomers accumulation plays a central role in amyloidoses 46   conclusion It is well known that amyloidal aggregation is the hallmark of amyloidoses. In order to quantitatively explore the process of early amyloidal aggregation, we have developed a stochastic mathematical model about oligomers aggregation from monomers using rate law in chemical kinetics. We have adopted a typical moment method based on lognormal closure to capture the statistical moments, and massive calculations show very good agreement between the semi-analytic method and the Gillespie's SSA for the stochastic model. Our results show that the aggregation of monomers into oligomers is highly dependent on the number of monomers and the misfolding rate, and decreasing the number of monomers and the misfolding rate can inhibit the formation of the toxic oligomers. Our research may be helpful in explaining the individual variation in suffering from transthyretin amyloid disease and emphasizes the importance of controlling the misfolding of protein in preventing this type of disease. Additionally, we remark that the lognormal moment method is successful for stochastic master equation relevant to the transthyretin with an upper limit of oligomer size N = 6 , and we expect this method would also be effective with a different upper limit which might be relevant to other amyloidogenic systems. The method in the present work can provide essential quantitative insights into the mechanism of the early steps in the aggregation reactions. The mechanism behind the formation of amyloidoses is far from clear, and in the future we will continue to investigate the kinetic of amyloid nucleation and the mechanism of amyloid fibril growth.