Global parameter identification of stochastic
reaction networks from single trajectories

arXiv:1111.4785v1 [q-bio.MN] 21 Nov 2011

Christian L. Müller† , Rajesh Ramaswamy† , and Ivo F. Sbalzarini

Abstract We consider the problem of inferring the unknown parameters of a
stochastic biochemical network model from a single measured time-course of the
concentration of some of the involved species. Such measurements are available,
e.g., from live-cell fluorescence microscopy in image-based systems biology. In
addition, fluctuation time-courses from, e.g., fluorescence correlation spectroscopy
provide additional information about the system dynamics that can be used to more
robustly infer parameters than when considering only mean concentrations. Estimating model parameters from a single experimental trajectory enables single-cell
measurements and quantification of cell–cell variability. We propose a novel combination of an adaptive Monte Carlo sampler, called Gaussian Adaptation, and efficient exact stochastic simulation algorithms that allows parameter identification
from single stochastic trajectories. We benchmark the proposed method on a linear
and a non-linear reaction network at steady state and during transient phases. In addition, we demonstrate that the present method also provides an ellipsoidal volume
estimate of the viable part of parameter space and is able to estimate the physical
volume of the compartment in which the observed reactions take place.

1 Introduction
Systems biology implies a holistic research paradigm, complementing the reductionist approach to biological organization [16, 15]. This frequently has the goal of
mechanistically understanding the function of biological entities and processes in
interaction with the other entities and processes they are linked to or communicate
with. A formalism to express these links and connections is provided by network
Christian L. Müller, Rajesh Ramaswamy, and Ivo F. Sbalzarini:
Institute of Theoretical Computer Science and Swiss Institute of Bioinformatics, ETH Zurich, CH–
8092 Zurich, Switzerland, e-mail: christian.mueller@inf.ethz.ch, rajeshr@ethz.ch, ivos@ethz.ch.
† These authors contributed equally to this work.

1

2

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

models of biological processes [4, 1]. Using concepts from graph theory [26] and
dynamic systems theory [44], the organization, dynamics, and plasticity of these
networks can then be studied.
Systems biology models of molecular reaction networks contain a number of
parameters. These are the rate constants of the involved reactions and, if spatiotemporal processes are considered, the transport rates, e.g. diffusion constants, of the
chemical species. In order for the models to be predictive, these parameters need to
be inferred. The process of inferring them from experimental data is called parameter identification. If in addition also the network structure is to be inferred from
data, the problem is called systems identification. Here, we consider the problem of
identifying the parameters of a biochemical reaction network from a single, noisy
measurement of the concentration time-course of some of the involved species.
While this time series can be long, ensemble replicas are not possible, either because the measurements are destructive or one is interested in variations between
different specimens or cells. This is particularly important in molecular systems biology, where cell–cell variations are of interest or large numbers of experimental
replica are otherwise not feasible [45].
This problem is particularly challenging and traditional genomic and proteomic
techniques do not provide single-cell resolution. Moreover, in individual cells the
molecules and chemical reactions can only be observed indirectly. Frequently, fluorescence microscopy is used to observe biochemical processes in single cells. Fluorescently tagging some of the species in the network of interest allows measuring the
spatiotemporal evolution of their concentrations from video microscopy and fluorescence photometry. In addition, fluorescence correlation spectroscopy (FCS) allows
measuring fluctuation time-courses of molecule numbers [23].
Using only a single trajectory of the mean concentrations would hardly allow
identification of network parameters. There could be several combinations of network parameters that lead to the same mean dynamics. A stochastic network model,
however, additionally provides information about the fluctuations of the molecular
abundances. The hope is that there is then only a small region of parameter space that
produces the correct behavior of the mean and the correct spectrum of fluctuations
[31]. Experimentally, fluctuation spectra can be measured at single-cell resolution
using FCS.
The stochastic behavior of biochemical reaction networks can be due to low copy
numbers of the reacting molecules [39, 10]. In addition, biochemical networks may
exhibit stochasticity due to extrinsic noise. This can persist even at the continuum
scale, leading to continuous–stochastic models. Extrinsic noise can, e.g., arise from
environmental variations or variations in how the reactants are delivered into the system. Also measurement uncertainties can be accounted for in the model as extrinsic
noise, modeling our inability to precisely quantify the experimental observables.
We model stochastic chemical kinetics using the chemical master equation
(CME). Using a CME forward model in biological parameter identification amounts
to tracking the evolution of a probability distribution, rather than just a single
function. This prohibits predicting the state of the system and only allows statements about the probability for the system to be in a certain state, hence requiring

Global stochastic parameter identification

3

sampling-based parameter identification methods. In the stochastic–discrete context, a number of different approaches have been suggested. Boys et al. proposed
a fully Bayesian approach for parameter estimation using an explicit likelihood for
data/model comparison and a Markov Chain Monte Carlo (MCMC) scheme for
sampling [5]. Zechner et al. developed a recursive Bayesian estimation technique
[46, 45] to cope with cell–cell variability in experimental ensembles. Toni and coworkers used an approximate Bayesian computation (ABC) ansatz, as introduced
by Marjoram and co-workers [25], that does not require an explicit likelihood [43].
Instead, sampling is done in a sequential Monte Carlo (or particle filter) framework.
Reinker et al. used a hidden Markov model where the hidden states are the actual
molecule abundances and state transitions model chemical reactions [40]. Inspired
by Prediction Error Methods [24], Cinquemani et al. identified the parameters of a
hybrid deterministic–stochastic model of gene expression from multiple experimental time courses [7]. Randomized optimization algorithms have been used, e.g., by
Koutroumpas et al. who applied a Genetic Algorithm to a hybrid deterministic–
stochastic network model [21]. More recently, Poovathingal and Gunawan used
another global optimization heuristic, the Differential Evolution algorithm [32]. A
variational approach for stochastic two-state systems is proposed by Stock and coworkers based on Maximum Caliber [41], an extension of Jaynes’ Maximum Entropy principle [14] to non-equilibrium systems. If estimates are to be made based
on a single trajectory, the stochasticity of the measurements and of the model leads
to very noisy similarity measures, requiring optimization and sampling schemes that
are robust against noise in the data.
Here, we propose a novel combination of exact stochastic simulations for a
CME forward model and an adaptive Monte Carlo sampling technique, called
Gaussian Adaptation, to address the single-trajectory parameter estimation problem
for monostable stochastic biochemical reaction networks. Evaluations of the CME
model are done using exact partial-propensity stochastic simulation algorithms [35].
Parameter optimization uses Gaussian Adaptation. The method iteratively samples
model parameters from a multivariate normal distribution and evaluates a suitable
objective function that measures the distance between the dynamics of the forward
model output and the experimental measurements. In addition to estimates of the
kinetic parameters in the network, the present method also provides an ellipsoidal
volume estimate of the viable part of parameter space and is able to estimate the
physical volume of the intra-cellular compartment in which the reactions take place.
We assume that quantitative experimental time series of either a transient or the
steady state of the concentrations of some of the molecular species in the network
are available. This can, for example, be obtained from single-cell fluorescence microscopy by translating fluorescence intensities to estimated chemical concentrations. Accurate methods that account for the microscope’s point-spread function
and the camera noise model are available to this end [12, 13, 6]. Additionally, FCS
spectra can be analyzed in order to quantify molecule populations, their intrinsic
fluctuations, and lifetimes [23, 34, 39]. The present approach requires only a single
stochastic trajectory from each cell. Since the forward model is stochastic and only a
single experimental trajectory is used, the objective function needs to robustly mea-

