Is this adaption of the Gillespie algorithm using Michaelis constants justifiable?

Is this adaption of the Gillespie algorithm using Michaelis constants justifiable?

We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

I would like to run a discrete simulation of a biological system as it can be done, e.g., using the Gillespie algorithm. However, the Gillespie algorithm requires you to know the reaction rate constants of each involved reaction whereas I only know the Michaelis constant $K_M$ and the maximum rate $V_{max}$ of each reaction.

Since the Gillespie algorithm only uses the reaction rate constants to calculate the reaction rate $v$ (which, again, is used to determine the probability of a given reaction), I wonder wheter it is justifiable to simply replace this calculation of $v = k * [S]$ by $v = frac{V_{max}cdot[S]}{K_M + [S]}$ for each reaction $S ightarrow P$ and, apart from that, use the Gillespie algorithm as is? If not, is there any better way to run a discrete simulation of the system with the given parameters?

Note: I've already asked in a different question whether it is possible to get (an estimation of) the reaction rates knowing only the parameters $K_M$ and $V_{max}$. I hope my both questions are sufficiently distinct to justify two different posts.

The stochastic simulation algorithm by Gillespie makes the assumption that reaction propensity (probability) in an infinitesimally small interval of time is same as the macroscopic reaction rate. (For details, you can see my answer in Chemistry.SE).

You can model and simulate every step of enzyme catalysis or you can also use the Michaelis-Menten model. In any case, you can extend the same approximations to the stochastic model.

This is also applicable to other Michaelis-Menten like (or Hill) functions of the form $$frac{x^n}{k^n+x^n} qquad ext{or}qquad frac{k^n}{k^n+x^n}$$.

There are many papers that have done this. Note that the non-linear nature of these functions results in deviation of the stochastic behaviour from the deterministic behaviour. Read about stochastic focussing. Modelling each step exactly (binding, unbinding and catalysis) may result in different behaviour depending on your parameters. Whether you should use Michaelis-Menten function to describe reaction rate totally depends on how reasonable your approximation is.

Very interesting question. People use Gillespie's algorithm for several sort of kinetic functions such as Michaelis-Menten (MM), Hill etc., no matter what the original assumptions of Gillespie's algorithm are.

Concerning MM, There was a 2011 paper by Gillespie himself that shows that the approximation is applicable in discrete stochastic models and that the validity conditions are the same as in the deterministic regime.

Remember that, because MM is derived under a QSSA, you can break down the S->P transformation as its single-step components. See for instance a comparison here

Derivation of stationary distributions of biochemical reaction networks via structure transformation

Long-term behaviors of biochemical reaction networks (BRNs) are described by steady states in deterministic models and stationary distributions in stochastic models. Unlike deterministic steady states, stationary distributions capturing inherent fluctuations of reactions are extremely difficult to derive analytically due to the curse of dimensionality. Here, we develop a method to derive analytic stationary distributions from deterministic steady states by transforming BRNs to have a special dynamic property, called complex balancing. Specifically, we merge nodes and edges of BRNs to match in- and out-flows of each node. This allows us to derive the stationary distributions of a large class of BRNs, including autophosphorylation networks of EGFR, PAK1, and Aurora B kinase and a genetic toggle switch. This reveals the unique properties of their stochastic dynamics such as robustness, sensitivity, and multi-modality. Importantly, we provide a user-friendly computational package, CASTANET, that automatically derives symbolic expressions of the stationary distributions of BRNs to understand their long-term stochasticity.


Motivation: Simulation and modeling is becoming a standard approach to understand complex biochemical processes. Therefore, there is a big need for software tools that allow access to diverse simulation and modeling methods as well as support for the usage of these methods.

Results: Here, we present COPASI, a platform-independent and user-friendly biochemical simulator that offers several unique features. We discuss numerical issues with these features in particular, the criteria to switch between stochastic and deterministic simulation methods, hybrid deterministic–stochastic methods, and the importance of random number generator numerical resolution in stochastic simulation.

Availability: The complete software is available in binary (executable) for MS Windows, OS X, Linux (Intel) and Sun Solaris (SPARC), as well as the full source code under an open source license from Author Webpage.

1. Introduction

Non-genetic heterogeneity is a hallmark of cell physiology. Isogenic cells can display markedly different phenotypes as a result of the stochasticity of intracellular processes and fluctuations in environmental conditions. Gene expression variability, in particular, has received substantial attention thanks to robust experimental techniques for measuring transcripts and proteins at a single-cell resolution (Golding et al., 2005 Taniguchi et al., 2010). This progress has gone hand-in-hand with a large body of theoretical work on stochastic models to identify the molecular processes that affect expression heterogeneity (Swain et al., 2002 Raj and van Oudenaarden, 2008 Thomas et al., 2014 Dattani and Barahona, 2017 Tonn et al., 2019).

In contrast to gene expression, our understanding of stochastic phenomena in metabolism is still in its infancy. Traditionally, cellular metabolism has been regarded as a deterministic process on the basis that metabolites appear in large numbers that filter out stochastic phenomena (Heinemann and Zenobi, 2011). But this view is changing rapidly thanks to a growing number of single-cell measurements of metabolites and co-factors (Bennett et al., 2009 Imamura et al., 2009 Lemke and Schultz, 2011 Paige et al., 2012 Ibá༞z et al., 2013 Yaginuma et al., 2014 Esaki and Masujima, 2015 Xiao et al., 2016 Mannan et al., 2017) that suggest that cell-to-cell metabolite variation is much more pervasive than previously thought. The functional implications of this heterogeneity are largely unknown but likely to be substantial given the roles of metabolism in many cellular processes, including growth (Weisse et al., 2015), gene regulation (Lempp et al., 2019), epigenetic control (Loftus and Finlay, 2016), and immunity (Reid et al., 2017). For example, metabolic heterogeneity has been linked to bacterial persistence (Radzikowski et al., 2017 Shan et al., 2017), a dormant phenotype characterized by a low metabolic activity, as well as antibiotic resistance (Deris et al., 2013) and other functional effects (Vilhena et al., 2018). In biotechnology applications, metabolic heterogeneity is widely recognized as a limiting factor on metabolite production with genetically engineered microbes (Binder et al., 2017 Schmitz et al., 2017 Liu et al., 2018).

