A dynamic model for functional mapping of biological rhythms

Functional mapping is a statistical method for mapping quantitative trait loci (QTLs) that regulate the dynamic pattern of a biological trait. This method integrates mathematical aspects of biological complexity into a mixture model for genetic mapping and tests the genetic effects of QTLs by comparing genotype-specific curve parameters. As a way of quantitatively specifying the dynamic behavior of a system, differential equations have proven to be powerful for modeling and unraveling the biochemical, molecular, and cellular mechanisms of a biological process, such as biological rhythms. The equipment of functional mapping with biologically meaningful differential equations provides new insights into the genetic control of any dynamic processes. We formulate a new functional mapping framework for a dynamic biological rhythm by incorporating a group of ordinary differential equations (ODE). The Runge-Kutta fourth order algorithm was implemented to estimate the parameters that define the system of ODE. The new model will find its implications for understanding the interplay between gene interactions and developmental pathways in complex biological rhythms.     

Continuous modeling of metabolic networks with gene regulation in yeast and in vivo determination of rate parameters

A continuous model of a metabolic network including gene regulation to simulate metabolic fluxes during batch cultivation of yeast Saccharomyces cerevisiae was developed. The metabolic network includes reactions of glycolysis, gluconeogenesis, glycerol and ethanol synthesis and consumption, the tricarboxylic acid cycle, and protein synthesis. Carbon sources considered were glucose and then ethanol synthesized during growth on glucose. The metabolic network has 39 fluxes, which represent the action of 50 enzymes and 64 genes and it is coupled with a gene regulation network which defines enzyme synthesis (activities) and incorporates regulation by glucose (enzyme induction and repression), modeled using ordinary differential equations. The model includes enzyme kinetics, equations that follow both mass-action law and transport as well as inducible, repressible, and constitutive enzymes of metabolism. The model was able to simulate a fermentation of S. cerevisiae during the exponential growth phase on glucose and the exponential growth phase on ethanol using only one set of kinetic parameters. All fluxes in the continuous model followed the behavior shown by the metabolic flux analysis (MFA) obtained from experimental results. The differences obtained between the fluxes given by the model and the fluxes determined by the MFA do not exceed 25% in 75% of the cases during exponential growth on glucose, and 20% in 90% of the cases during exponential growth on ethanol. Furthermore, the adjustment of the fermentation profiles of biomass, glucose, and ethanol were 95%, 95%, and 79%, respectively. With these results the simulation was considered successful. A comparison between the simulation of the continuous model and the experimental data of the diauxic yeast fermentation for glucose, biomass, and ethanol, shows an extremely good match using the parameters found. The small discrepancies between the fluxes obtained through MFA and those predicted by the differential equations, as well as the good match between the profiles of glucose, biomass, and ethanol, and our simulation, show that this simple model, that does not rely on complex kinetic expressions, is able to capture the global behavior of the experimental data. Also, the determination of parameters using a straightforward minimization technique using data at only two points in time was sufficient to produce a relatively accurate model. Thus, even with a small amount of experimental data (rates and not concentrations) it was possible to estimate the parameters minimizing a simple objective function. The method proposed allows the obtention of reasonable parameters and concentrations in a system with a much larger number of unknowns than equations. Hence a contribution of this study is to present a convenient way to find in vivo rate parameters to model metabolic and genetic networks under different conditions.     

A variational approach to parameter estimation in ordinary differential equations

ABSTRACT: BACKGROUND: Ordinary differential equations are widely-used in the field of systems biology andchemical engineering to model chemical reaction networks. Numerous techniques havebeen developed to estimate parameters like rate constants, initial conditions or steady stateconcentrations from time-resolved data. In contrast to this countable set of parameters, theestimation of entire courses of network components corresponds to an innumerable set ofparameters. RESULTS: The approach presented in this work is able to deal with course estimation for extrinsicsystem inputs or intrinsic reactants, both not being constrained by the reaction networkitself. Our method is based on variational calculus which is carried out analytically toderive an augmented system of differential equations including the unconstrainedcomponents as ordinary state variables. Finally, conventional parameter estimation isapplied to the augmented system resulting in a combined estimation of courses andparameters. CONCLUSIONS: The combined estimation approach takes the uncertainty in input courses correctly intoaccount. This leads to precise parameter estimates and correct confidence intervals. Inparticular this implies that small motifs of large reaction networks can be analysedindependently of the rest. By the use of variational methods, elements from control theoryand statistics are combined allowing for future transfer of methods between the two fields.     

Stochastic simulations on a model of circadian rhythm generation

