Stochastic vs. deterministic uptake of dodecanedioic acid by isolated rat livers

Deterministic and stochastic differential equations models of the uptake of dodecanedioic acid (C12) are fitted to experimental data obtained on nine isolated, perfused rat livers. 11500 μg of C12 were injected as a bolus into the perfusing liver solution. The concentrations of C12 in perfusate samples taken over 2 h from the beginning of the experiments were analyzed by High Performance Liquid Chromatography (HPLC). A two-compartment deterministic model is studied. To include spontaneous erratic variations in the metabolic processes the parameter for the uptake rate is randomized to obtain a stochastic differential equations model. Parameters are estimated in a two-step procedure: first, parameters in the drift part are estimated by least squares; then, the diffusion parameter is estimated using Monte-Carlo simulations to approximate the unknown likelihood function. Parameter estimation is carried out over a wide range of reasonable measurement error variances to check robustness of estimates. It is concluded that the kinetics of dodecanedioic acid, in the experimental conditions discussed, is well approximated by a model including spontaneous erratic variations in the liver uptake rate.     

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.     

Calibration of dynamic models of biological systems with KInfer

Methods for parameter estimation that are robust to experimental uncertainties and to stochastic and biological noise and that require a minimum of a priori input knowledge are of key importance in computational systems biology. The new method presented in this paper aims to ensure an inference model that deduces the rate constants of a system of biochemical reactions from experimentally measured time courses of reactants. This new method was applied to some challenging parameter estimation problems of nonlinear dynamic biological systems and was tested both on synthetic and real data. The synthetic case studies are the 12-state model of the SERCA pump and a model of a genetic network containing feedback loops of interaction between regulator and effector genes. The real case studies consist of a model of the reaction between the inhibitor κB kinase enzyme and its substrate in the signal transduction pathway of NF-κB, and a stiff model of the fermentation pathway of Lactococcus lactis.     

Modelling and analysis of time-lags in some basic patterns of cell proliferation

 In this paper, we present a systematic approach for obtaining qualitatively and quantitatively correct mathematical models of some biological phenomena with time-lags. Features of our approach are the development of a hierarchy of related models and the estimation of parameter values, along with their non-linear biases and standard deviations, for sets of experimental data. We demonstrate our method of solving parameter estimation problems for neutral delay differential equations by analyzing some models of cell growth that incorporate a time-lag in the cell division phase. We show that these models are more consistent with certain reported data than the classic exponential growth model. Although the exponential growth model provides estimates of some of the growth characteristics, such as the population-doubling time, the time-lag growth models can additionally provide estimates of: (i) the fraction of cells that are dividing, (ii) the rate of commitment of cells to cell division, (iii) the initial distribution of cells in the cell cycle, and (iv) the degree of synchronization of cells in the (initial) cell population. Received: 15 September 1997/Revised version: 1 April 1998     

Global parameter estimation methods for stochastic biochemical systems

Background  

The importance of stochasticity in cellular processes having low number of molecules has resulted in the development of stochastic models such as chemical master equation. As in other modelling frameworks, the accompanying rate constants are important for the end-applications like analyzing system properties (e.g. robustness) or predicting the effects of genetic perturbations. Prior knowledge of kinetic constants is usually limited and the model identification routine typically includes parameter estimation from experimental data. Although the subject of parameter estimation is well-established for deterministic models, it is not yet routine for the chemical master equation. In addition, recent advances in measurement technology have made the quantification of genetic substrates possible to single molecular levels. Thus, the purpose of this work is to develop practical and effective methods for estimating kinetic model parameters in the chemical master equation and other stochastic models from single cell and cell population experimental data.     

Deterministic extinction effect of parasites on host populations

 Experimental studies have shown that parasites can reduce host density and even drive host population to extinction. Conventional mathematical models for parasite-host interactions, while can address the host density reduction scenario, fail to explain such deterministic extinction phenomena. In order to understand the parasite induced host extinction, Ebert et al. (2000) formulated a plausible but ad hoc epidemiological microparasite model and its stochastic variation. The deterministic model, resembles a simple SI type model, predicts the existence of a globally attractive positive steady state. Their simulation of the stochastic model indicates that extinction of host is a likely outcome in some parameter regions. A careful examination of their ad hoc model suggests an alternative and plausible model assumption. With this modification, we show that the revised parasite-host model can exhibit the observed parasite induced host extinction. This finding strengthens and complements that of Ebert et al. (2000), since all continuous models are likely break down when all population densities are small. This extinction dynamics resembles that of ratio-dependent predator-prey models. We report here a complete global study of the revised parasite-host model. Biological implications and limitations of our findings are also presented. Received: 30 October 2001 / Revised version: 11 February 2002 / Published online: 17 October 2002 Work is partially supported by NSF grant DMS-0077790 Mathematics Subject Classification (2000): 34C25, 34C35, 92D25. Keywords or phrases: Microparasite model – Ratio-dependent predator-prey model – Host extinction – Global stability – Biological control     