A key challenge for quantifying metabolic variability is the difficulty in measuring cellular metabolites at a single-cell resolution (Amantonico et al., 2010 Takhaveev and Heinemann, 2018 Wehrens et al., 2018). As a result, most studies use other phenotypes as a proxy for metabolic variation, e.g., enzyme expression levels (Kotte et al., 2014 van Heerden et al., 2014), metabolic fluxes (Schreiber et al., 2016), or growth rate (Kiviet et al., 2014 Şimᗾk and Kim, 2018). From a computational viewpoint, the key challenge is that metabolic processes operate on two timescales: a slow timescale for expression of metabolic enzymes, and a fast timescale for enzyme catalysis. Such multiscale structure results in stiff models that are infeasible to solve with standard algorithms for stochastic simulation (Gillespie, 2007). Other strategies to accelerate stochastic simulations, such as τ-leaping (Rathinam et al., 2003), also fail to produce accurate simulation results due to the disparity in molecule numbers between enzymes and metabolites (Tonn, 2020). These challenges have motivated a number of methods to optimize stochastic simulations of metabolism (Puchałka and Kierzek, 2004 Cao et al., 2005 Labhsetwar et al., 2013 Lugagne et al., 2013 Murabito et al., 2014). Most of these methods exploit the timescale separation to accelerate simulations at the expense of some approximation error. This progress has been accompanied by a number of theoretical results on the links between molecular processes and the shape of metabolite distributions (Levine and Hwa, 2007 Oyarzún et al., 2015 Gupta et al., 2017b Tonn et al., 2019). Yet to date there are no general methods for computing metabolite distributions that can handle inherent features of metabolic pathways such as feedback regulation, complex stoichiometries, and the high number of molecular species involved.

In this paper we present a widely applicable method for approximating single-cell metabolite distributions. Our method is founded on the timescale separation between enzyme expression and enzyme catalysis, which we employ to approximate the stationary solution of the chemical master equation. The approximate solution takes the form of mixture distributions with: (i) mixture weights that can be computed from models for gene expression or single-cell expression data, and (ii) mixture components that are directly computable from deterministic pathway models. The resulting mixture model can be employed to explore the impact of biochemical parameters on metabolite variability. We illustrate the power of the method in two exemplar systems that are core building blocks of large metabolic networks. Our theory provides a quantitative basis to draw testable hypotheses on the sources of metabolite heterogeneity, which together with the ongoing efforts in single-cell metabolite measurements, will help to re-evaluate the role of metabolism as an active source of phenotypic variation.

Realizing ‘integral control’ in living cells: how to overcome leaky integration due to dilution?

A major problem in the design of synthetic genetic circuits is robustness to perturbations and uncertainty. Because of this, there have been significant efforts in recent years in finding approaches to implement integral control in genetic circuits. Integral controllers have the unique ability to make the output of a process adapt perfectly to disturbances. However, implementing an integral controller is challenging in living cells. This is because a key aspect of any integral controller is a ‘memory’ element that stores the accumulation (integral) of the error between the output and its desired set-point. The ability to realize such a memory element in living cells is fundamentally challenged by the fact that all biomolecules dilute as cells grow, resulting in a ‘leaky’ memory that gradually fades away. As a consequence, the adaptation property is lost. Here, we propose a general principle for designing integral controllers such that the performance is practically unaffected by dilution. In particular, we mathematically prove that if the reactions implementing the integral controller are all much faster than dilution, then the adaptation error due to integration leakiness becomes negligible. We exemplify this design principle with two synthetic genetic circuits aimed at reaching adaptation of gene expression to fluctuations in cellular resources. Our results provide concrete guidance on the biomolecular processes that are most appropriate for implementing integral controllers in living cells.

1. Introduction

A long-standing challenge in synthetic biology is the creation of genetic circuits that perform predictably and reliably in living cells [1,2]. Variations in the environment [3], unforeseen interactions between a circuit and the host machinery [4], and unknown interference between circuit components [5] can all significantly alter the intended behaviour of a genetic circuit. Integral control, which is ubiquitous in traditional engineering systems, is a promising approach to mitigate the impact of such unknowns on genetic circuits. At the most basic level, an integral controller ensures that the output of a system can reach a desired set-point while perfectly rejecting (i.e. adapting to) constant disturbances [6–8]. Owing to this unique ability, synthetic biology research has witnessed increasing efforts towards establishing procedures to realize integral control in living cells [9–17].

At the core of an integral controller is an indispensable ‘memory’ element that accumulates (i.e. integrates) the error between the desired set-point and the measured output over time. While computing this integral is almost trivial using electronic components, realizing it through biomolecular reactions in living cells is problematic. This is because information about the error is usually stored in the form of molecular concentrations, and as host cell grows and divides, biomolecules in every single cell dilute [18,19], leading to a ‘leaky’ memory that fades away in time. As cell growth is unavoidable and even beneficial in a number of applications, such as in biofuel production [20], leaky integration is a fundamental physical limitation for the implementation of integral control through genetic circuits in living cells.

Prior theoretical studies have proposed integral control motifs (ICMs), but neglected the presence of molecule dilution in the analysis [10,12,15,17]. Consequently, while these motifs have theoretically guaranteed performance in cell-free systems, where dilution is non-existent, their performance in the presence of molecule dilution is not preserved. In fact, the deteriorating effect of leaky integration can be significant, as shown in this paper and in other recent studies (e.g. [21,22]). Therefore, possible solutions to leaky integration have been proposed [10,21]. Specifically, the approach proposed in [10] hinges on ‘cancelling’ the effect of dilution by producing the ‘memory species’ at the same rate as its dilution. This approach, however, relies on exact parameter matching and is hence difficult to realize in practice. A different approach is to use exhaustive numerical simulation to extract circuit parameters that minimize the effects of leaky integration [21], leading to a lengthy and ad hoc design process.

In this report, we propose general design principles for ICMs that perform well despite the presence of dilution. In particular, we mathematically demonstrate that the undesirable effect of dilution can be arbitrarily suppressed by increasing the rate of all controller reactions. This is true under mild conditions that are largely independent of specific parameter values. This design principle guides the choice of core biomolecular processes and circuit parameters that are most suitable for realizing an ICM in living cells. We illustrate this guided choice on two circuits that are designed to mitigate the effects of transcription/translation resource fluctuations on gene expression, a problem that has gained significant attention recently [4,23–25].

2. Quasi-integral control

The approach that we take in this paper is as follows. We first describe two types of ideal ICMs in §2.1, which were previously proposed abstract circuit motifs for adapting to constant disturbances in the absence of dilution. We then introduce leaky ICMs, which add dilution to these ideal ICMs, and demonstrate that the adaptation property is lost. Finally, we describe quasi-ICMs in §2.2, the main novelty of this paper, in which the controller reactions are engineered to be much faster than dilution, enabling the circuit to restore almost perfect adaptation to constant disturbances.

