IEEE/CAA Journal of Automatica Sinica  2018, Vol. 5 Issue(6): 1142-1149   PDF    
Mathematical Study of A Memory Induced Biochemical System
Mithun Kumar Ghosh1, Tridip Sardar2, Xianbing Cao3, Priti Kumar Roy1     
1. Centre for Mathematical Biology and Ecology, Department of Mathematics, Jadavpur University, Kolkata 700032, India;
2. Department of Mathematics, Dinabondhu Andrews College, Kolkata 700084, India;
3. College of Science, Beijing Technology and Business University, Beijing 100037, China
Abstract: In this work, to study the effect of memory on a bi-substrate enzyme kinetic reaction, we have introduced an approach to fractionalize the system, considering it as a threecompartmental model. Solutions of the fractionalized system are compared with the corresponding integer-order model. The equilibrium points of the fractionalized system are derived analytically. Their stability properties are discussed from numerical aspect. We determine the changes of the substances due to the changes of "memory effect". The effect is discussed critically from the perspective of product formation. We have also analyzed the memory induced system with a control measure in view of optimizing the product. Our numerical result reveals that the solutions of the fractionalized system, when it is free from memory, are in good agreement with the integer-order system. It is noticed that the effect of memory influences the reaction in the forward direction and assists in yielding the product more quickly. However, an extensive use of memory makes the system slower, but introduction of a control input makes the reaction faster. It is possible to overcome the slowness of the reaction due to the undue effect of memory by appropriate use of a control measure.
Key words: Bi-substrate enzymatic reaction     compartmental system     control theoretic approach     fractional-order differential equations     memory effect    

Enzymes are biological catalysts that are necessary in almost every biochemical reaction [1]. These enzymes are proteins synthesized by genes [2]. The main function of an enzyme is to catalyze the making and breaking of chemical bonds depending on an accurate sequence of amino acids and its complicated tertiary structure. The catalytic ability of enzymes increases the rate of a reaction. The enzyme is not used up in the reactions and, it does not change the equilibria of the processes [3]. This raises a new dimension of thinking towards various fields viz. physics [4], chemistry [5], biology [6], ecology [7], epidemiology [8], pharmacokinetics [9] etc. A lot of research has been done about enzymatic processes of different chemical and biochemical transformations. Enzyme kinetics is the study of rates of these reactions to optimize the velocity of reactions, rate of intermediate complexes and products.

For a better understanding of the reaction kinetics, many authors have implemented different techniques to obtain approximate analytical solutions of the enzymatic systems [10]-[13]. Modern day literature related to enzyme activity in enzymatic processes consist of mathematical approaches to study system dynamics for optimization and quantification of product [14]. Single substrate or double substrate biochemical reactions make the approaches more interesting, of which the latter is more reasonable and important [15]-[17].

Westerlund stated in [18] that every matter has memory. Although it is debatable, a large number of theoretical physicists considered the memory function as an embedded characteristic of molecular properties, which is discussed in various domains of science and engineering branches [19]-[21]. Toledo-Hernandez et al. mentioned in [22] that biochemical reactions involve the participation of living organisms viz. enzymes. The dynamic behavior of living microorganisms not only depends on their current state conditions (e.g., substrate concentration, medium condition, etc.), but also on their previous states. They have explained this phenomena as the dynamics of the reactions that involve memory effects. Now, it is to be noted that integer-order (IO) derivatives consider only local properties (at time $t$) while fractional derivatives take into account the history of a process i.e., their previous states [21]. An enzymatic reaction system with IO derivatives is in general memory-less [20], [23] and hence it is unable to reflect the effect of memory. The memory effect can be incorporated in a system by introducing fractional-order ($\alpha\in(0, 1]$) derivatives as an index of memory [24] i.e., $\alpha\rightarrow0$ indicates that the system has an ideal memory and $\alpha\rightarrow1$ represents that the system is free from memory.

The conception of fractional calculus is first projected by Leibniz [25] in $1695$. A fractional-order differential equation is considered as an alternative model to special nonlinear differential equations [26], [27]. In enzyme kinetics, Abdullah [28] employed FDEs in $2011$ for modeling the Michaelis-Menten reaction in a 2-d region containing obstacles. In $2013$, Alawneh [29] used the multistep generalized differential transform method to solve a time-fractional enzyme kinetics. They investigated dynamical behavior of various complex materials and systems for the benefits of more degrees of freedom and introduction of memory in the model. The drawback of both the studies was the way of introduction of memory. Both of them fractionalized the systems only by changing the order of the ordinary derivatives on the left-hand side of the ODEs. However, fractionalization of a system of two or more ODEs is not possible without violating mass balance and the system may suffer from unit inconsistencies [21], [30]. The problem is not limited to the units. The above difficulty can be removed by considering a common order of all the FDEs of the system, but it is a very particular case and makes its application restricted. Here, we present a more accurate model of a bi-substrate enzymatic reaction, where dynamics are influenced by memory.

In this article, we have introduced fractionalization of a two-substrate enzymatic reaction to study the effect of memory on it. Nonlinear FDEs cannot, in general, be solved analytically [27], but can be solved by numerical techniques [31]. The numerical solutions of the system have been studied here and compared with the integer-order system. We have also observed the dynamics of the different substances of the system by varying the order of the fractional derivatives (which signifies a measure of memory effect in a system [24]). We formulate a control based mathematical model involving the memory effect to conquer the negative effect of the extensive use of memory.

We have organized the rest of the paper as follows. In Section Ⅱ, we formulate the model of a bi-substrate enzymatic reaction involving the memory effect. Some basic theoretical properties, and the existence and stability of equilibrium points are studied in Section Ⅲ. In Section Ⅳ, a control theoretic approach is introduced towards the fractional-order model. The numerical results are illustrated in Section Ⅴ. Finally, we have completed our article with a discussion and conclusion of the study in Section Ⅵ.


The schematic diagram of a two-substrate enzyme kinetic reaction, as described by Roy et al. [15], is given by

