Parameter estimation in systems biology models using spline approximation

Background Mathematical models for revealing the dynamics and interactions properties of biological systems play an important role in computational systems biology. The inference of model parameter values from time-course data can be considered as a "reverse engineering" process and is still one of the most challenging tasks. Many parameter estimation methods have been developed but none of these methods is effective for all cases and can overwhelm all other approaches. Instead, various methods have their advantages and disadvantages. It is worth to develop parameter estimation methods which are robust against noise, efficient in computation and flexible enough to meet different constraints. Results Two parameter estimation methods of combining spline theory with Linear Programming (LP) and Nonlinear Programming (NLP) are developed. These methods remove the need for ODE solvers during the identification process. Our analysis shows that the augmented cost function surfaces used in the two proposed methods are smoother; which can ease the optima searching process and hence enhance the robustness and speed of the search algorithm. Moreover, the cores of our algorithms are LP and NLP based, which are flexible and consequently additional constraints can be embedded/removed easily. Eight system biology models are used for testing the proposed approaches. Our results confirm that the proposed methods are both efficient and robust. Conclusions The proposed approaches have general application to identify unknown parameter values of a wide range of systems biology models.


Background
In recent years, the rapid development of sophisticated experiment tools in molecular biology allows the acquisition of high qualitative time series data which can significantly improve the ability of revealing the complex dynamics and interactions of biological systems. Profiting from the rapid technological advances, more and more researchers from different disciplines can now utilize such observation data to establish mechanism-based models which can incorporate every possible detail and functioning of biological systems [1]. One common approach is to characterize the biological system with a set of Ordinary Differential Equations (ODEs) [2][3][4][5][6][7]. Generally, there are two major aspects of building an ODE model for a biological system from experimentally measured time series: (1) to determine the structure of the system through a set of suitable ODEs with unknown parameters; (2) to determine the unknown parameters of this ODE model. The identification of these unknown parameter with fixed model structure from observations is one of the central issues of computational systems biology [8]. This type of approach can be considered as a "reverse engineering process" [9][10][11].
The parameter estimation problem is generally formulated as an optimization problem that minimizes an objective function which represents the fitness of the model with respect to a set of experimental data [8,[12][13][14][15][16][17]. Two major optimization approaches are commonly adopted; the gradient-based nonlinear optimization method and the evolutionary based method. Also, simulations had shown that the simulated annealing (SA) method can offer promising results [18]. In [19], many deterministic and stochastic global optimization (GO) methods for parameter estimation were further compared using a three-step pathway model with noise free data assumption; the best result was given by the Stochastic Ranking Evolution Strategy (SRES) method.
It is worth mentioning that, due to its simplicity in implementation, evolutionary algorithms, such as genetic algorithm and their variants, are extensively utilized for identifying unknown parameters of systems biology models [1,11,[20][21][22]. However, most of these aforementioned approaches need a numerical ODE solver to perform the numerical integration for the underlining differential equations. Studies have revealed that more than 90% of the computation time is consumed in the ODE solver during the identification process [19]. In particular, for nonlinear dynamical systems with highparameter-dimension, one trial usually consumes tens of hours or even days [10,20,21,23]. Furthermore, the convergence property is aggravated by numerical integration failure, which is a major problem in the optimization process [11]. The computational burden can be relieved by reformulating the system involving differential equations into a system of algebraic equations [12,15,17,24], which can be classified as "decomposition approaches". These decomposition approaches are widely employed in the parameter estimation of S-systems [7,25]. The reliability of the decomposed methods depends on the accuracy of the "smooth" estimated derivatives and the states of the system. In practice, these data are subject to significant observation noise. Without proper preprocessing, the estimation faces the potential of the overfitting problem and hence the estimation can deviate badly from the "true" value [26,27]. Regularization can be considered as a mathematical pre-processing on the measured noisy data set and be used to control the trade-off between the "roughness" of the solution and the infidelity of the data [28]. Since we are dealing with a known structured bio-system, the system model itself possesses a physical inertia and can serve as physical constraints which limit the system states within a set of possible trajectories. In this paper, the over-fitting problem can be relieved by embedding the model dynamics, the mass and energy balance constraints into our constrained optimization algorithms. Owing to the nonlinearity of systems biology models, the cost function to be minimized is complex and has multiple local minima. Minimization algorithms face the high possibility of getting trapped at local optima. For these reasons, the parameter estimation problem is still a bottleneck and a challenging task of computational analysis of systems biology [1,11]. Until now, none of the parameter estimation methods is effective in all cases and can overwhelm all the other methods. Instead, various methods have their advantages and disadvantages. Consequently, it is worthy to develop acceptably "good enough" methods within a given tolerance and time frame.
For practical purpose, some essential issues should be taken into account when developing a parameter estimation method: first, the method must be "efficient" enough that a trial can be completed within a reasonable computation time; second, for biological systems, the observation data is often corrupted by high level of noise, which complicates the objective function surface and introduces unwanted additional local minima in the search space [29]. Hence, the approach should be robust subject to noise; third, it needs to be flexible enough for adding/removing physical constraints, such as model dynamics, the mass and energy balance constraints. Furthermore, the representative cost function should have less local minima so as to ease the optimization algorithm in converging to the global minima. In this paper, two parameter estimation methods of combining spline theory [28] with Linear Programming (LP) and Nonlinear Programming (NLP) are developed, respectively. These methods remove the need for an ODE solver. Our analysis exhibits that the cost function surfaces of the two proposed methods are smooth.
Moreover, the cores of our algorithms are LP and NLP based, which are very flexible and hence additional constraints can be embeded/removed easily. Eight systems biology models were used to test the proposed algorithms. Experimental results show that the proposed methods are both efficient and robust (see additional file 1 for details).
This paper is organized as follows: The preliminary problem formulation is given and the bottleneck of the problem is highlighted in the next section. Then, two parameter estimation methods surmounting those bottlenecks are presented. In section 3, two trials are given, a simple enzyme kinetic model and the mammalian G1/S transition network model, in order to illustrate the robustness and the effectiveness of these two proposed methods (more models and trial results are given in additional file 1 and 2). Finally, conclusions and discussions are given in section 4.

Parameter estimation problem of systems biology models
Biological pathway dynamics can be modelled by the following continuous ODEs: ) ( ), where x R n is the system's state vector (for example the concentrations of a process), θ R k is the system's parameter vector (for instance, the reaction rates), u(t) R p is system's input, y R m denotes the measured data subject to a Gaussian white noise h(t)~N(0, σ 2 ), and x 0 is the initial state. f(·) is a set of nonlinear transition functions describing the dynamical properties of a biological system. Here, g(·) represents a measurement function. If all the states can be measured, the observer g(·) becomes an identity matrix. Otherwise, g(·) usually is a rectangular zero-one matrix with corresponding rows deleted (represent the immeasurable states) from the identity matrix I n .
The parameter estimation problem of nonlinear dynamical systems described in (1) can be formulated as a nonlinear programming problem (NLP) P 0 with differential-algebraic constraints: P 0 minimizes a cost function that measures the fitness of the model with respect to a given set of experiment data subjecting to a set of constraints, where ∈ R k is the set of parameters to be estimated, ||·|| l denotes the l-norm with l > 0,x 0 is the estimated initial condition,x R k ∈ is the estimated system states (ˆ( |ˆ) x t j  represents the estimated variable at time t j with parameter and initial conditionx 0 ), w ij are the weighting coefficients,ŷ is the estimated measured data. In some applications, additional constraints are introduced to impose special structural properties of a given system; they can be implemented in the form of the equality and inequality constraints C eq and C ineq (for instance the system performance and the mass balance constraints). Finally, θ L and θ U are simple structural constraints such as the parameter's upper/lower bounds (they can be part of the C ineq ).
For the NLP-P 0 , the direct optimization methods, such as Newton type methods and many GO methods, require solving the nonlinear dynamic model (1) forx in order to compute the cost function. The common method to estimate ( ) x t i andˆ( ) x t i is using ODE solvers, which perform the numerical integration with fixed at each iteration [19]. During the process of identification, the integration has to be executed thousands, even millions of times. That is the main reason more than 90% of the time is consumed in the ODE solver [24] and the computation time spent on the P 0 can be hours even days [10,20]. Moreover, P 0 is a nonlinear optimization problem subjecting to a set of linear and non-linear differential equation constraints. Hence, P 0 is often multimodal (non-convex) and has many local minima. In a highnoise environment, the situation becomes more difficult. Consequently, P 0 requires further manipulation in order to reduce the complexity so as to relieve the computation burden and also to avoid being trapped in local minima.
Instead of using ODE solvers to estimate x(t) and  x t ( ) , one can utilize spline approximation. Given L real values τ i , called knots, with τ 0 ≤ τ 1 ≤ · · · ≤ τ L-1 . Using the Cox-de Boor recursion formula, the B-spline basis of degree n d = 0, 1, 2, · · ·, L -2 can be defined as follows: , , , of length L in d -2, be the chosen basis functions. Then, the estimated variablex can be expressed in terms of the basis function expansion [28] ( ) Similarly, the estimated  ( ) x t i can be approximated by where , , , is the set of the derivatives of the basis functions. There are various types of splines suitable for this application, such as cubic spline, B-spline, uniform spline, nonuniform spline and interpolating spline. For more detail information about spline approximation theory, please refer to chapter (IX, XI, and XIV) in [28]. As B-spline is simple in formation and efficient for computation, it is adopted here. Our extensive tests have shown that uniform B-spline basis with produces good results. Hence, in this paper, unless otherwise indicated, the uniform B-spline basis with L N i ≈ 3 was used in the parameter identification process. Next, two techniques based on spline for parameter estimation will be proposed: one is based on linear programming (LP) which is very efficient and can cover many special structured systems and the other one is based on NLP which is flexible and can cater for general system structures.

The LP Approach
In many bio-system models, f (x, θ) is autonomous system and linear in θ as follows: where Φ(x) R n ×k is a matrix and its elements are a function of the state x. Systems with structure (8) covers a large set of systems biology models, such as enzyme kinetic pathway model, RKIP pathway model, IB-NF-B model TNFa-Mediated NF-B-signaling pathway model, irreversible inhibition of HIV proteinase model, Laub and Loomis model [2][3][4]30]. In addition, these types of models are usually subject to the mass balance constraints which can be incorporated into the LP easily (It is demonstrated in the results section via the Enzyme kinetic model).
For noisy data, good smoothing approximation can be achieved by minimizing the following cost function [26] whereˆ( ) x i m represents the mth derivative ofx i , m Z + and λ ≥ 0 control the trade-off between the "roughness" of the solution and hence can be used to relieve the overfitting problem. If the equality C eq represents the mass balance constraint and the inequality constraint C ineq represents parameter values' lower/upper bounds, the B-coefficient vector p = [p 1 ; p 2 ; · ·· p n ] T can be computed by solving the following quadratic programming sub-problem A 1 : ·ˆ= stands for the equality constraints C x eq ( , ) = 0 . A eq is a constant matrix, b eq is a constant vector. It is found empirically that m = 2 and 0 ≤ l ≤ 0.05 produce relatively good results. Hence, the parameters m = 2, and l = {0, 0.01, 0.03} corresponding to a noise level {0%, 5%, 10%}, respectively, were used in this study. Then, the "smooth" estimated statex can be generated by the B-spline approximation (5).
Replace x(t) by the estimated stateˆ( ) x t and integrate (8) yield: represents the transition matrix. Then, P 0 can be reformulated into the following optimization problem: Here, w ij ≥ 0 is weighting factor. Note that the L 1norm of a variable × is equivalent to the following rela- Then P 1 can be transformed into the following augmented optimization problem by introducing the slack variables a as follows: It is a Linear Programming (LP) problem with variable {a, θ}, which is a convex problem with a wealth of fast and efficient routines available [31,32].

