Application of the Multistep Generalized Differential Transform Method to Solve a Time-Fractional Enzyme Kinetics

Abstract

The multistep differential transform method is first employed to solve a time-fractional enzyme kinetics. This enzyme-substrate reaction is formed by a system of nonlinear ordinary differential equations of fractional order. The fractional derivatives are described in the Caputo sense. A comparative study between the new algorithm and the classical Runge-Kutta method is presented in the case of integer-order derivatives. The results demonstrate reliability and efficiency of the algorithm developed.

1. Introduction

Mathematical modeling of complex processes is a major challenge for contemporary scientists. In contrast to simple classical systems, where the theory of integer-order differential equations is sufficient to describe their dynamics, complex systems are characterized by the variability of structures in them, multiscale behavior, and nonlinearity in the mathematical description of the mutual relationship between parameters [1]. Fractional derivatives provide an excellent instrument for the description of dynamical behavior of various complex materials and systems. There are two important benefits: (1) we have more degrees of freedom in the model considered and (2) a “memory” is enclosed in the model. Therefore, the number of scientific and engineering problems involving fractional derivatives is already very large and still growing and perhaps the fractional calculus (i.e., derivatives and integrals of any real or complex order) will be the calculus of the twenty-first century [1].

Recent investigations have shown that many complex biological systems can be represented more accurately through fractional derivative formulation. Magin [2] was the first who used fractional derivatives and fractional integrals in order to model stress-strain relationship in biomaterials. Craiem et al. [3] applied fractional calculus to model arterial viscoelasticity. Abdullah [4] used fractional differential equations to model the Michaelis-Menten reaction in a 2d region containing obstacles. He investigated the phenomenon of subdiffusion, which occurs when there is molecular crowding, by proposing a continuous spatial model involving fractional differential equations. The Michaelis-Menten mechanism is the simplest chemical network which models the formation of a product through an enzymatic catalysis of a substrate. In particular, an enzyme reacts with the substrate and reversibly forms an intermediate complex, which then decays into the product and original enzyme. This enzyme reaction model is given by (for more details on the enzyme dynamics, see [5–9] and the references therein) 𝐸+𝑆⟷^𝑘1_𝑘−1𝐸𝑆−−→^𝑘2𝐸+𝑃, (1) where 𝐸 is the enzyme, 𝑆 the substrate, 𝐸𝑆 the enzyme-substrate intermediate complex, and 𝑃 the product. The parameters 𝑘₁, 𝑘₋₁, and 𝑘₂ are positive rate constants and denote the rates of reaction of these three processes. Note that substrate binding is reversible but product release is not. The concentration of the reactants in (1) is denoted by lower case letters: 𝑠=[𝑆],𝑒=[𝐸],𝑐=[𝑆𝐸],𝑝=[𝑃]. (2) Applying the law of mass action, which states that reaction rates are proportional to the concentrations of the reactants, the time evolution of scheme (1) can be determined from the solution of the system of coupled nonlinear ordinary differential equations (ODEs) [6, 7]: 𝑑𝑠𝑑𝑡=−𝑘₁𝑒𝑠+𝑘₋₁𝑐,𝑑𝑒𝑑𝑡=−𝑘₁𝑒𝑠+𝑘₋₁+𝑘₂𝑐,𝑑𝑐𝑑𝑡=𝑘₁𝑒𝑠−𝑘₋₁+𝑘₂𝑐,𝑑𝑝𝑑𝑡=𝑘₂𝑐, (3) subject to the initial conditions 𝑠(0)=𝑠₀,𝑒(0)=𝑒₀,𝑐(0)=0,𝑝(0)=0. (4) Systems (3) can be reduced to only two equations for 𝑠 and 𝑐, namely [6], 𝑑𝑠𝑑𝑡=−𝑘₁𝑒₀𝑠+𝑘₁𝑠+𝑘₋₁𝑐,𝑑𝑐𝑑𝑡=𝑘₁𝑒₀𝑠−𝑘₁𝑠+𝑘₋₁+𝑘₂𝑐, (5) subject to the initial conditions 𝑠(0)=𝑠₀ and 𝑐(0)=0.By introducing the following parameters: 𝑘𝜏=¹𝑒₀𝑡𝜀(𝜏)=𝑠(𝑡),𝑢𝑠₀(𝜏)=𝑐(𝑡),𝑣𝑐₀,𝑤(𝜏)=𝑝(𝑡)𝑒₀𝑘,𝜆=²𝑘₁𝑠₀,𝑘₃=𝑘₋₁+𝑘₂𝑘₁𝑠₀𝑒,𝜀=⁰𝑠₀. (6)