4

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

sure closeness between the experimental and the simulated trajectories. We review
previously considered measures and present a new distance function in Sec. 5. First,
however, we set out the formal stochastic framework and problem description below.
We then describe Gaussian Adaptation and its applicability to the current estimation
task. The evaluation of the forward model is outlined in Sec. 4. We consider a cyclic
linear chain as well as a non-linear colloidal aggregation model as benchmark test
cases in Sec. 6 and conclude in Sec. 7.

2 Background and problem statement
We consider a network model of a biochemical system given by M coupled chemical
reactions
N

kj

N

∑ νi,−j Si −−−−−→ ∑ νi,+j Si

i=1

∀ j = 1, . . . , M

(1)

i=1

between N species, where ν − = [νi,−j ] and ν + = [νi,+j ] are the stoichiometry matrices
of the reactants and products, respectively, and Si is the ith species in the reaction
network. Let ni be the population (molecular copy number) of species Si . The reactions occur in a physical volume Ω and the macroscopic reaction rate of reaction j
is k j . This defines a dynamic system with state n(t) = [ni (t)] and M + 1 parameters
θ = [k1 , . . . , kM , Ω ].
The state of such a system can be interpreted as a realization of a random variable
n(t) that changes over time t. All one can know about the system is the probability
for it to be in a certain state at a certain time t j given the system’s state history, hence
P(n(t j ) | n(t j−1 ), . . . , n(t1 ), n(t0 )) dNn
= Prob{n(t j ) ∈ [n(t j ), n(t j ) + dn) | n(ti ), i = 0, . . . , j − 1 } .

(2)

A frequently made model assumption, substantiated by physical reasoning, is
that the probability of the current state depends solely on the previous state, i.e.,
P(n(t j ) | n(t j−1 ), . . . , n(t1 ), n(t0 )) = P(n(t j ) | n(t j−1 )) .

(3)

The system is then modeled as a first-order Markov chain where the state n evolves
as
n(t + ∆t) = n(t) + Ξ (∆t; n,t)
(4)
This is the equation of motion of the system. If n is real-valued, it defines a
continuous–stochastic model in the form of a continuous-state Markov chain. Discrete n, as is the case in chemical kinetics, amount to discrete–stochastic models expressed as discrete-state Markov chains. The Markov propagator Ξ is itself a random variable, distributed with probability distribution Π (ξ | ∆t; n,t) =
P(n + ξ ,t + ∆t | n,t) for the state change ξ . For continuous-state Markov chains,
Π is a continuous probability density function (PDF), for discrete-state Markov

Global stochastic parameter identification

5

chains a discrete probability distribution. If Π (ξ ) = δ (ξ − ξ 0 ), with δ the Dirac
delta distribution, then the system’s state evolution becomes deterministic with predictable discrete or continuous increments ξ 0 . Deterministic models can hence be
interpreted as a limit case of stochastic models [22].
In chemical kinetics, the probability distribution Π of the Markov propagator
is a linear combination of Poisson distributions with weights given by the reaction
stoichiometry. This leads to the equation of motion for the population n given by


ψ1
n(t + ∆t) = n(t) + (ν + − ν − )  . . .  ,
(5)
ψM
where ψi ∼ P(ai (n(t))∆t) is a random variable from the Poisson distribution with
rate λ = ai (n(t))∆t. The second term on the right-hand side of Eq. 5 follows a
probability distribution Π (ξ | ∆t; n,t) whose explicit form is analytically intractable
in the general case. The rates a j , j = 1, . . . M, are called the reaction propensities
and are defined as

N 
kj
ni
aj = ∏
.
(6)
−
ν
1+
ν−
∑N
i, j
i0 =1 i0 , j
i=1
Ω
They depend on the macroscopic reaction rates and the reaction volume and can be
interpreted as the probability rates of the respective reactions. Advancing Eq. 5 with
a ∆t such that more than one reaction event happens per time step yields an approximate simulation of the biochemical network as done in approximate stochastic
simulation algorithms [9, 3].
An alternative approach consists in considering the evolution of the state probability distribution P(n,t | n0 ,t0 ) of the Markov chain described by Eq. 5, hence:
!
M
N ν − −ν +
∂P
i, j
i, j
(7)
= ∑ ∏ Ei E j − 1 a j (n(t))P(n,t)
∂t
j=1 i=1
with the step operator Eip f (n) = f (n + pî ) for any function f where î is the Ndimensional unit vector along the ith dimension. This equation is called the chemical
master equation (CME). Directly solving it for P is analytically intractable, but
trajectories of the Markov chain governed by the unknown state probability P can
be sampled using exact stochastic simulation algorithms (SSA) [8]. Exact SSAs
are exact in the sense that they sample Markov chain realizations from the exact
solution P of the CME, without ever explicitly computing this solution. Since SSAs
are Monte Carlo algorithms, however, a sampling error remains.
Assuming that the population n increases with the volume Ω , n can be approximated as a continuous random variable in the limit of large volumes, and Eq. 5
becomes


η1
n(t + ∆t) = n(t) + (ν + − ν − )  . . .  ,
(8)
ηM

6

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

where ηi ∼ N (ai (n(t))∆t, ai (n(t))∆t) are normally distributed random variables.
The second term on the right-hand side of this equation is a random variable that is
distributed according to the corresponding Markov propagator Π (ξ | ∆t; n,t), which
is a Gaussian. Equation 8 is called the chemical Langevin equation with Π given by
1

T Σ −1 (ξ −µ)

Π (ξ | ∆t; n,t) = (2π)−N/2 |Σ |−1/2 e− 2 (ξ −µ)

,

(9)

where



a1 (n(t))
µ = ∆t (ν + − ν − )  . . .  and Σ = ∆t (ν + − ν − ) diag (n(t))(ν + − ν − )T .
aM (n(t))
The corresponding equation for the evolution of the state PDF is the non-linear
Fokker-Planck equation, given by