$ \begin{eqnarray}\label{model0} E+S_{1} {\large\overset{k_{1}}{\underset{k_{-1}}{\rightleftharpoons}}} ES_{1}+S_{2} {\large\overset{k_{2}}{\underset{k_{-2}}{\rightleftharpoons}}} ES_{1}S_{2} \xrightarrow{k_{3}} E+P \end{eqnarray} $ (1)

where $S_1$, $S_2$ are substrates, $E$ is the enzyme, $C_1$ i.e., $ES_1$ and $C_2$ i.e., $ES_1S_2$ are intermediate complexes and $P$ is the product. $k_1$, $k_2$ are the rate constants of formation of the complexes $C_1$ and $C_2$ respectively, and $k_3$ is the rate of product formation. The rates of dissociation of $C_1$ and $C_2$ are $k_{-1}$ and $k_{-2}$ respectively.

Let us denote the concentrations $[S_1]$, $[S_2]$, $[E]$, $[C_1]$, $[C_2]$ and $[P]$ by $s_1$, $s_2$, $e$, $c_1$, $c_2$ and $p$ respectively. From the law of mass action [10], the above enzymatic reaction $(1)$ can be described by the following set of differential equations:

$ \frac{ds_1}{dt} = -k_1es_1+k_{-1}c_1\nonumber\\ \frac{ds_2}{dt}=-k_2c_1s_2+k_{-2}c_2\nonumber\\ \frac {de}{dt} = \ -k_1es_1+k_{-1}c_1+k_3c_2 \\ \frac {dc_1}{dt}=k_1es_1-k_{-1}c_1-k_2c_1s_2+k_{-2}c_2\nonumber\\ \frac {dc_2}{dt}=k_2c_1s_2-k_{-2}c_2-k_3c_2\nonumber\\ \frac {dp}{dt}= \ k_3c_2 $ (2)

with the initial conditions

$ \begin{eqnarray}\label{inicon} s_1(0) = s_{10}, \; s_2(0) = s_{20}, \; e(0) = e_0\nonumber\\ c_1(0) = 0, \; c_2(0) = 0,\; \mbox{and} \; p(0) = 0. \end{eqnarray} $ (3)

From system $(2)$, we have the following relations:

$ \begin{eqnarray}\label{reln} \frac{ds_1}{dt}-\frac{ds_2}{dt}+\frac{dc_1}{dt}\!\!&\;=\;&\!\!0\nonumber\\ \frac{ds_2}{dt}+\frac{dc_2}{dt}+\frac{dp}{dt}\!\!&\;=\;&\!\!0\nonumber\\ \frac{ds_2}{dt}+\frac{dc_2}{dt}+\frac{de}{dt}-\frac{ds_1}{dt}\;\; &\;=\;&\;\; 0. \end{eqnarray} $ (4)

Using the initial conditions $(3)$, from $(4)$, we have

$ \begin{eqnarray}\label{cons} s_1-s_2+c_1\!\!&\;=\;&\!\! s_{10}-s_{20}\nonumber\\ s_2+c_2+p\!\!&\;=\;&\!\!s_{20}\nonumber\\ s_2+c_2+e-s_1\!\!&\;=\;&\!\!s_{20}+e_0-s_{10}. \end{eqnarray} $ (5)

With the help of the relations $(5)$, system $(2)$ can be reduced to the following three dimensional model consisting of substrates $s_1$, $s_2$ and complex $c_2$ as

$ \begin{eqnarray}\label{model1.1} \frac{ds_1}{dt}\!\!&\;=\;&\!\!-k_1(s_1-s_2-c_2+s_{20}-s_{10}+e_0)s_1\nonumber\\ &+k_{-1}(s_2-s_1+s_{10}-s_{20}), \nonumber\\ \frac{ds_2}{dt}\!\!&\;=\;&\!\!-k_2(s_2-s_1+s_{10}-s_{20})s_2+k_{-2}c_2, \nonumber\\ \frac {dc_2}{dt}\;\; &\;=\;&\;\; k_2(s_2-s_1+s_{10}-s_{20})s_2-k_{-2}c_2-k_3c_2 \end{eqnarray} $ (6)

with initial conditions,

$ \begin{eqnarray}\label{ini1} s_1(0)=s_{10}, \;\;\; s_2(0)=s_{20}, \;\;\;c_2(0)=0. \end{eqnarray} $ (7)
A. Fractionalization of the System (6)

The schematic diagram $(1)$ can be considered as a three-compartmental model (as shown in Fig. 1) [30]. The initial stage of the reaction, where substrate $S_1$ is reacting with enzyme $E$ to form the complex $ES_1$, is termed as Compartment $1$. Compartment $2$ describes the intermediate stage where substrate $S_2$ is combining with $ES_1$ to form $ES_1S_2$ complex. Compartment $3$ consists of the yielding of $ES_1S_2$, which may either convert to the product or decompose back to the previous stage of reaction. Here, the mass flux $k_{-2}c_2$ is transferred from Compartment $3$ to $2$ (Fig. 1) and is common between the second and third equations of $(6)$. We can fractionalize the system $(6)$ as given below [30], [31]:

$ \begin{eqnarray}\label{modelrl} \frac{ds_1}{dt}\!\!&\;=\;&\!\!-k_1(s_1-s_2-c_2+s_{20}-s_{10}+e_0)s_1\nonumber\\ &+k_{-1}(s_2-s_1+s_{10}-s_{20}), \nonumber\\ \frac{ds_2}{dt}\!\!&\;=\;&\!\!-k_2(s_2-s_1+s_{10}-s_{20})s_2+{k^{\alpha}_{-2}}^{R}_0D^{1-\alpha}_t(c_2), \nonumber\\ \frac {dc_2}{dt}\!\!&\;=\;&\!\!k_2(s_2-s_1+s_{10}-s_{20})s_2\nonumber\\ &-{k^{\alpha}_{-2}}^{R}_0D^{1-\alpha}_t(c_2)-k_3c_2 \end{eqnarray} $ (8)
Fig. 1 Three-compartment model corresponding to the schematic diagram $(1)$ of a bi-substrate enzymatic reaction. Rectangular boxes containing $E+S_1$, $ES_1+S_2$ and $ES_1S_2$ represent Compartments $1$, $2$ and $3$ respectively. $E(0)$, $S_1(0)$ and $S_2(0)$ are the respective initial values of $E$, $S_1$ and $S_2$. $k_1$, $k_{-1}$, $k_2$ and $k_{-2}$ are the rate at which the mass fluxes are transferred from the source compartment to the targeted one as directed, where $k_3$ is the rate of elimination of product $P$ and enzyme $E$ from Compartment $3$.

