Biochemical systems identification by a random drift particle swarm optimization approach

Abstract

Background

Finding an efficient method to solve the parameter estimation problem (inverse problem) for nonlinear biochemical dynamical systems could help promote the functional understanding at the system level for signalling pathways. The problem is stated as a data-driven nonlinear regression problem, which is converted into a nonlinear programming problem with many nonlinear differential and algebraic constraints. Due to the typical ill conditioning and multimodality nature of the problem, it is in general difficult for gradient-based local optimization methods to obtain satisfactory solutions. To surmount this limitation, many stochastic optimization methods have been employed to find the global solution of the problem.

Results

This paper presents an effective search strategy for a particle swarm optimization (PSO) algorithm that enhances the ability of the algorithm for estimating the parameters of complex dynamic biochemical pathways. The proposed algorithm is a new variant of random drift particle swarm optimization (RDPSO), which is used to solve the above mentioned inverse problem and compared with other well known stochastic optimization methods. Two case studies on estimating the parameters of two nonlinear biochemical dynamic models have been taken as benchmarks, under both the noise-free and noisy simulation data scenarios.

Conclusions

The experimental results show that the novel variant of RDPSO algorithm is able to successfully solve the problem and obtain solutions of better quality than other global optimization methods used for finding the solution to the inverse problems in this study.

Background

Evolutionary algorithms (EAs) have been widely used for data mining tasks in Bioinformatics and Computational Biology [1, 2]. They are random search methods inspired by natural mechanisms existing in the biological world [1, 2]. EAs originally comprised four types of paradigms, namely, genetic algorithms (GAs), genetic programming (GP), evolution strategies (ES), and evolutionary programming (EP), with GAs being the most popular one. Data analysis tools traditionally used in Bioinformatics were mainly based on statistical techniques, such as regression and estimation, and EAs played significant roles in handling large biological data sets in a robust and computationally efficient manner [2].

Currently, evolutionary computing techniques mostly comprise conventional EAs (GAs, GP, ES and EP), swarm intelligence algorithms, artificial immune systems, differential evolution, as the main representative classes of evolutionary computing approaches[3]. Swarm intelligence is a class of evolutionary computing techniques simulating natural systems composed of many individuals that coordinate one another using decentralized control and self-organization. Two most influential and classical examples of swarm intelligence approaches are particle swarm optimization (PSO) and ant colony optimization (ACO) algorithms, which have been widely used in many different fields [3–7]. Particularly, PSO algorithms have shown their effectiveness in data mining tasks in bioinformatics due to their performance in solving difficult optimisation tasks [8–10].

Biochemical modelling can be considered a generic data-driven regression problem on the given experimental data. The goal of biochemical modeling is to build the mathematical formulations that quantitatively describe the dynamical behaviour of biochemical processes. For example, metabolic reactions are formulated as rate laws and described as a system of differential equations, the kinetic parameters of which are identified from a set of experimental data. Finding the solution of the parameter estimation problem, thus, plays a key role in building a dynamic model for a biochemical process, which, in turn, can help understand the functionality of the signalling pathways at the system level [11, 12].

Since solving the inverse problem in biochemical process modelling involves a task of nonlinear programming, many numerical optimization methods have been used to determine the parameters of biochemical models. These methods can be generally classified into two categories, namely, local optimization methods and global optimization methods [13]. The widely used local optimization tools for inverse problems are those based on gradient descent methods, the most popular being the Newton method [14]. This type of approaches, however, cannot be applied to non-smooth problems, since the objective functions of the problems are discontinuous or have discontinuous derivatives. Direct search methods, such as the Hooke-Jeeves method, the Needler-Mead simplex algorithm and the Downhill simplex algorithm, are also a kind of local optimization techniques that could be used to find a local minimum without the information from derivatives [13]. Normally, most local optimization approaches are used as single shooting methods. For each of them, the path of its optimization process leading to a final solution is determined by the initial conditions for the state variables. Therefore, the algorithm will lead to a wrong minimum, particularly if the initial conditions depend on model parameters. To overcome this shortcoming, one can adopt multiple shooting methods in which the time interval is partitioned and new initial conditions are used at the start of each time interval part [15]. The methods can offer the possibility to circumvent local optima by enlarging the parameter space during the optimization process.

The aforementioned local search methods are generally less efficient for the inverse problems of biochemical models, which are multimodal and high-dimensional. In order to solve these hard inverse problems efficiently, one can turn to global optimization methods, most of which incorporate stochastic search strategies to prevent the search process from being stuck into the local optimal or suboptimal solutions. The branch-and-bound approach is a global optimization method that converts the inverse problem into a convex optimization problem so that a global optimal solution can be obtained [16]. This method requires a finite search space that can be divided into smaller subspaces. A remarkable disadvantage is that it is applicable only if the lower and upper bounds of the objective function can be computed. Simulated annealing (SA) can be effectively used for parameter estimation from time-course biochemical data as shown in [17]. However, it has a slow convergence speed and high computational cost, and is not easy to be parallelized. Genetic algorithms (GAs) represent a widely used global search technique that could be employed to predict the parameters of dynamic models [18]. Nevertheless, GAs are always complained of slow convergence speed and high computation cost. The evolutionary strategy (ES) approach showed its ability to successfully solve inverse problems in a performance comparison made by Moles et al. [19] among a number of global optimization techniques on biochemical system identification problems. In contrast to SA, evolutionary algorithms, including ES and GAs, can be implemented as self-tuning methods and can be parallelizable, with the stochastic ranking evolutionary strategy (SRES) method being a very successful example [19–21]. Scatter search (SS) is known as a population-based random search approach that was proposed to identify the appropriate parameters for nonlinear dynamic biological systems [22, 23]. As an evolutionary algorithm method, the SS method, as well as its hybrid with a local search step after the recombination operation, showed to be efficient in solving inverse problems. Particle swarm optimization (PSO), also a population-based optimization technique from swarm intelligence and evolutionary computation area, has demonstrated its better performance than GAs in solving inverse problems [24, 25]. Hybrids of PSO with other methods have also shown their effectiveness in modelling biochemical dynamic systems [26–28]. However, PSO shows to be sensitive to the neighbourhood topology of the swarm, as commented in [29].

