BACKGROUND
elucidation of biological systems has gained wider interest in the last decade.
despite recent advances, fundamental understanding of life processes still requires powerful theoretical tools from mathematics and physical sciences.
particularly, mathematical modelling of biological systems is nowadays becoming an essential partner of experimental work.
one of the most challenging tasks in computational modelling of biological systems is the estimation of the model parameters.
the aim here is to obtain the set of parameter values that make the model response consistent with the data observed.
parameter estimation can be formulated as an optimization problem in which the sum of squared residuals between the measured and simulated data is minimized.
the biological model dictates the type of optimization problem being faced.
many biological systems are described through nonlinear ordinary differential equations  that provide the concentration profiles of certain metabolites over time.
recent methodological developments have enabled the generation of some dynamic profiles of gene networks and protein expression data, although the latter are still very rare.
in this context, there is a strong motivation for developing systematic techniques for building dynamic biological models from experimental data.
the parameter estimation of these models gives rise to dynamic optimization problems which are hard to solve.
existing approaches to optimize dynamic models can be roughly classified as direct or indirect   <cit> .
direct methods make use of gradient-based nonlinear programming  solvers and can in turn be divided into sequential and simultaneous.
in sequential approaches, the optimization of the control variables, which are discretized, is performed by a nlp solver, whereas the ode is calculated externally, that is, both steps are executed in a sequential manner.
in contrast, in simultaneous strategies, both the control and state profiles are approximated using polynomials  and discretized in time by means of finite elements  <cit> .
in the latter strategy, the ode system is replaced by a system of algebraic equations that is optimized with a standard gradient-based nlp solver.
simultaneous approaches can handle dynamic systems with unstable modes and with path constraints  <cit> .
furthermore, they allow performing automatic differentiation with respect to the control and state variables, avoiding the need to calculate the derivatives numerically as is the case in the sequential approach.
unfortunately, the discretization step can lead to large scale nlps that are difficult to solve.
models of biological systems are typically highly nonlinear, which gives rise to nonconvex optimization problems with multiple local solutions .
because of this, traditional gradient-based methods used in the sequential and simultaneous approaches may fall in local optima.
in the context of parameter estimation, these local solutions should be avoided, since they may lead to inaccurate models that are unable to predict the system’s performance precisely.
global optimization  algorithms are a special class of techniques that attempt to identify the global optimum in nonconvex problems.
these methods can be classified as stochastic and deterministic.
stochastic go methods are based on probabilistic algorithms that provide near optimal solutions in short cpu times.
despite having shown great potential with large-scale problems like parameter estimation  <cit> , these methods have as major limitation that are unable to guarantee convergence to the global optimum in a finite number of iterations.
in other words, they provide solutions whose optimality  is unknown, and may or may not be globally optimal.
in contrast, deterministic global optimization methods ensure global optimality within a desired tolerance, but lead to larger computational burdens.
hence, in addition to the solution itself, these methods provide as output a rigorous interval within which the best possible solution  must fall.
despite recent advances in deterministic global optimization methods  <cit> , their application to parameter estimation has been quite scarce.
two main deterministic go methods exist: spatial branch-and-bound   <cit> , and outer approximation  <cit> .
both algorithms rely on computing valid lower and upper bounds on the global optimum.
these bounds tend to approach as iterations proceed, thus offering a theoretical guarantee of convergence to the global optimum.
a rigorous lower bound on the global optimum of the original nonconvex problem is obtained by solving a valid relaxation that contains its feasible space.
to construct this relaxed problem, the nonconvex terms in the original formulation are replaced by convex envelopes that overestimate its feasible region.
there are different types of convex envelopes that provide relaxations for a wide variety of nonconvexities.
these relaxations are the main ingredient of deterministic go methods and play a key role in their performance.
in general, tighter relaxations provide better bounds , thereby expediting the overall solution procedure.
to the best of our knowledge, esposito and floudas were the first to propose a deterministic method for the global solution of dynamic optimization problems with embedded odes  <cit> .
their approach relies on reformulating the problem as a nonconvex nlp using orthogonal collocation on finite elements.
this reformulated nlp was then solved by means of a sbb method.
to this end, they constructed a convex relaxation of the reformulated problem following the αbb approach previously proposed by the authors  <cit> .
despite being valid for twice continuous differentiable functions, these relaxations may provide weak bounds in some particular cases and therefore lead to large cpu times when used in the context of a spatial branch and bound framework  <cit> .
this work proposes a computational framework for the deterministic global optimization of parameter estimation problems of nonlinear dynamic biological systems.
the main contributions of our work are:  the application of deterministic global optimization methods to dynamic models of biological systems, and  the use of several known techniques employed in dynamic  and global optimization  in the context of an outer approximation algorithm.
the approach presented relies on discretizing the set of nonlinear odes using orthogonal collocation on finite elements, thereby transforming the dynamic system into an equivalent nonconvex nlp problem.
a customized outer approximation algorithm that relies on a mixed-integer linear programming  relaxation is used in an iterative scheme along with the aforementioned nlp to solve the nonconvex model to global optimality.
the milp relaxation is tightened using a special type of cutting plane that exploits the problem structure, thereby expediting the overall solution procedure.
the capabilities of our algorithm are tested through its application to two case studies: the isomerisation of α-pinene  and the inhibition of hiv proteinase .
the results obtained are compared with those produced by the state-of-art commercial global optimization package baron .
our algorithm is proved from these numerical examples to produce near optimal solutions in a fraction of the cpu time required by baron.
methods
problem statement
the problem addressed in this work can be stated as follows: given is a dynamic kinetic model describing the mechanism of a set of biochemical reactions.
the goal is to determine the appropriate values of the model coefficients , so as to minimize the sum-of-squares of the residuals between the simulated data provided by the model and the experimental observations.
mathematical formulation
we consider dynamic parameter estimation optimization problems of the following form:
minθ,z^u∑j∈jm∑u∈u <dig>
s.t.z˙j=g∀j∈j
zj=z0∀j∈j
t∈
z^u,j=zj∀u∈u;∀j∈jm
where z˙ represents the state variables , z0
their initial conditions, z^u,j represents the experimental data variables, z¯u,j are the experimental observations, j is the set of state variables whose derivatives explicitly appear in the model, θ are the parameters to be estimated and tu, is the time associated with the uth experimental data point in the set u.
our solution strategy relies on reformulating the nonlinear dynamic optimization problem as a finite-dimensional nlp by applying a complete discretization using orthogonal collocation on finite elements.
this nlp is next solved using an outer approximation algorithm .
in the sections that follow, we explain in detail the main steps of our algorithm.
orthogonal collocation approach
there is a considerable number of collocation-based discretizations for the solution of differential-algebraic systems  <cit> .
without loss of generality, we employ herein the so-called orthogonal collocation on finite elements method  <cit> .
consider the following set of ode’s defined as
z˙j=g∀j∈j
the state variables are first approximated using lagrange polynomials as follows:
znk+1=∑k=0nkξkϕkϕk=∏q= <dig> q≠knkt−tqtk−tq
these polynomials have the property that at the orthogonal collocation points their coefficients, ξk, take the value of the state profile at that point.
therefore, the collocation coefficients ξk acquire physical meaning which allows to generate bounds for these variables.
because state variables may present steep variations, the whole solution space is commonly divided into time intervals called finite elements.
hence, the time variable t is divided into ne elements of length Δ
ηe and rescaled as τ∈ <cit> .
within each finite element, nk +  <dig> orthogonal collocation points τ, τ, τ, ⋯ ,τ
are distributed at the shifted  roots of the orthogonal legendre polynomial of nk degree.
recall that the 0th orthogonal collocation point is located at the beginning of each element .
collocation points evaluated at the shifted orthogonal legendre polynomials.
following the collocation method  <cit> , the residual equations arising from the combination of eqs.
<dig> and  <dig>  are defined for each element e in the set e and state variable in the set j, giving rise to the following constraints:
∑k=0nkξe,k,jϕ˙e,k,j−Δηegj=0∀e∈ek′= <dig> …,nk;∀j∈j
the state variables have to be continuous between elements, so we enforce the following continuity constrains:
ξe, <dig> j−∑k=0nkξe− <dig> k,jϕk=0e= <dig> …,ne∀j∈j
these equations extrapolate the polynomial at element e- <dig>  providing an accurate initial condition for the next element e.
moreover, initial conditions are enforced for the beginning of the first element using the following equation:
ξ <dig> ,j−z <dig> j=0∀j∈j
recall that collocation points in which time has been discretized will not necessarily match the times at which experimental profiles were registered.
hence, variable z^u,j is added to determine the value of the model states profiles at times tu
making it possible to fit the model to the experimental points.
this is accomplished by adding the following equation:
−z^u,j+∑k=0nkξeu,k,jϕk=0∀u∈u;∀j∈jm
where τu
is calculated as follows:
τu=tu−ηeuΔηeu
here, the subscript eu refers to the element where tu
falls, that is, eu≡{e:ηe≤tu<ηe + 1}.
npl formulation
the dynamic optimization problem is finally reformulated into the following nlp:
minθ,ξ,z^u∑j∈jm∑u∈u <dig>
s.t.∑k=0nkξe,k,jϕ˙e,k,j−Δηegj=0∀e∈ek′= <dig> …,nk;∀j∈j
ξe, <dig> j−∑k=0nkξe− <dig> k,jϕk=0e= <dig> …,ne∀j∈j
ξ <dig> ,j−z <dig> j=0∀j∈j
−z^u,j+∑k=0nkξeu,k,jϕk=0∀u∈u;∀j∈jm
RESULTS
optimization approach
the method devised for globally optimizing the nlp that arises from the reformulation of the parameter estimation problem  is based on an outer approximation algorithm  <cit>  used by the authors in previous works  <cit> .
this approach relies on decomposing the original nlp into two subproblems at different hierarchical levels: a lower level milp problem and an upper level slave nlp problem.
the master problem is a relaxation of the original nlp  and hence provides a rigorous lower bound on its global optimum.
the slave nlp yields a valid upper bound when it is solved locally.
the algorithm iterates between these two levels until the optimality gap  is reduced below a given tolerance .
in the following subsections, we provide a detailed description of the algorithm.
lower level master problem
designing efficient and smart strategies for attaining tight bounds is a mayor challenge in deterministic global optimization.
both the quality of the bounds and the time required to generate them drastically influence the overall performance of a deterministic global optimization algorithm.
any feasible solution of the original nlp is a valid upper bound and can be obtained by means of a local nlp solver.
to obtain lower bounds, we require a rigorous convex  relaxation.
this relaxation is obtained by replacing the nonconvex terms by convex overestimators.
since the relaxed problem is convex, it is possible to solve it to global optimality using standard local optimizers.
furthermore, since its feasible region contains that of the original problem and its objective function rigorously underestimates the original one, it is guaranteed to provide a lower bound on the global optimum of the original nonconvex model  <cit> .
androulakis et al.
<cit>  proposed a convex quadratic relaxation for nonconvex functions named αbb underestimator which can be applied to general twice continuously differentiable functions.
this technique, which was used in parameter estimation by esposito and floudas  <cit> , might lead in some cases to weak relaxations and therefore poor numerical performance  <cit> .
to construct a valid milp relaxation, we apply the following approach.
we first reformulate the nlp using the symbolic reformulation method proposed by smith and pantelides  <cit> .
this technique reformulates any system of nonlinear equations into an equivalent canonical form with the only nonlinearities being bilinear products, linear fractional, simple exponentiation and univariate function terms with the following standard form:
minwwobj
s.t.aw=b
wl≤w≤wu
y∈
wk≡wiwj∀∈tbt
wk≡wiwj∀∈tlft
wk≡win∀∈tet
wk≡fn∀∈tuft
where vector w comprises continuous variables x as well as integers y, while the sets tbt, tlft, tet and tuft are the bilinear product, linear fractional, simple exponentiation and univariate function terms, respectively.
a rigorous relaxation of the original model is constructed by replacing the nonconvex terms in the reformulated model by convex estimators.
the solution of the convex relaxation provides a valid lower bound on the global optimum.
more precisely, the bilinear terms are replaced by piecewise mccormick relaxations.
the fractional terms can be convexified in two different manners.
the first is to replace them by tailored convex envelopes that exploit their structure  <cit> .
the second is to express them as bilinear terms by performing a simple algebraic transformation, and then use the mccormick envelopes to relax the associated bilinear term.
univariate functions commonly used in process engineering models  are purely convex or purely concave, and can be replaced by the exact function-secant pair estimators  <cit> .
the reader is referred to the work by smith and pantelides  <cit>  for further details on the symbolic reformulation.
we focus next on explaining how the bilinear terms are approximated in the reformulated nlp.
piecewise mccormick-based relaxation
the bilinear terms appearing in the reformulated model are approximated using mccormick’s envelopes  <cit> .
for bilinear terms, this relaxation is tighter than the αbb-based relaxations  <cit> .
each bilinear term xy can be replaced by an auxiliary variable z as follows:
z=xyxl≤x≤xuyl≤y≤yu
the best known relaxation for approximating a bilinear term is given by the mccormick envelopes, obtained by replacing eq.
<dig> by the following linear under , and overestimators :
z≥xyl+xly−xlyl
z≥xyu+xuy−xuyu
z≤xyl+xuy−xuyl
z≤xyu+xly−xlyu
in this work we further tighten the mccormick envelopes by adding binary variables  <cit> .
particularly, two additional sets of variables are defined in the piecewise formulation:
binary switch: λ∈{ <dig> }np
continuous switch: Δy∈np
the binary switch λ
is active =1) for the segment where x is located ≤x≤xl+anp) and is otherwise inactive.
therefore, the partitioning scheme activates exactly only one np∈{ <dig> …,np} so that the feasible region corresponding to the relaxation of xy is reduced from the parallelogram in figure  <dig> to a significantly smaller one depicted in figure  <dig>
eq.
<dig> enforces that only one binary variable is active:
∑np=1npλ= <dig>
the continuous switch Δy takes on any positive value between  <dig> and yu−yl when the binary switch corresponding to the npth piecewise λ is active =1) and  <dig> otherwise.
therefore:
y=yl+∑np=1npΔy
0≤Δy≤λnp= <dig> …,np
finally, the under and overestimators for the active segment are defined in algebraic terms as follows:
z≥xyl+∑np=1npΔy
z≥xyu+∑np=1np
z≤xyl+∑np=1npΔy
z≤xyu+∑np=1np
xl≤x≤xu;yl≤y≤yu
note that the discrete relaxation is tighter than the continuous one over the entire feasible region.
the introduction of the binary variables required in the piecewise mccormick reformulation gives rise to a mixed-integer nonlinear programming  problem, with the only nonlinearities appearing in the objective function.
while this minlp is convex and can be easily solved to global optimality with standard minlp solvers, it is more convenient to linearize it in order to obtain an milp formulation, for which more efficient software packages exist.
the section that follows explains how this is accomplished.
hyper-planes underestimation
the convex minlp can be further reformulated into an milp by replacing the objective function by a set of hyper-planes.
for this, we define two new variables as zu,j′=z^u,j−z¯u,j and α≥z′u,j <dig>  the quadratic terms are then approximated by 1st degree taylor series.
that is, the square terms are replaced by l hyper-planes uniformly distributed between the maximum and minimum desired values of zu,j′  so that the objective function is reduced to a summation of quadratic terms as follows:
minθ,ξ,z^u∑j∈jm∑u∈uαu,j
αu,j≥z′02u,j,l+2z0′u,j,l∀u∈u∀j∈jm∀l∈l
upper level slave problem
a valid upper bound on the global optimum is obtained by optimizing the original nlp locally.
this nlp is initialized using the solution provided by the milp as starting point.
the solution of this nlp is used to tighten the milp, so the lower and upper bounds tend to converge as iterations proceed.
algorithm steps
the proposed algorithm comprises the following steps:
<dig>  set iteration count it =  <dig>  ub = ∞, lb = −∞
and tolerance error = tol.
<dig>  set it=it +  <dig>  solve the master problem milp.
if the milp is infeasible, stop .
otherwise, update the current lb making lb=maxit, where lbit
is the value of the objective function of the milp in the itth
iteration.
<dig>  solve the slave problem nlp.
if the nlp is infeasible add one more piecewise term and hyper-plane to the master milp and go to step  <dig> of the algorithm.
otherwise, update the current ub making ub = minit , where ubit
is the value of the objective function of the nlp in the itth
iteration.
<dig>  calculate the optimality gap og as og=|ub−lb|ub.
if og≤tol, then stop.
the current ub is regarded as the global optimum within the desired tolerance.
otherwise, add one more piecewise section and hyper-plane to the master milp and go to step  <dig> of the algorithm.
remarks:
there are different methods to update the piecewise bilinear approximation.
one possible strategy is to update it by dividing the active piecewise  into two equal-length segments.
the new hyper-plane term z0′u,j,l is added at the optimal solution of the milp  in the previous iteration.
the univariate convex and concave terms in the reformulated problem can be either approximated by the secant or by a piecewise univariate function similarly as done with the mccormick envelopes.
our algorithm needs to be tuned prior to its application.
this is a common practice in any optimization algorithm.
in a previous publication  <cit> , we studied the issue of defining the number of piecewise intervals and supporting hyper-planes in an optimal manner.
in practice, however, the optimal number of piecewise terms and hyper-planes is highly dependent on the specific instance being solved, so it is difficult to provide general guidelines on this.
the approach presented might lead to large computational burdens in large-scale models of complex biological systems.
future work will focus on expediting our algorithm through the addition of cutting planes and the use of customized decomposition strategies.
case studies
we illustrate the performance of the proposed algorithm through its application to two challenging benchmark parameter estimation problems: the isomerisation of α-pinene  and the inhibition of hiv proteinase .
the objective in these problems is to obtain the set of values of the model parameters such that the model response is as close as possible to the experimental data.
for comparison purposes we used the global optimization package baron .
baron is a commercial software for solving nonconvex optimization problems to global optimality.
baron combines constraint propagation, interval analysis, duality, and enhanced "branch and bound" concepts for efficient range reduction with rigorous relaxations constructed by enlarging the feasible region and/or underestimating the objective function.
the interested readers have the possibility to evaluate this software on their own for free in this link: http://www.neos-server.org/neos/solvers/go:baron/gams.html.
our algorithm was implemented in gams  <dig> .
<dig> using cplex  <dig> .
<dig>  for the milps and snopt  <dig> for the nlps subproblems.
all the calculations were performed in a pc/amd athlon ii at  <dig>  ghz using a single core.
data about the size of the models can be found in table  <dig>
case study 1: isomerisation of α-pinene
in this first case study, five kinetic parameters describing the thermal isomerisation of α-pinene are estimated.
the proposed reaction scheme for this process is depicted in figure  <dig>  in this homogeneous chemical reaction, α-pinene  is thermally isomerised to dipentene  and allo-ocimene , which in turn yields α- and β-pyronene  and a dimer .
this process was originally studied by fuguitt and hawkins  <cit> , which carried out a single experiment reporting the experimental concentrations  of the reactant and the four products measured at eight time intervals.
hunter and mcgregor  <cit>  postulated first-order kinetics and proposed the following set of ode’s describing the dynamic process:
dγ1dt=−γ <dig>
dγ2dt=p1γ <dig>
dγ3dt=p2γ1−γ3+p5γ <dig>
dγ4dt=p3γ <dig>
dγ5dt=p4γ3−p5γ <dig>
γ0= <cit> t∈ <cit>
rodriguez-fernandez et al.
<cit>  addressed this problem by applying a metaheuristic based on the scatter search method.
this strategy does not offer any theoretical guarantee of convergence to the global optimum in a finite number of iterations.
following our approach, the state variables were approximated by lagrange polynomials using three collocation points evaluated at the shifted roots of orthogonal legendre polynomials and defining five finite elements of equal length.
the nonconvexities in the resulting residual equations are given by the bilinear terms θiξe,k, j which were relaxed using piecewise mccormick approximations as described previously.
the objective function was underestimated using supporting hyper-planes.
it is well known that the quality of the lower bound predicted by a relaxation strongly depends on the bounds imposed on its variables  <cit> .
hence bounds on collocation coefficients  were tightened by performing a bound contraction procedure  <cit> .
particularly, tight lower and upper bounds were estimated for each collocation coefficient by maximizing and minimizing its value while satisfying the constraints contained in the master problem.
this is a costly process .
for this reason, it was only performed recursively  <dig> times before the initialization of the algorithm.
the milp was further tightened by adding the following constraint:
∑j∈jm∑u∈u2≤ <dig>
which forces the model to find a solution better than the one obtained at the beginning of the search by locally minimizing the original nlp .
furthermore, the parameter θi was allowed to take any value within the  <cit>
interval.
the problem was solved with  <dig> initial hyper-planes.
an extra hyper-plane was added in each iteration, but the total number of piecewise terms was kept constant  in order to keep the milp in a manageable size.
a tolerance of 5% was set as termination criterion.
for comparison purposes, we solved the same problem with the standard global optimization package baron using its default settings.
baron was able to find the global optimum but failed at reducing the optimality gap below the specified tolerance after 12h of cpu time.
in contrast, our algorithm closed the gap in less than 3h .
as shown in table  <dig>  the results obtained agree with those reported in the literature.
case study 2: inhibition of hiv proteinase
in this second case study, we considered a much more complex biological dynamic system.
particularly, we studied the reaction mechanism of the irreversible inhibition of hiv proteinase, as originally examined by kuzmic  <cit>  .
note that this dynamic model has lack of practical identifiability, as reported in rodriguez-fernandez et al <cit> .
nevertheless, we think that this example is still useful for the purpose of our analysis, since the emphasis here is placed on globally optimizing dynamic models of biological systems rather than analyzing identifiability issues.
the model can be described mathematically through a set of  <dig> nonlinear ode’s with ten parameters:
ddt=−2k11+2k <dig>
ddt=k3−2k41+2k <dig>
ddt=−k21+k <dig>
ddt=−k51+k <dig>
ddt=k21−k22−k <dig>
ddt=k41−k <dig>
ddt=k11−k12−k21+k22+k3−k41+k42−k51+k <dig>
ddt=k51−k52−k <dig>
ddt=k <dig>
where the following initial conditions and parameters are known:
0=00=00=00=00=00=00=00= <dig> = <dig> = <dig> = <dig>
k11= <dig> k12= <dig> k41=100k21=100k51= <dig>
t∈ <cit>
a series of five experiments where the enzyme hiv proteinase   was added to a solution of an irreversible inhibitor  and a fluorogenic substrate   were considered.
the five experiments were carried out at four different concentrations of the inhibitor .
the fluorescence changes were monitored during one hour.
the measured signal is a linear function of the product  concentration, as expressed in the following equation:
signal=ϵ+offset
in this fit, the offset  of the fluorimeter was considered as a degree of freedom.
a certain degree of uncertainty  was assumed for the value of the initial concentrations of substrate and enzyme .
the calibration of a total of  <dig> adjustable parameters was addressed: five rate constants, five initial concentrations of enzyme and substrate and five offset values.
mendes and kell  <cit>  solved this problem using simulated annealing and reported its first known solution.
later, rodriguez-fernandez et al.
<cit>  improved that solution by means of a scatter search metaheuristic, which required a fraction of the time employed by mendes’ simulated annealing.
recall that, despite producing near optimal solutions in short cpu times, stochastic algorithms provide no information on the quality of the solutions found and are unable to guarantee convergence to the global optimum in a finite number of iterations.
on the contrary, the proposed methodology ensures the global optimality of the solution computed within a desired tolerance.
in our study, the state variables were approximated using five orthogonal collocation points and five equal-length finite elements.
in this case, the nonconvexities arise from the bilinear terms ξe, k, jξe, k, j and θiξe, k, j.
the parameter bounds θi were set to θi∈ <cit> .
the lower and upper limits for the collocation coefficients ξe, k, j, n
were fixed to ξe, k, j, n∈
except for ξe, k, e, n∈ and ξe, k, s, n∈.
the bounds for all the offsets were set to offsetn∈.
the master problem was further tightened by adding a special type of strengthening cuts.
these cuts are generated by temporally decomposing the original full space milp into a series of milps in each of which we fit only a subset of the original dataset, and remove the continuity equations corresponding to the extreme elements included in the sub-problem.
the cuts are expressed as inequalities added to the master problem that impose lower bounds on the error of a subset of elements for which the sub-milps are solved.
these bounds are hence obtained from the solution of a set of milp sub-problems that optimize the error of only a subset of elements.
this case study was solved with  <dig> initial piecewise intervals and  <dig> initial hyper-planes.
two strengthening cuts involving elements  <dig>   <dig>   <dig> and  <dig>  and  <dig>   <dig>   <dig> and  <dig>  respectively, were added as constrains.
a tolerance of 20% was used in the calculations.
hyper-planes and piecewise terms were updated at each iteration of the algorithm.
in this case, baron failed to identify any feasible solution after 12h of cpu time.
in contrast, our algorithm was able to obtain the global optimum  with a gap of  <dig> % in approximately  <dig>  cpu s .
remarkably, the solution found by our algorithm improves the best known solution reported by rodriguez-fernandez et al.
<cit> .
hence, our algorithm clearly outperformed other parameter estimation methods, improving the best known solution  <cit> , and providing a rigorous lower bound on the minimum error that can be attained.
CONCLUSIONS
in this work, we have proposed a novel strategy for globally optimizing parameter estimation problems with embedded nonlinear dynamic systems.
the method presented was tested through two challenging benchmark problems: the isomerisation of α-pinene  and the inhibition of hiv proteinase .
the proposed algorithm identified the best known solution, which was originally reported by rodriguez-fernandez et al.
<cit> , in the case of the α-pinene, and improved the best known one in the hiv proteinase case study.
in both cases, rigorous lower bounds were provided on the global optimum, making it possible to determine the optimality gap of the solutions found.
the method proposed produced promising results, surpassing the capabilities of baron.
our method requires some knowledge on optimization theory as well as skills using modelling systems.
our final goal is to develop a software to automate the calculations, so our approach can be easily used by a wider community.
this is a challenging task, since nonlinear models are hard to handle and typically require customized solution procedures.
particularly, nonlinear models must be initialized carefully to ensure convergence even to a local solution.
in this regard, the use of an outer approximation scheme that relies on a master milp formulation is quite appealing, since the outcome of this milp can be used to initialize the nlp in a robust manner.
another key point here is how to construct tight relaxations of the nonconvex terms.
an efficient algorithm must exploit the problem structure to obtain high quality relaxations and therefore good bounds close to the global optimum.
these relaxations can be further tightened through the addition of cutting planes or the use of customized decomposition methods.
as observed, there is still much work to be done in this area, but we strongly believe that such an effort is worthy.
furthermore, recent advances in global optimization theory and software applications are paving the way to develop systematic deterministic tools for the global optimization of parameter estimation problems of increasing size.
our future work will focus on making the approach more efficient through the use of tailored cutting planes and decomposition strategies and also through the hybridization of deterministic methods with stochastic approaches.
competing interest
the authors declare that they have no competing interests.
authors’ contributions
g. g-g suggested the need for the approach and j.a.
e provided the biological problems.
a. m, c. p, g. g-g and l. j developed the optimization algorithms and performed the numerical analysis.
all authors evaluated the results, wrote the paper and contributed to its final form.

 elucidation of biological systems has gained wider interest in the last decade. despite recent advances, fundamental understanding of life processes still requires powerful theoretical tools from mathematics and physical sciences. particularly, mathematical modelling of biological systems is nowadays becoming an essential partner of experimental work. one of the most challenging tasks in computational modelling of biological systems is the estimation of the model parameters. the aim here is to obtain the set of parameter values that make the model response consistent with the data observed. parameter estimation can be formulated as an optimization problem in which the sum of squared residuals between the measured and simulated data is minimized. the biological model dictates the type of optimization problem being faced. many biological systems are described through nonlinear ordinary differential equations that provide the concentration profiles of certain metabolites over time. recent methodological developments have enabled the generation of some dynamic profiles of gene networks and protein expression data, although the latter are still very rare. in this context, there is a strong motivation for developing systematic techniques for building dynamic biological models from experimental data. the parameter estimation of these models gives rise to dynamic optimization problems which are hard to solve. existing approaches to optimize dynamic models can be roughly classified as direct or indirect. direct methods make use of gradient-based nonlinear programming solvers and can in turn be divided into sequential and simultaneous. in sequential approaches, the optimization of the control variables, which are discretized, is performed by a nlp solver, whereas the ode is calculated externally, that is, both steps are executed in a sequential manner. in contrast, in simultaneous strategies, both the control and state profiles are approximated using polynomials and discretized in time by means of finite elements. in the latter strategy, the ode system is replaced by a system of algebraic equations that is optimized with a standard gradient-based nlp solver. simultaneous approaches can handle dynamic systems with unstable modes and with path constraints. furthermore, they allow performing automatic differentiation with respect to the control and state variables, avoiding the need to calculate the derivatives numerically as is the case in the sequential approach. unfortunately, the discretization step can lead to large scale nlps that are difficult to solve. models of biological systems are typically highly nonlinear, which gives rise to nonconvex optimization problems with multiple local solutions. because of this, traditional gradient-based methods used in the sequential and simultaneous approaches may fall in local optima. in the context of parameter estimation, these local solutions should be avoided, since they may lead to inaccurate models that are unable to predict the system’s performance precisely. global optimization algorithms are a special class of techniques that attempt to identify the global optimum in nonconvex problems. these methods can be classified as stochastic and deterministic. stochastic go methods are based on probabilistic algorithms that provide near optimal solutions in short cpu times. despite having shown great potential with large-scale problems like parameter estimation, these methods have as major limitation that are unable to guarantee convergence to the global optimum in a finite number of iterations. in other words, they provide solutions whose optimality is unknown, and may or may not be globally optimal. in contrast, deterministic global optimization methods ensure global optimality within a desired tolerance, but lead to larger computational burdens. hence, in addition to the solution itself, these methods provide as output a rigorous interval within which the best possible solution must fall. despite recent advances in deterministic global optimization methods, their application to parameter estimation has been quite scarce. two main deterministic go methods exist: spatial branch-and-bound, and outer approximation. both algorithms rely on computing valid lower and upper bounds on the global optimum. these bounds tend to approach as iterations proceed, thus offering a theoretical guarantee of convergence to the global optimum. a rigorous lower bound on the global optimum of the original nonconvex problem is obtained by solving a valid relaxation that contains its feasible space. to construct this relaxed problem, the nonconvex terms in the original formulation are replaced by convex envelopes that overestimate its feasible region. there are different types of convex envelopes that provide relaxations for a wide variety of nonconvexities. these relaxations are the main ingredient of deterministic go methods and play a key role in their performance. in general, tighter relaxations provide better bounds, thereby expediting the overall solution procedure. to the best of our knowledge, esposito and floudas were the first to propose a deterministic method for the global solution of dynamic optimization problems with embedded odes. their approach relies on reformulating the problem as a nonconvex nlp using orthogonal collocation on finite elements. this reformulated nlp was then solved by means of a sbb method. to this end, they constructed a convex relaxation of the reformulated problem following the αbb approach previously proposed by the authors. despite being valid for twice continuous differentiable functions, these relaxations may provide weak bounds in some particular cases and therefore lead to large cpu times when used in the context of a spatial branch and bound framework. this work proposes a computational framework for the deterministic global optimization of parameter estimation problems of nonlinear dynamic biological systems. the main contributions of our work are: the application of deterministic global optimization methods to dynamic models of biological systems, and the use of several known techniques employed in dynamic and global optimization in the context of an outer approximation algorithm. the approach presented relies on discretizing the set of nonlinear odes using orthogonal collocation on finite elements, thereby transforming the dynamic system into an equivalent nonconvex nlp problem. a customized outer approximation algorithm that relies on a mixed-integer linear programming relaxation is used in an iterative scheme along with the aforementioned nlp to solve the nonconvex model to global optimality. the milp relaxation is tightened using a special type of cutting plane that exploits the problem structure, thereby expediting the overall solution procedure. the capabilities of our algorithm are tested through its application to two case studies: the isomerisation of α-pinene and the inhibition of hiv proteinase. the results obtained are compared with those produced by the state-of-art commercial global optimization package baron. our algorithm is proved from these numerical examples to produce near optimal solutions in a fraction of the cpu time required by baron. methods problem statement the problem addressed in this work can be stated as follows: given is a dynamic kinetic model describing the mechanism of a set of biochemical reactions. the goal is to determine the appropriate values of the model coefficients, so as to minimize the sum-of-squares of the residuals between the simulated data provided by the model and the experimental observations. mathematical formulation we consider dynamic parameter estimation optimization problems of the following form: minθ,z^u∑j∈jm∑u∈u s.t.z˙j=g∀j∈j zj=z0∀j∈j t∈ z^u,j=zj∀u∈u;∀j∈jm where z˙ represents the state variables, z0 their initial conditions, z^u,j represents the experimental data variables, z¯u,j are the experimental observations, j is the set of state variables whose derivatives explicitly appear in the model, θ are the parameters to be estimated and tu, is the time associated with the uth experimental data point in the set u. our solution strategy relies on reformulating the nonlinear dynamic optimization problem as a finite-dimensional nlp by applying a complete discretization using orthogonal collocation on finite elements. this nlp is next solved using an outer approximation algorithm. in the sections that follow, we explain in detail the main steps of our algorithm. orthogonal collocation approach there is a considerable number of collocation-based discretizations for the solution of differential-algebraic systems. without loss of generality, we employ herein the so-called orthogonal collocation on finite elements method. consider the following set of ode’s defined as z˙j=g∀j∈j the state variables are first approximated using lagrange polynomials as follows: znk+1=∑k=0nkξkϕkϕk=∏q= q≠knkt−tqtk−tq these polynomials have the property that at the orthogonal collocation points their coefficients, ξk, take the value of the state profile at that point. therefore, the collocation coefficients ξk acquire physical meaning which allows to generate bounds for these variables. because state variables may present steep variations, the whole solution space is commonly divided into time intervals called finite elements. hence, the time variable t is divided into ne elements of length Δ ηe and rescaled as τ∈. within each finite element, nk + orthogonal collocation points τ, τ, τ, ⋯,τ are distributed at the shifted roots of the orthogonal legendre polynomial of nk degree. recall that the 0th orthogonal collocation point is located at the beginning of each element. collocation points evaluated at the shifted orthogonal legendre polynomials. following the collocation method, the residual equations arising from the combination of eqs. and are defined for each element e in the set e and state variable in the set j, giving rise to the following constraints: ∑k=0nkξe,k,jϕ˙e,k,j−Δηegj=0∀e∈ek′= …,nk;∀j∈j the state variables have to be continuous between elements, so we enforce the following continuity constrains: ξe, j−∑k=0nkξe− k,jϕk=0e= …,ne∀j∈j these equations extrapolate the polynomial at element e- providing an accurate initial condition for the next element e. moreover, initial conditions are enforced for the beginning of the first element using the following equation: ξ,j−z j=0∀j∈j recall that collocation points in which time has been discretized will not necessarily match the times at which experimental profiles were registered. hence, variable z^u,j is added to determine the value of the model states profiles at times tu making it possible to fit the model to the experimental points. this is accomplished by adding the following equation: −z^u,j+∑k=0nkξeu,k,jϕk=0∀u∈u;∀j∈jm where τu is calculated as follows: τu=tu−ηeuΔηeu here, the subscript eu refers to the element where tu falls, that is, eu≡{e:ηe≤tu<ηe + 1}. npl formulation the dynamic optimization problem is finally reformulated into the following nlp: minθ,ξ,z^u∑j∈jm∑u∈u s.t.∑k=0nkξe,k,jϕ˙e,k,j−Δηegj=0∀e∈ek′= …,nk;∀j∈j ξe, j−∑k=0nkξe− k,jϕk=0e= …,ne∀j∈j ξ,j−z j=0∀j∈j −z^u,j+∑k=0nkξeu,k,jϕk=0∀u∈u;∀j∈jm  optimization approach the method devised for globally optimizing the nlp that arises from the reformulation of the parameter estimation problem is based on an outer approximation algorithm used by the authors in previous works. this approach relies on decomposing the original nlp into two subproblems at different hierarchical levels: a lower level milp problem and an upper level slave nlp problem. the master problem is a relaxation of the original nlp and hence provides a rigorous lower bound on its global optimum. the slave nlp yields a valid upper bound when it is solved locally. the algorithm iterates between these two levels until the optimality gap is reduced below a given tolerance. in the following subsections, we provide a detailed description of the algorithm. lower level master problem designing efficient and smart strategies for attaining tight bounds is a mayor challenge in deterministic global optimization. both the quality of the bounds and the time required to generate them drastically influence the overall performance of a deterministic global optimization algorithm. any feasible solution of the original nlp is a valid upper bound and can be obtained by means of a local nlp solver. to obtain lower bounds, we require a rigorous convex relaxation. this relaxation is obtained by replacing the nonconvex terms by convex overestimators. since the relaxed problem is convex, it is possible to solve it to global optimality using standard local optimizers. furthermore, since its feasible region contains that of the original problem and its objective function rigorously underestimates the original one, it is guaranteed to provide a lower bound on the global optimum of the original nonconvex model. androulakis et al. proposed a convex quadratic relaxation for nonconvex functions named αbb underestimator which can be applied to general twice continuously differentiable functions. this technique, which was used in parameter estimation by esposito and floudas, might lead in some cases to weak relaxations and therefore poor numerical performance. to construct a valid milp relaxation, we apply the following approach. we first reformulate the nlp using the symbolic reformulation method proposed by smith and pantelides. this technique reformulates any system of nonlinear equations into an equivalent canonical form with the only nonlinearities being bilinear products, linear fractional, simple exponentiation and univariate function terms with the following standard form: minwwobj s.t.aw=b wl≤w≤wu y∈ wk≡wiwj∀∈tbt wk≡wiwj∀∈tlft wk≡win∀∈tet wk≡fn∀∈tuft where vector w comprises continuous variables x as well as integers y, while the sets tbt, tlft, tet and tuft are the bilinear product, linear fractional, simple exponentiation and univariate function terms, respectively. a rigorous relaxation of the original model is constructed by replacing the nonconvex terms in the reformulated model by convex estimators. the solution of the convex relaxation provides a valid lower bound on the global optimum. more precisely, the bilinear terms are replaced by piecewise mccormick relaxations. the fractional terms can be convexified in two different manners. the first is to replace them by tailored convex envelopes that exploit their structure. the second is to express them as bilinear terms by performing a simple algebraic transformation, and then use the mccormick envelopes to relax the associated bilinear term. univariate functions commonly used in process engineering models are purely convex or purely concave, and can be replaced by the exact function-secant pair estimators. the reader is referred to the work by smith and pantelides for further details on the symbolic reformulation. we focus next on explaining how the bilinear terms are approximated in the reformulated nlp. piecewise mccormick-based relaxation the bilinear terms appearing in the reformulated model are approximated using mccormick’s envelopes. for bilinear terms, this relaxation is tighter than the αbb-based relaxations. each bilinear term xy can be replaced by an auxiliary variable z as follows: z=xyxl≤x≤xuyl≤y≤yu the best known relaxation for approximating a bilinear term is given by the mccormick envelopes, obtained by replacing eq. by the following linear under, and overestimators: z≥xyl+xly−xlyl z≥xyu+xuy−xuyu z≤xyl+xuy−xuyl z≤xyu+xly−xlyu in this work we further tighten the mccormick envelopes by adding binary variables. particularly, two additional sets of variables are defined in the piecewise formulation: binary switch: λ∈{ }np continuous switch: Δy∈np the binary switch λ is active =1) for the segment where x is located ≤x≤xl+anp) and is otherwise inactive. therefore, the partitioning scheme activates exactly only one np∈{ …,np} so that the feasible region corresponding to the relaxation of xy is reduced from the parallelogram in figure to a significantly smaller one depicted in figure eq. enforces that only one binary variable is active: ∑np=1npλ= the continuous switch Δy takes on any positive value between and yu−yl when the binary switch corresponding to the npth piecewise λ is active =1) and otherwise. therefore: y=yl+∑np=1npΔy 0≤Δy≤λnp= …,np finally, the under and overestimators for the active segment are defined in algebraic terms as follows: z≥xyl+∑np=1npΔy z≥xyu+∑np=1np z≤xyl+∑np=1npΔy z≤xyu+∑np=1np xl≤x≤xu;yl≤y≤yu note that the discrete relaxation is tighter than the continuous one over the entire feasible region. the introduction of the binary variables required in the piecewise mccormick reformulation gives rise to a mixed-integer nonlinear programming problem, with the only nonlinearities appearing in the objective function. while this minlp is convex and can be easily solved to global optimality with standard minlp solvers, it is more convenient to linearize it in order to obtain an milp formulation, for which more efficient software packages exist. the section that follows explains how this is accomplished. hyper-planes underestimation the convex minlp can be further reformulated into an milp by replacing the objective function by a set of hyper-planes. for this, we define two new variables as zu,j′=z^u,j−z¯u,j and α≥z′u,j the quadratic terms are then approximated by 1st degree taylor series. that is, the square terms are replaced by l hyper-planes uniformly distributed between the maximum and minimum desired values of zu,j′ so that the objective function is reduced to a summation of quadratic terms as follows: minθ,ξ,z^u∑j∈jm∑u∈uαu,j αu,j≥z′02u,j,l+2z0′u,j,l∀u∈u∀j∈jm∀l∈l upper level slave problem a valid upper bound on the global optimum is obtained by optimizing the original nlp locally. this nlp is initialized using the solution provided by the milp as starting point. the solution of this nlp is used to tighten the milp, so the lower and upper bounds tend to converge as iterations proceed. algorithm steps the proposed algorithm comprises the following steps: set iteration count it = ub = ∞, lb = −∞ and tolerance error = tol. set it=it + solve the master problem milp. if the milp is infeasible, stop. otherwise, update the current lb making lb=maxit, where lbit is the value of the objective function of the milp in the itth iteration. solve the slave problem nlp. if the nlp is infeasible add one more piecewise term and hyper-plane to the master milp and go to step of the algorithm. otherwise, update the current ub making ub = minit, where ubit is the value of the objective function of the nlp in the itth iteration. calculate the optimality gap og as og=|ub−lb|ub. if og≤tol, then stop. the current ub is regarded as the global optimum within the desired tolerance. otherwise, add one more piecewise section and hyper-plane to the master milp and go to step of the algorithm. remarks: there are different methods to update the piecewise bilinear approximation. one possible strategy is to update it by dividing the active piecewise into two equal-length segments. the new hyper-plane term z0′u,j,l is added at the optimal solution of the milp in the previous iteration. the univariate convex and concave terms in the reformulated problem can be either approximated by the secant or by a piecewise univariate function similarly as done with the mccormick envelopes. our algorithm needs to be tuned prior to its application. this is a common practice in any optimization algorithm. in a previous publication, we studied the issue of defining the number of piecewise intervals and supporting hyper-planes in an optimal manner. in practice, however, the optimal number of piecewise terms and hyper-planes is highly dependent on the specific instance being solved, so it is difficult to provide general guidelines on this. the approach presented might lead to large computational burdens in large-scale models of complex biological systems. future work will focus on expediting our algorithm through the addition of cutting planes and the use of customized decomposition strategies. case studies we illustrate the performance of the proposed algorithm through its application to two challenging benchmark parameter estimation problems: the isomerisation of α-pinene and the inhibition of hiv proteinase. the objective in these problems is to obtain the set of values of the model parameters such that the model response is as close as possible to the experimental data. for comparison purposes we used the global optimization package baron. baron is a commercial software for solving nonconvex optimization problems to global optimality. baron combines constraint propagation, interval analysis, duality, and enhanced "branch and bound" concepts for efficient range reduction with rigorous relaxations constructed by enlarging the feasible region and/or underestimating the objective function. the interested readers have the possibility to evaluate this software on their own for free in this link: http://www.neos-server.org/neos/solvers/go:baron/gams.html. our algorithm was implemented in gams. using cplex. for the milps and snopt for the nlps subproblems. all the calculations were performed in a pc/amd athlon ii at ghz using a single core. data about the size of the models can be found in table case study 1: isomerisation of α-pinene in this first case study, five kinetic parameters describing the thermal isomerisation of α-pinene are estimated. the proposed reaction scheme for this process is depicted in figure in this homogeneous chemical reaction, α-pinene is thermally isomerised to dipentene and allo-ocimene, which in turn yields α- and β-pyronene and a dimer. this process was originally studied by fuguitt and hawkins, which carried out a single experiment reporting the experimental concentrations of the reactant and the four products measured at eight time intervals. hunter and mcgregor postulated first-order kinetics and proposed the following set of ode’s describing the dynamic process: dγ1dt=−γ dγ2dt=p1γ dγ3dt=p2γ1−γ3+p5γ dγ4dt=p3γ dγ5dt=p4γ3−p5γ γ0= t∈ rodriguez-fernandez et al. addressed this problem by applying a metaheuristic based on the scatter search method. this strategy does not offer any theoretical guarantee of convergence to the global optimum in a finite number of iterations. following our approach, the state variables were approximated by lagrange polynomials using three collocation points evaluated at the shifted roots of orthogonal legendre polynomials and defining five finite elements of equal length. the nonconvexities in the resulting residual equations are given by the bilinear terms θiξe,k, j which were relaxed using piecewise mccormick approximations as described previously. the objective function was underestimated using supporting hyper-planes. it is well known that the quality of the lower bound predicted by a relaxation strongly depends on the bounds imposed on its variables. hence bounds on collocation coefficients were tightened by performing a bound contraction procedure. particularly, tight lower and upper bounds were estimated for each collocation coefficient by maximizing and minimizing its value while satisfying the constraints contained in the master problem. this is a costly process. for this reason, it was only performed recursively times before the initialization of the algorithm. the milp was further tightened by adding the following constraint: ∑j∈jm∑u∈u2≤ which forces the model to find a solution better than the one obtained at the beginning of the search by locally minimizing the original nlp. furthermore, the parameter θi was allowed to take any value within the interval. the problem was solved with initial hyper-planes. an extra hyper-plane was added in each iteration, but the total number of piecewise terms was kept constant in order to keep the milp in a manageable size. a tolerance of 5% was set as termination criterion. for comparison purposes, we solved the same problem with the standard global optimization package baron using its default settings. baron was able to find the global optimum but failed at reducing the optimality gap below the specified tolerance after 12h of cpu time. in contrast, our algorithm closed the gap in less than 3h. as shown in table the results obtained agree with those reported in the literature. case study 2: inhibition of hiv proteinase in this second case study, we considered a much more complex biological dynamic system. particularly, we studied the reaction mechanism of the irreversible inhibition of hiv proteinase, as originally examined by kuzmic. note that this dynamic model has lack of practical identifiability, as reported in rodriguez-fernandez et al. nevertheless, we think that this example is still useful for the purpose of our analysis, since the emphasis here is placed on globally optimizing dynamic models of biological systems rather than analyzing identifiability issues. the model can be described mathematically through a set of nonlinear ode’s with ten parameters: ddt=−2k11+2k ddt=k3−2k41+2k ddt=−k21+k ddt=−k51+k ddt=k21−k22−k ddt=k41−k ddt=k11−k12−k21+k22+k3−k41+k42−k51+k ddt=k51−k52−k ddt=k where the following initial conditions and parameters are known: 0=00=00=00=00=00=00=00= = = = k11= k12= k41=100k21=100k51= t∈ a series of five experiments where the enzyme hiv proteinase was added to a solution of an irreversible inhibitor and a fluorogenic substrate were considered. the five experiments were carried out at four different concentrations of the inhibitor. the fluorescence changes were monitored during one hour. the measured signal is a linear function of the product concentration, as expressed in the following equation: signal=ϵ+offset in this fit, the offset of the fluorimeter was considered as a degree of freedom. a certain degree of uncertainty was assumed for the value of the initial concentrations of substrate and enzyme. the calibration of a total of adjustable parameters was addressed: five rate constants, five initial concentrations of enzyme and substrate and five offset values. mendes and kell solved this problem using simulated annealing and reported its first known solution. later, rodriguez-fernandez et al. improved that solution by means of a scatter search metaheuristic, which required a fraction of the time employed by mendes’ simulated annealing. recall that, despite producing near optimal solutions in short cpu times, stochastic algorithms provide no information on the quality of the solutions found and are unable to guarantee convergence to the global optimum in a finite number of iterations. on the contrary, the proposed methodology ensures the global optimality of the solution computed within a desired tolerance. in our study, the state variables were approximated using five orthogonal collocation points and five equal-length finite elements. in this case, the nonconvexities arise from the bilinear terms ξe, k, jξe, k, j and θiξe, k, j. the parameter bounds θi were set to θi∈. the lower and upper limits for the collocation coefficients ξe, k, j, n were fixed to ξe, k, j, n∈ except for ξe, k, e, n∈ and ξe, k, s, n∈. the bounds for all the offsets were set to offsetn∈. the master problem was further tightened by adding a special type of strengthening cuts. these cuts are generated by temporally decomposing the original full space milp into a series of milps in each of which we fit only a subset of the original dataset, and remove the continuity equations corresponding to the extreme elements included in the sub-problem. the cuts are expressed as inequalities added to the master problem that impose lower bounds on the error of a subset of elements for which the sub-milps are solved. these bounds are hence obtained from the solution of a set of milp sub-problems that optimize the error of only a subset of elements. this case study was solved with initial piecewise intervals and initial hyper-planes. two strengthening cuts involving elements and and and respectively, were added as constrains. a tolerance of 20% was used in the calculations. hyper-planes and piecewise terms were updated at each iteration of the algorithm. in this case, baron failed to identify any feasible solution after 12h of cpu time. in contrast, our algorithm was able to obtain the global optimum with a gap of % in approximately cpu s. remarkably, the solution found by our algorithm improves the best known solution reported by rodriguez-fernandez et al.. hence, our algorithm clearly outperformed other parameter estimation methods, improving the best known solution, and providing a rigorous lower bound on the minimum error that can be attained.  in this work, we have proposed a novel strategy for globally optimizing parameter estimation problems with embedded nonlinear dynamic systems. the method presented was tested through two challenging benchmark problems: the isomerisation of α-pinene and the inhibition of hiv proteinase. the proposed algorithm identified the best known solution, which was originally reported by rodriguez-fernandez et al., in the case of the α-pinene, and improved the best known one in the hiv proteinase case study. in both cases, rigorous lower bounds were provided on the global optimum, making it possible to determine the optimality gap of the solutions found. the method proposed produced promising results, surpassing the capabilities of baron. our method requires some knowledge on optimization theory as well as skills using modelling systems. our final goal is to develop a software to automate the calculations, so our approach can be easily used by a wider community. this is a challenging task, since nonlinear models are hard to handle and typically require customized solution procedures. particularly, nonlinear models must be initialized carefully to ensure convergence even to a local solution. in this regard, the use of an outer approximation scheme that relies on a master milp formulation is quite appealing, since the outcome of this milp can be used to initialize the nlp in a robust manner. another key point here is how to construct tight relaxations of the nonconvex terms. an efficient algorithm must exploit the problem structure to obtain high quality relaxations and therefore good bounds close to the global optimum. these relaxations can be further tightened through the addition of cutting planes or the use of customized decomposition methods. as observed, there is still much work to be done in this area, but we strongly believe that such an effort is worthy. furthermore, recent advances in global optimization theory and software applications are paving the way to develop systematic deterministic tools for the global optimization of parameter estimation problems of increasing size. our future work will focus on making the approach more efficient through the use of tailored cutting planes and decomposition strategies and also through the hybridization of deterministic methods with stochastic approaches. competing interest the authors declare that they have no competing interests. authors’ contributions g. g-g suggested the need for the approach and j.a. e provided the biological problems. a. m, c. p, g. g-g and l. j developed the optimization algorithms and performed the numerical analysis. all authors evaluated the results, wrote the paper and contributed to its final form. 
