Parameter estimation techniques for interaction and redistribution models: a predator-prey example

Summary The use of parameter estimation techniques for partial differential equations is illustrated using a predatorprey model. Whereas ecologists have often estimated parameters in models, they have not previously been able to do so for models that describe interactions in heterogeneous environments. The techniques we describe for partial differential equations will be generally useful for models of interacting species in spatially complex environments and for models that include the movement of organisms. We demonstrate our methods using field data from a ladybird beetle (Coccinella septempunctata) and aphid (Uroleucon nigrotuberculatum) interaction. Our parameter estimation algorithms can be employed to identify models that explain better than 80% of the observed variance in aphid and ladybird densities. Such parameter estimation techniques can bridge the gap between detail-rich experimental studies and abstract mathematical models. By relating the particular bestfit models identified from our experimental data to other information on Coccinella behavior, we conclude that a term describing local taxis of ladybirds towards prey (aphids in this case) is needed in the model.     

Parameter Estimation Methods for Chaotic Intercellular Networks

We have investigated simulation-based techniques for parameter estimation in chaotic intercellular networks. The proposed methodology combines a synchronization–based framework for parameter estimation in coupled chaotic systems with some state–of–the–art computational inference methods borrowed from the field of computational statistics. The first method is a stochastic optimization algorithm, known as accelerated random search method, and the other two techniques are based on approximate Bayesian computation. The latter is a general methodology for non–parametric inference that can be applied to practically any system of interest. The first method based on approximate Bayesian computation is a Markov Chain Monte Carlo scheme that generates a series of random parameter realizations for which a low synchronization error is guaranteed. We show that accurate parameter estimates can be obtained by averaging over these realizations. The second ABC–based technique is a Sequential Monte Carlo scheme. The algorithm generates a sequence of “populations”, i.e., sets of randomly generated parameter values, where the members of a certain population attain a synchronization error that is lesser than the error attained by members of the previous population. Again, we show that accurate estimates can be obtained by averaging over the parameter values in the last population of the sequence. We have analysed how effective these methods are from a computational perspective. For the numerical simulations we have considered a network that consists of two modified repressilators with identical parameters, coupled by the fast diffusion of the autoinducer across the cell membranes.     

Parameter inversion estimation in photosynthetic models: Impact of different simulation methods

When we apply ecological models in environmental management, we must assess the accuracy of parameter estimation and its impact on model predictions. Parameters estimated by conventional techniques tend to be nonrobust and require excessive computational resources. However, optimization algorithms are highly robust and generally exhibit convergence of parameter estimation by inversion with nonlinear models. They can simultaneously generate a large number of parameter estimates using an entire data set. In this study, we tested four inversion algorithms (simulated annealing, shuffled complex evolution, particle swarm optimization, and the genetic algorithm) to optimize parameters in photosynthetic models depending on different temperatures. We investigated if parameter boundary values and control variables influenced the accuracy and efficiency of the various algorithms and models. We obtained optimal solutions with all of the inversion algorithms tested if the parameter bounds and control variables were constrained properly. However, the efficiency of processing time use varied with the control variables obtained. In addition, we investigated if temperature dependence formalization impacted optimally the parameter estimation process. We found that the model with a peaked temperature response provided the best fit to the data.     

Bayesian Parameter Inference and Model Selection by Population Annealing in Systems Biology

Parameter inference and model selection are very important for mathematical modeling in systems biology. Bayesian statistics can be used to conduct both parameter inference and model selection. Especially, the framework named approximate Bayesian computation is often used for parameter inference and model selection in systems biology. However, Monte Carlo methods needs to be used to compute Bayesian posterior distributions. In addition, the posterior distributions of parameters are sometimes almost uniform or very similar to their prior distributions. In such cases, it is difficult to choose one specific value of parameter with high credibility as the representative value of the distribution. To overcome the problems, we introduced one of the population Monte Carlo algorithms, population annealing. Although population annealing is usually used in statistical mechanics, we showed that population annealing can be used to compute Bayesian posterior distributions in the approximate Bayesian computation framework. To deal with un-identifiability of the representative values of parameters, we proposed to run the simulations with the parameter ensemble sampled from the posterior distribution, named “posterior parameter ensemble”. We showed that population annealing is an efficient and convenient algorithm to generate posterior parameter ensemble. We also showed that the simulations with the posterior parameter ensemble can, not only reproduce the data used for parameter inference, but also capture and predict the data which was not used for parameter inference. Lastly, we introduced the marginal likelihood in the approximate Bayesian computation framework for Bayesian model selection. We showed that population annealing enables us to compute the marginal likelihood in the approximate Bayesian computation framework and conduct model selection depending on the Bayes factor.     