Other methods for parameter estimation include the Newton-flow analysis [30], the alternating regression technique [31], decoupling approaches [32, 33], the collocation method [20, 34], the decomposing method [35, 36]. These approximation techniques, when incorporated into an optimization algorithm, can help reduce the number of objective function evaluations, which are very computationally expensive. Additionally, radial basis function neural networks [37] and a quantitative inference method [38] have also been employed to solve inverse problems in biochemical process modelling.

In all of the above cases, the optimization approach is used to minimize to the residual error of an inferred model against experimental data. Smaller error means that the model describes the dynamic behaviour of the biochemical system better and has more justification to be accepted as a valid mathematical representation of the system. Theoretically, the prediction error diminishes with the accuracy of the model increasing. This study focuses on developing an efficient optimization method for parameter estimation of a given dynamic biochemical model. However, since parameter estimation problems of complex dynamic systems (generally with many parameters and many differential equations) are high-dimensional, multimodal and more challenging to solve, but which allow to depict more complex biochemical processes, our goal in this study is to develop an efficient global optimization method for solving such inverse problems of complex biochemical dynamic systems.

After extensive and in-depth study, we selected the PSO algorithm as a candidate to be modified in order to achieve our goal of solving complex inverse problems. The reason why PSO attracted us is that PSO has many advantages, such as faster convergence speed, lower computational need, as well as being easily parallelizable and having fewer parameters to adjust. However, PSO has the following shortcomings. First of all, it was theoretically proven that the PSO is not a global convergent algorithm, even not a local convergent one, against the convergence criteria given in [39]. Practically, the algorithm is more prone to be trapped into local optimal or suboptimal points for a high-dimensional problem, due to the weakening of its global search ability during the mid and later stages of the search process. Next, PSO is widely known to be sensitive to its search parameters including upper limits of the velocity, and even to the "swarm topology", so that users may feel awkward when selecting the parameters and the topologies when using the algorithm [40]. Finally, the performance of PSO appears to be very sensitive to the setting of upper and lower bounds of the search scope [40]. If the global optimal solution is located near the boundary of the search scope, the algorithm may have little chance to catch it. We have found that these shortcomings are mainly attributed to the velocity update equation, which is the essence of the PSO algorithm, and where it seems to be much room for improvement so as to boost the global search ability of the PSO.

In this study, inspired by the free electron model in metal conductors placed in an external electric field [41], we propose to use a variant of the PSO algorithm, called the random drift particle swarm optimization (RDPSO), in order to achieve our goal of effectively estimating the parameters of complex biochemical dynamical systems. The motivation of the RDPSO algorithm is to improve the search ability of the PSO algorithm by fundamentally modifying the update equation of the particle's velocity, instead of by revising the algorithm based on the original equation so as to probably increase the complexity of the algorithmic implementation as well as its computational cost. It is different from the drift particle swarm optimization (DPSO) proposed by us in [42, 43] in that it can make a better balance between the global search and the local search of the particle swarm.

The original and basic RDPSO version was recently introduced by us in [44], which was used for solving other problems in [45, 46]. A novel variant of RDPSO algorithm is being proposed in this work to solve the parameter identification problem for two biochemical systems. The novel variant proposed here is different from the original one in that it employs an exponential distribution for sampling the velocity of the particles, whilst the original one used the Gaussian distribution.

The novel RDPSO variant is used for estimating the parameters of two benchmark models, one of which describes the thermal isomerization of α-pinene with 5 parameters [22, 47], the other of which has a three-step pathway with 36 parameters [19]. The results of RDPSO and some other well-known global optimization algorithms are then compared and discussed. It should be noted that although this paper is focused on the parameter estimation for biochemical modelling, just as PSO and other EAs, the proposed RDPSO variant can be employed as a general-purpose tool for optimization problems in data miming tasks, such as clustering, classification, regression, and so forth, which widely exist in bioinformatics and computational biology [1, 2, 8, 48, 49].

Methods

Problem statement

The inverse problem of a nonlinear dynamic system involves finding proper parameters so as to minimize the cost function of the model with respect to an experimental data set, with some given differential equality constraints as well as other algebraic constraints. Such a data-driven regression problem can be approached with statistical techniques, using the given experimental data and the proposed models with unknown parameters. As stated by Moles et al. [19], the problem can be mathematically formulated as a nonlinear programming problem (NLP) with differential-algebraic constraints, whose goal is to find$\theta$ so as to minimize

$$J=\underset{0}{\overset{t_f}{\int }}{\left(y_{\text{msd}}\left(t\right)-y\left(\theta ,t\right)\right)}^{T}W\left(t\right)\left(y_{\text{msd}}\left(t\right)-y\left(\theta ,t\right)\right)dt$$

(1)

subject to

$$f\left(\frac{d\text{x}}{dt},\text{x},\text{y},\theta ,\text{v},t\right)=0$$

(2)

$$\text{x(}{t}_0)={\text{x}}_0$$

(3)