where $0 <\alpha \leq 1$ and $^{R}_0D^{\alpha}_t f(t)$ represents the Riemann-Lioville (RL) fractional derivative. In this context, it is to be noted that the memory index $\alpha$ was assigned between $0$ to $1$ in several classical systems and various engineering disciplines [27], [30], [32]-[34].

The unit on left-hand side of all the sub-equations of system $(8)$ is $\mbox{Hour}^{-1}$. Mol/l, the unit of the substances, is nothing but a number. The unit on right-hand side of the first sub-equation of $(8)$ is also $\mbox{Hour}^{-1}$. Now, the unit of $k_{-2}^{\alpha}$ is $\mbox{Hour}^{-\alpha}$ and that of $^{R}_0D^{1-\alpha}_t(c_2)$ i.e., of $\frac{d^{1-\alpha}c_2(t)}{dt^{1-\alpha}}$ is $\mbox{Hour}^{\alpha-1}$. Thus, the unit of the term ${k_{-2}^{\alpha}}^{R}_0D^{1-\alpha}_t(c_2)$ and consequently, the right-hand side of both the second and third sub-equations of $(8)$ is $\mbox{Hour}^{-1}$. Hence, units of all of the sub-equations of $(8)$ remain consistent under the fractionalization process we have considered.

In order to use standard initial conditions, the Riemann-Lioville derivatives must be re-defined as Caputo fractional derivatives [22]. The relation between RL and Caputo's derivatives is given by the following equation:

$ \begin{eqnarray}\label{rlcaputorlsn} ^{R}_0D^{1-\alpha}_t f(t)\!\!&\;=\;&\!\!^{C}_0D^{\alpha}_t f(t)+\frac{f(0)t^{\alpha-1}}{\Gamma(\alpha)}. \end{eqnarray} $ (9)

The system $(8)$ using Caputo derivative can be expressed as follows:

$ \begin{eqnarray}\label{modelcp} \frac{ds_1}{dt}\!\!&\;=\;&\!\!-k_1(s_1-s_2-c_2+s_{20}-s_{10}+e_0)s_1\nonumber\\ &+k_{-1}(s_2-s_1+s_{10}-s_{20}), \nonumber\\ \frac{ds_2}{dt}\!\!&\;=\;&\!\!-k_2(s_2-s_1+s_{10}-s_{20})s_2\nonumber\\ &+{k^{\alpha}_{-2}}^{C}_0D^{1-\alpha}_t(c_2)+\frac{A_1 t^{\alpha-1}}{\Gamma(\alpha)}, \nonumber\\ \frac {dc_2}{dt}\!\!&\;=\;&\!\!k_2(s_2-s_1+s_{10}-s_{20})s_2\nonumber\\ &-{k^{\alpha}_{-2}}^{C}_0D^{1-\alpha}_t(c_2)-k_3c_2-\frac{A_1 t^{\alpha-1}}{\Gamma(\alpha)} \end{eqnarray} $ (10)


$ \begin{eqnarray}\label{a1a2} A_1 = k^{\alpha}_{-2} c_2(0)=0\nonumber \end{eqnarray} $

Therefore, system $(10)$ becomes

$ \begin{eqnarray}\label{modelcpa1a20} \frac{ds_1}{dt}\!\!&\;=\;&\!\!-k_1(s_1-s_2-c_2+s_{20}-s_{10}+e_0)s_1\nonumber\\ &+k_{-1}(s_2-s_1+s_{10}-s_{20}), \nonumber\\ \frac{ds_2}{dt}\!\!&\;=\;&\!\!-k_2(s_2-s_1+s_{10}-s_{20})s_2+{k^{\alpha}_{-2}}^{C}_0D^{1-\alpha}_t(c_2), \nonumber\\ \frac {dc_2}{dt}\!\!&\;=\;&\!\!k_2(s_2-s_1+s_{10}-s_{20})s_2\nonumber\\ &-{k^{\alpha}_{-2}}^{C}_0D^{1-\alpha}_t(c_2)-k_3c_2. \end{eqnarray} $ (11)

In this section, we have determined the equilibrium points of model $(11)$ and discussed their stability from numerical point of view.

A. Existence of Equilibria and Stability

It is not possible to understand the stability of the equilibrium points of the system $(11)$ directly because the fractional derivative does not satisfy Leibniz rule [35]. We apply the following transformation:

$ \begin{eqnarray}\label{trns} X\!\!&\;=\;&\!\!D^{\alpha} s_1, \nonumber\\ Y\!\!&\;=\;&\!\!D^{\alpha} s_2-k^{\alpha}_{-2} c_2, \nonumber\\ Z\!\!&\;=\;&\!\!D^{\alpha} c_2+k^{\alpha}_{-2} c_2. \end{eqnarray} $ (12)

Using the above transformation $(12)$, system $(11)$ is thus equivalent to the following system:

$ \begin{eqnarray}\label{trnseq} D^{1-\alpha} X\!\!&\;=\;&\!\!-k_1(s_1-s_2-c_2+s_{20}-s_{10}+e_0)s_1\nonumber\\ &+k_{-1}(s_2-s_1+s_{10}-s_{20}), \nonumber\\ D^{1-\alpha} Y\!\!&\;=\;&\!\!-k_2(s_2-s_1+s_{10}-s_{20})s_2, \nonumber\\ D^{1-\alpha} Z\!\!&\;=\;&\!\!k_2(s_2-s_1+s_{10}-s_{20})s_2-k_3 c_2, \nonumber\\ D^{\alpha} s_1\!\!&\;=\;&\!\!X, \nonumber\\ D^{\alpha} s_2\!\!&\;=\;&\!\!Y+k^{\alpha}_{-2} c_2, \nonumber\\ D^{\alpha} c_2\!\!&\;=\;&\!\!Z-k^{\alpha}_{-2} c_2. \end{eqnarray} $ (13)

