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  <cit> .
one common approach is to characterize the biological system with a set of ordinary differential equations   <cit> .
generally, there are two major aspects of building an ode model for a biological system from experimentally measured time series:  to determine the structure of the system through a set of suitable odes with unknown parameters;  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  <cit> .
this type of approach can be considered as a "reverse engineering process"  <cit> .
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  <cit> .
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  method can offer promising results  <cit> .
in  <cit> , many deterministic and stochastic global optimization  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  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  <cit> .
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  <cit> .
in particular, for nonlinear dynamical systems with high-parameter-dimension, one trial usually consumes tens of hours or even days  <cit> .
furthermore, the convergence property is aggravated by numerical integration failure, which is a major problem in the optimization process  <cit> .
the computational burden can be relieved by reformulating the system involving differential equations into a system of algebraic equations  <cit> , which can be classified as "decomposition approaches".
these decomposition approaches are widely employed in the parameter estimation of s-systems  <cit> .
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 pre-processing, the estimation faces the potential of the overfitting problem and hence the estimation can deviate badly from the "true" value  <cit> .
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  <cit> .
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  <cit> .
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  <cit> .
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  <cit>  with linear programming  and nonlinear programming  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 .
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  <dig>  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 .
finally, conclusions and discussions are given in section  <dig>
methods
parameter estimation problem of systems biology models
biological pathway dynamics can be modelled by the following continuous odes:
x˙=f,u,θ),   x=x <dig> y=g)+η,
where x ϵ rn is the system's state vector , θ ϵ rk is the system's parameter vector , u ϵ rp is system's input, y ϵ rm denotes the measured data subject to a gaussian white noise η ~ n, and x <dig> 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  from the identity matrix in.
the parameter estimation problem of nonlinear dynamical systems described in  can be formulated as a nonlinear programming problem  p <dig> with differential-algebraic constraints:
p0:minθ^,x^0∑j=0n−1∑i=1nwij‖yi−y^i‖l,s.t.
{x˙^=f,u,θ^), x=x^ <dig> y^=g),  j= <dig> ,⋯,n− <dig> ceq= <dig> cineq≤ <dig> θl≤θ^≤θu.
p <dig> 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 θ^∈rk is the set of parameters to be estimated, ||·||l denotes the l-norm with l >  <dig>  x^ <dig> is the estimated initial condition, x^∈rk is the estimated system states  represents the estimated variable at time tj with parameter θ^ and initial condition x^0), wij are the weighting coefficients, y^ 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 ceq and cineq .
finally, θ l and θ u are simple structural constraints such as the parameter's upper/lower bounds .
for the nlp-p <dig> , the direct optimization methods, such as newton type methods and many go methods, require solving the nonlinear dynamic model  for x^ in order to compute the cost function.
the common method to estimate x˙^ and x^ is using ode solvers, which perform the numerical integration with θ^ fixed at each iteration  <cit> .
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  <cit>  and the computation time spent on the p <dig> can be hours even days  <cit> .
moreover, p <dig> is a nonlinear optimization problem subjecting to a set of linear and non-linear differential equation constraints.
hence, p <dig> is often multimodal  and has many local minima.
in a high-noise environment, the situation becomes more difficult.
consequently, p <dig> 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 and x˙, one can utilize spline approximation.
given l real values τi, called knots, with τ <dig> ≤ τ <dig> ≤ · · · ≤ τl- <dig>  using the cox-de boor recursion formula, the b-spline basis of degree nd =  <dig>   <dig>   <dig>  · · ·, l -  <dig> can be defined as follows:
bm,0≜{1ifτm≤t<τm+ <dig> m= <dig> ,⋯,l− <dig> otherwise,
bm,nd≜t−τmτm+nd−τmbm,nd−1+τm+nd+1−tτm+nd+1−τm+1bm+ <dig> nd− <dig>    m= <dig> , ⋯, l−nd− <dig>
let bi=t, a vector of length li - nd -  <dig>  be the chosen basis functions.
then, the estimated variable x^ can be expressed in terms of the basis function expansion  <cit>
x^i=∑m=1li−nb−2pi,m·bi,nb,   t∈
where x^i∈r is the estimation of the ith state of , pi, m is the weighting coefficient.
let pi=t,  can be rewritten in matrix form
x^i=bit·pi,    i= <dig> , ⋯, n.
similarly, the estimated x˙^i can be approximated by
x˙^i=b˙it·pi,
where b˙j=t 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  in  <cit> .
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 n3≤li≤n <dig> produces good results.
hence, in this paper, unless otherwise indicated, the uniform b-spline basis with li≈n <dig> 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  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  is autonomous system and linear in θ as follows:
x˙=Φ)θ,  x=x <dig>
where Φ ϵ rn ×k is a matrix and its elements are a function of the state x. systems with structure  covers a large set of systems biology models, such as enzyme kinetic pathway model, rkip pathway model, iκb-nf-κb model tnfα-mediated nf-κb-signaling pathway model, irreversible inhibition of hiv proteinase model, laub and loomis model  <cit> .
in addition, these types of models are usually subject to the mass balance constraints which can be incorporated into the lp easily .
for noisy data, good smoothing approximation can be achieved by minimizing the following cost function  <cit>
ci=1n+1∑j=0n−x^i)2+λ∑j=0n) <dig>
where x^i represents the mth derivative of x^i, m ϵ z+ and λ ≥  <dig> control the trade-off between the "roughness" of the solution and hence can be used to relieve the overfitting problem.
if the equality ceq represents the mass balance constraint and the inequality constraint cineq represents parameter values' lower/upper bounds, the b-coefficient vector p = t can be computed by solving the following quadratic programming sub-problem a <dig> :
a1: minp{1n+1∑j=0n||y−x^||22+λ∑j=0n||x^||22}s.t.
{aeq⋅x^=beq,x^i=bit⋅pi,x^i=)t⋅pi,i,= <dig> ,⋯n,j= <dig> ,2⋯n.
aeq·x^=beq stands for the equality constraints ceq= <dig>  aeq is a constant matrix, beq is a constant vector.
it is found empirically that m =  <dig> and  <dig> ≤ λ ≤  <dig>  produce relatively good results.
hence, the parameters m =  <dig>  and λ = { <dig>   <dig> ,  <dig> } corresponding to a noise level {0%, 5%, 10%}, respectively, were used in this study.
then, the "smooth" estimated state x^ can be generated by the b-spline approximation .
replace x by the estimated state x^ and integrate  yield:
x˜=)dt)  ·  θ^+x^0=Ψ^jθ^+x^ <dig>
where Ψ^j=∫t=t0tjΦ)dt represents the transition matrix.
then, p <dig> can be reformulated into the following optimization problem:
p1:minθ^∈rk∑j=0n∑i=1nwij|yi−x˜i|s.t.
{x˜−Ψ^jθ^+x^ <dig> θl≤θ^≤θu,i= <dig> ..,n; j= <dig> , <dig> ⋯,n.
here, wij ≥  <dig> is weighting factor.
note that the l1-norm of a variable × is equivalent to the following relation: |x|=minα≥0{α}; s.t.
- α ≤ x ≤ α. then p <dig> can be transformed into the following augmented optimization problem by introducing the slack variables α as follows:
p2:minθ^,α{∑i,jwij⋅αij}s.t.
{−αij≤yi−x˜i,αij≥yi−x˜i,x˜=Ψ^jθ^+x^0αij≥ <dig> θl≤θ^≤θu,i= <dig> ..,n; j= <dig> ,...,n.
it is a linear programming  problem with variable {α, θ}, which is a convex problem with a wealth of fast and efficient routines available  <cit> .
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  and , the estimation of x^ and x˙^ can be constructed by a set of basis functions.
we can replace x^ and x˙^ in  of p <dig> with  and .
with little change, p <dig> can be reformulated as
p3:minθ^,p∑j=0n−y^)twj−y^)s.t.
{x^i=bit⋅pi,x˙^i=bit⋅pi,‖x˙^−f,u,θ^)‖22= <dig> y^=g)ceq,x˙^,θ^)= <dig> cineq,x˙^,θ^)≤ <dig> θl≤θ^≤θu,i= <dig> ,⋯,n, j= <dig> ,⋯,n.
note that the constraint- of  has been replaced by constraints - of p <dig>  then, nlp-p0- with differential-algebraic constraints turns into nlp-p3- with only algebraic equation constraints.
hence, p <dig> does not require ode solvers, which eases the computation burden .
in contrast to the the decomposition methods  <cit> , which divide the estimation of the system states  and the parameter estimation into two separate steps, p <dig> computes the estimated states  and parameter values at the same time.
note that constraint  of p <dig> governs the estimated state  so as to ensure these estimates belong to the trajectory x^, which is a solution of system .
thus, the system model itself serves as a filter performing regularization.
hence, the overfitting problem can be relieved .
for a non-linear system, the lagrangian of , l, is an implicit function of {θ^,x^0} <cit> .
however, many traditional optimization algorithms require the derivative ∂l/∂θ^ during the optimization process.
as l is an implicit function of θ^, ∂l/∂θ^ cannot be obtained directly, but has to be computed via approximation methods  <cit> , which makes the algorithm unreliable.
for p <dig>  the lagrangian l 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 <dig> can be solved by a number of optimization approaches; either via evolution type algorithms, such as genetic algorithm , simulated annealing  and etc, or via traditional nlp algorithms, such as sequential quadratic programming, sequential penalty function, the trust region approach and etc  <cit> .
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 <dig> and p <dig> 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 <dig>  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 <dig>  or it can further proceed to form a product p with a rate constant k <dig>  it is assumed that none of the products reverts to the initial substrate.
these relations can be represented by the following set of odes.
dsdt=−k1⋅e⋅s+k2⋅es,dedt=−k1·e·s+.
es ,desdt=k1⋅ e ⋅s− ⋅es dpdt=k3⋅ es ,
where k <dig>  k <dig>  k <dig>  are the system unknown parameters.
let x <dig>  x <dig>  x3and x <dig> represent s,e, es and p respectively.
then, the above equation can be rewritten into
x˙1=−k1x1x2+k2x <dig> x˙2=−k1x1x2+k2x3+k3x <dig> x˙3=k1x1x2−k2x3−k3x <dig> x˙4=k3x <dig> ⇒  =  ⋅
then, the mass balance constraint becomes:
{x2+x3=x2+x3x1+x3+x4=x1+x3+x <dig>
or in matrix form, we have
aeq·x=beq
where aeq= <cit>   and beq=.
according to , we have
Φ)= , θ= ,and  Ψ^)=∫t=t0ti dt.
an artificial data set with four time courses was created.
a total of  <dig> sampling points were assigned on each time courses.
the observation data was perturbed by a zero mean gaussian white noise η ~ n in order to simulate the observation error.
the estimated state x^ was computed by solving the quadratic programming sub-problem a <dig>  the unknown parameter values were estimated using p <dig>  the searching region of the parameter values was [ <dig>  +∞).
all the computations were performed on a pentium dual core computer  with  <dig> gb ram.
the algorithm was implemented with matlab- <dig> using the interior point algorithm.
to quantify the fitness of the estimated model, the following relative squared error  measure j is employed:
j=1n⋅n∑i=1n∑j=0n−xixi) <dig>
where x^i is the estimated time-course at time tj of a state variable xi, and xi represents the "true" time-courses without noise at time tj.
note that smaller rse j reflects better estimation.
in order to obtain a statistical result on the quality of the estimation,  <dig>  trials were performed.
at each trial an estimated θ^ is computed using p <dig>  then, the mean estimation and standard deviation were deduced.
the computation was very efficient and only took a few seconds for one estimation 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 <dig>  in the mammalian g1/s transition network, prb and ap- <dig> are the tumor suppressor from the family of pocket proteins and the family of transcription factors that mediate mitogenic signals, e2f <dig> is the transcription factor targeting genes that regulate cell cycle progression, cyclin d/cdk <dig> , cyclin e/cdk <dig>  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 e2f <dig> 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  <dig> and more details can refer to swat et al.
<cit> .
definition of variables for g1/s transition model is shown in table  <dig>  the corresponding ode model is as follows
x
1
x
2
x
3
x
4
x
5
x
6
x
7
x
8
x
9
x˙1=k1x2km1+x2j11j11+x1j61j61+x6−k16x1x4+k61x6−ϕ1x <dig> x˙2=kp+k2a22x22km22+x22j12j12+x1j62j62+x6−ϕ2x <dig> x˙3=k3x5+k23x2j13j13+x1j63j63+x6+k43x4−k34x3x4km4+x4          −ϕ3x <dig> x˙4=k34x3x4km4+x4−k43x4−ϕ4x <dig> x˙5=fm+k25x2j15j15+x1j65j65+x6−ϕ5x <dig> x˙6=k16x1x4−k61x6−k67x6x9+k76x7−ϕ6x <dig> x˙7=k67x6x9−k76x7−ϕ7x <dig> x˙8=k28x2j18j18+x1j68j68+x6+k98x9−k89x8x9km9+x9−ϕ8x <dig> x˙9=k89x8x9km9+x9−k98x9−ϕ9x <dig>
where x is the set of state variables.
there are totally  <dig> states and  <dig> parameters.
the nominal parameter values are shown in table  <dig>
here, p <dig> was solved by the stochastic raking evolution strategy  algorithm  <cit> .
the searching region of the parameters was .
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  <dig> sets of time course, which consists of  <dig> sample points.
for trials with noise free data, the algorithm converged in  <dig> ~  <dig> hours after  <dig>  ~  <dig>  iterations.
the estimated parameter values, as shown in table  <dig> are almost identical to the nominal parameter values.
however, for k <dig>  k <dig> and j <dig>  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 <dig>  k <dig> and j <dig>  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  <cit> .
it is worth 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  <dig>  generations in about  <dig>  hours with the rse measure j is smaller than 1%.
figure  <dig> 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 time-series 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  <dig> .
trials were performed using matlab- <dig>  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  slower than equivalent compiled fortan or c codes  <cit> .
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 <dig> with differential algebraic constraints as cost function.
in this condition, after running  <dig> day, the algorithm failed to produce a set of parameters that can produce reasonable simulation result.
we further reduced the searching region to  and used noise free data, but the estimation result was still not good and the rse j is larger than  <dig>
here, we use the g1/s model to show the differences of the cost function surfaces between nlp-p <dig> and nlp-p3: in this case, the cost function of p <dig> is a highly irregular and complicated manifold with multiple local minima; the augmented cost function adopted in problem p <dig> 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 <dig> and k <dig> over the range k <dig> ϵ  and k <dig> ϵ  and fix all other parameters at their nominal values.
figure  <dig> displays the cost function surface of p <dig> , while figure  <dig> exhibits the same data as figure  <dig> on the expanded scale and figure  <dig> is the corresponding contour plots.
figure  <dig> shows that the cost function surface of is a ridge, which drops suddenly from  <dig> to  <dig>  however, figure  <dig> reveals the cost function surface of p <dig> 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 <dig> a severe challenge to every optimization algorithm.
with the same condition, figure  <dig> displays the cost function surface of p <dig>  while figure  <dig> shows the corresponding contour line.
compared with the cost function surface of p <dig>  the cost function surface of p <dig> is bowl-shaped, which is smoother.
similar results has also been observed when other combination of parameters served as variables.
obviously, if all  <dig> 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 <dig> is smoother than the cost function surface of p <dig>
furthermore, p <dig> only involves algebraic equations as objective function and constraints.
these properties make the nlp-p <dig> easier to solve.
discussion and
CONCLUSIONS
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  <dig> 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 <dig>  k <dig>  kp, j <dig>  j <dig>  km <dig>  km <dig>  ϕ <dig>  ϕ2} were well within the respective nominal values.
while the set Φb ≡ {j <dig>  j <dig>  j <dig>  j <dig>  j68} 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  <cit> .
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 prgp is not a key element of the total system.
however, to identify which parameter values or variables are important, a sensitivity analysis is needed  <cit> , which is another important topic in systems biology and deserves a more detailed study.
this sensitivity analysis is a pre-process 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:
<dig>  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.
<dig>  perform a sensitivity analysis and identifiability analysis before the identification phase  <cit> .
<dig>  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  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  <dig>  it is not general enough to catch every case; for method  <dig>  the price for the simplicity and generality is at the expansion of the optimization variable dimension.
under high noise condition, method  <dig> 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  <cit> .
authors' contributions
all authors contributed to the development of the theory and analysis of the results.
cj designed, analyzed and implemented the algorithm, performed the experiments and wrote the main body of manuscript.
ly contributed to the development of the method and writing the manuscript and providing the critical review of the manuscript.
all authors have read and approved the final version of the manuscript.
supplementary material
additional file 1
in this additional file, we tested the proposed methods on seven systems biology models were used to test: tnfα -mediated nf-κb-signaling pathway model, rkip regulated erk pathway model and the model of irreversible inhibition of hiv proteinase; yeast fermentation pathway model, large-scale target genetic network model, a three step pathway model and the mammalian g1/s transition network model.
click here for file
additional file 2
in this additional file, we use e2f/dp dimmer model to illustrate the differences between the three different type methods mentioned in the paper:  the direct optimization method;  decomposition methods;  methods combine spline theory and nlp.
click here for file
acknowledgements
the work is supported by cerg grant  <dig> .
we would like to express our appreciation for dr. robin s. bradbeer and dr. k. t. ko for their fruitful suggestions.

 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. one common approach is to characterize the biological system with a set of ordinary differential equations. generally, there are two major aspects of building an ode model for a biological system from experimentally measured time series: to determine the structure of the system through a set of suitable odes with unknown parameters; 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. this type of approach can be considered as a "reverse engineering process". 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. 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 method can offer promising results. in, many deterministic and stochastic global optimization 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 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. 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. in particular, for nonlinear dynamical systems with high-parameter-dimension, one trial usually consumes tens of hours or even days. furthermore, the convergence property is aggravated by numerical integration failure, which is a major problem in the optimization process. the computational burden can be relieved by reformulating the system involving differential equations into a system of algebraic equations, which can be classified as "decomposition approaches". these decomposition approaches are widely employed in the parameter estimation of s-systems. 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 pre-processing, the estimation faces the potential of the overfitting problem and hence the estimation can deviate badly from the "true" value. 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. 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. 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. 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 with linear programming and nonlinear programming 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. 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 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. finally, conclusions and discussions are given in section methods parameter estimation problem of systems biology models biological pathway dynamics can be modelled by the following continuous odes: x˙=f,u,θ),   x=x y=g)+η, where x ϵ rn is the system's state vector, θ ϵ rk is the system's parameter vector, u ϵ rp is system's input, y ϵ rm denotes the measured data subject to a gaussian white noise η ~ n, and x 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 from the identity matrix in. the parameter estimation problem of nonlinear dynamical systems described in can be formulated as a nonlinear programming problem p with differential-algebraic constraints: p0:minθ^,x^0∑j=0n−1∑i=1nwij‖yi−y^i‖l,s.t. {x˙^=f,u,θ^), x=x^ y^=g),  j=,⋯,n− ceq= cineq≤ θl≤θ^≤θu. p 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 θ^∈rk is the set of parameters to be estimated, ||·||l denotes the l-norm with l > x^ is the estimated initial condition, x^∈rk is the estimated system states represents the estimated variable at time tj with parameter θ^ and initial condition x^0), wij are the weighting coefficients, y^ 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 ceq and cineq. finally, θ l and θ u are simple structural constraints such as the parameter's upper/lower bounds. for the nlp-p, the direct optimization methods, such as newton type methods and many go methods, require solving the nonlinear dynamic model for x^ in order to compute the cost function. the common method to estimate x˙^ and x^ is using ode solvers, which perform the numerical integration with θ^ fixed at each iteration. 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 and the computation time spent on the p can be hours even days. moreover, p is a nonlinear optimization problem subjecting to a set of linear and non-linear differential equation constraints. hence, p is often multimodal and has many local minima. in a high-noise environment, the situation becomes more difficult. consequently, p 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 and x˙, one can utilize spline approximation. given l real values τi, called knots, with τ ≤ τ ≤ · · · ≤ τl- using the cox-de boor recursion formula, the b-spline basis of degree nd = · · ·, l - can be defined as follows: bm,0≜{1ifτm≤t<τm+ m=,⋯,l− otherwise, bm,nd≜t−τmτm+nd−τmbm,nd−1+τm+nd+1−tτm+nd+1−τm+1bm+ nd−    m=, ⋯, l−nd− let bi=t, a vector of length li - nd - be the chosen basis functions. then, the estimated variable x^ can be expressed in terms of the basis function expansion x^i=∑m=1li−nb−2pi,m·bi,nb,   t∈ where x^i∈r is the estimation of the ith state of, pi, m is the weighting coefficient. let pi=t, can be rewritten in matrix form x^i=bit·pi,    i=, ⋯, n. similarly, the estimated x˙^i can be approximated by x˙^i=b˙it·pi, where b˙j=t 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 in. 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 n3≤li≤n produces good results. hence, in this paper, unless otherwise indicated, the uniform b-spline basis with li≈n 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 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 is autonomous system and linear in θ as follows: x˙=Φ)θ,  x=x where Φ ϵ rn ×k is a matrix and its elements are a function of the state x. systems with structure covers a large set of systems biology models, such as enzyme kinetic pathway model, rkip pathway model, iκb-nf-κb model tnfα-mediated nf-κb-signaling pathway model, irreversible inhibition of hiv proteinase model, laub and loomis model. in addition, these types of models are usually subject to the mass balance constraints which can be incorporated into the lp easily. for noisy data, good smoothing approximation can be achieved by minimizing the following cost function ci=1n+1∑j=0n−x^i)2+λ∑j=0n) where x^i represents the mth derivative of x^i, m ϵ z+ and λ ≥ control the trade-off between the "roughness" of the solution and hence can be used to relieve the overfitting problem. if the equality ceq represents the mass balance constraint and the inequality constraint cineq represents parameter values' lower/upper bounds, the b-coefficient vector p = t can be computed by solving the following quadratic programming sub-problem a: a1: minp{1n+1∑j=0n||y−x^||22+λ∑j=0n||x^||22}s.t. {aeq⋅x^=beq,x^i=bit⋅pi,x^i=)t⋅pi,i,=,⋯n,j=,2⋯n. aeq·x^=beq stands for the equality constraints ceq= aeq is a constant matrix, beq is a constant vector. it is found empirically that m = and ≤ λ ≤ produce relatively good results. hence, the parameters m = and λ = {, } corresponding to a noise level {0%, 5%, 10%}, respectively, were used in this study. then, the "smooth" estimated state x^ can be generated by the b-spline approximation. replace x by the estimated state x^ and integrate yield: x˜=)dt)  ·  θ^+x^0=Ψ^jθ^+x^ where Ψ^j=∫t=t0tjΦ)dt represents the transition matrix. then, p can be reformulated into the following optimization problem: p1:minθ^∈rk∑j=0n∑i=1nwij|yi−x˜i|s.t. {x˜−Ψ^jθ^+x^ θl≤θ^≤θu,i=..,n; j=, ⋯,n. here, wij ≥ is weighting factor. note that the l1-norm of a variable × is equivalent to the following relation: |x|=minα≥0{α}; s.t. - α ≤ x ≤ α. then p can be transformed into the following augmented optimization problem by introducing the slack variables α as follows: p2:minθ^,α{∑i,jwij⋅αij}s.t. {−αij≤yi−x˜i,αij≥yi−x˜i,x˜=Ψ^jθ^+x^0αij≥ θl≤θ^≤θu,i=..,n; j=,...,n. it is a linear programming problem with variable {α, θ}, which is a convex problem with a wealth of fast and efficient routines available. 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 and, the estimation of x^ and x˙^ can be constructed by a set of basis functions. we can replace x^ and x˙^ in of p with and. with little change, p can be reformulated as p3:minθ^,p∑j=0n−y^)twj−y^)s.t. {x^i=bit⋅pi,x˙^i=bit⋅pi,‖x˙^−f,u,θ^)‖22= y^=g)ceq,x˙^,θ^)= cineq,x˙^,θ^)≤ θl≤θ^≤θu,i=,⋯,n, j=,⋯,n. note that the constraint- of has been replaced by constraints - of p then, nlp-p0- with differential-algebraic constraints turns into nlp-p3- with only algebraic equation constraints. hence, p does not require ode solvers, which eases the computation burden. in contrast to the the decomposition methods, which divide the estimation of the system states and the parameter estimation into two separate steps, p computes the estimated states and parameter values at the same time. note that constraint of p governs the estimated state so as to ensure these estimates belong to the trajectory x^, which is a solution of system. thus, the system model itself serves as a filter performing regularization. hence, the overfitting problem can be relieved. for a non-linear system, the lagrangian of, l, is an implicit function of {θ^,x^0}. however, many traditional optimization algorithms require the derivative ∂l/∂θ^ during the optimization process. as l is an implicit function of θ^, ∂l/∂θ^ cannot be obtained directly, but has to be computed via approximation methods, which makes the algorithm unreliable. for p the lagrangian l 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 can be solved by a number of optimization approaches; either via evolution type algorithms, such as genetic algorithm, simulated annealing and etc, or via traditional nlp algorithms, such as sequential quadratic programming, sequential penalty function, the trust region approach and etc.  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 and p 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 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 or it can further proceed to form a product p with a rate constant k it is assumed that none of the products reverts to the initial substrate. these relations can be represented by the following set of odes. dsdt=−k1⋅e⋅s+k2⋅es,dedt=−k1·e·s+. es,desdt=k1⋅ e ⋅s− ⋅es dpdt=k3⋅ es, where k k k are the system unknown parameters. let x x x3and x represent s,e, es and p respectively. then, the above equation can be rewritten into x˙1=−k1x1x2+k2x x˙2=−k1x1x2+k2x3+k3x x˙3=k1x1x2−k2x3−k3x x˙4=k3x ⇒  =  ⋅ then, the mass balance constraint becomes: {x2+x3=x2+x3x1+x3+x4=x1+x3+x or in matrix form, we have aeq·x=beq where aeq=   and beq=. according to, we have Φ)=, θ=,and  Ψ^)=∫t=t0ti dt. an artificial data set with four time courses was created. a total of sampling points were assigned on each time courses. the observation data was perturbed by a zero mean gaussian white noise η ~ n in order to simulate the observation error. the estimated state x^ was computed by solving the quadratic programming sub-problem a the unknown parameter values were estimated using p the searching region of the parameter values was [ +∞). all the computations were performed on a pentium dual core computer with gb ram. the algorithm was implemented with matlab- using the interior point algorithm. to quantify the fitness of the estimated model, the following relative squared error measure j is employed: j=1n⋅n∑i=1n∑j=0n−xixi) where x^i is the estimated time-course at time tj of a state variable xi, and xi represents the "true" time-courses without noise at time tj. note that smaller rse j reflects better estimation. in order to obtain a statistical result on the quality of the estimation, trials were performed. at each trial an estimated θ^ is computed using p then, the mean estimation and standard deviation were deduced. the computation was very efficient and only took a few seconds for one estimation 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 in the mammalian g1/s transition network, prb and ap- are the tumor suppressor from the family of pocket proteins and the family of transcription factors that mediate mitogenic signals, e2f is the transcription factor targeting genes that regulate cell cycle progression, cyclin d/cdk, cyclin e/cdk 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 e2f 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 and more details can refer to swat et al.. definition of variables for g1/s transition model is shown in table the corresponding ode model is as follows x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 x˙1=k1x2km1+x2j11j11+x1j61j61+x6−k16x1x4+k61x6−ϕ1x x˙2=kp+k2a22x22km22+x22j12j12+x1j62j62+x6−ϕ2x x˙3=k3x5+k23x2j13j13+x1j63j63+x6+k43x4−k34x3x4km4+x4          −ϕ3x x˙4=k34x3x4km4+x4−k43x4−ϕ4x x˙5=fm+k25x2j15j15+x1j65j65+x6−ϕ5x x˙6=k16x1x4−k61x6−k67x6x9+k76x7−ϕ6x x˙7=k67x6x9−k76x7−ϕ7x x˙8=k28x2j18j18+x1j68j68+x6+k98x9−k89x8x9km9+x9−ϕ8x x˙9=k89x8x9km9+x9−k98x9−ϕ9x where x is the set of state variables. there are totally states and parameters. the nominal parameter values are shown in table here, p was solved by the stochastic raking evolution strategy algorithm. the searching region of the parameters was. 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 sets of time course, which consists of sample points. for trials with noise free data, the algorithm converged in ~ hours after ~ iterations. the estimated parameter values, as shown in table are almost identical to the nominal parameter values. however, for k k and j 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 k and j 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. it is worth 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 generations in about hours with the rse measure j is smaller than 1%. figure 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 time-series 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. trials were performed using matlab- 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 slower than equivalent compiled fortan or c codes. 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 with differential algebraic constraints as cost function. in this condition, after running day, the algorithm failed to produce a set of parameters that can produce reasonable simulation result. we further reduced the searching region to and used noise free data, but the estimation result was still not good and the rse j is larger than here, we use the g1/s model to show the differences of the cost function surfaces between nlp-p and nlp-p3: in this case, the cost function of p is a highly irregular and complicated manifold with multiple local minima; the augmented cost function adopted in problem p 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 and k over the range k ϵ and k ϵ and fix all other parameters at their nominal values. figure displays the cost function surface of p, while figure exhibits the same data as figure on the expanded scale and figure is the corresponding contour plots. figure shows that the cost function surface of is a ridge, which drops suddenly from to however, figure reveals the cost function surface of p 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 a severe challenge to every optimization algorithm. with the same condition, figure displays the cost function surface of p while figure shows the corresponding contour line. compared with the cost function surface of p the cost function surface of p is bowl-shaped, which is smoother. similar results has also been observed when other combination of parameters served as variables. obviously, if all 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 is smoother than the cost function surface of p furthermore, p only involves algebraic equations as objective function and constraints. these properties make the nlp-p easier to solve. discussion and  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 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 k kp, j j km km ϕ ϕ2} were well within the respective nominal values. while the set Φb ≡ {j j j j j68} 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. 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 prgp is not a key element of the total system. however, to identify which parameter values or variables are important, a sensitivity analysis is needed, which is another important topic in systems biology and deserves a more detailed study. this sensitivity analysis is a pre-process 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: 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. perform a sensitivity analysis and identifiability analysis before the identification phase. 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 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 it is not general enough to catch every case; for method the price for the simplicity and generality is at the expansion of the optimization variable dimension. under high noise condition, method 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. authors' contributions all authors contributed to the development of the theory and analysis of the results. cj designed, analyzed and implemented the algorithm, performed the experiments and wrote the main body of manuscript. ly contributed to the development of the method and writing the manuscript and providing the critical review of the manuscript. all authors have read and approved the final version of the manuscript. supplementary material additional file 1 in this additional file, we tested the proposed methods on seven systems biology models were used to test: tnfα -mediated nf-κb-signaling pathway model, rkip regulated erk pathway model and the model of irreversible inhibition of hiv proteinase; yeast fermentation pathway model, large-scale target genetic network model, a three step pathway model and the mammalian g1/s transition network model. click here for file additional file 2 in this additional file, we use e2f/dp dimmer model to illustrate the differences between the three different type methods mentioned in the paper: the direct optimization method; decomposition methods; methods combine spline theory and nlp. click here for file acknowledgements the work is supported by cerg grant. we would like to express our appreciation for dr. robin s. bradbeer and dr. k. t. ko for their fruitful suggestions.
many parameter estimation methods have been developed but none of these methods is effective for all cases and can overwhelm all other approaches. it is worth to develop parameter estimation methods which are robust against noise, efficient in computation and flexible enough to meet different constraints. two parameter estimation methods of combining spline theory with linear programming and nonlinear programming 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. instead, various methods have their advantages and disadvantages. the proposed approaches have general application to identify unknown parameter values of a wide range of systems biology models. moreover, the cores of our algorithms are lp and nlp based, which are flexible and consequently additional constraints can be embedded/removed easily. our results confirm that the proposed methods are both efficient and robust. mathematical models for revealing the dynamics and interactions properties of biological systems play an important role in computational systems biology.
 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.  two parameter estimation methods of combining spline theory with linear programming and nonlinear programming 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.  the proposed approaches have general application to identify unknown parameter values of a wide range of systems biology models. 