BACKGROUND
the development of mechanistic  models in order to quantitatively describe the dynamics of biological phenomena is one of the core research themes in systems biology.
during the last decade, fostered by the greater availability of the necessary experimental data, the development of large  kinetic models has become one of the main objectives in the field, as well as in related areas such as synthetic biology, metabolic engineering, or industrial biotechnology .
more recently, the first steps towards comprehensive whole-cell models have been taken  <cit> , which has great potential for applications e.g.
in personalized medicine  <cit> .
however, the development of these large-scale integrated dynamic models poses severe challenges  <cit> .
those associated with model building are common to the more general problem of reverse engineering of biological systems  <cit> .
in this context, parameter estimation  is arguably one of the most studied , yet more challenging step in model building.
parameter estimation in nonlinear dynamic models can be an extremely hard problem mostly due to the following issues  <cit> : lack of identifiability, ill-conditioning, multimodality and over-fitting.
the latter three can be handled via global optimization and regularization methods, as reviewed and illustrated recently  <cit> .
the present paper begins by continuing the line of work in  <cit> , addressing these three issues.
to this end we introduce a combination of a global optimization metaheuristic, ess  <cit> , and an efficient local search method, the adaptive algorithm nl2sol  <cit> .
by using this optimization technique jointly with regularization it is possible to reduce the calibration times of large dynamic models and simultaneously avoid over-fitting.
we show this for models from the recently presented biopredyn benchmark collection  <cit> .
then we focus on the remaining issue, that is, identifiability analysis of large dynamic models.
our aim is to develop a methodology which  is able to characterize high-order relationships among parameters, and  scales up well with model size.
thus, our objective goes beyond finding the subset of identifiable parameters: we also aim to systematically characterize the space of non-identifiable parameters, and to facilitate the advanced analysis of the results with scalable visualization tools.
identifiability analysis aims at establishing whether it is possible to determine the values of the unknown model parameters  <cit> .
it is common to distinguish between structural and practical identifiability.
structural or a priori identifiability analysis decides whether the model parameters are uniquely determinable based on the model formulation, which includes the dynamic equations, observation functions and stimuli  <cit> .
a parameter θ of the model is structurally identifiable if y=y⇔θ=θ
′, where y denotes the model predictions, which are observable in the experiments.
a parameter θ is structurally locally identifiable if for almost any value θ
∗ there is a neighbourhood v in which the above relationship holds.
it is globally identifiable if the relationship holds in all the range of values of the parameter.
if there is some region with non-zero measure where the relationship does not hold, θ is structurally unidentifiable.
structural identifiability analyses usually involve a high computational burden, which makes them difficult to apply to large models .
furthermore, structural identifiability is only a necessary but not sufficient condition for identifiability.
very often a structurally identifiable parameter is practically unidentifiable, that is, its value cannot be determined with precision due to limitations in the available data.
this can be quantified using practical or a posteriori identifiability analysis, which provides confidence intervals of the parameter values.
the two main sources of practical non-identifiability are  lack of influence of a parameter on the observables, and  interdependence among the parameters.
obviously, if a parameter does not influence the observables  it is not possible to determine its value.
the second situation, in which the effect on the observables of a change in one parameter can be compensated by changes in other parameters, can also prevent parameter identification.
both problems are related to the sensitivities of the observables to changes in model parameters.
while  is related to the average sensitivity of the model outputs to a specific parameter,  can be investigated based on the collinearity of the parametric sensitivities  <cit> .
in this paper we combine global optimization and regularization techniques to calibrate medium and large scale biological models .
then we evaluate the practical identifiability of estimated model parameters using sensitivity analysis and collinearity measures.
we determine the largest identifiable subsets of parameters, characterize the interplay among non-identifiable groups of parameters, and visualize the results using cytoscape.
the visualization tool shows the identifiable and non-identifiable groups of parameters together with the model structure in the same graph.
in this way, modellers can detect parts that are problematic and need refinement or reformulation, and experimentalists obtain information that can be helpful in the design of new experiments.
the methods for identifiability analysis and visualization presented here have been implemented as a matlab toolbox called visid, which is available from github  and as additional file  <dig>
methods
parameter estimation with regularization and global optimization
mathematical model
we consider deterministic models of biological systems that can be described by nonlinear ordinary differential equations  in the following form:
<dig> dxdt=f,u,θ),
<dig> y=g,θ),
<dig> x=x <dig> t∈t <dig> tf.
here x∈ℝnx denotes the state vector , f describes the interactions among the state variables , and u denotes the input variables .
the parameter vector θ∈ℝnθ contains the  parameters, e.g.
reaction rate coefficients or hill exponents.
their values are often unknown and must be estimated from data.
the model variables x are mapped to the measurable output variables y∈ℝny, also known as observables or model predictions, by the observation function g. these y signals are the quantities that can be experimentally measured.
we will denote by y
ijk the model prediction for the j-th observed quantity in the k-th experiment at time t
i∈[t
<dig>  t
f].
the corresponding measured data is denoted by y~ijk.
parameter estimation
the goal of parameter estimation is to determine the values of the unknown parameter vector θ. this is usually done by minimizing a distance between model prediction y
ijk and measured data y~ijk.
one of the simplest, but yet general, choices of this distance is the weighted sum-of-squares
<dig> qls=∑k=1ne∑j=1ny,k∑i=1nt,k,jwijkyijkx,θ−y~ijk <dig>
where n
e is the number of experiments, n
y,k is the number of observed compounds in the k-th experiment, and n
t,k,j is the number of measurement time points of the j-th observed quantity in the k-th experiment, and the weights are denoted by w
ijk.
the total number of data in all experiments is denoted by nd=∑k=1ne∑j=1ny,k∑i=1nt,k,j <dig>  in order to simplify the index triplet, from now on we will use only one index, i.e.
the weights and observables are denoted by w
i and y
i for i= <dig>  2…n
d.
then the parameter estimation problem is formulated as an optimization problem in the following form:
<dig> minimizeθqls+αΓ
<dig> subject toθmin≤θ≤θmax,
<dig> dxdt=f,x,θ),
<dig> y=g,θ),
<dig> x=x <dig> t∈t <dig> tf.
here Γ is a a regularization term, which is described in the following subsection, and θ
min and θ
max are lower and upper bounds of the parameter values.
the parameter vector θ^ that solves this minimization problem is called the optimal parameter vector or the parameter estimates.
regularization
large scale dynamic models are often over-parametrized, turning the estimation of their parameters into an ill-posed problem  <cit> .
this means that the minimum of the least-squares cost function  is non-unique, or that even a very small perturbation of the data results in very different estimated parameters.
furthermore, due to the large number of degrees of freedom, these models tend to capture the artificial dynamics of measurement noise.
this is known as overfitting  <cit>  and it usually results in poor predictive capability of the calibrated model.
regularization techniques incorporate a priori knowledge about the parameter values to make the problem well-posed.
the regularization parameter α in  balances the strength of this knowledge; its value can be found by regularization tuning methods  <cit> .
here we followed the guidelines presented in  <cit>  and chose a small regularization parameter , since we assume that we do not have good a priori estimates of the parameters.
regarding the regularization function, Γ, we chose the tikhonov regularization framework to match the form of the penalty to the least squares formalism of the objective function.
in this case the penalty is a quadratic penalty function,
<dig> Γ=θ−θreftwtwθ−θref,
where w∈ℝnθ×nθ is a diagonal scaling matrix and θref∈ℝnθ is a reference parameter vector, which is problem dependent and determined by the available information about the model parameters.
global optimization
we solve the minimization problem defined by – using optimization.
since the cost function  is usually multi-modal  , it is necessary to use an efficient global optimization method.
deterministic global optimization methods  can guarantee global optimality of the solution.
however, their computational cost increases exponentially with the number of parameters, which makes them unsatisfactory for large scale models.
stochastic and metaheuristic methods  <cit> , on the other hand, do not provide such guarantees, but are often capable of finding adequate solutions in reasonable computation times.
for this reason we use a method called enhanced scatter search   <cit> , which is an advanced implementation of a population-based algorithm called scatter search.
the scatter search metaheuristic works by evolving a number of solutions , which constitute the reference set .
members of this set are selected due to their quality and diversity.
they are updated at every iteration by combining them with other refset members and, occasionally, by applying an improvement method.
this improvement consists of a local search to speed-up the convergence to optimal solutions.
in the present work we have chosen nl2sol  <cit>  as a local method.
nl2sol is a quasi-newton algorithm with trust region strategy that exploits the structure of the nonlinear least squares problem.
note that the combination of a global method  with a local one makes ess a hybrid algorithm.
practical identifiability analysis
the shape of the cost function  in the surroundings of its optima determines the local identifiability of the parameters.
we assess parametric identifiability in two consecutive steps:
first we calculate the sensitivity of the model outputs  with respect to changes in the parameters.
those parameters which have no effect  on the observed signals are classified as non-identifiable.
note that this label is assigned on an individual basis, that is, taking only into account the effect of each parameter individually.
even if a parameter influences the model output, it may still be unidentifiable if its effect can be compensated by changes in other parameter.
hence in the second step we consider the interplay among parameters, aiming at finding groups of parameters which are non-identifiable due to their collinearity.
note that, while it would be possible at least in principle to perform both steps simultaneously, in practice the curse of dimensionality hampers the application of such a global sensitivity approach to large models  <cit> .
sensitivity analysis
the analysis of parametric sensitivity of kinetic models has a long tradition in model analysis  <cit> .
for the dynamical system –, the parametric sensitivities of the observables can be accurately calculated by solving the forward sensitivity equations:
<dig> dxidt=∂f∂xxi+∂f∂θfori= <dig> …,nθ
<dig> si=∂g∂xxi+∂g∂θfori= <dig> …,nθ
<dig> si=0ifθiis a model parameter1ifθiis an initial conditionfori= <dig> …,nθ.
here xi=∂x∂θi denotes the sensitivity of the state vector with respect to the i-th parameter and the vector si=∂y∂θi is the sensitivity of the observables with respect to this parameter.
this calculation requires the solution of the n
x×n
θ ordinary differential eq.
with initial conditions  for each experiment.
the numerical solution is determined for the time points for which there are experimental data available, and then the algebraic eq.
are evaluated.
if the partial derivatives of the dynamic equations are not available, an alternative is to calculate the sensitivities using finite differences or automatic differentiation.
the sensitivities of the observables are scaled using the same weights as in eq.
, resulting in scaled sensitivities for an output j and a parameter i:
<dig> s~ij=wj∂yj∂θi.
for each parameter we calculate an overall scoring called root mean squared sensitivity, s~imsqr, to take into account changes in time or across experiments  <cit> :
<dig> s~imsqr=1nd∑j=1nds~ij2fori= <dig> …,nθ.
below a certain threshold the parameters are considered non-influential to the outputs.
we set the threshold to four orders of magnitude smaller than the maximum root mean square value .
parameters whose sensitivity falls below this cut-off value are considered practically non-identifiable and they are kept out of further analysis.
the procedure is summarized in algorithm  <dig>
we remark that the outcome of the sensitivity calculations depends not only on the parameters, but also on the choice of initial conditions and external stimuli, which can have a strong influence in the practical identifiability of a model.
if insufficiently excitatory stimuli or initial conditions result in poor practical identifiability, a solution – if it is possible to carry out additional measurements – is to design and perform a new experiment to generate maximally informative data  <cit> .
collinearity of parameters
interplay among influential parameters can result in an unidentifiable model, because a variation in the cost function value due to a change in a parameter can be compensated by changes in other parameters.
pairwise interplay can be detected by plotting contours of the cost function versus pairs of parameters.
largely eccentric contours or “valleys” show that the cost function is almost unchanged in one direction, and the two parameters are highly correlated.
this approach has two drawbacks: it involves a large computational effort and is limited to interplay between pairs of parameters.
to compute higher dimensional interactions we use a different measure: the collinearity of parametric sensitivities.
to calculate collinearity we first normalize the scaled sensitivities  as follows:
<dig> s¯i=ŝi∥ŝi∥fori= <dig> …,nθ.
this normalization avoids biases caused by differences in the absolute values of the individual sensitivity vectors.
let us consider a set k of k parameters and their corresponding sensitivity vectors.
the parameters are linearly dependent if there exist k constants α
i≠ <dig> such that
<dig> α1s¯k1+α2s¯k2+…αks¯kk= <dig>
if the above relation does not hold, the set is independent.
when the equality  holds only approximately, the parameters are nearly dependent or nearly collinear.
the degree of collinearity among a set of parameters can be measured by the collinearity index, cik, which is defined as  <cit> :
<dig> cik=1min∥α∥=1∥s¯kα∥=1λk,min.
where s¯k is the sensitivity matrix built from the k sensitivity vectors, s¯k=, and λ
k,min is the smallest eigenvalue of s¯kts¯k.
the larger the collinearity index is, the more dependent the corresponding parameters are.
brun and co-authors  <cit>  proposed to classify a subset of parameters as identifiable if their collinearity index is smaller than a threshold which they chose as cik< <dig>  roughly speaking, a value of  <dig> means that 95% of the variation in the model output caused by changing one of the parameters in the subset can be compensated by changing the other parameters in the set.
other approaches for finding parameter correlations using sensitivity-based measures have been previously proposed in the literature.
li and vu presented two methods  <cit>  that search for relationships among parameters in the context of a priori identifiability analysis .
the method in  <cit>  provides a necessary but not sufficient condition for identifiability of nonlinear systems, which need to be fully observed .
the method in  <cit>  removes the requirement of measuring all the system states, replacing it with the restriction that the model must be linear.
we remark that the method proposed in the present manuscript does not have these limitations: it can be applied to partially observed, nonlinear systems with noisy, discrete-time measurements.
largest identifiable subset
as explained in the previous subsections, a subset of parameters is considered identifiable if its elements are influential and their sensitivity vectors are not collinear.
we are interested in finding the largest set of parameters for which the collinearity of the corresponding sensitivity vectors is below the chosen threshold, cik< <dig>  such a set of parameters represents all the degrees of freedom in the model.
this means that perturbing a parameter not included in this set has an effect in the model predictions that can be compensated  by changing other parameters in the set.
however, a perturbation in a parameter belonging to the set cannot be compensated by changes in the remaining parameters.
several methods have been developed for finding the group of identifiable parameters  <cit> .
iterative selection methods apply a step-wise procedure to select one parameter at a time, until no more parameters can be added to the identifiable set.
in each step the parameter to be included is selected based on an optimality criteria.
for example, the modified gram-schmidt orthogonalization method  <cit>  projects all the remaining sensitivity vectors to the subspace spanned by the already selected sensitivity vectors, and includes the parameter corresponding to the one with the largest projection value.
this step is repeated until the largest projection value falls below a threshold, which means that the next parameter would significantly interplay with the parameters already included.
the computational cost of this method scales up well with the number of sensitivity vectors.
however, the drawback of iterative procedures such as this one is that the solution might not be the global optimum, that is, it might fail to find the largest identifiable subset.
alternatively, we propose to solve the problem of finding the largest identifiable subset of all the estimated parameters using combinatorial optimization.
to this end we formulate it as a  integer optimization problem, where the goal is to maximize the number of sensitivity vectors included in the set, with the constraint that the corresponding collinearity index is below a threshold ci∗.
this algorithm can be stated as
<dig> maximizei∈{ <dig> }nθ∑k=1nθik
<dig> subject tosi=cat{sk|ik= <dig> fork= <dig> …,nθ}
<dig> ci<ci∗
<dig> ikis a binary variable fork= <dig> …,nθ
where the binary variable i
k indicates if the k-th parameter is included  or not included  in the identifiable group of parameters.
the sensitivity matrix corresponding to the selected parameters is s
i, and ‘cat’ stands for the concatenation of the column vectors in the constraint .
the collinearity index of this matrix is ci and it is determined by computing the minimum eigenvalue as in .
this combinatorial optimization problem has an exponentially scaling computational cost, and thus its solution requires an efficient algorithm.
we chose the variable neighbourhood search  technique  <cit> , which is a heuristic global optimization method for integer optimization problems.
we used the version of vns included in the meigo toolbox  <cit> , which is implemented in matlab.
we modified this initial formulation of the problem described in eqs.
– after finding that its solution is often not unique: even after maximizing the number of parameters in the subset, there may be multiple subsets that yield a collinearity index below the threshold ci∗.
indeed, we found large variability in the solutions if no initial guess was specified.
therefore, we reformulated the optimization problem in two ways, as described in the following paragraphs.
as a first modification, we transformed the collinearity requirement  from a ‘hard’ to a ‘soft’ constraint .
the modified optimization problem reads as
<dig> maximizei∈{ <dig> }nθ∑k=1nθik−p1−p <dig>
<dig> subject tosi=catsk|ik= <dig> fork= <dig> …,nθ
<dig> p1=12ci/ci∗
<dig> p2=0ifci<ci∗αci−ci∗βotherwise
<dig> ikis a binary variable fork= <dig> …,nθ
as above, the binary variable i
k indicates if the k-th parameter is included  or not included  in the selected group of parameters.
the penalty p
<dig> is a monotone increasing  function of the collinearity index ci, such that p
<dig> is  <dig>  when the collinearity equals to the threshold.
due to this small value, p
<dig> does not influence the size of the largest subset below the threshold.
in this way, when multiple sets of the same size co-exist, the set with smaller collinearity index is always favoured.
this results in an unique solution of the optimization problem if there are no sets with identical collinearity index.
the second penalty function p
<dig> represents a soft constraint that is active when the collinearity exceeds the threshold.
the steepness of this constraint is tuned by the values of α and β, which we set to α= <dig> and β= <dig>
our second improvement of the formulation of the optimization problem consists in providing a good initial guess of the solution using qr decomposition.
the rank revealing qr decomposition algorithm, or rrqr  <cit> , rewrites a matrix s as
<dig> Πs=qr,
where q is an orthogonal matrix, r is an upper triangular matrix, and Π is a permutation matrix.
due to the properties of this decomposition, the permutation matrix defines a reordering of the columns of s. in this re-ordered matrix s
ro=Π
s, the most orthogonal columns are located in the left.
in other words, the first n columns of the reordered matrix define a linear subspace, and the -th column has the largest projection value on this subspace among the remaining n
θ−n columns located to the right of the n-th column.
the outcome of the rrqr technique is similar to that of the aforementioned gram-schmidt orthogonalization method, but its implementation is more efficient.
we applied rank revealing qr decomposition to the sensitivity matrix, following the procedure described in algorithm  <dig>  then, we used the resulted ordering of the sensitivity vectors to initialize the global optimizer.
in this way we improved the performance of the global optimizer, which often found larger sets with collinearity index below the threshold value.
the whole procedure for identifying the largest non-collinear subset of parameters is summarized in algorithm  <dig>
the procedure presented in this subsection has similarities with the one proposed by chu and hahn  <cit> .
one difference is that we maximize the subset size for a given collinearity threshold, whereas chu and hahn adopted the opposite approach, i.e., maximizing parametric identifiability for a pre-specified subset size.
additionally, both methods differ in the optimization technique: we use variable neighbourhood search, which has better scalability than the genetic algorithm chosen in  <cit> .
recently, nienałtowski et al.
<cit>  have proposed a method for finding clusters of correlated parameters using so-called canonical correlation analysis .
cca is an extension of pearson correlation for measuring multidimensional correlations between groups of parameters.
given two groups of parameters of sizes m and n, with m<n, calculation of the canonical correlations provides m measures, which are summarized in a single measure, called mi-cca.
this similarity measure represents the mutual information between the two groups, although it should be noted that average mutual information is equivalent to canonical correlation only if the random variables follow an elliptically symmetric probability model.
nienałtowski et al.
use mi-cca to cluster parameters until an identifiable subset is reached.
this approach is sequential and yields a single parameter subset, which is possibly not maximal.
in contrast, the methodology described here combines an initial sequential phase with a subsequent combinatorial optimization procedure.
the second phase yields several identifiable parameter subsets and usually improves the initial solution.
finding all largest subsets
as mentioned above, the largest non-collinear subset of parameters is not unique.
to realize this, imagine that we have a non-collinear set of the parameters, and consider an additional pair of highly collinear parameters.
since we may add either of these two parameters to the set, but not both of them, we have two potential solutions.
the optimization algorithm described above would choose the option with a lower collinearity index.
however, we may also be interested in enumerating all the possible sets, instead of only one.
finding all the largest subsets is a combinatorial problem too, which is computationally expensive.
a naive approach for solving it could be to generate all possible sets of parameters and compute the corresponding collinearity index.
however, note that if two parameters θ
<dig> and θ
<dig> are collinear, then any sets including the pair { θ
<dig> θ
2} are highly collinear.
using this fact, we developed an incremental procedure for the systematic determination of the sets.
we start by considering all possible pairs of parameters and determining their collinearity.
then we extend only those pairs which have a small collinearity index, by considering all possible combinations of a third parameter.
this procedure is repeated until either all the sets are highly collinear, or there is only one set containing all the parameters.
in this way, summarized in algorithm  <dig>  we can find all the largest subsets of non-collinear parameters.
partitioning the non-identifiable parameters
the two procedures presented above can be used for finding  the largest, least collinear subset of parameters, and  all the largest subsets; in both cases, restricted to those subsets whose collinearity falls below a threshold.
however, it is often important to understand why certain parameters are not identifiable.
for example, a parameter may be unidentifiable because the model output has very low sensitivity to changes in its value.
but it could also be because it is highly correlated with another parameter, even when both parameters have high sensitivities.
finding small groups of highly collinear parameters can be helpful in determining the exact source of unidentifiability.
the collinearity of a subset always increases when a new parameter is added to the set.
for example, considering three parameters, the collinearity of the triplet is always higher than the collinearity of any pairs.
therefore, if a larger set of parameters contains a collinear pair, then the collinearity index of the large set is also large.
if we are interested in finding the smallest groups of highly collinear parameters, we can proceed as follows.
first we generate all possible pairs of parameters, and compute the collinearity of the corresponding sensitivity vectors.
then we evaluate all possible triplets.
the procedure can be extended for the analysis of larger sets.
however, due to the combinatorial explosion of the computational cost, this method can be applied only to models of moderate size .
visualization of identifiable subsets
it can be useful to represent the identifiability results graphically, because such visualization can provide modellers with insight about how to reformulate their models and/or design new experiments in order to avoid non-identifiable parameters.
with this aim, we display the model structure in the natural network visualization technique.
an example is shown in fig.
1
c. the model structure is represented as a graph whose nodes are state variables, observables, stimuli, and model parameters.
the edges – which can be directed  or undirected – have the following meaning: an arrow from node a to node b indicates that node b appears in the equation of a. for example, if the dynamic equation of a state x
<dig> is x˙1=p1·x <dig>  the corresponding graph would show two arrows x
2→x
<dig> and p
1→x
<dig>
fig.
<dig> tgf- β model: panel a shows the model sensitivities with respect to the parameters in logarithmic scale.
in panel b the maximum number of identifiable parameters is depicted depending on the collinearity threshold.
the group-size does not change much by the threshold.
the red vertical line indicates our choice  for the further analysis.
in panel c the mathematical model is depicted.
nodes indicate states , identifiable  and not identifiable  parameters, and observables .
in panel d the interplay among the collinear parameters are indicated.
there are  <dig> groups of parameters: in the first  <dig> groups triplets of parameters show large collinearity, while in the fifth group  <dig> parameters are collinear
more formally, we determine how the state, input variables, stimuli and parameters are connected and influence each other through symbolic manipulation of the model eq.
.
for this purpose we compute:  the jacobian matrix with respect to the states: ji,jss=∂fi∂xj,  the jacobian of the observation functions with respect to the states: ji,jso=∂gi∂xj,  the jacobian of the systems dynamics with respect to the stimuli ji,jsi=∂fi∂uj, and  the jacobian with respect to the parameters ji,jsp=∂fi∂θj.
all these matrices are evaluated symbolically, and then the expressions are converted to a logical  <dig>  or  <dig> when the symbolic result is zero.
additionally, we can connect parameters by undirected edges if their collinearity is larger than the collinearity threshold.
implementation: the visid software tool
we implemented the techniques proposed in subsections “practical identifiability analysis” and “visualization of identifiable subsets” as a matlab software package called visid, which is provided as additional file  <dig> and can also be downloaded from github .
it is free software, made available under the terms of the gnu general public license version  <dig>  the visid toolbox relies on three other matlab toolboxes, which are also freely available: amigo <dig>  <cit>  , which is used to to store, simulate and calibrate the models; meigo  <cit>  , which implements the variable neighbouring search  optimization method; and  rrqr , which performs the rank revealing qr decomposition used to initialize the global optimizer.
network visualization is performed with cytoscape  <cit>  .
further details can be found in section  <dig> of additional file  <dig>
RESULTS
in this section we demonstrate the application of the methodology presented in the previous section using several dynamic systems biology models of different type and complexity.
their main characteristics are given in table  <dig>  first we present detailed results of identifiability analysis and visualization for a model of the tgf- β signalling pathway.
we also provide similar results for the genetic network that controls the circadian clock in arabidopsis thaliana.
due to their complexity and yet relatively moderate size, these models are well suited as case studies for illustrating the identifiability methodology in depth.
a. thaliana
then we study two large scale benchmark problems included in the biopredyn-bench collection  <cit> .
since the analysis of these latter models is more challenging due to their larger size, we start by demonstrating the performance improvements that can be achieved during parameter estimation using the model calibration procedure proposed in section “parameter estimation with regularization and global optimization”.
then we perform identifiability analysis and report the corresponding results, including the graphical representation of the identifiable subset using the natural network visualization.
tgf- β signalling pathway
the dynamic model of the tgf- β signaling pathway was presented in  <cit>  as a tutorial example for model calibration.
it has  <dig> dynamic states and  <dig> kinetic parameters , of which  <dig> need to be estimated.
following  <cit> , we assumed that all the concentrations, except the smad rnas , can be measured in the experiments.
the algebraic equations of the reaction kinetics and the dynamic equations are provided in the additional file  <dig>
for the purpose of testing the methodology we generated a training dataset by simulating the model equations using the nominal values of the parameters k
1– k
<dig> .
then we sampled the simulated trajectories at equidistant time points, and added normally distributed random numbers to the data to mimic measurement errors.
finally, we estimated the model parameters from the generated data set.
this approach is widely used for testing calibration methods and assessing the extent to which they recover the nominal parameters.
it should be noted that, as the amount of noise in the dataset increases, the information/signal ratio decreases, making the estimation problem more ill-conditioned.
this makes it more difficult to recover the correct value of the parameters, but has a small effect in computation times.
the numerical values of the estimated parameters are reported in additional file 2: table s <dig>
we started the identifiability analysis by computing the sensitivities of the observations with respect to the estimated model parameters, according to algorithm  <dig>  we found that all the parameters have a non-negligible influence on the model outputs, thus there are no individually non-identifiable parameters .
next, following algorithm  <dig>  we applied qr decomposition and ranked the parameters according to their orthogonality.
we then solved the optimization problem – by initializing the variable neighboring search method with the results of the qr decomposition .
setting the threshold level for the collinearity index to ci = <dig> yielded  <dig> identifiable parameters, which are shown as green nodes in the network in fig.
1
c. parameters not present in the identifiable subset are shown as red nodes.
parameters are connected by arrows to state variables  if they appear in the equation of the corresponding dynamic equation.
states which directly influence each other are also connected by directed edges in the same manner.
blue squares represent measurements; a state is connected to a blue square if it appears in the corresponding observation function.
to see how the size of the identifiable subset is influenced by the choice of the collinearity index threshold , we solved the optimization problem for a range of threshold values.
the results are depicted in fig.
1
b. as the collinearity index threshold decreases, less parameters are considered identifiable.
we can see that the identifiability results are quite robust to the choice of threshold level: the number of identifiable parameters is always between  <dig> and  <dig>  and it is constant  for a very wide range of ci, 15≤ci≤ <dig>
the results presented so far tell us that the  <dig> parameters are not correlated.
however, they do not inform of the relationships among identifiable and non-identifiable parameters.
to investigate this point, we computed the smallest correlated subsets as described in section “partitioning the non-identifiable parameters”, up to groups of  <dig> parameters.
figure 1
d shows such groups; parameters are depicted as blue circles connected with group identifying nodes .
these nodes are labeled as gx, where x indicates the number of parameters in the group and y is the group index for a given number of parameters  stands for the second group of three correlated parameters).
we found that the large pairwise collinearity between k
14−k
<dig> and k
17−k
<dig> explains the non-identifiability of the model parameters only partially.
there are  <dig> groups of triplets and a group of  <dig> parameters which are highly correlated.
the members of the groups and the corresponding collinearity index are reported in table  <dig> d
it is important to note that collinearity might arise among multiple parameters, even if they are pairwise independent.
for example, despite the fact that none of the pairs in the group of k
<dig>  k
<dig> and k
<dig> has a high pairwise collinearity, the collinearity index of the triplet is extremely large.
algorithm  <dig> found  <dig> different sets of identifiable parameters with collinearity index ranging between  <dig>  and  <dig>  less than the threshold .
the sets are reported with the corresponding collinearity index in additional file 2: table s <dig>  we can see that parameters {k
1−k
<dig>  k
<dig>  k
15} are members of all the groups, and they do not participate in any of the small correlated groups in fig.
1
d. from each correlated group of size k, only k− <dig> parameters can participate in the largest set of identifiable parameters.
the aforementioned identifiability procedures can be carried out in a few seconds.
detailed computational costs are shown in table s <dig> of the additional file  <dig> for all the case studies considered in this paper.
circadian clock in arabidopsis thaliana
locke and co-authors  <cit>  described the genetic network controlling the circadian clock in arabidopsis thaliana; the dynamic equations of this model are provided in the additional file  <dig>
we generated training data by simulating the model equations with the nominal parameters  in two experimental conditions.
in the first one, the model input was kept constant , representing continuous light stimulation of the plant.
in the second experiment the input was changed pulse-wise in 12-hour cycles, repeated  <dig> times.
as in the previous example, the trajectories were sampled at equidistant time-points and disturbed by pseudo-random noise.
only two states, c
t
m and c
l
m, were observed.
the estimated model parameters are collected in additional file 2: table s <dig>
although the model outputs showed sensitivity to all the parameters , i.e.
there were no zero sensitivity vectors, we found that most of the model parameters are non-identifiable due to heavy collinearities.
the largest identifiable subset contains only  <dig> of the  <dig> parameters, depicted in fig.
2
d by green nodes.
the enumeration of the largest sets of identifiable parameters by algorithm  <dig> identified  <dig> parameter sets.
fig.
<dig> circadian clock in arabidopsis thaliana: panel a shows the model sensitivities with respect to parameters in logarithmic scale.
in panel b the maximum number of identifiable parameters is depicted as a function of the collinearity threshold.
the red vertical line indicates our choice of threshold  used for further analysis.
in panel c the interplay among the collinear parameters is indicated up to groups of  <dig> parameters.
in panel d the mathematical model is depicted.
nodes indicate states , identifiable  and not identifiable  parameters, and observables .
highly correlated parameter pairs are connected by undirected edges
benchmarks b <dig> and b <dig> from the biopredyn-bench collection
in this subsection we analyze two large scale benchmark problems taken from the biopredyn-bench collection  <cit> : the metabolic models of escherichia coli  and chinese hamster ovary cells .
they are highly non-linear, partially observed systems with more than  <dig> unknown parameters, which pose serious challenges for parameter identification.
in b <dig> the calibration data was generated by model simulation and disturbed by random noise, while in b <dig> it was experimentally measured.
further details about the models and the parameter estimation challenge can be found in  <cit> .
first we use these benchmarks to illustrate the benefits of the parameter estimation strategy proposed in subsection “parameter estimation with regularization and global optimization”, comparing it with the one used in  <cit> .
both approaches use a hybrid method, ess  <cit> , which combines a global optimization algorithm  with a local search.
in  <cit>  the local method of choice was fmincon; here we compare that configuration with nl2sol .
global optimization algorithms use pseudo-random numbers.
hence their performance changes at every run, and the calibration problem should be solved several times to obtain more robust results.
since each optimization takes several hours we limited the number of runs to five for each problem.
we used the approximate computation time  reported in  <cit>  as the stopping criterion for the model calibration.
convergence curves depict the best objective function value found versus cpu time, and can be used to compare the performance of different algorithms.
an optimization method is preferred if it achieves a lower objective function value at earlier cpu time.
the best convergence curves  corresponding to b <dig> and b <dig> are shown in the additional file  <dig> for  <dig> algorithms:  ess-fmincon, as reported in  <cit> ;  ess-nl2sol; and  ess-nl2sol using regularization as recommended in section “parameter estimation with regularization and global optimization”.
from those curves we see that the algorithm  proposed here converged earlier than the others to the optimal objective function value .
we stress that the main purpose of regularization is to avoid overfitting: we do not wish to obtain an excessively good fit, which would indicate that we are reproducing noise instead of the true dynamics.
therefore, regularization should not achieve a smaller objective function value.
next, we apply the identifiability analysis procedures presented in subsections “practical identifiability analysis” and “visualization of identifiable subsets” to these two models.
for b <dig> they yield an identifiable subset of size  <dig>  and for b <dig> of size  <dig> .
the corresponding networks are shown in figs.
<dig> and  <dig>  respectively.
it is also possible to find small groups of highly correlated parameters for models of this size; e.g.
for b <dig> we obtained those depicted in fig.
<dig>
fig.
<dig> representation of the connections in the b <dig> model using the network diagram formalism.
nodes indicate states , identifiable  and not identifiable  parameters, observables , and inputs .
the source file of this figure is provided with the visid toolbox; using cytoscape, the user can navigate through it and zoom on different areas to improve the visibility
fig.
<dig> representation of the connections in the b <dig> model using the network diagram formalism.
nodes indicate states , identifiable  and not identifiable  parameters, observables , and inputs
fig.
<dig> visualization of the relationships among highly collinear parameters in the b <dig> model.
the figure shows small groups, whose sizes range between  <dig> and  <dig> parameters.
unidentifiable parameters are shown in red; identifiable parameters in green.
highly correlated pairs are connected by lines
the aforementioned results show that both models are poorly identifiable in practice for the considered datasets; more informative data would be needed in order to obtain accurate estimates of their parameters.
discussion and
CONCLUSIONS
in this paper we have presented a workflow to efficiently estimate the parameters of dynamic models and analyze their practical identifiability.
our approach combines an advanced optimization technique, which reduces computation times in parameter estimation, and several identifiability analysis procedures, which can find subsets of identifiable and unidentifiable parameters.
results are visualized using network diagrams, which provide an intuitive representation of the findings and facilitate their analysis and understanding.
many approaches have been applied to study identifiability of kinetic models, but they suffer from lack of scalability.
an advantage of the integrated method presented here is its moderate computational cost, which enables its application to large-scale models; complete results can be obtained in a few hours for models of more than a hundred parameters.
another important aspect is the integration of identifiability analysis with visualization, which presents the results in a way that is easily interpretable for modelers and experimentalists.
currently, its main limitation arises when trying to find all the different existing groups of highly correlated parameters: the combinatorial explosion of this particular task makes it feasible only for models of moderate size, i.e.
of a few dozens of parameters.
however, all the remaining steps of the workflow presented in this manuscript scale up well up to several hundred parameters.
the usefulness of the methodology and workflow presented here goes beyond basic parameter identifiability analysis.
the procedure not only  determines the largest subset of identifiable parameters, but also  informs about the characteristics of the space of non-identifiable parameters, reporting small groups of highly correlated parameters, and  presents all these results in a coherent and scalable way using visualization techniques, facilitating the understanding of the underlying complex interactions.
uncovering these higher order relationships helps in determining the causes of unidentifiability and provides guidelines for remedying them, e.g.
by reformulating the model or by collecting new data through a new experimental design.
all this information can be readily used to improve the iterative model-building cycle.
a matlab implementation of the identifiability and visualization methodology, which we have called the visid software package , is available from github  as free, open source software.
this distribution includes the case studies discussed above.
additional files
additional file  <dig> visid toolbox.
this compressed folder contains the visid matlab toolbox.
additional file  <dig> supplementary material.
this document contains detailed descriptions of the case studies and of the visid toolbox, as well as additional details about the results.
abbreviations
ccacanonical correlation analysis
cicollinearity index
cpucentral processing unit
essenhanced scatter search
mimutual information
odeordinary differential equation
qrdecomposition of a matrix into an orthogonal matrix  and an upper triangular matrix
rnaribonucleic acid
rrqrrank revealing qr decomposition
tgf- βtransforming growth factor beta
vnsvariable neighbourhood search

 the development of mechanistic models in order to quantitatively describe the dynamics of biological phenomena is one of the core research themes in systems biology. during the last decade, fostered by the greater availability of the necessary experimental data, the development of large kinetic models has become one of the main objectives in the field, as well as in related areas such as synthetic biology, metabolic engineering, or industrial biotechnology. more recently, the first steps towards comprehensive whole-cell models have been taken, which has great potential for applications e.g. in personalized medicine. however, the development of these large-scale integrated dynamic models poses severe challenges. those associated with model building are common to the more general problem of reverse engineering of biological systems. in this context, parameter estimation is arguably one of the most studied, yet more challenging step in model building. parameter estimation in nonlinear dynamic models can be an extremely hard problem mostly due to the following issues: lack of identifiability, ill-conditioning, multimodality and over-fitting. the latter three can be handled via global optimization and regularization methods, as reviewed and illustrated recently. the present paper begins by continuing the line of work in, addressing these three issues. to this end we introduce a combination of a global optimization metaheuristic, ess, and an efficient local search method, the adaptive algorithm nl2sol. by using this optimization technique jointly with regularization it is possible to reduce the calibration times of large dynamic models and simultaneously avoid over-fitting. we show this for models from the recently presented biopredyn benchmark collection. then we focus on the remaining issue, that is, identifiability analysis of large dynamic models. our aim is to develop a methodology which is able to characterize high-order relationships among parameters, and scales up well with model size. thus, our objective goes beyond finding the subset of identifiable parameters: we also aim to systematically characterize the space of non-identifiable parameters, and to facilitate the advanced analysis of the results with scalable visualization tools. identifiability analysis aims at establishing whether it is possible to determine the values of the unknown model parameters. it is common to distinguish between structural and practical identifiability. structural or a priori identifiability analysis decides whether the model parameters are uniquely determinable based on the model formulation, which includes the dynamic equations, observation functions and stimuli. a parameter θ of the model is structurally identifiable if y=y⇔θ=θ ′, where y denotes the model predictions, which are observable in the experiments. a parameter θ is structurally locally identifiable if for almost any value θ ∗ there is a neighbourhood v in which the above relationship holds. it is globally identifiable if the relationship holds in all the range of values of the parameter. if there is some region with non-zero measure where the relationship does not hold, θ is structurally unidentifiable. structural identifiability analyses usually involve a high computational burden, which makes them difficult to apply to large models. furthermore, structural identifiability is only a necessary but not sufficient condition for identifiability. very often a structurally identifiable parameter is practically unidentifiable, that is, its value cannot be determined with precision due to limitations in the available data. this can be quantified using practical or a posteriori identifiability analysis, which provides confidence intervals of the parameter values. the two main sources of practical non-identifiability are lack of influence of a parameter on the observables, and interdependence among the parameters. obviously, if a parameter does not influence the observables it is not possible to determine its value. the second situation, in which the effect on the observables of a change in one parameter can be compensated by changes in other parameters, can also prevent parameter identification. both problems are related to the sensitivities of the observables to changes in model parameters. while is related to the average sensitivity of the model outputs to a specific parameter, can be investigated based on the collinearity of the parametric sensitivities. in this paper we combine global optimization and regularization techniques to calibrate medium and large scale biological models. then we evaluate the practical identifiability of estimated model parameters using sensitivity analysis and collinearity measures. we determine the largest identifiable subsets of parameters, characterize the interplay among non-identifiable groups of parameters, and visualize the results using cytoscape. the visualization tool shows the identifiable and non-identifiable groups of parameters together with the model structure in the same graph. in this way, modellers can detect parts that are problematic and need refinement or reformulation, and experimentalists obtain information that can be helpful in the design of new experiments. the methods for identifiability analysis and visualization presented here have been implemented as a matlab toolbox called visid, which is available from github and as additional file methods parameter estimation with regularization and global optimization mathematical model we consider deterministic models of biological systems that can be described by nonlinear ordinary differential equations in the following form: dxdt=f,u,θ), y=g,θ), x=x t∈t tf. here x∈ℝnx denotes the state vector, f describes the interactions among the state variables, and u denotes the input variables. the parameter vector θ∈ℝnθ contains the parameters, e.g. reaction rate coefficients or hill exponents. their values are often unknown and must be estimated from data. the model variables x are mapped to the measurable output variables y∈ℝny, also known as observables or model predictions, by the observation function g. these y signals are the quantities that can be experimentally measured. we will denote by y ijk the model prediction for the j-th observed quantity in the k-th experiment at time t i∈[t t f]. the corresponding measured data is denoted by y~ijk. parameter estimation the goal of parameter estimation is to determine the values of the unknown parameter vector θ. this is usually done by minimizing a distance between model prediction y ijk and measured data y~ijk. one of the simplest, but yet general, choices of this distance is the weighted sum-of-squares qls=∑k=1ne∑j=1ny,k∑i=1nt,k,jwijkyijkx,θ−y~ijk where n e is the number of experiments, n y,k is the number of observed compounds in the k-th experiment, and n t,k,j is the number of measurement time points of the j-th observed quantity in the k-th experiment, and the weights are denoted by w ijk. the total number of data in all experiments is denoted by nd=∑k=1ne∑j=1ny,k∑i=1nt,k,j in order to simplify the index triplet, from now on we will use only one index, i.e. the weights and observables are denoted by w i and y i for i= 2…n d. then the parameter estimation problem is formulated as an optimization problem in the following form: minimizeθqls+αΓ subject toθmin≤θ≤θmax, dxdt=f,x,θ), y=g,θ), x=x t∈t tf. here Γ is a a regularization term, which is described in the following subsection, and θ min and θ max are lower and upper bounds of the parameter values. the parameter vector θ^ that solves this minimization problem is called the optimal parameter vector or the parameter estimates. regularization large scale dynamic models are often over-parametrized, turning the estimation of their parameters into an ill-posed problem. this means that the minimum of the least-squares cost function is non-unique, or that even a very small perturbation of the data results in very different estimated parameters. furthermore, due to the large number of degrees of freedom, these models tend to capture the artificial dynamics of measurement noise. this is known as overfitting and it usually results in poor predictive capability of the calibrated model. regularization techniques incorporate a priori knowledge about the parameter values to make the problem well-posed. the regularization parameter α in balances the strength of this knowledge; its value can be found by regularization tuning methods. here we followed the guidelines presented in and chose a small regularization parameter, since we assume that we do not have good a priori estimates of the parameters. regarding the regularization function, Γ, we chose the tikhonov regularization framework to match the form of the penalty to the least squares formalism of the objective function. in this case the penalty is a quadratic penalty function, Γ=θ−θreftwtwθ−θref, where w∈ℝnθ×nθ is a diagonal scaling matrix and θref∈ℝnθ is a reference parameter vector, which is problem dependent and determined by the available information about the model parameters. global optimization we solve the minimization problem defined by – using optimization. since the cost function is usually multi-modal, it is necessary to use an efficient global optimization method. deterministic global optimization methods can guarantee global optimality of the solution. however, their computational cost increases exponentially with the number of parameters, which makes them unsatisfactory for large scale models. stochastic and metaheuristic methods, on the other hand, do not provide such guarantees, but are often capable of finding adequate solutions in reasonable computation times. for this reason we use a method called enhanced scatter search, which is an advanced implementation of a population-based algorithm called scatter search. the scatter search metaheuristic works by evolving a number of solutions, which constitute the reference set. members of this set are selected due to their quality and diversity. they are updated at every iteration by combining them with other refset members and, occasionally, by applying an improvement method. this improvement consists of a local search to speed-up the convergence to optimal solutions. in the present work we have chosen nl2sol as a local method. nl2sol is a quasi-newton algorithm with trust region strategy that exploits the structure of the nonlinear least squares problem. note that the combination of a global method with a local one makes ess a hybrid algorithm. practical identifiability analysis the shape of the cost function in the surroundings of its optima determines the local identifiability of the parameters. we assess parametric identifiability in two consecutive steps: first we calculate the sensitivity of the model outputs with respect to changes in the parameters. those parameters which have no effect on the observed signals are classified as non-identifiable. note that this label is assigned on an individual basis, that is, taking only into account the effect of each parameter individually. even if a parameter influences the model output, it may still be unidentifiable if its effect can be compensated by changes in other parameter. hence in the second step we consider the interplay among parameters, aiming at finding groups of parameters which are non-identifiable due to their collinearity. note that, while it would be possible at least in principle to perform both steps simultaneously, in practice the curse of dimensionality hampers the application of such a global sensitivity approach to large models. sensitivity analysis the analysis of parametric sensitivity of kinetic models has a long tradition in model analysis. for the dynamical system –, the parametric sensitivities of the observables can be accurately calculated by solving the forward sensitivity equations: dxidt=∂f∂xxi+∂f∂θfori= …,nθ si=∂g∂xxi+∂g∂θfori= …,nθ si=0ifθiis a model parameter1ifθiis an initial conditionfori= …,nθ. here xi=∂x∂θi denotes the sensitivity of the state vector with respect to the i-th parameter and the vector si=∂y∂θi is the sensitivity of the observables with respect to this parameter. this calculation requires the solution of the n x×n θ ordinary differential eq. with initial conditions for each experiment. the numerical solution is determined for the time points for which there are experimental data available, and then the algebraic eq. are evaluated. if the partial derivatives of the dynamic equations are not available, an alternative is to calculate the sensitivities using finite differences or automatic differentiation. the sensitivities of the observables are scaled using the same weights as in eq., resulting in scaled sensitivities for an output j and a parameter i: s~ij=wj∂yj∂θi. for each parameter we calculate an overall scoring called root mean squared sensitivity, s~imsqr, to take into account changes in time or across experiments: s~imsqr=1nd∑j=1nds~ij2fori= …,nθ. below a certain threshold the parameters are considered non-influential to the outputs. we set the threshold to four orders of magnitude smaller than the maximum root mean square value. parameters whose sensitivity falls below this cut-off value are considered practically non-identifiable and they are kept out of further analysis. the procedure is summarized in algorithm we remark that the outcome of the sensitivity calculations depends not only on the parameters, but also on the choice of initial conditions and external stimuli, which can have a strong influence in the practical identifiability of a model. if insufficiently excitatory stimuli or initial conditions result in poor practical identifiability, a solution – if it is possible to carry out additional measurements – is to design and perform a new experiment to generate maximally informative data. collinearity of parameters interplay among influential parameters can result in an unidentifiable model, because a variation in the cost function value due to a change in a parameter can be compensated by changes in other parameters. pairwise interplay can be detected by plotting contours of the cost function versus pairs of parameters. largely eccentric contours or “valleys” show that the cost function is almost unchanged in one direction, and the two parameters are highly correlated. this approach has two drawbacks: it involves a large computational effort and is limited to interplay between pairs of parameters. to compute higher dimensional interactions we use a different measure: the collinearity of parametric sensitivities. to calculate collinearity we first normalize the scaled sensitivities as follows: s¯i=ŝi∥ŝi∥fori= …,nθ. this normalization avoids biases caused by differences in the absolute values of the individual sensitivity vectors. let us consider a set k of k parameters and their corresponding sensitivity vectors. the parameters are linearly dependent if there exist k constants α i≠ such that α1s¯k1+α2s¯k2+…αks¯kk= if the above relation does not hold, the set is independent. when the equality holds only approximately, the parameters are nearly dependent or nearly collinear. the degree of collinearity among a set of parameters can be measured by the collinearity index, cik, which is defined as: cik=1min∥α∥=1∥s¯kα∥=1λk,min. where s¯k is the sensitivity matrix built from the k sensitivity vectors, s¯k=, and λ k,min is the smallest eigenvalue of s¯kts¯k. the larger the collinearity index is, the more dependent the corresponding parameters are. brun and co-authors proposed to classify a subset of parameters as identifiable if their collinearity index is smaller than a threshold which they chose as cik< roughly speaking, a value of means that 95% of the variation in the model output caused by changing one of the parameters in the subset can be compensated by changing the other parameters in the set. other approaches for finding parameter correlations using sensitivity-based measures have been previously proposed in the literature. li and vu presented two methods that search for relationships among parameters in the context of a priori identifiability analysis. the method in provides a necessary but not sufficient condition for identifiability of nonlinear systems, which need to be fully observed. the method in removes the requirement of measuring all the system states, replacing it with the restriction that the model must be linear. we remark that the method proposed in the present manuscript does not have these limitations: it can be applied to partially observed, nonlinear systems with noisy, discrete-time measurements. largest identifiable subset as explained in the previous subsections, a subset of parameters is considered identifiable if its elements are influential and their sensitivity vectors are not collinear. we are interested in finding the largest set of parameters for which the collinearity of the corresponding sensitivity vectors is below the chosen threshold, cik< such a set of parameters represents all the degrees of freedom in the model. this means that perturbing a parameter not included in this set has an effect in the model predictions that can be compensated by changing other parameters in the set. however, a perturbation in a parameter belonging to the set cannot be compensated by changes in the remaining parameters. several methods have been developed for finding the group of identifiable parameters. iterative selection methods apply a step-wise procedure to select one parameter at a time, until no more parameters can be added to the identifiable set. in each step the parameter to be included is selected based on an optimality criteria. for example, the modified gram-schmidt orthogonalization method projects all the remaining sensitivity vectors to the subspace spanned by the already selected sensitivity vectors, and includes the parameter corresponding to the one with the largest projection value. this step is repeated until the largest projection value falls below a threshold, which means that the next parameter would significantly interplay with the parameters already included. the computational cost of this method scales up well with the number of sensitivity vectors. however, the drawback of iterative procedures such as this one is that the solution might not be the global optimum, that is, it might fail to find the largest identifiable subset. alternatively, we propose to solve the problem of finding the largest identifiable subset of all the estimated parameters using combinatorial optimization. to this end we formulate it as a integer optimization problem, where the goal is to maximize the number of sensitivity vectors included in the set, with the constraint that the corresponding collinearity index is below a threshold ci∗. this algorithm can be stated as maximizei∈{ }nθ∑k=1nθik subject tosi=cat{sk|ik= fork= …,nθ} ci<ci∗ ikis a binary variable fork= …,nθ where the binary variable i k indicates if the k-th parameter is included or not included in the identifiable group of parameters. the sensitivity matrix corresponding to the selected parameters is s i, and ‘cat’ stands for the concatenation of the column vectors in the constraint. the collinearity index of this matrix is ci and it is determined by computing the minimum eigenvalue as in. this combinatorial optimization problem has an exponentially scaling computational cost, and thus its solution requires an efficient algorithm. we chose the variable neighbourhood search technique, which is a heuristic global optimization method for integer optimization problems. we used the version of vns included in the meigo toolbox, which is implemented in matlab. we modified this initial formulation of the problem described in eqs. – after finding that its solution is often not unique: even after maximizing the number of parameters in the subset, there may be multiple subsets that yield a collinearity index below the threshold ci∗. indeed, we found large variability in the solutions if no initial guess was specified. therefore, we reformulated the optimization problem in two ways, as described in the following paragraphs. as a first modification, we transformed the collinearity requirement from a ‘hard’ to a ‘soft’ constraint. the modified optimization problem reads as maximizei∈{ }nθ∑k=1nθik−p1−p subject tosi=catsk|ik= fork= …,nθ p1=12ci/ci∗ p2=0ifci<ci∗αci−ci∗βotherwise ikis a binary variable fork= …,nθ as above, the binary variable i k indicates if the k-th parameter is included or not included in the selected group of parameters. the penalty p is a monotone increasing function of the collinearity index ci, such that p is when the collinearity equals to the threshold. due to this small value, p does not influence the size of the largest subset below the threshold. in this way, when multiple sets of the same size co-exist, the set with smaller collinearity index is always favoured. this results in an unique solution of the optimization problem if there are no sets with identical collinearity index. the second penalty function p represents a soft constraint that is active when the collinearity exceeds the threshold. the steepness of this constraint is tuned by the values of α and β, which we set to α= and β= our second improvement of the formulation of the optimization problem consists in providing a good initial guess of the solution using qr decomposition. the rank revealing qr decomposition algorithm, or rrqr, rewrites a matrix s as Πs=qr, where q is an orthogonal matrix, r is an upper triangular matrix, and Π is a permutation matrix. due to the properties of this decomposition, the permutation matrix defines a reordering of the columns of s. in this re-ordered matrix s ro=Π s, the most orthogonal columns are located in the left. in other words, the first n columns of the reordered matrix define a linear subspace, and the -th column has the largest projection value on this subspace among the remaining n θ−n columns located to the right of the n-th column. the outcome of the rrqr technique is similar to that of the aforementioned gram-schmidt orthogonalization method, but its implementation is more efficient. we applied rank revealing qr decomposition to the sensitivity matrix, following the procedure described in algorithm then, we used the resulted ordering of the sensitivity vectors to initialize the global optimizer. in this way we improved the performance of the global optimizer, which often found larger sets with collinearity index below the threshold value. the whole procedure for identifying the largest non-collinear subset of parameters is summarized in algorithm the procedure presented in this subsection has similarities with the one proposed by chu and hahn. one difference is that we maximize the subset size for a given collinearity threshold, whereas chu and hahn adopted the opposite approach, i.e., maximizing parametric identifiability for a pre-specified subset size. additionally, both methods differ in the optimization technique: we use variable neighbourhood search, which has better scalability than the genetic algorithm chosen in. recently, nienałtowski et al. have proposed a method for finding clusters of correlated parameters using so-called canonical correlation analysis. cca is an extension of pearson correlation for measuring multidimensional correlations between groups of parameters. given two groups of parameters of sizes m and n, with m<n, calculation of the canonical correlations provides m measures, which are summarized in a single measure, called mi-cca. this similarity measure represents the mutual information between the two groups, although it should be noted that average mutual information is equivalent to canonical correlation only if the random variables follow an elliptically symmetric probability model. nienałtowski et al. use mi-cca to cluster parameters until an identifiable subset is reached. this approach is sequential and yields a single parameter subset, which is possibly not maximal. in contrast, the methodology described here combines an initial sequential phase with a subsequent combinatorial optimization procedure. the second phase yields several identifiable parameter subsets and usually improves the initial solution. finding all largest subsets as mentioned above, the largest non-collinear subset of parameters is not unique. to realize this, imagine that we have a non-collinear set of the parameters, and consider an additional pair of highly collinear parameters. since we may add either of these two parameters to the set, but not both of them, we have two potential solutions. the optimization algorithm described above would choose the option with a lower collinearity index. however, we may also be interested in enumerating all the possible sets, instead of only one. finding all the largest subsets is a combinatorial problem too, which is computationally expensive. a naive approach for solving it could be to generate all possible sets of parameters and compute the corresponding collinearity index. however, note that if two parameters θ and θ are collinear, then any sets including the pair { θ θ 2} are highly collinear. using this fact, we developed an incremental procedure for the systematic determination of the sets. we start by considering all possible pairs of parameters and determining their collinearity. then we extend only those pairs which have a small collinearity index, by considering all possible combinations of a third parameter. this procedure is repeated until either all the sets are highly collinear, or there is only one set containing all the parameters. in this way, summarized in algorithm we can find all the largest subsets of non-collinear parameters. partitioning the non-identifiable parameters the two procedures presented above can be used for finding the largest, least collinear subset of parameters, and all the largest subsets; in both cases, restricted to those subsets whose collinearity falls below a threshold. however, it is often important to understand why certain parameters are not identifiable. for example, a parameter may be unidentifiable because the model output has very low sensitivity to changes in its value. but it could also be because it is highly correlated with another parameter, even when both parameters have high sensitivities. finding small groups of highly collinear parameters can be helpful in determining the exact source of unidentifiability. the collinearity of a subset always increases when a new parameter is added to the set. for example, considering three parameters, the collinearity of the triplet is always higher than the collinearity of any pairs. therefore, if a larger set of parameters contains a collinear pair, then the collinearity index of the large set is also large. if we are interested in finding the smallest groups of highly collinear parameters, we can proceed as follows. first we generate all possible pairs of parameters, and compute the collinearity of the corresponding sensitivity vectors. then we evaluate all possible triplets. the procedure can be extended for the analysis of larger sets. however, due to the combinatorial explosion of the computational cost, this method can be applied only to models of moderate size. visualization of identifiable subsets it can be useful to represent the identifiability results graphically, because such visualization can provide modellers with insight about how to reformulate their models and/or design new experiments in order to avoid non-identifiable parameters. with this aim, we display the model structure in the natural network visualization technique. an example is shown in fig. 1 c. the model structure is represented as a graph whose nodes are state variables, observables, stimuli, and model parameters. the edges – which can be directed or undirected – have the following meaning: an arrow from node a to node b indicates that node b appears in the equation of a. for example, if the dynamic equation of a state x is x˙1=p1·x the corresponding graph would show two arrows x 2→x and p 1→x fig. tgf- β model: panel a shows the model sensitivities with respect to the parameters in logarithmic scale. in panel b the maximum number of identifiable parameters is depicted depending on the collinearity threshold. the group-size does not change much by the threshold. the red vertical line indicates our choice for the further analysis. in panel c the mathematical model is depicted. nodes indicate states, identifiable and not identifiable parameters, and observables. in panel d the interplay among the collinear parameters are indicated. there are groups of parameters: in the first groups triplets of parameters show large collinearity, while in the fifth group parameters are collinear more formally, we determine how the state, input variables, stimuli and parameters are connected and influence each other through symbolic manipulation of the model eq.. for this purpose we compute: the jacobian matrix with respect to the states: ji,jss=∂fi∂xj, the jacobian of the observation functions with respect to the states: ji,jso=∂gi∂xj, the jacobian of the systems dynamics with respect to the stimuli ji,jsi=∂fi∂uj, and the jacobian with respect to the parameters ji,jsp=∂fi∂θj. all these matrices are evaluated symbolically, and then the expressions are converted to a logical or when the symbolic result is zero. additionally, we can connect parameters by undirected edges if their collinearity is larger than the collinearity threshold. implementation: the visid software tool we implemented the techniques proposed in subsections “practical identifiability analysis” and “visualization of identifiable subsets” as a matlab software package called visid, which is provided as additional file and can also be downloaded from github. it is free software, made available under the terms of the gnu general public license version the visid toolbox relies on three other matlab toolboxes, which are also freely available: amigo, which is used to to store, simulate and calibrate the models; meigo, which implements the variable neighbouring search optimization method; and rrqr, which performs the rank revealing qr decomposition used to initialize the global optimizer. network visualization is performed with cytoscape. further details can be found in section of additional file  in this section we demonstrate the application of the methodology presented in the previous section using several dynamic systems biology models of different type and complexity. their main characteristics are given in table first we present detailed results of identifiability analysis and visualization for a model of the tgf- β signalling pathway. we also provide similar results for the genetic network that controls the circadian clock in arabidopsis thaliana. due to their complexity and yet relatively moderate size, these models are well suited as case studies for illustrating the identifiability methodology in depth. a. thaliana then we study two large scale benchmark problems included in the biopredyn-bench collection. since the analysis of these latter models is more challenging due to their larger size, we start by demonstrating the performance improvements that can be achieved during parameter estimation using the model calibration procedure proposed in section “parameter estimation with regularization and global optimization”. then we perform identifiability analysis and report the corresponding results, including the graphical representation of the identifiable subset using the natural network visualization. tgf- β signalling pathway the dynamic model of the tgf- β signaling pathway was presented in as a tutorial example for model calibration. it has dynamic states and kinetic parameters, of which need to be estimated. following, we assumed that all the concentrations, except the smad rnas, can be measured in the experiments. the algebraic equations of the reaction kinetics and the dynamic equations are provided in the additional file for the purpose of testing the methodology we generated a training dataset by simulating the model equations using the nominal values of the parameters k 1– k. then we sampled the simulated trajectories at equidistant time points, and added normally distributed random numbers to the data to mimic measurement errors. finally, we estimated the model parameters from the generated data set. this approach is widely used for testing calibration methods and assessing the extent to which they recover the nominal parameters. it should be noted that, as the amount of noise in the dataset increases, the information/signal ratio decreases, making the estimation problem more ill-conditioned. this makes it more difficult to recover the correct value of the parameters, but has a small effect in computation times. the numerical values of the estimated parameters are reported in additional file 2: table s we started the identifiability analysis by computing the sensitivities of the observations with respect to the estimated model parameters, according to algorithm we found that all the parameters have a non-negligible influence on the model outputs, thus there are no individually non-identifiable parameters. next, following algorithm we applied qr decomposition and ranked the parameters according to their orthogonality. we then solved the optimization problem – by initializing the variable neighboring search method with the results of the qr decomposition. setting the threshold level for the collinearity index to ci = yielded identifiable parameters, which are shown as green nodes in the network in fig. 1 c. parameters not present in the identifiable subset are shown as red nodes. parameters are connected by arrows to state variables if they appear in the equation of the corresponding dynamic equation. states which directly influence each other are also connected by directed edges in the same manner. blue squares represent measurements; a state is connected to a blue square if it appears in the corresponding observation function. to see how the size of the identifiable subset is influenced by the choice of the collinearity index threshold, we solved the optimization problem for a range of threshold values. the results are depicted in fig. 1 b. as the collinearity index threshold decreases, less parameters are considered identifiable. we can see that the identifiability results are quite robust to the choice of threshold level: the number of identifiable parameters is always between and and it is constant for a very wide range of ci, 15≤ci≤ the results presented so far tell us that the parameters are not correlated. however, they do not inform of the relationships among identifiable and non-identifiable parameters. to investigate this point, we computed the smallest correlated subsets as described in section “partitioning the non-identifiable parameters”, up to groups of parameters. figure 1 d shows such groups; parameters are depicted as blue circles connected with group identifying nodes. these nodes are labeled as gx, where x indicates the number of parameters in the group and y is the group index for a given number of parameters stands for the second group of three correlated parameters). we found that the large pairwise collinearity between k 14−k and k 17−k explains the non-identifiability of the model parameters only partially. there are groups of triplets and a group of parameters which are highly correlated. the members of the groups and the corresponding collinearity index are reported in table d it is important to note that collinearity might arise among multiple parameters, even if they are pairwise independent. for example, despite the fact that none of the pairs in the group of k k and k has a high pairwise collinearity, the collinearity index of the triplet is extremely large. algorithm found different sets of identifiable parameters with collinearity index ranging between and less than the threshold. the sets are reported with the corresponding collinearity index in additional file 2: table s we can see that parameters {k 1−k k k 15} are members of all the groups, and they do not participate in any of the small correlated groups in fig. 1 d. from each correlated group of size k, only k− parameters can participate in the largest set of identifiable parameters. the aforementioned identifiability procedures can be carried out in a few seconds. detailed computational costs are shown in table s of the additional file for all the case studies considered in this paper. circadian clock in arabidopsis thaliana locke and co-authors described the genetic network controlling the circadian clock in arabidopsis thaliana; the dynamic equations of this model are provided in the additional file we generated training data by simulating the model equations with the nominal parameters in two experimental conditions. in the first one, the model input was kept constant, representing continuous light stimulation of the plant. in the second experiment the input was changed pulse-wise in 12-hour cycles, repeated times. as in the previous example, the trajectories were sampled at equidistant time-points and disturbed by pseudo-random noise. only two states, c t m and c l m, were observed. the estimated model parameters are collected in additional file 2: table s although the model outputs showed sensitivity to all the parameters, i.e. there were no zero sensitivity vectors, we found that most of the model parameters are non-identifiable due to heavy collinearities. the largest identifiable subset contains only of the parameters, depicted in fig. 2 d by green nodes. the enumeration of the largest sets of identifiable parameters by algorithm identified parameter sets. fig. circadian clock in arabidopsis thaliana: panel a shows the model sensitivities with respect to parameters in logarithmic scale. in panel b the maximum number of identifiable parameters is depicted as a function of the collinearity threshold. the red vertical line indicates our choice of threshold used for further analysis. in panel c the interplay among the collinear parameters is indicated up to groups of parameters. in panel d the mathematical model is depicted. nodes indicate states, identifiable and not identifiable parameters, and observables. highly correlated parameter pairs are connected by undirected edges benchmarks b and b from the biopredyn-bench collection in this subsection we analyze two large scale benchmark problems taken from the biopredyn-bench collection: the metabolic models of escherichia coli and chinese hamster ovary cells. they are highly non-linear, partially observed systems with more than unknown parameters, which pose serious challenges for parameter identification. in b the calibration data was generated by model simulation and disturbed by random noise, while in b it was experimentally measured. further details about the models and the parameter estimation challenge can be found in. first we use these benchmarks to illustrate the benefits of the parameter estimation strategy proposed in subsection “parameter estimation with regularization and global optimization”, comparing it with the one used in. both approaches use a hybrid method, ess, which combines a global optimization algorithm with a local search. in the local method of choice was fmincon; here we compare that configuration with nl2sol. global optimization algorithms use pseudo-random numbers. hence their performance changes at every run, and the calibration problem should be solved several times to obtain more robust results. since each optimization takes several hours we limited the number of runs to five for each problem. we used the approximate computation time reported in as the stopping criterion for the model calibration. convergence curves depict the best objective function value found versus cpu time, and can be used to compare the performance of different algorithms. an optimization method is preferred if it achieves a lower objective function value at earlier cpu time. the best convergence curves corresponding to b and b are shown in the additional file for algorithms: ess-fmincon, as reported in; ess-nl2sol; and ess-nl2sol using regularization as recommended in section “parameter estimation with regularization and global optimization”. from those curves we see that the algorithm proposed here converged earlier than the others to the optimal objective function value. we stress that the main purpose of regularization is to avoid overfitting: we do not wish to obtain an excessively good fit, which would indicate that we are reproducing noise instead of the true dynamics. therefore, regularization should not achieve a smaller objective function value. next, we apply the identifiability analysis procedures presented in subsections “practical identifiability analysis” and “visualization of identifiable subsets” to these two models. for b they yield an identifiable subset of size and for b of size. the corresponding networks are shown in figs. and respectively. it is also possible to find small groups of highly correlated parameters for models of this size; e.g. for b we obtained those depicted in fig. fig. representation of the connections in the b model using the network diagram formalism. nodes indicate states, identifiable and not identifiable parameters, observables, and inputs. the source file of this figure is provided with the visid toolbox; using cytoscape, the user can navigate through it and zoom on different areas to improve the visibility fig. representation of the connections in the b model using the network diagram formalism. nodes indicate states, identifiable and not identifiable parameters, observables, and inputs fig. visualization of the relationships among highly collinear parameters in the b model. the figure shows small groups, whose sizes range between and parameters. unidentifiable parameters are shown in red; identifiable parameters in green. highly correlated pairs are connected by lines the aforementioned results show that both models are poorly identifiable in practice for the considered datasets; more informative data would be needed in order to obtain accurate estimates of their parameters. discussion and  in this paper we have presented a workflow to efficiently estimate the parameters of dynamic models and analyze their practical identifiability. our approach combines an advanced optimization technique, which reduces computation times in parameter estimation, and several identifiability analysis procedures, which can find subsets of identifiable and unidentifiable parameters. results are visualized using network diagrams, which provide an intuitive representation of the findings and facilitate their analysis and understanding. many approaches have been applied to study identifiability of kinetic models, but they suffer from lack of scalability. an advantage of the integrated method presented here is its moderate computational cost, which enables its application to large-scale models; complete results can be obtained in a few hours for models of more than a hundred parameters. another important aspect is the integration of identifiability analysis with visualization, which presents the results in a way that is easily interpretable for modelers and experimentalists. currently, its main limitation arises when trying to find all the different existing groups of highly correlated parameters: the combinatorial explosion of this particular task makes it feasible only for models of moderate size, i.e. of a few dozens of parameters. however, all the remaining steps of the workflow presented in this manuscript scale up well up to several hundred parameters. the usefulness of the methodology and workflow presented here goes beyond basic parameter identifiability analysis. the procedure not only determines the largest subset of identifiable parameters, but also informs about the characteristics of the space of non-identifiable parameters, reporting small groups of highly correlated parameters, and presents all these results in a coherent and scalable way using visualization techniques, facilitating the understanding of the underlying complex interactions. uncovering these higher order relationships helps in determining the causes of unidentifiability and provides guidelines for remedying them, e.g. by reformulating the model or by collecting new data through a new experimental design. all this information can be readily used to improve the iterative model-building cycle. a matlab implementation of the identifiability and visualization methodology, which we have called the visid software package, is available from github as free, open source software. this distribution includes the case studies discussed above. additional files additional file visid toolbox. this compressed folder contains the visid matlab toolbox. additional file supplementary material. this document contains detailed descriptions of the case studies and of the visid toolbox, as well as additional details about the results. abbreviations ccacanonical correlation analysis cicollinearity index cpucentral processing unit essenhanced scatter search mimutual information odeordinary differential equation qrdecomposition of a matrix into an orthogonal matrix and an upper triangular matrix rnaribonucleic acid rrqrrank revealing qr decomposition tgf- βtransforming growth factor beta vnsvariable neighbourhood search 
