# Optimized Direct Padé and HPM for Solving Equation of Oxygen Diffusion in a Spherical Cell

Authors: M. A. Sandoval-Hernandez; H. Vazquez-Leal; A. Sarmiento-Reyes; U. Filobello-Nino; F. Castro-Gonzalez; A. L. Herrera-May; R. A. Callejas-Molina; R. Ruiz-Gomez
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

This work presents the application of homotopy perturbation method (HPM) and Optimized Direct Padé (ODP) to obtain a handy and easily computable approximate solution of the nonlinear differential equation to model the oxygen diffusion in a spherical cell with nonlinear oxygen uptake kinetics. On one hand, the obtained HPM solution is fully symbolic in terms of the coefficients of the equation, allowing us to use the same solution for different values of the maximum reaction rate, the Michaelis constant, and the permeability of the cell membrane. On the other hand, the numerical experiments show the high accuracy of the proposed ODP solution, resulting in3.58×10-4 as the lowest absolute relative error (A.R.E.) for a set of coefficients. In addition, a novel technique is proposed to reduce the number of algebraic operations during the process of application of ODP method through the use of the Taylor series, which help to simplify the algebraic expressions used. The powerful process to obtain the solution shows that the Optimized Direct Padé and homotopy perturbation method are suitable methods to use.


---

## Body

## 1. Introduction

Michaelis-Menten kinetics describes the rate of enzymatic reactions. This model is valid when the concentration of a substrate is higher than the concentration of the enzymes, and for steady state conditions, that is, when the concentration of the complex enzyme-substrate is constant [1]. There are several works about the Michaelis-Menten oxygen uptake kinetics. For example, in [2], the relation between Michaelis-Menten direct and inverse kinetics, chemical kinetics of approximation, and second-order kinetics is presented. Moreover, in [3], the subject of transport phenomena is approached by means of the study of transport of quantity of motion (viscous flow), transport of energy (heat conduction, convection, and radiation), and matter transport (Diffusion) [4, 5]. Consequently, transport phenomena are employed for solving different problems in the area of sciences such as chemistry, biochemistry, soil physics, meteorology, biology, and semiconductors, disciplines in which the use of Bessel functions, differential equations, and Laplace transform is required [6–10]. In [11], the behaviour of dopamine released from a small iontophoresis electrode and its voltammetric detection by a carbon fiber sensor 100pm away is presented as a basis for developing a new paradigm for measuring dopamine kinetics in intact rat neostriatum. In [11], the presented model was derived from the work of diffusion in a punctual iontophoretic source where the lineal term of absorption is replaced by a nonlinear expression that describes a Michaelis-Menten kinetic given by a constant Vm, the Michaelis-Menten constant Km.In [12], the equation that models oxygen diffusion in a spherical cell, including nonlinear absorption kinetics, is solved by transforming the Lane-Emden equation into its equivalent Volterra integral form, then Adomian decomposition was employed to solve the nonlinear Fredholm-Volterra integral. Furthermore, in [13] a fully symbolic solution was proposed to solve the aforementioned equation by means of a modified Taylor series obtaining an accuracy of 1.66×10-15 as the lowest mean square error without the use of complicated integrals.There are several methods of solutions for solving nonlinear differential equations such as inverse scattering transformation [14, 15], Darboux transformation [16, 17], bilinear method [18, 19], the tanh-function method [20–22], the variable separation approach [23–25], the symmetry method [26, 27], sine-cosine method [28–30], Adomian decomposition method (ADM) [12, 31–33], and homotopy perturbation method (HPM) [34–36]. HPM is based on the use of power series of p, which transforms a differential equation into a set of linear differential equations.In [37] was presented a procedure to apply Padé method to find approximate solutions for nonlinear differential equations, which consist in that the solution of a differential equation can be directly expressed as a rational power series of the independent variable as a Padé approximant. From (6) ODP employs a polynomial-like rational expression as the proposal of approximation of the nonlinear differential equation to be solved. In general terms, it works by means of substituting the rational expression in the differential equation and then regroups the powers of the independent variable. It is important to note that due to the rational expression a large amount of algebraic operations is generated. However, this work proposes a Direct Padé (DP) modification oriented to reduce this algebraic operation by means of the Taylor expansion of the rational function.This paper is organized as follows. In Section2, a brief description of Michaelis-Menten kinetics and its equation is given. In Section 3, a brief description of Lane-Emden equation is presented. The equation to be solved and its boundary conditions are introduced in Section 4. Section 5 describes in detail Optimized Direct Padé. In Section 6, the HPM solution is presented. Section 7 presents the analysis of solution for the case study by HPM and ODP. In Section 8, its numerical simulations, comparisons, and discussion are presented. Lastly, in Section 9 a brief conclusion of this work is presented.


