Stochastic Parameter Estimation in Biochemical Signalling Pathways

It is common when modelling biochemical networks to use qualitative information such as the general ODE model structure so as to proceed in parameter estimation while at the same time retaining the basic model structure the best represents the biochemical process governing the cell. This is not the case however when the population of the available molecules from each of the participating species is very small (small copy number) deeming necessary the introduction of complex stochastic modelling techniques that make use of chemical master equations to simulate the trajectories of the states (species concentration) of the system [7]. Gene expression is stochastic by nature [7][5] and as a consequence gene regulatory and signal transduction networks follow a similar behaviour. Most importantly, a large number of gene expression data sets examined in yeast, mouse and human cells follow a Pareto-like distribution model skewed by many low-abundance transcripts, covering a large variety of eukaryotic cells [2]. It is therefore apparent that a stochastic modelling strategy should be structure so as to accommodate the specific needs of the system. Previous studies have focused on capturing the stochastic nature of the system by assuming that the noise terms follow a Gaussian style distribution [6][4], as do the priors for the parameters used in the models [6]. It follows that the systems behaviour could be more accurately represented by a system whose state trajectories (due to noise) have deviations of arbitrary, non-Gaussian, distribution in nature. In fact, the state trajectories obtained by a stochastic simulation algorithm might be thought of as the stochastic version of the trajectories that would be obtained by solving the reaction rate equations of a traditional deterministic system [7]. The time step differs in that in the stochastic case is exact and not a finite approximation to some infinitesimal dt [7]. In this paper a novel stochastic parameter estimation technique is introduced where the same ODE model of the biochemical network is used as in the deterministic case, under the assumption that it is corrupted by additive noise of arbitrary, non-Gaussian distribution. This is now a bounded dynamic stochastic system and as such it is possible to evaluate the joint probability density function (PDF) of the output that is proven to be dynamically linked to the systems parameters[9]. The purpose of the method is to select those parameters that drive the joint output PDF to follow a target stochastic distribution as close as possible and generally direct the system into a state of lower randomness. The method could be thought of as a way of shaping the joint output PDF so as to approximate as close as possible the given target [9]. The joint model PDF is approximated using B-Splines with the weights dynamically linked to the original parameters of the ODE model. This is an optimization problem where the objective is to minimize the integral square error between the two PDFs. A key advantage is that the method does not depend on time and only the stochastic properties of the system are considered conveyed by the model output PDF. The objective function is formulated as: J( ) = (y, ) g(ŷ) ( ) 2 dy