A methodology for performing global uncertainty and sensitivity analysis in systems biology

Accuracy of results from mathematical and computer models of biological systems is often complicated by the presence of uncertainties in experimental data that are used to estimate parameter values. Current mathematical modeling approaches typically use either single-parameter or local sensitivity analyses. However, these methods do not accurately assess uncertainty and sensitivity in the system as, by default, they hold all other parameters fixed at baseline values. Using techniques described within we demonstrate how a multi-dimensional parameter space can be studied globally so all uncertainties can be identified. Further, uncertainty and sensitivity analysis techniques can help to identify and ultimately control uncertainties. In this work we develop methods for applying existing analytical tools to perform analyses on a variety of mathematical and computer models. We compare two specific types of global sensitivity analysis indexes that have proven to be among the most robust and efficient. Through familiar and new examples of mathematical and computer models, we provide a complete methodology for performing these analyses, in both deterministic and stochastic settings, and propose novel techniques to handle problems encountered during these types of analyses.     

What Can We Learn from Global Sensitivity Analysis of Biochemical Systems?

Most biological models of intermediate size, and probably all large models, need to cope with the fact that many of their parameter values are unknown. In addition, it may not be possible to identify these values unambiguously on the basis of experimental data. This poses the question how reliable predictions made using such models are. Sensitivity analysis is commonly used to measure the impact of each model parameter on its variables. However, the results of such analyses can be dependent on an exact set of parameter values due to nonlinearity. To mitigate this problem, global sensitivity analysis techniques are used to calculate parameter sensitivities in a wider parameter space. We applied global sensitivity analysis to a selection of five signalling and metabolic models, several of which incorporate experimentally well-determined parameters. Assuming these models represent physiological reality, we explored how the results could change under increasing amounts of parameter uncertainty. Our results show that parameter sensitivities calculated with the physiological parameter values are not necessarily the most frequently observed under random sampling, even in a small interval around the physiological values. Often multimodal distributions were observed. Unsurprisingly, the range of possible sensitivity coefficient values increased with the level of parameter uncertainty, though the amount of parameter uncertainty at which the pattern of control was able to change differed among the models analysed. We suggest that this level of uncertainty can be used as a global measure of model robustness. Finally a comparison of different global sensitivity analysis techniques shows that, if high-throughput computing resources are available, then random sampling may actually be the most suitable technique.     

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.     

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.

     

Meta-models as a straightforward approach to the sensitivity analysis of complex models

Complex simulation models are important tools in applied ecological and conservation research. However sensitivity analysis of this important class of models can be difficult to conduct. High level interactions and non-linear responses are common in complex simulations, and this necessitates a global sensitivity analysis, where each parameter is tested at a range of values, and in combination with changes in many other parameters. We reviewed the literature, searching for population viability analyses that used simulation models. We found only 9 out of the 122 simulation population viability analysis used global sensitivity analysis. This result is typical of other simulation models in applied ecology, where global sensitivity analysis is rare. We then demonstrate how to conduct a meta-modeling sensitivity analysis, where a simpler statistically fit function (the meta-model, also known as the surrogate model or emulator) is used to approximate the behavior of the complicated simulation. This simpler meta-model is interrogated to inform on the behavior of simulation model. We fit two example meta-models, a generalized linear model and a boosted regression tree, to exemplify the approach. Our hope is that by going through these techniques thoroughly they will become more widely adopted.     

A simulation study of the genetic regulatory hierarchy for butterfly eyespot focus determination

