Local Perturbation Analysis: A Computational Tool for Biophysical Reaction-Diffusion Models

Diffusion and interaction of molecular regulators in cells is often modeled using reaction-diffusion partial differential equations. Analysis of such models and exploration of their parameter space is challenging, particularly for systems of high dimensionality. Here, we present a relatively simple and straightforward analysis, the local perturbation analysis, that reveals how parameter variations affect model behavior. This computational tool, which greatly aids exploration of the behavior of a model, exploits a structural feature common to many cellular regulatory systems: regulators are typically either bound to a membrane or freely diffusing in the interior of the cell. Using well-documented, readily available bifurcation software, the local perturbation analysis tracks the approximate early evolution of an arbitrarily large perturbation of a homogeneous steady state. In doing so, it provides a bifurcation diagram that concisely describes various regimes of the model’s behavior, reducing the need for exhaustive simulations to explore parameter space. We explain the method and provide detailed step-by-step guides to its use and application.     

‘Glocal’ Robustness Analysis and Model Discrimination for Circadian Oscillators

To characterize the behavior and robustness of cellular circuits with many unknown parameters is a major challenge for systems biology. Its difficulty rises exponentially with the number of circuit components. We here propose a novel analysis method to meet this challenge. Our method identifies the region of a high-dimensional parameter space where a circuit displays an experimentally observed behavior. It does so via a Monte Carlo approach guided by principal component analysis, in order to allow efficient sampling of this space. This ‘global’ analysis is then supplemented by a ‘local’ analysis, in which circuit robustness is determined for each of the thousands of parameter sets sampled in the global analysis. We apply this method to two prominent, recent models of the cyanobacterial circadian oscillator, an autocatalytic model, and a model centered on consecutive phosphorylation at two sites of the KaiC protein, a key circadian regulator. For these models, we find that the two-sites architecture is much more robust than the autocatalytic one, both globally and locally, based on five different quantifiers of robustness, including robustness to parameter perturbations and to molecular noise. Our ‘glocal’ combination of global and local analyses can also identify key causes of high or low robustness. In doing so, our approach helps to unravel the architectural origin of robust circuit behavior. Complementarily, identifying fragile aspects of system behavior can aid in designing perturbation experiments that may discriminate between competing mechanisms and different parameter sets.     

Systemic properties of ensembles of metabolic networks: application of graphical and statistical methods to simple unbranched pathways

MOTIVATION: Mathematical models are the only realistic method for representing the integrated dynamic behavior of complex biochemical networks. However, it is difficult to obtain a consistent set of values for the parameters that characterize such a model. Even when a set of parameter values exists, the accuracy of the individual values is questionable. Therefore, we were motivated to explore statistical techniques for analyzing the properties of a given model when knowledge of the actual parameter values is lacking. RESULTS: The graphical and statistical methods presented in the previous paper are applied here to simple unbranched biosynthetic pathways subject to control by feedback inhibition. We represent these pathways within a canonical nonlinear formalism that provides a regular structure that is convenient for randomly sampling the parameter space. After constructing a large ensemble of randomly generated sets of parameter values, the structural and behavioral properties of the model with these parameter sets are examined statistically and classified. The results of our analysis demonstrate that certain properties of these systems are strongly correlated, thereby revealing aspects of organization that are highly probable independent of selection. Finally, we show how specification of a given behavior affects the distribution of acceptable parameter values.     

Mathematical identification of critical reactions in the interlocked feedback model

Dynamic simulations are necessary for understanding the mechanism of how biochemical networks generate robust properties to environmental stresses or genetic changes. Sensitivity analysis allows the linking of robustness to network structure. However, it yields only local properties regarding a particular choice of plausible parameter values, because it is hard to know the exact parameter values in vivo. Global and firm results are needed that do not depend on particular parameter values. We propose mathematical analysis for robustness (MAR) that consists of the novel evolutionary search that explores all possible solution vectors of kinetic parameters satisfying the target dynamics and robustness analysis. New criteria, parameter spectrum width and the variability of solution vectors for parameters, are introduced to determine whether the search is exhaustive. In robustness analysis, in addition to single parameter sensitivity analysis, robustness to multiple parameter perturbation is defined. Combining the sensitivity analysis and the robustness analysis to multiple parameter perturbation enables identifying critical reactions. Use of MAR clearly identified the critical reactions responsible for determining the circadian cycle in the Drosophila interlocked circadian clock model. In highly robust models, while the parameter vectors are greatly varied, the critical reactions with a high sensitivity are uniquely determined. Interestingly, not only the per-tim loop but also the dclk-cyc loop strongly affect the period of PER, although the dclk-cyc loop hardly changes its amplitude and it is not potentially influential. In conclusion, MAR is a powerful method to explore wide parameter space without human-biases and to link a robust property to network architectures without knowing the exact parameter values. MAR identifies the reactions critically responsible for determining the period and amplitude in the interlocked feedback model and suggests that the circadian clock intensively evolves or designs the kinetic parameters so that it creates a highly robust cycle.     

