Stochastic approaches for modelling in vivo reactions

Phenotype prediction in regulated metabolic networks

Background

Due to the growing amount of biological knowledge that is incorporated into metabolic network models, their analysis has become more and more challenging. Here, we examine the capabilities of the recently introduced chemical organization theory (OT) to ease this task. Considering only network stoichiometry, the theory allows the prediction of all potentially persistent species sets and therewith rigorously relates the structure of a network to its potential dynamics. By this, the phenotypes implied by a metabolic network can be predicted without the need for explicit knowledge of the detailed reaction kinetics.

Results

We propose an approach to deal with regulation – and especially inhibitory interactions – in chemical organization theory. One advantage of this approach is that the metabolic network and its regulation are represented in an integrated way as one reaction network. To demonstrate the feasibility of this approach we examine a model by Covert and Palsson (J Biol Chem, 277(31), 2002) of the central metabolism of E. coli that incorporates the regulation of all involved genes. Our method correctly predicts the known growth phenotypes on 16 different substrates. Without specific assumptions, organization theory correctly predicts the lethality of knockout experiments in 101 out of 116 cases. Taking into account the same model specific assumptions as in the regulatory flux balance analysis (rFBA) by Covert and Palsson, the same performance is achieved (106 correctly predicted cases). Two model specific assumptions had to be considered: first, we have to assume that secreted molecules do not influence the regulatory system, and second, that metabolites with increasing concentrations indicate a lethal state.

Conclusion

The introduced approach to model a metabolic network and its regulation in an integrated way as one reaction network makes organization analysis a universal technique to study the potential behavior of biological network models. Applying multiple methods like OT and rFBA is shown to be valuable to uncover critical assumptions and helps to improve model coherence.

A stochastic automaton shows how enzyme assemblies may contribute to metabolic efficiency

BACKGROUND

The advantages of grouping enzymes into metabolons and into higher order structures have long been debated. To quantify these advantages, we have developed a stochastic automaton that allows experiments to be performed in a virtual bacterium with both a membrane and a cytoplasm. We have investigated the general case of transport and metabolism as inspired by the phosphoenolpyruvate:sugar phosphotransferase system (PTS) for glucose importation and by glycolysis.

RESULTS

We show that PTS and glycolytic metabolons can increase production of pyruvate eightfold at low concentrations of phosphoenolpyruvate. A fourfold increase in the numbers of enzyme EI led to a 40% increase in pyruvate production, similar to that observed in vivo in the presence of glucose. Although little improvement resulted from the assembly of metabolons into a hyperstructure, such assembly can generate gradients of metabolites and signaling molecules.

CONCLUSION

in silico experiments may be performed successfully using stochastic automata such as HSIM (Hyperstructure Simulator) to help answer fundamental questions in metabolism about the properties of molecular assemblies and to devise strategies to modify such assemblies for biotechnological ends.

Discrete stochastic simulation of cell signaling: comparison of computational tools

Several stochastic simulation tools have been developed recently for cell signaling. A comparative evaluation of the stochastic simulation tools is needed to highlight the current state of the development. In our study, we have chosen to evaluate three stochastic simulation tools: Dizzy, Systems Biology Toolbox, and Copasi, using our own MATLAB implementation as a benchmark. The Gillespie stochastic simulation algorithm is used in all tests. With all the tools, we are able to simulate stochastically the behavior of the selected test case and to produce similar results as our own MATLAB implementation. However, it is not possible to use time-dependent inputs in stochastic simulations in Systems Biology Toolbox and Copasi. The present study is one of the first evaluations of stochastic simulation tools for realistic signal transduction pathways.

Effects of intrinsic and extrinsic noise can accelerate juxtacrine pattern formation

Epithelial pattern formation is an important phenomenon that, for example, has roles in embryogenesis, development and wound-healing. The ligand Epithelial Growth Factor (EGF) and its receptor EGF-R, constitute a system that forms lateral induction patterns by juxtacrine signalling-binding of membrane-bound ligands to receptors on neighbouring cells. Owen et al. developed a generic ordinary differential equation model of juxtacrine lateral induction that exhibits stable patterning under some conditions. The model predicts relatively slow pattern formation. We examine here the effects of both intrinsic and extrinsic cellular noise arising from the stochastic treatment of this model, and show that this noise could have an accelerating effect on the patterning process.

Stochastic binding of Ca2+ ions in the dyadic cleft; continuous versus random walk description of diffusion

Ca(2+) signaling in the dyadic cleft in ventricular myocytes is fundamentally discrete and stochastic. We study the stochastic binding of single Ca(2+) ions to receptors in the cleft using two different models of diffusion: a stochastic and discrete Random Walk (RW) model, and a deterministic continuous model. We investigate whether the latter model, together with a stochastic receptor model, can reproduce binding events registered in fully stochastic RW simulations. By evaluating the continuous model goodness-of-fit for a large range of parameters, we present evidence that it can. Further, we show that the large fluctuations in binding rate observed at the level of single time-steps are integrated and smoothed at the larger timescale of binding events, which explains the continuous model goodness-of-fit. With these results we demonstrate that the stochasticity and discreteness of the Ca(2+) signaling in the dyadic cleft, determined by single binding events, can be described using a deterministic model of Ca(2+) diffusion together with a stochastic model of the binding events, for a specific range of physiological relevant parameters. Time-consuming RW simulations can thus be avoided. We also present a new analytical model of bimolecular binding probabilities, which we use in the RW simulations and the statistical analysis.