∂P
1
T
=∇
D ∇ − F P(n,t) ,
(10)
∂t
2
where



∂
∂
,...,
∇ =
,
∂ n1
∂ nN
T

1
∆t→0 ∆t

Fi = lim

Z +∞
−∞

dξi ξi Π (ξ | ∆t; n,t) ,

(11)

(12)

and
1
Di j = lim
∆t→0 ∆t

Z +∞ Z +∞
−∞

−∞

dξi dξ j ξi ξ j Π (ξ | ∆t; n,t) .

(13)

At much larger Ω , when the population n is on the order of Avogadro’s number,
Eq. 8 can be further approximated as


φ1 (n(t))∆t
+
− 
,
...
n(t + ∆t) = n(t) + (ν − ν )
(14)
φM (n(t))∆t
1−∑N

ν−

ν−

i, j
i0 =1 i0, j
where φ j (n) = k j Ω
∏Ni=1 ni (νi,−j !)−1 . Note that the second term on the
right-hand side of this equation is a random variable whose probability distribution
is the Dirac delta



φ1 (n(t))∆t
 .
...
Π (ξ | ∆t; n,t) = δ ξ − (ν + − ν − ) 
(15)
φM (n(t))∆t

Global stochastic parameter identification

7

Equation 14 hence is a deterministic equation of motion. In the limit ∆t → 0 this
equation can be written as the ordinary differential equation


φ1 (x(t))
dx
= (ν + − ν − )  . . . 
(16)
dt
φM (x(t))
for the concentration x = nΩ −1 . This is the classical reaction rate equation for the
system in Eq. 1.
By choosing the appropriate probability distribution Π of the Markov propagator,
one can model reaction networks in different regimes: small population n (small Ω )
using SSA over Eq. 7, intermediate population (intermediate Ω ) using Eq. 8, and
large population (large Ω ) using Eq. 16. The complete model definition therefore is
M (θ ) = {ν − , ν + , Π }.
The problem considered here can then be formalized as follows: Given a forward
model M (θ ) and a single noisy trajectory of the population of the chemical species
n̂(t0 + (q − 1)∆texp ) at K discrete time points t = t0 + (q − 1)∆texp , q = 1, . . . , K, we
wish to infer θ = [k1 , . . . , kM , Ω ]. The time between two consecutive measurements
∆texp and the number of measurements K are given by the experimental technique
used. As a forward model we use the full CME as given in Eq. 7 and sample trajectories from it using the partial-propensity formulation of Gillespie’s exact SSA as
described in Sec. 4.

3 Gaussian Adaptation for global parameter optimization,
approximate Bayesian computation, and volume estimation
Gaussian Adaptation (GaA), introduced in the late 1960’s by Gregor Kjellström
[17, 19], is a Monte Carlo technique that has originally been developed to solve
design-centering and optimization problems in analog electric circuit design. Design centering solves the problem of determining the nominal values (resistances,
capacitances, etc.) of the components of a circuit such that the circuit output is
within specified design bounds and is maximally robust against random variations
in the circuit components with respect to a suitable criterion or objective function.
This problem is a superset of general optimization, where one is interested in finding a parameter vector that minimizes (or maximizes) the objective function without
any additional robustness criterion. GaA has been specifically designed for scenarios where the objective function f (θ ) is only available in a black-box (or oracle)
model that is defined on a real-valued domain A ⊆ Rn and returns real-valued output R. The black-box model assumes that gradients or higher-order derivatives of
the objective function may not exist or may not be available, hence including the
class of discontinuous and noisy functions. The specific black-box function used
here is presented in Sec. 5.

8

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

The principle idea behind GaA is the following: Starting from a user-defined
point in parameter space, GaA explores the space by iteratively sampling single
parameter vectors from a multivariate Gaussian distribution N (m, Σ ) whose mean
m ∈ Rn and covariance matrix Σ ∈ Rn×n are dynamically adapted based on the information from previously accepted samples. The acceptance criterion depends on the
specific mode of operation, i.e., whether GaA is used as an optimizer or as a sampler [28, 27]. Adaptation is performed such as to maximize the entropy of the search
distribution under the constraint that acceptable search points are found with a predefined, fixed hitting (success) probability p < 1 [19]. Using the

qdefinition of the
(2πe)n det(Σ )
entropy of a multivariate Gaussian distribution H (N ) = log
shows that this is equivalent to maximizing the determinant of the covariance matrix Σ . GaA thus follows Jaynes’ Maximum Entropy principle [14].
GaA starts by setting the mean m(0) of the multivariate Gaussian to an initial
acceptable point θ (0) and the Cholesky factor Q(0) of the covariance matrix to the
identity matrix I. At each iteration g > 0, the covariance Σ (g) is decomposed as:


T

T

Σ (g) = r · Q(g) r · Q(g) = r2 Q(g) Q(g) , where r is the scalar step size
that controls the scale of the search. The matrix Q(g) is the normalized square root
of Σ (g) , found by eigen- or Cholesky decomposition of Σ (g) . The candidate parameter vector in iteration g + 1 is sampled from a multivariate Gaussian according to
θ (g+1) = m(g) + r(g) Q(g) η (g) , where η (g) ∼ N (0, I). The parameter vector is then
evaluated by the objective function f (θ (g+1) ).
Only if the parameter vector is accepted, the following adaptation rules are applied: The step size r is increased as r(g+1) = fe · r(g) , where fe > 1 is termed the
expansion factor. The mean of the proposal distribution is updated as


1
1 (g+1)
(g+1)
m
m(g) +
= 1−
θ
.
(17)
Nm
Nm
Nm is a weighting factor that
rate of the method. The successful
 controls the learning

(g+1)
(g+1)
(g)
= θ
−m
is used to perform a rank-one update of
search direction d


the covariance matrix: Σ (g+1) = 1 − N1C Σ (g) + N1C d (g+1) d (g+1) T . NC weights the
influence of the accepted parameter vector on the covariance matrix. In order to
decouple the volume of the covariance (controlled by r(g+1) ) from its orientation,
Q(g+1) is normalized such that det(Q(g+1) ) = 1.
In case θ (g+1) is not accepted at the current iteration, only the step size is adapted
as r(g+1) = fc · r(g) , where fc < 1 is the contraction factor.
The behavior of GaA is controlled by several strategy parameters. Kjellström
analyzed the information-theoretic optimality of the acceptance probability p for
GaA in general regions [19]. He concluded that the efficiency E of the process and
p are related as E ∝ −p log p, leading to an optimal p = 1e ≈ 0.3679, where e is
Euler’s number. A proof is provided in [18]. Maintaining this optimal hitting prob-

Global stochastic parameter identification

9