The color patterns on the wings of butterflies have been an important model system in evolutionary developmental biology. Two types of models have been used to study these patterns. The first type of model employs computational techniques and generalized mechanisms of pattern formation to make predictions about how color patterns will vary as parameters of the model are changed. These generalized mechanisms include diffusion gradient, reaction-diffusion, lateral inhibition, and threshold responses. The second type of model uses known genetic interactions from Drosophila melanogaster and patterns of candidate gene expression in one of several butterfly species (most often Junonia (Precis) coenia or Bicyclus anynana) to propose specific genetic regulatory hierarchies that appear to be involved in color pattern formation. This study combines these two approaches using computational techniques to test proposed genetic regulatory hierarchies for the determination of butterfly eyespot foci (also known as border ocelli foci). Two computer programs, STELLA 8.1 and Delphi 2.0, were used to simulate the determination of eyespot foci. Both programs revealed weaknesses in a genetic model previously proposed for eyespot focus determination. On the basis of these simulations, we propose two revised models for eyespot focus determination and identify components of the genetic regulatory hierarchy that are particularly sensitive to changes in model parameter values. These components may play a key role in the evolution of butterfly eyespots. Simulations like these may be useful tools for the study of other evolutionary developmental model systems and reveal similar sensitive components of the relevant genetic regulatory hierarchies.     

灵芝胞外多糖分批发酵非结构动力学模型

在２Ｌ搅拌发酵罐上提出了描述了灵芝胞外多糖分批发酵过程中菌球生长、底物消耗和胞外多糖形成的非结构动力学模型。首先研究了灵芝分批发酵特性,结果表明该发酵过程属菌体生长和产物形成相偶联型。然后在总结文献的基础上,运用动力学模型,经过非线性回归,得到了模型中的参数值。通过计算机模拟,证明模型预测值与实际实验值具有良好的拟合性。     

Parameter estimation in Hodgkin-Huxley-Type equations for membrane action potentials in nerve and heart muscle

A method is presented which renders parameter estimation possible in systems of non-linear differential equations where normally no solution exists in terms of analytic functions and which have to be solved numerically. The method uses the concept of sensitivity equations. Two examples are given, taking mathematical models for membrane action potentials in nerve and heart muscle by Hodgkin and Huxley and by Beeler and Reuter. The model equations together with the corresponding system of sensitivity equations are given, which are necessary to estimate maximum conductivity coefficients defining the interactions of different ionic current components. A computer program is described and results of action potential numerical analysis are presented using simulated data. It can be seen, that even with superimposed simulated noise the real parameter values are estimated in an excellent manner. The method can be used to interpret observed changes in action potential time courses under physiological and pharmacological conditions.     

Learning-Induced Synchronization and Plasticity of a Developing Neural Network

Learning-induced synchronization of a neural network at various developing stages is studied by computer simulations using a pulse-coupled neural network model in which the neuronal activity is simulated by a one-dimensional map. Two types of Hebbian plasticity rules are investigated and their differences are compared. For both models, our simulations show a logarithmic increase in the synchronous firing frequency of the network with the culturing time of the neural network. This result is consistent with recent experimental observations. To investigate how to control the synchronization behavior of a neural network after learning, we compare the occurrence of synchronization for four networks with different designed patterns under the influence of an external signal. The effect of such a signal on the network activity highly depends on the number of connections between neurons. We discuss the synaptic plasticity and enhancement effects for a random network after learning at various developing stages.     

Conditions for the Emergence of Shared Norms in Populations with Incompatible Preferences

Understanding norms is a key challenge in sociology. Nevertheless, there is a lack of dynamical models explaining how one of several possible behaviors is established as a norm and under what conditions. Analysing an agent-based model, we identify interesting parameter dependencies that imply when two behaviors will coexist or when a shared norm will emerge in a heterogeneous society, where different populations have incompatible preferences. Our model highlights the importance of randomness, spatial interactions, non-linear dynamics, and self-organization. It can also explain the emergence of unpopular norms that do not maximize the collective benefit. Furthermore, we compare behavior-based with preference-based punishment and find interesting results concerning hypocritical punishment. Strikingly, pressuring others to perform the same public behavior as oneself is more effective in promoting norms than pressuring others to meet one’s own private preference. Finally, we show that adaptive group pressure exerted by randomly occuring, local majorities may create norms under conditions where different behaviors would normally coexist.     