Then, (5) and the initial conditions (4) can be represented in dimensionless form as follows: 𝑑𝑢𝑑𝑡=−𝜀𝑢+𝜀𝑢+𝑘₃−𝜆𝑣,𝑑𝑣𝑑𝑡=𝑢−𝑢+𝑘₃𝑣,𝑑𝑤𝑢(0)(0)(0)𝑑𝑡=𝜆𝑣,=1,𝑣=0,𝑤=0, (7) where 𝜀, 𝑘, and 𝜆 are dimensionless parameters.

Now, we introduce fractional order into enzyme reaction model. The new system is described by the following set of fractional differential equations: 𝐷^𝛼1𝑢=−𝜀𝑢+𝜀𝑢+𝑘₃𝐷−𝜆𝑣,_𝛼2𝑣=𝑢−𝑢+𝑘₃𝐷𝑣,_𝛼3𝑢(0)(0)(0)𝑤=𝜆𝑣,=1,𝑣=0,𝑤=0, (8)(9)(10)(11) where the fractional derivatives 𝐷^𝛼1𝑢, 𝐷^𝛼2𝑣, and 𝐷^𝛼3𝑤 are considered in the Caputo sense, and 0<𝛼_𝑖≤1, 𝑖=1,2,3, are parameters describing the order of the fractional time-derivatives in the Caputo sense. The general response expression contains parameters describing the order of the fractional derivatives that can be varied to obtain various responses. Obviously, the classical integer-order enzyme reaction model can be viewed as a special case from the fractional-order enzyme reaction model by setting 𝛼₁=𝛼₂=𝛼₃=1. In other words, the ultimate behavior of the fractional system response must converge to the response of the integer-order version of the equation [10].

Our motivation for this work is to obtain the approximate solution of the time-fractional enzyme kinetics (8)–(10) using the multistep generalized differential transform method (MSGDTM). This method is only a simple modification of the generalized differential transform method (GDTM) [11–14], in which it is treated as an algorithm in a sequence of small intervals (i.e., time step) for finding accurate approximate solutions to the corresponding systems. The approximate solutions obtained by using GDTM are valid only for a short time, while the ones obtained by using the MSGDTM [15, 16] are more valid and accurate during a long time and are in good agreement with the RK4-5 numerical solution when the order of the derivative 𝛼_𝑖=1,   𝑖=1,2,3.

This paper is organized as follows. In Section 2, we present some necessary definitions and notations related to fractional calculus and differential transform method. In Section 3, the proposed method is described. In Section 4, the method is applied to the problems (8)–(10), and numerical simulations are presented graphically. Finally, the conclusions are given in Section 5.

2. Preliminaries

In this section, we give some basic definitions and properties of the fractional calculus theory and differential transform method which are used further in this paper [1, 10, 17–25].

2.1. Fractional Calculus

Definition 1. A function 𝑓(𝑥)(𝑥>0) is said to be in the space 𝐶_𝛼(𝛼∈𝑅) if it can be written as 𝑓(𝑥)=𝑥^𝑝𝑓₁(𝑥) for some 𝑝>𝛼, where 𝑓₁(𝑥) is continuous in [0,∞), and it is said to be in the space 𝐶^𝑚_𝛼 if 𝑓^(𝑚)∈𝐶_𝛼, 𝑚∈𝐍.