Biological phenomena are often modeled by differential equations, where states of a model system are described by continuous real values. When we consider concentrations of molecules as dynamical variables for a set of biochemical reactions, we implicitly assume that numbers of the molecules are large enough so that their changes can be regarded as continuous and they are described deterministically. However, for a system with small numbers of molecules, changes in their numbers are apparently discrete and molecular noises become significant. In such cases, models with deterministic differential equations may be inappropriate, and the reactions must be described by stochastic equations. In this study, we focus a clock gene expression for a circadian rhythm generation, which is known as a system involving small numbers of molecules. Thus it is appropriate for the system to be modeled by stochastic equations and analyzed by methodologies of stochastic simulations. The interlocked feedback model proposed by Ueda et al. as a set of deterministic ordinary differential equations provides a basis of our analyses. We apply two stochastic simulation methods, namely Gillespie's direct method and the stochastic differential equation method also by Gillespie, to the interlocked feedback model. To this end, we first reformulated the original differential equations back to elementary chemical reactions. With those reactions, we simulate and analyze the dynamics of the model using two methods in order to compare them with the dynamics obtained from the original deterministic model and to characterize dynamics how they depend on the simulation methodologies.     

Estimating a predator-prey dynamical model with the parameter cascades method

Summary .   Ordinary differential equations (ODEs) are widely used in ecology to describe the dynamical behavior of systems of interacting populations. However, systems of ODEs rarely provide quantitative solutions that are close to real field observations or experimental data because natural systems are subject to environmental and demographic noise and ecologists are often uncertain about the correct parameterization. In this article we introduce "parameter cascades" as an improved method to estimate ODE parameters such that the corresponding ODE solutions fit the real data well. This method is based on the modified penalized smoothing with the penalty defined by ODEs and a generalization of profiled estimation, which leads to fast estimation and good precision for ODE parameters from noisy data. This method is applied to a set of ODEs originally developed to describe an experimental predator–prey system that undergoes oscillatory dynamics. The new parameterization considerably improves the fit of the ODE model to the experimental data sets. At the same time, our method reveals that important structural assumptions that underlie the original ODE model are essentially correct. The mathematical formulations of the two nonlinear interaction terms (functional responses) that link the ODEs in the predator–prey model are validated by estimating the functional responses nonparametrically from the real data. We suggest two major applications of "parameter cascades" to ecological modeling: It can be used to estimate parameters when original data are noisy, missing, or when no reliable priori estimates are available; it can help to validate the structural soundness of the mathematical modeling approach.     

A dynamic model for functional mapping of biological rhythms

Functional mapping is a statistical method for mapping quantitative trait loci (QTLs) that regulate the dynamic pattern of a biological trait. This method integrates mathematical aspects of biological complexity into a mixture model for genetic mapping and tests the genetic effects of QTLs by comparing genotype-specific curve parameters. As a way of quantitatively specifying the dynamic behaviour of a system, differential equations have proved to be powerful for modelling and unravelling the biochemical, molecular, and cellular mechanisms of a biological process, such as biological rhythms. The equipment of functional mapping with biologically meaningful differential equations provides new insights into the genetic control of any dynamic processes. We formulate a new functional mapping framework for a dynamic biological rhythm by incorporating a group of ordinary differential equations (ODE). The Runge–Kutta fourth-order algorithm was implemented to estimate the parameters that define the system of ODE. The new model will find its implications for understanding the interplay between gene interactions and developmental pathways in complex biological rhythms.     

Transport Effects on Multiple-Component Reactions in Optical Biosensors

Optical biosensors are often used to measure kinetic rate constants associated with chemical reactions. Such instruments operate in the surface–volume configuration, in which ligand molecules are convected through a fluid-filled volume over a surface to which receptors are confined. Currently, scientists are using optical biosensors to measure the kinetic rate constants associated with DNA translesion synthesis—a process critical to DNA damage repair. Biosensor experiments to study this process involve multiple interacting components on the sensor surface. This multiple-component biosensor experiment is modeled with a set of nonlinear integrodifferential equations (IDEs). It is shown that in physically relevant asymptotic limits these equations reduce to a much simpler set of ordinary differential equations (ODEs). To verify the validity of our ODE approximation, a numerical method for the IDE system is developed and studied. Results from the ODE model agree with simulations of the IDE model, rendering our ODE model useful for parameter estimation.     

Immune network behavior—I. From stationary states to limit cycle oscilations

