Mathematical Modeling and Analysis of the Kinetics of Thermal Inactivation of Enzyme

Abstract

A theoretical model of Illeova et al. (2003) thermal inactivation of urease is discussed. Analytical expressions pertaining to the molar concentrations of the native and denatured enzyme are obtained in terms of second-order reaction rate constant. Simple and closed form of theoretical expression pertains to the temperature are also derived. In this paper, homotopy analysis method (HAM) is used to obtain approximate solutions for a nonlinear ordinary differential equation. The obtained approximate result in comparison with the numerical ones is found to be in satisfactory agreement.

1. Introduction

Urease is a good catalyst for the hydrolysis of urea. Several excellent techniques are available to assess urease activity [1, 2]. In 1926, urease was isolated by Summner from the seeds of jack bean as a pure, crystalline enzyme [3]. These crystals, the first obtained for a known enzyme, played a decisive role in proving the protein nature of enzymes. Approximately 50 years later, jack bean urease was identified as the first nickel metalloenzyme [4]. A method for the determination of mercury (II) ions at trace levels is described. The method is based on the profound inhibitory effect of mercury on the enzyme urease [5]. For unknown reasons some seeds are particularly rich sources of urease, and this enzyme has been extensively studied in seeds of various Leguminosae, Cucurbitaceae, Asteraceae and Pinaceae [6, 7]. Jack bean urease, which is the most widely used plant urease, is a nickel containing oligomeric enzyme exhibiting a high degree of specificity to urea [8]. Numerous papers have been published on the applications of urease in biotechnology, including the determination of urea for analytical and biomedical purposes and analysis of heavy metal content in natural drinking water and surface water [9]. Hirai et al. [10] study the structural change of jack bean urease induced by addition of synchrotron radiation. Lencki et al. [11] discuss the effect of subunit dissociation, denaturation, aggregation, coagulation, and decomposition on enzyme inactivation kinetics. Omar and Beauregard [12] investigate the unfolding of jack bean urease by fluorescence emission spectroscopy.

To our knowledge, no rigorous analytical expressions of molar concentrations of the native enzyme, denatured enzyme, and temperature for thermal inactivation of urease for the parameters 𝑘₊₁, 𝑘^′₋₁, 𝑘^′₃, 𝑘^′₂, ℎ, 𝑇_𝐵, and 𝐾 have been reported. The purpose of this communication is to derive simple approximate analytical expression for the nonsteady-state concentrations for thermal inactivation of urease using Homotopy analysis method.

2. Mathematical Formulation of Analysis and Problems

Equation (1) illustrates a three-step mechanism of inactivation with a dissociation reaction of the native form of the enzyme, N, into a denatured form, D, and with two parallel association reactions of the native and denatured forms into irreversibly denatured enzyme forms I₁ and I₂, respectively: N⇄^𝑘+1_𝑘−12D,2D−−→^𝑘2I₁,2N−−→^𝑘3I₂, (1) where 𝑘₊₁, 𝑘₋₁, 𝑘₂, and 𝑘₃ represent the rate constants of individual reactions. The material balances equations for N and D and temperature are given as follows [1]: 𝑑𝑐_N𝑑𝑡=−𝑘₊₁𝑐_N+𝑘^′₋₁𝑐²_D−2𝑘^′₃𝑐²_N,𝑑𝑐_D𝑑𝑡=2𝑘₊₁𝑐_N−2𝑘^′₋₁+𝑘^′₂𝑐²_D,𝑑𝑇𝑇𝑑𝑡=𝐾_𝐵−𝑇, (2) where 𝑐_N and 𝑐_D are the molar concentrations of the native and denatured enzyme form and 𝑇 is the temperature. The initial boundary condition is 𝑐_N=1,𝑐_D=0,𝑇=30+𝑇_𝐵at𝑡=0. (3)

3. Concentrations of Native Enzyme Form, Denatured Enzyme Form, and Temperature under Nonsteady State Condition

