We present a novel method for estimating rate coefficients from noisy observations of concentration levels at discrete time points. This is traditionally done by computing the least-squares estimator. However, estimation of the error function generally requires solving the reaction rate equations, which can become computationally unfeasible. Here we present an alternative approach based on a probabilistic, generative model of the variations in reactant concentration. Our method returns the rate coefficients, the level of noise and an error range on the estimates of rate constants. Its probabilistic formulation is key to a principled handling of the noise inherent in biological data, and it allows a number of further extensions

Lecca, Paola

Palmisano, Alida

Priami, Corrado

Sanguinetti, Guido

English

University of Trento : UNITN-eprints

Technical Report CoSBi 20/2007
A new model for kinetic parameter estimation
in biochemical reactions
Paola Lecca
The Microsoft Research - University of Trento
Centre for Computational and Systems Biology
lecca@cosbi.eu
Guido Sanguinetti
Dept. of Computer Science, University of Sheffield, UK
guido@dcs.shef.ac.uk
Alida Palmisano
The Microsoft Research - University of Trento
Centre for Computational and Systems Biology
DISI - University of Trento
palmisano@cosbi.eu
Corrado Priami
The Microsoft Research - University of Trento
Centre for Computational and Systems Biology
DISI - University of Trento
priami@cosbi.eu
A new model for kinetic parameter estimation in
biochemical reactions.
Paola Lecca, Guido Sanguinetti, Alida Palmisano, and Corrado Priami
Abstract
We present a novel method for estimating rate coefficients from
noisy observations of concentration levels at discrete time points. This
is traditionally done by computing the least-squares estimator. How-
ever, estimation of the error function generally requires solving the re-
action rate equations, which can become computationally unfeasible.
Here we present an alternative approach based on a probabilistic, gen-
erative model of the variations in reactant concentration. Our method
returns the rate coefficients, the level of noise and an error range on
the estimates of rate constants. Its probabilistic formulation is key to
a principled handling of the noise inherent in biological data, and it
allows a number of further extensions.
1
1 Introduction
The ability to infer kinetic parameters of biochemical reactions is emerging
as a crucial problem in system biology. Estimates of reaction rates are often
crucial for modelling and simulating the biochemical dynamics, yet their
direct measurement is an outstanding experimental problem. The current
need for new techniques for integrative and predictive modelling [1, 7, 6] is
based on the realisation that traditional reductionist approaches have to be
complemented by the reconstruction of a model of how systems function as a
whole. The importance and the great interest that the scientific community
attributes to this problem are evidenced by the increasing number of recent
works concerning this research topic.
In this paper we present a novel method for estimating rate coefficients
from noisy observations of concentration levels at discrete time points. Given
a number of N reactant species, we observe time series concentrations for
each of the species, gathered in N state vectors ~X1, . . . , ~XN . The relation
between the instantaneous rate of reaction and the concentrations of the
reactants at any moment is given by the law of mass action: i.e. the rate at
which a substance takes part in a reaction is proportional to its concentra-
tion raised to a power which represents the number of molecules taking part
in the reaction. Such formulation is made for simultaneous as well as iso-
lated reactions, and for heterogeneous as well as homogeneous systems. The
goal is the estimation of the constant of proportionality. This is tradition-
ally done by computing the least-squares estimators (LSE), resulting from
the minimisation of the sum of squared differences between the observed
data and the model. However, estimation of the error function generally
requires solving the reaction rate equations, which can easily become com-
putationally unfeasible. Here we present an alternative approach based on a
probabilistic, generative model of the variations in reactant concentration.
Our method discretizes the law of mass action and provides a tool to predict
the values of the variables ~Xi at time t, conditioned on their values at the
previous time point. The variations of the species concentration at different
time points are conditionally independent by the Markovian nature of the
discrete model of the law of mass action. Assuming the observation noise
to be Gaussian with variance σ2, the probability of observing a variation
Di for the concentration [X]i of species i between time tk−1 and tk is a
Gaussian with variance depending on σ and mean the expectation value of
the law mass action function under the noise distribution. The discretiza-
tion of the law of mass action provides a model for the variations of the
species concentration, rather than a model for the time-trajectory of the
species concentrations. This makes the evaluation of the expectation value
of law mass action function simpler and analytically tractable. The rate
coefficients and the level of noise σ2 are then obtained by maximising the
likelihood function defined by the observed variations. Our method returns
2
the rate coefficients, the level of noise and an error range on the estimates of
rate constants. Its probabilistic formulation is key to a principled handling
of the noise inherent in biological data, and it allows a number of further
extensions, such as a fully Bayesian treatment of the parameter inference
and automated model selection strategies based on the comparison between
marginal likelihoods of different models. Finally, the implementation of this
method may be used as an interface tool, connecting the outcomes of the
wet-lab activity for the concentration measurements and the softwares for
the simulation of chemical kinetics. We tested our algorithm on some real
case-study, whose rate coefficients were experimentally determined and well
documented in literature for simple chemical reactions and for more complex
biochemical interactions.
2 The model for inference
ConsiderN reactant species, S1, S2, . . . , SN , with concentrationsX1, X2, . . . , XN ,
that evolve according to a system of rate equations
dXi
dt
= fi(X(i)(t); θi) (1)
where θi , i = 1, 2, . . . , N , is the vector of the rate coefficients, which are
present in the expression of the function fi . We wish to estimate the set
of parameters Θ = ∪θi (i = 1, 2, . . . , N), whose element θi is the set of rate
coefficients appearing in the rate equations of i-th species, therefore
θ1 = {θ11, θ12, . . . , θ1N1}, . . . , θN = {θN1, θN2, . . . , θNNN }
X(i) is the vector of concentrations of chemicals that are present in the
expression of the function fi for the species i. According to the law of mass
action, the functions fi have the general form
fi(X(i)(t); θi) =
= θi1
∏
w∈S1⊆[1,N ]
Xαww + . . .+ θiNi
∏
w∈SNi⊆[1,N ]
Xαww =
Ni∑
h=1
(
θih
∏
w∈Sh
Xαww
)
(2)
where αw ∈ R, and Ni is the number of parameter in the fi rate equation
The rate equations in (2) form the so-called Generalized Mass Action law.
We assume we have noisy observations Xˆi = Xi +  at times t0, . . . , tM ,
where  ∼ N (0, σ2) is a Gaussian noise term with mean zero and variance σ.
With this choice we are assuming that the concentration measurements are
not significantly affected by systematic errors, but by uncontrolled random
3
errors and that an error is equally likely to occur in either positive or negative
direction with respect to the symmetry axis of the distribution.
We also assume a number M of concentration measurements for each
considered species. Approximating the rate equation (1) as a finite difference
equation between the observation times, gives
Xi(tk) = Xi(tk−1) + (tk − tk−1)fi(X(i)(tk−1); θi) (3)
where k = 1, . . . ,M . In Eq. (3) the rate equation is viewed as a model
of increments/decrements of reactant concentrations; i.e., given a value of
the variables at time tk−1 , the model can be used to predict the value
at the next time point tk. Increments/decrements between different time
points are conditionally independent by the Markov nature of the model
(3). Therefore, given the Gaussian model for the noise, it is possible to
estimate the probability to observe the value Xˆi(tk) given the model at time
tk−1, Xi(tk−1), and the set of parameters θi, as
p
(
Xˆi(tk−1)|Xi(tk−1)
)
= N
(
Xi(tk−1) + (tk − tk−1)fi(Xi(tk−1, θi)), σ2
)
(4)
We then also have that the true value of Xi(tk) is normally distributed
around the observed value Xˆi(tk), so that
p
(
Xi(tk−1)|Xˆi(tk−1))
)
= N
(
Xˆi(tk−1), σ2
)
= (5)
=
1√
2piσ
exp
[
− (Xi(tk−1)− Xˆi(tk−1))
2
2σ2
]
Therefore, the probability to observe a variation Di(tk) = Xi(tk)−Xi(tk−1)
for the concentration of the i-th species between the time tk−1 and tk, given
the parameter vector θi is
p(Di(tk)|θi, σ) = N
(
E[fi(X(i)(tk−1), θi)], 2σ2
)
(6)
and
E[fi(X(i)(tk−1, θi))] =
∫
Ω
X(i)
fi(X(i)(tk−1), θi)
Ki∏
i=1
[
pi
(
Xi(tk−1)|Xˆi(tk−1)
)]
dX(i)
(7)
where ΩX(i) is the sample space of X
(i), and Ki is the number of chemical
species in the expression for fi. While the increments/decrements are condi-
tionally independent given the starting point Xi (tk), the random variables
Di(tk) are not independent of each other. Intuitively, if Xi(tk) happens to
be below its expected value because of random fluctuations, then the follow-
ing increment Di(tk+1) can be expected to be bigger as a result, while the
4
previous one Di(tk) will be smaller. A simple calculation allows us to obtain
the covariance matrix of the vector of increments for the i-th species. This
is a banded matrix Ci ≡ C = Cov(Di) with diagonal elements given by
E
[
D2i (tk)− E[D2i (tk)]
]
= 2σ2
and a non-zero band above and below the diagonal given by
E
[
(Di(tk)− E[Di(tk)])(Di(tk−1)− E[Di(tk−1)])
]
= −σ2
with all other entries zero. The likelihood for the observed increments/decrements
therefore will be
p(D|Θ) =
N∏
i=1
N (Di|mi(Θ),C) =
(
1√
2pi det(C)
)N
e
∑N
i=1− 12 (Di−mi)TC−1(Di−mi)
(8)
where D = {D1, . . . ,DN}, Di = Di(t1), Di(t2), . . . Di(tM ) (i = 1, 2, . . . , N),
and mi(tk−1) ≡ E
[
fi(X(tk−1), θi)
]
.
The Eq. (8) can be optimized w. r. t. the parameters Θ = (θ1, θ2, . . . , θN )
of the model to yield estimates of the parameters themselves and of the noise
level. The chief numerical problem of this approach is the computation of
the expectations of the rate functions given by equation (7). Non-integer
values of the coefficients α can make estimating the integral analytically
difficult. We propose an approximate method in which the Gaussian noise
is replaced by an approximate uniform (white) noise, with the amplitude of
the uniform noise being obtained as a sample from the Gaussian cumulative
distribution function. At the first order, for small σ, we can approximate
the Gaussian with zero mean and variance σ with an uniform distribution
defined on the interval [−
√
2piσ
4 ,
√
2piσ
4 ], so that
Ki∏
i=1
pi =
Ki∏
i=1
χi (9)
where
χi(Xi) =
{
2√
2piσ
if −
√
2piσ
4 ≤ Xi ≤
√
2piσ
4
0 otherwise.
This approximation makes the calculation of the expectation value of the
rate equation (Eq. (7)) simpler and reduces the computational time of the
procedure. Moreover, experiments not illustrated in this paper demonstrate
that it does not influence the accuracy of the parameter estimates until σ is
less that 30% of the concetration measurement.
Substituting Eq. (9) in Eq, (7) gives
5
E[fi(X(i)(tk−1), θ)] =
(
2√
2piσ
)Ki ∫ Xˆ+√2piσ
4
Xˆ−
√
2piσ
4
fi(X(i)(tk−1), θi)dX(i) (10)
Now, substituting Eq. (2) in Eq. (10) leads to
E[fi(X(i)(tk−1), θi)] =
=
(
2√
2piσ
)Ki{ Ni∑
h=1
θih
[(√
2piσ
2
)#(S−Sh)
×
×
∏
w∈Sh
1
αw + 1
((
Xˆw +
√
2piσ
4
)αw+1 − (Xˆw − √2piσ4 )αw+1
)]}
(11)
where S is the set containing the indexes referring to all the Ki species
appearing in fi, and αw 6= −1. In case some orders are equal to -1 Eq. (11)
takes the following form
E[fi(X(i)(tk−1), θi)] =
(
2√
2piσ
)Ki Ni∑
h=1
θih
{(√
2piσ
2
)#(S−Sh)
×
×
[ ∏
w∈S′h
1
αw + 1
((
Xˆw +
√
2piσ
4
)αw+1 − (Xˆw − √2piσ4 )αw+1
)]
×
×
[ ∏
w∈S′′h
ln
Xˆw +
√
2piσ
4
Xˆw −
√
2piσ
4
]}
(12)
where S′h is the set of indexes {h′1, h′2, . . . , h′s} such that αh′ 6= −1 ∀h′ ∈ S′h,
and S′′h is the set of indexes {h′′1, h′′2, . . . , h′′s} such that αh′′ = −1 ∀h′′ ∈ S′′h.
If in the Eq. (8), mi is substituted with the expression (11) or (12),
Eq. (8) becomes more tractable and can be optimized w. r. t. the param-
eters Θ = (θ1, θ2, . . . , θN ) and σ. The values of the model’s parameters for
which p(D|Θ) has a maximum are the most likely values giving the observed
kinetics.
3 The estimation errors
For a zero-th order reaction
∅ θ−→ A
A
θ−→ ∅
6
∆θ =
{( ∂θ
∂A0
A0
)2
+
( ∂θ
∂Af
Af
)2} 1
2 (13)
∂k
∂Af
=
(Af −A0)
tf − t0 (14)
∂k
∂A0
=
(A0 −Af )
tf − t0 (15)
For a first order reaction
A
θ−→ P
if A0 (Af ) is the initial (final) concentration of the molecules of the
species A, the error that the algorithm makes in estimating θ is given by
the following
∆θ =
{( ∂θ
∂A0
A0
)2
+
( ∂θ
∂Af
Af
)2} 1
2 (16)
Since
θ =
1
tf − t0 (lnA0 − lnAf ) (17)
Eq. 16 becomes
∆θ =
{( ∂θ
A0
A0
)2
+
( ∂θ
Af
Af
)2} 1
2 (18)
For a second order reaction
A+B θ−→ P
if A0 (Af ) and B0 (Bf ) are the initial (final) concentrations of the molecules
of species A and B, respectively, the error that the algorithm makes in
estimating θ is given by the following
∆θ =
{( ∂θ
∂A0
A0
)2
+
( ∂θ
∂B0
B0
)2
+
( ∂θ
∂Af
Af
)2
+
( ∂θ
∂Bf
Bf
)2} 1
2 (19)
where
∂θ
∂A0
=
1
tf − t0
{ 1
(A0 −B0)2 ln
AfB0
A0Bf
− 1
A0 −B0
Bf
A0
}
(20)
∂θ
∂B0
=
1
tf − t0
{
− 1
(A0 −B0)2 ln
AfB0
A0Bf
+
1
A0 −B0
Af
B0
}
(21)
7
∂θ
∂Af
=
1
tf − t0
1
A0 −B0
B0
Af
(22)
∂θ
∂Bf
=
1
tf − t0
1
B0 −A0
A0
Bf
(23)
4 Case studies
All the mathematical operations described above has been automated through
the implementation of an algorithm in C++. We tested our algorithm on
some real case-study, whose rate coefficients were experimentally determined
and well documented in literature for simple chemical reactions and for more
complex biochemical interactions. We show the ability of our algorithm of
obtaining reasonable estimates for the rate coefficients in the following cases
(see Table 1):
R1 Bromine molecule formation 2Br −→ Br2. The data referred to this
reaction are the bromine concentration time-points after a flash pho-
tolysis of a mixture of bromine and SF6 with [Br2]/[SF6] = 3.2×10−2
(experiment of Graff and Lang [11]).
R2. In presence of an acid solution of phenol, the ion IO−3 is reduced to
IO−2 by Br
− accordingly to the reaction
IO−3 + 2Br
− + 2H+ −→ IO−2 +Br2 +H2O
With [IO−3 ] = 5×10−3M and [Br−] = 1×10−2M , Sharma and Gupta
obtained the following data for a solution at 35 C with [C6H5OH] =
2× 10−2M [11].
R3. the first-order decomposition in azomethane CH3N2CH3 −→ CH3CH3+
N2 ([13]);
R4. the alkaline hydrolysis of ethyl nitrobenzoathe ([13]);
R5. the gluthathione S-transferase-catalysed dehalogenation of haloaro-
matic compounds with aromatic substrates ([15, 10]);
R6. degradation of IkBa-UB in breast cancer cells (refer to [3] for a com-
parison);
R7. Nuclear factor (NF) kB translocation from cytoplasm into the nucleus
in breast cancer cells (r. t. [3] for a comparison);
R8. the glucose consumption rate in astrocytes [9].
8
Reaction k k σ
(experimental) (estimate) (estimate)
R1 (2.75± 1.00)× 108M−1s−1 (1.93± 1.28)× 108M−1s−1 3.2× 10−6
R2 (13.8± 1.00)M−2s−2 (15.92± 6.98)M−2s−1 3.5× 10−1
R3 (3.60± 2.16)× 10−4s−1 (0.80± 0.45)× 10−4s−1 5.36× 10−3
R4 (8.1± 2.31)× 10−2s−1 (5.9± 1.5)× 10−2s−1 3.3× 10−3
R5 (7.4± 0.5)× 10−4s−1 (2.60± 1.17)× 10−4s−1 4.1× 10−3
R6 (5.56± 0.11)× 10−4s−1 (3.1± 2.2)× 10−4s−1 5× 10−3
R7 (3.85± 2.67)× 10−4s−1 (1.3± 0.4)× 10−4s−1 8.5× 10−3
R8 (1.33± 0.04)× 10−2s−1 (1.24± 0.78)× 10−2s−1 6× 10−4
Table 1: The estimates of the rate coefficients compared with the experimental
values.
5 Conclusion
The ability of this algorithm to estimate the expected value of the kinetic
rate constant significantly depends on the experimental error on the data. If
the time-course measurements of the reactant concentration are affected by a
high level of noise, this technique estimate a rate constant very close to zero,
indicating that there is no dynamic in the systems. Therefore the observed
variations in the concentration values are attributed exclusively to the noise.
For an accurate estimate of the rate constants the procedure requires high-
quality and high-resolution data, that at the time of writing are relatively
difficult to obtain, In the meantime experimental biology evolves to make
such data more readily available, the development of inference techniques
allow identification of experimental design issues that need to be addressed
to get more accurate data for more reliable estimate of rate constants.
References
[1] J. M. Bower, H. Bolouri (Eds.), Computational Modeling of Genetic and
Biochemical Networks, Cambridge, Massachusetts, MIT Press, 2001.
[2] E. . Crampin, P. E. McSharry, and S. Schnell, Extracting Biochemical
Reaction Kinetics from Time Series Data, M. GH. Negoita et al. (Eds):
KES 2004, LNAI 3214, p. 329-336, Springer-Verlag Berlin Heidelberg
2004.
[3] J. W. Fuseler, D. M. Merril, J. A. Rogers, M. B. Grisham and R. E.
Wolf, Analysis and quatitation of NF-kB nuclear translocation in tu-
mor necrosis factor alpha (TNF-α) activated vascular endothelial cells,
Microsc. Microanal. 12, pp. 269-276, 2006.
9
[4] A. Golightly and D. J. Wilkinson, Bayesian sequential inference for
stochastic kinetic biochemical network inference, Journal of Computa-
tional Biology vol. 13, no. 3 pp. 838-851, 2006.
[5] A. Golightly and D. J. Wilkinson, Bayesian inference for stochastic
kinetic models using a diffusion approximation, Biometrics, vol. 61, no.
3, pp. 781-788, 2005.
[6] T. B. L. Kirkwood, R. J. Boys, C. S. Gillespie, C. J. Proctor, D. P.
Shanley, and D. J. Wilkinson, Toward an e-biology of ageing: inte-
grating theory and data, Nature Reviews Molecular Cell Biology 4, pp.
243-249, 2003.
[7] H. Kitano (Ed.), Foundations of System Biology, Cambridge, Mas-
sachusetts, MITOress, 2001.
[8] N. Lalam and C. A. J. Klassen, Pseudo maximum likelihood estimation
for differential equations,
[9] P. Lecca and M. Lecca, Molecular mechanism of glutamate-triggered
brain glucose metabolism: a parametric model from FDG PET-scans,
In Procs of BVAI 2007, To appear in LNCS, 2007.
[10] I. Mancini, G. Guella, G. Chiasera and F. Pietra, Glutathione S-
transferase catalysed dehalogenation of haloaromatic compounds which
lack nitro groups near the reaction centre, Tetrahedron Letters, 39, pp.
1611-1614, 1998.
[11] C. R. Metz, Chimica fisica, McGraw-Hill, 1994.
[12] I. Nachman, A. Regev and N. Friedman, Inferring quantitative models
of regulatory networks from expression data, Bioinformatics, vol. 20, no.
1, pp. 248-256, 2004.
[13] V. Nummert and M. Piirsalu, Kinetic study of hydrolysis of benzoates,
J. Chem. Soc., Perkin Trans. 2, pp. 583-594, 2000.
[14] S. Reinker, R. M. Altman and J. Timmer, Parameter estimation in
stochastic biochemical reactions, IEE Proc. Syst. Biol, vol. 153, no. 4,
July 2006.
[15] B. Rossi, I. Mancini and G. Guella, Kinetic studies of glutathione S-
transferase-catalysed process through on- line liquid chromatography-
electrospray mass spectrometry, Philosophical Magazine, 84(13-16), pp.
1373-1382, 2004.
[16] J. Yu, V. A. Smith, P. P. Wang, a. J. Hartemink and E. D. Jarvis,
Advances to Bayesian network inference for generating causal networks
10
from observational biological data, Bioinformatics vol 20, no. 18, pp.
3594 - 3603, 2004.
[17] R. J. Boys, D. J. Wilkinson and T. B. Kirkwood, Bayesian inference for
a discretely observed stochastic kinetic model, currently in submission
(available at http://www.mas.ncl.ac.uk/∼ndjw1/docs/BWK04.pdf).
11


A new model for kinetic parameter estimation in biochemical reactions.

Unitn-eprints Research

Technical Report CoSBi 20/2007A new model for kinetic parameter estimationin biochemical reactionsPaola LeccaThe Microsoft Research - University of TrentoCentre for Computational and Systems Biologylecca@cosbi.euGuido SanguinettiDept. of Computer Science, University of Sheffield, UKguido@dcs.shef.ac.ukAlida PalmisanoThe Microsoft Research - University of TrentoCentre for Computational and Systems BiologyDISI - University of Trentopalmisano@cosbi.euCorrado PriamiThe Microsoft Research - University of TrentoCentre for Computational and Systems BiologyDISI - University of Trentopriami@cosbi.euA new model for kinetic parameter estimation inbiochemical reactions.Paola Lecca, Guido Sanguinetti, Alida Palmisano, and Corrado PriamiAbstractWe present a novel method for estimating rate coefficients fromnoisy observations of concentration levels at discrete time points. Thisis traditionally done by computing the least-squares estimator. How-ever, estimation of the error function generally requires solving the re-action rate equations, which can become computationally unfeasible.Here we present an alternative approach based on a probabilistic, gen-erative model of the variations in reactant concentration. Our methodreturns the rate coefficients, the level of noise and an error range onthe estimates of rate constants. Its probabilistic formulation is key toa principled handling of the noise inherent in biological data, and itallows a number of further extensions.11 IntroductionThe ability to infer kinetic parameters of biochemical reactions is emergingas a crucial problem in system biology. Estimates of reaction rates are oftencrucial for modelling and simulating the biochemical dynamics, yet theirdirect measurement is an outstanding experimental problem. The currentneed for new techniques for integrative and predictive modelling [1, 7, 6] isbased on the realisation that traditional reductionist approaches have to becomplemented by the reconstruction of a model of how systems function as awhole. The importance and the great interest that the scientific communityattributes to this problem are evidenced by the increasing number of recentworks concerning this research topic.In this paper we present a novel method for estimating rate coefficientsfrom noisy observations of concentration levels at discrete time points. Givena number of N reactant species, we observe time series concentrations foreach of the species, gathered in N state vectors ~X1, . . . , ~XN . The relationbetween the instantaneous rate of reaction and the concentrations of thereactants at any moment is given by the law of mass action: i.e. the rate atwhich a substance takes part in a reaction is proportional to its concentra-tion raised to a power which represents the number of molecules taking partin the reaction. Such formulation is made for simultaneous as well as iso-lated reactions, and for heterogeneous as well as homogeneous systems. Thegoal is the estimation of the constant of proportionality. This is tradition-ally done by computing the least-squares estimators (LSE), resulting fromthe minimisation of the sum of squared differences between the observeddata and the model. However, estimation of the error function generallyrequires solving the reaction rate equations, which can easily become com-putationally unfeasible. Here we present an alternative approach based on aprobabilistic, generative model of the variations in reactant concentration.Our method discretizes the law of mass action and provides a tool to predictthe values of the variables ~Xi at time t, conditioned on their values at theprevious time point. The variations of the species concentration at differenttime points are conditionally independent by the Markovian nature of thediscrete model of the law of mass action. Assuming the observation noiseto be Gaussian with variance σ2, the probability of observing a variationDi for the concentration [X]i of species i between time tk−1 and tk is aGaussian with variance depending on σ and mean the expectation value ofthe law mass action function under the noise distribution. The discretiza-tion of the law of mass action provides a model for the variations of thespecies concentration, rather than a model for the time-trajectory of thespecies concentrations. This makes the evaluation of the expectation valueof law mass action function simpler and analytically tractable. The ratecoefficients and the level of noise σ2 are then obtained by maximising thelikelihood function defined by the observed variations. Our method returns2the rate coefficients, the level of noise and an error range on the estimates ofrate constants. Its probabilistic formulation is key to a principled handlingof the noise inherent in biological data, and it allows a number of furtherextensions, such as a fully Bayesian treatment of the parameter inferenceand automated model selection strategies based on the comparison betweenmarginal likelihoods of different models. Finally, the implementation of thismethod may be used as an interface tool, connecting the outcomes of thewet-lab activity for the concentration measurements and the softwares forthe simulation of chemical kinetics. We tested our algorithm on some realcase-study, whose rate coefficients were experimentally determined and welldocumented in literature for simple chemical reactions and for more complexbiochemical interactions.2 The model for inferenceConsiderN reactant species, S1, S2, . . . , SN , with concentrationsX1, X2, . . . , XN ,that evolve according to a system of rate equationsdXidt= fi(X(i)(t); θi) (1)where θi , i = 1, 2, . . . , N , is the vector of the rate coefficients, which arepresent in the expression of the function fi . We wish to estimate the setof parameters Θ = ∪θi (i = 1, 2, . . . , N), whose element θi is the set of ratecoefficients appearing in the rate equations of i-th species, thereforeθ1 = {θ11, θ12, . . . , θ1N1}, . . . , θN = {θN1, θN2, . . . , θNNN }X(i) is the vector of concentrations of chemicals that are present in theexpression of the function fi for the species i. According to the law of massaction, the functions fi have the general formfi(X(i)(t); θi) == θi1∏w∈S1⊆[1,N ]Xαww + . . .+ θiNi∏w∈SNi⊆[1,N ]Xαww =Ni∑h=1(θih∏w∈ShXαww)(2)where αw ∈ R, and Ni is the number of parameter in the fi rate equationThe rate equations in (2) form the so-called Generalized Mass Action law.We assume we have noisy observations Xˆi = Xi +  at times t0, . . . , tM ,where  ∼ N (0, σ2) is a Gaussian noise term with mean zero and variance σ.With this choice we are assuming that the concentration measurements arenot significantly affected by systematic errors, but by uncontrolled random3errors and that an error is equally likely to occur in either positive or negativedirection with respect to the symmetry axis of the distribution.We also assume a number M of concentration measurements for eachconsidered species. Approximating the rate equation (1) as a finite differenceequation between the observation times, givesXi(tk) = Xi(tk−1) + (tk − tk−1)fi(X(i)(tk−1); θi) (3)where k = 1, . . . ,M . In Eq. (3) the rate equation is viewed as a modelof increments/decrements of reactant concentrations; i.e., given a value ofthe variables at time tk−1 , the model can be used to predict the valueat the next time point tk. Increments/decrements between different timepoints are conditionally independent by the Markov nature of the model(3). Therefore, given the Gaussian model for the noise, it is possible toestimate the probability to observe the value Xˆi(tk) given the model at timetk−1, Xi(tk−1), and the set of parameters θi, asp(Xˆi(tk−1)|Xi(tk−1))= N(Xi(tk−1) + (tk − tk−1)fi(Xi(tk−1, θi)), σ2)(4)We then also have that the true value of Xi(tk) is normally distributedaround the observed value Xˆi(tk), so thatp(Xi(tk−1)|Xˆi(tk−1)))= N(Xˆi(tk−1), σ2)= (5)=1√2piσexp[− (Xi(tk−1)− Xˆi(tk−1))22σ2]Therefore, the probability to observe a variation Di(tk) = Xi(tk)−Xi(tk−1)for the concentration of the i-th species between the time tk−1 and tk, giventhe parameter vector θi isp(Di(tk)|θi, σ) = N(E[fi(X(i)(tk−1), θi)], 2σ2)(6)andE[fi(X(i)(tk−1, θi))] =∫ΩX(i)fi(X(i)(tk−1), θi)Ki∏i=1[pi(Xi(tk−1)|Xˆi(tk−1))]dX(i)(7)where ΩX(i) is the sample space of X(i), and Ki is the number of chemicalspecies in the expression for fi. While the increments/decrements are condi-tionally independent given the starting point Xi (tk), the random variablesDi(tk) are not independent of each other. Intuitively, if Xi(tk) happens tobe below its expected value because of random fluctuations, then the follow-ing increment Di(tk+1) can be expected to be bigger as a result, while the4previous one Di(tk) will be smaller. A simple calculation allows us to obtainthe covariance matrix of the vector of increments for the i-th species. Thisis a banded matrix Ci ≡ C = Cov(Di) with diagonal elements given byE[D2i (tk)− E[D2i (tk)]]= 2σ2and a non-zero band above and below the diagonal given byE[(Di(tk)− E[Di(tk)])(Di(tk−1)− E[Di(tk−1)])]= −σ2with all other entries zero. The likelihood for the observed increments/decrementstherefore will bep(D|Θ) =N∏i=1N (Di|mi(Θ),C) =(1√2pi det(C))Ne∑Ni=1− 12 (Di−mi)TC−1(Di−mi)(8)where D = {D1, . . . ,DN}, Di = Di(t1), Di(t2), . . . Di(tM ) (i = 1, 2, . . . , N),and mi(tk−1) ≡ E[fi(X(tk−1), θi)].The Eq. (8) can be optimized w. r. t. the parameters Θ = (θ1, θ2, . . . , θN )of the model to yield estimates of the parameters themselves and of the noiselevel. The chief numerical problem of this approach is the computation ofthe expectations of the rate functions given by equation (7). Non-integervalues of the coefficients α can make estimating the integral analyticallydifficult. We propose an approximate method in which the Gaussian noiseis replaced by an approximate uniform (white) noise, with the amplitude ofthe uniform noise being obtained as a sample from the Gaussian cumulativedistribution function. At the first order, for small σ, we can approximatethe Gaussian with zero mean and variance σ with an uniform distributiondefined on the interval [−√2piσ4 ,√2piσ4 ], so thatKi∏i=1pi =Ki∏i=1χi (9)whereχi(Xi) ={2√2piσif −√2piσ4 ≤ Xi ≤√2piσ40 otherwise.This approximation makes the calculation of the expectation value of therate equation (Eq. (7)) simpler and reduces the computational time of theprocedure. Moreover, experiments not illustrated in this paper demonstratethat it does not influence the accuracy of the parameter estimates until σ isless that 30% of the concetration measurement.Substituting Eq. (9) in Eq, (7) gives5E[fi(X(i)(tk−1), θ)] =(2√2piσ)Ki ∫ Xˆ+√2piσ4Xˆ−√2piσ4fi(X(i)(tk−1), θi)dX(i) (10)Now, substituting Eq. (2) in Eq. (10) leads toE[fi(X(i)(tk−1), θi)] ==(2√2piσ)Ki{ Ni∑h=1θih[(√2piσ2)#(S−Sh)××∏w∈Sh1αw + 1((Xˆw +√2piσ4)αw+1 − (Xˆw − √2piσ4 )αw+1)]}(11)where S is the set containing the indexes referring to all the Ki speciesappearing in fi, and αw 6= −1. In case some orders are equal to -1 Eq. (11)takes the following formE[fi(X(i)(tk−1), θi)] =(2√2piσ)Ki Ni∑h=1θih{(√2piσ2)#(S−Sh)××[ ∏w∈S′h1αw + 1((Xˆw +√2piσ4)αw+1 − (Xˆw − √2piσ4 )αw+1)]××[ ∏w∈S′′hlnXˆw +√2piσ4Xˆw −√2piσ4]}(12)where S′h is the set of indexes {h′1, h′2, . . . , h′s} such that αh′ 6= −1 ∀h′ ∈ S′h,and S′′h is the set of indexes {h′′1, h′′2, . . . , h′′s} such that αh′′ = −1 ∀h′′ ∈ S′′h.If in the Eq. (8), mi is substituted with the expression (11) or (12),Eq. (8) becomes more tractable and can be optimized w. r. t. the param-eters Θ = (θ1, θ2, . . . , θN ) and σ. The values of the model’s parameters forwhich p(D|Θ) has a maximum are the most likely values giving the observedkinetics.3 The estimation errorsFor a zero-th order reaction∅ θ−→ AAθ−→ ∅6∆θ ={( ∂θ∂A0A0)2+( ∂θ∂AfAf)2} 12 (13)∂k∂Af=(Af −A0)tf − t0 (14)∂k∂A0=(A0 −Af )tf − t0 (15)For a first order reactionAθ−→ Pif A0 (Af ) is the initial (final) concentration of the molecules of thespecies A, the error that the algorithm makes in estimating θ is given bythe following∆θ ={( ∂θ∂A0A0)2+( ∂θ∂AfAf)2} 12 (16)Sinceθ =1tf − t0 (lnA0 − lnAf ) (17)Eq. 16 becomes∆θ ={( ∂θA0A0)2+( ∂θAfAf)2} 12 (18)For a second order reactionA+B θ−→ Pif A0 (Af ) and B0 (Bf ) are the initial (final) concentrations of the moleculesof species A and B, respectively, the error that the algorithm makes inestimating θ is given by the following∆θ ={( ∂θ∂A0A0)2+( ∂θ∂B0B0)2+( ∂θ∂AfAf)2+( ∂θ∂BfBf)2} 12 (19)where∂θ∂A0=1tf − t0{ 1(A0 −B0)2 lnAfB0A0Bf− 1A0 −B0BfA0}(20)∂θ∂B0=1tf − t0{− 1(A0 −B0)2 lnAfB0A0Bf+1A0 −B0AfB0}(21)7∂θ∂Af=1tf − t01A0 −B0B0Af(22)∂θ∂Bf=1tf − t01B0 −A0A0Bf(23)4 Case studiesAll the mathematical operations described above has been automated throughthe implementation of an algorithm in C++. We tested our algorithm onsome real case-study, whose rate coefficients were experimentally determinedand well documented in literature for simple chemical reactions and for morecomplex biochemical interactions. We show the ability of our algorithm ofobtaining reasonable estimates for the rate coefficients in the following cases(see Table 1):R1 Bromine molecule formation 2Br −→ Br2. The data referred to thisreaction are the bromine concentration time-points after a flash pho-tolysis of a mixture of bromine and SF6 with [Br2]/[SF6] = 3.2×10−2(experiment of Graff and Lang [11]).R2. In presence of an acid solution of phenol, the ion IO−3 is reduced toIO−2 by Br− accordingly to the reactionIO−3 + 2Br− + 2H+ −→ IO−2 +Br2 +H2OWith [IO−3 ] = 5×10−3M and [Br−] = 1×10−2M , Sharma and Guptaobtained the following data for a solution at 35 C with [C6H5OH] =2× 10−2M [11].R3. the first-order decomposition in azomethane CH3N2CH3 −→ CH3CH3+N2 ([13]);R4. the alkaline hydrolysis of ethyl nitrobenzoathe ([13]);R5. the gluthathione S-transferase-catalysed dehalogenation of haloaro-matic compounds with aromatic substrates ([15, 10]);R6. degradation of IkBa-UB in breast cancer cells (refer to [3] for a com-parison);R7. Nuclear factor (NF) kB translocation from cytoplasm into the nucleusin breast cancer cells (r. t. [3] for a comparison);R8. the glucose consumption rate in astrocytes [9].8Reaction k k σ(experimental) (estimate) (estimate)R1 (2.75± 1.00)× 108M−1s−1 (1.93± 1.28)× 108M−1s−1 3.2× 10−6R2 (13.8± 1.00)M−2s−2 (15.92± 6.98)M−2s−1 3.5× 10−1R3 (3.60± 2.16)× 10−4s−1 (0.80± 0.45)× 10−4s−1 5.36× 10−3R4 (8.1± 2.31)× 10−2s−1 (5.9± 1.5)× 10−2s−1 3.3× 10−3R5 (7.4± 0.5)× 10−4s−1 (2.60± 1.17)× 10−4s−1 4.1× 10−3R6 (5.56± 0.11)× 10−4s−1 (3.1± 2.2)× 10−4s−1 5× 10−3R7 (3.85± 2.67)× 10−4s−1 (1.3± 0.4)× 10−4s−1 8.5× 10−3R8 (1.33± 0.04)× 10−2s−1 (1.24± 0.78)× 10−2s−1 6× 10−4Table 1: The estimates of the rate coefficients compared with the experimentalvalues.5 ConclusionThe ability of this algorithm to estimate the expected value of the kineticrate constant significantly depends on the experimental error on the data. Ifthe time-course measurements of the reactant concentration are affected by ahigh level of noise, this technique estimate a rate constant very close to zero,indicating that there is no dynamic in the systems. Therefore the observedvariations in the concentration values are attributed exclusively to the noise.For an accurate estimate of the rate constants the procedure requires high-quality and high-resolution data, that at the time of writing are relativelydifficult to obtain, In the meantime experimental biology evolves to makesuch data more readily available, the development of inference techniquesallow identification of experimental design issues that need to be addressedto get more accurate data for more reliable estimate of rate constants.References[1] J. M. Bower, H. Bolouri (Eds.), Computational Modeling of Genetic andBiochemical Networks, Cambridge, Massachusetts, MIT Press, 2001.[2] E. . Crampin, P. E. McSharry, and S. Schnell, Extracting BiochemicalReaction Kinetics from Time Series Data, M. GH. Negoita et al. (Eds):KES 2004, LNAI 3214, p. 329-336, Springer-Verlag Berlin Heidelberg2004.[3] J. W. Fuseler, D. M. Merril, J. A. Rogers, M. B. Grisham and R. E.Wolf, Analysis and quatitation of NF-kB nuclear translocation in tu-mor necrosis factor alpha (TNF-α) activated vascular endothelial cells,Microsc. Microanal. 12, pp. 269-276, 2006.9[4] A. Golightly and D. J. Wilkinson, Bayesian sequential inference forstochastic kinetic biochemical network inference, Journal of Computa-tional Biology vol. 13, no. 3 pp. 838-851, 2006.[5] A. Golightly and D. J. Wilkinson, Bayesian inference for stochastickinetic models using a diffusion approximation, Biometrics, vol. 61, no.3, pp. 781-788, 2005.[6] T. B. L. Kirkwood, R. J. Boys, C. S. Gillespie, C. J. Proctor, D. P.Shanley, and D. J. Wilkinson, Toward an e-biology of ageing: inte-grating theory and data, Nature Reviews Molecular Cell Biology 4, pp.243-249, 2003.[7] H. Kitano (Ed.), Foundations of System Biology, Cambridge, Mas-sachusetts, MITOress, 2001.[8] N. Lalam and C. A. J. Klassen, Pseudo maximum likelihood estimationfor differential equations,[9] P. Lecca and M. Lecca, Molecular mechanism of glutamate-triggeredbrain glucose metabolism: a parametric model from FDG PET-scans,In Procs of BVAI 2007, To appear in LNCS, 2007.[10] I. Mancini, G. Guella, G. Chiasera and F. Pietra, Glutathione S-transferase catalysed dehalogenation of haloaromatic compounds whichlack nitro groups near the reaction centre, Tetrahedron Letters, 39, pp.1611-1614, 1998.[11] C. R. Metz, Chimica fisica, McGraw-Hill, 1994.[12] I. Nachman, A. Regev and N. Friedman, Inferring quantitative modelsof regulatory networks from expression data, Bioinformatics, vol. 20, no.1, pp. 248-256, 2004.[13] V. Nummert and M. Piirsalu, Kinetic study of hydrolysis of benzoates,J. Chem. Soc., Perkin Trans. 2, pp. 583-594, 2000.[14] S. Reinker, R. M. Altman and J. Timmer, Parameter estimation instochastic biochemical reactions, IEE Proc. Syst. Biol, vol. 153, no. 4,July 2006.[15] B. Rossi, I. Mancini and G. Guella, Kinetic studies of glutathione S-transferase-catalysed process through on- line liquid chromatography-electrospray mass spectrometry, Philosophical Magazine, 84(13-16), pp.1373-1382, 2004.[16] J. Yu, V. A. Smith, P. P. Wang, a. J. Hartemink and E. D. Jarvis,Advances to Bayesian network inference for generating causal networks10from observational biological data, Bioinformatics vol 20, no. 18, pp.3594 - 3603, 2004.[17] R. J. Boys, D. J. Wilkinson and T. B. Kirkwood, Bayesian inference fora discretely observed stochastic kinetic model, currently in submission(available at http://www.mas.ncl.ac.uk/∼ndjw1/docs/BWK04.pdf).11

http://eprints.biblio.unitn.it/1724/1/TR%2D20%2D2007.pdf

A new model for kinetic parameter estimation in biochemical reactions.

Abstract

Similar works

Full text

Available Versions

University of Trento : UNITN-eprints

Unitn-eprints Research