2.1. Ideal integral control motifs and leaky integral control motifs

We illustrate in figure 1a two different types of ideal ICMs that abstract the two main mechanisms for biomolecular integral control proposed in the literature. In both types of motifs, we denote by x the species whose concentration needs to be kept at a set-point u, despite that its production rate is affected by a constant disturbance d. Therefore, the output of the process to be controlled is the concentration of species x, which is denoted by x (italic). The purpose of implementing an integral controller is thus to have the steady state of x, denoted by , to satisfy regardless of d (i.e. x adapts perfectly to d).

Figure 1. Quasi-integral control mitigates the effect of leaky integration due to dilution. (a) Two types of ideal ICMs. The controller reactions are boxed in pink, and the rest of the circuit belongs to the process to be controlled. Dilution of the controller species is neglected in ideal ICMs. (b) When dilution of the controller species is considered, ideal ICMs become leaky ICMs that cannot carry out integration. The adaptation error, defined as the extent to which the output is affected by the disturbance, can be arbitrarily large. (c) All controller reactions in a quasi-ICM are 1/ε times faster than those in the corresponding leaky ICM, with ε ≪ 1. (d) Simulations of the type I and type II ideal, leaky and quasi-ICMs. Simulation parameters for both motifs: α = γ = δ = k = 1 h −1 , θ = 1 nM −1 h −1 , ε = 0.02, u = 10 nM and d = 5 nM h −1 . Disturbance input d is applied at 15 h. (e) A general ε-quasi-integral control system. Output of the process y becomes an input to the quasi-integral controller. Variable z2 is the leaky memory variable, and z1 represents the remaining controller states, if any. All controller reactions are much faster than dilution, as characterized by parameter ε. Output of the controller v drives the process to set-point u and adapts to d. (Online version in colour.)

In these ideal ICMs, selected controller species realize the memory function essential for integral control. Specifically, a type I ideal ICM (figure 1a(i)) regulates the expression of x using a single controller species z. This motif arises from saturating certain Hill-type or Michaelis–Menten-type kinetics such that the production and the removal rates of z become roughly proportional to u and x, respectively, resulting in the following approximate mass action kinetics model for the controller species z [10,15,17]:

With reference to figure 1a(ii), a type II ideal ICM, arising from what was called the antithetic integral controller in [12], realizes integral action using two controller species z1 and z2. Their production rates are engineered to be proportional to the set-point (u) and to the concentration of the output species (x). The two controller species can bind together to form a complex that is removed with rate constant θ: . These biomolecular processes can be described by the following mass action kinetics model:

When dilution of the controller species due to host cell growth is taken into account, the key integral structure of the memory variables is disrupted (figure 1b). In fact, in a standard mass action kinetics model describing reactions in exponentially growing cells, the average effect of cell growth and division on the dynamics of any species z can be modelled by a first-order additive term −γz, where γ is the specific growth rate of the host cell [18,19,26,27] (see the electronic supplementary material, section S1.1, for details). Hence, with reference to figure 1b(i), the dynamics of the memory variable z in a type I ICM become

2.2. Quasi-integral control motifs: almost perfect adaptation to disturbances

In this section, we introduce quasi-ICMs, the main contribution of this paper. While leaky ICMs cannot achieve perfect adaptation to disturbances, and their performance can be arbitrarily deteriorated by the presence of dilution, we propose that one can achieve almost perfect adaptation to disturbances despite dilution by increasing the rate of all controller reactions in the leaky ICMs. Therefore, we call a leaky ICM whose controller reactions are much faster than the slower dilution process a quasi-ICM. These motifs are shown in figure 1c(i)(ii) for the two types of motifs. In particular, with reference to figure 1c(i), we use a small dimensionless positive parameter 0 < ε≪1 to capture the fact that the controller reaction rates in a type I quasi-ICM are 1/ε times faster than those in a type I leaky ICM (figure 1b(i)), resulting in the memory variable dynamics to become

We consider a biomolecular process to be controlled, whose dynamics can be written as

Definition 2.1

The closed loop system (2.9)–(2.10) realizes ε-quasi-integral control in an admissible input set if for all , the system's steady-state output satisfies

Equation (2.11) implies that, by increasing the rate of controller reactions (i.e. decreasing ε), the steady-state output of an ε-quasi-integral control system can be made arbitrarily close to the set-point u in the presence of disturbance d. The following claim provides an easy-to-check condition for ε-quasi-integral control.

Claim 2.2.

The closed loop system (2.9)–(2.10) realizes ε-quasi-integral control if the auxiliary ideal integral control system consisting of the process (2.9) and the ideal integral controller

The proof of this claim can be found in the electronic supplementary material, section S1.2. It establishes that if one can use the ideal integral controller (2.12) to obtain perfect adaptation in the absence of dilution, then when dilution is taken into account, one can always speed up all controller reactions to restore adaptation. From an engineering perspective, one should therefore select controller biomolecular processes whose rates are significantly faster with respect to dilution rates. Such processes include, for example, enzymatic reactions (i.e. phosphorylation, methylation and ubiquitination), protein–protein interaction, and RNA-based interactions (i.e. sRNA-or microRNA-enabled degradation or sequestration). With these processes, the designer can then realize type I or type II ideal ICMs as well as additional ideal ICMs that are in the form of system (2.12). It is worth to note that even in cases where only part of the controller reactions can be sped up, it may still be possible to reach ε-quasi-integral control, although the algebraic conditions to check become slightly more involved (see the electronic supplementary material, section S1.7). In the electronic supplementary material, section S4.3, we further demonstrate that when dilution rate constant varies slowly in time, the adaptation error is still guaranteed to decrease as ε is decreased.

Claim 2.2 provides the theoretical underpinning to the type I and type II quasi-ICM designs of figure 1c. As the corresponding type I and II ideal ICMs are in the form of (2.12) and have unique locally exponentially stable steady states (see proof in the electronic supplementary material, section S1.3–1.4), increasing all controller reaction rates will quench the steady-state effect of leaky integration. Furthermore, we also show in the electronic supplementary material, section S1.3–1.4, that the steady states of the type I and II quasi-ICMs are unique and stable.

3. Two biomolecular realizations of quasi-integral control motifs

In this section, we propose physical implementations of a phosphorylation-based quasi-integral controller and a small RNA-based quasi-integral controller. They are realizations of the type I and type II quasi-ICMs, respectively, which satisfy the general design principle of fast controller reactions. For illustration purposes, we use these controllers to mitigate the effect of transcriptional and translational resource fluctuations on gene expression, a problem that has received considerable attention in synthetic biology in recent years [4,23–25].