Statistical model building and model criticism for human circadian data

Mathematical models have played an important role in the analysis of circadian systems. The models include simulation of differential equation systems to assess the dynamic properties of a circadian system and the use of statistical models, primarily harmonic regression methods, to assess the static properties of the system. The dynamical behaviors characterized by the simulation studies are the response of the circadian pacemaker to light, its rate of decay to its limit cycle, and its response to the rest-activity cycle. The static properties are phase, amplitude, and period of the intrinsic oscillator. Formal statistical methods are not routinely employed in simulation studies, and therefore the uncertainty in inferences based on the differential equation models and their sensitivity to model specification and parameter estimation error cannot be evaluated. The harmonic regression models allow formal statistical analysis of static but not dynamical features of the circadian pacemaker. The authors present a paradigm for analyzing circadian data based on the Box iterative scheme for statistical model building. The paradigm unifies the differential equation-based simulations (direct problem) and the model fitting approach using harmonic regression techniques (inverse problem) under a single schema. The framework is illustrated with the analysis of a core-temperature data series collected under a forced desynchrony protocol. The Box iterative paradigm provides a framework for systematically constructing and analyzing models of circadian data.     

三参数增长模型拟合:以季风常绿阔叶林中两个优势乔木种群为例

Logistic、Mitscherlich、Gompertz方程是一类三参数饱和增长曲线模型,广泛地应用于许多学科领域．本文基于logistic方程饱和值K估计的三点法、四点法,推导出Mitscherlich、Gompertz方程K值的三点法、四点法估计公式,并以南亚热带季风常绿阔叶林中两种优势乔木厚壳桂、黄果厚壳桂种群为例,先用三点法或四点法估计出K值,再通过线性回归与非线性回归相结合的方法,可获得三个增长模型中三个参数的最优无偏估计．实例研究表明,两个优势种群增长数据均符合三个增长模型,但更符合增长曲线呈S形的logistic、Gompertz方程,且以logistic方程最适合于观察;黄果厚壳桂种群增长快于厚壳桂种群．     

Estimating parameters for generalized mass action models using constraint propagation

As modern molecular biology moves towards the analysis of biological systems as opposed to their individual components, the need for appropriate mathematical and computational techniques for understanding the dynamics and structure of such systems is becoming more pressing. For example, the modeling of biochemical systems using ordinary differential equations (ODEs) based on high-throughput, time-dense profiles is becoming more common-place, which is necessitating the development of improved techniques to estimate model parameters from such data. Due to the high dimensionality of this estimation problem, straight-forward optimization strategies rarely produce correct parameter values, and hence current methods tend to utilize genetic/evolutionary algorithms to perform non-linear parameter fitting. Here, we describe a completely deterministic approach, which is based on interval analysis. This allows us to examine entire sets of parameters, and thus to exhaust the global search within a finite number of steps. In particular, we show how our method may be applied to a generic class of ODEs used for modeling biochemical systems called Generalized Mass Action Models (GMAs). In addition, we show that for GMAs our method is amenable to the technique in interval arithmetic called constraint propagation, which allows great improvement of its efficiency. To illustrate the applicability of our method we apply it to some networks of biochemical reactions appearing in the literature, showing in particular that, in addition to estimating system parameters in the absence of noise, our method may also be used to recover the topology of these networks.     

Asymmetry in neural fields: a spatiotemporal encoding mechanism

Neural field models have been successfully applied to model diverse brain mechanisms like visual attention, motor control, and memory. Most theoretical and modeling works have focused on the study of the dynamics of such systems under variations in neural connectivity, mainly symmetric connectivity among neurons. However, less attention has been given to the emerging properties of neuron populations when neural connectivity is asymmetric, although asymmetric activity propagation has been observed in cortical tissue. Here we explore the dynamics of neural fields with asymmetric connectivity and show, in the case of front propagation, that it can bias the population to follow a certain trajectory with higher activation. We find that asymmetry relates linearly to the input speed when the input is spatially localized, and this relation holds for different kernels and input shapes. To illustrate the behavior of asymmetric connectivity, we present an application: standard video sequences of human motion were encoded using the asymmetric neural field and compared to computer vision techniques. Overall, our results indicate that asymmetric neural fields are a competitive approach for spatiotemporal encoding with two main advantages: online classification and distributed operation.     

