A Convex Approach to Steady State Moment Analysis for Stochastic Chemical Reactions

TL;DR

This paper introduces a convex optimization method to accurately compute steady state moments of molecular counts in stochastic chemical reactions, aiding in biomolecular circuit design.

Contribution

It presents a novel convex semi-algebraic set framework for bounding moments, enabling precise predictions with semidefinite programming.

Findings

Accurately predicts mean and variance of protein copy numbers.

Provides a rigorous convex optimization approach for stochastic moment analysis.

Demonstrates effectiveness on a protein dimerization example.

Abstract

Model-based prediction of stochastic noise in biomolecular reactions often resorts to approximation with unknown precision. As a result, unexpected stochastic fluctuation causes a headache for the designers of biomolecular circuits. This paper proposes a convex optimization approach to quantifying the steady state moments of molecular copy counts with theoretical rigor. We show that the stochastic moments lie in a convex semi-algebraic set specified by linear matrix inequalities. Thus, the upper and the lower bounds of some moments can be computed by a semidefinite program. Using a protein dimerization process as an example, we demonstrate that the proposed method can precisely predict the mean and the variance of the copy number of the monomer protein.

Tables

Table 1. Table 1: Transition rates w i ​ ( x ) subscript 𝑤 𝑖 𝑥 w_{i}(x) in Fig.  2

Reaction $i$	Transition rate $w_i\left(x\right)$
$1$	$k_1{D}_T$
$2$	$k_{2}x$
$3$	$k_{3}x\left(x-1\right)$

Table 2. Table 2: Parameters of chemical reactions in Fig.  2

​​​Parameter​​​	Value	Meaning
$k_1$	$0.7{\min }^{-1}$​​	​​​​​​​Transcription and translation rate​​
$k_2$	$\ln \left(2\right)/20{\mathrm{copy}}^{-1}\cdot {\min }^{-1}$	​​​​Protein (monomer) degradation​​
$k_3$	$0.1{\mathrm{copy}}^{-1}\cdot {\min }^{-1}$	Dimerization rate
$D_T$	$50\mathrm{copy}$	Total copy of DNA

Equations

$$\frac{dP{}_{\mathbf{x}}(t)}{dt}=\sum _{i=1}^r\left\{w_i(\mathbf{x}-{\mathbf{s}}_i)P{}_{\mathbf{x}-{\mathbf{s}}_i}(t)-w_i(\mathbf{x})P{}_{\mathbf{x}}(t)\right\},$$

$$w_i(\mathbf{x})=\underset{\Delta t\to 0}{\lim }\frac{P_{\mathbf{x}+{\mathbf{s}}_i}(t+\Delta t)-P_{\mathbf{x}}(t)}{\Delta t}(i=1,2,\cdots ,r).$$

$$m_{\mathbf{\alpha }}:=\mathbb{E}\left[\prod _{j=1}^n{x}_j^{{\alpha }_j}\right]=\sum _{x_1=0}^{\infty }\sum _{x_2=0}^{\infty }\cdots \sum _{x_n=0}^{\infty }\prod _{j=1}^n{x}_j^{{\alpha }_j}{P}_{\mathbf{x}}(t),$$

$$\frac{d}{dt}\left\{\sum _{\mathbf{x}}\left(\prod _{j=1}^n{x}_j^{{\alpha }_j}\right)P_{\mathbf{x}}(t)\right\}$$

$$=\sum _{\mathbf{x}}\prod _{j=1}^n{x}_j^{{\alpha }_j}\sum _{i=1}^r\left\{w_i(\mathbf{x}-{\mathbf{s}}_i)P_{\mathbf{x}-{\mathbf{s}}_i}(t)-w_i(\mathbf{x})P_{\mathbf{x}}(t)\right\},$$

$$\sum _{\mathbf{x}}:=\sum _{x_1=0}^{\infty }\sum _{x_2=0}^{\infty }\cdots \sum _{x_n=0}^{\infty }.$$

$$\frac{d}{dt}{m}_{\mathbf{\alpha }}=\sum _{{\beta }_1}\sum _{{\beta }_2}\cdots \sum _{{\beta }_n}{A}_{\mathbf{\alpha },\mathbf{\beta }}{m}_{\mathbf{\beta }}(\mathbf{\alpha }\in {\mathbb{N}}_0^n),$$