3.1. Type I quasi-integral control motif: phosphorylation cycle

A diagram of the phosphorylation-based quasi-integral control system is shown in figure 2a. The controller is intended to regulate the production of protein p to adapt to a disturbance d, which models a reduction in protein production rate due to, for example, the depletion of transcriptional and/or translational resources in the host cell [4,23]. In this system, the intended integral control is accomplished by a phosphorylation cycle. The regulated protein p is co-expressed with a phosphatase and the concentration of the kinase acts as the set-point u. The substrate b is expressed constitutively. When b is phosphorylated by kinase u to become active substrate b * , it transcriptionally activates production of protein p. A simplified mathematical model of this system is (see the electronic supplementary material, section S2.1, for detailed derivation from chemical reactions)

where k1 and k2 are the catalytic rate constants of the phosphorylation and dephosphorylation reactions, respectively, K1 and K2 are the Michaelis–Menten constants, R is the protein production rate constant, and λ is the dissociation constant between b * and the promoter of the regulated gene.

Figure 2. Two physical realizations of quasi-ICMs. (a) Genetic circuit diagram of the phosphorylation-based quasi-integral controller. Chemical reactions realizing the controller are boxed in pink. (b) Simulation of the circuit's response according to (3.1). A set-point input u = 20 nM is applied at time 0 and a disturbance input d = 0.5 is applied at 20 h. The vertical axis represents the ratio between output y, defined to be proportional to p (y = σp), and set-point u. The dashed black line is the response of the phosphorylation-based control system assuming no dilution of the active substrate b * . The dotted blue line, the thin green line with square markers and the solid red line represent circuit's response in the presence of nonzero substrate dilution (γ = 1 h −1 ) and decreasing ε, which corresponds to increasing catalytic rates (ki, i = 1, 2). (c) Genetic circuit diagram of the sRNA-based quasi-integral controller. (d) Simulation of the circuit's response according to (3.4). A set-point input u = 1 is applied at time 0 and a disturbance input d = 0.5 is applied at 30 h. The vertical axis represents the ratio between output y, defined in (3.6), and set-point input u. The dashed black line represents response of an ideal integral control system, where RNA decay rate δ is set to 0. The dotted blue line, the thin green line with square markers and the solid red line represent circuit's responses in the presence of nonzero RNA decay rate (δ = 3 h −1 , corresponding to half-life of approx. 13 min) and decreasing ε. Parameter ε is decreased by increasing the mRNA–sRNA removal rate (θ/β). The DNA copy numbers of the regulated gene and the sRNA are increased simultaneously by a factor of 1/ε as ε decreases. Simulation parameters are listed in electronic supplementary material, section S5, table S2. (Online version in colour.)

System (3.1) can be taken to the form of a type I leaky-ICM (2.3) if both the kinase and the phosphatase are saturated (i.e. bK1 and b * ≫K2). These design constraints are present in any type I ICM, which relies on saturation of Hill-type or Michaelis–Menten kinetics [10,15,17], and can be satisfied in practice when the substrate is over-expressed and the kinase concentration is not too small (see the electronic supplementary material, section S2.3, for details). Under these assumptions and by setting x: = p, z: = b * and σ: = k2/k1, we can approximate (3.1) by

3.2. Type II quasi-integral control motif: small RNA interference

The sRNA-based quasi-integral controller is a realization of the type II quasi-ICM. It is intended to regulate translation of protein p to adapt to a disturbance d, which models uncertainty in translation rate constant due to fluctuation in the amount of available ribosomes [4,23,24]. With reference to figure 2c, the controller consists of protein p transcriptionally activating production of an sRNA (s) that is complementary to the mRNA (m) of the regulated protein p. The sRNA and mRNA can bind and degrade together rapidly [28–30]. The mRNA concentration m is the control input to the translation process. A constant upstream transcription factor regulates mRNA production as a set-point input u to the controller. Based on the chemical reactions in the electronic supplementary material, section S3.1, the ODE model of this system is

Setting z1: = m, z2: = s and x: = p, system (3.4) can be taken to the form of a type II leaky-ICM:

4. Discussion

Owing to integral controllers' unique ability to drive a process to reach a set-point regardless of disturbances/uncertainties, they are especially appealing to synthetic biology, where various forms of disturbances and uncertainties limit the robustness and predictability of synthetic circuits. A unique and fundamental challenge to realize integral control in living cells is that the concentration of biomolecules dilutes as cells grow, resulting in the impossibility to create a non-leaky biomolecular memory element, which is indispensable to implement integral control.

In this report, we propose a general design principle to overcome this obstacle. We mathematically demonstrate that one can design an ideal integral controller to obtain perfect adaptation in the absence of dilution and then speed up all controller reactions with respect to dilution to restore adaptation in the presence of dilution. This implies that the designer should select controller biomolecular processes whose rates are significantly faster than dilution rates. Our result is equally applicable to situations where controller species are degraded in addition to being diluted. Owing to the asymptotic nature of our results, they are largely insensitive to system parameters. To solidify this point, we simulated the phosphorylation-based and the sRNA-based quasi-integral controllers with randomly generated parameters, and demonstrated that, regardless of the rest of the parameters, the adaptation error always decreases, as the rates of controller reactions are increased (electronic supplementary material, section S5.2).

While our mathematical analyses are carried out using ordinary differential equation models derived from mass action kinetics, we numerically exemplify their effectiveness using a time-varying Gillespie simulation [36], in which cell volume expansion and cell division as well as stochasticity in biomolecular reactions are explicitly taken into account (see the electronic supplementary material, section S5.3, for details). In addition, we have assumed that the host cells grow at constant rate (constant γ). In practice, however, their growth rate may change with, for example, temperature, nutrients and metabolic load of the circuit [37]. We have therefore studied the performance of quasi-integral controllers in changing growth conditions (see the electronic supplementary material, section S4.3), and found that the circuit's output can still adapt to disturbances as long as growth rate changes slowly.

Nature may have already used the time-scale separation strategy between process and controller that we propose here for realizing adaptation and homeostasis. For example, it has been suggested that the bacterial chemotaxis system [38] and the yeast osmoregulation system [39] contain integral controllers realized through methylation and phosphorylation, respectively. As these processes are much faster than dilution in their respective host cells, it is justifiable to neglect dilution in these systems, and almost perfect adaptation can be achieved. Therefore, experimental study of synthetic realizations of the quasi-integral controllers proposed here may enhance our understanding of their natural counterparts, which are often embedded in more complex biomolecular networks and are therefore more difficult to analyse.