It is sufficient to study the stability properties of system $(13)$.

System (13) has the equilibrium points $E^{*}_{1}(0, 0, 0, 0, s_{20}-s_{10}, 0)$ for $\delta=s_{20}-s_{10}>0$ and $E^{*}_{2}(0, 0, 0, s^{*}_1, 0, 0)$ for $\delta\leq0$ where $s^{*}_1$ is given by the following equation:

$ \begin{eqnarray}\label{s1qdeq} k_1 s^{*^{2}}_1+[k_1(s_{20}-s_{10}+e_0)+k_{-1}]s^{*}_1\nonumber\\ +k_{-1}(s_{20}-s_{10})=0, \end{eqnarray} $ (14)

i.e., $s^{*}_1= \frac{-A+\sqrt{A^{2}+B}}{2k_1}, $ where

$ \begin{eqnarray}\label{s1ab} A\!\!&\;=\;&\!\!k_1(s_{20}-s_{10}+e_0)+k_{-1}, \nonumber\\ B\!\!&\;=\;&\!\!4 k_1 k_{-1} (s_{10}-s_{20}). \nonumber \end{eqnarray} $

The Jacobian matrix $J_{E^{*}_{1}}$ of system $(13)$ at $E^{*}_{1}$ is given by

$ \begin{eqnarray}\label{jacobe1} \left( \begin{array}{cccccc} 0&0&0&-(k_{-1}+k_1 e_0)&k_{-1}&0 \\ 0&0&0&k_2\delta&-k_2\delta&0 \\ 0&0&0&-k_2\delta&k_2\delta&-k_3 \\ 1&0&0&0&0&0 \\ 0&1&0&0&0&k^{\alpha}_{-2} \\ 0&0&1&0&0&-k^{\alpha}_{-2} \\ \end{array} \right) \end{eqnarray} $ (15)

where $ \delta = s_{20}-s_{10}.$

Since $\alpha$ is real, it can also be an irrational number. However, there is no existing method for studying such a system with an irrational order of fractional derivatives. Therefore, we assume that $\alpha=\frac{M}{N}$ is rational, where $N>M>0$ and $\mbox{gcd}(M, N)=1.$

Therefore, the characteristic equation of the matrix $J_{E^{*}_{1}}$ is given by,

$ \begin{eqnarray}\label{jacobe1charlame1} \triangle(J_{E^{*}_{1}}-\mbox{diag}([\lambda^{N-M}, \lambda^{N-M}, \lambda^{N-M}, \nonumber\\ \lambda^{M}, \lambda^{M}, \lambda^{M}]))=0, \end{eqnarray} $ (16)

where "$\triangle$" and "diag" represent the determinant and the diagonal matrix respectively [36].

Expanding the characteristic equation $(16)$, we have

$ \begin{eqnarray}\label{jacobe1chareqe1} \lambda^{3N}+\Lambda_{11} \lambda^{2N}+\Lambda_{12} \lambda^{N}+\Lambda_{13} \lambda^{N(3-\alpha)}\nonumber\\ +\Lambda_{14} \lambda^{N(2-\alpha)} + \Lambda_{15}=0 \end{eqnarray} $ (17)


$ \begin{eqnarray}\label{jacobe1charcoeffe1} \Lambda_{11}\!\!&\;=\;&\!\!k_3+k_2\delta+k_{-1}+k_1e_0, \nonumber\\ \Lambda_{12}\!\!&\;=\;&\!\!k_2k_3\delta+k_3(k_{-1}+k_1e_0)+k_1k_2e_0\delta, \nonumber\\ \Lambda_{13}\!\!&\;=\;&\!\!k^{\alpha}_{-2}, \Lambda_{14} = (k_{-1}+k_1e_0)k^{\alpha}_{-2}\; \mbox{and}\nonumber\\ \Lambda_{15}\!\!&\;=\;&\!\!k_1k_2k_3e_0\delta. \nonumber \end{eqnarray} $

For $\alpha = 1$, from $(17)$, we have

$ \begin{eqnarray}\label{intordrchareqne1} \lambda^{3}+\eta_1\lambda^{2}+\eta_2\lambda+\eta_3=0 \end{eqnarray} $ (18)


$ \begin{eqnarray}\label{coeffalphaeq1e1} \eta_1\!\!&\;=\;&\!\!k_{-2}+k_{-1}+k_3+k_1e_0+k_2\delta, \nonumber\\ \eta_2\!\!&\;=\;&\!\!(k_3+k_{-2})(k_{-1}+k_1e_0)+k_2\delta(k_3+k_1e_0), \nonumber\\ \eta_3\!\!&\;=\;&\!\!k_1k_2k_3e_0\delta.\nonumber \end{eqnarray} $

Equation $(18)$ is same as the characteristic equation of the integer-order system $(6)$ for the equilibrium point $(0, \delta, 0)$.

The Jacobian matrix $J_{E^{*}_{2}}$ of system $(13)$ at $E^{*}_{2}$ is given by,

$ \begin{eqnarray}\label{jacobe2} \left( \begin{array}{cccccc} 0&0&0&-\sqrt{A^{2}+B}&k_1s^*_1+k_{-1}&k_1s^*_1 \\ 0&0&0&0&k_2(\delta+s^*_1)&0 \\ 0&0&0&0&-k_2(\delta+s^*_1)&-k_3 \\ 1&0&0&0&0&0 \\ 0&1&0&0&0&k^{\alpha}_{-2} \\ 0&0&1&0&0&-k^{\alpha}_{-2} \\ \end{array} \right).\nonumber \end{eqnarray} $

Proceeding as above, we have the characteristic equation of $J_{E^{*}_{2}}$ as follows:

$ \begin{eqnarray}\label{jacobe1chareqe2} (\lambda^{N}+\sqrt{A^{2}+B})(\lambda^{2N}+\Lambda_{21} \lambda^{N}\nonumber\\ +\Lambda_{22} \lambda^{N(2-\alpha)}+\Lambda_{23})=0 \end{eqnarray} $ (19)