Flux module decomposition for parameter estimation in a multiple-feedback loop model of biochemical networks

Computer simulation is an important technique to capture the dynamics of biochemical networks. Since few quantitative values are measured in vivo, the values for unmeasured parameters should be estimated so that the simulation agrees with the experimental data. Considering the sparsity and error rates of experimentally measured data, the first thing is not to find a numerically exact and global solution but to explore a variety of the plausible parameter solutions. To find many plausible parameter solutions without any biases, we developed the two-phase search (TPS) method. However, calculation complexity makes it hard for TPS to optimize a large-scale dynamic model. In this study divide-and-conquer methods are used to solve this problem. The flux module decomposition (FMD) is first proposed that separates a complex, large-scale dynamic model into multiple flux modules without deteriorating its basic control architectures. FMD is combined with TPS, named FMD-TPS, to find many plausible parameter solutions for a dynamic model. To demonstrate the feasibility of FMD-TPS, it is applied to the E. coli ammonia assimilation system that consists of multiple-feedback loops. The variability of the solutions is verified by measuring the space distribution of the parameter solution vectors and by defining the binary vectors checking the consistency with biological behaviors. Compared with non-decomposition methods, FMD-TPS efficiently explored a variety of plausible parameter solutions that reproduce the dynamic behaviors in vivo.     

A data-informed mean-field approach to mapping of cortical parameter landscapes

Constraining the many biological parameters that govern cortical dynamics is computationally and conceptually difficult because of the curse of dimensionality. This paper addresses these challenges by proposing (1) a novel data-informed mean-field (MF) approach to efficiently map the parameter space of network models; and (2) an organizing principle for studying parameter space that enables the extraction biologically meaningful relations from this high-dimensional data. We illustrate these ideas using a large-scale network model of the Macaque primary visual cortex. Of the 10-20 model parameters, we identify 7 that are especially poorly constrained, and use the MF algorithm in (1) to discover the firing rate contours in this 7D parameter cube. Defining a “biologically plausible” region to consist of parameters that exhibit spontaneous Excitatory and Inhibitory firing rates compatible with experimental values, we find that this region is a slightly thickened codimension-1 submanifold. An implication of this finding is that while plausible regimes depend sensitively on parameters, they are also robust and flexible provided one compensates appropriately when parameters are varied. Our organizing principle for conceptualizing parameter dependence is to focus on certain 2D parameter planes that govern lateral inhibition: Intersecting these planes with the biologically plausible region leads to very simple geometric structures which, when suitably scaled, have a universal character independent of where the intersections are taken. In addition to elucidating the geometry of the plausible region, this invariance suggests useful approximate scaling relations. Our study offers, for the first time, a complete characterization of the set of all biologically plausible parameters for a detailed cortical model, which has been out of reach due to the high dimensionality of parameter space.     

Efficient parametric analysis of the chemical master equation through model order reduction

ABSTRACT: BACKGROUND: Stochastic biochemical reaction networks are commonly modelled by the chemical master equation, and can be simulated as first order linear differential equations through a finite state projection. Due to the very high state space dimension of these equations, numerical simulations are computationally expensive. This is a particular problem for analysis tasks requiring repeated simulations for different parameter values. Such tasks are computationally expensive to the point of infeasibility with the chemical master equation. RESULTS: In this article, we apply parametric model order reduction techniques in order to construct accurate low-dimensional parametric models of the chemical master equation. These surrogate models can be used in various parametric analysis task such as identifiability analysis, parameter estimation, or sensitivity analysis. As biological examples, we consider two models for gene regulation networks, a bistable switch and a network displaying stochastic oscillations. CONCLUSIONS: The results show that the parametric model reduction yields efficient models of stochastic biochemical reaction networks, and that these models can be useful for systems biology applications involving parametric analysis problems such as parameter exploration, optimization, estimation or sensitivity analysis.     

A model of late long-term potentiation simulates aspects of memory maintenance