Definition 2. The Riemann-Liouville integral operator of order 𝛼>0 with 𝑎≥0 is defined as 𝐽^𝛼_𝑎𝑓(𝑥)=1Γ(𝛼)^𝑥_𝑎(𝑥−𝑡)^𝛼−1𝑓(𝑡)𝐽𝑑𝑡,𝑥>𝑎,_0𝑎𝑓(𝑥)=𝑓(𝑥). (12) Properties of the operator can be found in [1, 10]. We only need here the following: for 𝑓∈𝐶_𝛼,𝛼,𝛽>0,  𝑎≥0,  𝑐∈𝑅 and 𝛾>−1, we have 𝐽^𝛼_𝑎𝐽^𝛽_𝑎𝑓(𝑥)=𝐽^𝛽_𝑎𝐽^𝛼_𝑎𝑓(𝑥)=𝐽_𝑎^𝛼+𝛽𝑓(𝑥),𝐽^𝛼_𝑎𝑥^𝛾=𝑥^𝛾+𝛼Γ(𝛼)𝐵_{(𝑥−𝑎)/𝑥}𝛼,𝛾+1, (13) where 𝐵_𝜏(𝛼,𝛾+1) is the incomplete beta function which is defined as 𝐵_𝜏𝛼,𝛾+1=^𝜏₀𝑡^𝛼−1(1−𝑡)^𝛾𝐽𝑑𝑡,_𝛼𝑎𝑒^𝑐𝑥=𝑒^𝑎𝑐(𝑥−𝑎)^𝛼∞_𝑘=0[𝑐(𝑥−𝑎)]^𝑘Γ().𝛼+𝑘+1 (14) The Riemann-Liouville derivative has certain disadvantages when trying to model real-world phenomena with fractional differential equations. Therefore, we will introduce a modified fractional differential operator 𝐷^𝛼_𝑎 proposed by Caputo in his work on the theory of viscoelasticity.

Definition 3. The Caputo fractional derivative of 𝑓(𝑥) of order 𝛼>0 with 𝑎≥0 is defined as 𝐷^𝛼_𝑎𝑓(𝑥)=𝐽_𝑎^𝑚−𝛼𝑓^(𝑚)(𝑥)=1Γ(𝑚−𝛼)^𝑥_𝑎𝑓^(𝑚)(𝑡)(𝑥−𝑡)^{𝛼+1−𝑚}𝑑𝑡, (15) for 𝑚−1<𝛼≤𝑚,  𝑚∈𝐍,𝑥≥𝑎,𝑓(𝑥)∈𝐶^𝑚₋₁.

The Caputo fractional derivative was investigated by many authors; for 𝑚−1<𝛼≤𝑚, 𝑓(𝑥)∈𝐶^𝑚_𝛼 and 𝛼≥−1, we have 𝐽^𝛼_𝑎𝐷^𝛼_𝑎𝑓(𝑥)=𝐽^𝑚𝐷^𝑚𝑓(𝑥)=𝑓(𝑥)−^𝑚−1_𝑘=0𝑓^(𝑘)(𝑎)(𝑥−𝑎)^𝑘.𝑘! (16) For mathematical properties of fractional derivatives and integrals, one can consult the mentioned references.

2.2. The Definitions and Operations of Differential Transform