Correcting the Bias of Empirical Frequency Parameter Estimators in Codon Models

Markov models of codon substitution are powerful inferential tools for studying biological processes such as natural selection and preferences in amino acid substitution. The equilibrium character distributions of these models are almost always estimated using nucleotide frequencies observed in a sequence alignment, primarily as a matter of historical convention. In this note, we demonstrate that a popular class of such estimators are biased, and that this bias has an adverse effect on goodness of fit and estimates of substitution rates. We propose a “corrected” empirical estimator that begins with observed nucleotide counts, but accounts for the nucleotide composition of stop codons. We show via simulation that the corrected estimates outperform the de facto standard estimates not just by providing better estimates of the frequencies themselves, but also by leading to improved estimation of other parameters in the evolutionary models. On a curated collection of sequence alignments, our estimators show a significant improvement in goodness of fit compared to the approach. Maximum likelihood estimation of the frequency parameters appears to be warranted in many cases, albeit at a greater computational cost. Our results demonstrate that there is little justification, either statistical or computational, for continued use of the -style estimators.     

Parameter optimization by using differential elimination: a general approach for introducing constraints into objective functions

Background

The investigation of network dynamics is a major issue in systems and synthetic biology. One of the essential steps in a dynamics investigation is the parameter estimation in the model that expresses biological phenomena. Indeed, various techniques for parameter optimization have been devised and implemented in both free and commercial software. While the computational time for parameter estimation has been greatly reduced, due to improvements in calculation algorithms and the advent of high performance computers, the accuracy of parameter estimation has not been addressed.

Results

We propose a new approach for parameter optimization by using differential elimination, to estimate kinetic parameter values with a high degree of accuracy. First, we utilize differential elimination, which is an algebraic approach for rewriting a system of differential equations into another equivalent system, to derive the constraints between kinetic parameters from differential equations. Second, we estimate the kinetic parameters introducing these constraints into an objective function, in addition to the error function of the square difference between the measured and estimated data, in the standard parameter optimization method. To evaluate the ability of our method, we performed a simulation study by using the objective function with and without the newly developed constraints: the parameters in two models of linear and non-linear equations, under the assumption that only one molecule in each model can be measured, were estimated by using a genetic algorithm (GA) and particle swarm optimization (PSO). As a result, the introduction of new constraints was dramatically effective: the GA and PSO with new constraints could successfully estimate the kinetic parameters in the simulated models, with a high degree of accuracy, while the conventional GA and PSO methods without them frequently failed.

Conclusions

The introduction of new constraints in an objective function by using differential elimination resulted in the drastic improvement of the estimation accuracy in parameter optimization methods. The performance of our approach was illustrated by simulations of the parameter optimization for two models of linear and non-linear equations, which included unmeasured molecules, by two types of optimization techniques. As a result, our method is a promising development in parameter optimization.

     

A High-Throughput Screening Approach to Discovering Good Forms of Biologically Inspired Visual Representation

While many models of biological object recognition share a common set of “broad-stroke” properties, the performance of any one model depends strongly on the choice of parameters in a particular instantiation of that model—e.g., the number of units per layer, the size of pooling kernels, exponents in normalization operations, etc. Since the number of such parameters (explicit or implicit) is typically large and the computational cost of evaluating one particular parameter set is high, the space of possible model instantiations goes largely unexplored. Thus, when a model fails to approach the abilities of biological visual systems, we are left uncertain whether this failure is because we are missing a fundamental idea or because the correct “parts” have not been tuned correctly, assembled at sufficient scale, or provided with enough training. Here, we present a high-throughput approach to the exploration of such parameter sets, leveraging recent advances in stream processing hardware (high-end NVIDIA graphic cards and the PlayStation 3''s IBM Cell Processor). In analogy to high-throughput screening approaches in molecular biology and genetics, we explored thousands of potential network architectures and parameter instantiations, screening those that show promising object recognition performance for further analysis. We show that this approach can yield significant, reproducible gains in performance across an array of basic object recognition tasks, consistently outperforming a variety of state-of-the-art purpose-built vision systems from the literature. As the scale of available computational power continues to expand, we argue that this approach has the potential to greatly accelerate progress in both artificial vision and our understanding of the computational underpinning of biological vision.     