A powerful, easy-to-use analytic tool for nonlinear problems in general, namely, the homotopy analysis method. This analytic solution agrees well with numerical results and can be regarded as a definition of the solution of the considered nonlinear problem [13]. The homotopy analysis method (HAM) provides an analytical solution in terms of an infinite power series. In order to investigate the accuracy of the solution with a finite number of terms, the system of differential equations was solved. This method is a good technique when comparing to perturbation method and the numerical simulation. Homotopy analysis method has many advantages since it doesn’t rely on a special small parameter [14]. The results reveal that the method is very effective and simple. The method contains the auxiliary parameter ℏ, which provides us with a simple way to adjust and control the convergence region of solution series [15, 16]. Homotopy analysis method provides a powerful tool to strongly analyze linear and nonlinear problems [17, 18]. In this paper, homotopy analysis method (see Appendix A) is used to solve nonlinear differential equations. The analytical expression of molar concentrations of N, D, and temperature (see Appendices A and B) are 𝑐_N(𝑡)=−4ℎ𝑘^′₋₁𝑘₊₁+8ℎ𝑘^′₋₁𝑡𝑒^{−𝑘+1𝑡}+1+2ℎ𝑘^′₃𝑘₊₁𝑒^{−𝑘+1𝑡}+4ℎ𝑘^′₋₁𝑘₊₁−2ℎ𝑘^′₃𝑘₊₁𝑒^{−2𝑘+1𝑡},𝑐_D(𝑡)=2+2ℎ𝑘^′₃𝑘₊₁+20ℎ𝑘^′₋₁𝑘₊₁−12ℎ𝑘^′₋₁+𝑘^′₂𝑘₊₁−16ℎ𝑘^′₋₁𝑒^{−𝑘+1𝑡}𝑡+16ℎ𝑘^′₋₁+𝑘^′₂𝑘₊₁−4ℎ𝑘^′₃𝑘₊₁−16ℎ𝑘^′₋₁𝑘₊₁−2𝑒^{−𝑘+1𝑡}+2ℎ𝑘^′₃𝑘₊₁−4ℎ𝑘^′₋₁𝑘₊₁−4ℎ𝑘^′₋₁+𝑘^′₂𝑘₊₁𝑒^{−2𝑘+1𝑡},𝑇(𝑡)=𝑇_𝐵+30𝑒^−𝐾𝑡. (4)(5)(6) These equations represent the new analytical expression of the molar concentrations of the native enzyme, denatured enzyme, and temperature for various values of the parameters 𝑘₊₁, 𝑘^′₋₁, 𝑘^′₃, 𝑘^′₂, ℎ, 𝑇_𝐵, and 𝐾.

4. Numerical Simulation

To show the efficiency of the present method, our nonsteady-state result is compared with numerical solution. Equations (4) and (5) are also solved by numerical method using Matlab/Scilab program. The Scilab/Matlab program [19] is also given in Algorithm 1. The numerical values of parameters employed in Illeova et al. [1] and in this study are given in Table 1. The numerical solution is compared with our analytical result in Figures 1(a)–1(e).

Table 1 

Range of parameters used in Illeova et al. [1] and this work.

function main1 options= odeset (‘RelTol’,1e-6,‘Stats’,‘on’); %initial conditions T=100; Xo = [1; 0; T+30]; tspan = [0,10]; tic [t,X] = ode45(@TestFunction,tspan,Xo,options); toc figure hold on plot(t, X(:,1)) plot(t, X(:,2)) plot(t, X(:,3),‘.’) legend(‘x1’,‘x2’,‘x3’) ylabel(‘x’) xlabel(‘t’) return function [dx_dt]= TestFunction(t,x) k1=2;k2=0.1;k3=10∧(-3);k4=2.66*10∧(-3);K=4.44*10∧(-2);T=100; dx_dt(1) =-k1*x(1)+k2*x(2)∧2-2*k3*x(1)∧2; dx_dt(2) =2*k1*x(1)-2*(k2+k4)*x(2)∧2; dx_dt(3)=K*(T-x(3)); dx_dt = dx_dt’; return