$ \begin{eqnarray}\label{jacobe1charcoeffe2} \Lambda_{21}\!\!&\;=\;&\!\!k_3-k_2(\delta+s^{*}_1), \nonumber\\ \Lambda_{22}\!\!&\;=\;&\!\!k^{\alpha}_{-2}, \nonumber\\ \Lambda_{23}\!\!&\;=\;&\!\!-k_2k_3(\delta+s^{*}_1). \nonumber \end{eqnarray} $

Arguments of the roots of the first factor of $(19)$ are of the form $\mbox{arg}(\lambda_k)=\frac{\pi}{N}+\frac{2k\pi}{N}, k=0, 1, 2, ..., N-1$ and hence $|\mbox{arg}(\lambda_k)|>\frac{\pi}{2N}$ for $k=0, 1, 2, ..., N-1$.

Therefore, stability of $E^{*}_{2}$ of system $(13)$ depends upon the nature of the roots of the remaining factor of $(19)$ as given by [36],

$ \begin{eqnarray}\label{jacobe1dtrmnstceqe2} \lambda^{2N}+\Lambda_{21} \lambda^{N}+\Lambda_{22} \lambda^{N(2-\alpha)}+\Lambda_{23}=0. \end{eqnarray} $ (20)

For $\alpha = 1$, we have, from $(19)$

$ \begin{eqnarray}\label{intordrchareqne2} (\lambda+\sqrt{A^{2}+B})[\lambda^{2}+[k_{-2}+k_3-k_2(\delta+s^{*}_1)]\lambda\nonumber\\ -k_2k_3(\delta+s^{*}_1)]=0 \end{eqnarray} $ (21)

which is same as the characteristic equation of the integer-order system $(6)$ for the equilibrium point $(s^{*}_1, 0, 0)$. Here $s^{*}_1$ is defined exactly as in $(14)$.


To determine the effect of memory towards the system, we introduce control parameter $u(t)$ into the model $(11)$. Our aim is to get an optimum amount of the product as quick as possible. The control input $u(t)$ is used to reduce the rate of reverse reaction at the second stage satisfying $0 \leq u(t) \leq 1$. $u(t)=1$ and $0$ represent maximum and minimum use of the control measure respectively. With these assumptions, model $(11)$ becomes:

$ \begin{eqnarray}\label{modelcpa1a20con} \frac{ds_1}{dt}\!\!&\;=\;&\!\!-k_1(s_1-s_2-c_2+s_{20}-s_{10}+e_0)s_1\nonumber\\ &+k_{-1}(s_2-s_1+s_{10}-s_{20}), \nonumber\\ \frac{ds_2}{dt}\!\!&\;=\;&\!\!-k_2(s_2-s_1+s_{10}-s_{20})s_2\nonumber\\ &+(1-u(t)){k^{\alpha}_{-2}}^{C}_0D^{1-\alpha}_t(c_2), \nonumber\\ \frac {dc_2}{dt}\!\!&\;=\;&\!\!k_2(s_2-s_1+s_{10}-s_{20})s_2\nonumber\\ &-(1-u(t)){k^{\alpha}_{-2}}^{C}_0D^{1-\alpha}_t(c_2)-k_3c_2 \end{eqnarray} $ (22)

where $s_1(0)=s_{10}$, $s_2(0)=s_{20}$ and $c_2(0)=0$.

Here the control measure basically stands for temperature, pressure, concentrations of the substances [7], [14] etc.. We study the effect of the control input on the system $(22)$ from numerical point of view.


In this section, dynamics of reaction kinetics have been analyzed with the help of numerical methods. There are various methods to solve a system of fractional-order differential equations. We have used the numerical scheme given in [31] and solved our system of equations using the Matlab subroutine "lsqnonlin" and called this method as NS-lsq. Here, we have observed the solutions of the fractional-order system for $\alpha=1$ i.e., when the system is free from the "memory effect" and compared them with the integer-order system. The stability region of equilibrium points of the system have also been studied numerically. We have compared the concentration of the substances, particularly the product, for various values of $\alpha$. Consequently, we determine how the rate of formation of the substances is influenced by the "memory effect" of the system. Here the parameter values are taken from [37]-[39]. The units and recommended values of the kinetic parameters used for numerical simulation are as given in Table Ⅰ.

Table 1
A. Comparison of the Substance Profiles Obtained From the Fractional System for $\alpha=1$ and From the Integer-order System

Fig. 2 represents the behavioral pattern of the substances for the integer-order (IO) model $(6)$ and the fractional-order (FO) model $(13)$ for $\alpha=1$ simultaneously. The solutions obtained from both the systems are in good agreement with each other. Concentration of the two substrates ($s_1$ and $s_2$) decreases with the progression of the reaction. Consumption of $s_1$ is relatively quicker than $s_2$ due to faster reaction between enzyme and the primary substrate. Initially, the concentration of the enzyme decreases due to the formation of enzyme-substrate complexes $c_1$ and $c_2$. It is recovered as the reaction progresses. Concentration of $c_1$ increases gradually from its initial value and it decreases with time, as it binds with $s_2$ while forming the second complex. Moreover, the concentration of $c_2$ increases as soon as the first complex is formed, and then is decreased as time progresses due to its transformation to the product. Fig. 2 displays continuous formation of the product with time until it becomes steady.

Fig. 2 Concentration profiles of the substances for integer-order system $(6)$ (dashed line) and the fractional-order system $(13)$ for $\alpha=1$ (circle). The parameter values are $k_1=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=1\mbox{h}^{-1}$, $k_2=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-2}=1\mbox{h}^{-1}$, $k_3=5\mbox{h}^{-1}$ and the initial values are $s_{10}=5\mbox{M}$, $s_{20}=5\mbox{M}$, $e_0=4.5\mbox{M}$ where $\mbox{M}$ stands for $\mbox{mol/l}$.
B. Stability Region of the Equilibrium Points