Spatial capture–recapture with random thinning for unidentified encounters

1. Spatial capture–recapture (SCR) models have increasingly been used as a basis for combining capture–recapture data types with variable levels of individual identity information to estimate population density and other demographic parameters. Recent examples are the unmarked SCR (or spatial count model), where no individual identities are available and spatial mark–resight (SMR) where individual identities are available for only a marked subset of the population. Currently lacking, though, is a model that allows unidentified samples to be combined with identified samples when there are no separate classes of “marked” and “unmarked” individuals and when the two sample types cannot be considered as arising from two independent observation models. This is a common scenario when using noninvasive sampling methods, for example, when analyzing data on identified and unidentified photographs or scats from the same sites.
2. Here we describe a “random thinning” SCR model that utilizes encounters of both known and unknown identity samples using a natural mechanistic dependence between samples arising from a single observation model. Our model was fitted in a Bayesian framework using NIMBLE.
3. We investigate the improvement in parameter estimates by including the unknown identity samples, which was notable (up to 79% more precise) in low‐density populations with a low rate of identified encounters. We then applied the random thinning SCR model to a noninvasive genetic sampling study of brown bear (Ursus arctos) density in Oriental Cantabrian Mountains (North Spain).
4. Our model can improve density estimation for noninvasive sampling studies for low‐density populations with low rates of individual identification, by making use of available data that might otherwise be discarded.

     

Efficient Characterization of Parametric Uncertainty of Complex (Bio)chemical Networks

Parametric uncertainty is a particularly challenging and relevant aspect of systems analysis in domains such as systems biology where, both for inference and for assessing prediction uncertainties, it is essential to characterize the system behavior globally in the parameter space. However, current methods based on local approximations or on Monte-Carlo sampling cope only insufficiently with high-dimensional parameter spaces associated with complex network models. Here, we propose an alternative deterministic methodology that relies on sparse polynomial approximations. We propose a deterministic computational interpolation scheme which identifies most significant expansion coefficients adaptively. We present its performance in kinetic model equations from computational systems biology with several hundred parameters and state variables, leading to numerical approximations of the parametric solution on the entire parameter space. The scheme is based on adaptive Smolyak interpolation of the parametric solution at judiciously and adaptively chosen points in parameter space. As Monte-Carlo sampling, it is “non-intrusive” and well-suited for massively parallel implementation, but affords higher convergence rates. This opens up new avenues for large-scale dynamic network analysis by enabling scaling for many applications, including parameter estimation, uncertainty quantification, and systems design.     

Universally sloppy parameter sensitivities in systems biology models

Quantitative computational models play an increasingly important role in modern biology. Such models typically involve many free parameters, and assigning their values is often a substantial obstacle to model development. Directly measuring in vivo biochemical parameters is difficult, and collectively fitting them to other experimental data often yields large parameter uncertainties. Nevertheless, in earlier work we showed in a growth-factor-signaling model that collective fitting could yield well-constrained predictions, even when it left individual parameters very poorly constrained. We also showed that the model had a “sloppy” spectrum of parameter sensitivities, with eigenvalues roughly evenly distributed over many decades. Here we use a collection of models from the literature to test whether such sloppy spectra are common in systems biology. Strikingly, we find that every model we examine has a sloppy spectrum of sensitivities. We also test several consequences of this sloppiness for building predictive models. In particular, sloppiness suggests that collective fits to even large amounts of ideal time-series data will often leave many parameters poorly constrained. Tests over our model collection are consistent with this suggestion. This difficulty with collective fits may seem to argue for direct parameter measurements, but sloppiness also implies that such measurements must be formidably precise and complete to usefully constrain many model predictions. We confirm this implication in our growth-factor-signaling model. Our results suggest that sloppy sensitivity spectra are universal in systems biology models. The prevalence of sloppiness highlights the power of collective fits and suggests that modelers should focus on predictions rather than on parameters.     