$$\sum _{i=1}^r$$

$$\cdots (x_n+s_{in}{)}^{{\alpha }_n}-x_1^{{\alpha }_1}{x}_2^{{\alpha }_2}\cdots x_n^{{\alpha }_n}\}{w}_i(\mathbf{x}).$$

$$0\le \mathrm{deg}(m_{\mathbf{\beta }})\le \mathrm{deg}(m_{\mathbf{\alpha }}),$$

$$0\le \mathrm{deg}(m_{\mathbf{\beta }})\le \mathrm{deg}(m_{\mathbf{\alpha }})+1.$$

$$0=\sum _{{\beta }_1}\sum _{{\beta }_2}\cdots \sum _{{\beta }_n}{A}_{\mathbf{\alpha },\mathbf{\beta }}{m}_{\mathbf{\beta }}\left(\mathbf{\alpha }\in {\mathcal{A}}_{\mu }\subset {\mathbb{N}}_0^n\right),$$

$$\mathcal{A}_{\mu}\mathrel{\mathop{:}}=\left\{\bm{\alpha}\in\mathbb{N}_{0}^{n}~{}\bigg{|}~{}0\leq\sum_{i}^{n}\alpha_{i}\leq\mu\right\}.$$

$$(H{)}_{p,q}:=\mathbb{E}\left[{\mathbf{x}}^p{\left({\mathbf{x}}^q\right)}^{\mathrm{T}}\right],$$

$$(H_k{)}_{p,q}:=\mathbb{E}\left[x_k{\mathbf{x}}^p{\left({\mathbf{x}}^q\right)}^{\mathrm{T}}\right](k=1,2,\cdots ,n).$$