The stability regions of the equilibrium points of model system $(13)$ are studied numerically with the help of Theorem1 of Sardar et al. [31] for reasonable values of the model parameters, where $\alpha\in(0, 1)$. We have observed that both equilibrium points are stable for such parameter values, and differentiating only in the time allows it to reach the steady state. Fig. 3 represents two stability regions corresponding to the equilibrium points $E^{*}_{1}$ and $E^{*}_{2}$ of the system $(13)$.

Fig. 3 Stability regions of the equilibrium points $E^{*}_{1}$ (left panel) and $E^{*}_{2}$ (right panel) with respect to $\alpha$ and $k_{-2}$. Other parameter values corresponding to $E^{*}_{1}$ are taken as $k_1=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=1\mbox{h}^{-1}$, $k_2=4\mbox{M}^{-1}\mbox{h}^{-1}$, $k_3=5\mbox{h}^{-1}$, $e_0=2\mbox{M}$, $s_{10}=5\mbox{M}$, $s_{20}=6\mbox{M}$ and for $E^{*}_{2}$ as $k_1=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=1\mbox{h}^{-1}$, $k_2=4\mbox{M}^{-1}\mbox{h}^{-1}$, $k_3=5\mbox{h}^{-1}$, $e_0=1\mbox{M}$, $s_{10}=6\mbox{M}$, $s_{20}=5\mbox{M}$.
C. Comparison Among the Concentration Profiles for Different Values of $\alpha$

We have compared the dynamic profiles of the substances obtained by decreasing $\alpha$ values gradually to $1$ and $0.7$. It is to be noted that the solutions of $(13)$ for $\alpha=1$ correspond to the ODE system $(6)$. Fig. 4 displays the dynamic profiles of the substances obtained from the IO system $(\alpha=1)$ and the FO system $(\alpha=0.7)$.

Fig. 4 Concentration profiles of the substances for the fractional-order system $(13)$ for $\alpha=1$ (solid line) and $0.7$ (dashed line) where the parameter values are $k_1=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=3\mbox{h}^{-1}$, $k_2=6\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-2}=3\mbox{h}^{-1}$, $k_3=5\mbox{h}^{-1}$ and the initial values are $s_{10}=5\mbox{M}$, $s_{20}=5\mbox{M}$, $e_0=4.5\mbox{M}$.

As the value of $\alpha$ reduces to $0.7$, the concentrations of both the substrates ($s_1$ and $s_2$) decrease gradually. It is observed that with and without memory operator has no significant changes in the first substrate. This may be due to the fact that we do not consider memory in the first backward reaction step (see Fig. 1). The consumption of the second substrate is faster in comparison to the integer-order system. The profiles of the enzyme concentration show a faster recovery for the lower value of $\alpha$. Variations in concentration of the first complex ($c_1$) is observed under varying $\alpha$ values. It is found that, the complex concentration $c_1$ is lower for $\alpha=0.7$ than for $\alpha=1$. It implies that for the lesser value of $\alpha$, there exists a lower accumulation of $c_1$ due to its quicker conversion to the second complex by binding with the second substrate. A relatively faster accumulation of the concentration of second complex ($c_2$) is observed for the smaller value of $\alpha$. This indicates the possibility of higher conversion of it to the product. Yielding of product is relatively faster for a lower value of $\alpha$ and consequently, its concentration reaches the steady state more quickly. The effect of memory assists in comparatively quicker accumulation of $c_2$ and consequently, a rapid formation of the product.

To study the effect of changes of the time taken for formation of the product due to the changes in $\alpha$ values, we decrease gradually the values of $\alpha$ as $1$, $0.7$ and $0.25$. Fig. 5 represents the dynamic profiles of the product $(p)$ for the aforesaid $\alpha$ values. Concentrations of the product are observed as 4.771 mol/l, $4.948$ mol/l and 4.927 mol/l for $\alpha=1$, $0.7$ and $0.25$ respectively. It is to be noted that as the $\alpha$ values decrease from $1$ to $0.7$, the formation of product becomes relatively faster than the classic case. However, it is also to be noted that, if the values of $\alpha$ decreased again, the time taken for the formation of product is relatively greater (Fig. 5) and consequently, the system slows down.

Fig. 5 Concentration profiles of the product $p$ of system $(13)$ for $\alpha=1$, $0.7$ and $0.25$. Other parameter values are taken as $k_1=7\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=0.1\mbox{h}^{-1}$, $k_2=13\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-2}=3\mbox{h}^{-1}$, $k_3=12\mbox{h}^{-1}$, $e_0=4.5\mbox{M}$, $s_{10}=5\mbox{M}$, $s_{20}=5\mbox{M}$.
D. Comparison Among the Concentration Profiles for Different Values of $u(t)$

Here, we investigate the effect due to the changes of control parameter $u(t)$ to the control induced FO model $(22)$. Fig. 6 represents the variation in the substances for $u(t)=0$, $0.4$ and $0.9$. It is observed that the concentrations of both the substrates ($s_1$ and $s_2$) are decreasing more quickly for higher values of the control parameter. Accumulation of the first complex concentration $(c_1)$ is lower for upper values of $u(t)$ due to its fast conversion into the second complex. Higher values of the control input corresponds to a more accumulation of $c_2$ which ultimately leads to fast formation of the product $(p)$.

Fig. 6 Concentration profiles of the substances of $(22)$ for $u(t)=0$, $0.4$ and $0.9$. Other parameter values are taken as $k_1=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=2\mbox{h}^{-1}$, $k_2=5\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-2}=2\mbox{h}^{-1}$, $k_3=5\mbox{h}^{-1}$, $e_0=4.5\mbox{M}$, $s_{10}=5\mbox{M}$, $s_{20}=5\mbox{M}$ and $\alpha=0.7$.

We vary the values of the control parameter as $0$ and $0.6$ to rise above the negative effect of extensive use of memory. Fig. 7 represents the concentration profiles of the product for the above values of $u(t)$. The rest of the parameter values are taken exactly as in Fig. 5 with $\alpha=0.25$. Concentration of the product $(p)$ for $u(t)=0$ is observed as $4.927$ mol/l, which is decreasing from the value $4.948$ mol/l for $\alpha=0.7$ (see Fig. 5). While, in Fig. 7, concentration of the product for $u(t)=0.6$ is noted as 4.973 mol/l. Thus, with proper control measures, it is possible to overcome the above mentioned negative effect.