The differential transform technique is one of the seminumerical analytical methods for ordinary and partial differential equations that uses the form of polynomials as approximations of the exact solutions that are sufficiently differentiable. The basic definition and the fundamental theorems of the differential transform method (DTM) and its applicability for various kinds of differential equations are given in [21–23]. For convenience of the reader, we present a review of the DTM. The differential transform of the 𝑘th derivative of function 𝑓(𝑡) is defined as follows: 𝐹(𝑘)=1𝑑𝑘!^𝑘𝑓(𝑡)𝑑𝑡^𝑘_𝑡=𝑡0, (17) where 𝑓(𝑡) is the original function and 𝐹(𝑘) is the transformed function. The differential inverse transform of 𝐹(𝑘) is defined as: 𝑓(𝑡)=^∞_𝑘=0𝐹(𝑘)𝑡−𝑡₀^𝑘. (18) From (17) and (18), we get 𝑓(𝑧)=^∞_𝑘=0𝑧−𝑧₀^𝑘𝑑𝑘!^𝑘𝑓(𝑧)𝑑𝑧^𝑘_𝑧=𝑧0, (19) which implies that the concept of differential transform is derived from Taylor series expansion, but the method does not evaluate the derivatives symbolically. However, relative derivatives are calculated by an iterative way which is described by the transformed equations of the original function. For implementation purposes, the function 𝑓(𝑡) is expressed by a finite series, and (18) can be written as 𝑓(𝑡)≈^𝑁_𝑘=0𝐹(𝑘)𝑡−𝑡₀^𝑘. (20) Here, 𝑁 is decided by the convergence of natural frequency. The fundamental operations performed by differential transform can readily be obtained and are listed in Table 1. The main steps of the DTM, as a tool for solving different classes of nonlinear problems, are the following. First, we apply the differential transform (17) to the given problem (integral equation, ordinary differential equation, or partial differential equations), and then the result is a recurrence relation. Second, solving this relation and using the differential inverse transform (18), we can obtain the solution of the problem.

Table 1 

Operations of differential transformation.

3. Solving the Systems (8)–(10) Using the Multistep Generalized Differential Transform Method (MSGDTM)

It has been shown that the approximated solutions obtained using DTM are not valid for large 𝑡 for some systems [11–14]. Therefore, we use the multistep generalized differential transform method (MSGDTM) to solve the system of equations (8)–(10). The method is only a simple modification of the generalized differential transform method, in which it is treated as an algorithm in a sequence of small intervals (i.e., time step) for finding accurate approximate solutions to the corresponding systems [15, 16]. This MSGDTM offers accurate solutions over a longer time frame (more stable) compared to the standard GDTM.