Combine spline theory and NLP
To deal with systems biology models, of which the states and parameters are separable, the LP approach is suitable and efficient. In contrast, if the model does not belong to this category, such as the mammalian G1/S transition model and S-system model, the aforementioned approach cannot apply. Thus, a more general approach will be introduced. Recalling (6) and (7), the estimation ofˆ( ) x t j and ( ) x t j can be constructed by a set of basis functions. We can replaceˆ( ) x t j and ( ) x t j in (i-iv) of P 0 with (6) and (7). With little change, P 0 can be reformulated as . .
Note that the constraint-(i) of (2) has been replaced by constraints (i)-(iii) of P 3 . Then, NLP-P 0 -(2) with differential-algebraic constraints turns into NLP-P 3 -(14) with only algebraic equation constraints. Hence, P 3 does not require ODE solvers, which eases the computation burden (as shown in Examples). In contrast to the the decomposition methods [12,15], which divide the estimation of the system states (and its derivative) and the parameter estimation into two separate steps, P 3 computes the estimated states (and its derivative) and parameter values at the same time. Note that constraint (iii) of P 3 governs the estimated state (and its derivative) so as to ensure these estimates belong to the trajectorŷ ( |ˆ) x t  , which is a solution of system (1). Thus, the system model itself serves as a filter performing regularization. Hence, the overfitting problem can be relieved (see additional file 2 for details). For a non-linear system, the Lagrangian of (2), . However, many traditional optimization algorithms require the derivative ∂ ∂ L / during the optimization process. As L x ( , ) 0 is an implicit function of , ∂ ∂ L / cannot be obtained directly, but has to be computed via approximation methods [16], which makes the algorithm unreliable. For P 3 , the Lagrangian L p ( , )  now consists of simple algebraic constraints. Thus, ∂ ∂ L / and ∂L/ ∂p are explicit functions of and p. In conclusion, many of the aforementioned difficulties can be reduced. P 3 can be solved by a number of optimization approaches; either via evolution type algorithms, such as genetic algorithm (GA), simulated annealing (SA) and etc, or via traditional NLP algorithms, such as sequential quadratic programming(SQP), sequential penalty function, the trust region approach and etc [33,34].