## 2. Michaelis-Menten Equation

General principles of kinetics in chemical reactions are applicable to catalysed reactions by enzymes in living things. However, this shows a characteristic side which is not observed in nonenzymatic catalyst, the substrate saturation, in terms of enzyme molecules active sites occupation [1]. The study of the effect that substrate concentration has over enzyme activity is not simple task; logically thinking, substrate concentration lowers as the reaction increases. A simplification on kinetics experiments consists in measuring initial velocity Vo. If time is short enough, substrate lowering will be minimal and this could be considered, thus, constant. This behaviour is characteristic of most enzymes and was studied by Michaelis and Menten in 1913 [38]. Figure 1 presents the three phases of enzyme kinetics:(i)
For a low substrate concentration, velocity reaction is directly proportional to substrate concentration (linear relation), first-order kinetics.(ii)
For a high substrate concentration, velocity reaction is practically constant and independent of substrate concentration; kinetics is considered zero order.(iii)
For medium concentrations of substrate, velocity of the process becomes nonlinear and this phase is called mixed kinetics.Figure 1
Enzyme kinetics: Michaelis-Menten model.The curve which expresses the relationship between substrate concentration and initial velocity has the same shape for most enzymes, a rectangular hyperbola, whose algebraic expression is given by Michaelis-Menten equation (1), which describes the relationship between initial velocity Vo and substrate concentration, S [1]:(1)Vo=VmaxSKm+S,V m a x, maximum velocity, and Km, Michaelis constant, are two characteristic kinetics parameters of every enzyme that can be experimentally determined.Maximum velocityVmax is obtained when velocity of reaction becomes independent of substrate concentration; when this occurs, velocity reaches a maximum value. This value depends on the available amount of enzyme. Michaelis constant Km indicates substrate concentration in which reaction rate is half of maximum velocity. This parameter is independent of enzyme concentration and it is characteristic of every enzyme according to the employed substrate (or substrates). Constant Km indicates the affinity that the enzyme has for the substrate, being this higher at a lower value of Km. For a lower value Km, the lower the substrate needed to reach half of maximum velocity, thus, the higher the affinity an enzyme has for the substrate.


## 3. Lane-Emden Equation

Lane-Emden equation (2) is a basic nonlinear equation arising in the study of stellar structure theory and has(2)y′′x+mxy′x+fx,y=gx,0<x≤1,m≥0,wherem can be chosen to describe a Cartesian geometry, m=0; cylindrical, m=1; or spherical, m=2 [39]. f(x,y) is a real function and g(x) is an analytical function. Initial conditions are(3)y0=A,y′0=B,whereA and B are constants.Equation (2) is used to model physical and astrophysical phenomena such as theory of stellar structure, spherical cloud of gas thermal behaviour, isothermal gas sphere [40], and transport phenomena [5] as it occurs for the case of particle diffusion in biophysical processes [4] as well in physiological processes.A model of oxygen diffusion and nonlinear absorption in a sphere was originally proposed and resolved in [41]. Later, the same model was proposed, analysed, and solved in [42]. The equation of the model is a steady state diffusion-reaction equation that represents oxygen transport by means of linear diffusion in a sphere.


