www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE The Gompertz model revisited and modified using reaction networks: Mathematical analysis Svetoslav M. Markov Institute of Mathematics and Informatics Bulgarian Academy of Sciences smarkov@bio.bas.bg The paper is dedicated to Prof. Dr. Roumen Anguelov on the occasion of his 65th birthday. Received: 22 August 2021, accepted: 2 October 2021, published: 4 October 2021 Abstract— In the present work we discuss the usage of the framework of chemical reaction net- works for the construction of dynamical models and their mathematical analysis. To this end, the process of construction of reaction-network-based models via mass action kinetics is introduced and illustrated on several familiar examples, such as the exponential (radioactive) decay, the logistic and the Gompertz models. Our final goal is to modify the reaction network of the classic Gompertz model in a natural way using certain features of the exponential decay and the logistic models. The growth function of the obtained new Gompertz-type hybrid model possesses an additional degree of freedom (one more rate parameter) and is thus more flexible when applied to numerical simulation of measurement and experimental data sets. More specifically, the ordinate (height) of the inflection point of the new generalized Gompertz model can vary in the interval (0, 1/e], whereas the respective height of the classic Gompertz model is fixed at 1/e (assuming the height of the upper asymptote is one). It is shown that the new model is a generalization of both the classic Gompertz model and the one-step exponential decay model. Historically the Gompertz function has been first used for statistical/insurance purposes, much later this function has been applied to simulate biological growth data sets coming from various fields of science, the reaction network approach explains and unifies the two approaches. Keywords-Systems of ordinary differential equa- tions, Reaction networks, Chemical reaction net- works, Evolutionary growth-decay models, Relative growth rate, Exponential (radioactive) decay, Log- arithmic change rate, Logistic model, Gompertz model. I. INTRODUCTION: REACTION NETWORKS AND EVOLUTIONARY GROWTH-DECAY MODELS The present study is devoted to mathemati- cal models induced by chemical reaction net- works describing evolutionary changes of biolog- Copyright: © 2021 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 M. Markov, The Gompertz model revisited and modified using reaction networks: Mathematical analysis, Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 1 of 21 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.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... ical species. We are interested in smooth changes that are described either as monotone growth, or monotone decay, in certain (large) time intervals, oscillation processes will not be considered in this work. We briefly call such changes, resp. functions, growth-decay changes/functions. The exponential radioactive decay, the Gompertz and the logistic models are three textbook examples, where we can consider the state variables involved as species that react between each others or are catalysts of certain reaction(s). The main purpose of this work is to propose a new modification of the classic Gompertz model possessing more flexibility and functionalities. On the way to achieve this goal we offer a brief introduction into the method of chemical reaction networks, which turns to be essential in the con- struction of new meaningful mathematical models, in particular when it comes to biological growth and decay processes. In the preliminary part of this work we present several examples demonstrating the role of the chemical reaction network theory (CRNT) in trac- ing the characteristics of elementary familiar math- ematical models for the numerical simulation of complex phenomenological (biological) processes. The example with the Gompertz model demon- strates the need of a detailed knowledge on the basics of elementary reaction networks. Biological growth-decay functions, describing evolutionary processes, are often presented in the literature by means of explicit algebraic expres- sions. Such a presentation offers little or no infor- mation on the physico-chemical mechanism of the studied process. More information in this direction is provided when the growth-decay functions are defined as solutions to systems of ordinary differ- ential equations. In the latter case we may look for a possible (chemical) reaction network, which implies the particular dynamical system via mass action kinetics [20]. If such a network does exist, we say that the differential system has a realization (formulation) in terms of a (chemical) reaction network [6]. Models formulated in terms of reaction net- works offer additional knowledge for the partic- ular biological process, possibly leading to further modifications and improvements of the particular model. An essential step towards the generalizion of purely chemical reaction network, involving just particular chemical substances towards a more generalized chemical objects, such as enzymes and substrates, seem to be done first by the prominent scientist Victor Henri, who proposed the enzyme kinetic reaction network, see example 7 below. Later on scientists working in fields, such as pop- ulation dynamics and epidemiological modelling, began to realize that many of their models can be based on reaction networks, wherein the chemical substances are considered as more generalized biological objects often denoted as species [6]. The logistic model is an instructive example of a growth-decay model possessing a chemical reaction network [6], [13], [21]. The reaction net- work presentation enables an easy identification of the logistic model as a constituent part of other (more complex) growth-decay models and to expect similar behaviour of the growth and decay functions, such as sigmoidal growth and exponential decay. Section 2 is intended for readers who are not familiar with the Chemical Reaction Network The- ory (CRNT) and its application to mathematical modelling in biology. On several examples we give a brief introduction of method of “translation” of a reaction network into a system of ordinary differ- ential equations (ODE’s) using the simple “mass action kinetic” principle. Such a translation turns the reaction network into an unique mathematical problem for the time evolution of the masses (concentrations , densities) of the species. Readers who are interested in the implementation of the more involved “power law kinetic” postulates may consult some textbooks on CRNT [8], [20]. In Section 3 we consider the classical Gompertz model from the perspective of CRNT. For this pur- pose we formulate the Gompertz model in terms of a reaction network. Readers already equipped with the technique of translation will be able to easily Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 2 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... trace the relation between the “chemical” reaction network and the classical ODE’s involved. This approach allows us to enlighten certain character- istics of the Gompertz growth function and the closely related Gompetz decay function involved (known as mortality law). In Section 4 we propose a new modification of the reaction network of the classic Gompertz model obtaining thus an original Gompertz-like model possessing one more degree of freedom. The proposed new model is mathematically anal- ysed in the spirit of the reaction network approach. It is shown that the model is a generalization of both the classic Gompertz model and the one-step radioactive exponential decay model, forming thus a hybrid between these two familiar models. II. PRELIMINARIES: REACTION NETWORKS AND THEIR TRANSLATION INTO ODE’S We briefly recall some features of growth-decay models based on reaction networks as well as some appropriate terminology and notation. A. Reaction networks. Canonical forms. Systems of chemical reac- tions, briefly: reaction networks, are symbolically presented as systems of elementary reactions of the following canonical form: S + Q k−→ P + R, (1) showing that one, two or more species on the left side of the reaction arrow, called reactants or reagents (in this example species S and Q), react, and, as result of the reaction, one, two or more species, named products (here P and R), are produced. Note that the arrow should point to the right and the sign “+” has different meaning when placed on the left or on the right side of the reaction arrow: on the left side the “+” means reaction between the enlisted reactants, whereas on the right no reaction is assumed between the product speciess; if such a reaction exists, then it should be described by a separate elementary reaction. Just one species on each side of the arrow is also possible, in fact, as we shall see below, a reaction of the form S k−→ P is a basic one. Note that a presentation, such as S k−→ P k1−→ Q, is not canonical, the corresponding canonical presentation for this reaction network is: S k−→ P, P k1−→ Q. As another example, the often used non-standard presentation of a reverse reaction network: S k−→←− k1 P has the following canonical form: S k−→ P, P k1−→ S. As one more example of a non-canonical famil- iar reaction network let us mention the enzyme kinetic reaction scheme between an enzyme E with a single active site and a substrate S, forming an enzyme-substrate complex C, which then yields product P , known as Henri-Michaelis-Menten re- action [7]: S + E k1−→←− k−1 C k2−→ P + E. In canonical form the above reaction network should be presented as: S + E k1−→ C, C k−1−→ S + E, C k2−→ P + E. As we shall see below canonical forms are useful for an easier transform of a reaction network into a system of differential equations via mass action principle. Notation. All species (reactants and products) partaking in a reaction are denoted by uppercase letters, such as P,Q,R,S,X,Y . A positive num- ber called “rate parameter” is written over the reaction arrow and indicates the velocity of the reaction. The reactants on the left side of a reaction arrow either decay or remain constant, whereas the product species on the right side of the arrow are growing or constant. In some cases a species may appear two or more times at one side of the arrow, such as A + A, briefly denoted as 2A. Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 3 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... For example, in the reaction (1) the reactants S,Q on the left side decay with time, whereas the product species P,R on the right side are growing. In a biological context we can say that in this case the product species grow at the expenses of the reactant species, or that the growing species consume the species on the left side of the arrow (as food resources). In population dynamics (some of) the species on the left side could be considered as “parent” species that give birth (reproduce) into the outcome species on the right side of the arrow. Catalyst species. One more special case should be mentioned, namely when a species X appears on both sides of the arrow, e.g. S +X k−→ P +X. In this case species X does not change in time, it is called a catalyst. Catalyst species enable a reaction to perform, e. g. in this example without X, the reaction S k−→ P cannot practically happen. When studying biological growth/decay pro- cesses, it is important to identify the species with catalytic action, the catalysts. In this work catalysts will be usually denoted by some of the letters X,Y,Z. Once again, by definition a catalyst species appears on both sides of a reaction arrow: on the left side as a reactant and on the right side as a product. Note also that some species may partake as catalysts in a particular reaction, but can also be involved as reagents in other reaction(s) as part of the same network; there they may change (grow or decay). In such situations a catalyst species may change as result of its total participation in a particular reaction network [20]. In the “logistic” reaction S + X k−→ 2X species X is a catalyst which catalyses the reaction S k−→ X, that is X catalyzes its own production (growing). Such species are called auto-catalysts. As a (total) result of the logistic reaction, the catalyst species X is growing in time. B. Differential systems induced by reaction net- works via mass action kinetics Law of Mass action: The rate of a reaction is proportional to the (mathematical) product of the concentrations of the reactants. Using mass action kinetics principle every par- ticular reaction network can be uniquely trans- formed (translated) into a system of ordinary dif- ferential equations (ODE’s), briefly: system of rate equations, or dynamical system [20], [23]. Such a transformation (“translation”) is performed in the following way. Firstly, we assume homogenous distribution of the species involved, say P,Q,S, in a certain volume/area/compartment. Then the quantitative (numeric) values assigned to species P,Q,S, such as masses, concentrations, densities, number of entities (individuals, molecules, cells, etc.), are considered as smooth functions of time denoted respectively: p = p(t),q = q(t),s = s(t), so that their derivatives, resp. p′ = dp(t)/dt,q′ = dq(t)/dt,s′ = ds(t)/dt, exist up to a certain order. Functions p,q,s are briefly called concentrations or masses, and the first derivatives of the masses p′,q′,s′ are called rates of change (growth or decay or both) of the respective species. Under these assumptions, the mass action principle says that the rate of change of each species is proportional to the product of the masses of the reacting species, thereby the coefficient of proportionality is negative if the species decays (which is the case when it appears on the left hand-side of the reaction arrow) and is positive if the species grows (when appearing on the right hand-side of the arrow). The proportionality coefficient is called the rate parameter and is usually written over the reaction arrow. This procedure is performed for every elementary reaction of the system of reactions, that is the reaction network. The “translation” of a reaction network into a system of ODE’s via mass action kinetics is illustrated on the following examples. Example 1. Consider the reaction network S + R k−→ P, k > 0. (2) Under mass action kinetics principle, reaction (2) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 4 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... induces a system of three differential equations, one for each of the three species, also known as reaction equations or rate equations: s′ = −ksr, r′ = −ksr, p′ = ksr, k > 0. (3) System (3) demonstrates the application of the mass action law when two reacting species S,R produce a new species P . In this reaction species S,R decay, whereas species P grows; so the signs of the rate parameters for S and R are negative, and the sign of the rate parameter for P is positive. The rates (of change) of all three species are proportional to the product of the masses of the reacting species S and R, so the absolute values of all rates are of the form ksr, with k > 0. Dynamical system (3) implies the identities s′+ r′ = 0, s′+p′ = 0, resp. s+r = c1 = const, s+p = c2 = const. Such identities are often known as conservation relations (laws). Example 2. Consider the reaction network: S + X k−→ P + X, k > 0. (4) The induced system of ODEs is: s′ = −ksx, x′ = 0, p′ = ksx, k > 0. (5) Note that equation x′ = 0 is obtained as x′ = −ksx + ksx = 0. Species X is a catalyst. The masses of species S and P satisfy the identity s+ p = const. Example 3. Consider the following reaction network involving two reactions and three species: S + X k−→ P + X, X α−→ P, (6) wherein k > 0,α > 0. The first reaction S + X k−→ P + X of reaction network (6) does not cause changes in catalyst species X, whereas the second reaction X α−→ P causes exponential decay of X. The other declining species is S; species P is growing. The induced dynamical system is: s′ = −ksx, x′ = −αx, p′ = ksx + αx. (7) Note that the rate p′ of species P is obtained as the sum ksx + αx of the rates of P from the two reactions involving P . Note also that the catalyst species X changes (decays) due to reaction x′ = −αx. Dynamical system (7) induces the identity p + x + s = const . Example 4. For the Henri-Michaelis-Menten reaction network S + E k1−→←− k−1 C k2−→ P + E a correct translation produces the following system of ODE’s for the concentrations s,e,c,p of the resp. species S,E,C,P : s′ = −k1es + k−1c, e′ = −k1es + (k−1 + k2)c, c′ = k1es− (k−1 + k2)c, p′ = k2c. In this last example the “most difficult” rates formulation seem to be the rates e′ = de/dt and c′ = dc/dt as they correspond to three distinct reaction arrows: e.g. for c′ one incoming in species C and two outgoing arrows from C. The example also demonstrates one more practically useful property of reaction networks: they are more obvious than the resp. systems of ODE’s and more easy to memorize. The above examples illustrate the process of translation of a chemical reaction networks into systems of ordinary differential equations. They also illustrate the derivation of an identity relation connecting the state variables in the system of ODE’s. We are going next to illustrate how the induced dynamical systems can be further mathe- matically analysed. C. Growth-decay models based on reaction net- works Below we present three case studies of famil- iar growth/decay functions generated by reaction networks. The goal is to demonstrates the use of the reaction networks methodology in several aspects: i) parallel to the growth function other useful functions appear (such as decay and wave- like functions) that should be analysed; ii) the original reaction networks offers meaningful inter- pretations of the resulting model solutions; iii) in Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 5 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... the process of modelling and numerical simulation of particular data sets the modeler can modify the basic reaction network by introducing meaningful changes in (some of) the elementary reactions. This latter possibility will be demonstrated in Section 4, where we modify the classic Gompertz model. Case study 1. Consider the reaction (network): S k−→ P, k > 0, (8) known in chemistry as a “first-order (FO)” reaction and in nuclear physics as “one-step exponential radioactive decay (1-SERD)”. This elementary re- action is known under several additional names due to its application to various processes such as radioactive nuclear decay, fluid dynamics, enzyme kinetics, marine ecology, physico-chemistry, etc. By definition, a first-order reaction proceeds at a rate that depends linearly on only one reactant concentration. Indeed, reaction (8) induces the following dynamical system for the change rates of the concentrations s = s(t), p = p(t) of species S,P : s′ = −ks, p′ = ks, k > 0. (9) Dynamical system (9) illustrates how the expres- sion “product of masses (concentrations)” should be interpreted in the definition of the mass action principle when just one species appears on the left hand-side of the reaction arrow. In such a situation the “product” consists of only one state variable, in the case of system (9)—concentration s. System (9) implies the relation s′ + p′ = 0, which after integration gives the identity (conser- vation) relation s + p = c = const. (10) Identity (10) says that at any time moment variable p gains as muuch as s loses. In certain real life situations this could be interpreted either as: i) species P consumes species S as a food resource, or: ii) compartment S migrates (flows, transforms) into compartment X. Thus reaction network (8) exhibits a specific mechanism for the time evolution of the two species S and P . Species P grows for the expense of species S, which proportionally decays. When equipped with initial conditions s(0) = s0 > 0, p(0) = p0 ≥ 0, (11) such that s0 + p0 = c, dynamical system (9) turns into an initial value problem (IVP) (9)–(11) and relation (10) becomes s + p = s0 + p0 = c, hence s = s0 + p0 −p = c−p. Let us briefly analyse the IVP for ODEs (9)– (11). Substituting s = c − p in equation p′ = ks we obtain an autonomous ordinary differential equation for the growth function p of the form: p′ = k(c−p), (12) with initial condition p(0) = p0. The differential equations (9), (12) for functions s and p under initial conditions (11) admit explicit solutions as functions of t ≥ 0. The solution for s is the familiar first-order exponential (radioactive) decay: s(t) = s0e −kt. (13) Function s has an asymptote s(∞) = 0, that is concentration s vanishes at infinity. Species with such a property are known as “limiting reagents” in chemistry. The solution for p can be obtained using identity (10) when substituting s by c−p in (13): p(t) = c−s = c− (c−p0)e−kt, c = s0 + p0. (14) Function p solving (14) has an upper asymptote p(∞) = c; it is known as exponential (also saturation) growth model. To obtain the absolute change rates of functions s,p we can differentiate expressions (13), (14), or insert the expression for s,p in the resp. differen- tial equations. For the absolute decay rate of s we obtain: s′ = −ks0e−kt. (15) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 6 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... For the absolute growth rate of function p we have p′ = (c−s)′ = ks0e−kt = k(c−p0)e−kt, c = s0 + p0. (16) The logarithmic (relative) decay rate rs of func- tion s is constant, namely: rs = (ln s) ′ = s′ s = −ks s = −k. (17) The logarithmic (relative) growth rate rp of p can be obtained as follows: rp = (ln p) ′ = p′ p = ks c−s = k ( c−p p ) = k ( (c−p0)e−kt c−(c−p0)e−kt ) . (18) The second derivatives of functions s,p are: s′′ = (−ks)′ = k2s > 0, (19) resp.: p′′ = k(c−p)′ = −kp′ = −k2(c−p0)e−kt < 0. (20) Relations s′′ > 0 and p′′ < 0 show that function s is convex, whereas function p is concave for all t ≥ 0. Remarks. 1. In numerical simulation studies the upper asymptote c of the growth function p, also known as environmental capacity, is usually set to one: p(∞) = c = 1. 2. Reaction (8) is the first chain-link of a multi-step exponential radioactive decay chain. In the field of radioactive decay and some population studies one is only interested in the decay process and ignores the evolution of the growth species. In such situations one often presents reaction (8) in the form S k−→ Ø, k > 0. The symbol Ø indicates that the reaction equation for the growth species in dynamical system (9), in our case equation p′ = ks, is suppressed. Case study 2. This example is an extension of the previous reaction network (8). An exponential mechanism involving two sequential first order steps in the transformation of three species S,P,Q is given in the reaction network: S k1−→ P, P k2−→ Q, (21) where k1,k2 are positive rate parameters. (As already mentioned, reaction network (21) is often written in the concise non-canonical form S k1−→ P k2−→ Q,). In nuclear physics reaction (21) is known as rwo-step exponential radioactive decay (2-SERD). Denoting the concentrations (densities, masses) of species S,P,Q as functions of time t by s = s(t),p = p(t),q = q(t) and their derivatives resp. by s′,p′,q′, we arrive at the following dynamical system: s′ = −k1s, p′ = k1s−k2p, q′ = k2p. (22) Dynamical system (22) induces the following conservation identity: s + p + q = c = const. (23) System (22) shows that s′ < 0 and q′ > 0, so function s decays, whereas function q grows. It can be proved that function p first increases until a certain time moment t∗ and then decreases in [t∗,∞). Such functions are called unimodal, their graphs are wave-like; such functions will also be considered as growth-decay functions. In chem- istry, species like P , having zero concentration at the beginning and at the end of the process, are called “intermediate”. A detailed discussion of reaction network (21) and the solutions s,p,q to system (22) are given in [5]. For the solution to general n-step exponential radioactive decay system of differential equations the reader may consult [4]. Case study 3. Let us discuss the familiar logis- tic model as induced by a reaction network. The logistic (Verhulst) growth function is originally introduced in [28] as the solution of a differential equation of the form x′ = kx(c−x). The solution of this equation is a sigmoidal growth function x = x(t), t ∈ R. One usually ignores the related decay function which is implicitly involved in the right-hand side of the differential equation as c−x. Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 7 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... In contrast, the growth-decay presentation of the logistic model based on reaction networks involves simultaneously the two functions—growing and decaying—as a single tuple (pair). The logistic growth-decay pair is generated by the following reaction network involving two reacting species S,X: S + X k−→ 2X, (24) wherein k is a positive rate parameter. As already mentioned, the symbol 2X in (24) is an abbrevi- ation for X + X. Reaction network (24) shows that S is a de- caying species, and X is a growing species that catalyses the 1-SERD reaction S k−→ X, hence, X is an auto-catalyst species. Under the assumption of mass action kinetics re- action network (24) induces the following dynam- ical system of two differential reaction equations for the masses (concentrations, densities) s = s(t), x = x(t) of species S,X, resp.: s′ = −ksx, x′ = kxs, k > 0. (25) Due to s′+x′ = 0, after integration, system (25) induces the conservation identity relation: s + x = const = c. (26) Assume initial value conditions s(0) = s0 > 0, x(0) = x0 > 0, (27) satisfying relation (26), so that s0 + x0 = c. (28) The initial value problem (25)–(27) implies the following autonomous differential equations for the growth function x and the decay function s: x′ = kx(c−x), x(0) = x0, s′ = −ks(c−s), s(0) = s0 = c−x0. (29) Differential equations (29) show that function x is monotonically increasing and bounded in R+ with values the interval [x0,c), where the value c = s0 + x0 is known as (environmental) carrying capacity. More specifically, function x approaches asymptotically c: x(∞) = x∞ = c. Function s is monotonically decreases approaching zero: s(∞) = s∞ = 0. As traditionally accepted in the literature, we shall assume c = 1, thus relation (26) becomes s + x = 1. (30) Equations (29) posses explicit algebraic solu- tions for t ∈ R. To find solution x we have to solve: x′ x(1 −x) = k, which can be written as x′ x + x′ 1 −x = k. (31) Integrating (31) we obtain ln x− ln (1 −x) = kt + ln C, or x 1 −x = Cekt, C = x0 1 −x0 , which can be presented as x = Cekt 1 + Cekt = x0 (1 −x0)e−kt + x0 . (32) For the boundary values of x at t = 0, t = ∞, expression (32) gives resp. x(0) = x0, x(∞) = 1, as expected. Using expression (32) for the growth function x, the decay function s is readily obtained from identity relation (30) as follows: s = 1 −x = (1−x0)e −kt (1−x0)e−kt+x0 = s0e −kt s0e−kt+(1−s0) = s0 s0+(1−s0)ekt . (33) Absolute and logarithmic (relative) change rates. To obtain the absolute rate of change of the growing species X, also called absolute growth rate (AGR), we can differentiate expression (32). Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 8 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... Alternatively, we can substitute the obtained expressions for s and x in the equation x′ = ksx from (29) to obtain: x′ = k s0e −kt s0e−kt+(1−s0) · x0 (1−x0)e−kt+x0 = kx0(1−x0)e−kt [(1−x0)e−kt+x0]2 . (34) For the boundary values of function x′ we have x′(0) = kx0s0 = kx0(1 −x0), x′(∞) = 0. The logarithmic change rate of function x, known also as relative growth rate (RGR), is defined as: rx = (lnx) ′ = d(lnx)/dt = x′ x . (35) The RGR (35) of x can be obtained by substi- tuting expression (33) for function s in differential equation x′/x = ks to get: rx = x′ x = ks = k(1 −x0)e−kt (1 −x0)e−kt + x0 . (36) For the boundary values of function rx = x′/x we have rx(0) = ks0 = k(1 −x0), rx(∞) = 0. To obtain the absolute change (decay) rate of the species S, we can proceed as follows: s′ = (1 −x)′ = −x′ = − kx0(1−x0)e −kt [(1−x0)e−kt+x0]2 = − ks0(1−s0)e −kt [s0e−kt+(1−s0)]2 . (37) The boundary values of function s′ are s′(0) = −ks0(1 −s0), s′(∞) = 0. The logarithmic change rate (relative decay rate) of species S is: rs = (lnx) ′ = d(lnx)/dt = s ′ s = −kx = −k x0 (1−x0)e−kt+x0 = − k(1−s0) s0e−kt+(1−s0) . (38) For the boundary values of the relative decay rate we have rs(0) = −ks0, rs(∞) = 0. Inflection point of the growth function. To look for inflection points of growth function x we need an expression for function x′′ = x′′(t): x′′ = (x′)′ = (kxs)′ = k(x′s + s′x) = k ((kxs)s + (−kxs)x) = k2xs(s−x). (39) Expression (39) reduces the solution of equation x′′(t∗) = 0 for an inflection point t∗ to equation s(t∗) = x(t∗), (40) showing that the values of the decay function s and the growth function x at t∗ are identical. The two equations (40): s(t∗) = x(t∗), and s(t∗) + x(t∗) = 1, due to (30), imply s(t∗) = x(t∗) = 1/2. Thus we have: x(t∗) = x0 (1 −x0)e−kt ∗ + x0 = 1 2 , (41) equivalently e−kt ∗ = x0 1 −x0 = x0 s0 , (42) or t∗ = − 1 k ln x0 1 −x0 = ln( x0 (1 −x0) )− 1 k . (43) Expression (43) for t∗ shows that for t∗ > 0, that is for the existence of inflection of the growth function, it is necessary that the logarithm in (43) is positive, that is x0 (1 −x0) = x0 s0 < 1, that is x0 < s0, hence x0 < 1/2. Consequently, when x0 ≥ 1/2 growth function x has no inflection. In this case we have x′′(t) < 0 for all t ≥ 0, hence growth function x is concave on R+. In the special case x0 = s0 = 1/2, we obtain the simple expressions: x = 1/(1 + e−kt), s = e−kt/(1 + e−kt) = 1/(1 + ekt). Lag time (phase). Let us find an expression for the slope of function x at the inflection point. Using (34) and (42), we obtain x′(t∗) = kx0(1−x0)e−kt ∗ [(1−x0)e−kt∗+x0]2 = kx0 2 (x0+x0) 2 = k 4 . (44) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 9 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... Using (44), we can compute the lag time y, which satisfies the relation x(t∗) y = x′(t∗) = k 4 , hence y = 2/k. Applications of the logistic and the one-step exponential radioactive decay (1-SERD) mod- els. The logistic model finds numerous applica- tions. A popular application is the Lotka-Volterra “prey-predator” model in population dynamics. Denoting the prey species by S and the predator by X, in its simplest form the Lotka-Volterra model can be written as a reaction network [6]: S + X k−→ 2X, S ν−→ 2S, X µ −→ Ø, (45) wherein k,ν,µ are positive rate parameters. Reac- tion network (45) induces the dynamical system: s′ = −ksx + νs, x′ = ksx−µx. (46) The logistic reaction S + X k−→ 2X describes the natural reproduction of the predator population. Reaction X µ −→ Ø represents the mortality of the predator. Reaction S ν−→ 2S describes the natural reproduction of the prey population. The last two reactions vary in different versions of the Lotka-Volterra (45), however, the logistic reaction remains usually the same. Another familiar application of the logistic model is the epidemiological SI model, where S stays for susceptible population and I for infective one: S + I k−→ 2I. (47) As we see, the basic epidemiological reaction network (47) coincides with the logistic reaction (24). Again, the epidemiological SI model (47) is the backbone of various modifications, such as the popular SIR model, where R means “removed” (or “recovered”) population: S + I k−→ 2I, I ν−→ R, (48) where k > 0, ν > 0 are positive rate parameters. As a further extension to (48), the “vital” SIR model includes additionally newborn (B) and dead (D) population compartments; in the simple case of equal birth and death rates the reaction-network- formulation of the vital SIR model reads: S + I k−→ 2I, I ν−→ R, D µ −→ B, (49) where k > 0, ν > 0, µ > 0 are positive rate parameters. The last reaction: D µ −→ B looks somewhat strange; however, it describes adequately the situation in stable populations. Models (45), (48), (49) demonstrate an useful property of the reaction-network-formulation of models. Namely, such a formulation allows a mod- eller to construct easily various combinations of existing familiar elementary models with already established characteristics. We shall demonstrate this property in Section 4 with a modification of the Gompertz model implementing in it certain features of the 1-SERD and the logistic models. The two models considered next in the present work—the classic and modified Gompertz models—can serve as further examples for our proposed methodology of treating growth-decay models induced by reaction networks. III. THE CLASSIC GOMPERTZ MODEL REVISITED FROM THE PERSPECTIVE OF REACTION NETWORKS THEORY The Gompertz growth function has been ini- tially designed for insurance purposes [9], and later used more generally as a modelling growth function in life sciences, like the logistic one [33]. It is usually presented as the explicit algebraic solution x = x(t) to an autonomous differential equation of the form: x′ = νx ln(1/x), ν > 0. In the sequel we deduce the Gompertz function x = x(t) starting from a reaction network us- ing the terminology of CRNT. This allows us to obtain and analyse the Gompertz growth function together with the related decay function, giving Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 10 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... us a general view on the Gompertzian growth- decay process, as well as a meaningful physico- chemical interpretation of the state variables and rate parameters involved. The reaction network approach explanes and unifies the two approaches. Consider the following reaction network involv- ing three species S,X,Q and two reactions: S ν−→ Q, S + X k−→ 2X + S, (50) wherein ν,k are positive rate parameters [5], [21]. Remarks. Reaction S + X k−→ 2X + S of network (50) says that both species X and S act as catalysts. More specifically species X catalyses the reaction S−→X + S, turning it into reaction: S + X−→X + X + S. So, X is a growing species autocatalysing itself. On the other hand, species S is also a catalyst in this reaction, it catalyses the reaction: X−→X + X. As result in this reaction, species S does not change in time; however, globally S changes (declines) as result of the first-order exponential decay reaction S−→Q. The latter reaction shows that S flows (migrates) into species Q, that is, outside the system of the two compartments of our interest (S,X). As mentioned, species Q can be replaced by the symbol Ø: S ν−→ Ø, meaning that we shall ignore the time evolution of species Q. Assuming homogeneity, denoting the mass- related numerical characteristics (such as concen- trations, masses, densities, etc.) of species S,X, resp., by lowercase letters s,x, reaction network (50) induces the following system of two ODE’s for the state variables s = s(t), x = x(t), t ∈ R+ = [0,∞) [21]: s′ = −νs, x′ = kxs, (51) where ν > 0, k > 0 are rate parameters. System (51) belongs to the class of biochemical systems (S-systems), cf. [25], [27] [26], [29], [31]. From system (51) we see that function s satisfies the uncoupled autonomous first order differential equation: s′ = −νs, ν > 0. As mentioned in Section 2, Case study 1., solution s = s(t) is given by: s(t) = s0e −νt, t ∈ R+, (52) for any initial value s(0) = s0 > 0. Hence, function s is monotone decreasing, exponentially approaching zero at t −→∞. Proposition 1. Let functions s = s(t), x = x(t) satisfy the system of ODE’s (51) for t ∈ R+. Then the following identity relation holds true in R+: γs + ln x = ln x(∞), γ = k/ν. (53) wherein x(∞) = x∞ = x(t)|t−→∞ is the ordinate of the horizontal asymptote of growth function x. Proof: Dynamical system (51) implies the identity: x′ = kxs = −kx(s′/ν), which can be written as: γs′ + x′/x = 0, γ = k/ν. (54) The integration of (54) yields γs + ln x = const = c. This equation shows that, while func- tion s decreases with time, function x increases, however, the latter increase is bounded by the constant c in the equation. The constant c has an important geometric meaning. Indeed, boundary values s(∞),x(∞) satisfy identity (53), so that γs(∞) + ln x(∞) = c, γ = k/ν. (55) Using that for t −→ ∞ function s = s0e−νt approaches zero for any positive s0, ν, symbol- ically s(∞) = s∞ = s|t−→∞ = 0, expression (52) implies ln x(∞) = c. This proves identity (53). The asymptote of the growth function. The constant c = ln x(∞) from identity (53) deter- mines the value of the horizontal asymptote x = x(∞) of growth function x(t). As traditionally done in the literature on Gompertz model, we fix the value for the asymptote as x(∞) = x|t−→∞ = 1. This choice of the asymptote leads to the value of c in expression (53) as c = ln x(∞) = ln 1 = 0. (56) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 11 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... Fixing c = 0, identity (53) becomes: γs + ln x = 0, γ = k/ν. (57) In what follows we shall consider the solution to (51) as an ordered 2-tuple (pair) (s,x), satisfying identity (57) in R+. Under the choice c = 0 relation (57) guarantees that the growth solution x approaches the asymptote x = 1 at t −→∞. Initial value problem. We shall next consider the system of ODE’s (51) as initial value problem involving an initial tuple (s(0),x(0)) = (s0,x0) for the solutions. Identity (57) is satisfied by solution (s,x) for all t ≥ 0 including t = 0 and t = ∞. Hence, when considering system (51) as an initial value problem, we shall naturally assume that the initial tuple (s0,x0) satisfies identity (57), i.e.: γs0 + ln x0 = 0, γ = k/ν. (58) Relation (58) restricts the range of x0 in the interval x0 ∈ (0, 1). Indeed, if x0 ≥ 1, then (58) implies s0 ≤ 0, which makes no practical sense. So, the choice c = 0 scales the total evolution of the (monotonically increasing) growth function x in the range x ∈ [x0, 1). In contrast, the monotonically decreasing decay function s ranges in the interval (0,s0], thereby the value s0 can be greater than one, s0 > 1. Identity (57) implies the following practically useful relations: ks + ν ln x = 0, (59) or, equivalently, using notation δ = 1/γ = ν/k: s = −δ ln x = ln x−δ, x = e−γs, (60) in particular. at t = 0; s0 = −δ ln x0 = ln x0−δ, x0 = e −γs0, (61) to be used in the calculations to follow. In particu- lar, for s0 = 1 we need to have, according to (61): x0 = e −γ. We now formulate the following: Proposition 2. Let initial value pair (s0,x0) sat- isfy 0 < x0 < 1, s0 = −δ ln x0, δ = ν/k, (62) then: 1) solution (s,x) to initial value problem (51)– (62) satisfies in R+ = [0,∞) relation (57): γs + ln x = 0; 2) solution x to system (51) satisfies the au- tonomous differential equation: x′ = νx(− ln x); (63) 3) solution x to initial value problem (51)–(62), can be presented in the form x = x0 e−νt. (64) Proof: 1) Using initial values (62) the integration of relation (54) under the choice c = 0 yields (57) together with x∞ = 1. 2) Substituting (59): ks = −ν ln x, in differen- tial equation x′ = kxs yields: x′ = kxs = x(ks) = x(−ν ln x) = νx(− ln x), = νx ln(1/x), (65) which is the familiar Gompertz differential equa- tion (63). 3) Using expressions (60), (52), solution x = x(t) can be obtained from relation (57), as follows: ln x = −γs = −γ(s0e−νt) = (−γs0)e−νt = (ln x0)e−νt = ln x0 e−νt, (66) resp. for x we obtain (64). This proves the propo- sition. Using part 3 of Proposition 2 we can present the solution tuple (s,x) to initial value problem (51)–(62) in the form (s,x) = ( s0e −νt, x0 e−νt ) . (67) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 12 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... Change rates. To obtain an explicit algebraic expression for the absolute growth rate of species X we write: x′ = kxs = (ks)x = (−ν ln x)x = −νln x0 e−νt exp(ln x0 e−νt), (68) which is positive due to ln x0 < 0. For the boundary values of function x′ we have: x′(0) = kx0s0 = x0(−ν ln x0) > 0, x′(∞) = kx(∞)s(∞) = 0. For the logarithmic (relative) growth rate rx = rx(t) of Gompertz growth function x we obtain: rx = (ln x) ′ = x′/x = −ν ln x0 e−νt = ln x0−νe −νt . (69) For the boundary values of rx = x′/x we have: rx(0) = −ν ln x0 e0 = ln(1/x0)ν, rx(∞) = −ν ln x∞ e−∞ = 0. (70) To obtain the absolute change (decay) rate of species S we write s′ = −νs = −νs0e−νt. (71) The boundary values of function s′ are s′(0) = −νs0, s′(∞) = 0. The logarithmic (relative) change rate rs = rs(t) of species S is constant (and so are the boundary values of rs): rs = s′ s = −ν. (72) Inflection points. Consider next the existence of a possible inflection point t∗ for the Gompertz growth function x. For the second derivative x′′ = x′′(t) of growth function x we have: x′′ = (x′)′ = (kxs)′ = k(x′s + s′x) = k ((kxs)s + (−νs)x) = kxs(ks−ν) = k2xs(s−ν/k). (73) Expression (73) reduces the solution of equation x′′(t∗) = 0 for t∗ to equation s(t∗) = ν/k = δ, (74) showing that the value of the decay function s at the inflection point t∗ is equal to the rate parameter ratio δ = ν/k. Using (52), equation (74) reads: s(t∗) = s0e −νt∗ = δ, hence e−νt ∗ = δ/s0 = 1/(γs0), (75) thus we obtain t∗ = (1/ν) ln(γs0) = ln(γs0) 1 ν . (76) Expressed via x0, the inflection time moment t∗ can be obtained when substituting γs0 in (76) by ln(1/x0): t∗ = ln(ln 1 x0 ) 1 ν . (77) To compute the value x(t∗) we can use relations (74): x(t∗) = e−γs(t ∗) = e−γδ = e−1 = 1/e. (78) Expression (76) implies: in order to have t∗ > 0, that is to exist an inflection point for growth function x on R+, it is necessary relation 1 < γs0 to take place. In terms of x0 this reads (using γs0 = − ln x0): 1 < ln(1/x0), equivalently: x0 < 1/e. (79) The condition for existence of inflection point s(t∗) = δ < s0 is equivalent to 0 < x0 < 1/e = x(t∗). Lag time (phase). To compute the so-called lag time interval of growth function x for the classic Gompertz model, we need the slope of function x at the inflection point, that is x′(t∗). Denote the intersection of the tangent line to the graph of x through the inflection point t∗,x(t∗) with the abscissa and the asymptote x = x∞ = 1, resp. by (ta, 0) and (tb, 1). The length of interval [ta, t∗] on the abscissa is the lag time, whereas the length of the interval [t∗, tb] is the log time. Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 13 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... Substituting x(t∗) = 1/e, resp. ln x′(t∗) = −1, in expression x′ = ksx = −νx ln x, we obtain x′(t∗) = ν/e. (80) The lag time L = t∗ − ta is equal to the ratio L = x(t∗)/x′(t∗) = 1/ν, (81) observing the triangle below its vertex point (t∗,x(t∗)). We summarize the results obtained on the ex- pressions for solution rates, inflection points and lag/log times in the following Proposition 3. 1. Solution tuple (s,x) to Gom- pertz initial value problem (51)–(62) is character- ized by the following properties: 1a. Solution (s,x) is given by (67): (s,x) = ( s0e −νt, x0 e−νt ) , thereby γs = − ln x, in particular: γs0 = − ln x0, γ = k/ν = 1/δ. The boundary values of Gompertz growth/decay functions s,x are: s(0) = ln x0 −δ, s(∞) = 0; x(0) = e−γs0 , x(∞) = 1. 1b. The absolute change rates of species S,X are given by expressions (71), (68): s′ = −νs0e−νt; x′ = −ν ln x0 e−νt x0e −νt . The boundary values of functions s′, x′, are: s′(0) = −νs0, s′(∞) = 0; x′(0) = −νx0 ln x0 > 0, x′(∞) = 0. 1c. The logarithmic change rates of functions s,x are given by (72), (69): rs = (ln s) ′ = s ′ s = −ν, rx = (ln x) ′ = x′/x = ln x0 −νe−νt. For the boundary values of the logarithmic change rates of Gompertz growth/decay functions s,x we have: rs(0) = rs(∞) = −ν; rx(0) = ln(1/x0) ν, rx(∞) == 0. 2a. The inflection point t∗ of Gompertz growth function x is given by (77): t∗ = ln ( ln 1 x0 ) 1 ν . The values of the growth/decay functions at inflection point t∗ are, cf. (74), (78): s(t∗) = δ = ν/k, x(t∗) = e−1. For the existence of inflection point in [0,∞), the necessary and sufficient condition is s0 > δ, resp.: 0 < x0 < 1/e. 2b. The lag time L is given by the ratio (81): L = x(t∗)/x′(t∗) = 1/ν. Remarks on the logistic and Gompertz mod- els. 1) The inflection point of the growing species X in the Gompertz model is lower than those in the logistic model, 1/e < 1/2. As a consequence, the Gompertzian growth curve has a shorter lag time, resp. longer lag (ageing, mortality) time, than the logistic growth curve. In both models the growing species X reproduces by a doubling mechanism, being constrained by species S which declines with time until vanishing. The inhibiting decay mechanism of species S is different in the two growth-decay models. In the logistic case species S is consumed by X as nutritional (food) resource (S charges X); thereby X is the solely species using S. In contrast, in the Gompertz model species S serves as a catalyst for X; thereby S charges some “other” species as well. The catalytic vs. the resource-charging role of species S turns out to be decisive in the distinction of the two models. 2) Both the logistic and the Gompertz models make use of just one rate parameter, which is not so obvious in the Gompertz model. The parameter k in the Gompertz model participates only in the identity relation and its role there is to determine the value of s0, resp. the limit value (one) of the upper asymptote of the growth function. Without loss of generality the parameter k can be set to one, se e.g. [32]. The decisive role of the rate parameter ν is noticed by many authors, Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 14 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... using for ν names such as “relative maturity rate”, “mortality rate”, etc. IV. A NEW MODEL BASED ON A MODIFIED (HYBRID) GOMPERTZ-LIKE REACTION NETWORK In this section we propose and mathematically analyse a growth-decay model induced by a reac- tion network that is close to Gompertz network (50) but borrows some features of the one-step exponential decay 1-SERD model. Consider the following reaction network involv- ing two species S and X: S ν−→ X, S + X k−→ 2X + S, (82) wherein ν,k are positive rate parameters. Denoting the mass-related quantitative (numer- ical) characteristics of species S,X, resp. by s,x, under the assumption of mass action kinetics, reac- tion network (82) induces the following dynamical system of two reaction equations: s′ = −νs, x′ = kxs + νs, (83) where ν,k are positive rate parameters. Proposition 4. If functions s = s(t), x = x(t) satisfy the system of ODE’s (83) on R+ = [0,∞), then the following identity relation holds true on R+: γs + ln(x + 1/γ) = ln(x∞ + 1/γ), (84) wherein γ = k/ν and x∞ = x(∞) = x(t)|t−→∞. Proof: System (83) implies the relation: s′ ν + x′ kx + ν = 0, or γs′ + x′ x + ν/k = 0. After integration, the above relation leads to the following identity γs + ln(x + 1/γ) = const = c, γ = k/ν. (85) As in the classic Gompertz model, solution s to system (83) satisfies the autonomous ordinary differential equation s′ = −νs, with solution (13) (or (52)). Hence function s is monotone decreas- ing, approaching zero: s(∞) = s∞ = 0. Passing to limit t −→ ∞ in identity (85), using s(∞) = 0, we obtain const = c = ln(x(∞) + 1/γ), hence (84). Identity (84) suggests that while function s monotonically decays, function x monotonically grows remaining bounded from above by x(∞), so the line x = x(∞) is a horizontal asymptote for the growth function x = x(t). We have the freedom to choose the boundary value x(∞) for x at t = ∞; so, as done tradition- ally, we set x∞ = x(∞) = 1. Using boundary values s∞ = 0, x∞ = 1, we obtain relation (84) in the form γs + ln(x + 1/γ) = ln(1 + 1/γ), or, using notation δ == 1/γ = ν/k: γs + ln(x + δ) = ln(1 + δ), equivalently γs + ln x + δ 1 + δ = 0. (86) Remark. Introducing the “deviated” growth function xδ: xδ = x + δ 1 + δ , in relation (86), we obtain γs + ln xδ = 0, which formally matches the corresponding identity (57) for the classic Gompertz model: γs + ln x = 0. This similarity takes place for a number of results to follow. In fact it is possible to rewrite most of the classical Gompertz results from Section 3 replacing function x by xδ, then performing a reverse transformation: x = xδ(1 + δ) −δ. From relation (86) we can obtain expressions for s in terms of x and for x in terms of s. Here are given some practically useful relations: s = δ ln 1 + δ x + δ = ln ( x + δ 1 + δ )−δ , (87) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 15 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... γs = ln ( 1 + δ x + δ ) , (88) eγs = 1 + δ x + δ , (89) x = (1 + δ)e−γs −δ. (90) In the special case t = 0 we have the following relations for the initial pair (s0,x0), assuring the limit condition x∞ = 1: s0 = ln ( 1 + δ x0 + δ )δ = ln ( x0 + δ 1 + δ )−δ ; (91) γs0 = ln ( 1 + δ x0 + δ ) ; (92) eγs0 = 1 + δ x0 + δ ; (93) x0 = (1 + δ)e −γs0 −δ. (94) Substituting s from (87) in the differential equa- tion for growth function x in dynamical system (83), leads to the following autonomous differen- tial equation: x′ = kxs + νs = k(x + δ)s = k(x + δ) ln ( 1+δ x+δ )δ . (95) To deduce an explicit solution for growth func- tion x, we first use relation (89) to write: 1 + δ x + δ = eγs = eγs0e −νt = (eγs0 ) e−νt . (96) We then substitute the term eγs0 in (96), using the expression (93), to get 1 + δ x + δ = (eγs0 ) e−νt = ( 1 + δ x0 + δ )e−νt . (97) Relation (97) implies an explicit expression for growth function x = x(t): x(t) = (1 + δ) ( x0 + δ 1 + δ )e−νt −δ. (98) Based on the above considerations, we formu- late the following Proposition 5. Let initial value tuple (s0,x0) be such that 0 < x0 < 1, s0 = ln ( 1+δ x0+δ )δ > 0, δ = 1/γ = ν/k, (99) then i) solution (s,x) to initial value problem (83)– (99) satisfy on R+ = [0,∞) relation (84); in particular relations (87), (90): s = ln ( x + δ 1 + δ )−δ ; x = (1 + δ)e−γs −δ. ii) the growth function x satisfies the au- tonomous ordinary differential equation (95); x′ = k(x + δ) ln ( 1 + δ x + δ )δ ; iii) the solution x to equation (95) with initial value x(0) = x0, resp. system (83)–(99) can be presented in the explicit form (98): x(t) = (1 + δ) ( x0 + δ 1 + δ )e−νt −δ. Change rates. To obtain an explicit algebraic expression for the absolute growth rate of species X we use expressions (95) and (98) to obtain: x′ = k(x + δ) ln ( 1+δ x+δ )δ = k(1 + δ) ( x0+δ 1+δ )e−νt ln ( 1+δ x0+δ )e−νt . (100) For the boundary values of function x′ we have: x′(0) = kx0s0 = x0(−ν ln x0) > 0 x′(∞) = kx(∞)s(∞) = 0. (101) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 16 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... For the relative (logarithmic) growth rate rx = rx(t) of growth function x we obtain: rx = (ln x) ′ = x′/x = −ν ln x0 e−νt = ln x0−νe −νt . (102) For the boundary values of rx = x′/x we have: rx(0) = −ν ln x0 e0 = ln(1/x0)ν, rx(∞) = −ν ln x∞ e−∞ = 0. (103) The expressions for the absolute and logarithmic change (decay) rates of decay species S are the same as those for the classic Gompertz model, cf. (71), (72). Inflection points. To calculate the inflection points of the growth function x (if any) we need to obtain an expression for the second derivative x′′ of x: x′′ = (x′) ′ = (ksx + νs)′ = (ksx)′ + (νs)′ = k(s′x + sx′) + νs′ = k(−νsx + s(ksx + νs)) −ν2s = ks(−νx + ksx + νs−ν2/k) = ks (ks(x + ν/k) −ν(x + ν/k)) = ks(x + δ)(ks−ν) = k2s(x + δ)(s−δ). (104) According to expression (104) equation x′′(t) = 0 is equivalent to equation s(t)−δ = 0, or s(t) = δ. Let time instant t∗ solve the latter equations, then s0 > s(t ∗) = δ (105) is a necessary condition for the existence of an inflection point. Indeed, if (105): s0 > δ, then the monotone decreasing function s(t) equals to δ at time instant t∗: s(t∗) = δ. In other words, for s0 > δ then there exists time moment t∗, such that the pair t∗,x(t∗) is an inflection point for growth function x, such that s(t∗) = δ, resp. x′′(t∗) = 0. Expression (104) reduces the solution of equa- tion x′′(t∗) = 0 for t∗ to equation s(t∗) = ν/k = δ, (106) saying that the value of the decay function s at inflection time instant t∗ is equal to the rate parameter ratio δ = ν/k. Using (52), we have s(t∗) = s0e−νt ∗ = δ, hence e−νt ∗ = δ/s0 = 1/(γs0), (107) thus we obtain t∗ = (1/ν) ln(γs0) = ln(γs0) 1 ν . (108) Let us now “translate” formulae (105), (108) in terms of growth function x. Expressed via x0, the inflection time instant (108) can be obtained when substituting γs0 in (108) by ln((1 + δ)/(x0 + δ)): t∗ = ln(ln 1 + δ x0 + δ ) 1 ν . (109) Knowing the s-value s(t∗) = δ, we compute the x-value x(t∗) using expression (90): x(t∗) = (1 + δ)e−γs(t ∗) −δ = (1 + δ)e−γδ −δ = (1 + δ)e−1 −δ = (1 − (e− 1)δ)/e, thus finally we have: x(t∗) = 1 − (e− 1)δ e . (110) From (110) we obtain a necessary and sufficient condition for the existence of inflection: 1 e > x(t∗) = 1 − (e− 1)δ e > x0 > 0. (111) Relation (111) implies a necessary condition for the existence of inflection: x(t∗) > 0, namely: δ < 1 e− 1 = e ≈ 0.58197671, (112) resp. ν < e k. (113) Using (111) we obtain a second necessary condition for the existence of inflection: 1/e > x(t∗) > x0, namely: 1 − (e− 1)δ > ex0, (114) Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 17 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... resp. δ < 1 −ex0 e− 1 = e(1 −ex0), (115) resp. ν < e(1 −ex0)k. (116) Practically, the necessary and sufficient condi- tion (111) for the existence of inflection point can be tested by verifying the two necessary conditions (112), (115), resp. (113), (116). Inequality (112) implies the following restric- tion on the initial value x0 for existence of inflec- tion point: x0 < 1/e. (117) Remarks. i) Restriction (117) says that for the existence of inflection initial value x0 should be below the inflection value for the classical Gompertz model, i.e. 1/e. ii) Note that in the classic Gompertz case the growth function may have no inflection only when x0 > 1/e. In contrast, the hybrid Gompertz growth function may have no inflection for any initial values x0 ∈ (0, 1), even for initial values satisfying (117). iii) Depending on the values of the rate param- eters ν,k, the inflection point can be arbitrarily close to the classic Gompertz value 1/e, as well as to initial value x0 no matter how small x0 is. In the latter case the inflection point can be arbitrarily close to zero (providing x0 itself is sufficiently small). This possibility makes the shape of the graph of x extremely flexible, which makes a considerable difference with the classic Gompertz case. Under suitable choice of the initial conditions and rate parameters the hybrid model can be close to the one-step exponential decal model. iv) Inequality (112) implies δ < 1 −ex0 e− 1 < 1 e− 1 . (118) The “rough” inequality (118) can be used when x0 is close to 0. Lag time (lag phase). To compute the lag time of growth function x for the hybrid Gompertz model, we need the value of slope of function x at inflection time moment t∗, that is x′(t∗). Denote the intersection of the tangent line through the inflection point (t∗,x(t∗)) with the abscissa and the asymptote x = x∞, resp. by (ta, 0) and (tb, 1). The width (length) of interval [ta, t∗] is by definition the lag time. To compute the slope x′(x∗) of growth function x at inflection time moment t∗, we substitute the value: x∗ = x(t∗) from (110), resp. x∗ + δ = (1 + δ)/e, in the expression for the slope x′ to obtain: x′(t∗) = k(x∗ + δ) ln ( 1+δ x∗+δ )δ = kδ e (1 −δ) ln ( e1+δ 1+δ ) = ν e (1 −δ). (119) As in the classical Gompertz case, we define the lag time (phase) L as the length of the segment on the abscissa between inflection moment t∗ and the intersection point of the abscissa and the tangent with slope x′(t∗). Hence, for the lag time L we obtain: L = x∗ x′(t∗) = 1 ν ( 1 − eδ 1 + δ ) . (120) We summarize the obtained results as follows. Proposition 6. 1. Solution pair (s,x) to initial value hybrid Gompertz problem (83), (s(0) = s0,x(0) = x0), is characterized by the following properties: 1a. The absolute change rate of species X is given by: (100): x′ = k(1 + δ) ( x0 + δ 1 + δ )e−νt ln ( 1 + δ x0 + δ )e−νt . 1b. For the boundary values of function x′ we have expression (101): x′(0) = kx0s0 = x0(−ν ln x0) > 0 x′(∞) = kx(∞)s(∞) = 0. 1c. The logarithmic change rate rx = rx(t) of the hybrid Gompertz growth function x is given by expression (102): rx = (ln x) ′ = x′/x = −ν ln x0 e−νt. Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 18 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... For the boundary values of the logarithmic change rates of the hybrid Gompertz growth func- tion x we have (103): rx(0) = −ν ln x0 e0 = ln(1/x0)ν, rx(∞) = −ν ln x∞ e−∞ = 0. (121) 2a. The inflection time moment t∗ of the hybrid Gompertz growth function x is given by (109): t∗ = ln(ln 1 + δ x0 + δ ) 1 ν . The values of the growth/decay functions at inflection point t∗ are (74): s(t∗) = δ = ν/k, resp. (110): x(t∗) = 1 − (e− 1)δ e . The slope of the tangent line through the inflec- tion point is given by (119): x′(t∗) = ν e (1 −δ). For the existence of inflection point in [0,∞), the necessary and sufficient condition is (111): 1 e > x(t∗) = 1 − (e− 1)δ e > x0 > 0. 2b. The lag time L is given by the ratio (120): L = x∗ x′(t∗) = 1 ν ( 1 − eδ 1 + δ ) . Finally, the following proposition holds true: Proposition 7. The hybrid Gompertz function (98) is a generalization of the classical Gompertz func- tion (64). Proof: The classical Gompertz function (64) is obtained from the hybrid Gompertz function (98) for the special case k −→ ∞, resp. δ = ν/k −→ 0, while keeping the rate parameter ν fixed. V. CONCLUDING REMARKS Biological growth functions are usually pre- sented in the mathematical literature by means of their explicit expressions or as solutions to differential equations [11]–[18]. However, biolog- ical growth models are usually related to de- cay processes/functions, which becomes especially transparent when the models are based on reaction equations. Using chemical reaction network the- ory, one can easily observe close relations between various growth/decay processes, as well as be- tween existing growth-decay models, e.g. classes of biochemical systems [25]–[27], [29]–[31]. In the present work we propose an elemen- tary introduction in the reaction network approach based on mass action kinetics. To this end we dis- cuss in some detail several familiar examples, such as the one- and two-step exponential (radioactive) decay, the logistic and the Gompertz models. We focus on the simultaneous analysis of the growth and the decay functions using the iden- tity relation between the two functions naturally induced by the reaction equations. The power of the reaction network approach is fully revealed in Section 3 when applied to the analysis of the classical Gompertz model. There we propose a revision of the model based on the reaction network inducing the original Gompertz model, which we call “Gompertz reaction net- work” in honour of the author of the well-known growth model and his seminal paper [9]. Our final goal in this work is the modification of the Gompertz reaction network in a natural way, using fully the dynamical features of the one-step (first-order) exponential decay reaction. In this way we obtain a hybrid of the one-step exponential and the classic Gompertz model in a natural way, performing a small modification in the Gompertz reaction network. The growth function of the obtained new hybrid Gompertz- like model possesses one additional degree of free- dom (one more rate parameter) and is thus more flexible when applied to modelling and numeri- cal simulation of measurement and experimental Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 19 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... data sets. More specifically, the ordinate (height) of the inflection point of the hybrid model can largely vary, whereas the resp. height of the classic Gompertz model is fixed (at 1/e). The presented generalization of the Gompertz model possesses some common features with the Richards model in direction of improved flexibility when simulating measurement data sets [24], [34]. Acknowledgments. This paper is supported by the National Scientific Program “Information and Communication Technologies for a Single Dig- ital Market in Science, Education and Security (ICTinSES)”, contract No. DO1–205/23.11.2018, financed by the Ministry of Education and Science in Bulgaria. The author is indebted to his col- leagues the professors R. Anguelov, M. Krastanov, N. Kyurkchiev, for their constant encouragements and useful suggestions. Special thanks also go to Dr. A. Vassilev for his expert numerical simula- tions (to be published jointly in a forthcoming article). REFERENCES [1] Anguelov, R., Borisov M., Iliev A., Kyurkchiev N., S. Markov, On the chemical meaning of some growth mod- els possessing Gompertzian-type property. Math. Meth. Appl. Sci. 2017, 1-2, https://doi.org/10.1002/mma.4539 [2] Asadi, M., Di Crescenzo, A., Sajadi, F.A. et al. A generalized Gompertz growth model with applications and related birth-death processes. Ricerche mat (2020). https://doi.org/10.1007/s11587-020-00548-y [3] Bajzer, Z., Vuk-Pavlovic, S., New Dimensions in Gom- pertzian Growth. Journal of Theoretical Medicine 1997; 2:307–315. [4] Bateman, H., The solution of a system of differential equations occurring in the theory of radio-active trans- formations, Proc. Cambridge Phil. Soc. 15 (1910), 423– 427. [5] Borisov, M., S. Markov, The two-step exponential de- cay reaction network: analysis of the solutions and relation to epidemiological SIR models with logistic and Gompertz type infection contact patterns, Journal of Mathematical Chemistry, 59(5), 1283–1315, DOI: 10.1007/s10910-021-01240-8 [6] Chellaboina, V., Bhat, S. P., Haddat W. M., D. S. Bern- stein, Modeling and Analysis of Mass-Action Kinetics. IEEE Control Systems Magazine 2009; 60–78. [7] Deichmann, U., S. Schuster, J.-P. Mazat, A. Cornish- Bowden, Commemorating the 1913 Michaelis–Menten paper Die Kinetik der Invertinwirkung: three perspec- tives, FEBS J. 281 (2014) 435—463. J.-P. Mazat, Part 3: 452—463. [8] Feinberg, M., Foundations of Chemical Reaction Network Theory, Applied Mathematical Sciences 202, Springer, 2019. https://doi.org/10.1007/978-3-030- 03858-8 [9] 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. [10] Hethcote, H. W., The Mathematics of Infectuous Dis- eases, SIAM REVIEW, 42, 4, 599–653 (2000). [11] Iliev, A., N. Kyurkchiev, S. Markov, A Note on the New Activation Function of Gompertz Type, Biomath Communications, 4, No. 2 (2017), 51–64. [12] Iliev A, Kyurkchiev N, S. Markov, On the approxi- mation of the cut and step functions by logistic and Gompertz functions. Biomath 2015; 4: 2–13. [13] Kyurkchiev, N., S. Markov, On the Hausdorff distance between the Heaviside step function and Verhulst logis- tic function. J. Math. Chem. 2016; 54:109–119. [14] Kyurkchiev N., V. Kyurkchiev, A. Iliev, A. Rahnev, A new modification of the SIR/SEIR models with “inter- vention polynomial factor”. Methodological aspects, Int. J. of Differential Equations and Applications 20, No. 1 (2021), 15–30. [15] Kyurkchiev, N., The new trasmuted C.D.F. based on Gompertz function, Biomath Communications 5, 2018. [16] Kyurkchiev, N., Some new classes of growth functions generated by reaction networks and based on “correct- ing amendments of Bateman-Gompertz and Bateman- Gompertz-Makeham type. I., Communications in Ap- plied Analysis, 24, No. 1 (2020), 13–29. [17] Laird, A K, Tyler, S A, Barton, A D., Dynamics of normal growth. Growth 1965; 29:233–248. [18] Laird, A K., Dynamics of growth in tumors in normal organisms. National Cancer Institute Monograph 1969; 30:15–28. [19] Lazarova, M., Markov, S., A. Vassilev, On Some Classes of Growth Functions and Their Links to Reaction Net- work Theory, AMiTaNS’20, AIP CP 2302 (American Institute of Physics, Melville, NY), 2020, paper 080004 [20] Lente, G., Deterministic Kinetics in Chemistry and Systems Biology. Briefs in Molecular Science, Springer 2016. [21] Markov, S., Reaction networks reveal new links between Gompertz and Verhulst growth functions, Biomath 8 (2019), 1904167. [22] Markov, S., On a class of generalized Gompertz– Bateman growth–decay models, Biomath Communica- tions 6 (1) (2019) 52–59. [23] Murray, J D., Mathematical Biology: I. An Introduction, Third Edition. Springer 2002. [24] Richards, F. J. (1959) A flexible growth function for empirical use. J Exp Bot 10, 290–300. [25] Savageau, M A. Allometric morphogenesis of complex Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 20 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 S M Markov, The Gompertz model revisited and modified using reaction networks: Mathematical ... systems: A derivation of the basic equations from first principles. Proceedings of the National Academy of Sciences 1979; 76:5413–5417. [26] Savageau M A, Voit E O. Recasting nonlinear differen- tial equations as S-systems: A canonical nonlinear form. Math. Biosci. 1987; 87:83–115. [27] Savageau M.A. (1990) Biochemical Systems Theory: Alternative Views of Metabolic Control. In: Cornish- Bowden A., Cardenas M.L. (eds) Control of Metabolic Processes. NATO ASI Series (Series A: Life Sciences), vol. 190. Springer, Boston, MA. [28] Verhulst, P.-F., Notice sur la loi que la population poursuit dans son accroissement. Correspondance math- ematique et physique 1838; 10:113–121. [29] Voit E. O., Biochemical Systems Theory: A Review, International Scholarly Research Notices, vol. 2013, Article ID 897658, 53 pages, 2013. [30] Voit E.O., Biochemical Systems Theory (BST). In: Dubitzky W., Wolkenhauer O., Cho KH., Yokota H. (eds) Encyclopedia of Systems Biology. Springer, 2013. [31] Voit, E.O. The best models of metabolism. Wiley Interdiscip Rev Syst Biol Med. 2017;9(6):10.1002/wsbm.1391. [32] West, J. et al., An Evolutionary Model of Tumor Cell Kinetics and the Emergence of Molecular Heterogene- ity Driving Gompertzian Growth, SIAM review 58 4 (2016), 716–736. [33] Winsor, C., Gompertz curve as a growth equation. Proc. Nat. Acad. Sciences, 18, 1–8, 1932. [34] Zwietering, M. H., Jongenburger, I., Rombout, F. M., K. van’t Riet, Modeling of the Bacterial Growth Curve, Applied and Environmental Microbiology 56 (6): 1875– 1881, 1990. doi:10.1128/AEM.56.6.1875-1881.1990 Biomath 10 (2021), 2110023, http://dx.doi.org/10.11145/j.biomath.2021.10.023 Page 21 of 21 http://dx.doi.org/10.11145/j.biomath.2021.10.023 Introduction: reaction networks and evolutionary growth-decay models Preliminaries: reaction networks and their translation into ODE's Reaction networks. Differential systems induced by reaction networks via mass action kinetics Growth-decay models based on reaction networks The classic Gompertz model revisited from the perspective of reaction networks theory A new model based on a modified (hybrid) Gompertz-like reaction network Concluding remarks References