Results
Two biological system models, a simple enzyme kinetic model and the mammalian G1/S transition network model, are chosen as benchmarks for evaluating the performance of P 2 and P 3 respectively.

Enzyme kinetic model
Consider the well-known simplified enzyme kinetic model. E is the concentration of an enzyme that combines with a substrate S to form an enzyme-substrate complex ES with a rate constant k 1 . The complex ES holds two possible out comes in the next step. It can be dissociated into E and S with a rate constant k 2 , or it can further proceed to form a product P with a rate constant k 3 . It is assumed that none of the products reverts to the initial substrate. These relations can be represented by the following set of ODEs.
where k 1 , k 2 , k 3 , are the system unknown parameters. Let x 1 (t), x 2 (t), x 3 (t)and x 4 (t) represent S(t),E(t), ES(t) and P(t) respectively. Then, the above equation can be rewritten into Then, the mass balance constraint becomes: x t x t x t x t x t x t x t x t x t x Or in matrix form, we have where According to (16), we have An artificial data set with four time courses was created. A total of 40 sampling points were assigned on each time courses. The observation data was perturbed by a zero mean Gaussian white noise h(t)~N(0, σ 2 ) in order to simulate the observation error. The estimated statex was computed by solving the quadratic programming sub-problem A 1 . The unknown parameter values were estimated using P 2 . The searching region of the parameter values was [0, +∞). All the computations were performed on a Pentium Dual Core computer (2.13 GHz ×2) with 2 GB RAM. The algorithm was implemented with Matlab-7 using the interior point algorithm. To quantify the fitness of the estimated model, the following relative squared error (RSE) measure J is employed: whereˆ( ) x t i j is the estimated time-course at time t j of a state variable x i , and x i (t j ) represents the "true" timecourses without noise at time t j . Note that smaller RSE J reflects better estimation. In order to obtain a statistical result on the quality of the estimation, 5,000 trials were performed. At each trial an estimated is computed using P 2 . Then, the mean estimation and standard deviation were deduced. The computation was very efficient and only took a few seconds for one estimation trial. Table 1 shows the statistical results of the estimation. It reveals that all the system parameter values are estimated successfully with a relative error of around 0.01% in noise free condition. When the system is subjected to a 10% observation noise level, all the mean estimated parameter values are within a relative tolerance, better than 2%. The RSE between the time-courses, produced by inferred model, and the given time-series data, averaged smaller than 1% even subject to 10% observation noise level condition. For this reason, we have confidence that the proposed method is robust within ±10% noise ratio (more complicated models and trials on P 2 can be found in part I of additional file 1). Figure 1 shows the responses of one trial.