## 4. Equation of Oxygen Diffusion in a Spherical Cell

The differential equation that models the oxygen diffusion in a spherical cell with Michaelis-Menten oxygen uptake kinetics (4) [41–46] is obtained by substituting (1) into (2) and given by(4)C′′R+2RC′R-αCRK+CR=0,0<R≤1,with boundary conditions(5)C′0=0,C′1+HC1=H.VariablesC and R represent the oxygen concentration and the radial distance, respectively. Parameters α, K, and H are constants that represent the maximum reaction rate. Michaelis constant is the concentration at half saturation and cell membrane permeability, respectively.When boundary conditionR is equal to 0, oxygen distribution is symmetrical to the center of the sphere. Boundary condition at R=1 determines that oxygen flow through the cell membrane is proportional to C′(1)=1-C(1) differential, which is less than the cell membrane normalized concentration [44, 45].


## 5. Optimized Direct Padé

The method ODP [37] considers that a nonlinear differential equation can be expressed as follows:(6)Lu+Nu-fr=0,wherer∈Ωwith boundary conditions(7)Bu,∂u∂η=0,wherer∈Γ,whereL and N are the linear and nonlinear operators, respectively; f(r) is a known analytical function; B is a boundary operator; Γ is the boundary of the domain Ω [47].u 0 is the initial condition for (6) that satisfies the boundary condition and p is the homotopy perturbation parameter. The solution of (6) can be written as(8)ux=a0x1+a1x2+a2x3+⋯+anxn1+b1x1+b2x2+⋯+bmxm,where a0,a1,⋯ and b0,b1,⋯ are unknowns determined by the method; n,m are the numerator and denominator order.There are no systematic criteria to choose the optimal Padé order[n/m] for a given problem. However, generally a finite number of terms are required to obtain highly precise Padé approximation. The basic procedure for Direct Padé can be described as follows:(1)
Boundary conditions are substituted into (8) to generate an equation for each boundary condition. It is important to note that there exists an algebraic equation for each boundary condition.(2)
u(x) from (8) is substituted into (6), then regrouping the resulting equation in terms of x powers. Afterwards, the regrouping procedure includes eliminating denominator terms arising from Padé approximant (8). In this way, the resulting expression is a power series that represent the residual error from (6).(3)
In order to reduce the residual error, from the lowest order, each coefficient fromx powers in the resulting power series is equalized to zero to obtain an algebraic equation in terms of the unknown coefficients of (8).(4)
All the previous steps generate NAES in terms of the unknowns from (8). Finally, NAES are solved to obtain a0,a1,⋯ and b0,b1,⋯.The disadvantage of the method is that it becomes as difficult as the number of terms that the expression has in the numerator as well in the denominator because it is a rational expression. Even more, the substitution into (8) gets more complicated if higher order derivatives are encountered, mainly because the first derivative of a rational expression lies in terms of another rational expression with a squared denominator and if in the nonlinear equation exist several different x powers, another rational expression with even more terms is generated. For this reason, a common denominator has to be found, to formulate the system of equations that allows finding ai and bi. All this process may be too cumbersome and computationally expensive. In this matter, it is concluded that if the nonlinear equation possesses higher order derivatives of x (the independent variable) the number of algebraic operations is notoriously increased when x powers are regrouped to form the equation system to find ai and bi. Thus, in order to avoid the algebraic complexity that implies the substitution of (8) into (6), in this paper we propose a novel strategy to circumvent such issue by obtaining the Taylor expansion of (8).(9)Tux=f0+f11!x+f22!x2+f33!x3+⋯+fqn!xq,whereT(·) is the Taylor expansion of (8) with respect to x and fi(a1,a2,⋯,b0,b1,⋯) considering i=[0,1,2,…,n], and q=m+n-1 is the resulting expansion of (8). The replacement of (8) by (9) decreases enormously the number of algebraic steps when it is substituted into the differential equation obtaining a simpler algebraic procedure. Once (9) is replaced into (6), with regrouping in terms of x powers, a system of equations is constructed with the aim of finding ai (i=[1,2,⋯,n]) and bk(k=[1,2,⋯,m]).