Algorithm 1 

Scilab/Matlab program to find the numerical solution of (2) and (3).

(a)

(b)

(c)

(d)

(e)

Figure 1 

(a–d) Molar concentrations of N and D versus time 𝑡 for various values of the parameters 𝑘+1, 𝑘′−1, 𝑘′3, and 𝑘′2. Solid lines represent the analytical solution presented in this work ((4) and (5)) and dotted lines represent the numerical solution. (e) Temperature 𝑇 versus time 𝑡 for various values of 𝑇𝐵 when 𝐾=4.44×10−2 s−1.

5. Discussions

Equations (4) and (5) represent the closed form of an analytical expression of concentrations of native and denatured enzyme for all values of parameters. In Figures 1(a)–1(e), the analytical results are compared with numerical results for various values of parameters and satisfactory agreement is noted. Figures 1(a)–1(d) represent the molar concentrations for various values of rate constant. From the figure, it is inferred that, when time increases the concentration of 𝑐_N decreases whereas the concentration of 𝑐_D increases. The time taken to reach the maximum value of 𝑐_D is the same as the time taken to reach the minimum value of 𝑐_N. The steady value of 𝑐_D and 𝑐_N depends upon the rate constants.

Figures 2(a)–2(c) represent the molar concentration of N versus time 𝑡 for various values of rate constant 𝑘₊₁ and 𝑘^′₋₁. From that figure, it is observed that, 𝑐_N increases when 𝑘₊₁ and 𝑘^′₋₁ increase. Figure 3 represents the molar concentration of D versus time 𝑡 for various values of parameter 𝑘₊₁. From this figure, it is found that, the value of concentration 𝑐_D initially increases and reaches the steady state value when 𝑡≥5. The concentration 𝑐_D increases when 𝑘₊₁ increases. The concentration becomes zero when 𝑘₊₁≤0.01 s⁻¹. Figure 1(e) denotes the representation temperature 𝑇 versus time 𝑡 for various values of 𝑇_𝐵. The value of temperature increases when 𝑇_𝐵 increases. It is almost linear with respect to time 𝑡.

(a)

(b)

(c)

Figure 2 

(a) Molar concentration of N versus time 𝑡 using (4) for various values of 𝑘+1 when 𝑘′−1=0.01 s−1, 𝑘′3=10−3 s−1, and ℎ=−1. (b) Molar concentration of N versus time 𝑡 using (4) for various values of 𝑘′−1 when 𝑘+1=8.80×10−1 s−1, 𝑘′3=0.28×10−3 s−1, and ℎ=−1. (c) Molar concentration of N versus time 𝑡 using (4) for various values of 𝑘+1 when 𝑘′−1=8.64×101 s−1, 𝑘′3=0.28×10−3 s−1, and ℎ=−1.

Figure 3 

Molar concentration of D versus time 𝑡 using (5) for various values of 𝑘+1 when 𝑘′−1=0.1 s−1, 𝑘′2=2.66×10−4 s−1, 𝑘′3=10−3 s−1, and ℎ=−1.

6. Conclusion

In this work, we obtained the analytical expression of concentrations in terms of rate constants 𝑘₊₁, 𝑘^′₋₁, 𝑘^′₃, 𝑘^′₂, bath temperature 𝑇_𝐵, and coefficient in the enthalpy balance 𝐾. The nonlinear ordinary differential equations have been solved analytically. The closed analytical expressions of molar concentrations of 𝑐_N, 𝑐_D, and temperature 𝑇 are obtained using the homotopy analysis method. An agreement with the numerical result is noted. The information gained from this theoretical model can be useful for the kinetic analysis of the experimental results over handling rate constants and molar concentrations.

Appendices

A. Basic Concept of Homotopy Analysis Method

