Kinetic parameters estimation of protease production using penalty function method with hybrid genetic algorithm and particle swarm optimization

ABSTRACT

Almost all optimization techniques are restricted by the problems' dimensions and large search spaces. This research focuses on a special hybrid method combining two meta-heuristic techniques, genetic algorithms (GA) and particle swarm optimization (PSO) that aims at overcoming this issue. This method investigates the potential impact of constraints (or the feasible regions) on the search pattern of GA and PSO. The proposed algorithm was applied for parameter estimation of batch fermentation process for alkaline protease production by Bacillus licheniformis in submerged culture. Furthermore, a comparison of proposed hybrid GA/PSO with pure GA and pure PSO was carried out. The results revealed that combination of these two meta-heuristic algorithms speeds up the search (about two-fold faster) in comparison to pure algorithms, since it benefits from synergy. Hence, the proposed method can be considered as an applicable method for parameter estimation of biological models in particular for large search space problems. Also, it was concluded that PSO has a slightly better performance and possesses better convergence and computational time than GA.

KEYWORDS:

• Kinetic parameters
• parameter estimation
• protease
• genetic algorithm
• particle swarm optimization
• penalty function

Introduction

Accurate estimation of the kinetic parameters of enzymes is an important part in the analysis of biochemical problems.[1] Meta-heuristics are generic methods which offer good solutions, even global optima, within reasonable computing time. Thus, the use of meta-heuristics has received more and more attention.[2] Different meta-heuristic methods have been applied to surmount the parameter estimation difficulties. Genetic algorithms (GAs), originally developed by Holland in 1975,[3] are quite promising as a stochastic global optimization method.[4] GA belongs to the larger class of evolutionary algorithms, which use an approach inspired by natural evolution, e.g. mutation, selection and crossover, to solve optimization problems.[5] GA is identical to biological evolution in the principles of natural selection and the survival of the fittest.[6] It is perhaps the most popular method used in parameter estimation problems.[7] Particle swarm optimization (PSO), proposed by Kennedy and Eberhart in 1995, [8] is also a meta-heuristic approach inspired by social models and swarming theory to find good solutions to optimization problems.[9]

These optimization approaches can find global optima more quickly through cooperation and competition among the population of potential solutions of the search space even for complex optimization problems such as fermentation processes.[10]

In recent years, both the GA and the PSO, as well as hybrid algorithms combining the two, have been applied to parameter estimation problems in the context of kinetic models of microorganism growth.[7,10–17]

GA has an advantage in exploration search and PSO has an advantage in sharing movement information between particles.[18] Therefore in order to keep the advantages of both algorithms, combining between GA and PSO seems to be a logical method.

Also, many optimization problems are constrained with different conditions. Due to the existence of constraints, the optimization problems are more difficult than the unconstrained ones. Therefore, research on how to make the objective function find the optimal solution in the feasible region has great significance.[19]

Penalty methods were developed in order to eliminate some or all of the constraints and add to the objective function a penalty term which prescribes a high cost to infeasible points.[20] Therefore, in this work, growth kinetic parameter estimation has been carried out by a hybrid GA/PSO algorithm which investigates the potential impact of constraints (or the feasible regions) on the search pattern of GA and PSO. Finally, the result of the proposed hybrid algorithm has been compared to pure algorithms.

Materials and methods

Micro-organisms

The bacterial strain used in this study was isolated from soils collected at different locations in north of Iran. It was identified as Bacillus licheniformis (strain MG5), as described previously.[21] The strain was maintained on nutrient agar slant and used throughout the study.

Fermentation conditions

The pre-culture medium was nutrient broth containing 2 g/L yeast extract, 5 g/L peptone, 5 g/L NaCl and 1 g/L beef extract, sterilized at 121 °C for 15 min.[22]