the estimation of parameter values for mathematical models of biological systems is an optimization problem that is particularly challenging due to the nonlinearities involved. in contrast, deterministic methods provide solutions of a given quality, but tend to lead to large computational burdens. the algorithm iterates between these two levels until a termination criterion is satisfied. this nonconvex nlp is decomposed into two hierarchical levels: a master mixed-integer linear programming problem that provides a rigorous lower bound on the optimal solution, and a reduced-space slave nlp that yields an upper bound. the proposed strategy produced near optimal solutions in a fraction of the cpu time required by baron. this work presents a deterministic outer approximation-based algorithm for the global optimization of dynamic problems arising in the parameter estimation of models of biological systems. the former typically lead to lower cpu times but offer no guarantee of convergence to the global minimum in a finite number of iterations. our approach, which offers a theoretical guarantee of convergence to global minimum, is based on reformulating the set of ordinary differential equations into an equivalent set of algebraic equations through the use of orthogonal collocation methods, giving rise to a nonconvex nonlinear programming problem. deterministic global optimization methods overcome this limitation, ensuring convergence to the global optimum within a desired tolerance. the capabilities of our approach were tested in two benchmark problems, in which the performance of our algorithm was compared with that of the commercial global optimization package baron. global optimization techniques are usually classified into stochastic and deterministic.
 the estimation of parameter values for mathematical models of biological systems is an optimization problem that is particularly challenging due to the nonlinearities involved. one major difficulty is the existence of multiple minima in which standard optimization methods may fall during the search. deterministic global optimization methods overcome this limitation, ensuring convergence to the global optimum within a desired tolerance. global optimization techniques are usually classified into stochastic and deterministic. the former typically lead to lower cpu times but offer no guarantee of convergence to the global minimum in a finite number of iterations. in contrast, deterministic methods provide solutions of a given quality, but tend to lead to large computational burdens.  this work presents a deterministic outer approximation-based algorithm for the global optimization of dynamic problems arising in the parameter estimation of models of biological systems. our approach, which offers a theoretical guarantee of convergence to global minimum, is based on reformulating the set of ordinary differential equations into an equivalent set of algebraic equations through the use of orthogonal collocation methods, giving rise to a nonconvex nonlinear programming problem. this nonconvex nlp is decomposed into two hierarchical levels: a master mixed-integer linear programming problem that provides a rigorous lower bound on the optimal solution, and a reduced-space slave nlp that yields an upper bound. the algorithm iterates between these two levels until a termination criterion is satisfied.  the capabilities of our approach were tested in two benchmark problems, in which the performance of our algorithm was compared with that of the commercial global optimization package baron. the proposed strategy produced near optimal solutions in a fraction of the cpu time required by baron. 