Multisensory fusion and the stochastic structure of postural sway

 We analyze the stochastic structure of postural sway and demonstrate that this structure imposes important constraints on models of postural control. Linear stochastic models of various orders were fit to the center-of-mass trajectories of subjects during quiet stance in four sensory conditions: (i) light touch and vision, (ii) light touch, (iii) vision, and (iv) neither touch nor vision. For each subject and condition, the model of appropriate order was determined, and this model was characterized by the eigenvalues and coefficients of its autocovariance function. In most cases, postural-sway trajectories were similar to those produced by a third-order model with eigenvalues corresponding to a slow first-order decay plus a faster-decaying damped oscillation. The slow-decay fraction, which we define as the slow-decay autocovariance coefficient divided by the total variance, was usually near 1. We compare the stochastic structure of our data to two linear control-theory models: (i) a proportional–integral–derivative control model in which the postural system's state is assumed to be known, and (ii) an optimal-control model in which the system's state is estimated based on noisy multisensory information using a Kalman filter. Under certain assumptions, both models have eigenvalues consistent with our results. However, the slow-decay fraction predicted by both models is less than we observe. We show that our results are more consistent with a modification of the optimal-control model in which noise is added to the computations performed by the state estimator. This modified model has a slow-decay fraction near 1 in a parameter regime in which sensory information related to the body's velocity is more accurate than sensory information related to position and acceleration. These findings suggest that: (i) computation noise is responsible for much of the variance observed in postural sway, and (ii) the postural control system under the conditions tested resides in the regime of accurate velocity information. Received: 20 March 2001 / Accepted: 17 April 2002 Acknowledgements. We thank Tjeerd Dijkstra for bringing the slow-decay component of postural sway to our attention. Funding for this research was provided by National Institutes of Health grant R29 N35070–01A2, John J. Jeka, PI. Correspondence to: T. Kiemel (Tel.: +1-301-4056176, Fax: +1-301-3149358 e-mail: kiemel@glue.umd.edu)     

A PDMP model of the epithelial cell turn-over in the intestinal crypt including microbiota-derived regulations

Mechanistic models are a powerful tool to gain insights into biological processes. The parameters of such models, e.g. kinetic rate constants, usually cannot be measured directly but need to be inferred from experimental data. In this article, we study dynamical models of the translation kinetics after mRNA transfection and analyze their parameter identifiability. That is, whether parameters can be uniquely determined from perfect or realistic data in theory and practice. Previous studies have considered ordinary differential equation (ODE) models of the process, and here we formulate a stochastic differential equation (SDE) model. For both model types, we consider structural identifiability based on the model equations and practical identifiability based on simulated as well as experimental data and find that the SDE model provides better parameter identifiability than the ODE model. Moreover, our analysis shows that even for those parameters of the ODE model that are considered to be identifiable, the obtained estimates are sometimes unreliable. Overall, our study clearly demonstrates the relevance of considering different modeling approaches and that stochastic models can provide more reliable and informative results.

     

A Quasistationary Analysis of a Stochastic Chemical Reaction: Keizer’s Paradox

For a system of biochemical reactions, it is known from the work of T.G. Kurtz [J. Appl. Prob. 8, 344 (1971)] that the chemical master equation model based on a stochastic formulation approaches the deterministic model based on the Law of Mass Action in the infinite system-size limit in finite time. The two models, however, often show distinctly different steady-state behavior. To further investigate this “paradox,” a comparative study of the deterministic and stochastic models of a simple autocatalytic biochemical reaction, taken from a text by the late J. Keizer, is carried out. We compute the expected time to extinction, the true stochastic steady state, and a quasistationary probability distribution in the stochastic model. We show that the stochastic model predicts the deterministic behavior on a reasonable time scale, which can be consistently obtained from both models. The transition time to the extinction, however, grows exponentially with the system size. Mathematically, we identify that exchanging the limits of infinite system size and infinite time is problematic. The appropriate system size that can be considered sufficiently large, an important parameter in numerical computation, is also discussed.     