Late long-term potentiation (L-LTP) denotes long-lasting strengthening of synapses between neurons. L-LTP appears essential for the formation of long-term memory, with memories at least partly encoded by patterns of strengthened synapses. How memories are preserved for months or years, despite molecular turnover, is not well understood. Ongoing recurrent neuronal activity, during memory recall or during sleep, has been hypothesized to preferentially potentiate strong synapses, preserving memories. This hypothesis has not been evaluated in the context of a mathematical model representing ongoing activity and biochemical pathways important for L-LTP. In this study, ongoing activity was incorporated into two such models - a reduced model that represents some of the essential biochemical processes, and a more detailed published model. The reduced model represents synaptic tagging and gene induction simply and intuitively, and the detailed model adds activation of essential kinases by Ca(2+). Ongoing activity was modeled as continual brief elevations of Ca(2+). In each model, two stable states of synaptic strength/weight resulted. Positive feedback between synaptic weight and the amplitude of ongoing Ca(2+) transients underlies this bistability. A tetanic or theta-burst stimulus switches a model synapse from a low basal weight to a high weight that is stabilized by ongoing activity. Bistability was robust to parameter variations in both models. Simulations illustrated that prolonged periods of decreased activity reset synaptic strengths to low values, suggesting a plausible forgetting mechanism. However, episodic activity with shorter inactive intervals maintained strong synapses. Both models support experimental predictions. Tests of these predictions are expected to further understanding of how neuronal activity is coupled to maintenance of synaptic strength. Further investigations that examine the dynamics of activity and synaptic maintenance can be expected to help in understanding how memories are preserved for up to a lifetime in animals including humans.     

Scaling in biochemical kinetics: dissection of a relaxation oscillator

In a case study of a class of biochemical models, methods of scaling are employed to determine the expected limits of validity of singular perturbation approaches to a relaxation oscillator. This work complements earlier analysis of the use of the quasi-steady-state approximation for the Michaelis-Menten approximation. These studies present an advantage over more conventional approaches in which attention is concentrated on a single parameter, in that the range of convergence is delineated more precisely in the full parameter space of the problem.     

From global to local: exploring the relationship between parameters and behaviors in models of electrical excitability

Models of electrical activity in excitable cells involve nonlinear interactions between many ionic currents. Changing parameters in these models can produce a variety of activity patterns with sometimes unexpected effects. Further more, introducing new currents will have different effects depending on the initial parameter set. In this study we combined global sampling of parameter space and local analysis of representative parameter sets in a pituitary cell model to understand the effects of adding K ⁺ conductances, which mediate some effects of hormone action on these cells. Global sampling ensured that the effects of introducing K ⁺ conductances were captured across a wide variety of contexts of model parameters. For each type of K ⁺ conductance we determined the types of behavioral transition that it evoked. Some transitions were counterintuitive, and may have been missed without the use of global sampling. In general, the wide range of transitions that occurred when the same current was applied to the model cell at different locations in parameter space highlight the challenge of making accurate model predictions in light of cell-to-cell heterogeneity. Finally, we used bifurcation analysis and fast/slow analysis to investigate why specific transitions occur in representative individual models. This approach relies on the use of a graphics processing unit (GPU) to quickly map parameter space to model behavior and identify parameter sets for further analysis. Acceleration with modern low-cost GPUs is particularly well suited to exploring the moderate-sized (5-20) parameter spaces of excitable cell and signaling models.     

The evolution of aggressive losers

We examine the question of when aggressive behavior of likely losers should be part of an evolutionarily stable strategy. We modified an earlier model by the authors that found situations where likely losers initiate aggressive interactions more often than likely winners. The modifications allowed us to examine the robustness of the previous study by including an unusually high number of possible strategies (n=81) and to examine a wide range of parameter settings. First, we show that restricting attention to only a few most plausible strategies may change the overall results. Second, within the space where escalation is predicted, for a large percentage of the parameter settings (85%), an ESS exists that leads to a somewhat counterintuitive situation where escalation is more often initiated by the likely loser than by the likely winner of the contest. In contrast, an ESS that favors escalation by likely winners was found only for about 3% of parameter settings. Furthermore, we use simulations of evolution in a finite population to verify for certain parameter settings that the analytically predicted ESS's could in fact evolve. Our results suggest that ESSs in which the likely loser rather than the likely winner is expected to initiate escalation are generic and ESSs in which the opposite is true need to be explained by incorporating specific features of the biology of a given species into more detailed models.     

Comprehensive,Population-Based Sensitivity Analysis of a Two-Mass Vocal Fold Model

