An algorithm for finding globally identifiable parameter combinations of nonlinear ODE models using Gröbner Bases

The parameter identifiability problem for dynamic system ODE models has been extensively studied. Nevertheless, except for linear ODE models, the question of establishing identifiable combinations of parameters when the model is unidentifiable has not received as much attention and the problem is not fully resolved for nonlinear ODEs. Identifiable combinations are useful, for example, for the reparameterization of an unidentifiable ODE model into an identifiable one. We extend an existing algorithm for finding globally identifiable parameters of nonlinear ODE models to generate the ‘simplest’ globally identifiable parameter combinations using Gröbner Bases. We also provide sufficient conditions for the method to work, demonstrate our algorithm and find associated identifiable reparameterizations for several linear and nonlinear unidentifiable biomodels.     

Parameter Identifiability and Redundancy in a General Class of Stochastic Carcinogenesis Models

Background

Heidenreich et al. (Risk Anal 1997 17 391–399) considered parameter identifiability in the context of the two-mutation cancer model and demonstrated that combinations of all but two of the model parameters are identifiable. We consider the problem of identifiability in the recently developed carcinogenesis models of Little and Wright (Math Biosci 2003 183 111–134) and Little et al. (J Theoret Biol 2008 254 229–238). These models, which incorporate genomic instability, generalize a large number of other quasi-biological cancer models, in particular those of Armitage and Doll (Br J Cancer 1954 8 1–12), the two-mutation model (Moolgavkar et al. Math Biosci 1979 47 55–77), the generalized multistage model of Little (Biometrics 1995 51 1278–1291), and a recently developed cancer model of Nowak et al. (PNAS 2002 99 16226–16231).

Methodology/Principal Findings

We show that in the simpler model proposed by Little and Wright (Math Biosci 2003 183 111–134) the number of identifiable combinations of parameters is at most two less than the number of biological parameters, thereby generalizing previous results of Heidenreich et al. (Risk Anal 1997 17 391–399) for the two-mutation model. For the more general model of Little et al. (J Theoret Biol 2008 254 229–238) the number of identifiable combinations of parameters is at most less than the number of biological parameters, where is the number of destabilization types, thereby also generalizing all these results. Numerical evaluations suggest that these bounds are sharp. We also identify particular combinations of identifiable parameters.

Conclusions/Significance

We have shown that the previous results on parameter identifiability can be generalized to much larger classes of quasi-biological carcinogenesis model, and also identify particular combinations of identifiable parameters. These results are of theoretical interest, but also of practical significance to anyone attempting to estimate parameters for this large class of cancer models.     

Determination of parameter identifiability in nonlinear biophysical models: A Bayesian approach

A major goal of biophysics is to understand the physical mechanisms of biological molecules and systems. Mechanistic models are evaluated based on their ability to explain carefully controlled experiments. By fitting models to data, biophysical parameters that cannot be measured directly can be estimated from experimentation. However, it might be the case that many different combinations of model parameters can explain the observations equally well. In these cases, the model parameters are not identifiable: the experimentation has not provided sufficient constraining power to enable unique estimation of their true values. We demonstrate that this pitfall is present even in simple biophysical models. We investigate the underlying causes of parameter non-identifiability and discuss straightforward methods for determining when parameters of simple models can be inferred accurately. However, for models of even modest complexity, more general tools are required to diagnose parameter non-identifiability. We present a method based in Bayesian inference that can be used to establish the reliability of parameter estimates, as well as yield accurate quantification of parameter confidence.     

Parameter Identifiability and Redundancy: Theoretical Considerations

Background

Models for complex biological systems may involve a large number of parameters. It may well be that some of these parameters cannot be derived from observed data via regression techniques. Such parameters are said to be unidentifiable, the remaining parameters being identifiable. Closely related to this idea is that of redundancy, that a set of parameters can be expressed in terms of some smaller set. Before data is analysed it is critical to determine which model parameters are identifiable or redundant to avoid ill-defined and poorly convergent regression.