Consider the following differential equation [20–22]: 𝑁[𝑢(𝑡)]=0, (A.1) where 𝑁 is a nonlinear operator, 𝑡 denotes an independent variable, and 𝑢(𝑡) is an unknown function. For simplicity, we ignore all boundary or initial conditions, which can be treated in a similar way. By means of generalizing the conventional homotopy method, Liao [19] constructed the so-called zero-order deformation equation as follows: 1−𝑝𝐿𝜑𝑡;𝑝−𝑢₀(𝑡)=𝑝ℎ𝐻(𝑡)𝑁𝜑𝑡;𝑝, (A.2) where 𝑝∈[0,1] is the embedding parameter, ℎ≠0 is a nonzero auxiliary parameter, 𝐻(𝑡)≠0 is an auxiliary function, 𝐿 is an auxiliary linear operator, 𝑢₀(𝑡) is an initial guess of 𝑢(𝑡), and 𝜑(𝑡∶𝑝) is an unknown function. It is important that one has great freedom to choose auxiliary unknowns in HAM. Obviously, when 𝑝=0 and 𝑝=1, it holds the following: 𝜑(𝑡;0)=𝑢₀(𝑡),𝜑(𝑡;1)=𝑢(𝑡), (A.3) respectively. Thus, as 𝑝 increases from 0 to 1, the solution 𝜑(𝑡;𝑝) varies from the initial guess 𝑢₀(𝑡) to the solution 𝑢(𝑡). Expanding 𝜑(𝑡;𝑝) in Taylor series with respect to 𝑝, we have 𝜑𝑡;𝑝=𝑢₀(𝑡)+^+∞_𝑚=1𝑢_𝑚(𝑡)𝑝^𝑚, (A.4) where 𝑢_𝑚(𝑡)=1𝜕𝑚!^𝑚𝜑𝑡;𝑝𝜕𝑝^𝑚_𝑝=0. (A.5) If the auxiliary linear operator, the initial guess, the auxiliary parameter ℎ, and the auxiliary function are so properly chosen and the series (A.4) converges at 𝑝=1, then we have 𝑢(𝑡)=𝑢₀(𝑡)+^+∞_𝑚=1𝑢_𝑚(𝑡). (A.6) Differentiating (A.2) for 𝑚 times with respect to the embedding parameter 𝑝, and then setting 𝑝=0 and finally dividing them by 𝑚!, we will have the so-called 𝑚th-order deformation equation as follows: 𝐿𝑢_𝑚−𝜒_𝑚𝑢_𝑚−1=ℎ𝐻(𝑡)ℜ_𝑚𝑢_𝑚−1, (A.7) where ℜ_𝑚𝑢_𝑚−1=1()!𝜕𝑚−1^𝑚−1𝑁𝜑𝑡;𝑝𝜕𝑝^𝑚−1,𝜒_𝑚=0,𝑚≤1,1,𝑚>1. (A.8) Applying 𝐿⁻¹ on both sides of (A.7), we get 𝑢_𝑚(𝑡)=𝜒_𝑚𝑢_𝑚−1(𝑡)+ℎ𝐿⁻¹𝐻(𝑡)ℜ_𝑚𝑢_𝑚−1. (A.9) In this way, it is easily to obtain 𝑢_𝑚 for 𝑚≥1, at 𝑀th order, we have 𝑢(𝑡)=^𝑀_𝑚=0𝑢_𝑚(𝑡). (A.10) When 𝑀→+∞, we get an accurate approximation of the original equation (A.1). For the convergence of the above method, we refer the reader to Liao [13]. If (A.1) admits unique solution, then this method will produce the unique solution. If (A.1) does not possess unique solution, the HAM will give a solution among many other (possible) solutions.

B. Approximate Analytical Solutions of the System of (2) Using Homotopy Analysis Method

