Building Reaction Kinetic Models for Amiloid Fibril Growth

In this work we  discuss some methodological aspects of the creation and formulation of mathematical  models describing the growth of species from the point of view of reaction kinetics. Our discussion is based on familiar examples of growth models such as logistic growth and enzyme kinetics. We   propose several reaction network  models  for  the amiloid fibrillation processes in the citoplasm. The solutions of the models are sigmoidal functions graphically visualized using  the computer algebra system Mathematica .


I. INTRODUCTION
A field of considerable interest is the study of various biological growth processes and the application of mathematical models that facilitates the understanding of these processes.Growth processes usually evolve in time as sigmoidal functions.There exists a vast literature on sigmoidal functions.The field is characterized by a huge number of studies on real world growth phenomena and attempts to explain the intrinsic mechanisms of these phenomena using various mathematical methods [7], [8], [21], [33].
Sigmoidal growth functions are usually introduced in three main ways.Often growth functions are defined by an explicite arithmetic expression.Another way is to define them as a solution of a problem formulated in terms of a differential equation or a system of (integro-)differential equations.A third way is to formulate a chemical reaction network that induces (via the mass action law) a dynamical system, that in turn imply sigmoidal solution(s).This approach makes use of the reaction network theory-a well established field of applied mathematics (mathematical chemistry) that studies the behavior of real world chemical systems.In many situations the reaction network-approach can be applied to biological phenomena and has the advantage of suggesting possible (bio)-chemical mechanisms of the processes and phenomena under investigation.
In this work we are going to illustrate the above mentioned approaches for the formulation and study of growth functions on several familiar examples such as the Verhulst logistic function and Henri enzyme kinetic reaction network.We shall then formulate some models that can be possibly used to explain the growth of amiloid fibrils in the cell citoplasm.
Sigmoidal growth curves typically have three parts (phases): lag, log and stationary parts.It is a challenging question to characterize mathematically these phases.The lag phase is practically important in many medical and biotechnological applications as this phase is responsible for the acceleration or inhibition of the process and the possibility of controlling the lag phase depends on the understanding of the hidden mechanisms of the corresponding process.The reaction networkapproach to growth process modelling provides namely certain knowledge about the specific intrinsic mechanisms of the process.

II. BIOLOGICAL GROWTH/TRANSITION MODELS: THREE CASE STUDIES
In this section we present three case studies of familiar growth models in order to illustrate possible mechanisms of growth formulated in terms of reaction networks.Growth and transition are related processes, as a species (reactant, population) A grows for the expenses of another species B. In some situations this can be also expressed as "species B transits (transforms) into species A".The transition can be in both directions (reversible).The process can be catalyzed by a third species C, which in particular can coincide with some of the species A, B (autocatalysis).
Three familiar forms for presentation of the models will be illustrated by means of case studies.Explicitly formulated growth curves (models), also known as "empirical" models, are briefly denoted as E-type models.Models formulated in terms of systems of differential equations (dynamical systems) are denoted as D-type models.Finally, models formulated in terms of a (chemical) reaction network of certain reactants (species, populations) are classified as R-type models.Note that a model can have several types of formulation.An R-type model can be reformulated into a D-type model for the concentrations of the reactants by means of the mass action law [6], [32].For some growth models all three formulation types are available, as shown in case studies 1 and 2 below.