The mammalian G1/S transition network model
Next, the mammalian G1/S transition network model, which includes a set of proteins and regulatory gene network, is used to test P 3 . In the mammalian G1/S transition network, pRB and AP-1 are the tumor suppressor from the family of pocket proteins and the family of transcription factors that mediate mitogenic signals, E2F1 is the transcription factor targeting genes that regulate cell cycle progression, Cyclin D/cdk4,6, cyclin E/cdk2, complexes characterizing the G1-and S-phases. There are various positive and negative feedback loops in the network controlling the G1/S transition. The positive feedback regulation of E2F1 and a double activator-inhibitor module can lead to bistability. The double activator-inhibitor module of the antagonistic plays E2F/DP on pRB make up the key unit of this phase transition. The graph representation of the mammalian G1/S transition network model can be found in additional file 1 and more details can refer to Swat et al. [5]. Definition of Variables for G1/S Transition Model is shown in Table 2. The corresponding ODE model is as follows where x is the set of state variables. There are totally 9 states and 39 parameters. The nominal parameter values are shown in Table 3.
Here, P 3 was solved by the Stochastic Raking Evolution Strategy (SRES) algorithm [35]. The searching region of the parameters was [0, 50θ ]. SRES uses stochastic ranking as the constraint handling technique, which adjusts the balance between the objective and penalty functions automatically during the evolutionary search. The observation data include 4 sets of time course, which consists of 40 sample points. For trials with noise free data, the algorithm converged in 8~9 hours after 250,000~300,000 iterations. The estimated parameter values, as shown in Table 3 are almost identical to the nominal parameter values. However, for k 23 , k 25 and J 15 , the estimated values are far from the nominal values, but the RSE measure is almost zero, which possibly implies that the system is insensitive with the changes of k 23 , k 25 and J 15 . This phenomenon reveals that the G1/S transition model either has some parameters that are insensitive to the chosen observation, or they are non-identifiable parameters [36,37]. It is worth reaction kinetics of enzyme Figure 1 The dynamic profiles of a trial with observation data subject to 10% random noises. Solid lines represent the "true" time-series data without noise, dots represent the measured time-series data with added artificial noise, and diamonds represent the estimated time-series data produced by the model.