## 6. HPM Solution

Multiplying Michaelis-Menten kinetic equation (4) by [K+C(R)], the following is obtained:(10)Kd2dR2CR+CRd2dR2CR+2KRddRCR+2CRRddRCR-αCR=0.The proposed homotopy perturbation formulation is(11)H=1-pd2dR2CR+2RddRCR+pKd2dR2CR+CRd2dR2CR+2KRddRCR+2CRRddRCR-αCR.Assuming the solution for (11) is given by a power series of p(12)CR=C0R+C1Rp+C2Rp2+⋯,substituting (12) into (11), and regrouping in terms of p,(13)p0:d2dR2C0R=0C00=cC0′0=0p1:RK+C0R-1d2dR2C0R+Rd2dR2C1R+2K+2C0R-2ddRC0R-RαC0R+2ddRC1R=0,C10=0C1′0=0,⋮Substituting the first two system solutions (13) into (12), we obtain(14)CR=A+B+16AαR2.This is the solution employing HPM, whereA and B are constants to be determined. By substituting (14) into the second boundary condition of (5), i.e., C′(1)+HC(1)=H, we obtain(15)13Aα+HA+16HAα+HB-1=0.Solving forA in (15),(16)A=-6HB-16H+Hα+2αSubstituting in (14),(17)CR=B-6HB-16H+Hα+2α-HB-1αR26H+Hα+2α.With the end of reducing the residual error, (17) is substituted into (4), at middle point R=0.5, between boundaries in order to cancel the error at that point, to obtain for B(18)B=-112Hα3H+824KH2α+6H2α+144KH2+3H2α2+48KHα+144H2+14Hα2+16α2--4032H3α2+448Hα4-828H4α2-36H4α3+1728H4α-168H3α3-192H2α3+20736K2H4+9H4α4-84H4α4+292H2α4+2112KH2α3+960KH3α3+19008KH4α+2304K2H3α2+144KH4α3+6912K2H4α+2880KH4α2+576K2H4α2+12672KH3α2+1536KHα3+41472KH3α+2304K2H2α2+13824KH2α2+13824K2H3α+20736H4+256α4-4608H2α21/2.This strategy seeks to reduce the approximation error. Substituting (18) into (17), we obtain the solution in function of α, K, and H.(19)CR=A1+A2R2,whereA1 and A2 are given for(20)A1=B-6HB-16H+Hα+2α,A2=-HB-1α6H+Hα+2α.


## 7. ODP Solution Procedure