Fig. 7 Concentration profiles of the product $p$ of system $(22)$ for $u(t)=0$ and $0.6$. Other parameter values are taken as $k_1=7\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-1}=0.1\mbox{h}^{-1}$, $k_2=13\mbox{M}^{-1}\mbox{h}^{-1}$, $k_{-2}=3\mbox{h}^{-1}$, $k_3=12\mbox{h}^{-1}$, $e_0=4.5\mbox{M}$, $s_{10}=5\mbox{M}$, $s_{20}=5\mbox{M}$ and $\alpha=0.25$.

In this study, we have presented an approach of fractionalizing a bi-substrate enzyme kinetic reaction. The fractional-order system is solved numerically, as the system is unlikely to have analytical solutions. Our numerical results reveal that the solutions of the fractional-order system for $\alpha=1$, and the solutions of the corresponding integer-order system are overlapping. Benefits of the fractional-order model are observed from the solutions, mainly in the formation of the product.

In Section Ⅱ, we have calculated the equilibrium points of the FO system and discussed their stability regions. Our study shows that, similar to the integer-order system, equilibrium points remain stable for a fractional-order system with a realistic range of parameters.

We have studied the changes of the concentration profiles of the substances due to the changes in $\alpha$ values $(\alpha=1 \; \mbox{and} \; 0.7)$. Lower values of $\alpha$ signifies a faster reaction up to certain threshold values.

We have focused on the changes of concentration of the product due to the change in $\alpha$ values $(1, 0.7 \; \mbox{and} \;0.25)$. The system is highly sensitive to the $\alpha$ values. Formation of the product is perceived relatively faster due to "memory effect." However, extensive effect of memory makes the system slower. The results of our study can predict system dynamics with respect to optimization and quantification of the product.

The dynamical behavior of the substances is observed by varying the control input. Presence of the control parameter corresponds to a quicker reaction. The negative effect of the extensive use of memory can be recovered by proper use of a control measure.

The model can be extended by considering memory in both the backward reaction steps. One can consider the mass $k_{-1}c_1$, which is transferred from Compartment $2$ to $1$ similar to the way mass $k_{-2}c_2$ is transferred from Compartment $3$ to $2$ (see Fig. 1). In this context, with the help of relations $(5)$, system $(2)$ can be transformed to a three dimensional model consisting of the substrate $s_1$ and the complexes $c_1$, $c_2$. Proceeding as in Section Ⅱ, one can fractionalize the model. The fractionalized model would consist of terms ${k^{\alpha}_{-1}}^{C}_0D^{1-\alpha}_t(c_1)$ and ${k^{\beta}_{-2}}^{C}_0D^{1-\beta}_t(c_2)$ of different orders without violating mass balance. Hence in this study, we may summarily conclude that the presence of the mixing parameter may show complex dynamics.