Methodology/Principal Findings

In this paper we outline general considerations on parameter identifiability, and introduce the notion of weak local identifiability and gradient weak local identifiability. These are based on local properties of the likelihood, in particular the rank of the Hessian matrix. We relate these to the notions of parameter identifiability and redundancy previously introduced by Rothenberg (Econometrica 39 (1971) 577–591) and Catchpole and Morgan (Biometrika 84 (1997) 187–196). Within the widely used exponential family, parameter irredundancy, local identifiability, gradient weak local identifiability and weak local identifiability are shown to be largely equivalent. We consider applications to a recently developed class of cancer models of Little and Wright (Math Biosciences 183 (2003) 111–134) and Little et al. (J Theoret Biol 254 (2008) 229–238) that generalize a large number of other recently used quasi-biological cancer models.

Conclusions/Significance

We have shown that the previously developed concepts of parameter local identifiability and redundancy are closely related to the apparently weaker properties of weak local identifiability and gradient weak local identifiability—within the widely used exponential family these concepts largely coincide.     

Parameter identification in N-compartment mamillary models

We consider a general mamillary model with a central compartment (compartment 1) and n?1 peripheral compartments, each bidirectionally connected to the first. Elimination is allowed from any compartment and effectively occurs from the system. With input introduced into an arbitrary compartment and measurement performed in an arbitrary compartment, explicit equations are given to derive the parameters of the model from the input-output procedure. The calculations include essentially the determination of the roots of a polynomial plus some elementary algebra. If input and measurement are performed in the same compartment, then a set of 2n elementary combinations of the model parameters can be uniquely determined. However, the model parameters themselves can only be localized, each within an interval. These intervals are explicitly calculated and their width discussed.     

On Finding and Using Identifiable Parameter Combinations in Nonlinear Dynamic Systems Biology Models and COMBOS: A Novel Web Implementation

Parameter identifiability problems can plague biomodelers when they reach the quantification stage of development, even for relatively simple models. Structural identifiability (SI) is the primary question, usually understood as knowing which of P unknown biomodel parameters p ₁,…, p_i,…, p_P are-and which are not-quantifiable in principle from particular input-output (I-O) biodata. It is not widely appreciated that the same database also can provide quantitative information about the structurally unidentifiable (not quantifiable) subset, in the form of explicit algebraic relationships among unidentifiable p_i. Importantly, this is a first step toward finding what else is needed to quantify particular unidentifiable parameters of interest from new I–O experiments. We further develop, implement and exemplify novel algorithms that address and solve the SI problem for a practical class of ordinary differential equation (ODE) systems biology models, as a user-friendly and universally-accessible web application (app)–COMBOS. Users provide the structural ODE and output measurement models in one of two standard forms to a remote server via their web browser. COMBOS provides a list of uniquely and non-uniquely SI model parameters, and–importantly-the combinations of parameters not individually SI. If non-uniquely SI, it also provides the maximum number of different solutions, with important practical implications. The behind-the-scenes symbolic differential algebra algorithms are based on computing Gröbner bases of model attributes established after some algebraic transformations, using the computer-algebra system Maxima. COMBOS was developed for facile instructional and research use as well as modeling. We use it in the classroom to illustrate SI analysis; and have simplified complex models of tumor suppressor p53 and hormone regulation, based on explicit computation of parameter combinations. It’s illustrated and validated here for models of moderate complexity, with and without initial conditions. Built-in examples include unidentifiable 2 to 4-compartment and HIV dynamics models.     

Construction of a large scale integrated map of macrophage pathogen recognition and effector systems

Background

Mathematical models provide abstract representations of the information gained from experimental observations on the structure and function of a particular biological system. Conferring a predictive character on a given mathematical formulation often relies on determining a number of non-measurable parameters that largely condition the model's response. These parameters can be identified by fitting the model to experimental data. However, this fit can only be accomplished when identifiability can be guaranteed.

Results