A. Case study 1: saturated growth
The saturated growth is not sigmoidal one, as no lag phase is present.Nevertheless this growth model is practically important and is a good illustration of the three model types.
A (chemical) reaction network comprises a set of reactants, a set of products (often intersecting the set of reactants), and a set of reactions [6], [30], [32].
Consider the simple reaction network consisting of the two species (reactants) S, P and a transition reaction of S transforming into P with rate k, symbolically: In a situation when it is meaningful to speak of concentrations of the reagents (reacrants), e. g. when reaction (1) takes place in a liquid medium, then the concentrations s, p of the species S, P , resp., obey the mass action law and we obtain the dynamical system ds/dt = −ks, dp/dt = ks. ( Let us formulate an initial value problem (IVP) assuming initial conditions s(0) = s 0 = 1, p(0) = p 0 = 0. Then the first equation of (2) has as solution the exponential decay function: To find an expression for the concentration p of the product species P , we note that the conservation law for system (2): s + p = 0, induces s(t) + p(t) = s 0 .Substituting s = s 0 − p = 1 − p in the second equation of (2) we have Equation ( 4) with p(0) = 0 has as solution In this case study all three formulation types for the growth function p are present: the R-type (1), the D-type (2) and the E-type formulae for the saturation growth of p (5) as well for the decay of s (3).
Consider now reaction (1) slightly modified by adding a reaction in the reverse direction.The Rtype formulation is then: which is a shortcut for the following reaction network: Applying the mass action law to (6) we obtain the following D-type model for the evolution of the concentrations x, y of the reactants X, Y , resp.: Again an IVP can be considered by specifying the initial conditions x(0) = x 0 , y(0) = y 0 .The solutions x, y are graphically visualized on Fig. 1 for x 0 = 1.0, y 0 = 0.0.Assuming k ≥ k −1 means that x decays and y is growing.The saturation growth process has no lag phase.
Remark.Finding an R-type formulation for a certain D-model is called realization of the Dmodel [6].It is instructive to look for a realization of the Malthusian D-type model x = kx.A possible realization is: X k −→ X + X.According to this R-type model species X reproduces without using any resources, which may be a rather rough approximation of a real process (in the long term).
A generalization of the decay-saturation mechanism for many (more than two) species is discussed in [24].There we study an extension of the reaction network (6).Let us rewrite the latter in the form: Reaction network ( 7) can be extended for n species X 1 , X 2 , ..., X n , so that each pair of species interacts, that is: For example, for n = 3 we have six reactions From the R-type formulations (8), ( 9) one can straightforward obtain the corresponding D-type models applying the mass action law, cf.[24], [40].

B. Case study 2: Verhulst logistic growth
The logistic function is a smooth sigmoidal function finding numerous applications in biochemical, population and cell growth phenomena [26], [41]- [43].Consider the following autocatalytic reaction network: A possible "biological" interpretation of reaction network (10) in the context of population dynamics can be the following: the nutrient substrate (or species) U is utilized (consumed) by species (population) X leading to a reproduction of species X, thereby k is a specific growth rate of the process.
Assuming that U and X are uniformly spaced in a certain volume (or area), then we may denote the biomass of population X by The solution x to IVP (12) passing through the .point (0, x(0) = x * = 1/2) is the (basic) logistic sigmoid function: We see that the logistic model admits all three formulation types: the E-type (13), the D-type (11) or (12) and the R-type (10).
In analogy to the generalisation of the decaysaturation R-type mechanism for many (more than two) species (8) we can generalize the logistic Rtype mechanism i) introducing a reversible reaction, and ii) introducing many (more than two) species.
For the case of a reversible reaction we have: For the case of a (irreversible, closed) food chain of three species we have Reaction networks ( 14), ( 15) can be easily generalized for n species X 1 , X 2 , ..., X n .This demonstrates again the notational power of the R-type model formulations.

C. Lag time
The so-called lag time is a mathematical characteristic of the concept of lag phase which is the low rate time period of a sigmoidal process.We are going to define and calculate the lag time of the logistic model.To this end consider the shifted logistic function on R: Function ( 16) has an inflection point (γ, 1/2) and its slope κ at t = γ is κ = k/4.Let (γ − δ, 0), δ = 2/k be the point where the tangent through the inflection point intersects the abscissa.Assume that γ ≥ δ.The value is called the lag time of the logistic model ( 16), cf.[39].Expression (17) has been obtained from the E-type logistic model (16), however ( 17) can be derived from the D-type model as well.Indeed, from (12) we obtain showing that, if the initial condition x(0) is less than 1/2, then (since x is positive) at some point γ > 0 with x(γ) = 1/2 there exists an inflection (satisfying 1 − 2x(γ) = 0).Thus, at γ we have x(γ) = 1/2, x (γ) = κ, x (γ) = 0. Substituting the time t = γ in (12) we obtain Let us examine the solution x 0 of the IVP problem (12) with initial condition x 0 (0) = 1/2.
Proof.Let x be an arbitrary solution to (12), then we have Due to x > 0, x can be zero only for values of t satisfying 1 − 2x = 0, that is x(t) = 1/2.This is true for all solutions of (12), in particular for solution x 0 , satisfying initial condition x 0 (0) = 1/2.Hence, due to x 0 (0) = 0, we have that x 0 has an inflection at 0. Similarly, any shifted solution x γ (t) = x 0 (t−γ) has an inflection point at t = γ, as Problem.Find the value of γ so that the shifted function x γ (t) = x 0 (t − γ) satisfies the initial condition x γ (0) = x 0 (−γ) = x * for a given x * , 0 ≤ x * ≤ 1/2.
To solve the above problem we make use of the E-type formulation of the logistic model.Let x 0 (−γ) = x * < 1/2, γ > 0, then x * is the initial value for the D-type model (12) having as solution the shifted function x γ , that is x γ (0) = x * .Using the E-type presentation (13), the latter gives (1 + e kγ ) −1 = x * , hence In order to have t lag ≥ 0, the restriction x * ≤ (e 2 + 1) −1 should be satisfied.We summarize the above calculations in the following Proposition II.2.If the initial value x * of the IVP ( 12) is such that x * ≤ (e 2 + 1) −1 , then the sigmoidal solution x = x(t) has a positive lag time.
Proposition II.2 shows that the D-type model (12) should have a sufficiently small initial condition in order to possess a positive lag time.Growth processes with positive lag time are typical for bio-chemical reactions, that is reactions involving functional proteins, such as enzymes, receptors, ligands, as well as population models.
Finally we shall point out one more characteristic property of the lag time of the logistic growth function.To this end let us first note that any sigmoidal function induces two simple (non-smooth) functions-a so-called cut (or ramp) function and a step function.More specifically consider the shifted logistic function (16) with inflection point (γ, 1/2) and slope κ = k/4 at t = γ.The tangent through the inflection point hits the abscissa t the point A = (γ − δ, 0), δ = 2/k and hits the horizontal line with ordinate 1 in the point B = (γ + δ, 1).The segment of the tangent between the two points A, B together with the parts of the abscissa before A and the part of the horizontal line with ordinate 1 after point B is the graph of the cut function c γ induced by the logistic function x γ .Similarly one defines the step function through the inflection point of the logistic function [2].We can now formulate the following Proposition II.3.[12] The uniform distance between any logistic function and its induced cut function is (1 − e −2 ) −1 .
Noticing that the uniform distance mentioned in the above proposition is the value of the logistic function at the point γ−δ, we see that this distance does not depend on the slope κ of the logistic function (and the induced cut function as well).For a situation when κ is small it seems not natural to consider the point γ − δ as a definition of the lag time.That is why in a series of papers we propose another definition of lag time, namely γ − δ wherein δ is the Hausdorff distance between the sigmoidal function and the induced step function [17]- [25].