$$h\text{(x,y},\theta ,t)=0$$

(4)

$$g\text{(x,y},\theta ,t)\le 0$$

(5)

$${\theta }^L\le \theta \le {\theta }^U$$

(6)

where J is the cost function of the model, $\theta$ is a vector of model parameters to be estimated, $y_{\text{msd}}\left(t\right)$ is the experimental measure of a subset of the output state variables, $y\left(\theta ,t\right)$ is the prediction of those outputs by the model, x is the differential state variables and v is a vector of other (usually time-invariant) parameters that are not to be estimated. In Equation (1), $W\left(t\right)$ is the weighting (or scaling) matrix, and the equation can be discretized into a weighted least-square estimator. In Equation (2), f is the set of differential and algebraic equality constraints describing the system dynamics (i.e., the nonlinear process model). Equation (3) gives the initial value of x. In Equations (4) and (5), h and g are equality and inequality path and point constraints on system performance. In addition, $\theta$ is subject to upper and lower bounds, which are described by inequality constraints (6).

The above defined inverse problem is generally a multimodal (non-convex) optimization problem with multiple local optima due to the nonlinearity and constraints of the system dynamics. Even though many local and global optimization methods have been proposed to solve the problem as mentioned in Introduction, it is still challenging and very necessary to develop efficient optimization algorithms to deal with the parameter estimation problems, especially those for the dynamic systems with many parameters and many equations. Therefore this study focuses on the optimization approach for the inverse problem using the proposed variant of random drift particle swarm optimization (RDPSO) and other global optimization methods.

Particle swarm optimization

The original PSO algorithm was introduced by Kennedy and Eberhart in [50]. The algorithm was inspired by the observed social behavior of bird flocks or fish schooling, and it roots its methodology both in evolutionary computing and artificial life. It shares many similarities with EAs, in that both the PSO and the EAs are initialized randomly with a population of candidate solutions and then update the population iteratively, in order to approximate the global optimal solution to the given problem. However, unlike EAs, PSO has no evolution operators such as crossover and mutation, but perform optimization tasks by updating the particles' position (potential solutions) according to a set of discrete differential equations. It was shown that the PSO algorithm has comparable and even better performance than GAs [51]

In the PSO with m particles, each particle i ($1\le i\le m$), representing a potential solution of the given problem in a D-dimensional space, has three vectors at the k^th iteration, namely, the current position $X_i^k=\left(X_{i,1}^k,X_{i,2}^k,\cdots ,X_{i,D}^k\right)$, the velocity $V_i^k=\left(V_{i,1}^k,V_{i,2}^k,\cdots ,V_{i,D}^k\right)$ and its personal best (pbest) position $P_i^k=\left(P_{i,1}^k,P_{i,2}^k,\cdots ,P_{i,D}^k\right)$, which is defined as the position with the best objective function value found by the particle since initialization. A vector$G^k=\left(G_1^k,G_2^k,\cdots ,G_D^k\right)$, called the global best (gbest) position, is used to record the position with the best objective function value found by the all the particles in the particle swarm since initialization. With the above specification, the update equations for each particle's velocity and current position are given by:

$$V_{i,j}^{k+1}=w\cdot V_{i,j}^k+c_1{r}_{i,j}^k\left(P_{i,j}^k-X_{i,j}^k\right)+c_2{R}_{i,j}^k\left(G_j^k-X_{i,j}^k\right)$$

(7)

$$X_{i,j}^{k+1}=X_{i,j}^k+V_{i,j}^{k+1}$$

(8)

for$i=1,2,\cdots m;j=1,2\cdots ,D$, where $c_1$ and $c_2$ are known as acceleration coefficients and w is called the inertia weight, which can be adjusted to balance the exploration and exploitation ability of each particle [52]. Without loss of generality, we assume that the PSO is used to solve the following minimization problem:

$$\begin{array}{cc}\hfill Minimize\hfill & \hfill f\left(X\right),s.t.\hfill \end{array}X\in S\subseteq R^D$$

(9)

where $f\left(X\right)$is an objective function (or fitness function) and S denotes the feasible space. Consequently, $P_i^k$ can be updated by:

$$P_i^k=\{\begin{array}{c}\begin{array}{cc}\begin{array}{cc}{X}_i^k\hfill & \text{ if}\hfill \end{array}\hfill & f(X_i^k)<f(P_i^{k-1})\hfill \end{array}\hfill \\ \begin{array}{ccc}{P}_i^{k-1}\hfill & \text{if}\hfill & f(X_i^k)\ge f(P_i^{k-1})\hfill \end{array}\hfill \end{array}$$

(10)

and $G^k$ can be found by:

$$G^k=P_g^k,\text{ where }g=\underset{1\le i\le m}{\text{arg min}}\left[f\left(P_i^k\right)\right]$$

(11)

In Equation (7), $r_{i,j}^k$ and $R_{i,j}^k$ are the sequences of two different random numbers with uniform distribution on the interval (0, 1), namely, $r_{i,j}^k,R_{i,j}^k\text{~}U\left(0,1\right)$. In order to prevent the particle from flying away out of the search scope, $V_{i,j}^k$ is restricted on the interval $\left[-V_{\text{max}},V_{\text{max}}\right]$, where $V_{\text{max}}$ is also a user-specified algorithmic parameter.

Many researchers have proposed different variants of PSO in order to improve the search performance of the algorithm and proved this through empirical simulation [3–7, 52–56]

The proposed random drift particle swarm optimization

In [57], it was demonstrated that if the acceleration coefficients are properly valued, each particle converges to its local attractor,$p_i^k=\left(p_{i,1}^k,p_{i,2}^k,\cdots p_{i,D}^k\right)$, so that the convergence of the whole particle swarm can be achieved. Each coordinate of $P_i^k$ is given by:

$$p_{i,j}^k=\frac{c_1{r}_{i,j}^k{P}_{i,j}^k+c_2{R}_{i,j}^k{G}_j^k}{c_1{r}_{i,j}^k+c_2{R}_{i,j}^k},1\le j\le D$$

(12)

which can be restated as

$$p_{i,j}^k={\varphi }_{i,j}^k{P}_{i,j}^k+\left(1-{\varphi }_{i,j}^k\right),1\le j\le D$$

(13)

where

$${\varphi }_{i,j}^k=\frac{c_1{r}_{i,j}^k}{c_1{r}_{i,j}^k+c_2{R}_{i,j}^k}$$

(14)

In the PSO algorithm, c₁ and c₂ are set to be equal, and thus ${\varphi }_{i,j}^k$ is a random number with uniform distribution on the interval (0, 1), i.e. ${\varphi }_{i,j}^k~U\left(0,1\right)$.

During the search process of the PSO, as particles' current position are converging to their own local attractor, their current positions, pbest positions, local attractors and the gbest positions are all converging to one single point. The directional movement of each particle i towards $P_i^k$ resembles the drift motion of an electron in metal conductors placed in an external electric field. According to the free electron model [37], the electron has not only drift motion caused by the external electric field, but also a thermal motion, which appears to be a random movement. The superposition of the drift thermal motions makes the electron careen towards the location of the minimum potential energy. Thus, if the position of an electron in the metal is regarded as a candidate solution and the potential energy function as the objective function to be minimized, the movement of the electron resembles the process finding the minimum solution of the minimization problem.

The above facts can lead to a novel variant of PSO if the particle in PSO is assumed to behave like an electron moving in a metal conductor in an external electric field. More specifically, it can be assumed that at the k th iteration, each particle i has drift motion towards $P_i^k$ as well as a thermal motion, with their velocities in each dimension j denoted as $V{1}_{i,j}^{k+1}$ and $V{2}_{i,j}^{k+1}$, respectively. As a result, the velocity of the particle is given by $V_{i,j}^{k+1}=V{1}_{i,j}^{k+1}+V{2}_{i,j}^{k+1}$. In the drift particle swarm optimization proposed in [42, 43], we assume that $V{1}_{i,j}^{k+1}$ follows a Maxwell distribution, say a Gaussian probability distribution, and the $V{2}_{i,j}^{k+1}$ is given by the social part plus the cognitive part in Equation (7). This velocity update equation appears to add some effectiveness to the search performance of the particle swarm but has some shortcomings. Firstly, the Gaussian distribution has a thin tail so that it has less opportunity to generate outliers. As a result, the thermal motion of the particle has less randomness and cannot significantly improve the particle's global search ability. Secondly, although the update equation of $V{2}_{i,j}^{k+1}$:

$$V{2}_{i,j}^{k+1}=c_1{r}_{i,j}^k\left(P_{i,j}^k-X_{i,j}^k\right)+c_2{R}_{i,j}^k\left(G_j^k-X_{i,j}^k\right)$$

(15)

can guarantee the particle to converge towards its local attractor, the two random scaling coefficients add randomness to its motion, which means that the particle's position is sampled at uniformly random positions within the hyper-rectangle around the gbest position and its personal best position. It is not able to enhance the particle's global search ability because of the finite scope of the hyper-rectangle, but it may weaken its local search ability, which is the responsibility of the directional motion. Therefore, the velocity of the particle, which is given by the sum of $V{1}_{i,j}^{k+1}$ and$V{2}_{i,j}^{k+1}$, may not be able to make a good balance between the global search and the local search of the particle. In the present study, we employ a new way of determining $V{1}_{i,j}^{k+1}$ and$V{2}_{i,j}^{k+1}$.

Here, we assume that the velocity of the thermal motion $V{1}_{i,j}^{k+1}$ follows a double exponential distribution, whose probability density function and probability distribution function are

$$f_{V{1}_{i,j}^{k+1}}\left(v\right)=\frac{1}{{\sigma }_{i,j}^k}{e}^{\frac{-2\left|v\right|}{{\sigma }_{i,j}^k}}$$

(16)

and

$$F_{V{1}_{i,j}^{k+1}}\left(v\right)=1-e^{\frac{-2\left|v\right|}{{\sigma }_{i,j}^k}}$$

(17)

respectively, where v represents the value of the random variable $V{1}_{i,j}^{k+1}$ and ${\sigma }_{i,j}^k$ is the standard deviation of the distribution. By employing a stochastic simulation method, we can express $V{1}_{i,j}^{k+1}$ as

$$V{1}_{i,j}^{k+1}=\frac{{\sigma }_{i,j}^k}{2}{\varphi }_{i,j}^k$$

(18)

$${\varphi }_{i,j}^k=\{\begin{array}{c}\begin{array}{cc}+\ln \left(1/u_{i,j}^k\right)& \begin{array}{cc}if& s>0.5\end{array}\end{array}\\ \begin{array}{ccc}-\ln \left(1/u_{i,j}^k\right)& if& s\le 0.5\end{array}\end{array},$$

(19)

where

where s and $u_{i,j}^k$ are two different random numbers uniformly distributed on the interval (0,1), i.e. $\begin{array}{cc}\hfill s,\hfill & \hfill u_{i,j}^k\hfill \end{array}~U\left(0,1\right)$. As for the value of ${\sigma }_{i,j}^k$, an adaptive strategy is adopted to determine ${\sigma }_{i,j}^k$ by ${\sigma }_{i,j}^k=2\alpha |C_j^k-X_{i,j}^k|$, where $C^k=\left(C_1^k,C_2^k,\cdots ,C_D^k\right)$ is called the mean best (mbest) position, defined as the mean of the pbest positions of all the particles, i.e. $\begin{array}{cc}{C}_j^k=(1/m){\displaystyle {\sum }_{i=1}^m{P}_{i,j}^k}& (1\le j\le D)\end{array}$.

The velocity of the drift motion $V{2}_{i,j}^{k+1}$ may have many possible forms. However, the following simple linear expression is adopted in this study:

$$V{2}_{i,j}^{k+1}=\beta \left(p_{i,j}^k-X_{i,j}^k\right)$$

(20)

where $p_{i,j}^k$ is determined by

$$p_{i,j}^k={\varphi }_{i,j}^k{P}_{i,j}^k+(1-{\varphi }_{i,j}^k),\begin{array}{cc}\begin{array}{c}{\varphi }_{i,j}^k~U(0,1),\end{array}& 1\le j\le D\end{array}$$

(21)

It can be immediately proven that if $V_{i,j}^{k+1}=V{2}_{i,j}^{k+1}$, when $k\to \infty$, $X_{i,j}^k\to p_{i,j}^k$. Therefore the expression of $V{2}_{i,j}^{k+1}$ in Equation (20) can indeed guarantee that the particle move directionally to $P_i^k$ as an overall result.

With the definitions of the thermal motion and the drift motion of the particle, we can obtain a novel set of update equations for the particle:

$$V_{i,j}^{k+1}=\alpha |C_j^k-X_{i,j}^k|{\varphi }_{i,j}^k+\beta \left(p_{i,j}^k-X_{i,j}^k\right)$$

(22)

$$X_{i,j}^{k+1}=X_{i,j}^k+V_{i,j}^{k+1}$$

(23)

where $\alpha$ is called the thermal coefficient and $\beta$ is called the drift coefficient. The PSO with Equations (22) and (23) is a novel variant of RDPSO, which employs a double exponential distribution instead of a Gaussian one. The procedure of this RDPSO variant is outlined below.

Step 0: Randomly initialize the current positions and the pbest position of all the particles;

Step 1: Set k = 0;

Step 2: While the termination condition is not met, do the following steps;

Step 3: Set k=k+1 and compute the mbest position $C^k$, which is the centroid of the pbest positions of all the particles at iteration k;

Step 4: From $i=1$, carry out the following steps;

Step 5: Evaluate the objective function value $f\left(X_i^k\right)$, and update $P_i^k$ and $G^k$ according to Equation (10) and Equation (11), respectively;

Step 6: Update the components of the velocity and current position of particle i in each dimension, respectively, according to Equations (21), (22) and (23);

Step 7: Set i=i+1, and return to Step 5 until $i=m$;

Step 8: Return to Step 2;

In the RDPSO algorithm, in addition to the population size m, $\alpha$ and $\beta$ are two very important user-specified algorithmic parameters, which play the same roles as the inertia weight w and acceleration coefficients in the basic PSO algorithm. That is, the can be tuned to balance the exploration and exploitation ability of the particle. How to select the values of these parameters to prevent the particles from explosion is an open problem. Here, we performed the stochastic simulations for the one dimensional case, in which the local attractor was fixed at the origin and the mbest position was at x = 0.1. The results of two simulations are visualized in Figure 1 and Figure 2, in which the logarithmic value of the absolute of $X_k$ was recorded as ordinate, and the iteration number was the abscissa. Figure 1 shows that the particle's position was bounded when $\alpha =1$ and $\beta =1.5$. However, when $\alpha =1.8$ and $\beta =1.5$, the particle diverged to infinity. To obtain the sufficient and necessary condition for the particle to be bounded, we will focus our attention on a theoretical analysis in terms of probability measure in future.

Figure 1

The result of the stochastic simulation with$\alpha =1$and $\beta =1.5$.It is shown that, with this parameter setting, the particle's position is bounded.

Figure 2

The result of the stochastic simulation with$\alpha =1.8$ and $\beta =1.5$.It is shown that, with this parameter setting, the particle's position diverges as the iteration number increases.

Setting large values for $\alpha$ and $\beta$ implies better global search ability of the algorithm, while setting small values means better local search. When the RDPSO is used for solving a problem, a good balance between the global search and the local search of the algorithm is crucial for the algorithmic performance. However, in order to find out how to tune the parameters to generate generally good algorithmic performance we need a large number of experiments on benchmark functions, which will be performed in our future tasks. Here, we recommend that when the RDPSO is used, $\alpha$ should be set to be no larger than 1.0 and $\beta$ to be no larger than 1.5. More specifically, when the problem at hand is complex, the values of the two parameters should be set to be relatively large in order to make the particles search more globally, and on the other hand, when the problem is simple, relatively smaller values should be selected for the parameters, for the purpose of faster convergence speed of the algorithm. In the present study, the value of $\alpha$ and $\beta$ were set to be 0.75 and 1.0, respectively.

In addition, the population size and the maximum number of iterations (MaxIter) can also affect the performance of a population-based technique. Just as for other PSO variants, it is suggested that the population size should be larger than 20 for the RDPSO as well. The value of the MaxIter depends on the complexity of the problem. Generally, a smaller MaxIter value is used for simple problems, while a larger one is used for complex problems.

Moreover, $V_{i,j}^k$ is also restricted within the interval $\left[-V_{\text{max}},V_{\text{max}}\right]$during the search process of the RDPSO algorithm, just as in the original PSO algorithm.

The optimization methods compared

Besides the PSO and RDPSO algorithms, the Differential Evolution (DE), Scatter Search (SS) method and two versions of Evolution Strategies were also used to solve the selected inverse problems, for performance comparison purposes. The DE method, as presented by Storn and Price [58], is an evolutionary computing method, which has a faster convergence speed than GAs and can find the global optimal solution of a multidimensional and multimodal function effectively [58].

The SS method is also a population-based search techniques originally developed by Glover [59]. In [22], a novel meta-heuristic method, which is the combination of the original SS method with a local search technique, was proposed to solve inverse problems. It was shown that the local search operator can accelerate the convergence speed significantly. Thus, in our experiments, we used this novel SS method for performance comparison.

Evolutionary strategy (ES) is an important paradigm of EAs, which imitates the effects that genetics produces on the phenotype, rather than the genotype as in GAs [60]. The two canonical versions of ES we used in this study are denoted by (μ, λ)-ES and (μ + λ)-ES, where μ denotes the number of parents and λ the number of offspring. In the (μ, λ)-ES, the parents are deterministically selected from offspring (μ <λ must hold), while in (μ + λ)-ES, the parents are selected from both the parents and offspring.

In addition, the performances of the above mentioned algorithms, including the RDPSO, are also compared with those of the SRES method. The SRES is a version of (μ, λ)-ES that uses stochastic ranking to handle the constraints, by adjusting the balance between the objective function and the penalty function on the course of the search [19], [61].

Case studies

Two case studies involving two benchmark systems were carried out. For each system, we performed two groups of numerical experiments, one with noise-free simulation data, and the other with noisy simulation data.

Case study 1

The goal of this case study is to estimate the five rate constants of the homogeneous biochemical reaction describing the thermal isomerization of α-pinene, which is an organic compound of the terpene class, one of two isomers of pinene [22], [23]. The mathematical model of this process is formulated with the following linear equations:

$$\frac{d{y}_1}{dt}=-\left(p_1+p_2\right)y_1$$

(24)

$$\frac{d{y}_2}{dt}=p_1{y}_1$$

(25)

$$\frac{d{y}_3}{dt}=p_2{y}_1-\left(p_3+p_4\right)y_3+p_5{y}_5$$

(26)

$$\frac{d{y}_4}{dt}=p_3{y}_3$$

(27)

$$\frac{d{y}_5}{dt}=p_4{y}_3-p_5{y}_5$$

(28)

where $\left(p_1,p_2\cdots ,p_5\right)$ is the vector of unknown coefficients to be estimated, $y_1$, $y_2$, $y_3$, $y_4$ and $y_5$ denote the concentrations of the α-pinene, dipentene, alloocimen, β-pyronene and a dimer, respectively.

Case study 2

This case study involves the inverse problem to identify 36 kinetic parameters of a nonlinear biochemical dynamic model formed by the following 8 ordinary differential equations that describe the variation of the metabolite concentrates with time 1 [19].

$$\frac{d{G}_1}{dt}=\frac{V_1}{1+{\left(\frac{P}{K{i}_1}\right)}^{n{i}_1}+{\left(\frac{K{a}_1}{S}\right)}^{n{a}_1}}-k_1\cdot G_1$$

(29)

$$\frac{d{G}_2}{dt}=\frac{V_2}{1+{\left(\frac{P}{K{i}_2}\right)}^{n{i}_2}+{\left(\frac{K{a}_2}{M_1}\right)}^{n{a}_2}}-k_2\cdot G_2$$

(30)

$$\frac{d{G}_3}{dt}=\frac{V_3}{1+{\left(\frac{P}{K{i}_3}\right)}^{n{i}_3}+{\left(\frac{K{a}_3}{M_2}\right)}^{n{a}_3}}-k_3\cdot G_3$$

(31)

$$\frac{d{E}_1}{dt}=\frac{V_4\cdot G_1}{K_4+G_1}-k_4\cdot E_1$$

(32)

$$\frac{d{E}_2}{dt}=\frac{V_5\cdot G_2}{K_5+G_2}-k_5\cdot E_2$$

(33)

$$\frac{d{E}_3}{dt}=\frac{V_6\cdot G_3}{K_6+G_3}-k_6\cdot E_3$$

(34)

$$\begin{array}{c}\frac{d{M}_1}{dt}=\frac{kca{t}_1\cdot E_1\cdot \left(\frac{1}{K{m}_1}\right)\cdot \left(S-M_1\right)}{1+\frac{S}{K{m}_1}+\frac{M_1}{K{m}_2}}\\ \begin{array}{c}\hfill -\hfill \end{array}\frac{kca{t}_2\cdot E_2\cdot \left(\frac{1}{K{m}_3}\right)\cdot \left(M_1-M_2\right)}{1+\frac{M_1}{K{m}_3}+\frac{M_2}{K{m}_4}}\end{array}$$

(35)

$$\begin{array}{c}\frac{d{M}_2}{dt}=\frac{kca{t}_2\cdot E_2\cdot \left(\frac{1}{K{m}_3}\right)\cdot \left(M_1-M_2\right)}{1+\frac{M_1}{K{m}_3}+\frac{M_2}{K{m}_4}}\\ \begin{array}{c}\hfill -\hfill \end{array}\frac{kca{t}_3\cdot E_3\cdot \left(\frac{1}{K{m}_5}\right)\cdot \left(M_2-P\right)}{1+\frac{M_2}{K{m}_5}+\frac{P}{K{m}_6}}\end{array}$$

(36)

where$M_1$, $M_2$, $E_1$, $E_2$, $E_3$, $G_1$, $G_2$ and $G_3$ are the state variables representing the concentrations of the species involved in different biochemical reactions, and S and P are controlling parameters which are kept fixed at the initial values for each experiment. The inverse problem is then reduced to the optimization problem that fits the remaining 36 parameters represented by$\theta =\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_{36}\right)$.

Objective functions

The objective function (or fitness function) for the inverse problem in either of the two case studies is the discretization of Equation (1), which is formulated as the weighted sum of squares of the differences between the experimental and the predicted values of the state variables:

$$J={\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^l}{w}_{ij}{\left\{{\left[y_{\text{pred}}\left(i\right)-y_{\text{exp}}\left(i\right)\right]}_j\right\}}^2$$

(37)

where n is the number of data for each experiment, l is the number of experiments, $y_{\text{exp}}$ is the vector of experimental values of the state variables, and $y_{\text{pred}}$ is the vector of the values of state variables predicted by the model with a given set of parameters. In Case Study 1, each $w_{ij}$ was set to be 1 [22], while in Case Study 2, $w_{ij}$ was set as $w_{ij}={\left\{1/\text{max}{\left[y_{\text{pred}}\left(i\right)\right]}_j\right\}}^2$, which was used to normalize the contributions of each term [19].

Obtaining simulation data

In order to evaluate the performances of the global optimization methods in finding the solution of the inverse problems, we chose a set of parameters for each model, which are considered as the true or nominal values. For Case Study 1, the true values of the parameters are p₁ = 5.93e-5, p₂ = 2.96e-5, p₃ = 2.05e-5, p₄ = 27.5e-5 and p₅ = 4.00e-5. For Case Study 2, the nominal values of the model parameters are shown in Table 1.

Table 1 Experiment generation (simulation) and the nominal value.

The pseudo-experimental data (essentially the simulation data) in either case were generated by substituting the chosen parameters into the dynamic model and performing fourth order Runge-Kutta method on the corresponding differential equations. For Case Study 2, the pseudo-measurements of the concentrations of metabolites, proteins, and messenger RNA species were the results of 16 different pseudo-experiments, in which, with the given nominal values for the parameters, the initial concentrates of the pathway substrate S and product P were varied for each experiment (simulation) as shown in Table 1. These simulated data represent the exact experimental results devoid of measurement noise and they were used as noise-free data for the first group of numerical experiments in each case study. In order to test the optimization methods for noisy data, we added a white noise to each of the original noise-free data:

$$z^{\prime }=z+\sigma \epsilon$$

(38)

where z and $z^{\prime }$ represents the original noise-free data and the resulting noisy data, respectively, $\epsilon$ is an random number with standard normal distribution, namely, $\epsilon ~N\left(0,1\right)$, and $\sigma$ is the standard deviation of the white noise. In our case studies, $\sigma$ was set to 0.01 for both systems.

Initial problem solver used

During the search of each global optimization algorithm, each potential solution (i.e., a set of estimation values for the parameters) was substituted into the dynamic model. Then, the fourth order Runge-Kutta method was performed on the corresponding system of differential equations to generate a set of predicted values of the output state variables, from which the objective function value (or fitness value) of the potential solution could be evaluated according to Equation (37) with the obtained pseudo-experimental (simulation) data. This process is known as the solution to the forward problem, which was embedded in the iterations of the search during the solving of the inverse problem with the algorithm.

Experimental settings

For the sake of performance comparison, all the tested global optimization methods except the SS method (i.e., PSO, RDPSO, DE, (μ, λ)-ES, and (μ + λ)-ES) were programmed in C++ on a VC++6.0 platform in Windows XP environment, and implemented on Intel Pentium Dual-Core E5300 2.6GHz PCs, each with 2 MB cache and 4 GB main memory. The SS method was implemented in Matlab 7.0, on the same platform, for the purpose of calling the local solver SQP in Matlab during the search process. The software for SS for inverse problems can be found on http://www.iim.csic.es/~gingproc/ssmGO.html.

The configuration of the algorithm parameters including the population sizes are listed in Table 2. In Case Study 1, each optimization algorithm ran 20 times with each run executed for 500 iterations; that is, the maximum number of iterations (MaxIter) is 500. In Case Study 2, each algorithm also ran 20 times with each run executed for 2250000 function evaluations, which is the same as that for the DE in [19]. Since the population size of each algorithm was 100, the value of MaxIter was 22500 in Case Study 2. For the SS method, the initial population size was 100, and 10 individuals were selected to perform the iterative search after initialization. Other parameters were selected according to recommendations from the corresponding references and/or our preliminary runs. For all the algorithms tested on the inverse problems, the statistical values of J were figured out and the results with best values of J were selected, processed and visualized with Matlab 7.0.

Table 2 Configuration of search parameters in different algorithms (F is an algorithmic control parameter used to control the size of the perturbation in the mutation operator for DE, CR is the crossover constant in DE, varphi is a parameter determining the standard deviation of the mutation in the evolution strategy).

Results

For Case Study 1, the statistical values of J from 20 search runs with 500 iterations by each algorithm are listed in Table 3. The best value of J (J = 1.3740e-014) for the numerical experiments with noise-free simulation data was obtained by using our proposed RDPSO algorithm after running for 0.0348h (about 2 minutes). For the experiment with noisy data, the RDPSO generated the best J value (J = 0.2023) as well. The proposed algorithm also showed the best performance on average among all tested methods, as shown by the mean value of J over 20 runs. In this case study, the basic PSO algorithm showed good performance on low-dimensional inverse problems.

Table 3 J values and computational time for the global optimization methods in Case Study 1

The convergence process of each tested algorithm averaged over 20 runs in the numerical experiment with noisy data in Case Study 1 is shown by the convergence curve in Figure 3, which is plotted in the log-log scale with objective function values versus the iteration number. Evidently, the SS method showed a better convergence property than other algorithms. The best solution vector corresponding to the best value of J (J = 1.3740e-014) obtained by the RDPSO in the numerical experiment with noise-free data was p₁ = 5.930000e-005, p₂ = 2.960000e-005, p₃ = 2.750002e-004, p₄ = 2.750002e-004, p₅ = 4.000561e-005, extremely close to the real values of the parameters, and the best solution vector obtained by the RDPSO for the numerical experiment with noisy data was (when J = 0.2023) p₁ = 5.928818e-005, p₂ = 2.959918e-005, p₃ = 1.712580e-005, p₄ = 3.147252e-004, p₅ = 1.537512e-003. Figure 4and Figure 5 show good fits between the experimental data (simulation data) and the predicted data, with the best obtained parameters in the experiments with both noise-free and noisy data.

Figure 3

The figure visualizes the convergence curves of objective function values of all the algorithms averaged over 20 runs in the numerical experiments with noisy data for Case Study 1. It is shown that the RDPSO, PSO and SS methods had better convergence properties than other methods.

Figure 4

The figure shows the noise-free experimental data (marker) versus the predicted model (continuous line) for Case Study 1. It is shown that the predicted model obtained by RDPSO fits the experimental data well

Figure 5

The figure shows the noisy experimental data (marker) versus the predicted model (continuous line) for Case Study 1. It is also shown that the predicted model obtained by RDPSO fits the experimental data well.

For Case Study 2, the obtained values of J resulted from 20 runs with 22500 iterations by each algorithm are listed in Table 4. For the numerical experiments with noise-free data, the best result of J (J = 0.7358e-06) was obtained by using the SS method, which also had the better average performance than any other compared algorithm, as shown by the mean value of J over the 20 runs of the algorithm. The second best method was the RDPSO algorithm, which could converge to a value of J = 0.009124 and had a mean value of J = 0.178881 over 20 runs. Table 5 lists the estimated values of the model parameters corresponding to the best J value (J= 0.009124) found by the RDPSO algorithm. The results also show that (μ, λ)-ES is the winner in this inverse problem compared to (μ + λ)-ES, whose best and mean results of J were 0.123209 and 2.141820, respectively. The PSO and DE are two well-known efficient population-based optimization methods, which, however, could not arrive at the vicinity of the aforementioned solutions. When the experimental data (simulation data) was noisy, the SS method and RPSO obtained similar results for the best J value over 20 runs, but the former had a better average algorithmic performance. Table 6 shows the identified model parameters corresponding to the best J value (J = 0.2313) obtained by the RDPSO for the noisy data.

Table 4 J values and computational time of the global optimization methods for Case Study 2, including the results for the noise-free and noisy data.

Table 5 Decision vector for the solution found by RDPSO in the experiment with noise-free data for Case Study 2.

Table 6 Decision vector for the solution found by RDPSO in the experiment with noisy data for Case Study 2

We plotted in Figure 6 the convergence curve of each method averaged over 20 runs. It is shown that the SS method had a remarkably better convergence rate than others, probably due to its local solver that can enhance the local search ability of the algorithm significantly. Figures 7a and Figure 8 show comparisons between the predicted data and the experimental (simulation) data for the decision vectors found by the RDPSO in both groups of numerical experiments (with the noise-free and noisy data, respectively). It can be observed that there is good correlation between the experimental and predicted data.

Figure 6

The figure visualizes the convergence curves of objective function values of all the algorithms averaged over 20 runs of the numerical experiments with noisy data in Case Study 2. It is shown that the SS method had the fastest convergence speed and the RDPSO had the second fastest one.

Figure 7

The figure shows the predicted concentration provided by the RDPSO (continuous line) and noise-free experimental (marker) data in Case Study 2 when P = 1.0, S = 2.1544. It is shown that the predicted concentration had good correlation with the experimental data.

Figure 8

The figure shows the predicted concentration provided by the RDPSO (continuous line) and noisy experimental (marker) data in Case Study 2 when P = 1.0, S = 2.1544. It is shown that the predicted concentration had good correlation with the experimental data.

Conclusions

In this paper, a variant of RDPSO algorithm was proposed and showed to be able to successfully solve two inverse problems associated with the thermal isomerization of α-pinene and a three-step pathway, respectively. The results indicate that the proposed RDPSO algorithm outperformed its PSO predecessors and some other competitors in the first problem, and also had the second best algorithmic performance among all the compared algorithms.

Like other stochastic optimization methods, a possible drawback of the RDPSO method is the computational effort required. This is mainly because most of the computational time was spent on solving the forward problem. One measure that can be taken is to incorporate the local search technique into the algorithm in order to accelerate its convergence speed. Another is to develop a parallelized RDPSO implementation to solve inverse problems on computer clusters to reduce the computational cost to a reasonable level. Our future tasks will focus on these two ways of improving the algorithmic effectiveness of the RDPSO algorithm.

Availability and requirements

In the additional file 1 the source codes of five of the tested algorithms on the two benchmark systems are provided. It includes two file folds, one for benchmark system 1 and the other one for benchmark system 2. All the algorithms are programmed with C++ in Microsoft Visual C++ 6.0.

In additional file 2 the data files for the five algorithms used in the two case studies are provided. For each case, the data corresponding to the best results generated by 20 runs of each algorithm are provided in a .txt file.