From a classical control-theoretic perspective, integral controllers often operate with low integral gains (i.e. slow integration). This is because with a perfect integrator, steady-state adaptation can be achieved regardless of the magnitude of the integral gain, and a high integral gain often leads to undesirable effects such as higher energy consumption and tendency to instability [6,22]. However, a potential advantage of increasing the leaky integral gain is that it can improve a system's ability to track (reject) time-varying references (disturbances). While we demonstrate this capability for the two example systems through linearization and simulations in the electronic supplementary material, section S4.1–4.2, further theoretical assessment of this property requires deeper study [40,41].

Data accessibility

All simulation codes are provided as the electronic supplementary material.


In this paper we describe a Bayesian framework for estimating free parameters in ODE-based biochemical models, making probabilistic predictions about dynamical variables and discriminating between competing models having different topologies. We illustrate the use of this approach with a previously validated and non-identifiable model of receptor-mediated apoptosis in human cells (EARM1.3). Rather than return a single set of ‘best fit’ parameters, Bayesian estimation provides a statistically complete set of parameter vectors k that samples the posterior parameter distribution given a set of experimental observations (time-lapse data from live cells in the current case) and a value for experimental error. Estimation starts with a best-guess initial distribution (the prior), which is then modulated by a sum of squares log-likelihood function that scores the difference between model and data (Box 1). Recovery of the posterior parameter distribution makes it possible to compute confidence intervals for biologically interesting properties of the model (time and rapidity of apoptosis in the current case). These confidence intervals correctly account for measurement noise and parametric uncertainty, and can be remarkably precise in the face of substantial non-identifiablility ( Klinke, 2009 ). Simulations that include confidence intervals represent an advance on the prevailing practice of relying on error-free trajectories computed using a single set of maximum likelihood parameter values.

We also used Bayesian procedures to discriminate between competing direct and indirect models of MOMP, a key step in apoptosis. Discrimination involves estimating the evidence for the indirect model P(MI|data) divided by the evidence for the direct model P(MD|data), a ratio known as the Bayes factor ( Kass and Raftery, 1993 , 1995 Xu et al, 2010 ). Bayesian approaches to model discrimination account not only for uncertainty in parameter values but also for differences in the numbers of parameters. This is important because models that instantiate different hypotheses about biochemical mechanism usually have different numbers of parameters even when the number of unique model species is the same. The complexity penalty embedded in the Bayes factor represents a generalization of the Akaike or Bayesian Information Criteria (AIC and BIC) commonly used to score model likelihoods ( Kass and Raftery, 1995 ).

In the case of MOMP, we observed that both direct and indirect models fit experimental data equally well and, thus, cannot be distinguished on a maximum likelihood basis. However, computation of the Bayes factor revealed the direct model to be ∼16–24 times more probable than the indirect model at 90% confidence, reflecting the greater range of parameter values compatible with the direct model. This formalization of a ‘robustness’ criterion for preferring the direct model is consistent with recent experiments obtained from engineered cell lines ( Chipuk et al, 2010 ). With respect to the biological significance of this finding, however, it is important to note that published indirect and direct ‘word models’ are compatible with many different ODE networks. Thus, it will ultimately be necessary to distinguish among extended sets of competing models, not just the two presented here. With improvements in computational speed, methods for calculating the Bayes factor using ‘thermodynamic integration’ are well suited to this task.

Properties of the posterior distribution and implications for reporting parameter values

With EARM1.3, Bayesian estimation reveals substantial differences from one parameter to the next in the degree of identifiability (as reflected in the widths of the parameter distributions). This is expected given previous work showing that biochemical models are ‘sloppy’ ( Gutenkunst et al, 2007 ) even when calibrated against ‘complete’ data on all dynamic variables (which is not the case in the current work). A basic property of Bayesian estimation is that the posterior will resemble the prior when data provide little or no additional information. Conversely, when the data are informative, the shape of the posterior will differ substantially from that of the prior. In the case of EARM1.3 modal values for posterior distributions differed from the priors for about one-third of all parameters while still falling within a biophysically plausible range (with rate constants below the diffusion limit, for example). The exact shape of the prior did not appear to be critical in achieving convergent sampling, a fortunate outcome since we used general-purpose priors applicable to all cellular functions rather than priors derived from specific analysis of apoptosis proteins. One mistake we learned to avoid was constraining the MCMC walk to a fixed interval around nominal parameter values such hard limits result in artificially truncated marginal distributions. Gaussian priors also had the benefit of improving the fraction of parameter sets for which the EARM1.3 ODE system could be integrated.

It is incorrect to judge the impact of parameter estimation (i.e., what we learn by comparing models to data) simply by examining the shapes of individual parameter distributions: several lines of evidence show that marginal distributions contain only part of the information. Non-linear covariation among parameters accounts for the rest it is necessary for accurate model-based simulation and cannot be approximated by a covariance matrix. The reasons for this are evident from inspection of the landscape of the objective function. The landscape is steep (the eigenvalues of the Hessian matrix are high) in directions that do not point directly along raw parameter axes ( Gutenkunst et al, 2007 ). Thus, identifiable features of the systems correspond to ratios of rate constants (this is the basis of parameter covariation) but the value of the ratio varies through parameter space (this gives rise to curved high-probability valleys that are not well approximated by lines). By direct analogy, the identifiable parameter in a Michaelis–Menten treatment of a simple catalytic reaction is KM, a ratio of rate constants, rather than kf or kr themselves (( Chen et al, 2010 ) see also Box 2). When parameters in a model are highly covariant, it is almost always the case that the system can be described with a simpler model involving a smaller set of more identifiable parameters. In many applications in engineering it is desirable to use such reduced models, but in the case of biochemistry, parameter non-identifiability and high covariance appear to be the cost of representing systems as sets of reversible mass-action reactions. Under the assumption that mass-action kinetics (and also stochastic kinetics obtained by solving the chemical master equation) are uniquely appropriate as a means to describe the physics of biochemical systems, we are forced to use models such as EARM1.3. However, there is no reason, from a physical perspective, to believe that proteins in a network that do not actually bind to each other alter each other's rate constants. Thus, the presence of highly covariant parameters in best-fit vectors is not a property of the underlying biochemistry: instead, it represents a limitation on our ability to infer the properties of complex reaction networks based on the time-course data we typically acquire.