Production of protease was carried out in 500-mL shake-flasks with 100 mL medium containing the following constituents: 20 g/L glucose, 7 g/L K₂HPO₄, 2 g/L KH₂PO₄, 0.2 g/L MgSO₄.7H2O, 0.5 g/L NaCl, 1 g/L (NH₄)₂SO₄ and 0.07 g/L CaCl₂. Separately sterilized glucose was added into the medium as a carbon source just before inoculation. The above-mentioned pre-cultures were inoculated and the submerged cultures were incubated at 32 °C and 150 r/min in a shaker incubator.[11]

Analytical procedures

Periodically, samples were withdrawn and the biomass concentration, residual glucose concentration and protease enzyme activity were determined. Bacterial growth was monitored by using a spectrophotometer at an optical density of 600 nm (OD₆₀₀), which was then converted to cell dry weight (CDW) using a calibration curve. A standard graph was plotted for CDW versus absorbance for further estimation of CDW.[22] Dinitrosallicylic acid was used for residual sugar assessment.[23] Protease activity was determined as released tyrosine from the supernatants according to a modified Lowry method.[24]

All experiments were conducted in triplicate and the experimental data given in this paper are arithmetic mean values of the triplicate experiments.

Problem formulation

Batch fermentation refers to a partially closed system in which most of the required components for micro-organism growth and maintenance are loaded before the process starts.[15] The mathematical model of the batch fermentation process is presented by the following nonlinear differential equation system [25]: (1) $\frac{dX(t)}{dt}={\mu }_{\max }\left[1-\frac{X(t)}{X_{\max }}\right]X(t)$(1) (2) $\frac{dP(t)}{dt}=Y_{P/X}\frac{dX(t)}{dt}+kX(t)$(2) (3) $-\frac{dS(t)}{dt}=\frac{1}{Y_{X/S}}\frac{dX(t)}{dt}+m_{s}X(t)$(3)

In the above equations, X, S and P are the biomass concentration, substrate concentration and product concentration, respectively (g/L). X_max is the maximum biomass concentration (g/L). X₀, S₀ and P₀ are the initial concentration of biomass, substrate and product, respectively (g/L), t is the current fermentation time (h), μ_max is the maximum specific growth rate (h⁻¹), Y_P/X is the product yield in terms of biomass (gP/gX) and k is the secondary coefficient of product formation or destruction (gP/g.h), Y_X/S is the biomass yield in terms of substrate (gX/gS) and m_s is the maintenance coefficient (gS/gX.h).

Equation (1)(1) $\frac{dX(t)}{dt}={\mu }_{\max }\left[1-\frac{X(t)}{X_{\max }}\right]X(t)$(1) is logistic equation that describes the dynamics of cell mass concentration whose growth is dependent on the availability of substrate and also on the total cell population. Equation (2)(2) $\frac{dP(t)}{dt}=Y_{P/X}\frac{dX(t)}{dt}+kX(t)$(2) is the Luedeking–Piret equation that describes the product formation and Equation (3)(3) $-\frac{dS(t)}{dt}=\frac{1}{Y_{X/S}}\frac{dX(t)}{dt}+m_{s}X(t)$(3) is the modified Luedeking–Piret equation that describes the dynamics of substrate concentration.

Parameter estimation

In this study, the kinetic parameters estimated were {${\mu }_{\max },k,Y_{P/X},m_s,Y_{X/S}$ } as defined previously in the differential Equations (1)(1) $\frac{dX(t)}{dt}={\mu }_{\max }\left[1-\frac{X(t)}{X_{\max }}\right]X(t)$(1) –(3(3) $-\frac{dS(t)}{dt}=\frac{1}{Y_{X/S}}\frac{dX(t)}{dt}+m_{s}X(t)$(3) ). The following bounds were used for the five kinetic parameters (n = 5): $0.01\le {\mu }_{\max }\le 0.7,0.001\le k\le 0.5,0.1\le Y_{P/X}\le 1,0.01\le m_s\le 1,0.1\le Y_{X/S}\le 1$. Each bound is divided to m sections; m is the same for the all kinetic parameters (at first m was set to 2) (4) $x_d^l\le x_d\le x_d^{u}d=1,2,\mathrm{...},n(x_1={\mu }_{\max },x_2=k,x_3=Y_{P/X},x_4=m_s,x_5=Y_{X/S})$(4) (5) $x_d^l+(i-1)\frac{x_d^u-x_d^l}{m}\le x_{d_i}\le x_d^l+i\frac{x_d^u-x_d^l}{m}i=1,2,\mathrm{...},m$(5) where $x_d^{}$ is the dth parameter, $x_d^l$ and $x_d^u$ are the lower and the upper bounds of the dth parameter, respectively. Also $x_{d_i}^{}$ denote the dth parameter bounded to the ith sub-bound of the parameter bound.