Decreasing-Rate Pruning Optimizes the Construction of Efficient and Robust Distributed Networks

Robust, efficient, and low-cost networks are advantageous in both biological and engineered systems. During neural network development in the brain, synapses are massively over-produced and then pruned-back over time. This strategy is not commonly used when designing engineered networks, since adding connections that will soon be removed is considered wasteful. Here, we show that for large distributed routing networks, network function is markedly enhanced by hyper-connectivity followed by aggressive pruning and that the global rate of pruning, a developmental parameter not previously studied by experimentalists, plays a critical role in optimizing network structure. We first used high-throughput image analysis techniques to quantify the rate of pruning in the mammalian neocortex across a broad developmental time window and found that the rate is decreasing over time. Based on these results, we analyzed a model of computational routing networks and show using both theoretical analysis and simulations that decreasing rates lead to more robust and efficient networks compared to other rates. We also present an application of this strategy to improve the distributed design of airline networks. Thus, inspiration from neural network formation suggests effective ways to design distributed networks across several domains.     

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.     

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.     

Computer simulation of biological systems

The current status of mathematical models of biological systems is reviewed. Advances in supercomputer hardware allows more complex models to be constructed. The new generation of microcomputers are quite adequate for many computer simulations of biological systems. A theory of modeling is being developed to improve the relationship between the real biological system and the model. Deterministic models, stochastic models and applications of control theory and optimization methods are discussed. Examples given include models of molecular structure, of experimental techniques, and of biochemical reactions. It is recommended that experimental biologists consider the use of microcomputers to model the system under study as a part of their research program.     

Modeling fluctuations in default-mode brain network using a spiking neural network

Recently, numerous attempts have been made to understand the dynamic behavior of complex brain systems using neural network models. The fluctuations in blood-oxygen-level-dependent (BOLD) brain signals at less than 0.1 Hz have been observed by functional magnetic resonance imaging (fMRI) for subjects in a resting state. This phenomenon is referred to as a "default-mode brain network." In this study, we model the default-mode brain network by functionally connecting neural communities composed of spiking neurons in a complex network. Through computational simulations of the model, including transmission delays and complex connectivity, the network dynamics of the neural system and its behavior are discussed. The results show that the power spectrum of the modeled fluctuations in the neuron firing patterns is consistent with the default-mode brain network's BOLD signals when transmission delays, a characteristic property of the brain, have finite values in a given range.     

Robustness: predicting the effects of life history perturbations on stage-structured population dynamics

Matrix-based models lie at the core of many applications across the physical, engineering and life sciences. In ecology, matrix models arise naturally via population projection matrices (PPM). The eigendata of PPMs provide detailed quantitative and qualitative information on the dynamic behaviour of model populations, especially their asymptotic rates of growth or decline. A fundamental task in modern ecology is to assess the effect that perturbations to life-cycle transition rates of individuals have on such eigendata. The prevailing assessment tools in ecological applications of PPMs are direct matrix simulations of eigendata and linearised extrapolations to the typically non-linear relationship between perturbation magnitude and the resulting matrix eigenvalues. In recent years, mathematical systems theory has developed an analytical framework, called 'Robustness Analysis and Robust Control', encompassing also algorithms and numerical tools. This framework provides a systematic and precise approach to studying perturbations and uncertainty in systems represented by matrices. Here we lay down the foundations and concepts for a 'robustness' inspired approach to predictive analyses in population ecology. We treat a number of application-specific perturbation problems and show how they can be formulated and analysed using these robustness methodologies.