we also use the collinearity index to identify small groups of highly correlated parameters. therefore, before attempting to perform parameter estimation it is important to characterize the subset of identifiable parameters and their interplay. once this is achieved, it is still necessary to perform parameter estimation, which poses additional challenges. our contributions alleviate the difficulties that appear at different stages of the identifiability analysis and parameter estimation process. parameter estimationdynamic modelsidentifiabilityglobal optimizationregularizationoverfittinghttp://dx.doi.org/ /501100007601horizon 2020686282http://dx.doi.org/ /501100004963seventh framework programme289384minecodpi2014-55276-c5-2-rissue-copyright-statement© the author 2017 the results files can be visualized using cytoscape, showing the identifiable and non-identifiable groups of parameters together with the model structure in the same graph. we present a methodology that detects high-order relationships among parameters, and visualizes the results to facilitate further analysis. we use a collinearity index to quantify the correlation between parameters in a group in a computationally efficient way. it enables the practical identifiability analysis of dynamic models of large size, and accelerates their calibration. the identifiability analysis techniques are implemented as a matlab toolbox called visid, which is freely available as open source from github. we show how to combine global optimization and regularization techniques for calibrating medium and large scale biological models with moderate computation times. the visualization tool allows modellers to detect parts that are problematic and need refinement or reformulation, and provides experimentalists with information that can be helpful in the design of new experiments. kinetic models of biochemical systems usually consist of ordinary differential equations that have many unknown parameters. then we evaluate the practical identifiability of the estimated parameters using the proposed methodology. some of these parameters are often practically unidentifiable, that is, their values cannot be uniquely determined from the available data.
 kinetic models of biochemical systems usually consist of ordinary differential equations that have many unknown parameters. some of these parameters are often practically unidentifiable, that is, their values cannot be uniquely determined from the available data. possible causes are lack of influence on the measured outputs, interdependence among parameters, and poor data quality. uncorrelated parameters can be seen as the key tuning knobs of a predictive model. therefore, before attempting to perform parameter estimation it is important to characterize the subset of identifiable parameters and their interplay. once this is achieved, it is still necessary to perform parameter estimation, which poses additional challenges. methods we present a methodology that detects high-order relationships among parameters, and visualizes the results to facilitate further analysis. we use a collinearity index to quantify the correlation between parameters in a group in a computationally efficient way. then we apply integer optimization to find the largest groups of uncorrelated parameters. we also use the collinearity index to identify small groups of highly correlated parameters. the results files can be visualized using cytoscape, showing the identifiable and non-identifiable groups of parameters together with the model structure in the same graph.  our contributions alleviate the difficulties that appear at different stages of the identifiability analysis and parameter estimation process. we show how to combine global optimization and regularization techniques for calibrating medium and large scale biological models with moderate computation times. then we evaluate the practical identifiability of the estimated parameters using the proposed methodology. the identifiability analysis techniques are implemented as a matlab toolbox called visid, which is freely available as open source from github.  our approach is geared towards scalability. it enables the practical identifiability analysis of dynamic models of large size, and accelerates their calibration. the visualization tool allows modellers to detect parts that are problematic and need refinement or reformulation, and provides experimentalists with information that can be helpful in the design of new experiments. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords parameter estimationdynamic modelsidentifiabilityglobal optimizationregularizationoverfittinghttp://dx.doi.org/ /501100007601horizon 2020686282http://dx.doi.org/ /501100004963seventh framework programme289384minecodpi2014-55276-c5-2-rissue-copyright-statement© the author 2017 