Statistical Modeling of Agatston Score in Multi-Ethnic Study of Atherosclerosis (MESA)

The MESA (Multi-Ethnic Study of Atherosclerosis) is an ongoing study of the prevalence, risk factors, and progression of subclinical cardiovascular disease in a multi-ethnic cohort. It provides a valuable opportunity to examine the development and progression of CAC (coronary artery calcium), which is an important risk factor for the development of coronary heart disease. In MESA, about half of the CAC scores are zero and the rest are continuously distributed. Such data has been referred to as “zero-inflated data” and may be described using two-part models. Existing two-part model studies have limitations in that they usually consider parametric models only, make the assumption of known forms of the covariate effects, and focus only on the estimation property of the models. In this article, we investigate statistical modeling of CAC in MESA. Building on existing studies, we focus on two-part models. We investigate both parametric and semiparametric, and both proportional and nonproportional models. For various models, we study their estimation as well as prediction properties. We show that, to fully describe the relationship between covariates and CAC development, the semiparametric model with nonproportional covariate effects is needed. In contrast, for the purpose of prediction, the parametric model with proportional covariate effects is sufficient. This study provides a statistical basis for describing the behaviors of CAC and insights into its biological mechanisms.     

Three algorithms for interpreting models consisting of ordinary differential equations: Sensitivity coefficients,sensitivity functions,global optimization

This paper reports the development and application of three powerful algorithms for the analysis and simulation of mathematical models consisting of ordinary differential equations. First, we describe an extended parameter sensitivity analysis: we measure the relative sensitivities of many dynamical behaviors of the model to perturbations of each parameter. We check sensitivities to parameter variation over both small and large ranges. These two extensions of a common technique have applications in parameter estimation and in experimental design. Second, we compute sensitivity functions, using an efficient algorithm requiring just one model simulation to obtain all sensitivities of state variables to all parameters as functions of time. We extend the analysis to a behavior which is not a state variable. Third, we present an unconstrained global optimization algorithm, and apply it in a novel way: we determine the input to the model, given an optimality criterion and typical outputs. The algorithm itself is an efficient one for high-order problems, and does not get stuck at local extrema. We apply the sensitivity analysis, sensitivity functions, and optimization algorithm to a sixth-order nonlinear ordinary differential equation model for human eye movements. This application shows that the algorithms are not only practicable for high-order models, but also useful as conceptual tools.     

Parameter Estimation and Model Selection in Computational Biology

A central challenge in computational modeling of biological systems is the determination of the model parameters. Typically, only a fraction of the parameters (such as kinetic rate constants) are experimentally measured, while the rest are often fitted. The fitting process is usually based on experimental time course measurements of observables, which are used to assign parameter values that minimize some measure of the error between these measurements and the corresponding model prediction. The measurements, which can come from immunoblotting assays, fluorescent markers, etc., tend to be very noisy and taken at a limited number of time points. In this work we present a new approach to the problem of parameter selection of biological models. We show how one can use a dynamic recursive estimator, known as extended Kalman filter, to arrive at estimates of the model parameters. The proposed method follows. First, we use a variation of the Kalman filter that is particularly well suited to biological applications to obtain a first guess for the unknown parameters. Secondly, we employ an a posteriori identifiability test to check the reliability of the estimates. Finally, we solve an optimization problem to refine the first guess in case it should not be accurate enough. The final estimates are guaranteed to be statistically consistent with the measurements. Furthermore, we show how the same tools can be used to discriminate among alternate models of the same biological process. We demonstrate these ideas by applying our methods to two examples, namely a model of the heat shock response in E. coli, and a model of a synthetic gene regulation system. The methods presented are quite general and may be applied to a wide class of biological systems where noisy measurements are used for parameter estimation or model selection.     

Regression Analysis for Constraining Free Parameters in Electrophysiological Models of Cardiac Cells