ability corresponds to leaving the volume of the distribution, measured by det(Σ ),
constant under stationary conditions. Since det(Σ ) = r2n det(Q QT ), the expansion
and contraction factors fe and fc expand or contract the volume by a factor of fe2n
and fc2n , respectively. After S accepted and F rejected samples, a necessary condiS
, and
tion for constant volume thus is: ∏Si=1 ( fe )2n ∏Fi=1 ( fc )2n = 1. Using p = S+F
introducing a small β > 0, the choice fe = 1 + β (1 − p) and fc = 1 − β p satisfies
the constant-volume condition to first order. The scalar rate β is coupled to NC . NC
influences the update of Σ ∈ Rn×n , which contains n2 entries. Hence, NC should be
related to n2 . We suggested using NC = (n + 1)2 / log(n + 1) as a standard value, and
coupling β = N1C [29]. A similar reasoning is also applied to Nm . Since Nm influences the update of m ∈ Rn , it is reasonable to set Nm ∝ n. We propose Nm = en as
a standard value.
Depending on the specific acceptance rule used, GaA can be turned into a global
optimizer [29], an adaptive MCMC sampler [28, 27], or a volume estimation method
[30], as described next.

3.1 GaA for global black-box optimization
In a minimization scenario, GaA uses an adaptive-threshold acceptance mechanism.
(0)
Given an initial scalar cutoff threshold cT , we accept a parameter vector θ (g+1) at
(g)
iteration 
g + 1 if f (θ (g+1) ) < cT . Upon acceptance, the threshold cT is lowered as
(g+1)

(g)

cT
= 1 − N1T cT + N1T f (θ (g+1) ), where NT controls the weighting between the
old threshold and the objective-function value of the accepted sample. This sampledependent threshold update renders the algorithm invariant to linear transformations
of the objective function. The standard strategy parameter value is NT = en [28].
We refer to [28] for further information about convergence criteria and constraint
handling techniques in GaA.

3.2 GaA for approximate Bayesian computation and viable volume
estimation
Replacing the threshold acceptance-criterion by a probabilistic Metropolis criterion,
and setting Nm = 1, turns GaA into an adaptive MCMC sampler with global adaptive scaling [2]. We termed this method Metropolis-GaA [28, 27]. Its strength is that
GaA can automatically adapt to the covariance of the target probability distribution
while maintaining the fixed hitting probability. For standard MCMC, this cannot be
achieved without fine-tuning the proposal using multiple MCMC runs. We hypothesize that GaA might also be an effective tool for approximate Bayesian computation
(ABC) [43]. In essence, the ABC ansatz is MCMC without an explicit likelihood

10

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

function [25]. The likelihood is replaced by a distance function — which plays the
same role as our objective function — that measures closeness between a parameterized model simulation and empirical data D, or summary statistics thereof. When a
uniform prior over the parameters and a symmetric proposal are assumed, a parameter vector in ABC is unconditionally accepted if its corresponding distance function
value f (θ (g+1) ) < cT [25]. The threshold cT is a problem-dependent constant that
is fixed prior to the actual computation. Marjoram and co-workers have shown that
samples obtained in this manner are approximately drawn from the posterior parameter distribution given the data D. While Pritchard et al. used a simple rejection
sampler [33], Marjoram and co-workers proposed a standard MCMC scheme [25].
Toni and co-workers used sequential MC for sample generation [43]. To the best of
our knowledge, however, the present work presents the first application of an adaptive MCMC scheme for ABC in biochemical network parameter inference. Finally,
we emphasize that when GaA’s mean, covariance matrix, and hitting probability p
stabilize during ABC, they provide direct access to an ellipsoidal estimation of the
volume of the viable parameter space as defined by the threshold cT [30]. Hafner
and co-workers have shown how to use such viable volume estimates for model
discrimination [11].

4 Evaluation of the forward model
In each iteration of the GaA algorithm, the forward model of the network needs
to be evaluated for the proposed parameter vector θ . This requires an efficient and
exact SSA for the chemical kinetics of the reaction network, used to generate trajectories n(t) from M (θ ). Since GaA could well propose parameter vectors that lead
to low copy numbers for some species, it is important that the SSA be exact since
approximate algorithms are not appropriate at low copy number.
In its original formulation, Gillespie’s SSA has a computational cost that is
linearly proportional to the total number M of reactions in the network. If many
model evaluations are required, as in the present application, this computational cost
quickly becomes prohibitive. While more efficient formulations of SSA have been
developed for weakly coupled reaction networks, their computational cost remains
proportional to M for strongly coupled reaction networks [35]. A reaction network is
weakly coupled if the number of reactions that are influenced by any other reaction
is bounded by a constant. If a network contains at least one reaction whose firing
influences the propensities of a fixed proportion (in the worst case all) of the other
reactions, then the network is strongly coupled [35]. Scale-free networks as seem to
be characteristic for systems biology models [1, 42] are by definition strongly coupled. This is due to the existence of hubs that have a higher connection probability
than other nodes. These hubs frequently correspond to chemical reactions that produce or consume species that also participate in the majority of the other reactions,
such as water, ATP, or CO2 in metabolic networks.

Global stochastic parameter identification

11

We use partial-propensity methods [35, 36] to simulate trajectories according to
the solution of the chemical master Eq. 7 of the forward model. Partial-propensity
methods are exact SSAs whose computational cost scales at most linearly with the
number N of species in the network [35]. For large networks, this number is usually much smaller than the number of reactions. Depending on the network model
at hand, different partial-propensity methods are available for its efficient simulation. Strongly coupled networks where the rate constants span only a limited spectrum of values are best simulated with the partial-propensity direct method (PDM)
[35]. Multi-scale networks where the rate constants span many orders of magnitude are most efficiently simulated using the sorting partial-propensity direct method
(SPDM) [35]. Weakly coupled reaction networks can be simulated at constant computational cost using the partial-propensity SSA with composition-rejection sampling (PSSA-CR) [37]. Lastly, reaction networks that include time delays can be
exactly simulated using the delay partial-propensity direct method (dPDM) [38].
Different combinations of the algorithmic modules of partial-propensity methods
can be used to constitute all members of this family of SSAs [36]. We refer to the
original publications for algorithmic details, benchmarks of the computational cost,
and a proof of exactness of partial-propensity methods.

5 Objective function
In the context of parameter identification of stochastic biochemical networks, a
number of distance or objective functions have previously been suggested. Reinker
et al. proposed an approximate maximum-likelihood measure under the assumption
that only a small number of reactions fire between two experimental measurement
points, and a likelihood based on singular value decomposition that works when
many reactions occur per time interval [40]. Koutroumpas et al. compared objective functions based on least squares, normalized cross-correlations, and conditional
probabilities using a Genetic Algorithm [21]. Koeppl and co-workers proposed the
Kantorovich distance to compare experimental and model-based probability distributions [20]. Alternative distance measures include the Earth Mover’s distance or
the Kolomogorov-Smirnov distance [32]. These distance measures, however, can
only be used when many experimental trajectories are available. In order to measure the distance between a single experimental trajectory n̂(t) and a single model
output n(t), we propose a novel cost function f (θ ) = f (M (θ ), n̂) that reasonably
captures the kinetics of a monostable system. We define a compound objective function f (θ ) = f1 (θ ) + f2 (θ ) with
4