In the proposed hybrid GA/PSO algorithm, the search space is divided to mⁿ subspaces (Equation (5)(5) $x_d^l+(i-1)\frac{x_d^u-x_d^l}{m}\le x_{d_i}\le x_d^l+i\frac{x_d^u-x_d^l}{m}i=1,2,\mathrm{...},m$(5) ). Each subspace consists of a vector of solutions named x {$x={(x_1{}_{{}_{i1}},x_2{}_{{}_{i2}},\mathrm{...},x_d{}_{{}_{id}})}^T$ , where each one of i1, i2, … and id are an integer in the range [1,10]}. Then a group of GAs $\left\{G{A}_1,G{A}_2,\mathrm{...},G{A}_{m^n}\right\}$ searches the best solution in each subspace. This means that a set of populations is generated by each GA in its own subspace randomly; each population introduces a set of potential solutions to the problem. The potential solution (individual) is coded as a binary vector, called a “chromosome”, the elements of which are called ‘genes’. In the GA algorithm, the population has chromosomes that represent candidate solutions; each chromosome is an n-dimensional real value vector where n is the number of optimized parameters. Therefore, each optimized parameter represents a dimension of the problem space.

One chromosome consists of five genes, coding for the parameters {$x_1{}_{{}_{i1}},x_2{}_{{}_{i2}},\mathrm{...},x_d{}_{{}_{id}}$ }, respectively. One gene is coded by a binary string of 20 bits. Therefore, every possible solution of the parameters is simply represented by a binary string of 100 bits. Then, via decoding, the binary number of each gene in the chromosome is transformed to a decimal system and for obtaining simulated values of biomass, substrate and product, a set of nonlinear differential equations (Equation (1)(1) $\frac{dX(t)}{dt}={\mu }_{\max }\left[1-\frac{X(t)}{X_{\max }}\right]X(t)$(1) –(3(3) $-\frac{dS(t)}{dt}=\frac{1}{Y_{X/S}}\frac{dX(t)}{dt}+m_{s}X(t)$(3) )) is solved by the fourth-order Runge–Kutta method (equation kinetic parameters are replaced with the ones estimated by GAs). Next, the value of the penalized objective function f(x) to be minimized is calculated: (6) $Minimizef\left(x\right)x={(x_1{}_{{}_{i1}},x_2{}_{{}_{i2}},\mathrm{...},x_d{}_{{}_{id}})}^T\in F\subseteq W\subseteq R^n$(6)

Subject to (7) $L_1(x)=x_{d_i}-(x_d^l+i\frac{x_d^u-x_d^l}{m})\le 0$(7) (8) $L_2(x)=-(x_{d_i}-(x_d^l+(i-1)\frac{x_d^u-x_d^l}{m}))\le 0$(8)