In this appendix, we indicate how (4) and (5) in this paper are derived. The homotopy analysis method was constructed to determine the solution of (2) 𝑑𝑐_N𝑑𝑡=−𝑘₊₁𝑐_N+𝑘^′₋₁𝑐²_D−2𝑘^′₃𝑐²_N,𝑑𝑐_D𝑑𝑡=2𝑘₊₁𝑐_N−2𝑘^′₋₁𝑐²_D−2𝑘^′₂𝑐²_D. (B.1) In order to solve (B.1) by means of the HAM, we first construct the zeroth-order deformation equation by taking 𝐻(𝑡)=1, 1−𝑝𝑑𝜑𝑡;𝑝𝑑𝑡+𝑘₊₁𝜙𝑡;𝑝=𝑝ℎ𝑑𝜑𝑡;𝑝𝑑𝑡+𝑘₊₁𝜙𝑡;𝑝−𝑘^′₋₁𝜓²𝑡;𝑝+2𝑘^′₃𝜑²𝑡;𝑝,1−𝑝𝑑𝜓𝑡;𝑝𝑑𝑡−2𝑘₊₁𝜙𝑡;𝑝=𝑝ℎ𝑑𝜓𝑡;𝑝𝑑𝑡−2𝑘₊₁𝜙𝑡;𝑝+2𝑘^′₋₁+𝑘^′₂𝜓²𝑡;𝑝, (B.2) subject to the following initial conditions: 𝜑₀0;𝑝=1,𝜓₀𝜑0;𝑝=0,_𝑖0;𝑝=0,𝜓_𝑖0;𝑝=0,∀𝑖=1,2,…, (B.3)(B.4) where 𝑝∈[0,1] is an embedding parameter and ℎ≠0 is the so-called convergence control parameter. When 𝑝=0𝑑𝜑(𝑡;0)𝑑𝑡+𝑘₊₁𝜙()()𝑡;0=0,𝑑𝜓𝑡;0𝑑𝑡−2𝑘₊₁𝜙(𝑡;0)=0. (B.5) From (B.5) and initial condition (B.3), we get 𝜑(𝑡;0)=𝑒^{−𝑘+1𝑡},𝜓(𝑡;0)=21−𝑒^{−𝑘+1𝑡}. (B.6)(B.7) When 𝑝=1, (B.2) are equivalent to (B.1), thus it holds the following: 𝜑(𝑡;1)=𝑐_N(𝑡),𝜓(𝑡;1)=𝑐_D(𝑡). (B.8) Expanding 𝜑(𝑡;𝑝) and 𝜓(𝑡;𝑝) in Taylor series with respect to the embedding parameter 𝑝, we have 𝜑𝑡;𝑝=𝑐_N(𝑡)+^∞_𝑚=1𝑐_N𝑚(𝑡)𝑝^𝑚,𝜓𝑡;𝑝=𝑐_D(𝑡)+^∞_𝑚=1𝑐_D𝑚(𝑡)𝑝^𝑚, (B.9) where 𝑐_N0(𝑡)=𝜑(𝑡;0), (B.10) and 𝑐_N𝑚(𝑡)[𝑚=1,2,…] will be determined later. Note that the above series contains the convergence control parameter ℎ. Assuming that ℎ is chosen so properly that the above series is convergent at 𝑝=1, we have the solution series as follows: 𝑐_N(𝑡)=𝑐_N0(𝑡)+^+∞_𝑚=1𝑐_N𝑚(𝑡)𝑝^𝑚,𝑐_D(𝑡)=𝑐_D0(𝑡)+^+∞_𝑚=1𝑐_D𝑚(𝑡)𝑝^𝑚, (B.11) where 𝑐_N𝑚(𝑡)=1𝜕𝑚!^𝑚𝑐_N𝑡;𝑝𝜕𝑝^𝑚_𝑝=0,𝑐_D𝑚(𝑡)=1𝜕𝑚!^𝑚𝑐_D𝑡;𝑝𝜕𝑝^𝑚_𝑝=0. (B.12) Equating the like coefficients of the like powers of 𝑝 we have 𝑝¹∶𝑑𝜑₁𝑑𝑡+𝑘₊₁𝜑₁=(ℎ+1)𝑑𝜑₀𝑑𝑡+𝑘₊₁𝜑₀−ℎ𝑘^′₋₁𝜓²₀+2ℎ𝑘^′₃𝜑²₀,𝑝²∶𝑑𝜑₂𝑑𝑡+𝑘₊₁𝜑₂=(ℎ+1)𝑑𝜑₁𝑑𝑡+𝑘₊₁𝜑₁−2ℎ𝑘^′₋₁𝜓₀𝜓₁+4ℎ𝑘^′₃𝜑₀𝜑₁,𝑝¹∶𝑑𝜓₁𝑑𝑡−2𝑘₊₁𝜑₁=(ℎ+1)𝑑𝜓₀𝑑𝑡−2𝑘₊₁𝜑₀+2ℎ𝑘^′₋₁+𝑘^′₂𝜓²₀,𝑝²∶𝑑𝜓₂𝑑𝑡−2𝑘₊₁𝜑₂=(ℎ+1)𝑑𝜓₁𝑑𝑡−2𝑘₊₁𝜑₁+4ℎ𝑘^′₋₁+𝑘^′₂𝜓₀𝜓₁, (B.13) and so on. Solving (B.13) using the initial condition (B.3) and (B.4), we obtain the following result: 𝜑₁(𝑡)=2ℎ𝑘^′₃𝑘₊₁𝑒^{−𝑘+1𝑡}−4ℎ𝑘^′₋₁𝑘₊₁+4ℎ𝑘^′₋₁𝑘₊₁𝑒^{−2𝑘+1𝑡}+8ℎ𝑘^′₋₁𝑡𝑒^{−𝑘+1𝑡}−2ℎ𝑘^′₃𝑘₊₁𝑒^{−2𝑘+1𝑡},𝜓₁(𝑡)=2ℎ𝑘^′₃𝑘₊₁+20ℎ𝑘^′₋₁𝑘₊₁−12ℎ𝑘^′₋₁+𝑘^′₂𝑘₊₁−4ℎ𝑘^′₃𝑘₊₁𝑒^{−𝑘+1𝑡}−4ℎ𝑘^′₋₁𝑘₊₁𝑒^{−2𝑘+1𝑡}−16ℎ𝑘^′₋₁𝑒^{−𝑘+1𝑡}1𝑡+𝑘₊₁+2ℎ𝑘^′₃𝑘₊₁𝑒^{−2𝑘+1𝑡}−4ℎ𝑘^′₋₁+𝑘^′₂𝑘₊₁𝑒^{−2𝑘+1𝑡}+16ℎ𝑘^′₋₁+𝑘^′₂𝑘₊₁𝑒^{−𝑘+1𝑡}. (B.14)(B.15) Adding (B.6), (B.14), (B.7), and (B.15), we get (4) and (5) in the text. Here we can obtain the solution of the equation using only two iterations which leads to the high accuracy.

Nomenclature

𝑐N:	Molar concentration of the native enzyme form (mole/cm3)
𝑐D:	Molar concentration of the denatured enzyme form (mole/cm3)
𝑘+1, 𝑘−1, 𝑘2, 𝑘3:	Rate constants of individual reaction (s−1)
𝑘′−1, 𝑘′2, 𝑘′3:	Modified rate constants (s−1)
𝐾:	Coefficient in the enthalpy balance (s−1)
𝑇𝐵:	Bath temperature (K)
𝑇:	Temperature (K)
𝑡:	Time (s).

Acknowledgments

This work was supported by the Council of Scientific and Industrial Research (CSIR no. 01(2442)/10/EMR-II), the Government of India. The authors are thankful to Dr. R. Murali, the Principal, The Madura College, Madurai, and Mr. S. Natanagopal, the Secretary, Madura College Board, Madurai for their encouragement. It is a pleasure to thank the referees for their valuable suggestions.