D. Case study 3: growth models using Henri's reaction scheme
Biochemical processes involve functional proteins.The simplest reaction network involving a protein has been first formulated by Victor Henri [9], [14], [15], [16].The Henri reaction network involves two fractions of the enzyme (free and bound) denoted E and C, resp.: It is assumed that the rate parameters k 1 , k −1 , k 2 are positive constants such that k 1 > k −1 .Reaction scheme (18) describes the reaction mechanism between an enzyme E with a single active site and a substrate S, forming reversibly an enzyme-substrate complex C, which then yields irreversibly a product P .Reaction scheme (18) says that during the transition of the substrate S into product P the enzyme E bounds the substrate into a complex C having specific properties different than the properties of the free enzyme and thus being necessarily considered as a separate substance.
Proof.The solutions of ( 19)-( 20) are smooth functions.The first equation implies that s is monotonically decreasing function tending to zero with t −→ ∞ (note that k 1 > k −1 and e ≤ e 0 , c ≤ c 0 are bounded).Equation p = k 2 c ≥ 0 implies that p is a monotone increasing function.Function c increases from c(0) = 0 up to a maximum value c(t * ) achieved at some time t * (not necessarily unique) and then decreases tending to 0 with the exhaustion of s, see the second equation in (21), so c (t * ) = 0. Hence p (t * ) = k 2 c (t * ) = 0. meaning that p has an inflection at t * .From s +c +p = 0 we have s+c+p = s 0 showing that p(t) t−→∞ = s 0 .Therefore function p is sigmoidal.
The solutions to (19) with initial conditions s(0) = 0.5, e(0) = 1.0, c(0) = 0.0 p(0) = 0.2 are visualized on Fig. 4. In practice the rate constants k 1 , k −1 , k 2 are often not known and have to be determined for every specific enzyme-substrate pair.The contemporary approach to this task is to consider the rate constants as parameters in the dynamic system (19) and to compute them by fitting the solutions of the system to time course experimentally measured data [11].
The Verhulst and Henri reaction network models (10), (18) present two useful mechanisms for studying biochemical and biological growth.Various combinations of these two models have been proposed for the study of the EPS production [27], [28], [34].

FIBRILLATION
Recent intensive research into the physicochemical properties of amyloid and its formation into fibrils in the citoplasm points attention to growth models [3], [29], [38], [31].In [39] the authors consider the growth of aminoid fibrils and look for a mechanistic explanation of the process in terms of a biochemical reaction network.Fibril is an olygomer composed by monomers, thus Shoffner-Schnell model [39] involves two reactants: fibril F and monomer M , and additionally the intermediate reactant C. Reacrant C is the fibril "in action", that is the fibril that at the given time moment is in the process of storing the monomer molecule (adding it to self in a compact form).
The Shoffner-Schnell model describes the mechanism of the fibril growth in details and leads to interesting results.For educational purposes we present below several simple R-type models that may be helpful in illuminating certain particular steps of the fibril growth process and certain issues of interest (such as the lag phase).All presented models make use of the Verhulst and Henri reaction networks.