We develop a model for the idiotypic interaction between two B cell clones. This model takes into account B cell proliferation, B cell maturation, antibody production, the formation and subsequent elimination of antibody-antibody complexes and recirculation of antibodies between the spleen and the blood. Here we investigate, by means of stability and bifurcation analysis, how each of the processes influences the model's behavior. After appropriate nondimensinalization, the model consists of eight ordinary differential equations and a number of parameters. We estimate the parameters from experimental sources. Using a coordinate system that exploits the pairwise symmetry of the interactions between two clones, we analyse two simplified forms of the model and obtain bifurcation diagrams showing how their five equilibrium states are related. We show that the so-called immune states lose stability if B cell and antibody concentrations change on different time scales. Additionally, we derive the structure of stable and unstable manifolds of saddle-tye equilibria, pinpoint their (global) bifurcations and show that these bifurcations play a crucial role in determining the parameter regimes in which the model exhibits oscillatory behavior.     

Development of an Agent-Based Model (ABM) to Simulate the Immune System and Integration of a Regression Method to Estimate the Key ABM Parameters by Fitting the Experimental Data

Agent-based models (ABM) and differential equations (DE) are two commonly used methods for immune system simulation. However, it is difficult for ABM to estimate key parameters of the model by incorporating experimental data, whereas the differential equation model is incapable of describing the complicated immune system in detail. To overcome these problems, we developed an integrated ABM regression model (IABMR). It can combine the advantages of ABM and DE by employing ABM to mimic the multi-scale immune system with various phenotypes and types of cells as well as using the input and output of ABM to build up the Loess regression for key parameter estimation. Next, we employed the greedy algorithm to estimate the key parameters of the ABM with respect to the same experimental data set and used ABM to describe a 3D immune system similar to previous studies that employed the DE model. These results indicate that IABMR not only has the potential to simulate the immune system at various scales, phenotypes and cell types, but can also accurately infer the key parameters like DE model. Therefore, this study innovatively developed a complex system development mechanism that could simulate the complicated immune system in detail like ABM and validate the reliability and efficiency of model like DE by fitting the experimental data.     

Dynamical estimation of neuron and network properties II: path integral Monte Carlo methods

Hodgkin–Huxley (HH) models of neuronal membrane dynamics consist of a set of nonlinear differential equations that describe the time-varying conductance of various ion channels. Using observations of voltage alone we show how to estimate the unknown parameters and unobserved state variables of an HH model in the expected circumstance that the measurements are noisy, the model has errors, and the state of the neuron is not known when observations commence. The joint probability distribution of the observed membrane voltage and the unobserved state variables and parameters of these models is a path integral through the model state space. The solution to this integral allows estimation of the parameters and thus a characterization of many biological properties of interest, including channel complement and density, that give rise to a neuron’s electrophysiological behavior. This paper describes a method for directly evaluating the path integral using a Monte Carlo numerical approach. This provides estimates not only of the expected values of model parameters but also of their posterior uncertainty. Using test data simulated from neuronal models comprising several common channels, we show that short (<50 ms) intracellular recordings from neurons stimulated with a complex time-varying current yield accurate and precise estimates of the model parameters as well as accurate predictions of the future behavior of the neuron. We also show that this method is robust to errors in model specification, supporting model development for biological preparations in which the channel expression and other biophysical properties of the neurons are not fully known.     

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.

     

Function approximation approach to the inference of reduced NGnet models of genetic networks

Background  

The inference of a genetic network is a problem in which mutual interactions among genes are deduced using time-series of gene expression patterns. While a number of models have been proposed to describe genetic regulatory networks, this study focuses on a set of differential equations since it has the ability to model dynamic behavior of gene expression. When we use a set of differential equations to describe genetic networks, the inference problem can be defined as a function approximation problem. On the basis of this problem definition, we propose in this study a new method to infer reduced NGnet models of genetic networks.     

A new method of identifying a systems based on the minimal quadratic discrepancy criterion for biophysics problems

A wide range of biophysical systems are described by nonlinear dynamic models mathematically presented as a set of ordinary differential equations in the Cauchy explicit form: [formula: see text] Fij(X1(t),..,XN(t),t), (i = 1,...,N, j = 1,...,M), where Fij (X1(t), ..., XN(t), t) is a set of basis functions satisfying the Lipschitz condition. We investigate the problem of evaluation of model constants aij (the system identification) using experimental data about the time dependence of the dynamic parameters of the system Xi(t). A new method of system identification for the class of similar nonlinear dynamic models is proposed. It is shown that the problem of identifying an initial nonlinear model can be reduced to the solution of a system of linear equations for the matrix of the dynamic model constants [aj]i. It is proposed to determine the set of dynamic model constants aij using the criterion of minimal quadratic discrepancy for the time dependence of the set of dynamic parameters Xi(t). An important special case of the nonlinear model, the quadratic model, is considered. Test problems of identification using this method are presented for two nonlinear systems: the Van der Pol type multiparametric nonlinear oscillator and the strange attractor of Ressler, a widely known example of dynamic systems showing the stochastic behavior.     