We propose a novel iterative identification procedure for detecting and dealing with the lack of identifiability. The procedure involves the following steps: 1) performing a structural identifiability analysis to detect identifiable parameters; 2) globally ranking the parameters to assist in the selection of the most relevant parameters; 3) calibrating the model using global optimization methods; 4) conducting a practical identifiability analysis consisting of two (a priori and a posteriori) phases aimed at evaluating the quality of given experimental designs and of the parameter estimates, respectively and 5) optimal experimental design so as to compute the scheme of experiments that maximizes the quality and quantity of information for fitting the model.

Conclusions

The presented procedure was used to iteratively identify a mathematical model that describes the NF-κB regulatory module involving several unknown parameters. We demonstrated the lack of identifiability of the model under typical experimental conditions and computed optimal dynamic experiments that largely improved identifiability properties.     

Assessment of uncertainty in functional-structural plant models

Background and Aims

Constructing functional–structural plant models (FSPMs) is a valuable method for examining how physiology and morphology interact in determining plant processes. However, such models always have uncertainty concerned with whether model components have been selected and represented effectively, with the number of model outputs simulated and with the quality of data used in assessment. We provide a procedure for defining uncertainty of an FSPM and how this uncertainty can be reduced.

Methods

An important characteristic of FSPMs is that typically they calculate many variables. These can be variables that the model is designed to predict and also variables that give indications of how the model functions. Together these variables are used as criteria in a method of multi-criteria assessment. Expected ranges are defined and an evolutionary computation algorithm searches for model parameters that achieve criteria within these ranges. Typically, different combinations of model parameter values provide solutions achieving different combinations of variables within their specified ranges. We show how these solutions define a Pareto Frontier that can inform about the functioning of the model.

Key Results

The method of multi-criteria assessment is applied to development of BRANCHPRO, an FSPM for foliage reiteration on old-growth branches of Pseudotsuga menziesii. A geometric model utilizing probabilities for bud growth is developed into a causal explanation for the pattern of reiteration found on these branches and how this pattern may contribute to the longevity of this species.

Conclusions

FSPMs should be assessed by their ability to simulate multiple criteria simultaneously. When different combinations of parameter values achieve different groups of assessment criteria effectively a Pareto Frontier can be calculated and used to define the sources of model uncertainty.     

Identifiability of parameters in MCMC Bayesian inference of phylogeny

Methods for Bayesian inference of phylogeny using DNA sequences based on Markov chain Monte Carlo (MCMC) techniques allow the incorporation of arbitrarily complex models of the DNA substitution process, and other aspects of evolution. This has increased the realism of models, potentially improving the accuracy of the methods, and is largely responsible for their recent popularity. Another consequence of the increased complexity of models in Bayesian phylogenetics is that these models have, in several cases, become overparameterized. In such cases, some parameters of the model are not identifiable; different combinations of nonidentifiable parameters lead to the same likelihood, making it impossible to decide among the potential parameter values based on the data. Overparameterized models can also slow the rate of convergence of MCMC algorithms due to large negative correlations among parameters in the posterior probability distribution. Functions of parameters can sometimes be found, in overparameterized models, that are identifiable, and inferences based on these functions are legitimate. Examples are presented of overparameterized models that have been proposed in the context of several Bayesian methods for inferring the relative ages of nodes in a phylogeny when the substitution rate evolves over time.     

Estimating trees from filtered data: Identifiability of models for morphological phylogenetics

As an alternative to parsimony analyses, stochastic models have been proposed ( [Lewis, 2001] and [Nylander et al., 2004]) for morphological characters, so that maximum likelihood or Bayesian analyses may be used for phylogenetic inference. A key feature of these models is that they account for ascertainment bias, in that only varying, or parsimony-informative characters are observed. However, statistical consistency of such model-based inference requires that the model parameters be identifiable from the joint distribution they entail, and this issue has not been addressed.Here we prove that parameters for several such models, with finite state spaces of arbitrary size, are identifiable, provided the tree has at least eight leaves. If the tree topology is already known, then seven leaves suffice for identifiability of the numerical parameters. The method of proof involves first inferring a full distribution of both parsimony-informative and non-informative pattern joint probabilities from the parsimony-informative ones, using phylogenetic invariants. The failure of identifiability of the tree parameter for four-taxon trees is also investigated.     