[1] R. Dutta, Fundamentals of Biochemical Engineering. Berlin: Springer, 2008.
[2] I. Belgacem and J. L. Gouzé, "Global Stability of full open reversible michaelis-menten reactions, " IFAC Proc., vol. 45, no. 15, pp. 591-596, 2012.
[3] D. L. Nelson and M. M. Cox, Lehninger Principles of Biochemistry. 6th ed. Basingstoke: Macmillan Education, 2013.
[4] S. Sirin, D. A. Pearlman, and W. Sherman, "Physics-based enzyme design: predicting binding affinity and catalytic activity, " Proteins: Struct. Funct. Bioinform., vol. 82, pp. 3397-3409, Dec. 2014.
[5] D. Vasic-Racki, U. Kragl, and A. Liese, "Benefits of enzyme kinetics modelling, " Chem. Biochem. Eng. Quart., vol. 17, no. 1, pp. 7-18, Mar. 2003.
[6] R. Roskoski Jr, "The ErbB/HER family of protein-tyrosine kinases and cancer, " Pharmacol. Res., vol. 79, pp. 34-74, Jan. 2014.
[7] P. K. Roy, S. Datta, S. Nandi, and F. Al Basir, "Effect of mass transfer kinetics for maximum production of biodiesel from Jatropha Curcas oil: a mathematical approach, " Fuel, vol. 134, pp. 39-44, Oct. 2014.
[8] S. D. Thiberville, N. Moyen, L. Dupuis-Maguiraga, A. Nougairede, E. A. Gould, P. Roques, and X. De Lamballerie, "Chikungunya fever: epidemiology, clinical syndrome, pathogenesis and therapy, " Antiv. Res., vol. 99, no. 3, pp. 345-370, Sep. 2013.
[9] Y. L. Qi, D. G. Musson, B. Schweighardt, T. Tompkins, L. Jesaitis, A. J. Shaywitz, K. Yang, and C. A. O'Neill, "Pharmacokinetic and pharmacodynamic evaluation of Elosulfase Alfa, an enzyme replacement therapy in patients with morquio a syndrome, " Clin. Pharmacokinet., vol. 53, no. 12, pp. 1137-1147, Dec. 2014.
[10] J. D. Murray, Mathematical Biology:Ⅰ. An Introduction. 3rd ed. New York: Springer, 2002.
[11] L. A. Segel, Mathematical Models in Molecular and Cellular Biology. Cambridge: Cambridge University Press, 1980.
[12] G. Varadharajan and L. Rajendran, "Analytical solution of coupled non-linear second order reaction differential equations in enzyme kinetics, " Nat. Sci., vol. 3, no. 6, pp. 459-465, May 2011.
[13] A. Meena, A. Eswari, and L. Rajendran, "Mathematical modelling of enzyme kinetics reaction mechanisms and analytical solutions of non-linear reaction equations, " J. Math. Chem., vol. 48, no. 2, pp. 179-186, Aug. 2010.
[14] P. T. Benavides and U. Diwekar, "Optimal control of biodiesel production in a batch reactor: Part Ⅰ: deterministic control, " Fuel, vol. 94, pp. 211-217, Apr. 2012.
[15] P. K. Roy, S. Nandi, and M. K. Ghosh, "Modeling of a control induced system for product formation in enzyme kinetics, " J. Math. Chem., vol. 51, pp. 2704-2717, Nov. 2013.
[16] F. A. Basir, R. Bhattacharyya, and P. K. Roy, "Delay induced oscillation in a biochemical model and its control, " Nonlin. Stud., vol. 22, no. 3, pp. 453-472, Aug. 2015.
[17] R. A. Azizyan, A. E. Gevorgyan, V. B. Arakelyan, and E. S. Gevorgyan, "Mathematical modeling of uncompetitive inhibition of Bi-substrate enzymatic reactions". Int. Schol. Sci. Res. Innov. , vol.7, no.10, pp.974–977, 2013, 2013.
[18] S. Westerlund, "Dead matter has memory!". Phys. Scrip. , vol.43, no.2, pp.174–179, 1991, 1991. DOI:10.1088/0031-8949/43/2/011
[19] V. E. Tarasov, "Review of some promising fractional physical models, " Int. J. Mod. Phys. B, vol. 27, no. 9, pp. Article No. 1330005, Mar. 2013.
[20] A. A. Stanislavsky, "Memory effects and macroscopic manifestation of randomness, " Phys. Rev. E, vol. 61, no. 5, pp. 4752-4759, May 2000.
[21] J. K. Popović, M. T. Atanacković, A. S. Pilipović, M. R. Rapaić, S. Pilipović, and T. M. Atanacković, "A new approach to the compartmental analysis in pharmacokinetics: fractional time evolution of diclofenac, " J. Pharmacokinet. Pharmacodyn., vol. 37, no. 2, pp. 119-134, Apr. 2010.
[22] R. Toledo-Hernandez, V. Rico-Ramirez, G. A. Iglesias-Silva, and U. M. Diwekar, "A fractional calculus approach to the dynamic optimization of biological reactive systems. Part Ⅰ: fractional models for biological reactions, " Chem. Eng. Sci., vol. 117, pp. 217-228, Sep. 2014.
[23] E. Ahmed and A. S. Elgazzar, "On fractional order differential equations model for nonlocal epidemics, " Phys. A: Statist. Mechan. Appl., vol. 379, no. 2, pp. 607-614, Jun. 2007.
[24] M. L. Du, Z. H. Wang, and H. Y. Hu, "Measuring memory with the order of fractional derivative, " Sci. Rep., vol. 3, pp. Article No. 3431, Dec. 2013.
[25] M. El-Shahed and A. Alsaedi, "The fractional SIRC model and influenza A, " Math. Probl. Eng., vol. 2011, pp. Article No. 480378, Aug. 2011.
[26] S. Abbas, M. Benchohra, G. M. NǴuérékata, and B. A. Silmani, "Darboux problem for fractional-order discontinuous hyperbolic partial differential equations in Banach algebras". Compl. Variabl. Elliptic Equat.:Int. J. , vol.57, no.2-4, pp.337–350, 2012, 2012. DOI:10.1080/17476933.2011.555542
[27] I. Podlubny, Fractional Differential Equations. San Diego: Academic Press, 1999.
[28] F. A. Abdullah, "Using fractional differential equations to model the Michaelis-Menten reaction in a 2-d region containing obstacles". Science Asia , vol.37, no.1, pp.75–78, 2011, 2011. DOI:10.2306/scienceasia1513-1874.2011.37.075
[29] A. Alawneh, "Application of the multistep generalized differential transform method to solve a time-fractional enzyme kinetics". Discr. Dyn. Nat. Soc., vol. 2013, pp. Article No. 592938 , vol.2013, no.592938, 2013.
[30] A. Dokoumetzidis, R. Magin, and P. Macheras, "Fractional kinetics in multi-compartmental systems, " J. Pharmacokinet. Pharmacodyn., vol. 37, no. 5, pp. 507-524, Oct. 2010.
[31] T. Sardar, S. Rana, and J. Chattopadhyay, "A mathematical model of dengue transmission with memory, " Commun. Nonlin. Sci. Numer. Simulat., vol. 22, no. 1-3, pp. 511-525, May 2015.
[32] B. S. Chen, C. Y. Li, B. Wilson, and Y. J. Huang, "Fractional modeling and analysis of coupled MR damping system, " IEEE/CAA J. of Autom. Sinica, vol. 3, no. 3, pp. 288-294, Jul. 2016.
[33] S. Rana, S. Bhattacharya, J. Pal, G. M. NǴuérékata, and J. Chattopadhyay, "Paradox of enrichment: a fractional differential approach with memory, " Phys. A: Statist. Mechan. Appl., vol. 392, no. 17, pp. 3610-3621, Sep. 2013.
[34] M. K. Ghosh, J. Pal, and P. K. Roy, "How memory regulates drug resistant pathogenic bacteria? a mathematical study, " Int. J. Appl. Comput. Math., vol. 3, no. S1, pp. 747-773, Dec. 2017.
[35] V. E. Tarasov, "No violation of the Leibniz rule. No fractional derivative, " Commun. Nonlin. Sci. Numer. Simulat., vol. 18, no. 11, pp. 2945-2948, Nov. 2013.
[36] M. S. Tavazoei, and M. Haeri, "Chaotic attractors in incommensurate fractional order systems, " Phys. D, vol. 237, no. 20, pp. 2628-2637, Oct. 2008.
[37] S. Nandi, M. K. Ghosh, R. Bhattacharya, and P. K. Roy, "Mathematical modeling to optimize the product in enzyme kinetics". Control Cybern. , vol.42, no.2, pp.431–442, 2013, 2013.
[38] S. Schnell and P. K. Maini, "Enzyme kinetics at high enzyme concentration, " Bull. Math. Biol., vol. 62, no. 3, pp. 483-499, May 2000.
[39] S. Schnell and P. K. Maini, "A century of enzyme kinetics: reliability of the KM and vmax estimates, " Comm. Theoret. Biol., vol. 8, no. 2-3, pp. 169-187, Mar-Jan. 2003.