Abstract
Background
Stochastic effects can be important for the behavior of processes involving small population numbers, so the study of stochastic models has become an important topic in the burgeoning field of computational systems biology. However analysis techniques for stochastic models have tended to lag behind their deterministic cousins due to the heavier computational demands of the statistical approaches for fitting the models to experimental data. There is a continuing need for more effective and efficient algorithms. In this article we focus on the parameter inference problem for stochastic kinetic models of biochemical reactions given discrete timecourse observations of either some or all of the molecular species.
Results
We propose an algorithm for inference of kinetic rate parameters based upon maximum likelihood using stochastic gradient descent (SGD). We derive a general formula for the gradient of the likelihood function given discrete timecourse observations. The formula applies to any explicit functional form of the kinetic rate laws such as massaction, MichaelisMenten, etc. Our algorithm estimates the gradient of the likelihood function by reversible jump Markov chain Monte Carlo sampling (RJMCMC), and then gradient descent method is employed to obtain the maximum likelihood estimation of parameter values. Furthermore, we utilize flux balance analysis and show how to automatically construct reversible jump samplers for arbitrary biochemical reaction models. We provide RJMCMC sampling algorithms for both fully observed and partially observed timecourse observation data. Our methods are illustrated with two examples: a birthdeath model and an autoregulatory gene network. We find good agreement of the inferred parameters with the actual parameters in both models.
Conclusions
The SGD method proposed in the paper presents a general framework of inferring parameters for stochastic kinetic models. The method is computationally efficient and is effective for both partially and fully observed systems. Automatic construction of reversible jump samplers and general formulation of the likelihood gradient function makes our method applicable to a wide range of stochastic models. Furthermore our derivations can be useful for other purposes such as using the gradient information for parametric sensitivity analysis or using the reversible jump samplers for full Bayesian inference. The software implementing the algorithms is publicly available at http://cbcl.ics.uci.edu/sgd webcite
Background
It is becoming increasingly apparent that stochasticity, whether intrinsic or extrinsic, plays an important role in the dynamics and behavior of biological systems. In systems biology and the study of gene expression [13], the consequences of stochasticity can manifest in numerous ways such as slow promoter kinetics leading to gene transcription bursting [4,5], finitenumber effects and mRNA translation bursting [69], propagation of noise in gene regulatory cascades [4,10], and phenotypic switching [11,12]. In some cases, biological systems evolve to minimize the effects of noise such as through negative feedback loops [1315], but there is also evidence that biology exploits randomness such as to create phenotypic diversity in populations thus allowing better adaptation to changing environments [1618]. With the growing awareness of stochasticity in biology and the increasing use of stochastic models in computational systems biology, there is a need to develop new analysis and computational techniques for studying, understanding and designing these stochastic models.
One particular analysis technique and challenge in computational systems biology is the inference of rate parameters from experimental data for a specified biochemical system [19]. Parameter inference for continuous deterministic models has a considerable body of research literature and can often be converted into an optimization problem for which many computational methods are available [20]. The strategies of these methods can be classified as either deterministic or stochastic. Deterministic strategies are generally only applicable for specific mathematical formulations of the model where a statement about the existence of the global optimum can be guaranteed along with a constructive algorithm to find it. Many problems are not that welldefined so stochastic strategies are popular including stochastic gradient descent [21], simulated annealing [2224], evolutionary computation [25], and other heuristics. Regardless, considerable computational effort is required for all of these methods as many simulations of the continuous deterministic model are performed. A discrete stochastic model is essentially a more adequate description for a biochemical system, but it has the disadvantage of being computationally expensive to simulate as well as requiring numerous independent simulations to be performed in order to calculate expectation values of various model outputs [2628]. These computational challenges mean that approximation techniques are frequently used for parameter inference including simplification of the stochastic model [29] and approximate inference such as using the chemical Langevin equation [30] in place of the Markov jump process [31,32]. Recent research has shown that parameter inference for stochastic models is feasible given time course observations of the system, even if only a partial set of molecular species are observed [32,33]. However the current algorithms, based on the Bayesian framework, are typically timeconsuming due to the need of sampling high dimensional space. Therefore there are significant challenges in applying the method to real systems, such as gene regulatory networks [34].
Most proposed methods for parameter inference in stochastic biochemical models consider how to calculate the maximum likelihood for the rate parameter values given a stochastic model and observational data. Except for the simple models, the likelihood function is computationally intractable, so these methods either perform exact inference on an approximated model where the likelihood computation is tractable, or they approximate the likelihood with a more tractable function, or some combination of the two. Tian et al. [35] considered the simulated maximum likelihood (SML) method that estimates likelihood by generating samples from many simulations of the stochastic model. The ratio of samples matching observations to the total number of samples is used to estimate the transitional density and the loglikelihood. Then a genetic algorithm is used to obtain the optimal rate parameter values that minimize the loglikelihood function. While the SML approach is straightforward, it is computationally expensive because it requires a large number of simulations of the stochastic model. Similarly, approximate Bayesian computation also requires the stochastic model to be simulated, but it avoids calculating the likelihood function by comparing simulated data with observations using a rejection sampler [36,37]. In a similar framework, Yosiphon et al. [38] used a simulated annealing procedure in an MCMC algorithm to estimate the parameters in stochastic models of reaction networks. Reinker et al. [39] proposed a method utilizing a hidden Markov model to approximate the stochastic model that takes observational error into account. Boys et al. [33] showed how full Bayesian inference can be performed on the stochastic LotkaVolterra model along with performance of various Markov chain Monte Carlo (MCMC) algorithms. Interestingly they showed that with partially observed data, i.e., only one of the two species in the model, they can still make inferences about all three rate parameters in the model; though it is unclear how well this would work on larger models with many parameters. Wilkinson and colleagues have investigated additional methods including using diffusion approximations [29,31] and incorporating multiple data sources [40].
In this paper, we describe an alternative method for parameter inference in discretely observed stochastic kinetic models. Instead of calculating and approximating the likelihood function as in the previous methods, we focus on estimating the gradients of the likelihood function with respect to the parameters. In particular, we propose a general methodology for efficiently estimating the gradients using reversible jump Markov chain Monte Carlo (RJMCMC). RJMCMC is an extension of the standard MCMC method that allows for generating samples on spaces of varying dimensions [41]. An implementation challenge for RJMCMC is the lack of a general way to construct the jump proposals such that detailed balance is preserved [42]. For stochastic kinetic models, the jump proposal corresponds to moves that change the number and the time of reactions that occur between two observations of the system. For most models, there is an infinite set of possible reaction processes (constrained by the observation data) that can occur between two time points, and the probabilities of different reaction paths depend upon the rate parameter values. Utilizing the research in flux balance analysis for metabolic networks [4345], we provide an algorithm so that jump proposals can be automatically constructed from any standard biochemical model, thus allowing RJMCMC to be used without requiring any manual analysis by the modeler.
The availability of the gradient information allows for inference of the rate parameters of stochastic kinetic models using gradient descentbased methods. We implement a steepest gradient descent method for parameter inference using the estimated gradient information in a MATLAB software package http://cbcl.ics.uci.edu/sgd webcite. We demonstrate the utility of our algorithms using two example stochastic models, including a birthdeath process and a gene autoregulation model.
Methods
Stochastic kinetic model of reaction systems with discrete states
Consider a general reaction system involving M reactions R_{1}, R_{2},..., R_{M }and K species S_{1}, S_{2},..., S_{K}. We denote the state of the reaction system by X = (x_{1},⋯,x_{K}) where x_{a }is the number of species S_{a}. Each reaction R_{i }has an associated rate law, represented by a hazard function h_{i}(X, Θ) (also called rate function), where Θ ≡ {θ_{r}} is a set of parameters associated with the reactions. Suppose the reaction system takes the following form
where u_{ra }and v_{ra }are the positive integer stoichiometries associated with reaction R_{r }for reactant S_{a}, representing the amount of species S_{a }that decrease and increase respectively when reaction R_{r }occurs. Eq. (1.1) can be represented more compactly as US → VS, where U = [u_{ij}] and V = [v_{ij}] are M × K matrices. We define the net effect reaction matrix A = V  U, which reflects the net change of species numbers associated with reactions.
Asumming the reaction system is in wellstirred condition with a fixed volume, we can introduce the master equation model, also known as "chemical master equation (CME)" in the biochemical modeling field [26], which describe the time evolution of the state probability using a set of ordinary differential equations. The CME can be derived for any biochemical reaction system using the standard continuous time Markov process theory. Denote P(X;t) the probability of the system in state X at time t. For an infinitesimal time increment Δt, P(X;t + Δt) can be written as the sum of probabilities of the number of ways in which the system can reach or leave the current state:
where A_{i }denotes the ith row of the net effect matrix A, and h_{i}(X, Θ) is the hazard function, determining the rate of probability transition out of state X due to reaction type i. In the limit of Δt → 0, Eq. (1.2) adopts the standard master equation form,
with
Suppose all possible system states (usually countably infinite) are ordered and represented by indices 1, 2,..., etc. Then Eq. (1.3) can be rewritten as , where is a row vector with P_{i}(t) representing the probability of the ith state at time t.
For reactions that obey massaction law kinetics, one rate parameter θ_{r }is associated with each reaction type i, and consequently the hazard function has the form of
where u_{ra }is the stoichiometry coefficient of reactants a in reaction R_{r}. Forms of other rate laws for chemical kinetics, e.g. the MichaelisMenten model, can be found in [46]. Although we will focus our discussion on the hazard function in the form of Eq.(1.5), the following analysis can handle more general cases as long as the explicit functional form of the hazard function is known.
Gradient of the likelihood function with discrete observations
Our goal is to estimate the rate parameters of a stochastic model based on the observations at a set of discrete time points. Suppose we have observations {X_{Γ}(t_{1}), X_{Γ}(t_{2}),⋯,X_{Γ}(t_{m})} of the system at m discrete time points {t_{1}, t_{2},...,t_{m}} for a subset of species Γ ⊆ {1,⋯,K}. We say the system is fully observed if Γ = {1,⋯,K}, and partially observed otherwise. Denoting the likelihood of the observations for a given set of rate parameters by L(X_{Γ}(t_{1}), X_{Γ}(t_{2}),⋯,X_{Γ}(t_{m}); Θ), we estimate the rate parameters by maximizing the likelihood function.
For simplicity of discussion, consider first a single time interval [t_{s}, t_{s+1}] with full observations available at the start and the end of the interval, denoted by X(t_{s}) and X(t_{s+1}) respectively. Let L(X(t_{s}), X(t_{s+1});Θ) denote the likelihood of observing X(t_{s}) and X(t_{s+1}) under a model with parameters Θ. In Appendix, we show that the gradient of the likelihood function with respective to parameters can be calculated using the following formula, for any stochastic system with a master equation Eq. (1.3)
where T_{k }is the time duration of the system at state k, and N_{k, k' }is the number of transitions from state k to k' occurred during the interval. Both T_{k }and N_{k, k' }are random variables, and can be viewed as the sufficient statistics of the model. E[·] represents the expectation of the random variables. The formula suggests that we can calculate the gradient of the likelihood function by estimating the expectations of the two sufficient statistics.
For the biochemical reaction system in Eq. (1.1), suppose J reactions have occurred during the time interval [t_{s}, t_{s+1}] with the types and the corresponding times of the reactions denoted by . Then by Eq. (2.1), the gradient of the likelihood function can be rewritten as
where , which is fully specified by Ξ, denotes the state of the system between , and . Eq. (2.2) can also be written in an alternative form
where
is the likelihood of the reaction process Ξ. If all the reactions follow massaction law in Eq.(1.5), the gradient formula can be further written as
Now return to the general case where the observations are available at multiple time points from a subset of the species. The above formula for calculating gradient can still hold if we view the entire duration of the observations as a single time interval. However, the expectation in Eq. (2.3) is now taken on the systems states whose distribution is conditioned on the observations at the intermediate time points.
In general, the expectation in Eq. (2.3) cannot be calculated exactly. Instead we utilize a sampling method to approximate the expectation. More specifically, we sample the latent path conditioned on the parameters and the observations, and then calculate the quantity in Eq. (2.3) by averaging over the sampled paths to obtain the gradient. The same strategy also applies to the partially observed case, as long as the reaction paths are sampled conditioned on the partial observation data.
Reversible jump Markov chain Monte Carlo sampling
To calculate the gradient, we need to find an efficient way to sample the latent reaction processes conditioned on the observations. One commonly used sampling method is the stochastic simulation algorithm (SSA) [26], which can be used as a rejection method to discard samples that do that match the end state. The SSA method is computationally inefficient for generating samples between two measurements when the total number of possible states is high (as in the case of the biochemical reactions), because the chance of a sampled trajectory matching the end state is typically small and consequently most of the samples will likely be rejected.
Here we use the framework of RJMCMC [41] to sample the latent process. RJMCMC is a generalized MCMC method that can construct a sampler between models of different dimensions, which in our case corresponds to reaction paths with different number of reactions. To sample latent paths in biochemical reaction systems, the RJMCMC method [33] first generates an initial reaction path that is consistent with the observations. Then RJMCMCM constructs a Markov chain by a) proposing a new sample path by adding or deleting a specific set of reactions from the current path, and b) determining whether to accept the new sample or keep the previous one according to an acceptance probability.
Therefore, to construct a RJMCMC sampler, we will need to consider three issues: 1) how to generate the initial path; 2) how to propose a set of reactions for addition or deletion; and 3) how to determine the acceptance probability of a new path. Note that both the initial path and the proposed path have to match the observations at the start and the end of the interval, implying that only a subset of the reactions can be used for either initialization or addition/deletion. While the RJMCMC sampler exists for some specific reaction systems [33], usually taking advantage of the domainspecific knowledge, the challenge, however, is to find a general method that can work for any arbitrary reaction system.
Next we address the three issues mentioned above, and describe a general method to automatically construct a RJMCMC sampler for an arbitrary reaction system.
1) Generating an initial reaction path using integer programming
The first issue of generating the initial path is relatively easy to address. Let r be a vector representing the number of each reaction type occurred within the initial path. To match the observations at the start and the end of each interval, r has to satisfy certain constraints. Fortunately, all these constraints are linear, and thus we can use linear integer programming to find a solution. In practice, we used the GNU Linear Programming Kit (GLPK library) [47], which is incorporated into our MATLAB package using the interface GLPKMEX [48].
2) Proposing a new sample by adding or removing reactions
After an initial path is generated, our next step is to use proposal moves to add or remove reactions. Before describing our method, we first introduce two concepts that are used in studying biochemical reaction systems.
Definition 1: Elementary Mode
An elementary mode (EM) of a biochemical reaction network is a set of reactions that does not alter the observed number of molecular species. Formally, an elementary mode is a column vector of nonnegative integers that satisfies , where Ã = A ( the neteffect reaction matrix) when all species are observable, and is a submatrix of A with columns corresponding to the observed species when only a subset of species are observable.
Definition 2: Null Set
The null set is a set consisting of all independent elementary modes, denoted by .
Note that the null set is usually different between the fully and the partially observed case because of the different Ã matrix used.
Elementary modes analysis is well studied in metabolic networks theory and is used to find the flux distribution of the metabolic network at a steady state [49]. Various tools have been developed to identify EMs [4345]. In this work, we used the metatool package [44] to calculate the null set of any specific reaction network, which has been shown to be efficient for large networks.
Provided with a valid reaction path and the null set, we then proceed to generate a new sample by taking one of the following three move types. After randomly choosing an elementary mode from the null set,
1. With probability α_{1}, add the set of reactions in with random reaction times uniformly distributed within the interval.
2. With probability α_{2}, remove one set of randomly selected reactions in from the current path within the interval.
3. With probability 1α_{1}α_{2}, randomly move the time of all reactions.
Using ensures that the proposed reaction path is always consistent with the observations. However, there are two additional conditions for a new sample path to be valid: 1) the number of any reaction type must be positive after the move, and 2) the population numbers for all species remain positive throughout the whole process. If either of the two conditions is violated, we set the likelihood of the new sample path to be zero and reject the new sampled path. The proposal probability in RJMCMC for different moves is set be to α_{1 }= α_{2 }= 0.25 in practice. Note that the initial path and the null set only need to be calculated once, and thus they only impose a modest computational burden on the sampling algorithm.
3) Determining acceptance probability
Next we address the third issue on how to determine the acceptance probability of a proposed sample. We discuss the fully observed case first, and then the partially observed case.
Fully observed case
The observations at m discrete time points break the entire observation window into m1 subintervals. Because all species are observed, the reaction path at each subinterval is completely independently of each other conditioned on the observations. The reaction path at each subinterval can therefore be sampled independently using RJMCMC. Let Ξ denote the current reaction path and Ξ denote the proposed reaction path. The probability of accepting the new path is specified by min(1, AR_{p}), with p = 1, 2, or 3 denoting the type of the move
where π(ΞΘ), defined in Eq. (2.3), is the likelihood of sample path Ξ, r_{j }is the number of type j reaction in the current sample path, q_{k, j }denotes the number of reaction type j in the elementary mode , and τ is the time length of the subinterval. Appendix, Algorithm 1 provides the pseudocode for the fully observed case.
Partially observed case
In the partially observed case, observations are only available for a subset of the species. Different from the fully observed case, the reaction paths at different subintervals are now correlated, caused by unobserved species. Consequently, RJMCMC can no longer be applied independently for each subinterval.
To account for the correlation, we use a new strategy in which the reaction paths at two consecutive subintervals are sampled together at each sampling step using correlated moves. Let {q_{k}'}, k ∈ (1, K') be the null set corresponding to the partially observed case. Note that adding/deleting the set of reactions in only ensures that the observed species' numbers remain unchanged, but not the unobserved species. Suppose we are to update the reaction path following the time point t_{i}. We first generate a new sample path in the ith interval [t_{i}, t_{i+1}] using the same reversible jump moves as described for the fully observed case, with a randomly chosen elementary mode. If the move changes the unobserved species numbers at time t_{i+1}, we subsequently update the (i+1)th interval using a complementary move that keeps the system state at the end of the second interval unchanged. For example, if move type 1 (or 2) is chosen to update the first interval with an elementary mode q'_{k}, move type 2 (or 1) will be applied to the second interval to remove (or add) the same elementary mode q'_{k}. The complementary moves guarantee that the new reaction paths proposed for the two subintervals do not change the species numbers, including those of the unobserved species, at the end of the second interval. As with the fully observed case, the two conditions of a valid path (positive reaction type number and positive species number) must be satisfied, otherwise the proposal move will be rejected. The acceptance probability is calculated as , where p' denotes the complementary move type of p. In this way, the state of unobserved species at time t_{i }(i = 2,⋯,m) can be updated sequentially. An additional step is used to update the state at the first observation time point t_{1}, which is done by keeping the species number at the end of the first interval fixed and changing the start state according to the proposed move. Appendix, Algorithm 2 provides the pseudocode of using RJMCMC for the partially observed case.
Stochastic gradient descent algorithm
Given the estimated gradient of the likelihood function, we use the method of steepest descent to find an optimal solution of the parameters. At each step of the algorithm, we first generate sample paths using the RJMCMC algorithm at current parameter values. After burnin, we calculate the gradient of the likelihood function using the formula in (2.4). The estimated gradient is then used to update the parameter values until convergence. A simple strategy for choosing the step size is to set it to be a constant. Although this works well for simple systems, it sometimes induces overshooting of the parameter values or slow convergence during the gradient descent. When this happens, we adaptively adjust the step size within a certain range according to the gradient value. An overview of the stochastic gradient descent algorithm is given in Appendix, Algorithm 3.
Results
Next we illustrate the utility of our algorithm using two example reaction systems. In both cases, we simulated the reactions of the system using the stochastic simulation algorithm, and recorded the species numbers at a set of discrete time points, which were treated as observations of the system. Our method was then applied to infer the rate parameters for each system based on these observations.
Example 1: Birthdeath process
We first applied our algorithm to a wellstudied birthdeath process, which can be seen as a simplified model of production and degradation of a single molecular species [34]. The reactions are
We assume that R_{1 }and R_{2 }follow the firstorder and zerothorder massaction law respectively. Denote the number of A molecules by n_{A}, thus the hazard function is given by h_{1 }= k_{1}n_{A }and h_{2 }= k_{2}. The neteffect reaction matrix of the system is A = (1, 1)^{T}. Consequently, the null set of the system contains only one elementary mode , i.e. the combination of R_{1 }and R_{2}.
We generated observations by simulating the reaction process using SSA with different parameter sets (k_{1}, k_{2}) = (0.03,0.6), (0.06, 0.6), (0.1, 0.6), (0.1, 0.3) and (0.03, 0.2). For each parameter set, four observation datasets were generated that differ on the total observation time (T) and the observation interval (Δt) (see Table 1).
Table 1. Parameter inference result for the birthdeath model
We first examined the convergence property of the RJMCMC sampler with the different datasets generated with the first parameter set. Figure 1 shows the trace plots and autocorrelations of the total number of reactions in the sample paths. In all cases, we found the RJMCMC sampler is efficient and induces good mixing of sample paths, with convergence occurring typically within 100 samples. As expected, larger correlation length is observed for data with longer observation intervals (Figure 1d and 1f).
Figure 1. Property of the RJMCMC sampler for the birthdeath model. Panel a, c, and e plot the total number of reactions in each of the samples generated by RJMCMC, corresponding to three different observation datasets with Δt = 2, 5 and 10 respectively. Panel b, d, and f shows the autocorrelation on the number of reactions corresponding to panel a, c, and e respectively. The rate parameters used are (k_{1}, k_{2}) = (0.03, 0.6).
We applied the SGD algorithm to estimate the two rate parameters for each dataset. The convergence criterion is set to be that the relative changes of all parameter values are less than 0.005. We used 1000 samples after a burnin of 100 samples to estimate the gradient for a given set of parameter values. The estimated parameters for each dataset are summarized in Table 1. In all cases, the inferred parameters showed a good agreement with the true values, although the accuracy of the estimation clearly correlates with the number of observations and the observation time intervals. For datasets with larger observation interval and fewer data points, larger variation is the observed for the inferred value between different datasets, indicating the parameters are less constrained in these cases (results not shown). Additional file 1, Figure S1 shows a typical gradient descent run using the one of the datasets generated with (k_{1}, k_{2}) = (0.03,0.6), which consists of 21 data points with a total time period of T = 40. We observed that the parameters converge very quickly during the gradient descent, typically within 20 steps for our tested random start values.
Additional file 1. Supplementary figures and tables. This file contains Supplementary Figure S1S4, Table S1.
Format: DOCX Size: 977KB Download file
Example 2: Prokaryotic autoregulatory gene network
The second model we tested is a prokaryotic autoregulatory gene network in which dimmers of a protein repress its own gene transcription by binding to a regulatory region upstream of the gene. The system, involving both transcription and translation, can serve as a simple, yet illustrative, example of gene regulation [31,39]. The reactions in the network are given below:
Here DNA, P, P_{2 }and mRNA represent promoter sequences, proteins, protein dimmers and messenger RNA respectively. In this model, mRNAs and proteins are synthesized by transcription and translation processes (R_{3 }and R_{7}), and destroyed by degradation (R_{4 }and R_{8}). The proteins can form a dimmer P_{2 }(R_{5 }and R_{6}), which binds and unbinds to DNA (R_{1 }and R_{2}). When a protein dimmer binds to the promoter, it represses mRNA production. Overall, the network implements a selfregulatory mechanism to control the synthesis of the protein product, suppressing the transcription when the protein product is abundant. Note that DNA_{t }= DNA + DNA.P_{2 }is a conserved quantity in the system. The rate functions of reactions are assumed to follow massaction law with rate parameters k_{1 }to k_{8}, e.g. h_{1 }= k_{1}·P_{2}·DNA.
We applied our algorithm to both the fully and partially observed cases. We generated 10 datasets as observations within a time window of [0 50) with (k_{1},⋯,k_{8}) = (0.1, 0.7, 0.143, 0.35, 0.3, 0.1, 0.9, 0.11, 0.2, 0.1). Datasets D_{1 } D_{5 }have total copy number DNA_{t }to be 10 with the time interval between observations (Δt) from 1.0 to 0.1. The other five datasets D_{6 } D_{10 }are generated with DNA_{t }= 2. Detailed information of the datasets is shown in Table 2. For the partially observed case, we assume that only three of the species, mRNA, P and P_{2}, are observed. In addition, we assume that the total copy number DNA_{t }is known to avoid systematic bias in the sampling the system. While using the same datasets in the fully observed case, we only retain the observations corresponding to mRNA, P and P_{2}. Hereinafter, we denote the datasets by D_{1}*, D_{2}* etc.
Table 2. Parameter inference result for autoregulatory gene network model (Fully observed case)
For the fully observed case, the net effect reaction matrix is, shown with the corresponding reactions and species
The corresponding null set contains four elementary modes, consisting of the following four pairs of reactions: R_{1 } R_{2}, R_{3 } R_{4}, R_{5 } R_{6}, and R_{7 } R_{8}.
We focus our analysis on datasets D_{1 } D_{5}, of which the observation intervals range from 1.0 to 0.1. The results from datasets D_{6 } D_{10 }are similar. The convergence property of the RJMCMC sampler is shown in Figure 2. It shows the RJMCMC sampler is efficient and induces good mixing for all the datasets, with convergence occurring typically after 200 samples. The correlation lengths between the samples are smaller for the more densely observed dataset (D_{5 }with Δt = 0.1). The correlation length increases for the datasets with increasing Δt, suggesting the need of using larger sample size for sparse observed datasets.
Figure 2. Property of the RJMCMC sampler for the fully observed case of the autoregulatory gene network. Panel a, c, and e plot the total number of reactions in each of the samples generated by RJMCMC, corresponding to three different observation datasets with Δt = 0.1, 0.5 and 1 respectively. Panel b, d, and f shows the autocorrelation on the number of reactions corresponding to panel a, c, and e respectively. The rate parameters used are (k_{1},...k_{8}) = (0.1, 0.7, 0.143, 0.35, 0.3, 0.1, 0.9, 0.11, 0.2, 0.1).
We applied the stochastic gradient descent method to estimate the rate parameters given the observations. The initial parameter values were randomly chosen between 0.1 and 10. We used 5000 samples to calculate the gradient with a burnin size of 200. The estimated parameters for each dataset are summarized in Table 2. We observed a good agreement between the estimated and true values for most of the parameters. Also we observed that there is differences in the estimation accuracy for different parameters, with some (e.g. k_{2}, k_{3 }and k_{4}) showing consistently better results than others (Additional file 1, Figure S2). The estimation of k_{5 }and k_{6 }showed large deviation for the first two datasets with large observation interval (D_{1 }and D_{2}, Table 3), but improved with finersampled data. This is likely due to the faster dynamics of the two reactions (R_{5 }and R_{6}) than other reactions in the system.
Table 3. Parameter inference result for autoregulatory gene network model (Partially observed case)
Next we applied our algorithm to the partially observed case. The convergence property of the RJMCMC sampler for the partially observed case is shown in Figure 3. Compared with the fully observed case, the autocorrelation length in the partially observed case is typically longer, but the RJMCMC sampler can still induce good mixing for each dataset with adequately large sample size.
Figure 3. Property of the RJMCMC sampler for the partially observed case of the autoregulatory gene network. Panel a, c, and e plot the total number of reactions in each of the samples generated by RJMCMC, corresponding to three different observation datasets with Δt = 0.1, 0.5 and 1 respectively. Panel b, d, and f shows the autocorrelation on the number of reactions corresponding to panel a, c, and e respectively. The rate parameters used are (k_{1},...k_{8}) = (0.1, 0.7, 0.143, 0.35, 0.3, 0.1, 0.9, 0.11, 0.2, 0.1). Only three species, mRNA, P and P_{2}, are observed.
The parameter inference results are summarized in Table 3. We found that the accuracy of the inferred parameter varies for different datasets. For the densely observed dataset D_{5}*, the estimated values of all eight parameters are similar to those in the fully observed case and close to the true values. But for more sparsely observed datasets, the average percent of error of the inferred parameters increases significantly (compared to the fully observed case) for some of the datasets (D_{1}* and D_{4}*). The parameters k_{1 }and k_{2}, which are associated with the unobserved species, showed large variations between different datasets. In general, the results showed that parameter inference with partially observed data is more difficult than the one with fully observed data, and to achieve good estimation accuracy, more observations with small observation intervals will be needed.
Additional file 1, Figure S3 shows the changes of parameters and gradients during one gradient descent run for the most sparsely observed dataset D_{1 }and D_{1}* with the dataset (copy number of each species) shown in Additional file 1, Figure S4. Some of the parameters (e.g. k_{2 }and k_{6}) showed slow convergence during gradient descent in both fully and partially observed cases, which may reflect a flat likelihood surface in the corresponding parameter direction and an inherent difficulty in identifying these parameters.
Discussions
Recently there has been a growing interest in describing biological systems using stochastic models. However, most of the parameters in the stochastic models are unknown and difficult to measure. In this paper we described a maximum likelihood method to infer the parameters of a stochastic kinetic model directly from observations. Our method works by estimating the gradient of the likelihood function first, and then searching for an optimal solution by iteratively updating the parameters along the gradient descent direction. We developed a general RJMCMC algorithm to sample the latent reaction path in a constrained setting, where the reaction path has to match the observations given at the two ends of a time interval. The sampled reaction paths are used to calculate the gradient of the likelihood function using a formula that we derived. The availability of the gradient information makes it possible to develop other algorithms to solve the maximum likelihood estimation problem, in addition to the steepest descent method that we implemented. Furthermore, the availability of the gradient information also enables other possible applications such as parameter sensitivity analysis, which has already attracted considerable interest in deterministic modeling [50,51].
Our method is significantly faster than the SML method [34], which is also a maximum likelihood based parameter inference method. SML uses two steps to estimate parameters. First, it estimates the transition density on reaction species numbers after a given time interval, using a SSAbased sampling methods. The estimated transition density is then used to calculate the likelihood function. Because the gradient of the likelihood function is not directly available, SML uses a genetic algorithm to solve the maximum likelihood problem. Comparing the SML and our method for the birthdeath example, we tested the CPU time used to generate a new sample for both methods, eg. SSA for SML (unconstrained) and RJMCMC for SGD, which is approximately the same. However, SML uses 3 × 10^{4 }evaluations of transition density to reach a solution. By contrast, SGD typically requires less than 20 evaluations of the gradient before convergence. If we ignore the computational time of the gradient descent steps, overall SGD achieves a reduction of computational time by an order of 10 ^{3 }compared to SML.
In terms of accuracy, our approach, based on exact sampling, should be less biased than approximationbased methods. In this regard, we compared SGD to the method by Golightly et al. [31], who used a diffusion approximation to calculate the transition density. Comparing the results obtained by both methods on the same datasets (in courtesy of Dr. Golightly), we note that the estimated values for k_{1 }and k_{2 }by our method are closer to the true results in all three test datasets while the result from [31] are biased toward low values, although the estimates for other parameters from the two methods are similar (Additional file 1, Table S1). Interestingly, k_{1 }and k_{2 }are associated with low copy number of reaction species (DNA and DNA.P_{2}). We also tested the method in [31] with the datasets of DNA_{t }= 2 and found that the algorithm gives worse results, especially for the first two parameters (result not shown). This reflects the advantage of our method and possibly the limitation of the diffusion approximation, which assumes that the values of the hazard functions are approximately constants between two observation/latent states. This assumption is not valid if the copy numbers of species are small in the reactions. For example, in case of DNA_{t }= 1, reactions R_{1 }and R_{2 }can only happen alternatively and this clearly violates the approximation assumption.
Our method is closely related to the full Bayesian approach proposed by Boys et al. [33] as both methods use RJMCMC to sample the reaction process. Comparing to the method by Boys et al., our method offers two improvements. First, we provide a general method for RCMCMC sampling, which can be applied to an arbitrary biochemical reaction system, while the previous method is only tailored to a specific reaction system (more specifically, the LotkaVolterra system). Second, the gradientbased method is significantly faster than the full Bayesian method as sampling the parameter space is often computationally challenging. However, the Bayesian approach offers certain advantage over the maximum likelihood method in that it provides a posterior distribution of the parameters rather than just an optimal solution. In this regard, we note that the general RJMCMC sampling method we developed can be easily extended for Bayesian inference after introducing additional MetropolisHasting steps for sampling parameters.
Conclusion
In this paper, we proposed a new algorithm for inferring rate parameters in stochastic models and tested it using simulated data. Although few biological systems with measurements of species numbers across multiple time points are currently available, this type of data will likely become more common in the future, given rapid advances in single cell measurement technology [9,52,53]. The method could also be applied to cell colony data, e.g. in [54], which proposed some interesting models involving stem cell homeostasis process. As we observed, the current RJMCMC sampler can be inefficient in some cases with large observation intervals. One possible improvement of the current algorithm is to use more efficient sampling algorithm, for example, the blocking updating scheme in [33]. It is evident that significant challenges remain in dealing with true biological systems, including measurement noise, uncertainty in models, and sparsity of the data. However, studying stochastic systems with parameters inferred directly from data should be able to lead to a better understanding of the systems than the current approach of manually setting these parameters.
Authors' contributions
The study was initially conceived by EM and YW, and later extended by SC and XX. YW implemented the algorithm and carried out most of the computational analysis. YW, SC and XX wrote the paper. All authors read and approved the final manuscript.
Appendix
Derivation of the formula on calculating the gradient of the likelihood function
Consider the time interval [t_{s}, t_{s+1}] with full observations available at the start and the end of the interval, denoted by X(t_{s}) and X(t_{s+1}) respectively. To calculate the likelihood function L(X(t_{s}), X(t_{s+1});Θ), we discretize the time interval into N subintervals and denote the system states at these discrete points by {X^{i}i = 0,⋯,N}, where X^{0 }= X(t_{s}) and X^{N }= X(t_{s+1}) are two observations, and all other X^{i }s are intermediate states and not directly observable.
After the discretization, the likelihood function becomes, after using the Markov property of the process
For sufficiently large N, from the master equation Eq. (1.3), the conditional probability can be approximated by , where dt = (t_{s + 1 } t_{s})/N and δ_{X', X }is the Kronecker delta function.
We are interested in the gradient of the likelihood function instead of calculating the likelihood function explicitly. So we take the partial derivative of L(X(t_{s}), X(t_{s+1});Θ) w.r.t. the parameters,
Note that when N → ∞,
Therefore
where T_{k }is the time duration of the system at state k, and N_{k, k' }is the number of transitions from state k to k' occurred during the interval. Both T_{k }and N_{k, k' }are random variables, and can be viewed as the sufficient statistics of the model. E[·] represents the expectation of the random variables. For hazard functions of the form (1.5), and , thus equation (2.1) and (2.4) follow.
The formula can also be derived by the time ordered product expansion described in [55] without resorting to time discretization, as shown below. This result (Eq. A.3) suggests that the gradient of the likelihood function can be calculated by estimating the expectation of the two sufficient statistics. Note that the formula is quite general and holds for any stochastic system obeying the master equation (Eq. 1.3).
Derivation of the RJMCMC algorithm based on the time ordered product expansion of master equation
Representing the probability vector of system state as , the master equation (Eq. 1.3) can be written in a compact matrix form, using column vector notation (in this section),
where H is the time evolution matrix, of which the elements are uniquely determined by the stoichiometry matrix and the hazard functions.
The formal solution of master equation is
The probability of system evolving from a particular start state to an end state is given by the corresponding elements of the probability matrix . The time evolution matrix H is usually an infinite dimension matrix, for there are usually no upper bound for the species numbers.
The time ordered product expansion (TOPE) formula, which originates from quantum field theory, is useful to make series expansion of the matrix exponential. If we decompose the evolution matrix into two parts, H = H_{0 }+ H_{1}, the TOPE formula gives [55]
where τ_{0 }= t_{1}, τ_{1 }= t_{2 } t_{1 }etc. A proper choice is to decompose H into diagonal and offdiagonal matrices H = Ĥ  D, i.e. H_{0 }=  D and H_{1 }= Ĥ. This leads to the TOPE formula
where D represents the diagonal part (nonnegative) and Ĥ is the offdiagonal part of the matrix. The terms inside the integral, conditioned on a given Markov jump process, is the likelihood (or probability density) of the process. In case of reaction systems, a process corresponds to a set of reaction events. Thus the kth order integration gives the total probability of all reaction events with k reactions. We note that the TOPE formula provides a possible way to estimate the matrix exponential (probability matrix) by Monte Carlo integration by randomly casting reaction events and summing up the likelihood.
In the fully observed case, the likelihood function is the product of the likelihood of each subinterval,
The likelihood for each subinterval can be denoted as
where τ = t_{s+1 } t_{s}, and in the last step of the above equation we approximate the integration by a MonteCarlo integral with π(Ξ({r_{i}, t_{i}})θ_{r}) to be the likelihood of latent process Ξ({r_{i}, t_{i}}) (see Eq. 2.3) which is constrained by start/end observation (x_{a}(t_{s+1}) and x_{a}(t_{s})) and in which n_{r }is the number of reaction type r. V_{k }is the multiplication of two parts: the first part arises from the simplex integration of the time variables, which can be viewed as the measure of integration space when we convert the integration to summation; the second part is a combinatorial factor resulting from the permutation invariance of the same reaction type in a given reaction path.
Recalling that in the RJMCMC algorithm, we generate samples with different number and type of reactions via the MetropolisHasting steps. The ratio between π(·)/V_{k }of two samples gives the same acceptance probability as in Eq. (2.5).
Assuming all the reaction follows massaction law, we can derive the gradient of the likelihood function using TOPE formula. We can write Ĥ and D matrix in terms of the component of each reaction type, i.e. . Thus
We define , which is the branching ratio for reaction r in state X'. Then
Thus,
which gives the gradient formula in Eq. (2.1), (2.4) and (A.3), since the average of a frequency gives the probability.
Algorithm 1. Pseudocode of RJMCMC algorithm for fully observed case
Input observations {n(t_{s})} and generate initial path for each interval using GLPK;
Calculate the null set {q'_{k}} with the neteffect reaction matrix A;
for iter = 1: maxiteration
Randomly choose an elementary mode q_{k·};
for i = 1: number of time intervals
Randomly choose a move type p and update the reaction path in subinterval [t_{i}, t_{i+1}) according to ;
Calculate the number of each reaction type r_{m}
if min(r_{m}) == 0, AP = 0; break
else
for j = 1: J (total number of reactions within the interval)
if is negative for any a ∈ (1,2,...K),
AP = 0; break
else
Calculate the intermediate species number after the reaction:
endif
endfor
AP = min(1, AR_{ip});
endif
if AP > rand(0,1)
Accept the new path;
endif
endfor
endfor
Algorithm 2. Pseudocode of RJMCMC algorithm for partially observed case
Input observations {n_{Γ}(t_{s})} and randomly specify state for the unobserved species, generate initial path for each interval with GLPK;
Calculate the null set {} using the partial reaction matrix A_{p};
for iter = 1: maxiteration
for i = 1: number of time intervals
Randomly choose an elementary mode and a move type p; Update the reaction path in subinterval (t_{i }, t_{i+1}) according to ;
Calculate the number of each reaction type: , AP = 0; break
else
for j = 1: number of reactions within the ith interval
if is negative for any species a AP = 0; break
else
Calculate the intermediate species number after the reaction: ;
endif
endfor
endif
if x_{a, J + 1 }== x_{a}(t_{t+1})
else
Update the second interval via complementary move p';
Calculate the number of each reaction type:
if , AP = 0; break
else
for j' = 1: number of reactions within the (i + 1)th interval
if is negative for any species a AP = 0; break
else
Calculate the intermediate species number after the reaction: ;
endif
endfor
endif
Calculate for the new path and the acceptance probability ;
endif
if AP > rand(0,1)
Accept the new path;
endif
endfor
endfor
Algorithm 3: Stochastic gradient descent algorithm
Input: timecourse data
Output: set of inferred parameters {θ_{r}}
1. Initialize the reaction path using GLPK and set initial values of rate parameters;
2. Sample the latent paths for the entire observation interval with reversible jump MCMC
For fully observed case: sample latent paths for each interval s ∈ (0, m1) using Algorithm 1;
For partially observed case: sample latent paths for each paired intervals and separately for the first interval using Algorithm 2;
Calculate the gradient of each sample path according to Eq. (2.3) after burn in;
3. Calculate the gradient by averaging over sample paths;
4. Update parameter values by gradient descent step:
for all r, where η is the step size;
5. If convergence condition (to be specified) is not satisfied, return to step 2.
Acknowledgements
This work was partially supported by NSF grant DBI0846218 to XX. SC was supported by NIH grants R01GM75309 and P50GM76516. EM was supported by grants NSF EF0330786 and NIH R01GMO86883. We thank Jacob Biesinger, Qing Nie and GuiBo Ye for helpful discussions.
References

Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theories to phenotypes.
Nat Rev Genet 2005, 6:451464. PubMed Abstract  Publisher Full Text

Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell.
Science 2002, 297:11831186. PubMed Abstract  Publisher Full Text

Swain PS, Elowitz MB, Siggia ED: Intrinsic and extrinsic contributions to stochasticity in gene expression.
Proc Natl Acad Sci USA 2002, 99:1279512800. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Blake WJ, Kaern M, Cantor CR, Collins JJ: Noise in eukaryotic gene expression.
Nature 2003, 422:633637. PubMed Abstract  Publisher Full Text

Pirone JR, Elston TC: Fluctuations in transcription factor binding can explain the graded and binary responses observed in inducible gene expression.
J Theor Biol 2004, 226:111121. PubMed Abstract  Publisher Full Text

Thattai M, van Oudenaarden A: Intrinsic noise in gene regulatory networks.
Proc Natl Acad Sci USA 2001, 98:86148619. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

McAdams HH, Arkin A: Stochastic mechanisms in gene expression.
Proc Natl Acad Sci USA 1997, 94:814819. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ozbudak EM, Thattai M, Kurtser I, Grossman AD, van Oudenaarden A: Regulation of noise in the expression of a single gene.
Nat Genet 2002, 31:6973. PubMed Abstract  Publisher Full Text

BarEven A, Paulsson J, Maheshri N, Carmi M, O'Shea E, Pilpel Y, Barkai N: Noise in protein expression scales with natural protein abundance.
Nat Genet 2006, 38:636643. PubMed Abstract  Publisher Full Text

Pedraza JM, van Oudenaarden A: Noise propagation in gene networks.
Science 2005, 307:19651969. PubMed Abstract  Publisher Full Text

Choi PJ, Cai L, Frieda K, Xie XS: A stochastic singlemolecule event triggers phenotype switching of a bacterial cell.
Science 2008, 322:442446. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathway bifurcation in phage lambdainfected Escherichia coli cells.
Genetics 1998, 149:16331648. PubMed Abstract  PubMed Central Full Text

Fraser HB, Hirsh AE, Giaever G, Kumm J, Eisen MB: Noise minimization in eukaryotic gene expression.
PLoS Biol 2004, 2:e137. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Becskei A, Serrano L: Engineering stability in gene networks by autoregulation.
Nature 2000, 405:590593. PubMed Abstract  Publisher Full Text

Rao CV, Wolf DM, Arkin AP: Control, exploitation and tolerance of intracellular noise.
Nature 2002, 420:231237. PubMed Abstract  Publisher Full Text

Thattai M, van Oudenaarden A: Stochastic gene expression in fluctuating environments.
Genetics 2004, 167:523530. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schultz D, Ben Jacob E, Onuchic JN, Wolynes PG: Molecular level stochastic model for competence cycles in Bacillus subtilis.
Proc Natl Acad Sci USA 2007, 104:1758217587. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Beaumont HJE, Gallie J, Kost C, Ferguson GC, Rainey PB: Experimental evolution of bet hedging.
Nature 2009, 462:9093. PubMed Abstract  Publisher Full Text

Lawrence ND, Girolami M, Rattray M, Sanguinetti G, (Eds): Learning and Inference in Computational Systems Biology. Cambridge, MA, The MIT Press; 2010.

Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods.
Genome Res 2003, 13:24672474. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kirkpatrick S, Gelatt C, Vecchi M: Optimization by Simulated Annealing.
Science 1983, 220:671680. PubMed Abstract  Publisher Full Text

Janssens H, Hou S, Jaeger J, Kim AR, Myasnikova E, Sharp D, Reinitz J: Quantitative and predictive model of transcriptional control of the Drosophila melanogaster even skipped gene.
Nat Genet 2006, 38:11591165. PubMed Abstract  Publisher Full Text

Chen WW, Schoeberl B, Jasper PJ, Niepel M, Nielsen UB, Lauffenburger DA, Sorger PK: Inputoutput behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data.
Mol Syst Biol 2009, 5:239. PubMed Abstract  PubMed Central Full Text

Busch H, CamachoTrullio D, Rogon Z, Breuhahn K, Angel P, Eils R, Szabowski A: Gene network dynamics controlling keratinocyte migration.
Mol Syst Biol 2008, 4:199. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gillespie D: Exact stochastic simulation of coupled chemical reactions.
J Phys ChemUs 1977, 81:23402361. Publisher Full Text

Gillespie DT: Stochastic simulation of chemical kinetics.
Annual review of physical chemistry 2007, 58:3555. PubMed Abstract  Publisher Full Text

Li H, Cao Y, Petzold LR, Gillespie DT: Algorithms and software for stochastic simulation of biochemical reacting systems.
Biotechnol Prog 2008, 24:5661. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Henderson DA, Boys RJ, Krishnan KJ, Lawless C, Wilkinson DJ: Bayesian Emulation and Calibration of a Stochastic Computer Model of Mitochondrial DNA Deletions in Substantia Nigra Neurons.
J Am Stat Assoc 2009, 104:7687. Publisher Full Text

Gillespie D: The chemical Langevin equation.
J Chem Phys 2000, 113:297306. Publisher Full Text

Golightly A, Wilkinson DJ: Bayesian inference for stochastic kinetic models using a diffusion approximation.
Biometrics 2005, 61:781788. PubMed Abstract  Publisher Full Text

Golightly A, Wilkinson DJ: Bayesian sequential inference for stochastic kinetic biochemical network models.
J Comput Biol 2006, 13:838851. PubMed Abstract  Publisher Full Text

Boys RJ, Wilkinson DJ, Kirkwood TBL: Bayesian inference for a discretely observed stochastic kinetic model.
Stat Comput 2008, 18:125135. Publisher Full Text

Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems.
Nat Rev Genet 2009, 10:122133. PubMed Abstract  Publisher Full Text

Tian T, Xu S, Gao J, Burrage K: Simulated maximum likelihood method for estimating kinetic rates in gene expression.
Bioinformatics 2007, 23:8491. PubMed Abstract  Publisher Full Text

Sisson SA, Fan Y, Tanaka MM: Sequential Monte Carlo without likelihoods.
Proc Natl Acad Sci USA 2007, 104:17601765. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MPH: Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems.
Journal of The Royal Society Interface 2009, 6:187202. Publisher Full Text

Yosiphon G, Mjolsness E: Towards the inference of stochastic biochemical network and parameterized grammar models. In Learning and Inference in Computational Systems Biology. Edited by Lawrence ND, Girolami M, Rattray M, Sanguinetti G. MIT Press; 2009:297314.

Reinker S, Altman RM, Timmer J: Parameter estimation in stochastic biochemical reactions.
Systems biology 2006, 153:168178. PubMed Abstract

Henderson D, Boys R, Wilkinson D: Bayesian Calibration of a Stochastic Kinetic Computer Model Using Multiple Data Sources.
Biometrics 2009, 66:24956. PubMed Abstract  Publisher Full Text

Green P: Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.
Biometrika 1995, 82:711732. Publisher Full Text

Brooks S, Giudici P, Roberts G: Efficient construction of reversible jump Markov chain Monte Carlo proposal distributions.
J Roy Stat Soc B 2003, 65:339. Publisher Full Text

Urbanczik R, Wagner C: An improved algorithm for stoichiometric network analysis: theory and applications.
Bioinformatics 2005, 21:12031210. PubMed Abstract  Publisher Full Text

von Kamp A, Schuster S: Metatool 5.0: fast and flexible elementary modes analysis.
Bioinformatics 2006, 22:19301931. PubMed Abstract  Publisher Full Text

Wright J, Wagner A: Exhaustive identification of steady state cycles in large stoichiometric networks.
BMC systems biology 2008, 2:61. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

EdelsteinKeshet L: Mathematical Models in Biology. Society for Industrial and Applied Mathematics; 2005.

GLPK  GNU Linear Programming Toolkit [http://www.gnu.org/software/glpk/] webcite

GLPKMEX  a Matlab MEX interface for the GLPK library [http://sourceforge.net/projects/glpkmex/] webcite

Schuster S, Fell DA, Dandekar T: A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks.
Nat Biotechnol 2000, 18:326332. PubMed Abstract  Publisher Full Text

Rabitz H: Systems Analysis at the Molecular Scale.
Science 1989, 246:221226. PubMed Abstract  Publisher Full Text

Saltelli A, Ratto M, Tarantola S, Campolongo F: Sensitivity analysis for chemical models.
Chem Rev 2005, 105:28112828. PubMed Abstract  Publisher Full Text

Rosenfeld N, Perkins TJ, Alon U, Elowitz MB, Swain PS: A fluctuation method to quantify in vivo fluorescence data.
Biophys J 2006, 91:759766. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rosenfeld N, Young JW, Alon U, Swain PS, Elowitz MB: Gene regulation at the singlecell level.
Science 2005, 307:19621965. PubMed Abstract  Publisher Full Text

Clayton E, Doupe DP, Klein AM, Winton DJ, Simons BD, Jones PH: A single type of progenitor cell maintains normal epidermis.
Nature 2007, 446:185189. PubMed Abstract  Publisher Full Text

Mjolsness E, Yosiphon G: Stochastic process semantics for dynamical grammars.
Annals of Mathematics and Artificial Intelligence 2006, 47:329395. Publisher Full Text