The galvanotaxis response mechanism of keratinocytes can be modeled as a proportional controller

Human keratinocytes actively crawl in vitro when plated onto a collagen-coated glass substrate, and their direction of migration is totally random. In response to an imposed dc electric field, they migrate asymmetrically, moving mostly toward the negative pole of the field. The authors have analyzed experimental data reported by others to determine the basic characteristics of the cellular response machinery in these keratinocytes. This movement can be completely described mathematically using two independent variables: the speed, V, and the angle of migration, ϕ. The authors propose a model in which a steerer (controller without feedback) is responsible for determining the speed, and an automatic controller (controller with feedback) is responsible for determining the angle of migration. The torque to rotate is induced by a deterministic cellular signal and a stochastic cellular signal. The cellular machine characteristics are determined as follows: The angular dependence of the detection unit is sin ϕ; the detection unit detects the guiding field in a linear fashion; the cellular reaction unit can be described by a constant; the chemical amplifier, as well as the cellular motor work, is linear; the cellular characteristic time, which quantifies the cellular stochastic signal, is 50 min.     

Population extinction and quasi-stationary behavior in stochastic density-dependent structured models

Density-independent and density-dependent, stochastic and deterministic, discrete-time, structured models are formulated, analysed and numerically simulated. A special case of the deterministic, density-independent, structured model is the well-known Leslie age-structured model. The stochastic, density-independent model is a multitype branching process. A review of linear, density-independent models is given first, then nonlinear, density-dependent models are discussed. In the linear, density-independent structured models, transitions between states are independent of time and state. Population extinction is determined by the dominant eigenvalue λ of the transition matrix. If λ ≤ 1, then extinction occurs with probability one in the stochastic and deterministic models. However, if λ > 1, then the deterministic model has exponential growth, but in the stochastic model there is a positive probability of extinction which depends on the fixed point of the system of probability generating functions. The linear, density-independent, stochastic model is generalized to a nonlinear, density-dependent one. The dependence on state is in terms of a weighted total population size. It is shown for small initial population sizes that the density-dependent, stochastic model can be approximated by the density-independent, stochastic model and thus, the extinction behavior exhibited by the linear model occurs in the nonlinear model. In the deterministic models there is a unique stable equilibrium. Given the population does not go extinct, it is shown that the stochastic model has a quasi-stationary distribution with mean close to the stable equilibrium, provided the population size is sufficiently large. For small values of the population size, complete extinction can be observed in the simulations. However, the persistence time increases rapidly with the population size. This author received partial support by the National Science Foundation grant # DMS-9626417.     

Exciting chaos with noise: unexpected dynamics in epidemic outbreaks

 In this paper, we identify a mechanism for chaos in the presence of noise. In a study of the SEIR model, which predicts epidemic outbreaks in childhood diseases, we show how chaotic dynamics can be attained by adding stochastic perturbations at parameters where chaos does not exist apriori. Data recordings of epidemics in childhood diseases are still argued as deterministic chaos. There also exists noise due to uncertainties in the contact parameters between those who are susceptible and those who are infected, as well as random fluctuations in the population. Although chaos has been found in deterministic models, it only occurs in parameter regions that require a very large population base or other large seasonal forcing. Our work identifies the mechanism whereby chaos can be induced by noise for realistic parameter regions of the deterministic model where it does not naturally occur. Received: 13 October 2000 / Revised version: 15 May 2001 / Published online: 7 December 2001     

Consequences of chemosensory phenomena for leukocyte chemotactic orientation

The stochastic nature of cell surface receptor-ligand binding is known to limit the accuracy of detection of chemoattractant gradients by leukocytes (11, 12), thus limiting the orientation ability that is crucial to the chemotactic response in host defense. The probabilistic cell orientation model of Lauffenburger (11) is extended here to assess the consequences of recently discovered receptor phenomena: “down-regulation” of total surface receptor number, spatial asymmetry of surface receptors, and existence of a higher-affinity receptor subpopulation. In general, a reduction in orientation accuracy is predicted by inclusion of these phenomena. An orientation signal based on a simple model of chemosensory adaptation (i.e., a spatial difference inrelative receptor occupancy) is found to be functionally different from the signal suggested by an experimental correlation (i.e., a spatial difference inabsolute receptor occupancy). However, in the context of receptor “signal noise,” the signal based on adaptation yields predictions in better qualitative agreement with the experimental orientation data of Zigmond (10). From this cell orientation model we can estimate the effective timeaveraging period required for noise diminution to a level allowing orientation predictions to match observed levels. This time-averaging period presumably reflects the time constant for receptor signal transduction and locomotory response.     