A major challenge in computational biology is constraining free parameters in mathematical models. Adjusting a parameter to make a given model output more realistic sometimes has unexpected and undesirable effects on other model behaviors. Here, we extend a regression-based method for parameter sensitivity analysis and show that a straightforward procedure can uniquely define most ionic conductances in a well-known model of the human ventricular myocyte. The model''s parameter sensitivity was analyzed by randomizing ionic conductances, running repeated simulations to measure physiological outputs, then collecting the randomized parameters and simulation results as “input” and “output” matrices, respectively. Multivariable regression derived a matrix whose elements indicate how changes in conductances influence model outputs. We show here that if the number of linearly-independent outputs equals the number of inputs, the regression matrix can be inverted. This is significant, because it implies that the inverted matrix can specify the ionic conductances that are required to generate a particular combination of model outputs. Applying this idea to the myocyte model tested, we found that most ionic conductances could be specified with precision (R² > 0.77 for 12 out of 16 parameters). We also applied this method to a test case of changes in electrophysiology caused by heart failure and found that changes in most parameters could be well predicted. We complemented our findings using a Bayesian approach to demonstrate that model parameters cannot be specified using limited outputs, but they can be successfully constrained if multiple outputs are considered. Our results place on a solid mathematical footing the intuition-based procedure simultaneously matching a model''s output to several data sets. More generally, this method shows promise as a tool to define model parameters, in electrophysiology and in other biological fields.     

Evolutionary Triplet Models of Structured RNA

The reconstruction and synthesis of ancestral RNAs is a feasible goal for paleogenetics. This will require new bioinformatics methods, including a robust statistical framework for reconstructing histories of substitutions, indels and structural changes. We describe a “transducer composition” algorithm for extending pairwise probabilistic models of RNA structural evolution to models of multiple sequences related by a phylogenetic tree. This algorithm draws on formal models of computational linguistics as well as the 1985 protosequence algorithm of David Sankoff. The output of the composition algorithm is a multiple-sequence stochastic context-free grammar. We describe dynamic programming algorithms, which are robust to null cycles and empty bifurcations, for parsing this grammar. Example applications include structural alignment of non-coding RNAs, propagation of structural information from an experimentally-characterized sequence to its homologs, and inference of the ancestral structure of a set of diverged RNAs. We implemented the above algorithms for a simple model of pairwise RNA structural evolution; in particular, the algorithms for maximum likelihood (ML) alignment of three known RNA structures and a known phylogeny and inference of the common ancestral structure. We compared this ML algorithm to a variety of related, but simpler, techniques, including ML alignment algorithms for simpler models that omitted various aspects of the full model and also a posterior-decoding alignment algorithm for one of the simpler models. In our tests, incorporation of basepair structure was the most important factor for accurate alignment inference; appropriate use of posterior-decoding was next; and fine details of the model were least important. Posterior-decoding heuristics can be substantially faster than exact phylogenetic inference, so this motivates the use of sum-over-pairs heuristics where possible (and approximate sum-over-pairs). For more exact probabilistic inference, we discuss the use of transducer composition for ML (or MCMC) inference on phylogenies, including possible ways to make the core operations tractable.     

Parameter estimation for bursting neural models

This paper presents work on parameter estimation methods for bursting neural models. In our approach we use both geometrical features specific to bursting, as well as general features such as periodic orbits and their bifurcations. We use the geometry underlying bursting to introduce defining equations for burst initiation and termination, and restrict the estimation algorithms to the space of bursting periodic orbits when trying to fit periodic burst data. These geometrical ideas are combined with automatic differentiation to accurately compute parameter sensitivities for the burst timing and period. In addition to being of inherent interest, these sensitivities are used in standard gradient-based optimization algorithms to fit model burst duration and period to data. As an application, we fit Butera et al.'s (Journal of Neurophysiology 81, 382-397, 1999) model of preB?tzinger complex neurons to empirical data both in control conditions and when the neuromodulator norepinephrine is added (Viemari and Ramirez, Journal of Neurophysiology 95, 2070-2082, 2006). The results suggest possible modulatory mechanisms in the preB?tzinger complex, including modulation of the persistent sodium current.     

Improved free energy parameters for RNA pseudoknotted secondary structure prediction

Accurate prediction of RNA pseudoknotted secondary structures from the base sequence is a challenging computational problem. Since prediction algorithms rely on thermodynamic energy models to identify low-energy structures, prediction accuracy relies in large part on the quality of free energy change parameters. In this work, we use our earlier constraint generation and Boltzmann likelihood parameter estimation methods to obtain new energy parameters for two energy models for secondary structures with pseudoknots, namely, the Dirks–Pierce (DP) and the Cao–Chen (CC) models. To train our parameters, and also to test their accuracy, we create a large data set of both pseudoknotted and pseudoknot-free secondary structures. In addition to structural data our training data set also includes thermodynamic data, for which experimentally determined free energy changes are available for sequences and their reference structures. When incorporated into the HotKnots prediction algorithm, our new parameters result in significantly improved secondary structure prediction on our test data set. Specifically, the prediction accuracy when using our new parameters improves from 68% to 79% for the DP model, and from 70% to 77% for the CC model.     