The following rational solution is now proposed for (4):(21)CR=a0+a1R+a2R21+b1R+b2R2In order to reduce the amount of algebraic steps, it is proposed in this work, and for the first time, to expand the rational expression (21) and employ the Taylor series of order four:(22)CR=a0+a1-a0b1R+a2-a0b2+-a1+a0b1b1R2+-a1+a0b1b2+-a2+a0b2+b1a1-a0b12b1R3+-a2+a0b2+b1a1-a0b12b2+b-2a1-2b2a0b1+b1a2-b12a1+a0b13b1R4,differentiating with respect toR and to satisfy the first boundary with R=0 yields(23)a1+a0b1=0.With expansion (22) we have the advantage of eliminating the denominator and thus reducing the amount of algebraic terms when substituting (21) into (4). To satisfy the second boundary condition using (21) and (5)(24)a1+2a21+b1+b2+a0+a1+a2b1+2b21+b1+b22+Ha0+a1+a21+b1+b2-H=0.Substituting fourth-order Taylor expansion of (21) into (10),(25)32a2-2a0b2+2-a1+a0b1b1K+32a2-2a0b2+2-a1+a0b1b1a0+2a1-a0b12-αa0=0,6-a1+a0b1b2+6-a2+a0b2+b1a1-a0b12b1K+32a2-2a0b2+2-a1+a0b1b1a1-a0b1+6-a1+a0b1b2+6-a2+a0b2+b1a1-a0b12b1a0+23-a1+a0b1b2+3-a2+a0b2+b1a1-a0b12b1K+2a1-a0b1a2-a0b2+-a1+a0b1b1+23-a1+a0b1b2+3-a2+a0b2+b1a1-a0b12b1a0-αa1-a0b1=0,12-a2+a0b2+b1a1-a0b12b2+12b2a1-2b2a0b1+b1a2-b12a1+a0b13b1K+32a2-2a0b2+2-a1+a0b1b1a2-a0b2+-a1+a0b1b1+31,6-a1+a0b1b2+6-a2+a0b2+b1a1-a0b12b1a1-a0b1+12-a2+a0b2+b1a1-a0b12b2+12b2a1-2b2a0b1+b1a2-b12a1+a0b13b1a0+24-a2+a0b2+b1a1-a0b12b2+4b2a1-2b2a0b1+b1a2-b12a1+a0b13b1K+2a1-a0b1-a1+a0b1b2+-a2+a0b2+b1a1-a0b12b1+23-a1+a0b1b2+3-a2+a0b2+b1a1-a0b12b1a1-a0b1+24-a2+a0b2+b1a1-a0b12b2+4b2a1-2b2a0b1+b1a2-b12a1+a0b13b1a0-αa2-a0b2+-a1+a0b1b1=0.Substituting the values ofα, K, H and solving the simultaneous system of equations given by (23), (24), and (25), the numerical values of a0, a1, a2, b1, and b2 are obtained to satisfy (21).


## 8. Numerical Simulations and Discussion