The penalized objective function f(x) is defined as (9) $f(x)=\phi (x)+p(x)$(9) (10) $\phi (x)={{\displaystyle \sum _{j=1}^z\left(\frac{X_{\exp }(j)-X_{\text{pred}}(j)}{X_{\text{pred}}(j)}\right)}}^2+{{\displaystyle \sum _{j=1}^z\left(\frac{S_{\exp }(j)-S_{\text{pred}}(j)}{S_{\text{pred}}(j)}\right)}}^2+{{\displaystyle \sum _{j=1}^z\left(\frac{P_{\exp }(j)-P_{\text{pred}}(j)}{P_{\text{pred}}(j)}\right)}}^2$(10) (11) $p(x)={\displaystyle \sum _{d=1}^n{\displaystyle \sum _{i=1}^m\max \left\{0,L_1(x)\right\}}}+{\displaystyle \sum _{d=1}^n{\displaystyle \sum _{i=1}^m\max \left\{0,L_2(x)\right\}}}$(11)

If $L_1(x)\le 0$ and $L_2(x)\le 0$ , $\max \left\{0,L_1(x)\right\}$ and $\max \left\{0,L_2(x)\right\}$ will be zero. Therefore, the constraint does not affect f(x), but if the constraint is violated, that means $L_1(x)>0$ or $L_2(x)>0$ , f(x) will increase. where (12) $\{\begin{array}{ll}{L}_1(x)=x_{d_i}-(x_d^l-i\frac{x_d^u-x_d^l}{m})\&L_2(x)=-(x_{d_i}-(x_d^l+(i-1)\frac{x_d^u-x_d^l}{m}))\hfill & \text{i}\text{f}{x}_d{}_{{}_i}\in F\hfill \\ L_1(x)=L_2(x)=0\hfill & \text{otherwise}\hfill \end{array}$(12) where p(x) is the penalty term; z is the number of experimental data items; X_exp, S_exp and P_exp are the experimental data for biomass (B. licheniformis), substrate (glucose) and product (alkaline protease); X_pred, S_pred and P_pred are the predicted values for biomass, substrate and production with a given set of parameters; $x={(x_1{}_{{}_{i1}},x_2{}_{{}_{i2}},\mathrm{...},x_d{}_{{}_{id}})}^T$ is a vector of feasible solutions; F is the feasible region and W is the whole search space.

In each feasible region, after sorting the populations according to descending f(x) values, encoding is done and a genetic selection procedure is applied. Then two GA operators, mutation and crossover, are applied. When the final criterion (maximum number of generations: 100) is met, the GA operators are stopped.

Next, in the PSO algorithm, the population has mⁿ particles. PSO particles search for the best solution across the solutions of GAs (from the previous step) in the whole search space. The position and velocity of particles are updated by applying an operator until particles can be expected to change towards the better solution. The updating procedure is defined by the following equations: (13) ${\omega }^t={\omega }_{\max }-\frac{({\omega }_{\max }-{\omega }_{\min }).t}{M}$(13)

• New velocity (14) $V_i^{t+1}={\omega }^t.V_i^t+c_1.r_1.\left({\text{pbest}}_i^t-X_i^t\right)+c_2.r_2.\left({\text{gbest}}^t-X_i^t\right)$(14)

• New position (15) $X_i^{t+1}=X_i^t+V_i^{t+1}$(15)

where $V_i^t$ is the velocity of the ith particle in the tth iteration, $X_i^t$ is the position of the ith particle in the tth iteration, ${\text{pbest}}_i^t$ is the best previous position of the ith particle (personal best), ${\text{gbest}}^t$ is the best previous position among all particles in the tth iteration (global best); $c_1$ and $c_2$ are two positive constants, called ‘cognitive’ and ‘social’ parameter, respectively, which are normally taken as 2; $r_1$ and $r_2$ are random numbers in the range of [0,1]; their role is to keep the population diversity; ${\omega }^t$ is the inertia weight factor and controls the impact of the previous velocity of the particle. M is the maximum number of iterations (in this study M = 100) and t is the current number of iterations. Also, ${\omega }_{\max }$ and ${\omega }_{\min }$ are the maximum and the minimum value of the weighting factor, respectively (${\omega }_{\min }=0.4\begin{array}{cc}& \begin{array}{cc}\&& {\omega }_{\max }=0.9\end{array}\end{array}$ ).

It should be noted that after decoding, the binary number of particles' position and velocity is transformed to a decimal system and for obtaining simulated values of biomass, substrate and product, a set of nonlinear differential equations (Equation (1)(1) $\frac{dX(t)}{dt}={\mu }_{\max }\left[1-\frac{X(t)}{X_{\max }}\right]X(t)$(1) –(3(3) $-\frac{dS(t)}{dt}=\frac{1}{Y_{X/S}}\frac{dX(t)}{dt}+m_{s}X(t)$(3) )) is solved by the fourth-order Runge–Kutta method (the equation kinetic parameters are replaced with the ones estimated by particles, in each iteration). Next, in each PSO algorithm iteration, the personal best and global best are determined according to the values of a penalized objective function similar to the previous step except that the search is done in the whole search space: (16) $\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}ef(x)x={(x_1,x_2,\mathrm{...},x_d)}^T\in W\subseteq R^{n}d=1,2,\mathrm{...},n$(16)

Subject to (17) $L_1(x)=x_d-x_d^u\le 0$(17) (18) $L_2(x)=-(x_d-x_d^l)\le 0$(18) (19) $f(x)=\phi (x)+p(x)$(19) (20) $\phi (x)={{\displaystyle \sum _{j=1}^z\left(\frac{X_{\exp }(j)-X_{\text{pred}}(j)}{X_{\text{pred}}(j)}\right)}}^2+{{\displaystyle \sum _{j=1}^z\left(\frac{S_{\exp }(j)-S_{\text{pred}}(j)}{S_{\text{pred}}(j)}\right)}}^2+{{\displaystyle \sum _{j=1}^z\left(\frac{P_{\exp }(j)-P_{\text{pred}}(j)}{P_{\text{pred}}(j)}\right)}}^2$(20) (21) $p(x)={\displaystyle \sum _{d=1}^n\max \left\{0,L_1(x)\right\}}+{\displaystyle \sum _{d=1}^n\max \left\{0,L_2(x)\right\}}$(21) where (22) $\{\begin{array}{ll}{L}_1(x)=x_d-x_d^u\&L_2(x)=-(x_d-x_d^l)\hfill & if{x}_d\in W\hfill \\ L_1(x)=L_2(x)=0\hfill & \text{otherwise}\hfill \end{array}$(22)

When the final criterion (maximum number of iterations: 100) is met, the PSO algorithm is stopped. Next, the solution is investigated. If the value of the penalized objective function (Equation (16)(16) $\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}ef(x)x={(x_1,x_2,\mathrm{...},x_d)}^T\in W\subseteq R^{n}d=1,2,\mathrm{...},n$(16) ) is less than $ϵ$ ($f(x)\le ϵ\begin{array}{cc}& \begin{array}{cc}\&& ϵ=0.5\end{array}\end{array}$ ), values of parameters would be printed; otherwise, the number of feasible region in first step will be double (m = 2m) and the procedure is repeated. The flowchart of the hybrid GA/PSO is shown in Figure 1.

Figure 1. Flowchart of hybrid GA/PSO.

All calculations of hybrid GA/PSO were implemented in Matlab R2012b (MathWorks Inc., Massachusetts, USA).

Results and discussion

Figure 2 shows the comparison between simulated data by hybrid GA/PSO and experimental data in batch mode. The simulated data are obtained by solving Equations (1)(1) $\frac{dX(t)}{dt}={\mu }_{\max }\left[1-\frac{X(t)}{X_{\max }}\right]X(t)$(1) –(3(3) $-\frac{dS(t)}{dt}=\frac{1}{Y_{X/S}}\frac{dX(t)}{dt}+m_{s}X(t)$(3) ) by the fourth-order Runge–Kutta method (the equation parameters are replaced with the ones estimated by the hybrid GA/PSO method).

Figure 2. Experimental concentrations of X (٭), S (○) and P (◊) vs. hybrid GA/PSO simulated data (solid lines). Note: X, biomass (B. licheniformis); S, substrate (glucose); P, product (protease). Initial concentration of substrate 40 g/L.

Table 1 shows the parameter estimation results by the hybrid GA/PSO method. The performance of this method is evaluated by calculating the value of the penalized objective function (f(x)). The value of f(x) in the hybrid GA/PSO method was obtained $\phi =0.0561$ . This proves the hybrid GA/PSO method to have a good ability to estimate kinetic parameters of alkaline protease production in batch mode.

Table 1. Kinetic parameters estimated by hybrid GA/PSO algorithm.

Kinetic parameters	Value (estimated by hybrid GA/PSO)
µmax	0.5442
K	0.0146
YP/X	0.6013
ms	0.0137
YX/S	0.5291

For comparison of the performance of the hybrid algorithm, we used pure GA and the pure PSO. Figure 3 shows the comparison between experimental data and simulated ones by GA and PSO, respectively.

Figure 3. Comparison of the experimental and simulated data of batch fermentation by using GA (A) and PSO (B). Experimental concentrations of X (٭), S (○) and P (◊) vs. simulated data (solid lines).

Figures 2 and 3 show that all the three algorithms successfully model the dynamics of protease production in batch mode. Nevertheless, the values of the penalized objective function and the computation time that are presented in Table 2 reveal that the hybrid GA/PSO fits the real experimental data more precisely than conventional algorithms.

Table 2. Comparison of performance across GA, PSO and hybrid GA/PSO algorithm.

Performance
Penalized objective function $(f(x))$	AlgorithmComputation time (s)
0.0956	104.92	GA
0.0748	92.13	PSO
0.0561	48.34	Hybrid GA/PSO

We observed that our hybrid algorithm achieves results similar to the pure algorithms (the value of the objective function is close to these obtained by the pure GA and the pure PSO algorithms), but it uses less computation time (t = 48.34 s). The time used is about two times less compared to that used by the other ones.

Also the results showed that the computational effort required by PSO to arrive at such high quality solutions is less than the effort required to arrive at the same high quality solutions by the GA (t = 92.13 s for PSO and t = 104.92 s for GA). Furthermore, PSO has slightly better performance and possesses better convergence property than GA because PSO is easy to implement and there are few parameters to adjust. Thus, it must be concluded that the PSO approach is superior to the GA approach in this particular case.

Overall, the results indicate that both GA and PSO can be used in the parameter estimation of the biological models. In terms of computational effort, the PSO approach is faster, although it should be noted that neither algorithm takes what can be considered an unacceptably long time to determine their results. In our particular case, the hybrid GA/PSO technique gave successful results in parameter estimation because optimizers work jointly to efficiently locate quality design points better than what either of these could alone. The results prove that this hybrid algorithm can be considered a good and often superior technique for solving parameter estimation problem in biological models by fitting the simulation data with the corresponding experimental measurements.

Conclusions

Fermentation processes are very complex and thus it is often very difficult to draw a comprehensive picture of what actually takes place in a given fermentation process of interest. Therefore, the estimation of growth kinetic parameters in fermentation processes has become the primary objective for investigation of fermentation behaviour. In this study, a hybrid GA/PSO technique was successfully used for parameter estimation of simple kinetic models for protease production by using experimental data obtained from batch bacterial fermentation in submerged culture. The hybrid GA/PSO technique can be considered superior to pure GA or pure PSO because it optimizes the work jointly to efficiently locate quality design points better than either of these could alone. The proposed hybrid GA/PSO is also more effective in terms of execution time and solution quality, since it benefits from synergy. The results prove that this hybrid algorithm is a good and often superior technique for parameter estimation problem in biological models by fitting simulation data with the corresponding experimental measurements.

Acknowledgments

The authors would like to express their profound gratitude to Keivan Bolouri (Department of Industrial Engineering, Science and Research Branch, Islamic Azad University, Tehran, Iran) for his invaluable guidance and helpful comments.

Disclosure statement

No potential conflict of interest was reported by the authors.