A mathematical model for chronic myelogenous leukemia (CML) and T cell interaction

In this paper, we propose and analyse a mathematical model for chronic myelogenous leukemia (CML), a cancer of the blood. We model the interaction between naive T cells, effector T cells, and CML cancer cells in the body, using a system of ordinary differential equations which gives rates of change of the three cell populations. One of the difficulties in modeling CML is the scarcity of experimental data which can be used to estimate parameters values. To compensate for the resulting uncertainties, we use Latin hypercube sampling (LHS) on large ranges of possible parameter values in our analysis. A major goal of this work is the determination of parameters which play a critical role in remission or clearance of the cancer in the model. Our analysis examines 12 parameters, and identifies two of these, the growth and death rates of CML, as critical to the outcome of the system. Our results indicate that the most promising research avenues for treatments of CML should be those that affect these two significant parameters (CML growth and death rates), while altering the other parameters should have little effect on the outcome.     

Vulnerabilities in the tau network and the role of ultrasensitive points in tau pathophysiology

The multifactorial nature of disease motivates the use of systems-level analyses to understand their pathology. We used a systems biology approach to study tau aggregation, one of the hallmark features of Alzheimer's disease. A mathematical model was constructed to capture the current state of knowledge concerning tau's behavior and interactions in cells. The model was implemented in silico in the form of ordinary differential equations. The identifiability of the model was assessed and parameters were estimated to generate two cellular states: a population of solutions that corresponds to normal tau homeostasis and a population of solutions that displays aggregation-prone behavior. The model of normal tau homeostasis was robust to perturbations, and disturbances in multiple processes were required to achieve an aggregation-prone state. The aggregation-prone state was ultrasensitive to perturbations in diverse subsets of networks. Tau aggregation requires that multiple cellular parameters are set coordinately to a set of values that drive pathological assembly of tau. This model provides a foundation on which to build and increase our understanding of the series of events that lead to tau aggregation and may ultimately be used to identify critical intervention points that can direct the cell away from tau aggregation to aid in the treatment of tau-mediated (or related) aggregation diseases including Alzheimer's.     

Multistability in an age-structured model of hematopoiesis: Cyclical neutropenia

Cyclical neutropenia (CN) is a rare hematopoietic disorder in which the patient's neutrophil level drops to extremely low levels for a few days approximately every three weeks. CN is effectively treated with granulocyte colony stimulating factor (G-CSF), which is known to interfere with apoptosis in neutrophil precursors and to consequently increase the circulating neutrophil level. However, G-CSF treatment usually fails to eliminate the oscillation. In this study, we establish an age-structured model of hematopoiesis, which reduces to a set of four delay differential equations with specific forms of initial functions. We numerically investigate the possible stable solutions of the model equations with respect to changes in the parameters as well as the initial conditions. The results show that the hematopoietic system possesses multistability for parameters typical of the normal healthy state. From our numerical results, decreasing the proliferation rate of neutrophil precursors or increasing the stem cell death rate are two possible mechanisms to induce cyclical neutropenia, and the periods of the resulting oscillations are independent of the changing parameters. We also discuss the dependence of the model solution on the initial condition at normal parameter values corresponding to a healthy state. Using insight from our results we design a hybrid treatment method that is able to abolish the oscillations in CN.     

Mathematical description of a bursting pacemaker neuron by a modification of the Hodgkin-Huxley equations.

Modifications based on experimental results reported in the literature are made to the Hodgkin-Huxley equations to describe the electrophysiological behavior of the Aplysia abdominal ganglion R15 cell. The system is then further modified to describe the effects with the application of the drug tetrodotoxin (TTX) to the cells' bathing medium. Methods of the qualitative theory of differential equations are used to determine the conditions necessary for such a system of equations to have an oscillatory solution. A model satisfying these conditions is shown to preduct many experimental observations of R15 cell behavior. Numerical solutions are obtained for differential equations satisfying the conditions of the model. These solutions are shown to have a form similar to that of the bursting which is characteristic of this cell, and to preduct many results of experiments conducted on this cell. The physiological implications of the model are discussed.