Table 2 Definition of Variables for G1/S Transition Model
Acronym pRB E2F1 CycD i CycD a AP -1 pRB p pRB p p CycE i CycE a mentioning that the this large computational effort is the consequence of the very tight convergence criteria, an almost equal good result can be reached within 200,000 generations in about 6.5 hours with the RSE measure J is smaller than 1%. Figure 2(a) shows the "true" time-series data without noise and computed dynamic time-series data from one identified model. When 10% random noises are added, the convergence time increased and the relative estimation errors between estimated parameters and nominal parameters increased with the increase of noise. However, the timeseries produced by the estimated model is very similar to the original data, namely the RSE J is still small. This phenomenon may imply that there is no need to estimated every parameters accurately to achieve a model with equivalent dynamical properties with a good degree of accuracy. As the simulation time is long, performing thousands of simulations as the first method in order to evaluate the mean and variance of estimated parameters is impractical. Thus, due to the lack of space, results of just a few selected trial are shown in Table 3 (more trial results can be found in additional file 1). Trials were performed using Matlab-7. The main reason to use Matlab is that it is a convenient environment to visualize all the information arising from the optimization runs of the solver, evaluate new algorithms and modify existing algorithms. In contrast to the convenience, it is worth mentioning that Matlab programs usually are one order of magnitude (10 times or more) slower than equivalent compiled Fortan or C codes [19]. This is the major drawbacks of carrying programs out with Matlab. However, even in this situation, the performance of the proposed methods is acceptable.
For fair comparison, we also used the SRES algorithm to solve the same parameter estimation problem in the same searching region, but using NLP-P 0 with differential algebraic constraints as cost function. In this condition, after running 1 day, the algorithm failed to produce a set of parameters that can produce reasonable simulation result. We further reduced the searching region to [0, 3θ ] and used noise free data, but the estimation result was still not good and the RSE J is larger than 10.
Here, we use the G1/S model to show the differences of the cost function surfaces between NLP-P 0 and NLP-P 3 : in this case, the cost function of P 0 is a highly irregular and complicated manifold with multiple local minima; the augmented cost function adopted in problem P 3 is a much "smoother" function and hence it is easier for the NLP algorithm to converge to the solution. In order to simplify the analysis for exposition purpose, we only vary parameters k 1 and k 2 over the range  Figure 3(a) shows that the cost function surface of is a ridge, which drops suddenly from 10 9 to 0. However, Figure 3(b) reveals the cost function surface of P 0 are actually banana-shaped valley around the nominal value of the fixed parameters, this unfavorable profile can slow down the convergence rate of the algorithm. Furthermore, there are many local minima in the banana-shaped valley. Some algorithms, such as simulated annealing, genetic algorithm, have been proposed to overcome these problem. However, these algorithms are all computationally demanding. In conclusion, these cost function features make the problem P 0 a severe challenge to every optimization algorithm.
With the same condition, Figure 4(a) displays the cost function surface of P 3 , while Figure 4 Figure 2 The dynamic profiles of two trials. Solid lines represent the "true" time-series data without noise, dots represent the measured timeseries data with added artificial noise, and diamonds represent the estimated time-series data produced by the model: (a) noise free condition (b) 10% random noise condition. bowl-shaped, which is smoother. Similar results has also been observed when other combination of parameters served as variables. Obviously, if all 39 parameters vary at the same time, the surface of the cost function will be more "uneven" and more complicated. However, in this case, from the previous observations, the cost function surface of P 3 is smoother than the cost function surface of P 0. Furthermore, P 3 only involves algebraic equations as objective function and constraints. These properties make the NLP-P 3 easier to solve.