A global distribution of biodiversity inferred from climatic constraints: results from a process‐based modelling study

We investigated the connection between plant species diversity and climate by using a process‐based, generic plant model. Different ‘species' were simulated by different values for certain growth‐related model parameters. Subsequently, a wide range of values were tested in the framework of a ‘Monte Carlo' simulation for success; that is, the capability of each plant with these parameter combinations to reproduce itself during its lifetime. The range of successful parameter combinations approximated species diversity. This method was applied to a global grid, using daily atmospheric forcing from a climate model simulation. The computed distribution of plant ‘species' diversity compares very well with the observed, global‐scale distribution of species diversity, reproducing the majority of ‘hot spot' areas of biodiversity. A sensitivity analysis revealed that the predicted pattern is very robust against changes of fixed model parameters. Analysis of the climatic forcing and of two additional sensitivity simulations demonstrated that the crucial factor leading to this distribution of diversity is the early stage of a plant's life when water availability is highly coupled to the variability in precipitation because in this stage root‐zone storage of water is small. We used cluster analysis in order to extract common sets of species parameters, mean plant properties and biogeographic regions (biomes) from the model output. The successful ‘species' cannot be grouped into typical parameter combinations, which define the plant's functioning. However, the mean simulated plant properties, such as lifetime and growth, can be grouped into a few characteristic plant ‘prototypes', ranging from short‐lived, fast growing plants, similar to grasses, to long‐lived, slow growing plants, similar to trees. The classification of regions with respect to similar combinations of successful ‘species' yields a distribution of biomes similar to the observed distribution. Each biome has typical levels of climatic constraints, expressed for instance by the number of ‘rainy days' and ‘warm days'. The less the number of days favourable for growth, the greater the level of constraints and the less the ‘species' diversity. These results suggest that climate as a fundamental constraint can explain much of the global scale, observed distribution of plant species diversity.     

Minimal output sets for identifiability

Ordinary differential equation models in biology often contain a large number of parameters that must be determined from measurements by parameter estimation. For a parameter estimation procedure to be successful, there must be a unique set of parameters that can have produced the measured data. This is not the case if a model is not uniquely structurally identifiable with the given set of outputs selected as measurements. In designing an experiment for the purpose of parameter estimation, given a set of feasible but resource-consuming measurements, it is useful to know which ones must be included in order to obtain an identifiable system, or whether the system is unidentifiable from the feasible measurement set. We have developed an algorithm that, from a user-provided set of variables and parameters or functions of them assumed to be measurable or known, determines all subsets that when used as outputs give a locally structurally identifiable system and are such that any output set for which the system is structurally identifiable must contain at least one of the calculated subsets. The algorithm has been implemented in Mathematica and shown to be feasible and efficient. We have successfully applied it in the analysis of large signalling pathway models from the literature.     

Model‐based identifiable parameter determination applied to a simultaneous saccharification and fermentation process model for bio‐ethanol production