$$\mathrm{deg}(m_{\mathbf{\beta }})=\{\begin{array}{ll}p+q& \mathrm{for}(H{)}_{p,q}\\ p+q+1& \mathrm{for}(H_k{)}_{p,q}\end{array}.$$

$$H^{(\gamma_{1})}\!(\{m_{\bm{\beta}}\})\!\mathrel{\mathop{:}}=\!\left[\begin{array}[]{cccc}\!\!\!(H)_{0,0}\!&\!(H)_{0,1}\!&\!\cdots\!&\!(H)_{0,\gamma_{1}}\!\!\!\\ \!\!\!(H)_{1,0}\!&\!(H)_{1,1}\!&\!\cdots\!&\!(H)_{1,\gamma_{1}}\!\!\!\\ \!\!\!\vdots\!&\!\vdots\!&\!\ddots\!&\!\vdots\!\!\!\\ \!\!\!(H)_{\gamma_{1},0}\!&\!(H)_{\gamma_{1},1}\!&\!\cdots\!&\!(H)_{\gamma_{1},\gamma_{1}}\!\!\!\end{array}\right],$$

$$H_{k}^{(\gamma_{2})}\!(\{m_{\bm{\beta}}\})\!\mathrel{\mathop{:}}=\!\left[\begin{array}[]{cccc}\!\!\!(H_{k})_{0,0}\!\!&\!\!(H_{k})_{0,1}\!\!&\!\!\cdots\!\!&\!\!(H_{k})_{0,\gamma_{2}}\!\!\!\\ \!\!\!(H_{k})_{1,0}\!\!&\!\!(H_{k})_{1,1}\!\!&\!\!\cdots\!\!&\!\!(H_{k})_{1,\gamma_{2}}\!\!\!\\ \!\!\!\vdots\!\!&\!\!\vdots\!\!&\!\!\ddots\!\!&\!\!\vdots\!\!\!\\ \!\!\!(H_{k})_{\gamma_{2},0}\!\!&\!\!(H_{k})_{\gamma_{2},1}\!\!&\!\!\cdots\!\!&\!\!(H_{k})_{\gamma_{2},\gamma_{2}}\!\!\!\end{array}\right],$$

$$\gamma_{1}=\left\{\begin{array}[]{c}(\mu+1)/2~{}~{}~{}({\rm if}~{}\mu~{}{\rm is~{}odd})\\ ~{}~{}~{}~{}~{}\mu/2~{}~{}~{}~{}~{}~{}({\rm if}~{}\mu~{}{\rm is~{}even})\end{array},\right.$$

$$\gamma_{2}=\left\{\begin{array}[]{c}(\mu-1)/2~{}~{}~{}({\rm if}~{}\mu~{}{\rm is~{}odd})\\ ~{}~{}~{}~{}~{}\mu/2~{}~{}~{}~{}~{}~{}({\rm if}~{}\mu~{}{\rm is~{}even})\end{array}.\right.$$

$$H^{({\gamma }_1)}(\{m_{\mathbf{\beta }}\})$$

$$H_k^{({\gamma }_2)}(\{m_{\mathbf{\beta }}\})$$

$$0=\sum _{{\beta }_1}\sum _{{\beta }_2}\cdots \sum _{{\beta }_n}{A}_{\mathbf{\alpha },\mathbf{\beta }}{m}_{\beta }^{\ast }\left(\mathbf{\alpha }\in {\mathcal{A}}_{\mu }\subset {\mathbb{N}}_0^n\right),$$

$$H^{({\gamma }_1)}(\{m_{\beta }^{\ast }\})⪰0,$$

$$H_k^{({\gamma }_2)}(\{m_{\beta }^{\ast }\})⪰0(k=1,2,\cdots ,n),$$

$$\begin{array}[]{c}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm min}~{}~{}~{}~{}f(m_{\bm{\alpha}})\\ {\rm subject~{}to}~{}~{}~{}~{}~{}(\ref{eq:subject}),\end{array}$$

$$\begin{array}[]{c}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm max}~{}~{}~{}~{}f(m_{\bm{\alpha}})\\ {\rm subject~{}to}~{}~{}~{}~{}~{}(\ref{eq:subject}).\end{array}$$

$$f(m_{\mathbf{\alpha }})=\mathbb{E}[x_i]=m_{{\mathbf{e}}_i}$$

$$f(m_{\mathbf{\alpha }})$$

$$=m_{2{\mathbf{e}}_i}-m_{{\mathbf{e}}_i}^2.$$

$$\max$$

$$\mathrm{subject}\mathrm{to}$$

$$\displaystyle\left[\begin{array}[]{cc}m_{2\bm{e}_{i}}-V&m_{\bm{e}_{i}}\\ m_{\bm{e}_{i}}&1\end{array}\right]\succeq O.$$

$$m_{\alpha }:=\mathbb{E}\left[x^{\alpha }\right]=\sum _{x=0}^{\infty }{x}^{\alpha }{P}_x(t),$$

Full text

1 Introduction

A grand challenge of synthetic biology is to build and control layers of artificial biomolecular reaction networks that enable complex tasks by utilizing naturally existing reaction resources in micro-scale organisms such as E. coli. In control engineering community, many theoretical tools were developed over the last decade to enable model-guided design of biomolecular circuits including a fold-change detector [1], oscillators [2, 3, 4, 5], an event detector [6], and a light-controlled feedback circuit for setpoint control [7, 8]. Other successfully implemented examples include logic gates [9], a genetic memory [10], and a communication system between cells [11], to name a few. Given the elementary modules of artificial bimolecular networks, the next step stone is to assemble the circuit parts to build multiple layers of biomolecular networks that can provide practical functions.

As is the case with mechanical and electrical engineering, guaranteeing the performance of individual circuit modules at high precision is a key to successfully building large-scale and robust biomolecular circuits. Biosystems engineering, however, has a unique challenge, in contrast with other fields of engineering, that the signal-to-noise ratio is so low that the heterogeneity of circuit states between biological cells is almost impossible to avoid. A major source of the heterogeneity is the stochastic chemical reactions that are caused by the low copy nature of reactive molecules in a small-volume reactor, that is, a biological cell. In other words, the events of molecular collision that fire reactions are more accurately captured by stochastic process rather than continuous and deterministic process governed by ordinary differential equations (ODEs). This situation raises a strong need for the theoretical tools that can rigorously certify the performance of stochastic chemical reactions using statistical norms such as mean and variance of molecular copy numbers.

The random fluctuation of the copy number of molecules can be regarded as a Markov process that follows the chemical master equation (CME) [12]. Despite a linear ODE, the CME is often hard to solve analytically because of the infinite dimensionality of the equation. A typical solution to this problem is to use Monte Carlo simulations [13]. This approach, however, requires high computational time. Moreover, a strict error bound is hard to obtain. To resolve these issues, the finite state project [14] provides a systematic way to truncate the equation with an error bound. Although useful for analyzing transient dynamics, the FSP still suffers from high dimensionality if one wants to evaluate the probability distribution near the steady state, which is often the operation point of interest.

A more direct approach to evaluating the statistics of biomolecular circuits is to use the moment equation, an infinite dimensional ODE derived from the CME. Recently developed moment closure approaches [15, 16] allow for approximately solving the ODE by order reduction technique. The central idea is to approximate the high order moments in the equation using the lower order moments to derive a finite order closed ODE. Using this method, we can efficiently and directly compute the statistics of the molecular copy numbers of biomolecular circuits. Although appealing, this approach requires assuming the underlying probability distribution, which is not apriori in practice. As a result, biocircuit designers have to engineer biomolecular circuits based on approximated statistics with unknown precision.

The objective of this paper is to propose an algebraic approach to rigorously quantifying the steady state statistics of stochastic biomolecular reactions without assuming underlying probability distributions. Specifically, we formulate a semidefinite program that computes the upper and the lower bounds of statistical values of molecular copy numbers using the moment approach [Lasserre2009]. During the review process, the authors were informed that the same approach has been explored independently by two other research groups [Ghusinga2016, Kuntz2017]. In comparison with these works, this paper presents an optimization algorithm that directly computes the upper bound of the second order central moment, or variance, without running multiple optimizations for computing raw moments, allowing for direct and tighter quantification of the upper bound. We illustrate the proposed method using a stochastic protein dimerization process.

The organization of this paper is as follows. In Section II, we introduce the CME and derive the moment equation. Then, the optimization problem for moment computation is formulated in Section III. In section IV, we analyze the mean and the variance of a stochastic dimerization process. Finally, Section V summarizes our findings and concludes this paper.

The following notations are used in this paper. $\mathbb{N}_{0}:=\{0,1,2,\cdots\}$. $\mathbb{N}_{+}:=\{1,2,3,\cdots\}$. $\mathbb{R}_{+}:=\{x\in\mathbb{R}~{}|~{}x\geq 0\}$. ${\rm deg}(p(x))$ is the degree of the polynomial $p(x)$. For multivariate polynomials $p({\bm{x}})=\prod_{j=1}^{n}x_{j}^{p_{j}}$, ${\rm deg}(p({\bm{x}})):=\sum_{j=1}^{n}p_{j}$.

2 Moment Dynamics of Stochastic Chemical Reactions

2.1 Model of stochastic chemical reactions

We consider a chemical reaction network consisting of $n\in\mathbb{N}_{+}$ species of molecules and denote the copy number of each molecular species by $\bm{x}=[x_{1},x_{2},\cdots,x_{n}]^{\mathrm{T}}\in\mathbb{N}_{0}^{n}$. There are $r$ types of reactions in the reaction network, and let $\bm{s}_{i}=[s_{i1},s_{i2},\cdots,s_{in}]^{\mathrm{T}}~{}\in\mathbb{Z}^{n}$ be the increment of the molecular copy $\bm{x}$ by the $i$-th chemical reaction ($i=1,2,\cdots,r$).

Chemical reactions are caused by collisions of molecules. When the reaction volume is sufficiently large, the copy number of molecules is so large that the collision event occurs almost continuously in time. On the other hand, in small reactor systems such as microbes, the collision events tend to be stochastic as the copy number of molecules is small. This results in the stochastic fluctuation of $\bm{x}$. It is known that stochastic chemical reactions can be modeled by a Markov process. More specifically, we define $P_{\bm{x}}(t)$ as the probability that there are $\bm{x}$ molecules at time $t$. The stochastic chemical reactions can be modeled by the following Chemical Master Equation (CME) [12],

[TABLE]

where $w_{i}(\bm{x})$ is the transition rate associated with the $i$-th chemical reaction and is defined by

[TABLE]

The transition rate $w_{i}(\bm{x})$ is characterized by the rate of the reaction $i~{}(i=1,2,\cdots,r)$ We assume that all of the $r$ reactions are elementary as non-elementary reactions can always be decomposed into elementary reactions. For elementary reactions, the transition rate $w_{i}({\bm{x}})$ is a polynomial of $x_{i}$’s, and we utilize this fact in the following subsection to derive the equations of stochastic moments.

As ${\bm{x}}$ is the vector of non-negative positive integers $\mathbb{N}_{0}^{n}$, the CME (1) is a linear but an infinite dimensional ODE with respect to $P_{\bm{x}}(t)$. An analytic solution to the CME is, thus, hardly obtained except for some simple examples. In the next subsection, we derive the equations of moment dynamics to directly compute the statistical values without computing the distribution $P_{\bm{x}}(t)$.

2.2 The dynamics of moments

In this section, we derive the model of moment dynamics based on the CME (1). First, we define a stochastic moment by

[TABLE]

where $\bm{\alpha}\mathrel{\mathop{:}}=[\alpha_{1},\alpha_{2},\cdots,\alpha_{n}]\in\mathbb{N}_{0}^{n}$. We refer to the sum of the entries of $\bm{\alpha}$ as the order of the moment and denote by ${\rm deg}(m_{\bm{\alpha}}):=\sum_{j}\alpha_{j}$ with a little abuse of notation.

The stochastic moments carry all the information of the probability distribution, and some of the low order moments are useful to define the design specification of stochastic biomolecular circuits. Figure 1 illustrates an example of the distribution of a molecular copy number $x$ of some biomolecular reaction. In this case, the analytic form of the distribution is not known, and thus it is more appropriate to specify the design goal by statistical norms such mean $\mathbb{E}[x]=m_{1}$ and the standard deviation $\sqrt{\mathbb{E}[(x-\mathbb{E}[x])^{2}]}=\sqrt{m_{2}-m_{1}^{2}}$ of the distribution using the moments.

To derive the dynamics of stochastic moments, we multiply the model (1) by $\prod_{j}x_{j}^{\alpha_{j}}$ and calculate the expected value by taking sum over $\mathbb{N}_{0}^{n}$ to obtain the following differential equation.

[TABLE]

where we use $\sum_{\bm{x}}$ to abbreviate the sum over the positive orthant.

[TABLE]

As $w_{i}({\bm{x}})$ is a polynomial of $x_{i}~{}(i=1,2,\cdots,n)$, we can rewrite the equation (4) as

[TABLE]

where $A_{\bm{\alpha},\bm{\beta}}$ are the coefficients of the monomials of the polynomial

[TABLE]

Note that the range of the summation in (6) depends on the degree of the polynomial (7). In many cases, the higher order moments are necessary to characterize the low order moments, i.e., there are moments $m_{\bm{\beta}}$ satisfying ${\rm deg}(m_{\bm{\beta}})\!>\!{\rm deg}(m_{\bm{\alpha}})$ in the right-hand side of (6). More specifically, when the reactions are elementary, it is reasonable to assume that the polynomial order of the transition rates $w_{i}(\bm{x})$ are at most two since most practical reactions are unimolecular or bimolecular reactions [17]. When ${\rm deg}(w_{i}({\bm{x}}))=0$ or $1$, $m_{\bm{\beta}}$ in the right-hand side of equation (6) satisfies

[TABLE]

implying that the moment equation is closed, i.e, the dynamics of the $i$-th order moment can be written by the $i$-th or lower order moments. On the other hand, when one or more $w_{i}({\bm{x}})$’s are quadratic, i.e., ${\rm deg}(w_{i}({\bm{x}}))=2$,

[TABLE]

Thus, the moment equation (6) becomes an infinite dimensional coupled linear ODE, which is hard to solve analytically.

3 Moment Analysis of Stochastic Chemical Reactions

3.1 Truncation of the moment equation

When we design stochastic chemical reactions using a genetic circuit, the specification of the circuit is often given by a set of statistical constraints such as the mean and the variance of the copy numbers of the molecules. In this paper, we consider the following problem.

Problem 1. For a given CME (1), find the lower and the upper bounds of the mean and the variance of molecular copy number $\bm{x}$ at steady state.

To solve this problem, let $dm_{\bm{\alpha}}/dt=0$ in the equation (6), and consider the subset of linear equations by limiting the equations up to the $\mu$-th order moments. Specifically,

[TABLE]

where the set $\mathcal{A}_{\mu}$ is defined by

[TABLE]

In what follows, we refer to $\mu$ as the truncation order. When ${\rm deg}(w_{i}({\bm{x}}))=0$ or $1$, the number of equations is the same as the number of variables (see (8)), and thus we can compute the unique solution (if the equations are not degenerated). In the more general case where the transition rates $w_{i}(\bm{x})$ are quadratic, the right-hand side of (10) depends on the higher order moments, i.e., ${\rm deg}(m_{\bm{\beta}})=\mu+1$, as shown in (9) Consequently, the equations become underdetermined, and we can only conclude that the solution lies on a certain hyperplane unless we consider an infinite number of equations.

To further specify the solution space of the moments $m_{\bm{\beta}}$, we utilize the fact that $m_{\bm{\beta}}$ must be the moments of some (probability) measure defined on the positive orthant. In particular, we use a so-called moment condition to constrain $m_{\bm{\beta}}$ and formulate a semidefinite program solving Problem 1.

3.2 Necessary condition for $m_{\bm{\beta}}$ to be moments

Let $\bm{x}^{p}$ be the vector that consists of all of the monomial bases satisfying ${\rm deg}(\prod_{j}x_{j}^{p_{j}})=p$, and $(H)_{p,q}$ be a matrix of the form

[TABLE]

where $\mathbb{E}[\bm{x}^{p}(\bm{x}^{q})^{\mathrm{T}}]$ represents the entry-wise expected value of the matrix $\bm{x}^{p}(\bm{x}^{q})^{\mathrm{T}}$. In a similar manner, we define $(H_{k})_{p,q}$ by

[TABLE]

It should be noted that the entries of the matrices $(H)_{p,q}$ and $(H_{k})_{p,q}$ can be written with $m_{\bm{\beta}}$ satisfying

[TABLE]

Using the matrices $(H)_{p,q}$ and $(H_{k})_{p,q}$, we define the following block Hankel matrices $H^{(\gamma_{1})}(\{m_{\bm{\beta}}\})$ and $H_{k}^{(\gamma_{2})}(\{m_{\bm{\beta}}\})~{}(k=1,2,\cdots,n)$.

[TABLE]

[TABLE]

where $\gamma_{1}$ and $\gamma_{2}$ are determined as follows.

[TABLE]

[TABLE]

Using $H^{(\gamma_{1})}(\{m_{\bm{\beta}}\})$ and $H_{k}^{(\gamma_{2})}(\{m_{\bm{\beta}}\})$, the following proposition provides linear matrix inequality (LMI) conditions that the moments of some (probability) measure on $\mathbb{R}_{+}^{n}$ must satisfy.

Proposition 1. [18] Consider a sequence $(m_{\bm{\beta}})_{\bm{\beta}\in\mathcal{A}_{\mu+1}}$ with the set $\mathcal{A}_{\mu+1}$ defined by (11). The sequence $(m_{\bm{\beta}})_{\bm{\beta}\in\mathcal{A}_{\mu+1}}$ constitutes moments of some measure defined on the positive orthant $\mathbb{R}_{+}^{n}\mathrel{\mathop{:}}=\{\bm{x}\in\mathbb{R}^{n}~{}|~{}x_{k}\geq 0~{}(k=1,2,\cdots,n)\}$ only if

[TABLE]

The conditions (19) and (20) are called a moment condition and localized moment conditions, respectively. Combining Proposition 1 with (10), the following theorem specifies semi-algebraic conditions that the moments of the stochastic chemical reactions must satisfy.

Theorem 1.  Consider the stochastic chemical reaction modeled by the equation (1). Let $m_{\bm{\beta}}^{*}$ denote the steady state moments of the random variable ${\bm{x}}$. For a given truncation order $\mu$, the moments $m_{\bm{\beta}}^{*}$ satisfy the following conditions.

[TABLE]

where $A_{\bm{\alpha},\bm{\beta}}$ are the coefficients of the polynomial (7).

Theorem 1 implies that the moments of the stochastic chemical reaction governed by (1) lie in the semi-algebraic set (21). Thus, Problem 1 can be recast as a relaxation problem over the convex semi-algebraic set. In particular, we show, in the following section, that the bounds of the mean and the variance can be computed by semidefinite programming.

Remark 1.  In the special case of $n=1$ random variable, the conditions (19) and (20) with $\gamma_{1},\gamma_{2}\rightarrow\infty$ become both necessary and sufficient for the sequence $(m_{\beta})_{\beta\in\mathbb{N}_{0}}$ to be the moments of some measure on $\mathbb{R}_{+}$ [18]. The associated problem is called Stieltjes moment problem named after the renowned mathematician Thomas Stieltjes. Similar conditions are available for the moments of measures defined on various domains such as $[0,1]$, $\mathbb{R}$ and semi-algebraic sets (see [18, 19] for example).

3.3 Semidefinite programming for mean and variance computation

Theorem 1 asserts that the lower and the upper bounds of the steady state statistics, say $f(m_{\bm{\alpha}})$, can be calculated by solving the following optimization problems, respectively.

[TABLE]

[TABLE]

Note that the constrains are convex. In what follows, we show that the bounds of mean and variance, the two most important statistics of stochastic chemical reactions, can be formulated as semidefinite programs.

Mean The objective function $f(m_{\bm{\alpha}})$ is obviously linear, and thus the problem falls into the class of semidefinite programming. For example, if one wants to compute the mean of the copy number of the $i$-th molecule, the objective function should be

[TABLE]

with $\bm{e}_{i}$ representing the row vector whose $i$-th entry is 1 and 0 elsewhere.

In general, the gap between the upper and the lower bounds decreases monotonically with the increase of the truncation order $\mu$. Thus, there is a tradeoff between computational cost and the conservativeness of the bounds.

Upper bound of variance   By definition, the variance of the copy number of the $i$-th molecule is given by

[TABLE]

The second term of the objective function $f(m_{\bm{\alpha}})=m_{2\bm{e}_{i}}-m_{\bm{e}_{i}}^{2}$ is quadratic, and the optimization problem (23) is seemingly not semidefinite programming. However, using the Schur complement, we can convert the optimization (23) into an equivalent semidefinite programming as follows.

Theorem 2.  The solution of the optimization problem (23) with the objective function (25) is the same as that of the following optimization problem.

[TABLE]

The proof of this theorem can be found in Appendix A.

Lower bound of variance The lower bound of variance can be also computed by semidefinite programming with additional relaxation. Specifically, we first compute the upper bound of the mean value $m_{\bm{e}_{i}}$, and then substitute the result into $m_{\bm{e}_{i}}^{2}$ of the equation (25) to obtain a linear objective function. The estimated lower bound of variance monotonically increases with the truncation order $\mu$ since the upper bound of the mean value $m_{\bm{e}_{i}}$ monotonically decreases.

As the optimization problems shown in this subsection use convex relaxation, there may be the case where the truncation order $\mu$ needs to be unrealistically large to reduce the conservativeness of the bounds to a practically useful level. To the authors’ experience, however, one will only need to use $\mu=6$ to 12-th order moments to obtain a practically useful bounds for small reaction networks with $n=3$ to 5 molecules. The general quantitative analysis, however, is still an open question.

Remark 2. The moment closure approaches [15, 16] allow for approximating the high order moments, $m_{\bm{\beta}}$ with ${\rm deg}(m_{\bm{\beta}})=\mu+1$, using the low order moments, $m_{\bm{\beta}}$ with ${\rm deg}(m_{\bm{\beta}})\leq\mu$, to solve the equation (10). This approach, however, is based on assumptions on the underlying distribution $P_{\bm{x}}(t)$. Thus, the quantification of the approximation error is not easy. On the other hand, the proposed semidefinite programs are advantageous in that they can compute the upper and the lower bounds of the statistics of stochastic chemical reactions with theoretical rigor without assuming an underlying distribution.

4 Application to a Stochastic Genetic Circuit

To illustrate the proposed semidefinite programming approach, we analyze a protein dimerization process, which is one of the simplest examples of stochastic biomolecular reactions. Consider the chemical reactions shown in Fig. 2, where there are three species of molecules, DNA, monomer protein and dimer protein, and $r=3$ reactions. In Fig. 2, the protein monomer is first produced from DNA, and then it is either degraded or dimerized. As a result, the copy number of protein monomer reaches a steady state at $t\rightarrow\infty$.

In what follows, our goal is to analyze the steady state mean and variance of the copy number of the protein monomer. Let the copy number of monomer protein and its $\alpha$-th order moment be defined by $x$ and

[TABLE]

respectively. It should be noticed that $m_{0}=1$ by definition. Let $D_{T}$ denote the total copy number of DNA molecules and assume that $D_{T}$ is a constant. The reaction rates of monomer production, degradation and dimerization $w_{i}({\bm{x}})~{}(i=1,2,3)$ are then defined by the polynomials in Table 1. The stoichiometry of the reactions, or the increment of the copy number of monomer protein by each reaction, is given by

[TABLE]

Suppose the values of the parameters are given by Table 2. First, we derive the moment equation (10) by expanding the polynomial (7). Specifically, let the truncation order be $\mu=3$. Then, we have

[TABLE]

It should be noted that the low order moments $m_{1},m_{2}$ and $m_{3}$ depend on the higher order moment $m_{4}$ since the transition rate $w_{3}(x)$ is quadratic in $x$ (see Table 1). The equation (30) implies that the moments of the stochastic reactions are in the nullspace of the matrix. To further constrain the feasibility set, we use the moment condition and the localized moment condition, which are given by

[TABLE]

respectively. In practice, the truncation order $\mu$ needs to be determined so that the gap between the lower and the upper bounds is sufficiently small. This might require computing the optimization problems (22) and (23) multiple times.

To analyze the mean copy number of protein monomer, let the objective function be $f(m_{\bm{\alpha}})=\mathbb{E}[x]=m_{1}$ and solve the optimization problems (22) and (23). Figure 3 (A) illustrates the computed lower and the upper bounds of the mean monomer protein copy number for each truncation order $\mu$. The gap between the bounds monotonically decreases with increasing $\mu$, and with $\mu=5$, we can conclude

[TABLE]

which is a sufficient resolution in practice.

Next, to calculate the upper bound of the variance, let $f(m_{\bm{\alpha}})=\mathbb{E}[(x-\mathbb{E}[x])^{2}]=m_{2}-m_{1}^{2}$, and solve the optimization problem (26). The lower bound is also computed using the relaxation described in Section III. The optimization result in Fig. 3 (B) shows that the variance is bounded by

[TABLE]

It should be noticed that these results are mathematically strict and thus enable rigorous quantification of biocircuit performance using the statistical norms.

Remark 3.  The optimization problem was solved with SeDuMi 1.32 on MATLAB 2016b. To avoid numerical issues, we normalized the variables $m_{i}$ by constants $20^{i}$ in the above numerical example. It should be noted that optimization problems associated with the moment matrices tend to be numerically unstable as reported in [Kuntz2016].

5 Conclusion

In this paper, we have formulated the optimization problem to compute the bounds of the steady state statistics of stochastic chemical reactions in genetic circuits. First, we have introduced the steady state moment equation, which is a set of underdetermined linear equations. To restrict the solution space of the moment equation, our key idea is to use the LMI conditions, which the moments of some measure must satisfy. Consequently, we have obtained a convex relaxation problem that can be solved by semidefinite programming. A distinctive feature of the proposed approach is that it can provide mathematically rigorous upper and lower bounds of the statistics for any stochastic chemical reactions modeled by the CME. To demonstrate the method, we have analyzed the protein dimerization process and have obtained tight bounds of the mean and the variance.

Although we have only considered the steady state moments, we can extend the proposed approach to the analysis of transient moments. The authors are currently working to build a stochastic biocircuit design tool based on the specifications of transient statistics.

Appendix A Proof of Theorem 2.

Proof   The optimization problem (23) is equivalent to the following problem.

[TABLE]

As the objective function is $f(m_{\bm{\alpha}})=m_{2\bm{e}_{i}}-m_{\bm{e}_{i}}^{2}$, the inequality $f(m_{\bm{\alpha}})\geq V$ in the constraint can be equivalently written by

[TABLE]

using Schur compliment [20]. The optimal value of (23) is thus equal to that of (26). $\hfill\Box$