Taking the differential transform of (8)–(10) with respect to time 𝑡 gives 𝑈()=Γ()𝑘+1𝛼𝑘+1Γ(𝛼())×(𝑘)𝑘+1+1−𝜀𝑈+𝜀^𝑘_𝑙=0𝑈(𝑙)𝑉(𝑘−𝑙)+𝑘₃𝑉(𝑘),𝑉()=Γ()−𝜆𝑘+1𝛼𝑘+1Γ(𝛼())×𝑈(𝑘)−𝑘+1+1_𝑘𝑙=0𝑈(𝑙)𝑉(𝑘−𝑙)−𝑘₃𝑉(𝑘,𝑊()=Γ())𝑘+1𝛼𝑘+1Γ(𝛼())(𝑘),𝑘+1+1𝜆𝑉 (21) where 𝑈(𝑘), 𝑉(𝑘), and 𝑊(𝑘) are the differential transformation forms of 𝑢(𝑡), 𝑣(𝑡), and 𝑤(𝑡), respectively. The differential transform forms of the initial conditions are given by 𝑈(0)=1, 𝑉(0)=0, and 𝑊(0)=0. In view of the differential inverse transform, the differential transform series solution for the systems (8)–(10) can be obtained as 𝑢(𝑡)=^𝑁_𝑛=0𝑈(𝑛)𝑡^𝛼𝑛,𝑣(𝑡)=^𝑁_𝑛=0𝑉(𝑛)𝑡^𝛼𝑛,𝑤(𝑡)=^𝑁_𝑛=0𝑊(𝑛)𝑡^𝛼𝑛. (22) Now, according to the multistep generalized differential transform method, the series solution for the systems (8)–(10) is suggested by 𝑢(𝑡)=^𝐾_𝑛∑₌₀𝑈₁(𝑛)𝑡^𝛼𝑛,𝑡∈0,𝑡₁,_𝐾𝑛∑₌₀𝑈₂(𝑛)𝑡−𝑡₁^𝛼𝑛,𝑡∈𝑡₁,𝑡₂,⋮_𝐾𝑛∑₌₀𝑈_𝑀(𝑛)𝑡−𝑡_𝑀−1^𝛼𝑛,𝑡∈𝑡_𝑀−1,𝑡_𝑀,𝑣(𝑡)=^𝐾_𝑛∑₌₀𝑉₁(𝑛)𝑡^𝛼𝑛,𝑡∈0,𝑡₁,_𝐾𝑛∑₌₀𝑉₂(𝑛)𝑡−𝑡₁^𝛼𝑛,𝑡∈𝑡₁,𝑡₂,⋮_𝐾𝑛∑₌₀𝑉_𝑀(𝑛)𝑡−𝑡_𝑀−1^𝛼𝑛,𝑡∈𝑡_𝑀−1,𝑡_𝑀,𝑤(𝑡)=^𝐾_𝑛∑₌₀𝑊₁(𝑛)𝑡^𝛼𝑛,𝑡∈0,𝑡₁,_𝐾𝑛∑₌₀𝑊₂(𝑛)𝑡−𝑡₁^𝛼𝑛,𝑡∈𝑡₁,𝑡₂,⋮_𝐾𝑛∑₌₀𝑊_𝑀(𝑛)𝑡−𝑡_𝑀−1^𝛼𝑛,𝑡∈𝑡_𝑀−1,𝑡_𝑀, (23) where 𝑈_𝑖(𝑛), 𝑉_𝑖(𝑛), and 𝑊_𝑖(𝑛) for 𝑖=1,2,…,𝑀 satisfy the following recurrence relations: 𝑈_𝑖()=Γ()𝑘+1𝛼𝑘+1Γ(𝛼())×𝑘+1+1−𝜀𝑈_𝑖(𝑘)+𝜀^𝑘_𝑙=0𝑈_𝑖(𝑙)𝑉_𝑖(𝑘−𝑙)+𝑘₃−𝜆𝑉_𝑖(𝑘),𝑉_𝑖()=Γ()𝑘+1𝛼𝑘+1Γ(𝛼())×𝑈𝑘+1+1_𝑖(𝑘)−^𝑘_𝑙=0𝑈_𝑖(𝑙)𝑉_𝑖(𝑘−𝑙)−𝑘₃𝑉_𝑖(𝑘,𝑊_𝑖)()=Γ()𝑘+1𝛼𝑘+1Γ(𝛼(𝑘+1)+1)𝜆𝑉_𝑖(𝑘), (24) such that 𝑈_𝑖(0)=𝑈_𝑖−1(0), 𝑉_𝑖(0)=𝑉_𝑖−1(0), and 𝑊_𝑖(0)=𝑊_𝑖−1(0). Finally, if we start with 𝑈₀(0)=1, 𝑉₀(0)=0, and 𝑊₀(0)=0, using the recurrence relation given in (24), then we can obtain the multistep solution given in (23).

This multistage MSGDTM offers accurate solutions over a longer time frame (more stable) compared to the standard GDTM [11–14]. This distinctive strategy grants the iterative algorithm a time-marching scheme which significantly drives forward the convergence of the solutions precisely with great rapidity. The value of 𝑡 increases according to the designated time-step size in each iteration computation. Each amplified time step will produce a new approximation value for the desired iteration step, and the process is continued until the targeted time frame is achieved. Errors are therefore minimized and depend on the time step Δ𝑡, but the main disadvantage lies in having longer computational time. Generally, the convergence is achieved more quickly as the iteration steps increase.

4. Numerical Results

In this work, we carefully propose the MSGDTM, a reliable modification of the GDTM, that improves the convergence of the series solution. The method provides immediate and visible symbolic terms of analytic solutions, as well as numerical approximate solutions to both linear and nonlinear differential equations. To demonstrate the effectiveness of the proposed algorithm as an approximate tool for solving the nonlinear system of fractional differential equations (8)–(10) for larger 𝑡, we apply the proposed algorithm on the interval [0,10]. We divide this interval to subintervals with time step Δ𝑡=0.1 or 𝑀=100 and the tem-number of GDTM series solutions is fixed, 𝐾=10. All the results are calculated by using the computer algebra package Mathematica.