Previous vocal fold modeling studies have generally focused on generating detailed data regarding a narrow subset of possible model configurations. These studies can be interpreted to be the investigation of a single subject under one or more vocal conditions. In this study, a broad population-based sensitivity analysis is employed to examine the behavior of a virtual population of subjects and to identify trends between virtual individuals as opposed to investigating a single subject or model instance. Four different sensitivity analysis techniques were used in accomplishing this task. Influential relationships between model input parameters and model outputs were identified, and an exploration of the model’s parameter space was conducted. Results indicate that the behavior of the selected two-mass model is largely dominated by complex interactions, and that few input-output pairs have a consistent effect on the model. Results from the analysis can be used to increase the efficiency of optimization routines of reduced-order models used to investigate voice abnormalities. Results also demonstrate the types of challenges and difficulties to be expected when applying sensitivity analyses to more complex vocal fold models. Such challenges are discussed and recommendations are made for future studies.     

Improving evolutionary models of protein interaction networks

MOTIVATION: Theoretical models of biological networks are valuable tools in evolutionary inference. Theoretical models based on gene duplication and divergence provide biologically plausible evolutionary mechanics. Similarities found between empirical networks and their theoretically generated counterpart are considered evidence of the role modeled mechanics play in biological evolution. However, the method by which these models are parameterized can lead to questions about the validity of the inferences. Selecting parameter values in order to produce a particular topological value obfuscates the possibility that the model may produce a similar topology for a large range of parameter values. Alternately, a model may produce a large range of topologies, allowing (incorrect) parameter values to produce a valid topology from an otherwise flawed model. In order to lend biological credence to the modeled evolutionary mechanics, parameter values should be derived from the empirical data. Furthermore, recent work indicates that the timing and fate of gene duplications are critical to proper derivation of these parameters. RESULTS: We present a methodology for deriving evolutionary rates from empirical data that is used to parameterize duplication and divergence models of protein interaction network evolution. Our method avoids shortcomings of previous methods, which failed to consider the effect of subsequent duplications. From our parameter values, we find that concurrent and existing existing duplication and divergence models are insufficient for modeling protein interaction network evolution. We introduce a model enhancement based on heritable interaction sites on the surface of a protein and find that it more closely reflects the high clustering found in the empirical network.     

Bifurcation analysis of continuous biochemical reactor models

The validity of a biochemical reactor model often is evaluated by comparing transient responses to experimental data. Dynamic simulation can be a rather inefficient and ineffective tool for analyzing bioreactor models that exhibit complex nonlinear behavior. Bifurcation analysis is a powerful tool for obtaining a more efficient and complete characterization of the model behavior. To illustrate the power of bifurcation analysis, the steady-state and transient behavior of three continuous bioreactor models consisting of a small number of ordinary differential equations are investigated. Several important features, as well as potential limitations, that are difficult to ascertain via dynamic simulation are disclosed through the bifurcation analysis. The results motivate the use of dynamic simulation and bifurcation analysis as complementary tools for analyzing the nonlinear behavior of bioreactor models.     

RECOMMENDATIONS FOR USING MSBAYES TO INCORPORATE UNCERTAINTY IN SELECTING AN ABC MODEL PRIOR: A RESPONSE TO OAKS ET AL.

Prior specification is an essential component of parameter estimation and model comparison in Approximate Bayesian computation (ABC). Oaks et al. present a simulation‐based power analysis of msBayes and conclude that msBayes has low power to detect genuinely random divergence times across taxa, and suggest the cause is Lindley's paradox. Although the predictions are similar, we show that their findings are more fundamentally explained by insufficient prior sampling that arises with poorly chosen wide priors that critically undersample nonsimultaneous divergence histories of high likelihood. In a reanalysis of their data on Philippine Island vertebrates, we show how this problem can be circumvented by expanding upon a previously developed procedure that accommodates uncertainty in prior selection using Bayesian model averaging. When these procedures are used, msBayes supports recent divergences without support for synchronous divergence in the Oaks et al. data and we further present a simulation analysis that demonstrates that msBayes can have high power to detect asynchronous divergence under narrower priors for divergence time. Our findings highlight the need for exploration of plausible parameter space and prior sampling efficiency for ABC samplers in high dimensions. We discus potential improvements to msBayes and conclude that when used appropriately with model averaging, msBayes remains an effective and powerful tool.     

Accelerated maximum likelihood parameter estimation for stochastic biochemical systems