Mathematical models for somite formation

Somitogenesis is the process of division of the anterior-posterior vertebrate embryonic axis into similar morphological units known as somites. These segments generate the prepattern which guides formation of the vertebrae, ribs and other associated features of the body trunk. In this work, we review and discuss a series of mathematical models which account for different stages of somite formation. We begin by presenting current experimental information and mechanisms explaining somite formation, highlighting features which will be included in the models. For each model we outline the mathematical basis, show results of numerical simulations, discuss their successes and shortcomings and avenues for future exploration. We conclude with a brief discussion of the state of modeling in the field and current challenges which need to be overcome in order to further our understanding in this area.

Modelling reaction kinetics inside cells

In the past decade, advances in molecular biology such as the development of non-invasive single molecule imaging techniques have given us a window into the intricate biochemical activities that occur inside cells. In this chapter we review four distinct theoretical and simulation frameworks: (i) non-spatial and deterministic, (ii) spatial and deterministic, (iii) non-spatial and stochastic and (iv) spatial and stochastic. Each framework can be suited to modelling and interpreting intracellular reaction kinetics. By estimating the fundamental length scales, one can roughly determine which models are best suited for the particular reaction pathway under study. We discuss differences in prediction between the four modelling methodologies. In particular we show that taking into account noise and space does not simply add quantitative predictive accuracy but may also lead to qualitatively different physiological predictions, unaccounted for by classical deterministic models.

Mobile trap algorithm for zinc detection using protein sensors