An autoregressive method for the measurement of synchronization of interictal and ictal EEG signals

We propose a new measure of synchronization of multichannel ictal and interictal EEG signals. The measure is based on the residual covariance matrix of a multichannel autoregressive model. A major advantage of this measure is its ability to be interpreted both in the framework of stochastic and deterministic models. A preliminary analysis of EEG data from three patients using this measure documents the expected increased synchronization during ictal periods but also reveals that increased synchrony persists for prolonged periods (up to 2 h or more) in the postictal period. Received: 20 July 1997 / Accepted in revised form: 26 January 1999     

Mathematical modelling of the intravenous glucose tolerance test

 Several attempts at building a satisfactory model of the glucose-insulin system are recorded in the literature. The minimal model, which is the model currently mostly used in physiological research on the metabolism of glucose, was proposed in the early eighties for the interpretation of the glucose and insulin plasma concentrations following the intravenous glucose tolerance test. It is composed of two parts: the first consists of two differential equations and describes the glucose plasma concentration time-course treating insulin plasma concentration as a known forcing function; the second consists of a single equation and describes the time course of plasma insulin concentration treating glucose plasma concentration as a known forcing function. The two parts are to be separately estimated on the available data. In order to study glucose-insulin homeostasis as a single dynamical system, a unified model would be desirable. To this end, the simple coupling of the original two parts of the minimal model is not appropriate, since it can be shown that, for commonly observed combinations of parameter values, the coupled model would not admit an equilibrium and the concentration of active insulin in the “distant” compartment would be predicted to increase without bounds. For comparison, a simple delay-differential model is introduced, is demonstrated to be globally asymptotically stable around a unique equilibrium point corresponding to the pre-bolus conditions, and is shown to have positive and bounded solutions for all times. The results of fitting the delay-differential model to experimental data from ten healthy volunteers are also shown. It is concluded that a global unified model is both theoretically desirable and practically usable, and that any such model ought to undergo formal analysis to establish its appropriateness and to exclude conflicts with accepted physiological notions. Received: 22 June 1998 / Revised version: 24 February 1999     

A Stochastic Model Correctly Predicts Changes in Budding Yeast Cell Cycle Dynamics upon Periodic Expression of CLN2

In this study, we focus on a recent stochastic budding yeast cell cycle model. First, we estimate the model parameters using extensive data sets: phenotypes of 110 genetic strains, single cell statistics of wild type and cln3 strains. Optimization of stochastic model parameters is achieved by an automated algorithm we recently used for a deterministic cell cycle model. Next, in order to test the predictive ability of the stochastic model, we focus on a recent experimental study in which forced periodic expression of CLN2 cyclin (driven by MET3 promoter in cln3 background) has been used to synchronize budding yeast cell colonies. We demonstrate that the model correctly predicts the experimentally observed synchronization levels and cell cycle statistics of mother and daughter cells under various experimental conditions (numerical data that is not enforced in parameter optimization), in addition to correctly predicting the qualitative changes in size control due to forced CLN2 expression. Our model also generates a novel prediction: under frequent CLN2 expression pulses, G1 phase duration is bimodal among small-born cells. These cells originate from daughters with extended budded periods due to size control during the budded period. This novel prediction and the experimental trends captured by the model illustrate the interplay between cell cycle dynamics, synchronization of cell colonies, and size control in budding yeast.     