In this work, a methodology for the model‐based identifiable parameter determination (MBIPD) is presented. This systematic approach is proposed to be used for structure and parameter identification of nonlinear models of biological reaction networks. Usually, this kind of problems are over‐parameterized with large correlations between parameters. Hence, the related inverse problems for parameter determination and analysis are mathematically ill‐posed and numerically difficult to solve. The proposed MBIPD methodology comprises several tasks: (i) model selection, (ii) tracking of an adequate initial guess, and (iii) an iterative parameter estimation step which includes an identifiable parameter subset selection (SsS) algorithm and accuracy analysis of the estimated parameters. The SsS algorithm is based on the analysis of the sensitivity matrix by rank revealing factorization methods. Using this, a reduction of the parameter search space to a reasonable subset, which can be reliably and efficiently estimated from available measurements, is achieved. The simultaneous saccharification and fermentation (SSF) process for bio‐ethanol production from cellulosic material is used as case study for testing the methodology. The successful application of MBIPD to the SSF process demonstrates a relatively large reduction in the identified parameter space. It is shown by a cross‐validation that using the identified parameters (even though the reduction of the search space), the model is still able to predict the experimental data properly. Moreover, it is shown that the model is easily and efficiently adapted to new process conditions by solving reduced and well conditioned problems. © 2013 American Institute of Chemical Engineers Biotechnol. Prog., 29:1064–1082, 2013     

Probabilistic Prediction of Contacts in Protein-Ligand Complexes

We introduce a statistical method for evaluating atomic level 3D interaction patterns of protein-ligand contacts. Such patterns can be used for fast separation of likely ligand and ligand binding site combinations out of all those that are geometrically possible. The practical purpose of this probabilistic method is for molecular docking and scoring, as an essential part of a scoring function. Probabilities of interaction patterns are calculated conditional on structural x-ray data and predefined chemical classification of molecular fragment types. Spatial coordinates of atoms are modeled using a Bayesian statistical framework with parametric 3D probability densities. The parameters are given distributions a priori, which provides the possibility to update the densities of model parameters with new structural data and use the parameter estimates to create a contact hierarchy. The contact preferences can be defined for any spatial area around a specified type of fragment. We compared calculated contact point hierarchies with the number of contact atoms found near the contact point in a reference set of x-ray data, and found that these were in general in a close agreement. Additionally, using substrate binding site in cathechol-O-methyltransferase and 27 small potential binder molecules, it was demonstrated that these probabilities together with auxiliary parameters separate well ligands from decoys (true positive rate 0.75, false positive rate 0). A particularly useful feature of the proposed Bayesian framework is that it also characterizes predictive uncertainty in terms of probabilities, which have an intuitive interpretation from the applied perspective.     

Identifying evolutionary trees and substitution parameters for the general Markov model with invariable sites

The general Markov plus invariable sites (GM+I) model of biological sequence evolution is a two-class model in which an unknown proportion of sites are not allowed to change, while the remainder undergo substitutions according to a Markov process on a tree. For statistical use it is important to know if the model is identifiable; can both the tree topology and the numerical parameters be determined from a joint distribution describing sequences only at the leaves of the tree? We establish that for generic parameters both the tree and all numerical parameter values can be recovered, up to clearly understood issues of 'label swapping'. The method of analysis is algebraic, using phylogenetic invariants to study the variety defined by the model. Simple rational formulas, expressed in terms of determinantal ratios, are found for recovering numerical parameters describing the invariable sites.     

A database of computational models of a half-center oscillator for analyzing how neuronal parameters influence network activity

A half-center oscillator (HCO) is a common circuit building block of central pattern generator networks that produce rhythmic motor patterns in animals. Here we constructed an efficient relational database table with the resulting characteristics of the Hill et al.’s (J Comput Neurosci 10:281–302, 2001) HCO simple conductance-based model. The model consists of two reciprocally inhibitory neurons and replicates the electrical activity of the oscillator interneurons of the leech heartbeat central pattern generator under a variety of experimental conditions. Our long-range goal is to understand how this basic circuit building block produces functional activity under a variety of parameter regimes and how different parameter regimes influence stability and modulatability. By using the latest developments in computer technology, we simulated and stored large amounts of data (on the order of terabytes). We systematically explored the parameter space of the HCO and corresponding isolated neuron models using a brute-force approach. We varied a set of selected parameters (maximal conductance of intrinsic and synaptic currents) in all combinations, resulting in about 10 million simulations. We classified these HCO and isolated neuron model simulations by their activity characteristics into identifiable groups and quantified their prevalence. By querying the database, we compared the activity characteristics of the identified groups of our simulated HCO models with those of our simulated isolated neuron models and found that regularly bursting neurons compose only a small minority of functional HCO models; the vast majority was composed of spiking neurons.     