For all the simulations presented in this work, with the purpose of comparison, Maple17 numerical routines were employed. For example, numerical values of α=2.5, K=2, H=4 were used, with (19) and (21) expressed with the following values:(26)CR=0.8143254327+0.1237830449R2,(27)CR=0.8116575042+0.09461382624R21-0.03162388156R2.In order to verify the precision of the solutions (26) and ODP (27) obtained by HPM, their performance is compared with ADM solutions (φ2(R) and φ3(R)) reported in [12], as well against the numerical solution given by Maple 17 numerical routines.(28)φ2R=1+α-2-H+HR26HK+1,(29)φ3R=1+α-2-H+HR26HK+1+α2K40+28H+7H2+-20H-10H2R2+3HR4360H2K+13,and (13) given by MTSM [13] is(30)CTR=c+16αcK+cR2,(31)c=112HH+22α2+12HK-1H+2α+36H2K+12+-6K-α+6H-2α,H≠0.Equations (19) and (21) are second-order compact expressions that present a good accuracy at the interval R∈[0,1]. Figure 2(a) shows, in a graphical manner, that the obtained solutions with ODP and HPM present a better accuracy than the obtained for (28) and (30), (31), for values of α=2.5, K=2,H=4. Additionally, Figures 2(b) and 2(c) show the same comparisons with the values for α=3.5, K=2, H=4 and α=4.5, K=3, H=4.Figure 2
Comparison of HPM (19) and ODP (21) versus numerical solution, (28), (29), and (30).
(a)
Case forα = 2.5, K=2, H=4 (b)
Case forα = 3.5, K=2, H=4 (c)
Case forα = 4.5, K=3, H=4 (d)
Case forα = 8.5, K=5, H=10As for the previous cases, in Figure2(d), for higher values of α=8.5,K=5,H=10, the obtained solution employing ODP has better accuracy than other solutions. For this case, (32) and (33) describe the algebraic solutions for this value for HPM and ODP, respectively. For the values mentioned before, Figure 2(d) presents the graphics for the solutions employing ODP and HPM compared to the solutions proposed by MTSM and φ3.(32)CR=0.7624971978+0.1979190019R2,(33)CR=0.7581194350+0.1379309752R21-0.06409108102R2.In Table1 is presented and compared the absolute relative error (A.R.E.) for each approximation according to the values of α, K, H, for HPM and ODP solutions, including (28), (29), and (30).Table 1
Absolute relative error comparison.
Case ODP HPM MTSM [13] φ 3 ( R ) [12] φ 2 ( R ) [12] α = 2.5, K=2, H=4 3.58 × 1 0 - 4 3.21 × 1 0 - 3 6.99 × 1 0 - 3 1.11×10-2 4.54 × 1 0 - 1  α = 3.5, K=2, H=4 1.00 × 1 0 - 3 6.26 × 1 0 - 3 1.37 × 1 0 - 2 5.16×10-2 5.06 × 1 0 - 2  α = 4.5, K=3, H=4 8.37 × 1 0 - 4 6.18 × 1 0 - 3 1.34 × 1 0 - 2 4.04×10-2 2.90 × 1 0 - 1  α = 8.5, K=5, H=10 9.59 × 1 0 - 4 5.90 × 1 0 - 3 1.46 × 1 0 - 2 2.39×10-2 1.665×100The absolute relative error is defined as the absolute difference of the areas under the curve for the numerical and approximate solutions, divided by the area of the exact solution. The employed integrating method was the Rule of Simpson1/3.In all of the presented cases, the A.R.E. obtained with ODP is much lower compared to the Adomian and MTSM methods reported in [12, 13] and matches with all the graphics presented in this work.In the same manner, for the cases (α=4.5,K=3,H=5) and (α=8.5,K=5,H=10) for HPM it is observed that the relative absolute error is lower than the ones presented in Adomian and MTSM. Besides, the solution obtained by HPM has the advantage of being simpler because it only consists of two terms, compared against (28), (29), and (30) that possess more terms.Thus, for Optimized Direct Padé and HPM, the proposed results in this work coincide with [12, 13, 39, 48, 49]; that is, the oxygen concentration is higher at the substrate and diminishes as it approaches the center of spherical cell, which is due to the oxygen diffusing into the cell; it reacts to give rise to the boundary condition at the center. Likewise, results are congruent with the spherical geometry proposed by Lane-Emden (2).Finally, it is important to remark that although the precision for ODP tends to be the same as DP, it has the advantage of having a simpler algebraic procedure, avoiding the algebraic complexity of substituting a rational expression that generates more rational expressions with higher powers in the denominator o in its case, more numerator terms, making the formulation of the system of equations difficult for findingai and bi from the rational expression of the solution.


## 9. Conclusions

In this article were proposed two compact solutions for Michaelis-Menten oxygen uptake kinetics equation employing HPM and ODP. ODP presents a quotient of two polynomials, one of them of second order divided by a first-order one. On the other hand, HPM presents only two algebraic terms, the independent term and the second degree term. Likewise, the solutions found in this work were compared to those reported by other authors, showing that ODP has a better accuracy. HPM solution shows a good approximation in comparison to (29) obtained by Adomian; despite the fact that it only has two terms, it concluded that the solution is in good agreement. Furthermore, in this paper, we propose to replace the rational expression (8) by its Taylor series in order to reduce the amount of algebraic steps during the application of DP. The graphs presented in this work show that the solution given by ODP and HPM tends to converge as we approach R=1 unlike Adomian φ3(R) which tends to diverge by increasing the absolute error. For the case of α=8.5, K=5, and H=10 with ODP, the A.R.E. was 9.59×10-4. Also the HPM and ODP methods were compared against the numerical one obtaining very good approximations. In addition, ODP requires fewer iterations to obtain a more compact and accurate expression than the solution proposed by the Adomian method (29).



---

*Source: 9142124-2018-09-02.xml*