One consequence of parameter covariation in EARM1.3 is that parameter values in best-fit vectors do not correspond to the means of marginal parameter distributions, and sampling the means of marginal distributions does not result in a good fit. It is common practice in biochemical modeling to report parameters as a table of single values (with no allowance for non-identifiablility) or as a list of means and ranges. If these parameters are derived from calibration, critical information on covariation is lost. It is therefore necessary to report the actual vectors recovered by sampling the posterior parameter distribution. In principle, this is an array of size (C × M) × (N+1) where C is the number of MCMC chains, M the number of steps, and N the number of parameters (N+1 appears because we record a posterior value for each N-dimensional vector) corresponding to ∼1.5 × 10 8 entries for EARM1.3. However, steps in MCMC chains have characteristic ‘decorrelation lengths’ over which parameter values vary relatively little (∼10 2 –10 4 steps, depending on the parameter). Thinning by this amount yields an array of ∼10 4 –10 6 entries, still a much more complex representation of parameters than the simple tabular summary assumed by current standards such as SBML. It is also important to note that the posterior landscape needs to be revisited repeatedly when adding new data or extracting new hypotheses. In this sense, parameter estimates are best viewed as computational procedures and sets of possible values rather than fixed information.

Model discrimination in the Bayesian framework

A solid theoretical foundation and large body of literature speaks to the value of Bayesian frameworks for analyzing models having different numbers of uncertain parameters ( Kass and Raftery, 1995 MacKay, 2004 Xu et al, 2010 ). The Bayes factor described here, the odds ratio for competing models (i.e., the ratio of the evidence), is computed from an overlap integral between a likelihood and a prior (that is L(data|Θ) and π(Θ|M)). At first glance it might seem that making a model more complex would increase the number of dimensions and always increase the evidence, but a simple example shows that this is not the case. Consider a pair of one- and two-parameter models of the same hypothetical physical process and a function f that is the ratio of relevant likelihoods: f(k1,k2)=L(k1,k2)/L(k1). The evidence for the one-parameter model is the overlap integral between its likelihood and a normalized prior ∬dk1 L(k1)π(k1) and for the two-parameter model it is ∬dk1 dk2 L(k1,k2)π(k1)π(k2). In the case where f(k1,k2)<1 for all k1, k2 the likelihood of the two-parameter model is no better than that of the simpler one-parameter model (note that the evidence for the two-parameter model is ∬dk1 L(k1)π(k1)g(k1) where the function g(k1)=∬dk2 π(k2)f(k1,k2) must be less than 1 everywhere, as the priors are normalized to 1). The evidence for the two-parameter model will therefore be less than the evidence for the one-parameter model, meaning that it will lose out in a Bayes factor comparison, as it should.

When the function f(k1,k2)>1 then it must be true that introduction of a second parameter ‘rescues’ or ‘lifts’ the likelihood function by improving the fit to data. In this case, the more complex model will have greater evidence. In the special but interesting case where f(k1,k2)=1 for all k1, k2, model 2 is completely insensitive to the new parameter. The presence of a parameter with respect to which a model is completely insensitive has no impact in model assessment (the Bayes factor is one). Finally, in the general case where f(k1,k2) has values both above and below 1, explicit integration is needed to determine which model is favored, precisely what we do in this paper.

The Bayes factor is not unique as a means to balance goodness-of-fit and model complexity. The most commonly used metrics are the AIC and the BIC ( Akaike, 1974 Schwarz, 1978 ):

where n is the number of parameters, ML is the maximum likelihood value, and N is the number of data points (ML is simply the highest value achieved by the likelihood function). AIC and BIC do not explicitly account for parameter non-identifiability and the two metrics are therefore good for comparing models only in the ‘asymptotic’ limit where the number of experimental data points becomes ‘large’ and identifiability is ‘greater’ ( Akaike, 1977 , 1983 Kass and Raftery, 1995 ). It is rare in the field of biochemical modeling to explicitly check whether the conditions for computing the AIC and BIC are valid and, in our experience, they are frequently violated. In contrast, the Bayes factor is applicable even with limited data and reduces to the BIC and, under some conditions, to the AIC in the asymptotic limit ( Akaike, 1977 , 1983 Kass and Raftery, 1995 ). Moreover, whereas the AIC or BIC compare models solely on the basis of goodness-of-fit, Bayesian methods allows formal introduction of a prior degree of belief in each model. An arbitrary model (i.e., a physically impossible model) exhibiting a better fit to data might get a better AIC or BIC score than a more realistic mechanistic model, but in a Bayesian approach it would receive a low prior value. We therefore consider evaluation of the Bayes factor to be a better way than AIC or BIC to compare models when models have different numbers of non-identifiable parameters and data are limited.

Limitations of the approach

The computational approach described here has several practical and algorithmic limitations, albeit ones that can be mitigated with further work. A practical concern is that current methods for computing the Bayes factor are too slow to incorporate the full range of data we have collected from cells exposed to drugs, siRNA-mediated protein knockdown and ligand treatment. Using only a subset of available training data, computing the Bayes factor for EARM1.3 required ∼6 × 10 4 CPU-hr (4 weeks on a 100 core general-purpose computer cluster). It should be possible to improve this by optimizing the code (e.g., porting it from MatLab to C/C++) and performing multiple analyses in parallel. It also remains to be determined how inclusion of more calibration data will alter the topology of the posterior landscape. It may become more rugged, decreasing the requirement for Hessian guidance during MCMC walks but increasing the need for simulated annealing (SA) to move out of local minima. Readers interested in these developments should check our Github repository or contact the authors directly.

An additional conceptual concern with the current work involves the way in which MCMC walks sample the posterior landscape. To establish that sampling is correct, it is necessary to show that chains starting at independent positions converge. Convergent sampling is not an abstruse point because probability distributions can differ in shape and modal value, when sampling is convergent as opposed to non-convergent. With EARM1.3 we observed that convergence was not possible in a reasonable amount of time (e.g., a week-long cluster-based calculation) using a conventional MCMC walk. We therefore used an adaptive walk involving periodic recalculation of the local landscape as a means to guide MCMC walks and improve convergence. However, this approach may violate the detailed balance requirement of M–H sampling. With large models and existing methods, we are therefore in the position of having to choose between convergent Hessian-guided chains and partially non-convergent, conventional chains (( Klinke, 2009 ) chose the latter alternative). Moreover, using the Gelman–Rubin test to judge convergence has the weakness that it is one-sided: failing the test demonstrates non-convergence but passing the test does not guarantee it. Analysis of posterior distributions for EARM1.3 computed in different ways suggests that we are on relatively solid ground in the current paper (we did not observe significant differences in posterior distributions using different sampling approaches), but the development of methods for analyzing MCMC walks represents an active area of research in applied mathematics and it is necessary to be aware of future developments.