f1 (θ ) = ∑ γi ,
i=1

where

N

z

x
|ACFl (n̂i ) − ACFl (ni )|
∑l=0
,
zx
ACFl (n̂i )
∑l=0
i=1

f2 (θ ) = ∑

(18)

12

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini
N

γi =

s

∑

j=1

µi (n j ) − µi (n̂ j )
µi (n̂ j )

2
(19)

with the central moments given by
(
if i = 1
∑Kp=1 n j (t0 + (p − 1)∆texp )
i
µi (n j ) =
(| ∑Kq=1 n j (t0 + (q − 1)∆texp ) − µ1 (n j ) |)1/i otherwise

(20)

and the time-autocorrelation function (ACF) at lag l given by
ACFl (ni ) =

ni (t0 )ni (t0 + l ∆texp ) − (µ1 (ni ))2
.
µ2 (ni )

The variable zx is the lag at which the experimental ACF crosses 0 for the first time.
The function f1 (θ ) measures the difference between the first four moments of n and
n̂. This function alone would, however, not be enough to capture the kinetics since
it lacks information about correlations in time. This is taken into account by f2 (θ ),
measuring the difference in the lifetimes of all chemical species. These lifetimes are
systematically modulated by the volume Ω [39], hence enabling volumetric measurements of intra-cellular reaction compartments along with the identification of
the rate constants.
The present objective function allows inclusion of experimental readouts from
image-based systems biology. The moment-matching part is a typical readout from
fluorescence photometry, whereas the autocorrelation of the fluctuations can directly
be measured using, e.g., FCS.

a

b

c

100

300

300

250

250

80

n̂1 200

n̂1

200

n̂1

60

n̂2

n̂2

150

n̂3

40

n̂2

150
100

20
0
0

100

50
20

40

60

time [s]

80

100

0
0

50
20

40

60

time [s]

80

100

0
0

20

40

60

80

100

time [s]

Fig. 1 In silico data for all test cases. a. Time evolution of the populations of three species in the
cyclic chain model at steady state (starting at t0 = 2000). b. Time evolution of the populations of
two species in the aggregation model at steady state (starting at t0 = 5000). c. Same as b, but during
the transient phase (starting at t0 = 0).

Global stochastic parameter identification

13

6 Results
We estimate the unknown parameters θ for two reaction networks: a weakly coupled
cyclic chain and a strongly coupled non-linear colloidal aggregation network. For
the cyclic chain we estimate θ at steady state. For the aggregation model we estimate
θ both at steady state and in the transient phase. Every kinetic parameter is allowed
to vary in the interval [10−3 , 103 ] and the reaction volume Ω in [1, 500]. Each GaA
run starts from a point selected uniformly at random in logarithmic parameter space.

6.1 Weakly coupled reaction network: cyclic chain
The cyclic chain network is given by:
k

i
Si −
→
Si+1

kN

Si −→ S1

i = 1, . . . , N − 1 ,
i=N.

(21)

In this linear network, the number of reactions M is equal to the number of species
N. The maximum degree of coupling of this reaction network is 2, irrespective of the
size of the system (length of the chain), rendering it weakly coupled [35]. We hence
use PSSA-CR to evaluate the forward model with a computational complexity in
O(1) [37]. In the present test case, we limit ourselves to 3 species and 3 reactions,
i.e., N = M = 3. The parameter vector for this case is given by θ = [k1 , k2 , k3 ], since
the kinetics of linear reactions is independent of the volume Ω [39].
We simulate steady-state “experimental” data n̂ using PSSA-CR with ground
truth k1 = 2, k2 = 1.5, k3 = 3.2 (see Fig. 1a). We set the initial population of the
species to n1 (t = 0) = 50, n2 (t = 0) = 50, and n3 (t = 0) = 50 and sample a single
CME trajectory at equi-spaced time points with ∆texp = 0.1 between t = t0 and
t = t0 + (K − 1)∆texp with t0 = 2000 and K = 1001 for each of the 3 species S1 , S2 ,
and S3 . For the generated data we find zx = 7.
We generate trajectories from the forward model for every parameter vector θ
proposed by GaA using PSSA-CR between t = 0 and t = (K − 1)∆texp = 100, starting from the initial population ni (t = 0) = n̂i (t = t0 ).
Before turning to the actual parameter identification, we illustrate the topography
of the objective function landscape for the present example. We fix k3 = 3.2 to its
optimal value and perform a two-dimensional grid sampling for k1 and k2 over the
full search domain. We use 40 logarithmically spaced sample points per parameter,
resulting in 402 parameter combinations. For each combination we evaluate the objective function. The resulting landscapes of f1 (θ ), f2 (θ ), and f (θ ) are depicted in
Fig. 2a. Figure 2b shows refined versions around the global optimum. We see that
the moment-matching term f1 (θ ) is largely responsible for the global single-funnel
topology of the landscape. The autocorrelation term f2 (θ ) sharpens the objective
function near the global optimum and renders it locally more isotropic.

14

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

We perform both global optimization and ABC runs using GaA. In each of the 15
independent optimization runs, the number of function evaluations (FES) is limited
to MAX FES= 1000M = 3000. We set the initial step size to r(0) = 1 and perform all
searches in logarithmic scale of the parameters. Independent restarts from uniformly
random points are performed when the step size r drops below 10−4 [29]. For each
of the 15 independent runs, the 30 parameter vectors with the smallest objective
function value are collected and displayed in the box plot shown in the left panel of
Fig. 3a. All 450 collected parameter vectors have objective function values smaller
than 1.6. The resulting data suggest that the present method is able to accurately
determine the correct scale of the kinetic parameters from a single experimental
trajectory, although an overestimation of the rates is apparent.
We use the obtained optimization results for subsequent ABC runs. We conduct
15 independent ABC runs using cT = 2. The starting points for the ABC runs are
selected uniformly at random from the 450 collected parameter vectors in order to
ensure stable initialization. For each run, we again set MAX FES= 1000M = 3000.
The initial step size r(0) is set to 0.1, and the parameters are again explored in logarithmic scale. For all runs we observe rapid convergence of the empirical hitting

a
3

3

3

2

8

ï1

6
4

ï2

1
2.5

0

ï1

2

ï2

1.5

log10 k2

1

log10 k2

log10 k2

0

16

3

12
10

18

2

14
1

20

3.5

16
2

14
12

0
10

ï1

8
6

ï2
4

2
ï3
ï3

ï2

ï1

0

log10 k1

1

2

ï3
ï3

3

ï2

ï1

0

log10 k1

1

2

ï3
ï3

3

ï2

ï1

0

log10 k1

1

2

3

2

b
1

1

1

4.5

1.9

0.8

2.5

0.4

2
1.5

0.2

0
0

0.2

0.4

0.6

log10 k1

0.8

1

1.7

0.6

log10 k2

log10 k2

3

5.5
5

1.8

3.5
0.6

6
0.8

1.6
1.5

0.4

4.5

0.6

log10 k2

4

0.8

2

4
3.5

0.4

1.4
1.3

0.2

1

1.2

0.5

1.1

0
0

0.2

0.4

0.6

log10 k1

0.8

1

3
0.2

2.5
2

0
0

0.2

0.4

0.6

0.8

1

log10 k1

Fig. 2 a. Global objective function landscape for the cyclic chain over the complete search domain
for optimal k3 = 3.2. The three panels from left to right show f1 (θ ), f2 (θ ), and f (θ ), respectively.
b. A refined view of the global objective function landscape near the global optimum. The three
panels from left to right show f1 (θ ), f2 (θ ), and f (θ ), respectively. The white dots mark the true,
optimal parameters.

Global stochastic parameter identification

15

probability pemp to the optimal p = 1e (see Sec. 3). We collect the ABC samples
along with the means and covariances of GaA as soon as |pemp − p| < 0.05. As an
example we show the histograms of the posterior samples for a randomly selected
run in Fig. 3b. The means of the posterior distributions are again larger than the true
kinetic parameters. Using GaA’s means, covariance matrices, and the corresponding hitting probabilities that generated the posterior samples, we can construct an
ellipsoidal volume estimation [30]. This is done by multiplying each eigenvalue of
the average of the collected covariance matrices with c pemp = inv χn2 (pemp ), the ndimensional inverse Chi-square distribution evaluated at the empirical hitting probability. The product of these scaled eigenvalues and the volume of the n-dimensional
n

unit sphere, |S(n)| = Γ (πn 2+1) , then yields the ellipsoid volume with respect to a uni2
form distribution (see [30] for details). The resulting ellipsoid contains the optimal
kinetic parameter vector and is depicted in the right panel of Fig. 3a. It has a volume
of 0.045 in log-parameter space. This constitutes only 0.0208% of the initial search
space volume, indicating that GaA significantly narrows down the viable parameter
space around the true optimal parameters.

a

0.7
0.6

log10 k3

log10 ki

0.5
0.4
0.3
0.2
0.1

b

log10 k2
1

2
i

3

log10 k1

90

90

90

80

80

80

70

70

70

60

60

60

50

50

50

40

40

40

30

30

30

20

20

20

10

10

0
0.1

0.2

0.3

0.4

log10 k1

0.5

0.6

0
0

10
0.1

0.2

0.3

log10 k2

0.4

0.5

0
0.3

0.4

0.5

0.6

0.7

0.8

log10 k3

Fig. 3 a. Left panel: Box plot of the 30 best parameter vectors from each of the 15 independent
optimization runs. The blue dots mark the true parameter values. Right panel: Ellipsoidal volume
estimation of the parameter space below an objective-function threshold cT = 2 from a single ABC
run. b. Empirical posterior distributions of the kinetic parameters from the same single ABC run
with cT = 2. The red lines indicate the true parameters.

16

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

6.2 Strongly coupled reaction network: colloidal aggregation
The colloidal aggregation network is given by:
kon

1
→
0/ −−

S1

ki j

Si + S j −→ Si+ j
k̄i j

Si+ j −→ Si + S j
kioff

Si −−→

0/

i + j = 1, . . . , N
i + j = 1, . . . , N
i = 1, . . . , N .

(22)

j 2k
For this network of N species, the number of reactions is M = N2 + N + 1. The
maximum degree of coupling of this reaction network is proportional to N, rendering the network strongly coupled [35]. We hence use SPDM to evaluate the
forward model with a computational complexity in O(N) [35]. We use SPDM instead of PDM since the search path of GaA is unpredictable and could well generate parameters that lead to multi-scale networks. For this test case, we limit ourselves to two species, i.e., N = 2 and M = 5. The parameter vector for this case is
θ = [k11 , k̄11 , k1on , k1off , k2off , Ω ].
We perform GaA global optimization runs following the same protocol as for the
cyclic chain network with MAX FES = 1000(M + 1) = 6000.

6.2.1 At steady state
We simulate “experimental” data n̂ using SPDM with ground truth k11 = 0.1, k̄11 =
1.0, k1on = 2.1, k1off = 0.01, k2off = 0.1, and Ω = 15 (see Fig. 1b). We set the initial
population of the species to n1 (t = 0) = 0, n2 (t = 0) = 0, and n3 (t = 0) = 0 and
sample K = 1001 equi-spaced data points between t = t0 and t = t0 + (K − 1)∆texp
with t0 = 5000 and ∆texp = 0.1.
We generate trajectories from the forward model for every parameter vector θ
proposed by GaA using SPDM between t = 0 and t = (K − 1)∆texp = 100, stating
from the initial population ni (t = 0) = n̂i (t = t0 ).
The optimization results are summarized in the left panel of Fig. 4a. For each
of 15 independent runs, the 30 lowest-objective parameters are collected and shown
in the box plot. We observe that the true parameters corresponding to θ2 = k̄11 ,
θ3 = k1on , θ4 = k1off , and θ5 = k2off are between the 25th and 75th percentiles of the
identified parameters. Both the first parameter and the reaction volume are, on average, overestimated. Upon rescaling the kinetic rate constants with the estimated
volume, we find θ norm = [θ1 /θ6 , θ2 , θ3 θ6 , θ4 , θ5 ], which are the specific probability
rates of the reactions. The identified values are shown in the right panel of Fig. 4a.
The median of the identified θ3norm coincides with the true specific probability rate.
Likewise, θ1norm is closer to the 25th percentile of the parameter distribution. This

Global stochastic parameter identification

17

suggests a better estimation performance of GaA in the space of specific probability
rates, at the expense of not obtaining an estimate of the reaction volume.

6.2.2 In the transient phase
We simulate “experimental” data in the transient phase of the network dynamics
using the same parameters as above between t = t0 and t = (K − 1)∆texp with t0 =
0, ∆texp = 0.1, and K = 1001 (see Fig. 1c). We evaluate the forward model with
ni (t = 0) = n̂i (t = t0 ) to obtain trajectories between t = 0 and t = (K − 1)∆texp for
every proposed parameter vector θ .
The optimization results for the transient case are summarized in Fig. 4b. We
observe that the true parameters corresponding to θ3 = k1on , θ5 = k2off , and θ6 = Ω
are between the 25th and 75th percentiles of the identified parameters. The remaining
parameters are, on average, overestimated. In the space of rescaled parameters θ norm
we do not observe a significant improvement of the estimation.

3

3

2

2

log10 θinorm

a
log10 θi

1
0

1
0

ï1

ï1

ï2

ï2

ï3

ï3

1

2

3

4

5

6

1

2

3
i

4

5

1

2

3
i

4

5

i
3

3

2

2

log10 θinorm

b
log10 θi

1
0

1
0

ï1

ï1

ï2

ï2

ï3

ï3

1

2

3

4
i

5

6

Fig. 4 a. Left panel: Box plot of the 30 best parameter vectors from each of the 15 independent
optimization runs for the steady-state data set. Right panel: Box plots of the normalized parameters
(see main text for details). b. Left panel: Box plot of the 30 best parameter vectors from each of the
15 independent optimization runs for the transient data set. Right panel: Box plot of the normalized
parameters (see main text for details). The blue dots indicate the true parameter values.

18

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

7 Conclusions and Discussion
We have considered parameter estimation of stochastic biochemical networks from
single experimental trajectories. Parameter identification from single time series is
desirable in image-based systems biology, where per-cell estimates of the fluorescence evolution and its fluctuations are available. This enables quantifying cell–cell
variability on the level of network parameters. The histogram of the parameters
identified for different cells provides a biologically meaningful way of assessing
phenotypic variability beyond simple differences in the fluorescence levels.
We have proposed a novel combination of a flexible Monte Carlo method, the
Gaussian Adaptation (GaA) algorithm, and efficient exact stochastic simulation algorithms, the partial-propensity methods. The presented method can be used for
global parameter optimization, approximate Bayesian inference under a uniform
prior, and volume estimation of the viable parameter space. We have introduced an
objective function that measures closeness between a single experimental trajectory
and a single trajectory generated by the forward model. The objective function comprises a moment-matching and a time-autocorrelation part. This allows including
experimental readouts from, e.g., fluorescence photometry and fluorescence correlation spectroscopy.
We have applied the method to estimate the parameters of two monostable reaction networks from a single simulated temporal trajectory each, both at steady
state and during transient phases. We considered the linear cyclic chain network and
a non-linear colloidal aggregation network. For the linear model we were able to
robustly identify a small region of parameter space containing the true kinetic parameters. In the non-linear aggregation model, we could identify several parameter
vectors that fit the simulated experimental data well. There are two possible reasons
for this reduced parameter identifiability: either GaA cannot find the globally optimal region of parameter space due to high ruggedness and noise in the objective
function, or the non-linearity of the aggregation network modulates the kinetics in
a non-trivial way [39, 10]. Both cases are not accounted for in the current objective
function, thus leading to reduced performance for non-linear reaction networks.
We also used GaA as an adaptive MCMC method for approximate Bayesian inference of the posterior parameter distributions in the linear chain network. This enabled estimating the volume of the viable parameter space below a given objectivefunction value threshold. We found these volume estimates to be stable across independent runs. We thus believe that GaA might be a useful tool for exploring the
parameter spaces of stochastic systems.
Future work will include (i) alternative objective functions that include temporal cross-correlations between species and the derivative of the autocorrelation; (ii)
longer experimental trajectories; (iii) multi-stable and oscillatory systems; and (iv)
alternative global optimization schemes. Moreover, the applicability of the present
method to large-scale, non-linear biochemical networks and real-world experimental data will be tested in future work.

Global stochastic parameter identification

19

Acknowledgements RR was financed by a grant from the Swiss SystemsX.ch initiative (grant
WingX), evaluated by the Swiss National Science Foundation. This project was also supported
with a grant from the Swiss SystemsX.ch initiative, grant LipidX-2008/011, to IFS.

References
1. Albert, R.: Scale-free networks in cell biology. J. Cell Sci. 118(21), 4947–4957 (2005)
2. Andrieu, C., Thoms, J.: A tutorial on adaptive MCMC. Statistics and Computing 18(4), 343–
373 (2008)
3. Auger, A., Chatelain, P., Koumoutsakos, P.: R-leaping: Accelerating the stochastic simulation
algorithm by reaction leaps. J. Chem. Phys. 125, 084,103 (2006)
4. Barabási, A.L., Oltvai, Z.N.: Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5(2), 101–113 (2004)
5. Boys, R.J., Wilkinson, D.J., Kirkwood, T.B.L.: Bayesian inference for a discretely observed
stochastic kinetic model. Statistics and Computing 18(2), 125–135 (2008)
6. Cardinale, J., Rauch, A., Barral, Y., Székely, G., Sbalzarini, I.F.: Bayesian image analysis with
on-line confidence estimates and its application to microtubule tracking. In: Proc. IEEE Intl.
Symposium Biomedical Imaging (ISBI), pp. 1091–1094. IEEE, Boston, USA (2009)
7. Cinquemani, E., Milias-Argeitis, A., Summers, S., Lygeros, J.: Stochastic dynamics of genetic
networks: modelling and parameter identification. Bioinformatics 24(23), 2748–2754 (2008)
8. Gillespie, D.T.: A rigorous derivation of the chemical master equation. Physica A 188, 404–
425 (1992)
9. Gillespie, D.T.: Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys. 115(4), 1716–1733 (2001)
10. Grima, R.: Noise-induced breakdown of the Michaelis-Menten equation in steady-state conditions. Phys. Rev. Lett. 102(21), 218103 (2009). DOI 10.1103/PhysRevLett.102.218103. URL
http://link.aps.org/abstract/PRL/v102/e218103
11. Hafner, M., Koeppl, H., Hasler, M., Wagner, A.: ‘Glocal’ robustness analysis and model discrimination for circadian oscillators. PLoS Comput. Biol. 5(10), e1000,534 (2009)
12. Helmuth, J.A., Burckhardt, C.J., Greber, U.F., Sbalzarini, I.F.: Shape reconstruction of subcellular structures from live cell fluorescence microscopy images. J. Struct. Biol. 167, 1–10
(2009)
13. Helmuth, J.A., Sbalzarini, I.F.: Deconvolving active contours for fluorescence microscopy images. In: Proc. Intl. Symp. Visual Computing (ISVC), Lecture Notes in Computer Science,
vol. 5875, pp. 544–553. Springer, Las Vegas, USA (2009)
14. Jaynes, E.T.: Information theory and statistical mechanics. Phys. Rev. 106(4), 620–630 (1957).
DOI 10.1103/PhysRev.106.620
15. Kitano, H.: Computational systems biology. Nature 420(6912), 206–210 (2002)
16. Kitano, H.: Systems biology: A brief overview. Science 295(5560), 1662–1664 (2002)
17. Kjellström, G.: Network optimization by random variation of component values. Ericsson
Technics 25(3), 133–151 (1969)
18. Kjellström, G.: On the efficiency of Gaussian Adaptation. J. Optim. Theory Appl. 71(3),
589–597 (1991)
19. Kjellström, G., Taxen, L.: Stochastic optimization in system design. IEEE Trans. Circ. and
Syst. 28(7), 702–715 (1981)
20. Koeppl, H., Setti, G., Pelet, S., Mangia, M., Petrov, T., Peter, M.: Probability metrics to calibrate stochastic chemical kinetics. In: Proc. IEEE Intl. Symp. Circuits and Systems, pp.
541–544. Paris, France (2010)
21. Koutroumpas, K., Cinquemani, E., Kouretas, P., Lygeros, J.: Parameter identification for
stochastic hybrid systems using randomized optimization: A case study on subtilin production by Bacillus subtilis. Nonlinear Anal.: Hybrid Syst. 2(3), 786–802 (2008)

20

C. L. Müller, R. Ramaswamy, and I. F. Sbalzarini

22. Kurtz, T.G.: Relationship between stochastic and deterministic models for chemical reactions.
J. Chem. Phys. 57(7), 2976–2978 (1972)
23. Lakowicz, J.R.: Principles of Fluorescence Spectroscopy. Springer US (2006). DOI
10.1007/978-0-387-46312-4
24. Ljung, L.: Prediction error estimation methods. Circuits, systems, and signal processing 21(1),
11–21 (2002)
25. Marjoram, P., Molitor, J., Plagnol, V., Tavare, S.: Markov chain Monte Carlo without likelihoods. Proc. Natl. Acad. Sci. USA 100(26), 15,324–15,328 (2003)
26. Mason, O., Verwoerd, M.: Graph theory and networks in biology. Systems Biology, IET 1(2),
89–119 (2007)
27. Müller, C.L.: Exploring the common concepts of adaptive MCMC and Covariance Matrix Adaptation schemes. In: A. Auger, J.L. Shapiro, D. Whitley, C. Witt (eds.) Theory of Evolutionary Algorithms, no. 10361 in Dagstuhl Seminar Proceedings. Schloss
Dagstuhl – Leibniz-Zentrum für Informatik, Germany, Dagstuhl, Germany (2010). URL
http://drops.dagstuhl.de/opus/volltexte/2010/2813
28. Müller, C.L., Sbalzarini, I.F.: Gaussian Adaptation as a unifying framework for continuous
black-box optimization and adaptive Monte Carlo sampling. In: Proc. IEEE Congress on
Evolutionary Computation (CEC), pp. 2594–2601. Barcelona, Spain (2010)
29. Müller, C.L., Sbalzarini, I.F.: Gaussian Adaptation revisited — an entropic view on covariance matrix adaptation. In: Proc. EvoStar, Lect. Notes Comput. Sci., vol. 6024, pp. 432–441.
Springer, Istanbul, Turkey (2010)
30. Müller, C.L., Sbalzarini, I.F.: Gaussian Adaptation for robust design centering. In: Proc. EuroGen, Intl. Conf. Evolutionary and Deterministic Methods for Design, Optimization and Control. Capua, Italy (2011)
31. Munsky, B., Trinh, B., Khammash, M.: Listening to the noise: random fluctuations reveal gene
network parameters. Mol. Sys. Biol. 5(1), 318 (2009)
32. Poovathingal, S.K., Gunawan, R.: Global parameter estimation methods for stochastic biochemical systems. BMC Bioinformatics 11(1), 414 (2010)
33. Pritchard, J.K., Seielstad, M.T., Perez-Lezaun, A., Feldman, M.W.: Population growth of human Y chromosomes: A study of Y chromosome microsatellites. Molecular Biology and
Evolution 16(12), 1791–1798 (1999)
34. Qian, H., Elson, E.L.: Fluorescence correlation spectroscopy with high-order and dual-color
correlation to probe nonequilibrium steady states. Proc. Natl. Acad. Sci. USA 101(9), 2828–
2833 (2004)
35. Ramaswamy, R., González-Segredo, N., Sbalzarini, I.F.: A new class of highly efficient exact
stochastic simulation algorithms for chemical reaction networks. J. Chem. Phys. 130(24),
244,104 (2009)
36. Ramaswamy, R., Sbalzarini, I.F.: Fast exact stochastic simulation algorithms using partial
propensities. In: Proc. ICNAAM, Numerical Analysis and Applied Mathematics, International Conference, pp. 1338–1341. AIP (2010)
37. Ramaswamy, R., Sbalzarini, I.F.: A partial-propensity variant of the composition-rejection
stochastic simulation algorithm for chemical reaction networks. J. Chem. Phys. 132(4),
044,102 (2010)
38. Ramaswamy, R., Sbalzarini, I.F.: A partial-propensity formulation of the stochastic simulation
algorithm for chemical reaction networks with delays. J. Chem. Phys. 134, 014,106 (2011)
39. Ramaswamy, R., Sbalzarini, I.F., González-Segredo, N.: Noise-induced modulation of the relaxation kinetics around a non-equilibrium steady state of non-linear chemical reaction networks. PLoS ONE 6(1), e16,045 (2011)
40. Reinker, S., Altman, R.M., Timmer, J.: Parameter estimation in stochastic biochemical reactions. IEE Proc.-Syst. Biol 153(4), 168 (2006)
41. Stock, G., Ghosh, K., Dill, K.A.: Maximum Caliber: A variational approach applied to twostate dynamics. J. Chem. Phys. 128(19), 194,102 (2008)
42. Strogatz, S.H.: Exploring complex networks. Nature 410, 268–276 (2001)

Global stochastic parameter identification

21

43. Toni, T., Welch, D., Strelkowa, N., Ipsen, A., Stumpf, M.P.H.: Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J. R. Soc.
Interface 6(31), 187–202 (2009)
44. Wolkenhauer, O.: Systems biology: The reincarnation of systems theory applied in biology?
Briefings in Bioinformatics 2(3), 258 (2001)
45. Zechner, C., Ganguly, A., Pelet, S., Peter, M., Koeppl, H.: Accounting for extrinsic variability in the estimation of stochastic rate constants. Intl. J. Robust Nonlinear Control (2011,
submitted)
46. Zechner, C., Pelet, S., Peter, M., Koeppl, H.: Recursive Bayesian estimation of stochastic rate
constants from heterogeneous cell populations. In: Proc. IEEE CDC, Conference on Decision
and Control (2011, submitted)