Analysis of the reflexive feedback control loop during posture maintenance

 In previous work it has been shown in posture experiments of the human arm that reflexive dynamics were substantial for narrow-band stochastic force disturbances. The estimated reflex gains varied substantially with the frequency content of the disturbances. The present study analyses a simplified linear model of the reflexive feedback control loop, to provide an explanation for the observed behaviour. The model describes co-activation and reflexive feedback. The task instruction `minimize the displacements' is represented mathematically by a cost function that is minimized by adjusting the parameters of the model. Small-amplitude displacements allow the system to be analysed with a quasi-linear approach. The optimization results clarify the limited effectiveness of reflexive feedback on the system's closed-loop behaviour, which emanates from the time delay present in the reflex loops. For low-frequency inputs less than 5 Hz, boundary-stable solutions with high reflex gains are predicted to be optimal. Input frequencies near the system's eigenfrequency (about 5 Hz), however, would be amplified and result in oscillatory behaviour. As long as the disturbance does not excite these frequencies, boundary stability will be optimal. The predicted reflex gains show a striking similarity with the estimated reflex gains from the experimental study. The present model analysis also provides a clear explanation for the negative reflex gains, estimated for near-sinusoidal inputs beyond 1.5 Hz. Received: 24 January 2000 / Accepted in revised form: 7 July 2000     

Exterior exposure estimation using a one-compartment toxicokinetic model with blood sample measurements

Exposure assessment of individuals exposed to certain chemicals plays an important role in the analysis of occupational—as well as environmental-health problems. Biological monitoring, as an alternative to direct environmental measurements, may be applied to relate the exterior exposure with the amount of individual intake. In this paper, we estimate individuals’ (inhalation) exposure retrospectively from their blood concentrations via a simplified one-compartment toxicokinetic model. Considering stochastic variations to the toxicokinetic model, the solution to the resultant stochastic differential equation (SDE), together with measurement error, is transformed into a dynamic linear state-space model. The unknown model parameters and the mean inhalation concentration are then estimated via Markov Chain Monte Carlo (MCMC) simulations. The proposed method is used in the analysis of the styrene data (Wang et al. in Occup Environ Med 53:601–605, 1996) to backward estimate the inhalation concentration, assuming it is unknown. The data analysis showed that the internal stochastic variations, often ignored in toxicokinetic model analysis, outweighed in standard deviation almost twice that of the measurement error. Also, the simulation results showed that the method performed relatively well to the approach considering measurement error only.

---

Contract/grant sponsor: National Science Council of Taiwan (NSC 93-2118-M-032-004); National Health Research Institutes of Taiwan (BS-096-PP-11).     

An alternative stochastic model of generation of oligodendrocytes in cell culture

According to our previous model, oligodendrocyte – type 2 (O-2A) astrocyte progenitor cells become competent for differentiation in vitro after they complete a certain number of critical mitotic cycles. After attaining the competency to differentiate, progenitor cells divide with fixed probability p in subsequent cycles. The number of critical cycles is random; analysis of data suggests that it varies from zero to two. The present paper presents an alternative model in which there are no critical cycles, and the probability that a progenitor cell will divide again decreases gradually to a plateau value as the number of completed mitotic cycles increases. In particular all progenitor cells have the ability to differentiate from the time of plating. The Kiefer-Wolfowitz procedure is used to fit the new model to experimental data on the clonal growth of purified O-2A progenitor cells obtained from the optic nerves of 7 day old rats. The new model is shown to fit the experimental data well, indicating that it is not possible to determine whether critical cycles exist on the basis of these experimental data. In contrast to the fit of the previous model, which suggested that the addition of thyroid hormone increased the limiting probability of differentiation as the number of mitotic cycles increases, the fit of the new model suggests that the addition of thyroid hormone has almost no effect on the limiting probability of differentiation. Received: 6 March 2000 / Revised version: 18 September 2000 / Published online: 30 April 2001     

Stochastic resonance in psychophysics and in animal behavior

 A recent analysis of the energy detector model in sensory psychophysics concluded that stochastic resonance does not occur in a measure of signal detectability (d′), but can occur in a percent-correct measure of performance as an epiphenomenon of nonoptimal criterion placement [Tougaard (2000) Biol Cybern 83: 471–480]. When generalized to signal detection in sensory systems in general, this conclusion is a serious challenge to the idea that stochastic resonance could play a significant role in sensory processing in humans and other animals. It also seems to be inconsistent with recent demonstrations of stochastic resonance in sensory systems of both nonhuman animals and humans using measures of system performance such as signal-to-noise ratio of power spectral densities and percent-correct detections in a two-interval forced-choice paradigm, both closely related to d′. In this paper we address this apparent dilemma by discussing several models of how stochastic resonance can arise in signal detection systems, including especially those that implement a “soft threshold” at the input transform stage. One example involves redefining d′ for energy increments in terms of parameters of the spike-count distribution of FitzHugh–Nagumo neurons. Another involves a Poisson spike generator that receives an exponentially transformed noisy periodic signal. In this case it can be shown that the signal-to-noise ratio of the power spectral density at the signal frequency, which exhibits stochastic resonance, is proportional to d′. Finally, a variant of d′ is shown to exhibit stochastic resonance when calculated directly from the distributions of power spectral densities at the signal frequency resulting from transformation of noise alone and a noisy signal by a sufficiently steep nonlinear response function. All of these examples, and others from the literature, imply that stochastic resonance is more than an epiphenomenon, although significant limitations to the extent to which adding noise can aid detection do exist. Received: 22 January 2001 / Accepted in revised form: 8 March 2002