Discussion and Conclusion
In this paper, two parameter estimation methods based on spline theory are proposed. One aims at a narrower class of systems which is linear in parameters; however, it can cover many commonly found biological systems. The benefit is that the estimation problem can be transformed in an LP sub-algorithm which are fast and robust. Additional linear constraints can be embedded relative easily. For general systems, the problem is solved by an NLP with algebraic constraints, which is more computationally demanding.
A simple enzyme kinetic model and the mammalian G1/S transition network model were used as benchmarks to evaluate the performance of the two proposed methods. We illustrate the usefulness with more examples in additional files 1 but these do not remotely cover all the conditions. During the simulation of the mammalian G1/S transition network model, we found that the estimated parameter set Φ A ≡ {k 1 , k 2 , k p , J 11 , J 12 , K m1 , K m2 , j 1 , j 2 } were well within the respective nominal values. While the set Φ B ≡ {J 61 , J 62 , J 63 , J 65 , J 68 } were far from their nominal values. However, the time-series produced by the estimated model were very similar to the original data. This phenomena reveals that some parameter values are insensitive in the searching region. Interestingly, we find that the "sensitive" or "easily identified" parameters set Φ A are also the parameters of the double-activator-inhibitor module of the antagonistic players E2F/DP and pRB, which makes up the core unit of the G1/S transition model [5]. This phenomenon may imply that the parameter values of the core module are sensitive and easy to identify. In contrast, the parameters set Φ B seems to be insensitive, which may reflect that pRG p (x 6 ) is not a key element of the total system. However, to identify which parameter values or variables are important, a sensitivity analysis is needed [38], which is another important topic in systems biology and deserves a more detailed study. This sensitivity analysis is a preprocess for isolating those states and parameters which are sensitive in order to reduce the dimension of the system model and to improve the numerical stability for the core estimation problem.
For most biological systems, the ODE models are often high-dimensional and nonlinear. The problem of system parameter estimation is computationally expensive and can easily be trapped in local minima. We find that under noisy conditions, it is almost impossible to accurately estimate every parameter of the sloppy biological system model. However, in practice, a model with equivalent dynamical properties with a good degree of accuracy can be constructed based on dominant sensitive parameters and system states. The following are some of practical observations: 1. High quality experiment data is essential for identifying accurate biology systems. When the experiment data is corrupted with high level noise, it needs more experimental data. If an insufficient amount of time-series data is given as observed profiles, the high degree-freedom of systems biology models ensures that many candidate solutions will be found. 2. Perform a sensitivity analysis and identifiability analysis before the identification phase [36][37][38]. 3. For systems models with insensitive or non-identifiable parameters, the search may lead to a solution where some parameters can have large deviation, but still produce satisfactory system responses. This problem can be partly relieved by introducing auxiliary information (additional constraints such as shrinking the searching region) of the model into the algorithm. However, it remains difficult to be solved completely by improving parameter estimation strategy. It indicates that researchers should focus on predictions rather than on accurately estimating every parameter.
Although the proposed algorithms are fast and robust, there is certainly room for improvement: for method 1, it is not general enough to catch every case; for method 2, the price for the simplicity and generality is at the expansion of the optimization variable dimension. Under high noise condition, method 2 is still not robust enough. At the moment, the testing is based on Matlab which is much slower than native codes produced by C, Fortran, etc, however the conversion is straight forward. Currently, many high-speed computation engines are available that make use of parallelism, for instance multi-cluster engines, array-processing engines etc. Hence, one possible way is developed algorithm on these high-speed computation engines environment. Another possible way is developing hybrid algorithms to incorporate elements from evolution algorithms such as GA, SA and PSO. In this paper, we have considered the parameter estimation problem with known structure. However, it is easy to expand our method to structure identification by introducing an additional penalty term to the objective function [39].