A Global Sensitivity Test for Evaluating Statistical Hypotheses with Nonidentifiable Models

Summary We consider the problem of evaluating a statistical hypothesis when some model characteristics are nonidentifiable from observed data. Such a scenario is common in meta‐analysis for assessing publication bias and in longitudinal studies for evaluating a covariate effect when dropouts are likely to be nonignorable. One possible approach to this problem is to fix a minimal set of sensitivity parameters conditional upon which hypothesized parameters are identifiable. Here, we extend this idea and show how to evaluate the hypothesis of interest using an infimum statistic over the whole support of the sensitivity parameter. We characterize the limiting distribution of the statistic as a process in the sensitivity parameter, which involves a careful theoretical analysis of its behavior under model misspecification. In practice, we suggest a nonparametric bootstrap procedure to implement this infimum test as well as to construct confidence bands for simultaneous pointwise tests across all values of the sensitivity parameter, adjusting for multiple testing. The methodology's practical utility is illustrated in an analysis of a longitudinal psychiatric study.     

Identifiability and interval identifiability of mammillary and catenary compartmental models with some known rate constants

The identifiability problem is addressed for n-compartment linear mammillary and catenary models, for the common case of input and output in the first compartment and prior information about one or more model rate constants. We first define the concept of independent constraints and show that n-compartment linear mammillary or catenary models are uniquely identifiable under n-1 independent constraints. Closed-form algorithms for bounding the constrained parameter space are then developed algebraically, and their validity is confirmed using an independent approach, namely joint estimation of the parameters of all uniquely identifiable submodels of the original multicompartmental model. For the noise-free (deterministic) case, the major effects of additional parameter knowledge are to narrow the bounds of rate constants that remain unidentifiable, as well as to possibly render others identifiable. When noisy data are considered, the means of the bounds of rate constants that remain unidentifiable are also narrowed, but the variances of some of these bound estimates increase. This unexpected result was verified by Monte Carlo simulation of several different models, using both normally and lognormally distributed data assumptions. Extensions and some consequences of this analysis useful for model discrimination and experiment design applications are also noted.     

Structural identifiability of polynomial and rational systems

Since analysis and simulation of biological phenomena require the availability of their fully specified models, one needs to be able to estimate unknown parameter values of the models. In this paper we deal with identifiability of parametrizations which is the property of one-to-one correspondence of parameter values and the corresponding outputs of the models. Verification of identifiability of a parametrization precedes estimation of numerical values of parameters, and thus determination of a fully specified model of a considered phenomenon. We derive necessary and sufficient conditions for the parametrizations of polynomial and rational systems to be structurally or globally identifiable. The results are applied to investigate the identifiability properties of the system modeling a chain of two enzyme-catalyzed irreversible reactions. The other examples deal with the phenomena modeled by using Michaelis–Menten kinetics and the model of a peptide chain elongation.     

Modeling association among demographic parameters in analysis of open population capture-recapture data

We present a hierarchical extension of the Cormack-Jolly-Seber (CJS) model for open population capture-recapture data. In addition to recaptures of marked animals, we model first captures of animals and losses on capture. The parameter set includes capture probabilities, survival rates, and birth rates. The survival rates and birth rates are treated as a random sample from a bivariate distribution, thus the model explicitly incorporates correlation in these demographic rates. A key feature of the model is that the likelihood function, which includes a CJS model factor, is expressed entirely in terms of identifiable parameters; losses on capture can be factored out of the model. Since the computational complexity of classical likelihood methods is prohibitive, we use Markov chain Monte Carlo in a Bayesian analysis. We describe an efficient candidate-generation scheme for Metropolis-Hastings sampling of CJS models and extensions. The procedure is illustrated using mark-recapture data for the moth Gonodontis bidentata.