For reliable, probabilistic model-based simulation, we also need to consider the fact that the sufficiency of sampling is contingent not only on the model structure and available training data, but also on the types of predictions being sought. Assuming convergence, the posterior landscape sampled by multi-start MCMC walks represents a reasonable approximation to the true but unknown posterior distribution of the parameters, but the same is not necessarily true for predictions or simulated trajectories based on these parameters: the posterior landscape may be poorly sampled in regions of parameter space that have a significant impact on some simulations. In the current work we show that MCMC chains used to predict Ts and Td satisfy the Gelman–Rubin test, but this is a weak criterion and importance sampling using umbrella, non-Boltzmann or other methods ( Allen and Tildesley, 1987 ) will generally be necessary to revisit regions of the landscape that have low posterior values but contribute strongly to the distribution of a prediction. This suggests a workflow in which MCMC walks based on calibration data (as described here) are only the first step in model calibration. Inclusion of any new training data mandates a new round of estimation. Additional sampling should also be performed as required by importance sampling to reliably inform predictions. Finally, the use of formal methods for modeling experimental error ( Jaqaman and Danuser, 2006 ) should make it possible to distinguish errors arising from photon noise, infrequent sampling, incorrect normalization, etc, thereby improving the comparison between data and simulation.


The ubiquity of Bayesian methods in other scientific fields has motivated multiple, parallel efforts to apply the approach to biochemical models ( Flaherty et al, 2008 Calderhead and Girolami, 2009 Klinke, 2009 Xu et al, 2010 ), but significant challenges have remained with respect to development of widely available methods for discriminating between competing models. The algorithms and open-source code described in this paper enable a rigorous approach to estimating the parameters of non-identifiable biochemical models, making model-based predictions probabilistic, and using the Bayes factor to distinguish between models with different topologies and numbers of parameters. The generic priors we report are a good starting point for estimating forward, reverse and catalytic rates in any mass-action model of cellular biochemistry, but in some cases it makes more sense to substitute narrower and more specific priors derived from previous studies in vitro and in vivo such specific priors better capture pre-existing knowledge, improve parameter estimation and speed of computation. It is our opinion that application of rigorous probabilistic analysis of biochemical models will advance the long-term goal of understanding complex biochemical networks in diverse cell types and disease states ( Klinke, 2009 Xu et al, 2010 ). Preliminary application of Bayesian reasoning suggests that some long-standing disputes about cell signaling can be laid to rest, (e.g., direct versus indirect control of MOMP), whereas others truly cannot be discriminated based on available data.

5 Specifying Execution Settings

We have introduced the chemical perceptron models structurally as collections of species and reactions with catalysts and inhibitors. Now, we shall specify the setting of executions (or chemical experiments) in terms of the action series and the translation series (Section 2.2).

5.1 Binary-Function Modeling

Figure 6 presents simulations of the weight-loop perceptron and the weight-race perceptron, each computing the NAND function on four consecutive inputs.

Simulation of the NAND function on four different combinations of the input species by (a) the weight-loop perceptron with initial weight concentrations [W0 ⊕ ] = 6 M, [W1 ⊖ ] = 2 M, and [W2 ⊖ ] = 3 M, and (b) the weight-race perceptron with initial weight concentrations [W0 ⊕ ] = 0.3 M, [W1 ⊖ ] = 4 M, and [W2 ⊖ ] = 3.5 M. From top to bottom: concentrations of input species X1 0 and X1 1 , concentrations of input species X2 0 and X2 1 , and concentrations of output Y and weight species W. By applying the translation that compares the maximal concentrations of Y 1 and Y 0 in the four intervals 5,000 steps long, we obtain the NAND output sequence 1, 1, 1, 0.

Simulation of the NAND function on four different combinations of the input species by (a) the weight-loop perceptron with initial weight concentrations [W0 ⊕ ] = 6 M, [W1 ⊖ ] = 2 M, and [W2 ⊖ ] = 3 M, and (b) the weight-race perceptron with initial weight concentrations [W0 ⊕ ] = 0.3 M, [W1 ⊖ ] = 4 M, and [W2 ⊖ ] = 3.5 M. From top to bottom: concentrations of input species X1 0 and X1 1 , concentrations of input species X2 0 and X2 1 , and concentrations of output Y and weight species W. By applying the translation that compares the maximal concentrations of Y 1 and Y 0 in the four intervals 5,000 steps long, we obtain the NAND output sequence 1, 1, 1, 0.

5.2 Learning

In learning mode, the initial concentrations of weights are generated randomly in the interval 2–10 M, with equal probability of selecting positive and negative variants. A learning rate α, which is constant throughout the whole training, is incorporated into the concentration of the desired output D.