Figures 1 and 2 show the approximate solutions obtained using the MSGDTM and the fourth-order Runge-Kutta method of concentrations of substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for various values of dimensionless reaction parameters 𝑘₃, 𝜆, and 𝜀 when 𝛼₁=𝛼₂=𝛼₃=1. It can be seen that the results from the MSGDTM match the results of the Runge-Kutta method very well, which implies that the MSGDTM can predict the behavior of these variables accurately for the region under consideration. Therefore, the proposed method is a very efficient and accurate method that can be used to provide analytical solutions for nonlinear systems of differential equations.

Figure 1 

Profile of the normalized concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for 𝑘=1, 𝜆=0.5, 𝜀=0.6, and 𝛼1=𝛼2=𝛼3=1. MSDTM solution (solid line); Runge-Kutta method solution (dotted line).

Figure 2 

Profile of the normalized concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for 𝑘3=2, 𝜆=1.0, 𝜀=0.5, and 𝛼1=𝛼2=𝛼3=1. MSDTM solution (solid line); Runge-Kutta method solution (dotted line).

Next, we consider how the concentrations of substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 depend upon the magnitude of the order of fractional derivatives. We fix the dimensionless reaction parameters 𝑘₃, 𝜆, and 𝜀 and perform the numerical simulation for different values of  𝛼₁, 𝛼₂, and 𝛼₃. Simulation results are presented in Figures 3, 4, 5, and 6, and it is clear that the concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 continuously depend on the fractional derivatives. For example, the concentration of substrate 𝑢 decreases gradually from its initial value of the concentration (𝑢(0)=1) and converges to zero when 𝛼₁=𝛼₂=𝛼₃=1 and decreases slowly to zero when 𝛼₁=0.9, 𝛼₂=0.85, and 𝛼₃=0.95. However, the product 𝑤 increases rapidly for the integer-order derivatives and increases slowly in the case of fractional-order derivatives as shown in Figures 3–6. The key finding of these graphs is that the generalized enzyme kinetics model with fractional derivatives has more degrees of freedom and therefore can be varied to obtain various responses of the concentrations of substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤.

Figure 3 

Profile of the normalized concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for 𝑘3=1, 𝜆=0.5, and 𝜀=0.6, 𝛼1=0.98, 𝛼2=0.98, and 𝛼3=0.98 (solid line); 𝛼1=𝛼2=𝛼3=1 (dotted line).

Figure 4 

Profile of the normalized concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for 𝑘3=1, 𝜆=0.5, and 𝜀=0.6, 𝛼1=1.0, 𝛼2=0.95, and 𝛼3=0.95 (solid line); 𝛼1=𝛼2=𝛼3=1 (dotted line).

Figure 5 

Profile of the normalized concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for 𝑘3=1, 𝜆=0.5, and 𝜀=0.6, 𝛼1=0.95, 𝛼2=1.0, and 𝛼3=0.9 (solid line); 𝛼1=𝛼2=𝛼3=1 (dotted line).

Figure 6 

Profile of the normalized concentrations of the substrate 𝑢, enzyme-substrate complex 𝑣, and product 𝑤 for 𝑘3=1, 𝜆=0.5, and 𝜀=0.6, 𝛼1=0.9, 𝛼2=0.85, and 𝛼3=0.95 (solid line); 𝛼1=𝛼2=𝛼3=1 (dotted line).

5. Conclusions

In this paper, the application of MSGDTM was extended to obtain explicit and numerical solutions of a time-fractional enzyme kinetics. The MSGDTM was clearly a very efficient and powerful technique in finding the solutions of the proposed equations. The obtained results demonstrate the reliability of the algorithm and its wider applicability to fractional nonlinear evolution equations.

Acknowledgment

The author gratefully acknowledges the financial support received from The University of Jordan during the sabbatical leave from Department of Mathematics, The University of Jordan, Amman, Jordan.