We present a mobile trap algorithm to sense zinc ions using protein-based sensors such as carbonic anhydrase (CA). Zinc is an essential biometal required for mammalian cellular functions although its intracellular concentration is reported to be very low. Protein-based sensors like CA molecules are employed to sense rare species like zinc ions. In this study, the zinc ions are mobile targets, which are sought by the mobile traps in the form of sensors. Particle motions are modeled using random walk along with the first passage technique for efficient simulations. The association reaction between sensors and ions is incorporated using a probability (p1) upon an ion-sensor collision. The dissociation reaction of an ion-bound CA molecule is modeled using a second, independent probability (p2). The results of the algorithm are verified against the traditional simulation techniques (e.g., Gillespie's algorithm). This study demonstrates that individual sensor molecules can be characterized using the probability pair (p1,p2), which, in turn, is linked to the system level chemical kinetic constants, kon and koff. Further investigations of CA-Zn reaction using the mobile trap algorithm show that when the diffusivity of zinc ions approaches that of sensor molecules, the reaction data obtained using the static trap assumption differ from the reaction data obtained using the mobile trap formulation. This study also reveals similar behavior when the sensor molecule has higher dissociation constant. In both the cases, the reaction data obtained using the static trap formulation reach equilibrium at a higher number of complex molecules (ion-bound sensor molecules) compared to the reaction data from the mobile trap formulation. With practical limitations on the number sensors that can be inserted/expressed in a cell and stochastic nature of the intracellular ionic concentrations, fluorescence from the number of complex sensor molecules at equilibrium will be the measure of the intracellular ion concentration. For reliable detection of zinc ions, it is desirable that the sensors must not bind all the zinc ions tightly, but should rather bind and unbind. Thus for a given fluorescence and with association-dissociation reactions between ions and sensors, the static trap approach will underestimate the number of zinc ions present in the system.

Modeling and stochastic simulation of the Ras/cAMP/PKA pathway in the yeast Saccharomyces cerevisiae evidences a key regulatory function for intracellular guanine nucleotides pools

In the yeast Saccharomyces cerevisiae, the Ras/cAMP/PKA pathway is involved in the regulation of metabolism and cell cycle progression. The pathway is tightly regulated by several control mechanisms, as the feedback cycle ruled by the activity of phosphodiesterase. Here, we present a discrete mathematical model for the Ras/cAMP/PKA pathway that considers its principal cytoplasmic components and their mutual interactions. The tau-leaping algorithm is then used to perform stochastic simulations of the model. We investigate this system under various conditions, and we test how different values of several stochastic reaction constants affect the pathway behaviour. Finally, we show that the level of guanine nucleotides, GTP and GDP, could be relevant metabolic signals for the regulation of the whole pathway.

Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics

In cell biology, cell signaling pathway problems are often tackled with deterministic temporal models, well mixed stochastic simulators, and/or hybrid methods. But, in fact, three dimensional stochastic spatial modeling of reactions happening inside the cell is needed in order to fully understand these cell signaling pathways. This is because noise effects, low molecular concentrations, and spatial heterogeneity can all affect the cellular dynamics. However, there are ways in which important effects can be accounted without going to the extent of using highly resolved spatial simulators (such as single-particle software), hence reducing the overall computation time significantly. We present a new coarse grained modified version of the next subvolume method that allows the user to consider both diffusion and reaction events in relatively long simulation time spans as compared with the original method and other commonly used fully stochastic computational methods. Benchmarking of the simulation algorithm was performed through comparison with the next subvolume method and well mixed models (MATLAB), as well as stochastic particle reaction and transport simulations (CHEMCELL, Sandia National Laboratories). Additionally, we construct a model based on a set of chemical reactions in the epidermal growth factor receptor pathway. For this particular application and a bistable chemical system example, we analyze and outline the advantages of our presented binomial tau-leap spatial stochastic simulation algorithm, in terms of efficiency and accuracy, in scenarios of both molecular homogeneity and heterogeneity.

Open-system nonequilibrium steady state: statistical thermodynamics, fluctuations, and chemical oscillations

Gibbsian equilibrium statistical thermodynamics is the theoretical foundation for isothermal, closed chemical, and biochemical reaction systems. This theory, however, is not applicable to most biochemical reactions in living cells, which exhibit a range of interesting phenomena such as free energy transduction, temporal and spatial complexity, and kinetic proofreading. In this article, a nonequilibrium statistical thermodynamic theory based on stochastic kinetics is introduced, mainly through a series of examples: single-molecule enzyme kinetics, nonlinear chemical oscillation, molecular motor, biochemical switch, and specificity amplification. The case studies illustrate an emerging theory for the isothermal nonequilibrium steady state of open systems.

Family tree of Markov models in systems biology

Motivated by applications in systems biology, a probabilistic framework based on Markov processes is proposed to represent intracellular processes. The formal relationships between different stochastic models referred to in the systems biology literature are reviewed. As part of this review, a novel derivation of the differential Chapman-Kolmogorov equation for a general multidimensional Markov process made up of both continuous and jump processes, is presented. First, the definition of a time-derivative for a probability density is focused, but placing no restrictions on the probability distribution, in particular, it is not assumed to be to be confined to a region that has a surface (on which the probability is zero). In this derivation, the master equation gives the jump part of the Markov process and the Fokker-Planck equation gives the continuous part. As a result, a 'family tree' for stochastic models in systems biology is sketched, providing explicit derivations of their formal relationship and clarifying assumptions involved.

A systems biology approach to the blood–aluminium problem: The application and testing of a computational model

Transport and distribution of systemic aluminium are influenced by its interaction with blood. Current understanding is centred upon the role played by the iron transport protein transferrin which has been shown to bind up to 90% of serum total aluminium. We have coined what we have called the blood-aluminium problem which states that the proportion of serum aluminium which, at any one moment in time, is bound by transferrin is more heavily influenced by kinetic constraints than thermodynamic equilibria with the result that the role played by transferrin in the transport and distribution of aluminium is likely to have been over estimated. To begin to solve the blood-aluminium problem and therewith provide a numerical solution to the aforementioned kinetic constraints we have applied and tested a simple computational model of the time-dependency of a putative transferrin ligand (L) binding aluminium to form an Al-L complex with a probability of existence, K(E), between 0% (no complex) and 100% (complex will not dissociate). The model is based upon the principles of a lattice-gas automaton which when ran for K(E) in the range 0.1-98.0% demonstrated the emergence of complex behaviour which could be defined in the terms of a set of parameters (equilibrium value, E(V), equilibrium time, E(T), peak value, P(V), peak time, P(T), area under curve, AUC) the values of which varied in a predictable way with K(E). When K(E) was set to 98% the model predicted that ca. 90% of the total aluminium would be bound by transferrin within ca. 350 simulation timesteps. We have used a systems biology approach to develop a simple model of the time-dependency of the binding of aluminium by transferrin. To use this approach to begin to solve the blood-aluminium problem we shall need to increase the complexity of the model to better reflect the heterogeneity of a biological system such as the blood.

The interplay between discrete noise and nonlinear chemical kinetics in a signal amplification cascade

We used various analytical and numerical techniques to elucidate signal propagation in a small enzymatic cascade which is subjected to external and internal noises. The nonlinear character of catalytic reactions, which underlie protein signal transduction cascades, renders stochastic signaling dynamics in cytosol biochemical networks distinct from the usual description of stochastic dynamics in gene regulatory networks. For a simple two-step enzymatic cascade which underlies many important protein signaling pathways, we demonstrated that the commonly used techniques such as the linear noise approximation and the Langevin equation become inadequate when the number of proteins becomes too low. Consequently, we developed a new analytical approximation, based on mixing the generating function and distribution function approaches, to the solution of the master equation that describes nonlinear chemical signaling kinetics for this important class of biochemical reactions. Our techniques work in a much wider range of protein number fluctuations than the methods used previously. We found that under certain conditions the burst phase noise may be injected into the downstream signaling network dynamics, resulting possibly in unusually large macroscopic fluctuations. In addition to computing first and second moments, which is the goal of commonly used analytical techniques, our new approach provides the full time-dependent probability distributions of the colored non-Gaussian processes in a nonlinear signal transduction cascade.

Grand canonical Markov model: a stochastic theory for open nonequilibrium biochemical networks

In this paper we present the results of a stochastic model of reversible biochemical reaction networks that are being driven through an open boundary, such that the system is interacting with its surrounding environment with explicit material exchange. The stochastic model is based on the master equation approach and is intimately related to the grand canonical ensemble of statistical mechanics. We show that it is possible to analytically calculate the joint probability function of the random variables describing the number of molecules in each state of the system for general linear networks. Definitions of reaction chemical potentials and conductances follow from inherent properties of this model, providing a description of energy dissipation in the system. We are also able to suggest novel methods for experimentally determining reaction fluxes and biochemical affinities at nonequilibrium steady state as well as the overall network connectivity.

Living with noisy genes: how cells function reliably with inherent variability in gene expression

Within a population of genetically identical cells there can be significant variation, or noise, in gene expression. Yet even with this inherent variability, cells function reliably. This review focuses on our understanding of noise at the level of both single genes and genetic regulatory networks, emphasizing comparisons between theoretical models and experimental results whenever possible. To highlight the importance of noise, we particularly emphasize examples in which a stochastic description of gene expression leads to a qualitatively different outcome than a deterministic one.

Developing Itô stochastic differential equation models for neuronal signal transduction pathways

Mathematical modeling and simulation of dynamic biochemical systems are receiving considerable attention due to the increasing availability of experimental knowledge of complex intracellular functions. In addition to deterministic approaches, several stochastic approaches have been developed for simulating the time-series behavior of biochemical systems. The problem with stochastic approaches, however, is the larger computational time compared to deterministic approaches. It is therefore necessary to study alternative ways to incorporate stochasticity and to seek approaches that reduce the computational time needed for simulations, yet preserve the characteristic behavior of the system in question. In this work, we develop a computational framework based on the Itô stochastic differential equations for neuronal signal transduction networks. There are several different ways to incorporate stochasticity into deterministic differential equation models and to obtain Itô stochastic differential equations. Two of the developed models are found most suitable for stochastic modeling of neuronal signal transduction. The best models give stable responses which means that the variances of the responses with time are not increasing and negative concentrations are avoided. We also make a comparative analysis of different kinds of stochastic approaches, that is the Itô stochastic differential equations, the chemical Langevin equation, and the Gillespie stochastic simulation algorithm. Different kinds of stochastic approaches can be used to produce similar responses for the neuronal protein kinase C signal transduction pathway. The fine details of the responses vary slightly, depending on the approach and the parameter values. However, when simulating great numbers of chemical species, the Gillespie algorithm is computationally several orders of magnitude slower than the Itô stochastic differential equations and the chemical Langevin equation. Furthermore, the chemical Langevin equation produces negative concentrations. The Itô stochastic differential equations developed in this work are shown to overcome the problem of obtaining negative values.

Evolution of complex probability distributions in enzyme cascades

Unusual probability distribution profiles, including transient multi-peak distributions, have been observed in computer simulations of cell signaling dynamics. The emergence of these complex distributions cannot be explained using either deterministic chemical kinetics or simple Gaussian noise approximation. To develop physical insights into the origin of complex distributions in stochastic cell signaling, we compared our approximate analytical solutions of signaling dynamics with the exact numerical simulations. Our results are based on studying signaling in 2-step and 3-step enzyme amplification cascades that are among the most common building blocks of cellular protein signaling networks. We have found that while the multi-peak distributions are typically transient, and eventually evolve into single peak distributions, in certain cases these distributions may be stable in the limit of long times. We also have shown that introducing positive feedback loops results in diminution of the probability distribution complexity.

Phantom bursting is highly sensitive to noise and unlikely to account for slow bursting in beta-cells: considerations in favor of metabolically driven oscillations

Pancreatic beta-cells show bursting electrical activity with a wide range of burst periods ranging from a few seconds, often seen in isolated cells, over tens of seconds (medium bursting), usually observed in intact islets, to several minutes. The phantom burster model [Bertram, R., Previte, J., Sherman, A., Kinard, T.A., Satin, L.S., 2000. The phantom burster model for pancreatic beta-cells. Biophys. J. 79, 2880-2892] provided a framework, which covered this span, and gave an explanation of how to obtain medium bursting combining two processes operating on different time scales. However, single cells are subjected to stochastic fluctuations in plasma membrane currents, which are likely to disturb the bursting mechanism and transform medium bursters into spikers or very fast bursters. We present a polynomial, minimal, phantom burster model and show that noise modifies the plateau fraction and lowers the burst period dramatically in phantom bursters. It is therefore unlikely that slow bursting in single cells is driven by the slow phantom bursting mechanism, but could instead be driven by oscillations in glycolysis, which we show are stable to random ion channel fluctuations. Moreover, so-called compound bursting can be converted to apparent slow bursting by noise, which could explain why compound bursting and mixed Ca(2+) oscillations are seen mainly in intact islets.

Studying genetic regulatory networks at the molecular level: Delayed reaction stochastic models

Current advances in molecular biology enable us to access the rapidly increasing body of genetic information. It is still challenging to model gene systems at the molecular level. Here, we propose two types of reaction kinetic models for constructing genetic networks. Time delays involved in transcription and translation are explicitly considered to explore the effects of delays, which may be significant in genetic networks featured with feedback loops. One type of model is based on delayed effective reactions, each reaction modeling a biochemical process like transcription without involving intermediate reactions. The other is based on delayed virtual reactions, each reaction being converted from a mathematical function to model a biochemical function like gene inhibition. The latter stochastic models are derived from the corresponding mean-field models. The former ones are composed of single gene expression modules. We thus design a model of gene expression. This model is verified by our simulations using a delayed stochastic simulation algorithm, which accurately reproduces the stochastic kinetics in a recent experimental study. Various simplified versions of the model are given and evaluated. We then use the two methods to study the genetic toggle switch and the repressilator. We define the "on" and "off" states of genes and extract a binary code from the stochastic time series. The binary code can be described by the corresponding Boolean network models in certain conditions. We discuss these conditions, suggesting a method to connect Boolean models, mean-field models, and stochastic chemical models.

In vivo modeling of beta-glucan degradation in contrasting barley (Hordeum vulgare L.) genotypes

An important determinative of malt quality is the malt beta-glucan content, which in turn depends on the initial barley beta-glucan content as well as the beta-glucan depolymerization by beta-glucanase (EC 3.2.1.73) during malting. Another enzyme, named beta-glucan solubilase, has been suggested to act prior to beta-glucanase; its existence, however, has not been unequivocally proven. We monitored changes in beta-glucan levels and in the development of beta-glucan-degrading enzymes during malting of five lots of contrasting barley genotypes. Two models of in vivo kinetics for beta-glucan degradation were then compared as follows: (i) a biphasic model based on the sequential action of beta-glucan solubilase and beta-glucanase and (ii) a monophasic model assuming that all beta-glucans are depolymerized by beta-glucanase without the previous intervention of another enzyme. Confirmatory regression analysis was used to test the fit of the models to the observed data. Our results show that beta-glucan degradation is mostly monophasic, although some enzyme other than beta-glucanase seems to be required for the early solubilization of a small fraction of insoluble beta-glucans (on average, 7% of total beta-glucans). Furthermore, the genotype-dependent kinetic rate constant (indicating beta-glucan degradability), in addition to beta-glucanase activity, is suggested to play a major role in malting quality.

Simulating the temporal modulation of inducible DNA damage response in Escherichia coli

Living organisms make great efforts to maintain their genetic information integrity. However, DNA is vulnerable to many chemical or physical agents. To rescue the cell timely and effectively, the DNA damage response system must be well controlled. Recently, single cell experiments showing that after DNA damage, expression of the key DNA damage response regulatory protein oscillates with time. This phenomenon is observed both in eukaryotic and bacterial cells. We establish a model to simulate the DNA damage response (SOS response) in bacterial cell Escherichia coli. The simulation results are compared to the experimental data. Our simulation results suggest that the modulation observed in the experiment is due to the fluctuation of inducing signal, which is coupled with DNA replication. The inducing signal increases when replication is blocked by DNA damage and decreases when replication resumes.

Sustained oscillations via coherence resonance in SIR

Sustained oscillations in a stochastic SIR model are studied using a new multiple scale analysis. It captures the interaction of the deterministic and stochastic elements together with the separation of time scales inherent in the appearance of these dynamics. The nearly regular fluctuations in the infected and susceptible populations are described via an explicit construction of a stochastic amplitude equation. The agreement between the power spectral densities of the full model and the approximation verifies that coherence resonance is driving the behavior. The validity criteria for this asymptotic approximation give explicit expressions for the parameter ranges in which one expects to observe this phenomenon.

Process noise: an explanation for the fluctuations in the immune response during acute viral infection

The parameters of the immune response dynamics are usually estimated by the use of deterministic ordinary differential equations that relate data trends to parameter values. Since the physical basis of the response is stochastic, we are investigating the intensity of the data fluctuations resulting from the intrinsic response stochasticity, the so-called process noise. Dealing with the CD8+ T-cell responses of virus-infected mice, we find that the process noise influence cannot be neglected and we propose a parameter estimation approach that includes the process noise stochastic fluctuations. We show that the variations in data can be explained completely by the process noise. This explanation is an alternative to the one resulting from standard modeling approaches which say that the difference among individual immune responses is the consequence of the difference in parameter values.

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.

Extracting biochemical parameters for cellular modeling: A mean-field approach

Recent developments in molecular biology have made it feasible to carry out experimental verification of mathematical models for biochemical processes, offering the eventual prospect of creating a detailed, validated picture of gene expression. A persistent difficulty with this long-term goal is the incompleteness of the kinetic information available in the literature: Many rate constants cannot or have not yet been measured. Here, we present a method of filling in missing parameters using an approach conceptually analogous to mean-field approaches in statistical mechanics: When studying a particular gene, we extract key parameters by considering the averaged effect of all other genes in the system, analogously to considering the averaged magnetic field in a physical spin model. This methodology has been applied to account for the effect of the presence of the Escherichia coli genome on the availability of key enzymes involved in gene expression (RNA polymerases and ribosomes), yielding the number of free enzymes as a function of cellular growth rate. These conclusions have been obtained by deriving genome-wide averages and matching them to bulk literature values of E. coli K-12 and B/r. Average rate constants have been found for RNA polymerases and ribosomes binding to promoter and ribosome-binding sites, respectively; these results suggest that cells vary not only their production rates of RNA polymerase and ribosomes under different growth-rate conditions but also change their global level of transcriptional/translational activation and repression, thus altering the average binding rate constants for these enzymes. To test the mean-field method, the results from the genome-wide averages have been applied to the induced lac operon, where our derived on-rate for binding of RNA polymerase to the promoter is in good agreement with previous experimental results.

Sources of anomalous diffusion on cell membranes: a Monte Carlo study

A stochastic random walk model of protein molecule diffusion on a cell membrane was used to investigate the fundamental causes of anomalous diffusion in two-dimensional biological media. Three different interactions were considered: collisions with fixed obstacles, picket fence posts, and capture by, or exclusion from, lipid rafts. If motion is impeded by randomly placed, fixed obstacles, we find that diffusion can be highly anomalous, in agreement with previous studies. In contrast, collision with picket fence posts has a negligible effect on the anomalous exponent at realistic picket fence parameters. The effects of lipid rafts are more complex. If proteins partition into lipid rafts there is a small to moderate effect on the anomalous exponent, whereas if proteins are excluded from rafts there is a large effect on the anomalous exponent. In combination, these mechanisms can explain the level of anomaly in experimentally observed membrane diffusion, suggesting that anomalous diffusion is caused by multiple mechanisms whose effects are approximately additive. Finally, we show that the long-range diffusion rate, D(macro), estimated from fluorescence recovery after photobleaching studies, can be much smaller than D(micro), the small-scale diffusion rate, and is highly sensitive to obstacle densities and other impeding structures.

The role of compensatory mutations in the emergence of drug resistance

Pathogens that evolve resistance to drugs usually have reduced fitness. However, mutations that largely compensate for this reduction in fitness often arise. We investigate how these compensatory mutations affect population-wide resistance emergence as a function of drug treatment. Using a model of gonorrhea transmission dynamics, we obtain generally applicable, qualitative results that show how compensatory mutations lead to more likely and faster resistance emergence. We further show that resistance emergence depends on the level of drug use in a strongly nonlinear fashion. We also discuss what data need to be obtained to allow future quantitative predictions of resistance emergence.

Quasi-steady-state kinetics at enzyme and substrate concentrations in excess of the Michaelis-Menten constant

In vitro enzyme reactions are traditionally conducted under conditions of pronounced substrate excess since this guarantees that the bound enzyme is at quasi-steady-state (QSS) with respect to the free substrate, thereby justifying the Briggs-Haldane approximation (BHA). In contrast, intracellular reactions, amplification assays, allergen digestion assays and industrial applications span a range of enzyme-to-substrate ratios for which the BHA is invalid, including the extreme of enzyme excess. The quasi-equilibrium approximation (QEA) is valid for a subset of enzyme excess states. Previously, we showed that the total QSSA (tQSSA) overlaps and extends the validity of the BHA and the QEA, and that it is at least roughly valid for any total substrate and enzyme concentrations. The analysis of the tQSSA is hampered by square root nonlinearity. Previous simplifications of the tQSSA rate law are valid in a parameter domain that overlaps the validity domains of the BHA and the QEA and only slightly extends them. We now integrate the tQSSA rate equation in closed form, without resorting to further approximations. Moreover, we introduce a complimentary simplification of the tQSSA rate law that is valid in states of enzyme excess when the absolute difference between total enzyme and substrate concentrations greatly exceeds the Michaelis-Menten constant. This includes a wide range of enzyme and substrate concentrations where both the BHA and the QEA are invalid and allows us to define precisely the conditions for zero-order and first-order product formation. Remarkably, analytical approximations provided by the tQSSA closely match the expected stochastic kinetics for as few as 15 reactant molecules, suggesting that the conditions for the validity of the tQSSA and for its various simplifications are also of relevance at low molecule numbers.

Steady-state kinetic behaviour of functioning-dependent structures

A fundamental problem in biochemistry is that of the nature of the coordination between and within metabolic and signalling pathways. It is conceivable that this coordination might be assured by what we term functioning-dependent structures (FDSs), namely those assemblies of proteins that associate with one another when performing tasks and that disassociate when no longer performing them. To investigate a role in coordination for FDSs, we have studied numerically the steady-state kinetics of a model system of two sequential monomeric enzymes, E(1) and E(2). Our calculations show that such FDSs can display kinetic properties that the individual enzymes cannot. These include the full range of basic input/output characteristics found in electronic circuits such as linearity, invariance, pulsing and switching. Hence, FDSs can generate kinetics that might regulate and coordinate metabolism and signalling. Finally, we suggest that the occurrence of terms representative of the assembly and disassembly of FDSs in the classical expression of the density of entropy production are characteristic of living systems.

Parameter estimation in stochastic biochemical reactions

Gene regulatory, signal transduction and metabolic networks are major areas of interest in the newly emerging field of systems biology. In living cells, stochastic dynamics play an important role; however, the kinetic parameters of biochemical reactions necessary for modelling these processes are often not accessible directly through experiments. The problem of estimating stochastic reaction constants from molecule count data measured, with error, at discrete time points is considered. For modelling the system, a hidden Markov process is used, where the hidden states are the true molecule counts, and the transitions between those states correspond to reaction events following collisions of molecules. Two different algorithms are proposed for estimating the unknown model parameters. The first is an approximate maximum likelihood method that gives good estimates of the reaction parameters in systems with few possible reactions in each sampling interval. The second algorithm, treating the data as exact measurements, approximates the number of reactions in each sampling interval by solving a simple linear equation. Maximising the likelihood based on these approximations can provide good results, even in complex reaction systems.

A systematic investigation of the rate laws valid in intracellular environments

Recently there has been significant interest in deducing the form of the rate laws for chemical reactions occurring in the intracellular environment. This environment is typically characterized by low-dimensionality and a high macromolecular content; this leads to a spatial heterogeneity not typical of the well stirred in vitro environments. For this reason, the classical law of mass action has been presumed to be invalid for modeling intracellular reactions. Using lattice-gas automata models, it has recently been postulated [H. Berry, Monte Carlo simulations of enzyme reactions in two dimensions: Fractal kinetics and spatial segregation, Biophys. J. 83 (2002) 1891-1901; S. Schnell, T.E. Turner, Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws, Prog. Biophys. Mol. Biol. 85 (2004) 235-260] that the reaction kinetics is fractal-like. In this article we systematically investigate for the first time how the rate laws describing intracellular reactions vary as a function of: the geometry and size of the intracellular surface on which the reactions occur, the mobility of the macromolecules responsible for the crowding effects, the initial reactant concentrations and the probability of reaction between two reactant molecules. We also compare the rate laws valid in heterogeneous environments in which there is an underlying spatial lattice, for example crystalline alloys, with the rate laws valid in heterogeneous environments where there is no such natural lattice, for example in intracellular environments. Our simulations indicate that: (i) in intracellular environments both fractal kinetics and mass action can be valid, the major determinant being the probability of reaction, (ii) the geometry and size of the intracellular surface on which reactions are occurring does not significantly affect the rate law, (iii) there are considerable differences between the rate laws valid in heterogeneous non-living structures such as crystals and those valid in intracellular environments. Deviations from mass action are less pronounced in intracellular environments than in a crystalline material of similar heterogeneity.

Effects of quenched and annealed macromolecular crowding elements on a simple model for signaling in T lymphocytes

Biochemical reactions in cells occur in an environment that is crowded in the sense that various macromolecular species and organelles occupy much of the space. The effects of molecular crowding on biochemical reactions have usually been studied in the past in a spatially homogeneous environment. However, signal transduction in cells is often initiated by the binding of receptors and ligands in two apposed cell membranes, and the pertinent biochemical reactions occur in a spatially inhomogeneous environment. We have studied the effects of crowding on biochemical reactions that involve both membrane proteins and cytosolic molecules by investigating a simplified version of signaling in T lymphocytes using a Monte Carlo algorithm. We find that, if signal transduction occurs on time scales that are slow compared to the motility of the molecules and organelles that constitute the crowding elements, the effects of crowding are qualitatively the same as in a homogeneous three-dimensional (3D) medium. In contrast, if signal transduction occurs on a time scale that is much faster than the time over which the crowding elements move, then the effects of varying the extent of crowding are very different when reactions occur in both 2- and 3D space. We discuss these differences and their origin. Since many signaling reactions are fast, our results may be useful for diverse situations in cell biology.

Accommodating space, time and randomness in network simulation

Interest in the possibility of dynamically simulating complex cellular processes has escalated markedly in recent years. This interest has been fuelled by three factors: the generally accepted value in understanding living processes as integrated systems; the dramatic increase in computational capability; and the availability of new or improved technology for making the quantitative measurements that are needed to drive and validate cellular simulations. Between the extremes of atom-scale and organism-scale simulation is a vast middle-ground requiring simulation strategies that are capable of dealing with a range of spatial, temporal and molecular abundance scales that are crucial for a comprehensive understanding of integrative cell biology. Although at an early stage, methodological improvements and the development of computational platforms provide some hope that simulations will emerge that can bridge the gap between network models and the true operation of the cell as a complex machine.

Solving the chemical master equation for monomolecular reaction systems analytically

The stochastic dynamics of a well-stirred mixture of molecular species interacting through different biochemical reactions can be accurately modelled by the chemical master equation (CME). Research in the biology and scientific computing community has concentrated mostly on the development of numerical techniques to approximate the solution of the CME via many realizations of the associated Markov jump process. The domain of exact and/or efficient methods for directly solving the CME is still widely open, which is due to its large dimension that grows exponentially with the number of molecular species involved. In this article, we present an exact solution formula of the CME for arbitrary initial conditions in the case where the underlying system is governed by monomolecular reactions. The solution can be expressed in terms of the convolution of multinomial and product Poisson distributions with time-dependent parameters evolving according to the traditional reaction-rate equations. This very structured representation allows to deduce easily many properties of the solution. The model class includes many interesting examples. For more complex reaction systems, our results can be seen as a first step towards the construction of new numerical integrators, because solutions to the monomolecular case provide promising ansatz functions for Galerkin-type methods.

The Total Quasi-Steady-State Approximation for Fully Competitive Enzyme Reactions

The validity of the Michaelis-Menten-Briggs-Haldane approximation for single enzyme reactions has recently been improved by the formalism of the total quasi-steady-state approximation. This approach is here extended to fully competitive systems, and a criterion for its validity is provided. We show that it extends the Michaelis-Menten-Briggs-Haldane approximation for such systems for a wide range of parameters very convincingly, and investigate special cases. It is demonstrated that our method is at least roughly valid in the case of identical affinities. The results presented should be useful for numerical simulations of many in vivo reactions.

An optimal number of molecules for signal amplification and discrimination in a chemical cascade

Understanding the information processing ability of signal transduction pathways is of great importance because of their crucial roles in triggering various cellular responses. Despite continuing theoretical investigation, some important aspects of signal transduction such as a transient response and its connection to stochasticity originating from a small number of molecules have not yet been well understood. It is, however, through these aspects that unexpected and nontrivial properties of the information processing emerge. In this article, we analyze the transient behavior of a simple signaling cascade by taking into account the stochasticity originating from the small number of molecules. We identify several properties of the signaling cascade that emerge as a result of the interplay between the stochasticity and transient dynamics of the cascade. We specifically demonstrate that each step of the cascade has an optimal number of signaling molecules at which the average signal amplitude becomes maximal. We further investigate the connection between a finite number of molecules and the ability of the cascade to discriminate between true and error signals, which cannot be inferred from deterministic descriptions. The implications of our results are discussed from both biological and mathematical viewpoints.

Multifractality in intracellular enzymatic reactions

Enzymatic kinetics adjust well to the Michaelis-Menten paradigm in homogeneous media with dilute, perfectly mixed reactants. These conditions are quite different from the highly structured cell plasm, so applications of the classic kinetics theory to this environment are rather limited. Cytoplasmic structure produces molecular crowding and anomalous diffusion of substances, modifying the mass action kinetic laws. The reaction coefficients are no longer constant but time-variant, as stated in the fractal kinetics theory. Fractal kinetics assumes that enzymatic reactions on such heterogeneous media occur within a non-Euclidian space characterized by a certain fractal dimension, this fractal dimension gives the dependence on time of the kinetic coefficients. In this work, stochastic simulations of enzymatic reactions under molecular crowding have been completed, and kinetic coefficients for the reactions, including the Michaelis-Menten parameter KM, were calculated. The simulations results led us to confirm the time dependence of michaelian kinetic parameter for the enzymatic catalysis. Besides, other chaos related phenomena were pointed out from the obtained KM time series, such as the emergence of strange attractors and multifractality.

The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior

A key to advancing the understanding of molecular biology in the post-genomic age is the development of accurate predictive models for genetic regulation, protein interaction, metabolism, and other biochemical processes. To facilitate model development, simulation algorithms must provide an accurate representation of the system, while performing the simulation in a reasonable amount of time. Gillespie's stochastic simulation algorithm (SSA) accurately depicts spatially homogeneous models with small populations of chemical species and properly represents noise, but it is often abandoned when modeling larger systems because of its computational complexity. In this work, we examine the performance of different versions of the SSA when applied to several biochemical models. Through our analysis, we discover that transient changes in reaction execution frequencies, which are typical of biochemical models with gene induction and repression, can dramatically affect simulator performance. To account for these shifts, we propose a new algorithm called the sorting direct method that maintains a loosely sorted order of the reactions as the simulation executes. Our measurements show that the sorting direct method performs favorably when compared to other well-known exact stochastic simulation algorithms.

Stochastic kinetics of viral capsid assembly based on detailed protein structures

We present a generic computational framework for the simulation of viral capsid assembly which is quantitative and specific. Starting from PDB files containing atomic coordinates, the algorithm builds a coarse-grained description of protein oligomers based on graph rigidity. These reduced protein descriptions are used in an extended Gillespie algorithm to investigate the stochastic kinetics of the assembly process. The association rates are obtained from a diffusive Smoluchowski equation for rapid coagulation, modified to account for water shielding and protein structure. The dissociation rates are derived by interpreting the splitting of oligomers as a process of graph partitioning akin to the escape from a multidimensional well. This modular framework is quantitative yet computationally tractable, with a small number of physically motivated parameters. The methodology is illustrated using two different viruses which are shown to follow quantitatively different assembly pathways. We also show how in this model the quasi-stationary kinetics of assembly can be described as a Markovian cascading process, in which only a few intermediates and a small proportion of pathways are present. The observed pathways and intermediates can be related a posteriori to structural and energetic properties of the capsid oligomers.

Theoretical analysis of internal fluctuations and bistability in CO oxidation on nanoscale surfaces

The bistable CO oxidation on a nanoscale surface is characterized by a limited number of reacting molecules on the catalytic area. Internal fluctuations due to finite-size effects are studied by the master equation with a Langmuir-Hinshelwood mechanism for CO oxidation. Analytical solutions can be found in a reduced one-component model after the adiabatic elimination of one variable which in our case is the oxygen coverage. It is shown that near the critical point, with decreasing surface area, one cannot distinguish between two macroscopically stable stationary states. This is a consequence of the large fluctuations in the coverage which occur on a fast time scale. Under these conditions, the transition times between the macroscopic states also are no longer separated from the short-time scale of the coverage fluctuations as is the case for large surface areas and far away from the critical point. The corresponding stationary solutions of the probability distribution and the mean first passage times calculated in the reduced model are supported by numerics of the full two-component model.

About implementing a Monte Carlo simulation algorithm for enzymatic reactions in crowded media

In this paper, several aspects of implementing a Monte Carlo simulation algorithm for studying the Michaelis-Menten mechanism of enzymatic reactions in crowded media are presented. Using a two dimensional lattice with obstacles, it is shown how the initial distribution of the reactants and obstacles on the lattice affects the values of the rate coefficients and the concentration of the reactants. The influence of the number of considered nearest neighbors and of the obstacle concentration on the values of the rate coefficients is also demonstrated. The results strongly suggest fractal kinetics for enzyme reactions in crowded media.

Using Petri Net tools to study properties and dynamics of biological systems

Petri Nets (PNs) and their extensions are promising methods for modeling and simulating biological systems. We surveyed PN formalisms and tools and compared them based on their mathematical capabilities as well as by their appropriateness to represent typical biological processes. We measured the ability of these tools to model specific features of biological systems and answer a set of biological questions that we defined. We found that different tools are required to provide all capabilities that we assessed. We created software to translate a generic PN model into most of the formalisms and tools discussed. We have also made available three models and suggest that a library of such models would catalyze progress in qualitative modeling via PNs. Development and wide adoption of common formats would enable researchers to share models and use different tools to analyze them without the need to convert to proprietary formats.