Parameter estimation in a stochastic model of the tubuloglomerular feedback mechanism in a rat nephron

A key parameter in the understanding of renal hemodynamics is the gain of the feedback function in the tubuloglomerular feedback mechanism. A dynamic model of autoregulation of renal blood flow and glomerular filtration rate has been extended to include a stochastic differential equations model of one of the main parameters that determines feedback gain. The model reproduces fluctuations and irregularities in the tubular pressure oscillations that the former deterministic models failed to describe. This approach assumes that the gain exhibits spontaneous erratic variations that can be explained by a variety of influences, which change over time (blood pressure, hormone levels, etc.). To estimate the key parameters of the model we have developed a new estimation method based on the oscillatory behavior of the data. The dynamics is characterized by the spectral density, which has been estimated for the observed time series, and numerically approximated for the model. The parameters have then been estimated by the least squares distance between data and model spectral densities. To evaluate the estimation procedure measurements of the proximal tubular pressure from 35 nephrons in 16 rat kidneys have been analyzed, and the parameters characterizing the gain and the delay have been estimated. There was good agreement between the estimated values, and the values obtained for the same parameters in independent, previously published experiments.     

Improving Parameter Inference from FRAP Data: an Analysis Motivated by Pattern Formation in the <Emphasis Type="BoldItalic">Drosophila</Emphasis> Wing Disc

Fluorescence recovery after photobleaching (FRAP) is used to obtain quantitative information about molecular diffusion and binding kinetics at both cell and tissue levels of organization. FRAP models have been proposed to estimate the diffusion coefficients and binding kinetic parameters of species for a variety of biological systems and experimental settings. However, it is not clear what the connection among the diverse parameter estimates from different models of the same system is, whether the assumptions made in the model are appropriate, and what the qualities of the estimates are. Here we propose a new approach to investigate the discrepancies between parameters estimated from different models. We use a theoretical model to simulate the dynamics of a FRAP experiment and generate the data that are used in various recovery models to estimate the corresponding parameters. By postulating a recovery model identical to the theoretical model, we first establish that the appropriate choice of observation time can significantly improve the quality of estimates, especially when the diffusion and binding kinetics are not well balanced, in a sense made precise later. Secondly, we find that changing the balance between diffusion and binding kinetics by changing the size of the bleaching region, which gives rise to different FRAP curves, provides a priori knowledge of diffusion and binding kinetics, which is important for model formulation. We also show that the use of the spatial information in FRAP provides better parameter estimation. By varying the recovery model from a fixed theoretical model, we show that a simplified recovery model can adequately describe the FRAP process in some circumstances and establish the relationship between parameters in the theoretical model and those in the recovery model. We then analyze an example in which the data are generated with a model of intermediate complexity and the parameters are estimated using models of greater or less complexity, and show how sensitivity analysis can be used to improve FRAP model formulation. Lastly, we show how sophisticated global sensitivity analysis can be used to detect over-fitting when using a model that is too complex.     

Parameter redundancy in discrete state‐space and integrated models

Discrete state‐space models are used in ecology to describe the dynamics of wild animal populations, with parameters, such as the probability of survival, being of ecological interest. For a particular parametrization of a model it is not always clear which parameters can be estimated. This inability to estimate all parameters is known as parameter redundancy or a model is described as nonidentifiable. In this paper we develop methods that can be used to detect parameter redundancy in discrete state‐space models. An exhaustive summary is a combination of parameters that fully specify a model. To use general methods for detecting parameter redundancy a suitable exhaustive summary is required. This paper proposes two methods for the derivation of an exhaustive summary for discrete state‐space models using discrete analogues of methods for continuous state‐space models. We also demonstrate that combining multiple data sets, through the use of an integrated population model, may result in a model in which all parameters are estimable, even though models fitted to the separate data sets may be parameter redundant.