Numerical Solution to Nonlinear Biochemical Reaction Model Using Hybrid Polynomial Basis Differential Evolution Technique

In this paper, we present an approximate numerical solution to the well known Michaelis-Menten nonlinear biochemical reaction system using a stochastic technique based on hybrid polynomial basis evolutionary computing. The approximate solution is expanded as a linear combination of polynomial basis with unknown parameters. The system of nonlinear differential equation is transformed into an equivalent global error minimization problem. A trial solution is formulated using a fitness function with unknown parameters. Two popular evolutionary algorithms such as Genetic algorithm (GA) and Differential evolution (DE) are used to solve the minimization problem and to obtain the unknown parameters. The effectiveness of the proposed technique is demonstrated in contrast with fourth-order Runge Kutta method (RK-4) and some well known standard methods including homotopy perturbation method (HPM), variational iteration method (VIM), differential transform method (DTM), and modified Picard iteration method (Picard-Padé). The comparisons of numerical results validate the efficacy and viability of the suggested technique. The results are found to be in sharp agreement with RK-4 compared to some popular standard methods.


Introduction
A wide spectrum of scientific phenomena appearing in real life problems are governed by nonlinear ordinary differential equations (ODEs).The analytical treatment of such nonlinear ODEs is of much importance to get the insight of the systems behavior.Since most of nonlinear ODEs either do not have an exact solution or obtaining the same analytically is difficult, these problems must be tackled using approximate analytical and numerical methods.A great number of approximate analytical and numerical methods including homotopy perturbation method (HPM), variational iteration method (VIM), adomian decomposition method (ADM), differential transform method (DTM) etc. have been utilized for solving nonlinear ODEs [1][2][3][4][5].
In recent years many authors have used nature and biologically inspired based computation techniques as an alternative for solving nonlinear ODEs arising in diverse fields of engineering and science [6][7][8][9][10][11][12][13][14][15][16].Very recently Malik et al. [6,7] employed heuristic technique based on hybrid genetic algorithm for numerically solving the nonlinear singular boundary value problems in physiology and the Bratu problem.Arqub et al. [8] used genetic algorithm (GA) based method for solving linear and nonlinear singular boundary value problems (BVPs).Caetano et al. [9] used genetic algorithm (GA) based neural network (NN) for the solution Numerical solution to nonlinear biochemical reaction model 101 of nonlinear ODEs arising in atomic and molecular physics.Khan et al. [10] used particle swarm optimization (PSO) based NN technique for solving nonlinear ODEs including Wessinger's equation.Behrang et al. [11] employed particle swarm optimization (PSO) based neural network for solving the nonlinear ODE of Darcian fluid of vertical cone implanted in porous media.In this paper, we consider the well-known Michaelis-Menten biochemical reaction model [1][2][3][4][5] ⇌ → where is the enzyme, the substrate, the intermediate complex and the product.
The time evolution of scheme (1) can be determined from the solution of the system of coupled nonlinear ODEs [1][2][3][4][5] ( subject to the initial conditions where the parameters , and are positive rate constants for each reaction.Systems (2) -( 5) can be reduced to only two equations for and and in dimensionless form of concentrations of substrate, , and intermediate complex between enzyme and substrate, , [1][2][3][4][5] (7) where , , and are dimensionless parameters.The reader may refer [1] and references therein for a detailed mathematical formulation of equations ( 7) and (8).Many researchers have investigated the numerical solution of the system of coupled nonlinear ODEs (7) and ( 8) due to its practical importance to biochemists.Various different methods including homotopy perturbation method (HPM) [2], variational iteration method (VIM) [3], differential transform method (DTM) [4], and modified Picard iteration method (Picard-Padé) [5] have been utilized for numerically solving (7) and (8).
The main purpose of this paper is to obtain the numerical solution of the biochemical reaction model ( 7) and ( 8) using a novel technique based on nature inspired computation.The method employs two evolutionary algorithms including genetic algorithm (GA) and differential evolution (DE) for the minimization of a problem specific fitness function.To the best of our knowledge no body yet has solved the system of ( 7) and ( 8) using the approach proposed in this paper.Moreover the proposed technique is stochastic in nature as compared to the standard methods [1][2][3][4][5] utilized for solving (7) and (8), which are deterministic in nature.
To prove the efficiency and the viability of the presented technique comparisons of our numerical results are made with fourth-order Runge Kutta method (RK-4) and some well known deterministic standard techniques including HPM [2], VIM [3], DTM [4] and modified Picard-Padé [5].
The rest of the paper is organized as follows: In section 2, we provide a brief description of the proposed method, and evolutionary algorithms such as genetic algorithm (GA) and differential evolution (DE) are also introduced.In section 3, we present numerical results followed by the discussion on our findings.Finally in section 4 we give some concluding remarks.

Descript of Proposed Method for Biochemical Reaction Model
To obtain the numerical solution of biochemical reaction model ( 7) and (8), we may assume the approximate solution of , and y and their first derivates , are a linear combinations of some polynomial basis functions given as follows.
( ) in ( 10) -( 13) are determined by formulating a trial solution of the given problem using a fitness function ( ) given by (14) where j is the cycle index, and represent the mean square errors associated with the system of coupled nonlinear ODEs (7) and ( 8) respectively, given as follows The FF function given by ( 14) contains the unknown adaptable coefficients ( , .The minimization of fitness function ( is carried out by employing the evolutionary algorithms GA and DE for acquiring the unknown parameters.The optimal values of the unknown adaptable parameters corresponding to the minimum fitness function ( ) are achieved.Consequently the approximate numerical solutions and of the biochemical reaction model are obtained using (10) and (11).

Stochastic Global Search Evolutionary Algorithms: GA and DE
Evolutionary algorithms (EAs) are population based heuristic search methods that use the nature and biologically inspired process in an iterative manner to reach an optimal solution [17].
Genetic algorithm (GA) [17] is one of the well-known techniques in EAs that finds the optimal solution of a problem from a randomly generated population of individuals called chromosome.Each individual within a population is regarded as a possible solution to the problem.The individuals within a population are evaluated using a fitness function that is specific to the problem at hand.The algorithm evolves population iteratively by means of three primary operations: selection, crossover, and mutation to reach the optimal solution [17].The procedural steps of GA are given in algorithm 1 while the parameters settings for its implantation used in this study are given in Table 1.

Algorithm 1: Genetic Algorithm (GA)
Step 1: (Population Initialization) A population of N individuals or chromosomes (C1, C2, ….,CN) each of length M is generated using random number generator.The length of each chromosome represents the number of unknown coefficients.

Population creation function Step 2: (Fitness Evaluation)
A fitness function is used to compute the fitness of each chromosome.

Step 3: (Selection and Reproduction)
The chromosomes from the current population are chosen on the basis of their fitness value which acts as parents for new generation.These parents produce children (offsprings) with a probability to their fitness through crossover operation.
Step 4: Mutation Mutation operation introduces random alterations in the genes to maintain the genetic diversity to find a good solution.

Step 5: (Stoppage Criteria)
The algorithm terminates if a certain level of fitness has reached or the maximum number of cycles (generations) has exceeded.Else repeat steps 2 to 4.  [18][19] is another popular and powerful global search evolutionary algorithm (EA) that evolves populations to reach an optimal solution using the operations crossover, mutation, and selection like GA. Contrary to GA, the algorithm uses mutation operation as search mechanism, in which the weighted difference between two population vectors is added to create a new parameter vector.The mutated vectors are randomly mixed with the last target vector to generate a trial vector, this operation is called crossover.If the trial vector gives a smaller value of the objective function than the target vector, the trial vector replaces the target vector in the current population.The algorithm iteratively evolves populations by applying the aforementioned steps until the optimal results are achieved or the stopping criterion is reached.The steps of algorithm are given in algorithm 2 while the values of the parameters used for implementation in this work are given in Table 2.The algorithm stops if the number of generations or a desired fitness reaches, else the algorithm goes to step 2.

Results and Discussion
In this section proposed methodology is applied to the biochemical reaction model ( 7) -( 8) to obtain the numerical solutions in the interval 0 1 with the initial conditions, 0 1 and 0 0 , and the values of dimensionless parameters are taken 0.3752, 1.0 , 0.1, for a direct comparison with other classical methods used in [2-5] .
To apply the proposed method to biochemical reaction model we construct a trial solution in the form of a fitness function ( ) given by ( 17)- (19).
where x , , and their first derivates , , are given by ( 10) -( 13), and m (no. of basis functions) has been taken 7. Therefore the number of unknown parameters ( , , … , , , … that need to be tailored is 14.The values of these unknown parameters are obtained using GA and DE by performing the minimization of the fitness function ( ) given by (19).For implementation Matlab has been utilized.
The parameter values of the evolutionary algorithms GA and DE used for the execution are given in Table 2. To study the effectiveness of the proposed method, the approximate solution of the biochemical reaction model ( 7) and ( 8) was obtained employing three different number of basis functions (m = 3, 5, and 7).The optimal values of the unknown parameters corresponding to the minimum fitness achieved by GA and DE are provided in Table 3  The approximate numerical solutions and of the biochemical reaction model ( 7) and ( 8) are obtained by substituting the optimal values of unknown parameters from Tables 3 and Table 4 in (10) and (11).The numerical solutions for m = 7 obtained by the proposed method are presented in Table 5.In Table 6 and Table 7, the results of proposed method are compared against the forth-order Runge-Kutta (RK-4) and some classical method including HPM [2], VIM [3], DTM [4], and PIM-Padé [5].Furthermore Table 8 shows a comparison of absolute errors yielded by the proposed method and classical methods [2][3][4][5], in contrast with RK-4.The comparison of results and absolute errors in Tables 6 -8 evidently reveals that the proposed method provides the approximate solution of biochemical reaction model with a high degree of accuracy and in a sharp agreement with RK-4 as compared to some well known classical methods such as HPM [2], VIM [3], and DTM [4].Moreover it is observed from the comparison that the results obtained by the proposed method are fairly comparable to PIM-Padé [5/5] used in [5].[2], VIM [3], DTM [4], and Picard-Padé [5].The influence of different number of basis functions (i.e.change in m) on the performance of the proposed scheme and convergence of the evolutionary algorithms GA and DE are analyzed next.We used m = 3, 5 for evaluating the performance, therefore the number of unknown parameters to be tailored are 6 and 10 respectively.Without any other change in the parameter values except the chromosome size 6 and 10 for m= 3 and 7 respectively, GA and DE are executed to acquire the unknown parameters.The optimal values of the unknown parameters acquired by GA and DE are given in Tables 3 and 4 respectively.In Fig. 1, 2 we provide the approximate solutions x(t) and y(t) obtained using the proposed method with m = 3,5, also solutions using RK-4 and with m=7 by the proposed method are shown for the comparison.It is evident from Fig. 1, 2 that the accuracy of the results improves with an increase in m i.e number of basis functions.In Table 9 we also give a comparison of average absolute errors yielded by the proposed method for m = 3, 5, 7 and the number of generations (number of cycles) taken for achieving the desired minimum fitness by GA and DE.From Table 9 it could be concluded that the proposed method yields improved performance with increase in m but at the cost of large number of generations consequently high computational cost.Nonetheless the proposed method provides the solution of the biochemical reaction model with better accuracy even with m = 3 as compared to classical methods HPM [2], VIM [3], and DTM [4], which proves the effectiveness of our method.

Conclusions
A simple yet an efficient stochastic heuristic method based on hybrid approach of polynomial basis functions and evolutionary computation has been presented for numerically solving system of nonlinear ordinary coupled differential equations.The effectiveness of the proposed method has been demonstrated by numerically solving the nonlinear biochemical reaction system.On the basis of the numerical results and comparisons made with the fourth-order Runge Kutta method (RK-4) and some popular classical methods, it can be concluded that the proposed method is efficient and viable for solving the nonlinear biochemical reaction model.The findings have revealed that the proposed method shows significant supremacy on classical methods HPM, VIM, and DTM.Moreover the proposed method can give the approximate solution of the given problem with convenience and on continuous grid of time once the unknown parameters have been acquired.

( 13 )
Numerical solution to nonlinear biochemical reaction model103where are real valued unknown parameters, and m is the number of basis functions.The values of unknown parameters (

Algorithm 2 : 1
Differential Evolution(DE) Step 1: (Population Initialization) A population of N chromosomes , , … ., is generated randomly.Each chromosome has D number of genes representing the number of unknown coefficients.Step 2: (Updating): ∀ 1 Update the target vector, using the following operations.a) Mutation: choose 3 members randomly between 1 and N, all three different among themselves and also different from to generate a mutant vector as follows b) Crossover: The trial vector is formed as follows and 1 c) Selection: Fitness of parents and mutants are evaluated for next generation.Parents and mutants are sorted according to their fitness values.values of parent , and mutant , .Step 3: (Stopping Criterion)

Figure1: 4 Figure 2 :
Figure1: Effect of change in m on approximate solution x(t) and comparison with RK-4

Table 1 .
Parameter settings for GA

Table 2 .
and Table 4 respectively.Parameter values of GA and DE

Numerical solution to nonlinear biochemical reaction model 107Table 4 .
Optimal

Table 7 .
Comparison of Numerical Results for y(t) between the proposed method, RK-4, and Classical methods HPM

Numerical solution to nonlinear biochemical reaction model 109Table 8 .
Comparison of absolute errors between the proposed method and the classical methods.

Table 9 .
Comparison of average absolute errors and convergence speed for different values of m