The original definition of perceptron learning (Section 3 adjusts weight wi by Δwi = α(dy)xi for a given output y and desired output d. Assuming dy, each weight participating in the output production increments by Δw = α(dy). Since the sign of the weight sum fully determines output, weight adaptation is stronger for inputs with the higher number of ones. For instance, the weight sum is adjusted by Δw for input (0, 0), but by 3Δw for input (1, 1), as shown in Table 4a.

The adaptation of a weight sum during learning by a two-input binary perceptron for four inputs: (a) a uniform adaptation of individual weights, (b) a uniform adaptation of the weight sum.

(a) . (b) .
x1 . x2 . Adapted weight sum . x1 . x2 . Adapted weight sum .
0 0 w0 + Δw0 0 w0 + Δw
1 0 w0 + w1 + 2Δw1 0 w0 + w1 + Δw
0 1 w0 + w2 + 2Δw0 1 w0 + w2 + Δw
1 1 w0 + w1 + w2 + 3Δw1 1 w0 + w1 + w2 + Δw
(a) . (b) .
x1 . x2 . Adapted weight sum . x1 . x2 . Adapted weight sum .
0 0 w0 + Δw0 0 w0 + Δw
1 0 w0 + w1 + 2Δw1 0 w0 + w1 + Δw
0 1 w0 + w2 + 2Δw0 1 w0 + w2 + Δw
1 1 w0 + w1 + w2 + 3Δw1 1 w0 + w1 + w2 + Δw

In the chemical perceptron, the concentration of the desired output is divided between weights hence the weight adaptation of the original perceptron would correspond to the concentration profile of the desired output for four consecutive inputs: [D] = α, [D] = 2α, [D] = 2α, [D] = 3α. Our simulations proved that this unfairness results in a split of performance for function pairs, such as CONST0 , CONST1 or AND , NAND . Since these functions are the inverse of each other, the chemical perceptron should learn them with the same success.

Essentially, the uniform adaptation of individual weights causes a bias in the weight adaptation. To avoid that, we do not adapt individual weights, but the whole weight sum uniformly for all inputs, as presented in Table 4b. More specifically, the chemical perceptron divides Δw, the concentration of the desired output D among weights, so the whole weight sum is adjusted by Δw. Note that a small amount of the desired output disappears because of a decay.

The concentration of the desired-output species D is constant for all inputs in the chemical perceptron. By experiments we determined that the optimal concentration of D is 2 M. If the concentration was too high, the weights would oscillate and would not converge on a stable solution. Conversely, a low concentration of D prolongs the learning process and does not provide enough pressure to drive weights out of the zero region if their concentrations are very low.

As opposed to our chemical perceptron, the biased adaptation in the original formal perceptron does not cause substantial problems, because the weight sum is further processed by an activation function and the learning rate α decreases over time. As a result, small differences in the weight adaptation become unimportant.

Learning consists of a series of actions, each randomly chosen from the training set and performed every S steps. The total number of actions, L, per action series is Lf = 120 for the fitness evaluation (Section 6), and Lp = 200 for the learning performance and robustness analysis (Sections 7.1 and 7.2), so the perceptron runs for either S × Lf = 6 × 10 5 or S × Lp = 10 6 time steps.

If we injected inputs together with the desired output at the same time, the adaptation of the weights would start immediately, changing the actual output, so the actual output would differ from the one we would obtain by providing only input species. In the extreme case, the chemical perceptron could just copy the desired output to the actual output by having very low concentrations of weights. To prevent this, we inject an input, wait a certain number of steps, measure the output, and then provide the desired output. Note that a reaction DW requires both catalysts, the input species X and the output species Y, to have a sufficient concentration at the moment of weight adaptation. Therefore, we must allow enough time for the output production, but we cannot postpone the injection of D for too long if we did, the chemical perceptron would process or decay all input species. We found experimentally that this delay can be fairly short. More precisely, in our learning simulations we inject the desired output 100 steps after the input species.

Now, the only question is how to interpret the output, or in other words, how the concentrations of species Y 0 and Y 1 translate to the value of the variable y from the formal definition. Since y is a Boolean variable, the translation compares the concentrations of value-0 species and value-1 species just before a desired output is injected. Hence, the value of the variable y is defined as [Y 1 ] > [Y 0 ] at relative time step 99.


While there are many caveats for deploying any allometric-scaled equations to describe complex behavioral and ecological processes, SCEB, has clear advantages over RHt2, its multi-prey derivations, and other functions that make reference to half-saturation constants with units of external prey abundance. While SCEB has similarities, and hence common mechanistic grounds, with the classic Holling type 2 disk equation, the use of a single handling constant (here labeled for “satiation”) and a defined maximum ingestion rate makes SCEB more attractive for deployment in planktonic ecosystem models. In addition to providing solutions to the challenges associated with deployment of RHt2, the simplicity of the SCEB function offers potential for:

(a) An extended theoretical consideration of predator-prey interactions, including explorations of activity in patches, and the challenges of mean-field computations of predation dynamics (Grünbaum, 2012).

(b) A role as a front-end to more complex models considering multiple prey types of different quality and quantity, together with satiation and variable assimilation efficiency, and changes in predator size and motility with developmental stage.

With no significant additional computational cost, the SCEB function overcomes all the shortcomings of the RHt2 function (and allies) and provides a direct route to involving physics, changes in predator and/or prey size and motility, and prey palatability (thus allowing true de-selection of prey types). SCEB can also be directly incorporated within individual-based as well as biomass-based models. Importantly, SCEB requires (enforces) an acknowledgement of allometrics in biomass-based models, which typically make scant reference to size even when (as for predator-prey interactions, as we have seen) it is necessary to do so.

Zero-Information Closure Scheme

For the sake of brevity, we limit the discussion in this section to one-dimensional state spaces. In particular, for a single random variable that can attain a discrete set of values, , each with probability , the information entropy is defined as (23): We conjecture that a finite number of probability moments may capture all of the information needed to precisely construct the probability distribution. In consequence, we maximize the information entropy under the constraints of the first M moment definitions: where is the Lagrange multiplier of the jth moment constraint. These are readily computed with appropriate root-finding numerical methods, such as the simple Newton–Raphson method we present in the SI Appendix, section S1.

Trivially, the final form for the maximum-entropy probability distribution is (24): Now the –th order moment can be computed from the maximum entropy distribution: Eqs. 2 are now closed and can be integrated to evolve the probability moments by a time step . Given an initial condition for the value of the probability moments, the system may be propagated in time. With newly calculated moments up to order M, the information entropy is maximized again, generating new values for the Lagrange multipliers, and so on.

It is important to note that, in a separate calculation without simulations in time, the steady state of reaction networks can be computed by determining the most likely distribution moments that lie in the null space of the augmented matrix . This means that the steady-state moments are found so that they maximize the information entropy as well as satisfy: We note that in this calculation there is no need to provide an initial condition for the probability distribution in state space. Because only a single optimization step is used this method is far more efficient than traditional Monte Carlo approaches.

The algorithm that couples the optimization with the ODE solver, and the algorithm for determining steady-state probability distributions are detailed in the SI Appendix, sections S2 and S3.


A number of methods exist for measuring or inferring A over a wide range of scales, including leaf gas exchange (Long & Bernacchi 2003 ), canopy chamber methods (Leadley & Drake 1993 Johnson, Polley & Whitis 2000 ), eddy covariance (Baldocchi 2003 ) and remote sensing techniques (Field et al. 1995 Damm et al. 2010 Frankenberg et al. 2011 Serbin et al. 2012 ). The basis for the ability to compare A among these spatial scales is through modelling. The leaf A model (Farquhar et al. 1980 ) provides the foundation for scaling A among a variety of spatial and temporal ranges while also providing the basis for generating predictions and hypotheses. A key attribute of the leaf A model is the mechanistic basis upon which equations are derived. In the context of scaling A to the globe, the leaf A model provides the backbone for predictions of GPP that are being validated against an increasing number of remote sensing techniques. Recent advances in remote-sensing methods to acquire the parameters needed to model A (Vc,max and Jmax e.g. Serbin et al. 2012 ) will further help to constrain models using the most accurate parameterizations.

Watch the video: Εισαγωγή στην έννοια του αλγορίθμου και του προγραμματισμού (July 2022).


  1. Faolan

    Problem is, a quick answer :)

  2. Akibei

    I congratulate this idea just about

  3. Atteworthe

    the choice at home difficult

  4. Fitz

    It is not joke!

  5. Taucage

    I agree, this is a great opinion

Write a message