Regularized maximum likelihood estimation of sparse stochastic monomolecular biochemical reaction networks

Highlights

• A sparse parameter matrix estimation method is proposed for a stochastic monomolecular biochemical reaction network system.

• We considers stochastic monomolecular reaction systems for which an exact analytical solution of the CME can be derived.

• The estimation method incorporates the closed-form solution into a regularized MLE for which model complexity is penalized.

• Simulation results are provided to verify performance improvement of regularized MLE over LSE.

Abstract

A sparse parameter matrix estimation method is proposed for identifying a stochastic monomolecular biochemical reaction network system. Identification of a reaction network can be achieved by estimating a sparse parameter matrix containing the reaction network structure and kinetics information. Stochastic dynamics of a biochemical reaction network system is usually modeled by a chemical master equation (CME) describing the time evolution of probability distributions for all possible states. This paper considers closed monomolecular reaction systems for which an exact analytical solution of the corresponding chemical master equation can be derived. The estimation method presented in this paper incorporates the closed-form solution into a regularized maximum likelihood estimation (MLE) for which model complexity is penalized. A simulation result is provided to verify performance improvement of regularized MLE over least-square estimation (LSE), which is based on a deterministic mass-average model, in the case of a small population size.

Introduction

Stochastic dynamics of biological systems have lately received increased attention from researchers in the field of biological engineering (Raj and van Oudenaarden, 2009). In the past, such studies were greatly hampered by lack of qualitative measurement data. Nowadays, quantitative but noisy data can be obtained using the microarray technology. Recent developments in sensing techniques that can provide real-time observations of intrinsic stochastic dynamics at small length scales have motivated many scientific investigations (Raj et al., 2010, Xie et al., 2008). One such sensing technique for biological systems is bio-imaging using fluorescent proteins (García-Parajó et al., 2001). By grafting a fluorescent protein into gene expression, proteins and mRNA expressions originating from targeted DNA can be detected quantitatively in real time. Specifically, yellow fluorescent protein (yfp) (Elowitz et al., 2002, Li and Xie, 2011, Yu et al., 2006) has been widely used for detecting changes with single-macromolecule sensitivity in individual live cells.

In the real-time data reported in the aforementioned papers, strongly stochastic behavior has been observed. For example, bursts of transcribed protein molecules from the cell, controlled by an identical messenger RNA molecule, have different copy numbers (Elowitz et al., 2002). The total population number of species in the system within the detectable range is small, ranging from tens to thousands of copies. Stochastic dynamics of systems with discrete states can be modeled by the chemical master equation (CME) (Feinberg, 1979, Fichthorn and Weinberg, 1991),

$$\frac{\partial P\left(\sigma ,t\right)}{\partial t}={\displaystyle {\displaystyle \sum }_{{\sigma }^{\prime }}}W\left({\sigma }^{\prime },\sigma \right)P\left({\sigma }^{\prime },t\right)-{\displaystyle {\displaystyle \sum }_{{\sigma }^{\prime }}}W\left(\sigma ,{\sigma }^{\prime }\right)P\left(\sigma ,t\right)$$

where$P\left(\sigma ,t\right)$ is the probability of the system being in discrete state $\sigma$ at time t, and $W\left({\sigma }^{\prime },\sigma \right)$is the transition rate from state ${\sigma }^{\prime }$ to state $\sigma$. The CME describes the time evolution of the probability distribution among all possible configurations, and can be written as (MacNamara et al., 2008, Munsky and Khammash, 2006)

$$\frac{dP\left(t\right)}{dt}=A\left(t;\theta \right)P\left(t\right)$$

where $P\left(t\right)$ is the state vector containing all the state probability variables and $A\left(t;\theta \right)$ is a matrix containing all the transition rate constants with dependence on the model parameter vector $\theta$ related with physicochemical phenomena, e.g., biochemical reactions. Many numerical algorithms have been developed for solving the matrix ordinary differential Eq. (2), which can be divided into direct methods (MacNamara et al., 2008, Munsky and Khammash, 2006) and indirect methods (Gibson and Bruck, 2000, Gillespie, 1977). Direct methods, such as the finite state projection (FSP) algorithm, attempt to evaluate the matrix exponential directly. In practice, given the large size of the state space, indirect methods that use stochastic simulation algorithms (SSAs) to generate approximate probability distributions have been more popular.

A reaction network is one of fundamental models for describing biological mechanisms. Different types of biological reaction networks exist and interact to accomplish needed functions including gene regulatory networks, metabolic networks, and signaling networks. Network identification is an important aspect of studying biological network systems. Given time series experimental data from sensors, a parameter estimation method can be used to identify the parameters of a given reaction network model. Typically, the estimation is formulated to find parameter values minimizing the difference between the experimental data and their model predictions. Most of the literature has employed least-squares estimation (LSE) approaches, which fit stochastic data to a deterministic continuum model (Gennemark and Wedelin, 2009).

The LSE method based on deterministic models, however, can provide poor parameter estimates for highly stochastic systems (Tian et al., 2007). Several stochastic parameter estimation methods based on the stochastic differential equation of type (1) have shown improved estimation performance (Munsky et al., 2012). These methods attempt to solve for the probability density functions (PDFs) of the CME and use them for estimation. Previous studies employed the moment-based method (Munsky et al., 2009, Zechner et al., 2012), the Bayesian method (Boys et al., 2008, Golightly and Wilkinson, 2011, Lillacci and Khammash, 2010), the maximum likelihood estimation (MLE) method (Daigle et al., 2012, Tian et al., 2007), and the density function distance (DFD) method (Lillacci and Khammash, 2013, Poovathingal and Gunawan, 2010). However, these published methods are based on approximated PDFs of the CME. In many cases, PDFs are approximated using the SSA approach (Gillespie, 2007), which typically demands a very large number of simulations to be performed for an accurate estimation.

In many cases, the exact network structure of a biological system is not known. An important problem in biological systems is to use measurement data to identify the topology of reaction networks and estimate associated parameters at the same time. Craciun and Pantea, 2008 discussed the issue of identifiability of reaction networks with the fact that there might exist more than one parameter set that fits the measurement data with the same level of accuracy. To cope with such challenges, model discrimination/invalidation of reaction networks has been studied (Conradi et al., 2005, Kremling et al., 2004). Non-uniqueness of an estimation algorithm can be overcome by using the principle of parsimony and choosing simple models over complex models (Jefferys and Berger, 1991). An approach to reduce model complexity is to use penalization term for the number of parameters such as ridge and $\ell$₁ regularization (Hesterberg et al., 2008).

Taking the above issues into consideration altogether, this paper considers stochastic monomolecular reaction systems in a sparse parameter matrix estimation problem. Stochastic monomolecular reaction system, described by an exact probability distribution solution of the CME (Jahnke and Huisinga, 2007), can represent realistic gene regulatory networks or metabolic networks. The exact solution enables formulation of a regularized MLE method rather than LSE. Improved performance over the LSE method is verified with a simulation study involving a small scale reaction network system. Applicability of the proposed MLE method to a more stochastic and larger scale reaction network system is also tested.