www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Reaction networks reveal new links between Gompertz and Verhulst growth functions Svetoslav Markov Institute of Mathematics and Informatics, Bulgarian Academy of Sciences smarkov@bio.bas.bg Received: 8 January 2019, accepted: 16 April 2019, published: 20 April 2019 Abstract—New reaction network realizations of the Gompertz and logistic growth models are pro- posed. The proposed reaction networks involve an additional species interpreted as environmental re- source. Some natural generalizations and modifi- cations of the Gompertz and the logistic models, induced by the proposed networks, are formulated and discussed. In particular, it is shown that the induced dynamical systems can be reduced to one dimensional differential equations for the growth (resp. decay) species. The reaction network formu- lation of the proposed models suggest hints for the intrinsic mechanism of the modeled growth process and can be used for analyzing evolutionary measured data when testing various appropriate models, especially when studying growth processes in life sciences. Keywords-Dynamical growth model; logistic func- tion; Gompertz function; sigmoidal function; dy- namical system; reaction network, first integral; conservation equation. I. INTRODUCTION Sigmoidal functions are an useful tool for mod- eling measurenent data for the study of evolution- ary growth processes in life sciences [9], [10]. When studying the time evolution of biological growth processes we are often given a set of measured data of the form (ti,yi), i = 1, ...,n, where yi is an experimentally obtained value at time moment ti. We then have to choose a model function y = f(t) that fits the measured data. More specifically, function f is chosen from a family of functions depending on some parameters and the fitting process consists in finding a suitable parameter set so that a good approximation (fit) is achieved. The definition of the family of modeling functions is a major challenge. To achieve a good fit, we need to choose a family that indicates (incorporates, reflects) the “mechanism” (law) of the physical process generating the experimen- tally measured data set. In practice, the modeling function f is chosen either from a family of explicitly defined functions, or f is defined as a solution to a class of dynamical systems. In many situations the intrinsic mechanism of the growth process is little or not known. In such situations, a good idea is to look for a dynamical model that consists of a system of reaction equations induced Copyright: c© 2019 Markov. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Svetoslav Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions, Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 1 of 14 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions by some chemical reaction network via reaction kinetic principles, such as mass action kinetics [22], [25], [5]. In this way we readily have a physico-chemical interpretation of the dynamical model and its constituents such as reacting species, rate parameters, order of interactions, sigmoidal shape of the graph of the growth variable, lag time, etc. What then remains is to give a meaning of the chemical terms relative to the studied real- world biological process. Reaction networks are well known for a number of basic dynamical growth models, such as the saturation-decay mod- els, epidemiological models, population predator- prey type models, demographic models, etc. Un- doubtedly, a reaction network formulation of a dy- namical model, if possible, contributes to a better understanding of the mechanism of the specific physical process and to possible improvements of the existing mathematical model. In the present work we propose a reaction net- work realization of the Gompertz model [8]. The proposed reaction network suggests an interpreta- tion of the reaction mechanism of the Gompertz model in terms of population dynamics theory and reveals its relation to the Verhulst logistic growth model. The proposed reaction network conforms with some of the already suggested interpretations based on the differential formulation of the model and recent studies in cancer research [7], [34]. In particular, our reaction mechanism involves two species, one for the population size/volume plus an additional species presenting the nutrient resources of the environment. In addition, some natural generalizations and modifications of the Gompertz model induced by the reaction network are proposed, discussed and compared to classical logistic and Gompertz models. The next section provides a brief introduction to reaction network theory, in order to fix notations and for the exposition to be more self-contained. Section III introduces the Gompertz models, while its reaction network realization and generalizations are discussed in Section IV. Section V considers reaction networks and respective models with lo- gistic decay of the resource species. It is followed by some notes on applications and concluding remarks in the last two sections of the paper. II. GROWTH MODELS AND REACTION NETWORKS We focus our attention on growth functions (models) formulated as solutions to differential equations or systems of differential equations. In the latter case we speak of dynamical system that may have several variables apart from the one considered as growth function. Without loss of generality we shall consider growth functions f(t) defined in the interval T = [0,∞) with nonnegative values f(t) ≥ 0. In many situations the dynamical system suggests some “insight” for the “inner mechanism” that controls the behavior of the solutions and for the physical meaning of the parameters involved in the system. The “mechanism” of the process is especially well pre- sented when the dynamical model has a realization in the form of a reaction network [5], [22]. A reaction network is a set of (elementary) reactions. A reaction is defined by a set of species that are either reactants (reagents) or products or both. For example, let a reaction have three species S, P , X and let species S, X be reactants, whereas P and X be products, then the reaction is written symbolically in the form: S + X k−→ P + X, (1) k > 0 denoting the rate of the reaction. Applying mass action (MA) kinetics, reaction (1) is uniquely “translated” into a system of three differential equations, one for the mass (concentration) of each species involved. The differential equations in the dynamical system are then called “reaction equations” (although this term is often used in a much broader sense in the literature). In our case, the reaction network (1) is translated as a dynamical system of three reaction equations for the corresponding masses (concentrations) s = s(t), p = p(t), x = x(t) of the species S, P , X: s′ = −ksx, p′ = ksx, x′ = 0, Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 2 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions where s′ = ds/dt, p′ = dp/dt, x′ = dx/dt denote derivatives with respect to the time variable t. The solutions of this system are the familiar exponential decay s, the saturation growth p and the constant catalyst x. On the same example we can demonstrate the converse process — passing from a dynamical system to a reaction network. This process is not unique, if at all possible. Indeed, from the third equation x′ = 0 we have x = const = c, so that the dynamical system reduces to two equations of the form s′ = −acs, p′ = acs. A possible realization of this system in terms of reaction network is then S → P , with a rate parameter equal to ac, which is distinct from the reaction network (1) as the catalyst species X is now missing. Let us note that the presence of species X can be important, especially if X participates in some other reaction(s). A. The logistic model Let us consider the Verhulst logistic model [33]. The logistic growth function is usually defined as solution to the differential equation x′ = kx(1 − x K ), (2) wherein k, K are positive parameters, resp. called (intrinsic) reaction rate and carrying capacity. For our purposes it is more convenient to con- sider (2) in the following form: x′ = kx(c−x), (3) k and c being positive parameters. The solutions of equations (2) and (3) coincide up to an affine transformation of the time variable of the form t∗ = at, a = const. The two solutions are same (a = 1) whenever K = c = 1. We can say that the two forms (2) and (3) are equivalent in the above sense. Equation (3) can be “recast” into two differential equations by introducing a new variable s = c−x. Then we have x′ = kxs, s′ = −x′ = −ksx, so, we can write (3) as the dynamical system: s′ = −ksx, x′ = ksx. (4) System (4) is equivalent to equation (3) in the sense that the solutions for x in both dynamical systems coincide. The form of system (4) suggests that the new variable s can be interpreted as mass of some new “intermediate” species S. Looking for a suitable reaction network involv- ing two species S, X with masses/concentrations resp. s, x, satisfying the dynamical system (4), we may consider the following reaction network [14]: S + X k−→ 2X, (5) where k > 0 is the reaction rate and 2X is an abbreviation of X +X. Indeed, applying the mass action principle to reaction network (5) we obtain the dynamical reaction system (4). Conversely, due to s′+x′ = 0, and consequently s + x = const = c > 0, (6) we have s = c−x, (7) which substituted in the second equation of (4) gives the differential equation (3). Assuming in (4) initial conditions s(0) = a, x(0) = b, we have s + x = c = a + b. Due to s > 0 and 0 < x < c, we have x′ = kx(c−x) > 0, showing that the solution x is monotone increasing and tends asymptotically to c = a+b with t →∞, thus justifying the interpretation of the number c as environmental carrying capacity. Species S can be interpreted as the resource (food) consumed (used, uptaken) by species X in order to reproduce itself. Note that species X appears on both sides of reaction network (5), playing thus simultaneously the roles of a reactant and a product, so X is a catalyst. As the catalysts X reproduces itself, reaction network (5) is called autocatalytic. For the second derivative of x we obtain x′′ = kx′(c−x) + kx(−x′) = kcx′ − 2kxx′ = kx′(c− 2x). This shows that, for x(0) = b < c/2 = (a + b)/2, Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 3 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions that is b < a, an inflection occurs when x = (a + b)/2. Thus, function x has a sigmoidal shape when increasing from x(0) < c/2 up to approaching asymptotically the value c. We wish to note again that every reaction net- work induces a differential reaction system in a unique way, but there is no uniqueness in the op- posite direction, that is there can be many distinct reaction networks that induce the same differential systems. The logistic model offers an illustrative example. Indeed, consider the reaction network: X k1−→←− k2 X + X. (8) This so-called “reversible reaction” consists of two elementary reactions that can be written equiva- lently as X k1−→ X +X, X +X k2−→ X. Applying MA kinetics, we obtain the differential equation x′ = k1x−k2x2, which is of the same type as equation (3). In such situations it is up to the modeler to choose the reaction network that offers a more adequate inter- pretation of the particular modeling situation. Note that reaction network (5) involves two species, whereas (8) uses only one, so it is important to decide how many species are involved in the phys- ical process. A dynamical process may involve intermediate species whose mass is zero at the be- ginning and at the end of the process. Such species remain hidden as they cannot be easily measured. Reaction network (8) is reversible, something con- sidered normal for chemical reactions where the so-called “principle of microscopic reversibility” takes place. On the other hand, reaction network (5) is irreversible, which is a normal situation when modeling complex biological systems such as organs, organisms and populations. A (chemical) reaction network formulation of a dynamical model, if possible, contributes to a deeper understanding of the mechanism of the par- ticular real-word biological process and to possible improvements of the dynamical model whenever necessary. In Section V we propose one more reaction net- work that induces the logistic differential equation (3). B. Generalized Verhulst growth models The theory of reaction networks offers a pow- erful tool not only to understand the mechanism of classical models, but also to construct various modifications of existing models in order to adapt the model better to particular real-world phenom- ena. We next illustrate this idea on a modification and generalization of the logistic growth model. The Verhulst model is generalized by the reac- tion network formulated in the following proposi- tion. Proposition 1. The autocatalytic reaction network: X + n∑ i=1 Si k−→ X + mX, (9) where n,m ≥ 1 are integers, induces the following dynamical growth model x′ = kmx n∏ i=1 ( ci − x m ) , (10) where ci > 0 are constants. Proof: Applying the MA law to reaction network (9) we obtain: s′i = −kx n∏ j=1 sj, i = 1, ...,n, (11) x′ = mkx n∏ i=1 si. (12) From (11)–(12) we have for every i = 1, ...,n s′i + x ′/m = 0, implying si + x/m = ci, where ci > 0 are constants. Substituting si = ci−x/m, i = 1, ...,n in (12) we obtain (10). Special cases. For n = 1, m = 1 we obtain the special case of reaction network (5) inducing Verhulst logistic equation (3). Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 4 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions For n = m = 2 we obtain the reaction network X + S1 + S2 k−→ X + 2X. (13) For n = m and Si = S, i = 1, ...,n, we obtain reaction network nS + X k−→ X + nX. (14) The special cases (13) and (14) are proposed in [1]. Discussion. A possible “biochemical” interpre- tation of model (9) is as follows: the model takes into account the interaction between several species, such as various types of foods and other environmental resources (water, air, light, etc). In chemistry it is unlikely that more than three species interact simultaneously [22]. However, in models related to biology and social sciences, this restriction can be relaxed. III. THE GOMPERTZ GROWTH MODEL: GENERAL NOTES The Gompertz growth model is used in numer- ous applications. For a review and classification of various formulations of the Gompertz model see e.g. [32]. Next, we briefly recall some familiar facts related to the Gompertz model. The Gom- pertz growth function is often defined as a solution x = x(t) to the differential equation x′ = kx(c− ln x), (15) where k > 0 and c are parameters. Similarly to the case of the logistic differential equation (3), equation (15) can be recast into a system of two differential equations, e.g. s′ = −ks, x′ = ksx, (16) as discussed in works related to a special class of dynamical systems, called S-systems [4], [29], [30]. It is remarkable that, similarly to the logistic case, a new variable s = s(t) appears in system (16), apart from the growth function x = x(t) in the single equation (15). Instead of system (16), some authors use a slightly different model: s′ = −s, x′ = ksx, (17) see e.g. [7]; other authors make use of the form s′ = −ks, x′ = sx, (18) see, e.g, [34]. In order to discuss the three slightly different model forms (16), (17) and (18), let us introduce some notation that is independent on the notation used in the available literature. Each one of the three systems (16), (17) and (18), consists of two differential equations. One of these equations is uncoupled, involving only one unknown function, which is exponentially decreasing. This function will be called decay function or briefly d-function, in our case, this is function s = s(t). The uncoupled equation for the d-function will be called decay equation or d-equation. The other equation, named growth equation or g-equation, involves, apart from the d- function, one more function, named growth func- tion or g-function, which is increasing in time. In this work the g-function is denoted by x. Let us now discuss the differences between the three systems (16), (17) and (18). System (16) involves two equal rate parameters in both equations, system (17) has a fixed rate equal to 1 in the d-equation and a (free) rate parameter k > 0 in the g-equation. System (18) has a (free) rate parameter in the d-equation and a fixed rate 1 in the g-equation. Using Gompertz-type models of the forms (17) and (18) one accepts the possibility of distinct rate parameters. In the present work we also adopt this assumption and shall explicitly denote the two rate parameters by different symbols, namely k1,k2: s′ = −k1s, x′ = k2sx. (19) Dynamical system (19) induces the relation s′/k1 + x ′/(k2x) = 0, which can be written as s′/ρ+x′/x = 0, ρ = k1/k2. Integrating (obtaining a first integral) leads to the “conservation” relation: s ρ + ln x = const = c. (20) Relation (20) gives an expression for the variable s in terms of the variable x, namely s = ρ(c− ln x). (21) Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 5 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions Using relation (21) we can establish a relation between the parameters of the single Gompertz equation (15), and those of the 2D-system (19). Let functions s, x be solutions to system (19). Hence, (21) holds. Substituting s = ρ(c − ln x) in the equation for x in (19) we obtain x′ = k2xs = k2xρ(c− ln x) = k1x(c− ln x). We thus obtain Gompertz equation (15) with pa- rameters k1 and c. If initial conditions s(0) = a, x(0) = b are given, then c = a/ρ + ln b, so that the Gompertz equation resulting from the system (19) with rate constants k1,k2 is x′ = k1x(a/ρ + ln b− ln x) or, equivalently, x′ = k1x(a/ρ + ln b/x). Conversely, let g-function x be solution to (15) and x(0) = b. Define function s = ρ(c − ln x), where ρ > 0 and c are parameters to be deter- mined, so that: x′ = k ρ xs, and, s′ = ρ(c− ln x)′ = −ρ x′ x = −ρ kxs x = −ks. Hence, s,x satisfy (19) with k1 = k and k2 = k/ρ. Note that for fixed parameters c,x(0) the initial condition s(0) for the decay function is determined from s(0) = ρ(c− ln b). The above discussion on the relation between the Gompertz equation (15) and system (19) can be summarized in the following proposition: Proposition 2 [1]. Let functions s,x be solutions to system (19), with initial conditions s(0) = s0, x(0) = x0, then x is a solution to a differ- ential equation of the form (15), with parameters k = k1 and c = k2s(0)/k1 + ln(x(0)) and initial condition x(0) = x0. Conversely, if function x is a solution to (15), x(0) = x0, then for any ρ > 0 functions s,x, where s = ρ(c − ln x), satisfy a dynamical system of the form (19) with k1 = k and k2 = k/ρ. Remarks. 1) Relation (21) gives an expression for the variable s in terms of variable x. Relation (21) reminds of the analogous relation (7) for the logistic resource variable, resp. the interpretation of the parameter c as an envi- ronmental carrying capacity. 2) From the d-equation s′ = −k1s we have s(t) = s(0)e−k1t. Then from the second equation in (19) we have x′ x = k2s = k2s(0)e −k1t, an expression known as Gompertz law of mortality. Some authors call the expression of the form x′/x per capita population growth rate. 3) The second part of Proposition 2 shows that representing the Gompertz equation (15) as a system of two differential equations (19) introduces one free parameter. It can be taken as ρ or as a = s(0), either one being a function of the other through a = ρ(c−ln b). 4) The appearance of a new variable (the d- function s) in the Gompertz model (19) sug- gests that this variable is mass/concentration of some particular species S depending on the modelling situation. The biological meaning of species S in systems of the form (19) is much discussed in the literature, see e.g. [7], [34]. 5) System (19) belongs to the class of “S- systems” [2], [3], [4], [29], [30], it can be considered as “recast” form of equation (15). In the literature on S-systems one can find theorems that generalize Proposition 2. Our next aim is to consider the d- and g- functions s = s(t), x = x(t) as masses (con- centrations) of two reacting species S, X and to formulate reaction network(s) involving the two species S, X, that generate system (19) under the MA principle. As shown above, Gompertz model can be for- mulated equivalently in several forms. We shall next focus on the Gompertz model in the form (19) as expressing reactions between two species. Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 6 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions So let us look at the variables s,x as masses (con- centrations) of two species S, resp. X. The first equation (s′ = −k1s) of system (19) indicates no interaction between species S and X, but merely an exponential decay of S with a rate parameter k1. The second equation (x′ = k2sx) suggests that S and X interact as S+X. Due to the independent decay of S, the latter observation may lead to the conclusion that no realization of the Gompertz model as a reaction network is possible, as stated in [1]. However, there exists such a realization as formulated in the next section. IV. REACTION NETWORK REALIZATION OF THE GOMPERTZ MODEL AND GENERALIZATIONS A. Main result The Gompertz model can be formulated by means of a reaction network as follows. Proposition 3. The reaction network involving species S,X,P : S k1−→ P S + X k2−→ 2X + S (22) induces Gompertz reaction equations (19), resp. Gompertz differential equation (15) for the masses/concentrations s,x of species S,X. Proof: Applying the mass action law to reaction network (22) yields the dynamical system: s′ = −k1s, p′ = k1s, x′ = k2sx. (23) System (23) incorporates system (19) plus an additional reaction equation for the by-product P (p′ = k1s); the first two reaction equations for s and p are uncoupled representing a saturation- decay mechanism. Thus the reaction for p can be either ignored or used depending on the particular modeling situation. Hence, system (23) is equiv- alent to system (19) as far as species s and x are concerned. As we know, system (19) induces Gompertz differential equation (15) in the sense of Proposition 2. This proves the proposition. Remark. In the special case k1 = 1, we obtain dynamical system (17), if k2 = 1 we obtain system (18), and if k1 = k2 we obtain dynamical system (19). Discussion. The second reaction in network (22) is based on the logistic reaction network (5) with the modification that species S catalyzes the repro- duction of species X. In certain applications when the mass of S is much greater than the mass of X, the reaction network says that the exhaustion of S is due not only to species X, but also to other factors leading to transformation of S into a third species P . In other applications it may be the case that species X uses S as a resource (food) for its reproduction but takes care to sustain the mass of S. Depending on the modeling situation reaction S → P can be replaced by some other suitable d-reaction, e.g. S → Ø meaning an outflow of resource S from the modeled system. Other possibilities of the decay mode of S, such as a logistic decay (S +Z → 2Z) will be explored later. B. Generalized Gompertz-type models based on reaction networks The second reaction in (22) unifies two pro- cesses: a) reproduction of species X at the expense of an uptake of the resource S, and b) recover of species S by a reproduction process catalysed by species X. These two processes can be formulated separately as follows: S k1−→ P S + X k2−→ 2X S + X k2−→ 2S + X (24) Indeed, reaction network (24) induces the follow- ing dynamical system: s′ = −k1s−k2sx + k2sx = −k1s, p′ = k1s, x′ = k2sx, which is identical with (23). Discussion. Reaction network (24) involves the logistic reaction network (5): S + X k2−→ 2X. In Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 7 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions the logistic model species X consumes species S and reproduces itself. The consumption of S is compensated by the third reaction: S + X k2−→ 2S + X , indicating that species X catalyzes the reproduction of species S. Note that instead of reaction S → P , a reaction S → Ø can be used, depending on the real-word modeling situation. From chemical point of view the rate parameters in the reactions with S + X in the left-hand side should be equal, however in certain modeling situations the rate parameters may be considered distinct, which leads to a more flexible generalized model of Gompertz type. More specifically the following model generalizes model (22): S k1−→ P S + X k2−→ 2X S + X k3−→ 2S + X (25) Indeed, reaction network (25) coincides with (24) if k3 = k2. The above generalization (25) of reaction net- work (22) can be naturally extended by replacing the logistic reaction in (25) by the generalized logistic reaction (9) as follows: Proposition 4. The reaction network involving species S,S1,S2, ...,Sn,P,X: S k1−→ P X + n∑ i=1 Si k2−→ X + mX S + X k3−→ 2S + X (26) where n,m are positive integers, generalizes Gom- pertz model (22). Proof: Follows from Proposition 1 and assuming k3 = k2. Depending on the modeling situation the gener- alized logistic reaction network (9): X + n∑ i=1 Si k2−→ X + mX used in (26) can be replaced by some special case such as (13), (14). 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 t 0.0 0.2 0.4 0.6 0.8 1.0 x( t) c = 1.0, c1 = 0.5, c2 = 0.5, x(0) = 0.1 dx/dt = kGx( lnx), kG = 3.50 dx/dt = kMx(c1 x2 )(c2 x 2 ), kM = 19.90 dx/dt = kLx(c x), kL = 8.95 _c- Fig. 1. Graphs of the solutions to models (3), (13) and (15). The rate constants and initial conditions are equal in order to compare the different sigmoidal shapes of the growth functions. Fig. 1 represents the graphs of the solutions to models (3), (13) and (15). The rate constants and initial conditions are chosen to be equal in order to compare the shapes of the growth functions. C. A mixed Verhulst-Gompertz model with decay- type resource uptake Our approach to formulate Verhulst and Gom- pertz models in terms of reaction networks reveals an important link between these two classical growth models. In order to pass from the Verhulst logistic model (5) to Gompertz model (22) we need to perform two steps in the modification of the reaction networks: a) add a reaction S k1−→ P , and b) add a species S in the right-hand side of (5) in order to obtain the second reaction (S + X k2−→ 2X + S) of the Gompertz reaction network (22), providing thus a sustainability of the resource species S. This observation suggests to take a closer look at the first step a) in the above modification process, namely adding a decay reaction to the logistic one, leading to the reaction network: S k1−→ P, S + X k2−→ 2X. (27) Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 8 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions Proposition 5. Reaction network (27) involving species S,P,X induces the following “mixed Verhulst-Gompertz” dynamical system for the masses s,p,x: s′ = −k1s−k2xs, p′ = k1s, x′ = k2sx, (28) where k1,k2 are positive parameters. Dynami- cal system (28) generates the following mixed Verhulst-Gompertz differential equation for the g- function x: x′ = k2x(c−x− ln xρ), ρ = k1/k2, (29) where c = s(0) + x(0) + ln x(0)ρ. Proof: Reaction network (27) induces (28) accord- ing to mass action kinetics. From (28) we obtain s′ + x′ = −k1s = −k1 x′ k2x . Hence, s′ + x′ + ρ x′ x = 0, where ρ = k1/k2. Integrating yields the “conser- vation” relation s + x + ln xρ = const = c, (30) wherein c = s(0) + x(0) + ln x(0)ρ. Relation (30) allows us to express s in terms of x: s = c−x− ln xρ. Substituting this expression in equation x′ = k2xs we obtain (29). This proves the Proposition. Equation (29) suggests that the reaction mecha- nism (27) presents an intermediate step between the logistic and the Gompertz growth models. More specifically, the main steps in the construc- tion of the Gompertz reaction mechanism, starting from the logistic one, are as follows. A construction of the Gompertz model in three steps starting from the logistic model. The reaction network approach offers a simple presentation of the consecutive steps leading from the logistic to the Gompertz model. Here is a list of three “elementary” steps in the construction of the Gompertz model, starting from the logistic one. Step 1: S + X k−→ 2X. Logistic model: x′ = kx(c−x). Step 2: S k1−→ P, S + X k2−→ 2X. Mixed V-G model: x′ = k2x(c−x− ln xρ). Step 3: S k1−→ P, S + X k2−→ 2X + S. Gompertz model: x′ = k1x(c− ln x). The above steps suggest a variety of combina- tions of different mechanisms that can replace the “elementary” reactions, for example the exponen- tial decay reaction S k1−→ P can be replaced by logistic decay. Fig. 2 represents the graphs of the solutions to models (3), (27) and (15). The rate constants and initial conditions are chosen to be equal in order to compare the shapes of the growth functions. It is observed that the solution of the mixed Verhulst model is “between” the solutions of the logistic and the Gompertz models. Fig. 2. Graphs of the solutions to models (3), (27) and (15). The rate constants and initial conditions are equal in order to compare the different sigmoidal shapes of the growth functions. D. Another reaction network realizations of the Verhulst logistic model In Section 3 we considered the reaction network (5) generating the logistic model (3). Below we propose another reacrion network that generates the logistic model. Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 9 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions Proposition 6. The reaction network involving species S,X,P : S + X k1−→ P + X S + X k2−→ 2X + S (31) induces the following dynamical system for the masses/concentrations s,x of species, resp. S,X: s′ = −k1sx, x′ = k2sx, (32) where k1,k2 are positive rate parameters. Dynam- ical system (32) generates the Verhulst differential equation (2) for the growth function x. Proof: Applying the mass action law to reaction network (31) yields the dynamical system: s′ = −k1sx, p′ = k1sx, x′ = k2sx. The above system incorporates system (32) plus an additional reaction equation for a by-product P (p′ = k1sx). The latter equation is uncoupled and can be ignored obtaining thus system (32). Due to s′/k1 + x′/k2 = 0, and consequently s/k1 + x/k2 = const, we can write s + ρx = const = c > 0, ρ = k1/k2, or s = c−ρx. (33) The constant c can be determined from the initial condititions for s,x, namely c = s(0) + ρx(0). Substituting the expression (33) in x′ = k2xs gives the logistic differential equation: x′ = k2x(c−ρx). Formally the above equation obtains the form (2) for c = 1,K = 1/ρ. This proves the proposi- tion. Remark. From chemical point of view it is un- common that the two reactions in (31) have dis- tinct rate parameters, providing that the reactants involved (S,X) are the same. However from bi- ological perspective the two reactions in network (31) can be considered to be of different nature, hence they may have distinct rates. In addition, reaction network (31) has some methodological value, demonstrating the difference between the logistic and the Gompertz models. Namely, the difference consists in the decay equation, showing the independence of the decay reaction S k1−→ P , resp. of the resource species S, on X in the case of the Gompertz model towards the dependant (catalyzed) decay S + X k1−→ P + X of species S on species X in the case of the logistic model. V. REACTION NETWORK MODELS WITH LOGISTIC DECAY FUNCTION FOR THE RESOURCE SPECIES As shown above the logistic reaction network in the form (31) has two reactions: i) a growth reac- tion S + X−→2X + S inducing a logistic growth of the g-species X, and ii) a decay reaction of the form S + X −→ P + X representing a logistic decay (consumption, uptake, outflow) of the d- species (resource, food) S. For comparison, the mixed logistic reaction network (27) involves a d- reaction describing an exponential decay: S−→P , resp. s′ = −ks, inducing the exponential solution s(t) = s0 exp(−kt). A logistic decay process results from a logistic d-reaction of the form S + Z−→2Z. Such a process represents a consumption of the resource S by a competitive autocatalytic species Z inducing thereby a sigmoidal “logistic- type” decay function s(t). We next apply the logistic decay mode of the d-function to the mixed Verhulst model (27). A. Reaction network model using logistic reac- tions for the g- and d-species Consider the reaction network S + Z k1−→ 2Z, S + X k2−→ 2X. (34) Proposition 7. Reaction network (34) involv- ing species S,Z,X induces the following mixed Verhulst-Gompertz differential equation for the g- function x: x′ = k2x(c−x− c1xρ), (35) Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 10 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions where ρ = k1/k2, c,c1 = const, or, equivalently, x′ = k2x exp(−k1t)/(1 + exp(−k1t), (36) In (35) and (36, k1,k2 are the positive rate param- eters in (34). Proof: Reaction network (34) induces the follow- ing dynamical system: s = −k1sz −k2sx, z = k1sz, (37) x = k2sx. Adding the three equations of dynamical system (37) we obtain: s′ + z′ + x′ = 0. Hence, s+z +x = const = c = s(0) +z(0) +x(0). (38) From the equations for z′ and x′ in (37), we obtain the relation z′ k1z = x′ k2x (= s), Hence, ln z k1 − ln x k2 = const = c∗, yielding ln(zk2/xk1 ) = c, or zk2/xk1 = const = c∗∗. We thus obtain zk2 = c∗∗xk1 , or z = c∗∗xρ, ρ = k1/k2. Substituting z in (38), we obtain s = c−x−z = c−x− c∗∗xρ. Substituting this expression for s in the third equation of (37) we obtain differential equation (35). To prove the non-autonomous equation (36) it is enough to note that the sigmoidal logistic-type uptake of the resource S by species Z does not depend on x. Remark. Model (34) is symmetrical with re- spect to species X and Z in the sense that both species can exchange their places in the reactions. Due to this symmetry the equations for the g- function z are similar to those of g-function x. B. Reaction network model using Gompertz-type growth reaction and Verhulst-type decay reaction As shown above the Gompertz reaction network model has two important constituents: i) a growth reaction S + X−→2X + S inducing a logistic growth of the g-function, and ii) an additional reaction of the form S −→ P representing the decay (consumption, uptake, outflow) of the re- source (food) species S. In the Gompertz reac- tion network the form of the d-reaction describes an exponential decay: s(t) = s0 exp(−kt). We next suggest that the d-reaction mat represent a consumption of the resource S by a competitive species following a sigmoidal “logistic-mode”, that is the d-species S is consumed by an autocat- alyst Z according to the d-reaction S + Z−→2Z. We shall next apply this consumption mode to the mixed Verhulst-Gompertz model (27). Consider the reaction network S + Z k1−→ 2Z, S + X k2−→ 2X + S. (39) Proposition 8. Reaction network (39) involving species S,Z,X induces the following logistic dif- ferential equation for the g-function x: x′ = k2x(c−µxρ), (40) where k1,k2 are the positive rate parameters in (39), ρ = k1/k2 and c,µ = const. Proof: Reaction network (39) induces the follow- ing dynamical system: s = −k1sz, z = k1sz, (41) x = k2sx. Adding the first two equations of dynamical sys- tem (41) we obtain: s′ + z′ = 0. Hence, s + z = const = c. From the equations for z′ and x′ in (41) we obtain the relation z′ k1z = x′ k2x (= s), Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 11 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions Hence, ln z k1 = ln x k2 + ln µ, µ = const, yielding z(1/k1) = µx(1/k2)), or z = µxρ. Substituting z in s = c−z, we obtain s = c−z = c−µxρ, which substituted in third equation of (41), gives the differential equation (40). Remark. Model (39) is another more general for- mulation of the logistic model. Indeed, for k1 = k2 equation (40) is similar to the logistic differential equation (2). Another interesting conclusion is that all forms of the logistic reaction equation can be written in the form: s′ = −k1s(c−s), x′ = k2sx. (42) The above “recast” form of the logistic equation shows that the d-equation is uncoupled from the g- equation and its solution can be given in the form: s(t) = c 1 + ek1(t−γ) , where γ is determined by the initial condition. Then from the second equation in (42) we obtain x′ x = k2c 1 + ek1(t−γ) (43) Formula (43) can be interpreted as “logistic mor- tality law” or “logistic per capita population growth rate”. VI. NOTES ON APPLICATIONS Gompertz-type models have been widely used to simulate the kinetics of various natural phenomena such as the growth of various species (microorgan- isms, animals, plants) and bio-products formation, e.g. methane [31], hydrogen [6], [24], etc. Numer- ous research articles are devoted to application of Gompertz-type models in tumor growth [7], [26], [20], [21], [34]. Gompertz models find application in software reliability models (GSRM, GMSRM) [28], [27], [23]. Gompertz-type modeling and data fitting problems stimulate various mathematical and computational studies, such as applications of new cumulative distribution functions, transfor- mations (type I–type III) to construct families of sigmoidal functions and new activation functions using ’correcting amendments’ [17], [18], [19]. For other related works the reader may consult [11]–[14]. VII. CONCLUDING REMARKS Modeling growth processes in life sciences is an important scientific area [9], [10]. Dynami- cal growth processes are often described by sig- moidal growth functions, such as saturation, logis- tic, Gompertz, etc., many of them being solutions to dynamical systems. Reaction network theory ia an important tool for generating dynamical growth models providing thereby useful interpretation of the intrinsic mechanism of the biological process. In this work we focus our attention on logistic and Gompertzian-type growth models from the perspective of the reaction networks theory. There are well-known reaction network realizations for a number of dynamical growth models, such as saturation, logistic, epidemic, etc., however, to our knowledge no such realization is known for the Gompertz model. In this paper a reaction network inducing the Gompertz model is proposed. The proposed reaction network involves an additional reaction for the uptake of the resource species. We also propose several reaction networks inducing dynamical models that generalize the Gompertzian one. Discussed are important links between the Gompertz and the logistic model. Our method can be considered as an extension of the work of previous authors who “recast” the Gompertz differential equation into a dynamical system of two differential equations involving thereby an additional variable (species) that can be interpreted as “resource” or “food” consumed by the growth variable (species). Discussed is also how the in- duced dynamical systems can be reduced to one- dimensional differential equations for the growth (resp. decay) species, by finding a first integral leading to a conservation equation.. The proposed reaction networks are simple and may seem trivial, Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 12 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions but are of some importance to those who construct new models to study biological growth processes whose underlying mechanism is unknown. The proposed reaction network realization of Gompertz growth model can be interpreted from the per- spective of demographic and socio-economic sci- ences. It is remarkable that the Gompertz reaction network comprises a reaction equation describing biological activity that is characteristic for highly organized biological organs, organisms or popu- lations. This explains why using the Gompertz model in demographic studies and cancer research is so successful. The reaction network approach clearly explains the close links between the Gom- pertz model and the Verhulst logistic model. Acknowledgments. The author is grateful to Dr. N. Kyurkchiev and Dr. M. Borisov for stimulating discussions on the subject and numerous encour- agements. The author is grateful to the anonymous reviewer for his careful reading of the manuscript, critical remarks and suggested improvements. REFERENCES [1] Anguelov, R., Borisov M., Iliev A., Kyurkchiev N., Markov S., On the chemical meaning of some growth models possessing Gompertzian- type property. Math Meth Appl Sci. 2017;112 https://doi.org/10.1002/mma.4539 [2] Bajzer Z., Marusic M, Vuk-Pavlovic S., Concep- tual frameworks for mathematical modelling of tumor growth dynamics. Math Comput Model. 1996;23:31-46. [3] Bajzer, Z., Vuk-Pavlovic S., New dimensions in Gom- pertzian growth. J Theor Med. 1999;2:307-315. Journnl of Theorericnl Medicine 2, 307-315 (1999) [4] Bajzer Z., Gompertzian growth as a self-similar and allometric process. Growth Develop Aging. 1999;63:3- 11. [5] Chellaboina V., Bhat S.P., Haddat W.M., Bernstein D.S. Modeling and analysis of mass-action kinetics. IEEE Contr Syst Mag. 2009;29:60-78. [6] Costa J.C., Barbosa S.G., Alves M.M., Sousa D.Z. (2012) Thermo chemical pre- and biological co- treatments to improve hydrolysis and methane produc- tion from poultry litter. Bio resource Technology 111: 141-147. [7] Gerlee, P. (2013) The Model Muddle: In Search of Tu- mor Growth Laws. Cancer research. 73. 10.1158/0008- 5472.CAN-12-4355. [8] Gompertz B., On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philos Trans R Soc London. 1825;115: 513–585. [9] Garfinkel, A., Shevtsov, J., Guo, Y., Modeling Life. The Mathematics of Biological Systems, Springer, 2018. [10] Goriely, A. The Mathematics and Mechanics of Biolog- ical Growth, Springer, 2017. [11] Iliev A., Kyurkchiev N., Markov S., On the approxi- mation of the cut and step functions by logistic and Gompertz functions. Biomath. 2015;4:2-13. [12] Kyurkchiev N., Markov S., Sigmoidal functions: some computational and modelling aspects. Biomath Com- mun. 2014;1:30-48. [13] Kyurkchiev N., Markov S., On the Hausdorff distance between the Heaviside step function and Verhulst logis- tic function. J Math Chem. 2016;54:109-119. [14] Kyurkchiev N., Markov S., On the numerical solution of the general kinetic ”K-angle” reaction system. J Math Chem. 2016;54:792-805. [15] Kyurkchiev, N., The new transmuted C.D.F. based on Gompertz function, Biomath Communications, 5 (1), 2018. [16] Kyurkchiev, N. Investigations on a hyper-logistic model. Some applications, Dynamic Systems and Applications 28(2) (2019), 351–369. [17] Kyurkchiev, N., A new class of activation func- tion based on the correcting amendments, Int. J. for Sci., Res. and Developments, 6 (2), 2018, 565-568; http://www.ijsrd.com/articles/IJSRDV6I20203.pdf [18] Kyurkchiev, N., Iliev, A., Extension of Gompertz-type equation in modern science. 240 Anniversary of the birth of B. Gompertz, LAP LAMBERT Academic Publishing, 2018, ISBN: 978-613-9-90569-0; [19] Kyurkchiev, N., A. Iliev, 240-th Anniversary of the Birth of Benjamin Gompertz, International Journal of Pure and Applied Mathematics, 120 (2), 2018, 223-227, ISSN: 1311-8080, [20] Laird A.K., Tyler S.A., Barton A.D., Dynamics of normal growth. Growth. 1965;29:233-248. [21] Laird, A.K., Dynamics of growth in tumors in normal organisms. Natl Cancer I Monogr. 1969;30:15-28. [22] Lente G., Deterministic kinetics in chemistry and sys- tems biology. In: Briefs in Molecular Science. Cham Heidelberg New York Dordrecht London: Springer; 2015. [23] Markov, S., A. Iliev, A. Rahnev, N. Kyurkchiev, On the exponential–generalized extended Compertz cumulative sigmoid, International Journal of Pure and Applied Mathematics, 120, No. 4 (2018). [24] Mu Y., Yu H.Q., Wang G. (2007) A kinetic approach to anaerobic hydrogen-producing process. Water Research 41(5): 1152-1160. [25] Murray J.D., Mathematical Biology: I. An Introduction. 3rd edn. New York Berlin Heidelberg: Springer-Verlag; 2002. Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 13 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 S. Markov, Reaction networks reveal new links between Gompertz and Verhulst growth functions [26] Norton L., A Gompertzian model of human breast cancer growth. Cancer Res. 1988;48(24 Pt1):7067-7071. [27] Pavlov, N., A. Iliev, A. Rahnev, N. Kyurkchiev, Some software reliability models: Approximation and mod- eling aspects, LAP LAMBERT Academic Publishing, 2018, ISBN: 978-613-9-82805-0. [28] Pavlov, N., G. Spasov, A. Rahnev, N. Kyurkchiev, A new class of Gompertz–type software reliability models, International Electronic Journal of Pure and Applied Mathematics, 12 (1), 2018, 43-57. [29] Savageau M.A., Allometric morphogenesis of complex systems: a derivation of the basic equations from first principles. P Natl Acad Sci. 1979;76:5413-5417. [30] Savageau, M.A., Voit EO., Recasting nonlinear differen- tial equations as S-systems: a canonical nonlinear form. Math Biosci. 1987;87:83-115. [31] Shen, J., J. Zhu, Development of General Gompertz Models and Their Simplified Two-Parameter Forms Based on Specific Microbial Growth Rate for Microbial Growth, Bio-Products and Substrate Consumption, Adv Biotech & Micro 4(3): AIBM.MS.ID.555640 (2017). [32] Tjorve K.M.C., Tjorve E. (2017) The use of Gompertz models in growth analyses, and new Gompertz-model approach: An addition to the Unified-Richards family. PLoS ONE 12(6): e0178691. https://doi.org/10.1371/journal.pone.0178691 [33] Verhulst P.F., Notice sur la loi que la population poursuit dans son accroissement. Corresp Math Phys. 1838;10:113-121. [34] West, J., Z. Hasnain, P. Macklin, P. K. Newton, An Evolutionary Model of Tumor Cell Kinetics and the Emergence of Molecular Heterogeneity Driving Gompertian Growth, SIAM Rev Soc Ind Appl Math. 2016 ; 58(4): 716736. https://doi.org/10.1137/15M1044825. Biomath 8 (2019), 1904167, http://dx.doi.org/10.11145/j.biomath.2019.04.167 Page 14 of 14 http://dx.doi.org/10.11145/j.biomath.2019.04.167 Introduction Growth models and reaction networks The logistic model Generalized Verhulst growth models The Gompertz growth model: general notes Reaction network realization of the Gompertz model and generalizations Main result Generalized Gompertz-type models based on reaction networks A mixed Verhulst-Gompertz model with decay-type resource uptake Another reaction network realizations of the Verhulst logistic model Reaction network models with logistic decay function for the resource species Reaction network model using logistic reactions for the g- and d-species Reaction network model using Gompertz-type growth reaction and Verhulst-type decay reaction Notes on applications Concluding remarks References