A. A basic Verhulst-Henri model
Let M be the total amount (concentration) of monomer, F be the fibril and C = M -F be the monomer-fibril complex at the time t of aggregation [3], [29], [38], [39].After aggregation the complex C turns into fibril F , that is, the added monomer molecule M converts into (part of) the fibril F .A simple reaction network model of these processes is Reaction scheme (22) is almost same as (18) with product P substituted by fibril F (obtained as result of the aggregation of the monomer M ).On the other side reaction network (22) can be seen as an extension of Verhulst reaction network (10) by involving an intermediate catalyst C.
A conservation law for system ( 23) is m + f + 2c = 0 implying the relation m + f + 2c = const = m 0 + f 0 = a.We have m = a − f − 2c in system (23) reducing it to From (25) we have f + c = k c c. Assuming that the monomer m is abundant in a time interval ∆ we can apply the QSSA principle and assume that c is approximately equal to 0 in ∆, and accordingly, in ∆ we have approximately f = k c c [1], [4], [5], [13], [32], [35], [36], [37].Hence, for some t * ∈ ∆ we have f (t * ) = k c c (t * ) = 0, hence function f has inflection and has a sigmoidal form.In addition, from f + c = k c c and c (t * ) = 0, we obtain the slope at the inflection point: f (t * ) = k c c(t * ) > 0. We thus proved the following Proposition III.1.The solution f of the IVP ( 23)-( 24) is a sigmoidal function.The slope at the inflection point is f The solutions m(t), f (t) of model ( 23) are visualized on Figure 5 for the following input data k + = 0.1, k − = 0.05; k c = 0.2; m 0 = 10; f 0 = 0.1; c 0 = 0 (in the time interval [0,200]).It should be noted that the fibril growth sigmoid function possesses a clearly expressed lag phase.
Remark.The inverse problem-so-called parameter identification problem-is to find out val-Fig. 5. Graphs of the solutions m, f of model (23) ues for the rate parameters and the initial conditions by fitting (some of) the solutions to available time course experimental measurements.Computational tools for the solution of this problem are proposed in [10], [11].

B. Three variants of the basic model
Our basic model describes the aggregation of the monomer while the monomer is compressed and becomes part of the fibril.
Below we propose three more reaction networks.All they involve an intermediate product P representing the aggregated monomer before turning into fibril.The three models present possible mechanisms for the transition of the aggregated monomer into fibril particles.
Reaction network 1.An intermediate product P is added as follows: Here the transition of the aggregated monomer P into fibril follows the decay-saturation mechanism (1).Applying the mass action law we obtain the ODE system The transition of the aggregated monomer P into fibril according to reaction network (28 follows the Verhulst autocatalitic mechanism (10).Applying the mass action law we obtain the ODE system:

IV. CONCLUSION
The proposed reaction networks (R-models) ( 22), ( 27), ( 29), (31) together with the implied reaction dynamical systems of equations (D-models) describe possible biochemical mechanisms of the fibril elongation in the citoplasm.The proposed models can be used to fit time course data for real measurements of fibril growth using fitting simulation procedures for the identification of the parameters.The variants producing a good fit will be considered as candidate for a probable mechanism of the amiloid fibrillation process on the base of its generating reaction network.
The application of mathematical modelling in the field of amiloid fibrillation necessarily focuses attention to the lag phases of the growth curves that appear as model solutions.Therefore it is an open problem to study the above discussed models with respect to this practically important issue.In particular the various reaction networks can be compared with respect to the lag phases of their solutions under various sets of parameters and initial conditions of the related differential IVPs.

( 23 )
If the three rate constants k + , k − , k c are known, system (23) can be treated as a Cauchy ODE IVP with initial conditions

Fig. 6 .
Fig.6.Graph of the solutions m, f, c, p of model(27) for system(29) is m + f + p + 2c = 0. Reaction network 3.Here an intermediate product P is added as follows: + F.(30) In the above reaction network (30 the transition of the aggregated monomer P into fibril follows the Henri enzyme kinetic mechanism(18).Applying the mass action law we obtain:dm/dt = −k + f m + k − c, df /dt = −k + f m + (k − + k c )c + k c p − k f p + 2k p c p , dc/dt = k + f m − (k − + k c )c, dp/dt = k c c − k f p + k c p , dc p /dt = k f p − (k + k p )c p .(31)A conservation law for system(31) is m + f + p + 2c + 2c p = 0.