www.biomathforum.org/biomath/index.php/biomath REVIEW ARTICLE Building reaction kinetic models for amiloid fibril growth Svetoslav Markov Institute of Mathematics and Informatics Bulgarian Academy of Sciences Sofia, Bulgaria smarkov@bio.bas.bg Dedicated to the 180th anniversary of the birth of Nestor Markov, a pioneer of lexicography and phys-math education in Bulgaria Received: 30 June 2016, accepted: 31 July 2016, published: 8 August 2016 Abstract—In this work we discuss some method- ological 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 en- zyme 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. Keywords-reaction kinetics; reaction network; growth models; sigmoidal functions; dynamical sys- tems; ODE’s I. INTRODUCTION A field of considerable interest is the study of various biological growth processes and the application of mathematical models that facili- tates the understanding of these processes. Growth processes usually evolve in time as sigmoidal functions. There exists a vast literature on sig- moidal 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 intro- duced 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 equa- tion 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 so- lution(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 un- der investigation. In this work we are going to illustrate the above mentioned approaches for the formulation and study of growth functions on several familiar Citation: Svetoslav Markov, Building reaction kinetic models for amiloid fibril growth, Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 1 of 11 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth 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 mathemat- ically 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 network- approach to growth process modelling provides namely certain knowledge about the specific in- trinsic 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 pos- sible mechanisms of growth formulated in terms of reaction networks. Growth and transition are related processes, as a species (reactant, popula- tion) 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 (re- versible). 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 sys- tems) are denoted as D-type models. Finally, mod- els 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 illus- tration 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: S k−→ P. (1) 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. (2) Let us formulate an initial value problem (IVP) assuming initial conditions s(0) = s0 = 1, p(0) = p0 = 0. Then the first equation of (2) has as solution the exponential decay function: s(t) = exp(−kt). (3) To find an expression for the concentration p of the product species P , we note that the conser- vation law for system (2): s′ + p′ = 0, induces s(t) + p(t) = s0. Substituting s = s0 −p = 1 −p in the second equation of (2) we have dp/dt = k(1 −p). (4) Equation (4) with p(0) = 0 has as solution p(t) = 1 − exp(−kt). (5) Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 2 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth 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 R- type formulation is then: X k−→←− k−1 Y, which is a shortcut for the following reaction network: X k−→ Y, Y k−1−→ X. (6) 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.: dx/dt = −kx + k−1y, dy/dt = kx−k−1y. Again an IVP can be considered by specifying Fig. 1. Solutions to the saturation-decay model (6) the initial conditions x(0) = x0, y(0) = y0. The solutions x,y are graphically visualized on Fig. 1 for x0 = 1.0, y0 = 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 D- model [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 mech- anism for many (more than two) species is dis- cussed in [24]. There we study an extension of the reaction network (6). Let us rewrite the latter in the form: X1 k12−→ X2, X2 k21−→ X1. (7) Reaction network (7) can be extended for n species X1,X2, ...,Xn, so that each pair of species interacts, that is: Xi kij−→ Xj, i = 1, ...,n, j = 1, ...,n, i 6= j. (8) For example, for n = 3 we have six reactions X1 k12−→ X2, X2 k21−→ X1, X2 k23−→ X3, X3 k32−→ X2, X3 k31−→ X1, X1 k13−→ X3. (9) 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 bio- chemical, population and cell growth phenomena [26], [41]–[43]. Consider the following autocat- alytic reaction network: U + X k−→ X + X. (10) A possible “biological” interpretation of reaction network (10) in the context of population dynam- ics can be the following: the nutrient substrate (or species) U is utilized (consumed) by species (pop- ulation) 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 x and the mass (or concentration) of the (nutrient) substrate U by Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 3 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth u and apply the mass action law, to obtain the dynamical system du/dt = −kxu, dx/dt = kxu. (11) Adding initial conditions u(0) = u∗ > 0, x(0) = x∗ > 0 to (11), one obtains a corresponding IVP. The solution of (11) for u∗ = 1.0, x∗ = 0.1 is visualized on Fig. 2. The conservation relation du/dt + dx/dt = 0 implies u + x = x∗ + u∗ = const = a. Let us assume a = 1 meaning that the initial conditions are normalized by 1/a, setting u∗ := u∗/a < 1, x∗ := x∗/a < 1. We then substitute u = 1 − x in the differential equation for x to obtain the Verhulst differential equation: dx dt = kx (1 −x) , x(0) = x∗ < 1. (12) The solution x to IVP (12) passing through the Fig. 2. Solutions to (11) with u(0) = u∗ = 1.0, x(0) = x∗ = 0.1.. point (0,x(0) = x∗ = 1/2) is the (basic) logistic sigmoid function: x0(t) = 1 1 + e−kt . (13) 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 decay- saturation R-type mechanism for many (more than two) species (8) we can generalize the logistic R- type mechanism i) introducing a reversible reac- tion, and ii) introducing many (more than two) species. For the case of a reversible reaction we have: X1 + X2 k12−→ 2X2, 2X2 k21−→ X1. (14) For the case of a (irreversible, closed) food chain of three species we have X1 + X2 k2−→ 2X2, X2 + X3 k3−→ 2X3, X3 + X1 k1−→ 2X1. (15) Reaction networks (14), (15) can be easily gen- eralized for n species X1,X2, ...,Xn. This demon- strates again the notational power of the R-type model formulations. The basic logistic function (13) is a sigmoidal function with asymptotes limt→−∞x0(t) = 0, limt→∞x0(t) = 1. Due to x′′0(0) = 0, x0 has an inflection at 0, inflection point is (0, 1/2), cf. Fig. 3. Assume that the tangent to the graph of function (13) through the inflection point (0, 1/2) intersects the abscissa at the point (−δ, 0),δ > 0. The slope κ of the tangent through the inflection point (0, 1/2) is equal to κ = x′0(0) = k/4. Thus, we have (1/2)/δ = k/4, hence δ = 2/k. The value of δ may be called log time (or high rate time period). C. Lag time The so-called lag time is a mathematical char- acteristic 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: xγ(t) = xγ(k; t) = 1 1 + e−k(t−γ) . (16) Function (16) has an inflection point (γ, 1/2) and its slope κ at t = γ is κ = k/4. Let (γ − Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 4 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth Fig. 3. A logistic function (13) with γ = 0 and reaction rate k = 40. δ, 0), δ = 2/k be the point where the tangent through the inflection point intersects the abscissa. Assume that γ ≥ δ. The value tlag = γ −δ = γ − 2/k (17) 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 x′′ = (kx(1 −x))′ = kx′(1 − 2x), 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 in- flection (satisfying 1 − 2x(γ) = 0). Thus, at γ we have x(γ) = 1/2,x′(γ) = κ,x′′(γ) = 0. Substituting the time t = γ in (12) we obtain x′(γ) = kx(γ)(1−x(γ)) implying κ = k(1/2)(1− 1/2) = k(1/4), hence κ = k/4. Thus, again (1/2)/δ = k/4, hence δ = 2/k and (17) holds true. Let us examine the solution x0 of the IVP problem (12) with initial condition x0(0) = 1/2. Proposition II.1. The solution x0 of (12) with ini- tial condition x0(0) = 1/2 has an inflection point (0, 1/2). Any shifted solution xγ(t) = x0(t − γ) has an inflection point at (γ, 1/2). Proof. Let x be an arbitrary solution to (12), then we have (1/k)x′ = x(1 − x), hence (1/k)x′′ = x′(1 − x) − xx′ = x′(1 − 2x). 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 so- lution x0, satisfying initial condition x0(0) = 1/2. Hence, due to x′′0(0) = 0, we have that x0 has an inflection at 0. � Similarly, any shifted solution xγ(t) = x0(t−γ) has an inflection point at t = γ, as xγ(γ) = x0(t−γ)|t=γ = x0(0) = 1/2. Problem. Find the value of γ so that the shifted function xγ(t) = x0(t − γ) satisfies the initial condition xγ(0) = x0(−γ) = 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 x0(−γ) = 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 + ekγ)−1 = x∗, hence γ = 1 k ln 1 −x∗ x∗ . Therefore tlag = γ − 2 k = 1 k ( ln 1 −x∗ x∗ − 2 ) . In order to have tlag ≥ 0, the restriction x∗ ≤ (e2 + 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∗ ≤ (e2 + 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 condi- tion in order to possess a positive lag time. Growth processes with positive lag time are typical for bio-chemical reactions, that is reactions involving Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 5 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth functional proteins, such as enzymes, receptors, ligands, as well as population models. Finally we shall point out one more charac- teristic 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 be- tween 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 pro- teins. 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.: S + E k1−→←− k−1 C k2−→ P + E. (18) It is assumed that the rate parameters k1,k−1,k2 are positive constants such that k1 > k−1. Re- action scheme (18) describes the reaction mech- anism between an enzyme E with a single ac- tive 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 sub- strate into a complex C having specific properties different than the properties of the free enzyme and thus being necessarily considered as a separate substance. Denoting the concentrations s = [S], e = [E], c = [C], p = [P ] and applying the mass action law to Henri’s reaction scheme (18) we obtain the following system of ODEs: ds/dt = −k1es + k−1c, de/dt = −k1es + (k−1 + k2)c, dc/dt = k1es− (k−1 + k2)c, dp/dt = k2c, (19) If the three rate constants k1,k−1,k2 are known, system (19) can be treated as an IVP with initial conditions s(0) = s0, e(0) = e0, c(0) = c0, p(0) = p0, (20) usually s0 > 0, e0 > 0, c0 = 0, p0 = 0. System (19) has two conservation laws: e′+c′ = 0 and s′ + c′ + p′ = 0. These two relations can be used to reduce the system to two equations: ds/dt = −k1(e0 − c)s + k−1c, dc/dt = k1(e0 − c)s− (k−1 + k2)c. (21) Proposition II.4. The solution p of the IVP (19)– (20) is a sigmoidal function. 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 k1 > k−1 and e ≤ e0,c ≤ Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 6 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth Fig. 4. Solutions to (19) with initial conditions s(0) = 0.5, e(0) = 1.0, c(0) = 0.0 p(0) = 0.2. c0 are bounded). Equation p′ = k2c ≥ 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∗) = k2c′(t∗) = 0. meaning that p has an inflection at t∗. From s′+c′+p′ = 0 we have s+c+p = s0 showing that p(t)t−→∞ = s0. 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 con- stants k1,k−1,k2 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 mod- els (10), (18) present two useful mechanisms for studying biochemical and biological growth. Var- ious combinations of these two models have been proposed for the study of the EPS production [27], [28], [34]. III. R-TYPE MODELS FOR AMILOID FIBRILLATION Recent intensive research into the physico- chemical 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 reac- tants: 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 reac- tion 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 aggre- gation [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 M + F k+−→←− k− C kc−→ F + F. (22) 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. Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 7 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth Denoting the concentrations m = [M], f = [F], c = [C] and applying the mass action law to reaction scheme (22) we obtain the following system of three ODE’s: dm/dt = −k+fm + k−c, df/dt = −k+fm + (k− + 2kc)c, dc/dt = k+fm− (k− + kc)c. (23) If the three rate constants k+,k−,kc are known, system (23) can be treated as a Cauchy ODE IVP with initial conditions m(0) = m0, f(0) = f0, c(0) = c0, p(0) = p0, (24) thereby m0 > 0, f0 > 0, c0 = 0, p0 = 0. A conservation law for system (23) is m′ + f ′ + 2c′ = 0 implying the relation m + f + 2c = const = m0 + f0 = a. We have m = a−f − 2c in system (23) reducing it to f ′ = −k+f(a−f − 2c) + (k− + 2kc)c, c′ = k+f(a−f − 2c) − (k− + kc)c. (25) From (25) we have f ′ + c′ = kcc. 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 ′ = kcc [1], [4], [5], [13], [32], [35], [36], [37]. Hence, for some t∗ ∈ ∆ we have f ′′(t∗) = kcc′(t∗) = 0, hence function f has inflection and has a sig- moidal form. In addition, from f ′ + c′ = kcc and c′(t∗) = 0, we obtain the slope at the inflection point: f ′(t∗) = kcc(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 ′(t∗) = kcc(t∗) > 0. The solutions m(t),f(t) of model (23) are vi- sualized on Figure 5 for the following input data k+ = 0.1,k− = 0.05; kc = 0.2; m0 = 10; f0 = 0.1; c0 = 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 pa- rameter 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 condi- tions by fitting (some of) the solutions to available time course experimental measurements. Compu- tational 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 net- works. All they involve an intermediate product P representing the aggregated monomer before turn- ing 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: M + F k+−→←− k− C kc−→ F + P, P kp−→ F. (26) Here the transition of the aggregated monomer P into fibril follows the decay-saturation mecha- nism (1). Applying the mass action law we obtain the ODE system dm/dt = −k+fm + k−c, df/dt = −k+fm + (k− + kc)c + kpp, dc/dt = k+fm− (k− + kc)c, dp/dt = kcc−kpp. (27) Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 8 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth Fig. 6. Graph of the solutions m,f,c,p of model (27) A conservation law for system (27) is m′ + f ′ + p′ + 2c′ = 0. The solutions of the system are visualized on Fig. 6. for the following data: k1 = 2.62; k2 = 0.62; k3 = 2.62; k4 = 0.63; m0 = 1; f0 = 0.05; c0 = 0; p0 = 0. Reaction network 2. An intermediate product P is involved as follows: M + F k+−→←− k− C kc−→ F + P kp−→ F + F. (28) The transition of the aggregated monomer P into fibril according to reaction network (28 fol- lows the Verhulst autocatalitic mechanism (10). Applying the mass action law we obtain the ODE system: dm/dt = −k+fm + k−c, df/dt = −k+fm + (k− + kc)c + kpfp, dc/dt = k+fm− (k− + kc)c, dp/dt = kcc−kpfp. (29) A conservation law for system (29) is m′+f ′+ p′ + 2c′ = 0. Reaction network 3. Here an intermediate product P is added as follows: M + F k+−→←− k− C kc−→ F + P k′−→←− k′′ Cp kp−→ F + 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+fm + k−c, df/dt = −k+fm + (k− + kc)c + k′′cp −k′fp + 2kpcp, dc/dt = k+fm− (k− + kc)c, dp/dt = kcc−k′fp + k′′cp, dcp/dt = k ′fp− (k′′ + kp)cp. (31) A conservation law for system (31) is m′+f ′+ p′ + 2c′ + 2c′p = 0. IV. CONCLUSION The proposed reaction networks (R-models) (22), (27), (29), (31) together with the implied re- action 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. ACKNOWLEDGMENT The author is grateful to Dr. N. Kyurkchiev for his careful reading, useful suggestions, stimulating discussions and technical help. REFERENCES [1] Alt, R., S. Markov, Theoretical and computational studies of some bioreactor models, Computers and Mathematics with Applications 64 (2012), 350–360. Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 9 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth [2] Anguelov, R., Markov, S.: Hausdorff Continuous In- terval Functions and Approximations, In: Nehmeier, M. et al. (Eds), Scientific Computing, Computer Arithmetic, and Validated Numerics (SCAN 2014), LNCS 9553, Springer, 2016, 3–13. DOI: 10.1007/978-3-319-31769-4 [3] Arosio, P., T. P. J. Knowles, S. Linse, On the lag phase in amyloid fibril formation, Physical Chem- istry Chemical Physics 17, 2015, 7606–7618. [4] Bersani, A. M., G. Dell’Acqua, Is there anything left to say on enzyme kinetic constants and QSSA?, J. Math. Chem. 50 (2012) 335–344. [5] Bersani, A. M., E. Bersani, G. Dell’Acqua, M. G. Pedersen, New trends and perspectives in nonlinear intracellular dynamics: one century from Michaelis- Menten paper, Continuum Mech. Thermodyn. 27 (2014) 659–684. [6] Chellaboina, V., S. P. Bhat, W. M. Haddat, D. S. Bernstein, Modeling and Analysis of Mass-Action Kinetics, IEEE Control Systems Magazine, August 2009, 60–78. DOI: 10.1109/MCS.2009.932926 [7] Costarelli, D., Spigler, R.: Approximation results for neural network operators activated by sigmoidal functions, Neural Networks 44, 101–106 (2013). [8] Costarelli, D., R. Spigler, Constructive Approxima- tion by Superposition of Sigmoidal Functions, Anal. Theory Appl. 29, 169–196 (2013), DOI: 10.4208/ata.2013.v29.n2.8 [9] Deichmann, U., S. Schuster, J.-P. Mazat, A. Cornish-Bowden, Commemorating the 1913 Michaelis-Menten paper Die Kinetik der Invertinwirkung: three perspectives, FEBS Journal 281 (2014), 435–463. (see Part 3: before Michaelis and Menten: Victor Henris equation by Jean-Pierre Mazat) DOI: 10.1111/febs.12598 [10] Dimitrov, S., G. Velikova, V. Beschkov, S. Markov, On the Numerical Computation of Enzyme Kinetic Parameters, Biomath Communications 1 (2), 2014, 6–29. [11] Dimitrov, S., S. Markov, Metabolic Rate Constants: Some Computational Aspects, Math. Comput. in Simulation, (2015). DOI: 10.1016/j.matcom.2015.11.003 [12] J. Dombi, Z. Gera, The approximation of piecewise linear membership functions and Likasiewicz oper- ators, Fuzzy Sets and Systems 154, 275–286 (2005). DOI: 10.1016/j.fss.2005.02.016 [13] Hanson, S. M., S. Schnell, The reactant stationary approximation in enzyme kinetics, J. Phys. Chem. A 112 (2008) 8654–8658. [14] Henri, V., Ueber das Gesetz der Wirkung des In- vertins. Z. Phys. Chem., 39 (1901), 194–216. [15] Henri, V., Recherches sur la loi de l’action de la sucrase, C. R. Hebd. Acad. Sci. 133 (1901) 891– 899. [16] Henri, V., Lois générales de l’action des diastases, Paris: Librairie Scientifique A. Hermann. 1903. [17] Kyurkchiev, N., A note on the new geometric repre- sentation for the parameters in the fibril elongation process, Compt. rend. Acad. bulg. Sci., 69 (8) (2016), 961–970. [18] Iliev, A., N. Kyurkchiev, S. Markov, On the approximation of the cut and step functions by logistic and Gompertz functions, BIOMATH 4 (2), 2015. DOI: 10.11145/j.biomath.2015.10.101 [19] Iliev, A., N. Kyurkchiev, S. Markov, On the ap- proximation of the step function by some sigmoid functions, Math. Comput. Simulation, (2015). [20] Kyurkchiev, N., On the Approximation of the step function by some cumulative distribution functions, Compt. rend. Acad. bulg. Sci., 68, (12), 1475–1482. [21] Kyurkchiev N., Markov, S.: Sigmoidal Func- tions: some Computational and Modelling Aspects, Biomath Communications 1 (2), (2014), 30–48. [22] Kyurkchiev N., S. Markov, On the approximation of the generalized cut function of degree p+1 by smooth sigmoid functions, Serdica J. Computing 9, 1, 2015, 93–104. [23] Kyurkchiev, N., S. Markov, On the Hausdorff dis- tance between the Heaviside step function and Ver- hulst logistic function, J. Math. Chem. 54(1), 109– 119 (2016). DOI:10.1007/S10910-015-0552-0 [24] Kyurkchiev, N., S. Markov, On the numerical so- lution of the general kinetic “K-angle” reaction system, J. Math. Chem. 54(3), 792–805 (2016). DOI:10.1007/S10910-016-0592-0 [25] Kyurkchiev, N., S. Markov, Sigmoid Functions: Some Approximation and Modelling Aspects. Some Moduli in Programming Environment MATHE- MATICA, LAP (Lambert Academic Publ.), Saar- brucken, 120 pp., 2015; ISBN 978-3-659-76045-7 [26] McKendrick, A. G., M. Kesava Pai, The Rate of Multiplication of Micro-organisms: A Mathematical Study, Proc. of the Royal Society of Edinburgh 31, 649–653 (1912). DOI: 10.1017/S0370164600025426 [27] Markov, S., On the mathematical modelling of microbial growth: some computational aspects, Serdica J. Computing 5 (2), 153–168 (2011). http://serdica-comp.math.bas.bg/ index.php/serdicajcomputing/article/ view/125 [28] Markov, S.: Cell Growth Models Using Reaction Schemes: Batch Cultivation, Biomath 2/2 (2013), 1312301 DOI: 10.11145/j.biomath.2013.12.301 [29] Meisl, G., JB Kirkegaard, P Arosio, TCT Michaels, M Vendruscolo, CM Dobson, S Linse, TPJ Knowles, Molecular mechanisms of protein aggre- Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 10 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 S. Markov, Building reaction kinetic models for amiloid fibril growth gation from global fitting of kinetic models, Nature Protocols 11, 2016, 252–272. [30] P. Mendes, S. Hoops, S. Sahle, R. Gauges, J. Dada, U. Kummer, Computational modeling of biochem- ical networks using COPASI, Methods Mol. Biol. 500 (2009) 17–59. [31] Michaels, Thomas C. T., Tuomas P. J. Knowles, Kinetic theory of protein filament growth: Self- consistent methods and perturbative techniques, In- tern. J. of Modern Physics B 29, 2 (2015) 1530002. World Scienti?c. DOI: 10.1142/S0217979215300029 [32] Murray J. D., Mathematical Biology: I. An Intro- duction, Third Edition, Springer, 2002. [33] Putz, M. V., A. M. Putz, Logistic vs. W-Lambert in- formation in quantum modeling of enzyme kinetics, Int. J. Chemoinf. Chem. Eng. 1, 42–60 (2011). [34] Radchenkova, N., M. Kambourova, S. Vassilev, R. Alt, S. Markov, On the Mathematical Modelling of EPS Production by a Thermophilic Bacterium, Biomath 3/1, 1407121, (2014). DOI: 10.11145/j.biomath.2014.07.121 [35] Schnell, S., P. K. Maini, Enzyme kinetics at high enzyme concentration, Bull. Math. Biol. 62 (2000), 483–499. [36] Schnell, S., P. K. Maini, A century of enzyme kinetics: Reliability of the Km and vmax estimates, Comments Theor. Biol. 8 (2003), 169–187. [37] Schnell, S., Validity of the Michaelis-Menten equation—Steady-state, or reactant stationary as- sumption: that is the question, FEBS J. 281 (2014), 464–472. [38] Shammas, S. L., et al., A mechanistic model of tau amyloid aggregation based on direct observation of oligomers, Nature communications 6, 2015. [39] Shoffner, S. K., S. Schnell, Estimation of the lag time in a subsequent monomer addition model for fibril elongation, BioRxive preprint, 2016. DOI: 10.1101/034900 [40] R. Tobias, G. Tasi, Simple algebraic solutions to the kinetic problems of triangle, quandrangle and pentangle reactions, J. Math. Chem. 54 (1), 85–99 (2016). DOI: 10.1007/s10910-015-0550-2 [41] Verhulst, P.-F., Notice sur la loi que la population poursuit dans son accroissement, Correspondance mathematique et physique 10: 113–121 (1838). [42] Verhulst, P.-F., Recherches mathematiques sur la loi d’accroissement de la population (Mathematical Researches into the Law of Population Growth In- crease), Nouveaux Memoires de l’Academie Royale des Sciences et Belles-Lettres de Bruxelles 18: 1–42 (1845). [43] Verhulst, P.-F., Deuxieme memoire sur la loi d’accroissement de la population, Memoires de l’Academie Royale des Sciences, des Lettres et des Beaux-Arts de Belgique 20: 1–32 (1847). Biomath 5 (2016), 1607311, http://dx.doi.org/10.11145/j.biomath.2016.07.311 Page 11 of 11 http://dx.doi.org/10.11145/j.biomath.2016.07.311 Introduction Biological growth/transition models: three case studies Case study 1: saturated growth Case study 2: Verhulst logistic growth Lag time Case study 3: growth models using Henri's reaction scheme R-type models for amiloid fibrillation A basic Verhulst-Henri model Three variants of the basic model Conclusion References