ABSTRACT: BACKGROUND: A prerequisite for the mechanistic simulation of a biochemical system is detailed knowledge of its kinetic parameters. Despite recent experimental advances, the estimation of unknown parameter values from observed data is still a bottleneck for obtaining accurate simulation results. Many methods exist for parameter estimation in deterministic biochemical systems; methods for discrete stochastic systems are less well developed. Given the probabilistic nature of stochastic biochemical models, a natural approach is to choose parameter values that maximize the probability of the observed data with respect to the unknown parameters, a.k.a. the maximum likelihood parameter estimates (MLEs). MLE computation for all but the simplest models requires the simulation of many system trajectories that are consistent with experimental data. For models with unknown parameters, this presents a computational challenge, as the generation of consistent trajectories can be an extremely rare occurrence. RESULTS: We have developed Monte Carlo Expectation-Maximization with Modified Cross-Entropy Method (MCEM2): an accelerated method for calculating MLEs that combines advances in rare event simulation with a computationally efficient version of the Monte Carlo expectation-maximization (MCEM) algorithm. Our method requires no prior knowledge regarding parameter values, and it automatically provides a multivariate parameter uncertainty estimate. We applied the method to five stochastic systems of increasing complexity, progressing from an analytically tractable pure-birth model to a computationally demanding model of yeast-polarization. Our results demonstrate that MCEM2 substantially accelerates MLE computation on all tested models when compared to a stand-alone version of MCEM. Additionally, we show how our method identifies parameter values for certain classes of models more accurately than two recently proposed computationally efficient methods. CONCLUSIONS: This work provides a novel, accelerated version of a likelihood-based parameter estimation method that can be readily applied to stochastic biochemical systems. In addition, our results suggest opportunities for added efficiency improvements that will further enhance our ability to mechanistically simulate biological processes.     

A simple model can unify a broad range of phenomena in retinotectal map development

A paradigm model system for studying the development of patterned connections in the nervous system is the topographic map formed by retinal axons in the optic tectum/superior colliculus. Starting in the 1970s, a series of computational models have been proposed to explain map development in both normal conditions, and perturbed conditions where the retina and/or tectum/superior colliculus are altered. This stands in contrast to more recent models that have often been simpler than older ones, and tend to address more limited data sets, but include more recent genetic manipulations. The original exploration of many of the early models was one-dimensional and limited by the computational resources of the time. This leaves open the ability of these early models to explain both map development in two dimensions, and the genetic manipulation data that have only appeared more recently. In this article, we show that a two-dimensional and updated version of the XBAM model (eXtended Branch Arrow Model), first proposed in 1982, reproduces a range of surgical map manipulations not yet demonstrated by more modern models. A systematic exploration of the parameter space of this model in two dimensions also reveals richer behavior than that apparent from the original one-dimensional versions. Furthermore, we show that including a specific type of axon?Caxon interaction can account for the map collapse recently observed when particular receptor levels are genetically manipulated in a subset of retinal ganglion cells. Together these results demonstrate that balancing multiple influences on map development seems to be necessary to explain many biological phenomena in retinotectal map formation, and suggest important constraints on the underlying biological variables.     

Genetic noise control via protein oligomerization

Background

Biochemical equilibria are usually modeled iteratively: given one or a few fitted models, if there is a lack of fit or over fitting, a new model with additional or fewer parameters is then fitted, and the process is repeated. The problem with this approach is that different analysts can propose and select different models and thus extract different binding parameter estimates from the same data. An alternative is to first generate a comprehensive standardized list of plausible models, and to then fit them exhaustively, or semi-exhaustively.

Results

A framework is presented in which equilibriums are modeled as pairs (g, h) where g = 0 maps total reactant concentrations (system inputs) into free reactant concentrations (system states) which h then maps into expected values of measurements (system outputs). By letting dissociation constants K _d be either freely estimated, infinity, zero, or equal to other K _d , and by letting undamaged protein fractions be either freely estimated or 1, many g models are formed. A standard space of g models for ligand-induced protein dimerization equilibria is given. Coupled to an h model, the resulting (g, h) were fitted to dTTP induced R1 dimerization data (R1 is the large subunit of ribonucleotide reductase). Models with the fewest parameters were fitted first. Thereafter, upon fitting a batch, the next batch of models (with one more parameter) was fitted only if the current batch yielded a model that was better (based on the Akaike Information Criterion) than the best model in the previous batch (with one less parameter). Within batches models were fitted in parallel. This semi-exhaustive approach yielded the same best models as an exhaustive model space fit, but in approximately one-fifth the time.

Conclusion

Comprehensive model space based biochemical equilibrium model selection methods are realizable. Their significance to systems biology as mappings of data into mathematical models warrants their development.