Inferring gene regulatory networks from single-cell data: a mechanistic approach

TL;DR

This paper introduces a mechanistic, stochastic modeling approach for inferring gene regulatory networks from single-cell transcriptomics data, accounting for bursty mRNA synthesis and providing a tractable inference method.

Contribution

It develops a modular mechanistic model based on piecewise-deterministic Markov processes and derives a tractable hidden Markov model approximation for network inference.

Findings

The model accurately fits real single-cell data.

Successfully recovers known gene networks in simulations.

Demonstrates potential for analyzing large-scale single-cell data.

Abstract

The recent development of single-cell transcriptomics has enabled gene expression to be measured in individual cells instead of being population-averaged. Despite this considerable precision improvement, inferring regulatory networks remains challenging because stochasticity now proves to play a fundamental role in gene expression. In particular, mRNA synthesis is now acknowledged to occur in a highly bursty manner. We propose to view the inference problem as a fitting procedure for a mechanistic gene network model that is inherently stochastic and takes not only protein, but also mRNA levels into account. We first explain how to build and simulate this network model based upon the coupling of genes that are described as piecewise-deterministic Markov processes. Our model is modular and can be used to implement various biochemical hypotheses including causal interactions between genes. However, a naive fitting procedure would be intractable. By performing a relevant approximation of the stationary distribution, we derive a tractable procedure that corresponds to a statistical hidden Markov model with interpretable parameters. This approximation turns out to be extremely close to the theoretical distribution in the case of a simple toggle-switch, and we show that it can indeed fit real single-cell data. As a first step toward inference, our approach was applied to a number of simple two-gene networks simulated in silico from the mechanistic model and satisfactorily recovered the original networks. Our results demonstrate that functional interactions between genes can be inferred from the distribution of a mechanistic, dynamical stochastic model that is able to describe gene expression in individual cells. This approach seems promising in relation to the current explosion of single-cell expression data.

Figures

Tables

Table 1. Table 1: Chemical reactions defining the two-state model. The rate constants are usually abbreviated to rates as they correspond to actual reactions rates when only one molecule of reactant is present. In the stochastic setting, these rates are in fact propensities, i.e. probabilities per unit of time.

Reaction	Rate constant	Interpretation
$G\to G^{\ast }$	$k_{\text{on}}$	gene activation
$G^{\ast }\to G$	$k_{\text{off}}$	gene inactivation
$G^{\ast }\to G^{\ast }+M$	$s_0$	transcription
$M\to M+P$	$s_1$	translation
$M\to \varnothing$	$d_0$	mRNA degradation
$P\to \varnothing$	$d_1$	protein degradation

Table 2. Table 2: Successive dynamical models introduced in this article, recalling for each step the main feature and the form of the mRNA stationary distribution. The full network model (step 3) is used for simulations, while the simplified one (step 4) is used to derive the approximate statistical likelihood.

1	Single-gene, discrete [29]
	$\diamond$ All molecules are discrete
	$\diamond$ mRNA distribution: Beta-Poisson	$\downarrow$	Abundant species treated continuously
2	Single-gene, PDMP (1)
	$\diamond$ Only the promoter is discrete
	$\diamond$ mRNA distribution: Beta	$\downarrow$	Introduction of interactions via $k_{\text{on}}$, $k_{\text{off}}$
3	Network (2), normalized version (4)
	$\diamond$ Both accurate and fast to simulate
	$\diamond$ mRNA distribution: unknown	$\downarrow$	Timescale separation of Protein/mRNA ($d_0\gg d_1$)
4	Simplified network (7)
	$\diamond$ mRNA is removed from the network
	$\diamond$ Conditional mRNA distribution: Beta (8)

Table 3. Table 3: General parameters used in the examples. The s i , j subscript 𝑠 𝑖 𝑗 s_{i,j} correspond to the normalized model: counterparts in absolute protein numbers are s ¯ i , j = s i , j × ( s 0 ​ s 1 ) / ( d 0 ​ d 1 ) = 2 × 10 3 subscript ¯ 𝑠 𝑖 𝑗 subscript 𝑠 𝑖 𝑗 subscript 𝑠 0 subscript 𝑠 1 subscript 𝑑 0 subscript 𝑑 1 2E3 \overline{s}_{i,j}=s_{i,j}\times(s_{0}s_{1})/(d_{0}d_{1})=$2\text{\times}{10}^{3}$ for i ≠ j 𝑖 𝑗 i\neq j and s ¯ i , i = s i , i × ( s 0 ​ s 1 ) / ( d 0 ​ d 1 ) = 1.9 × 10 4 subscript ¯ 𝑠 𝑖 𝑖 subscript 𝑠 𝑖 𝑖 subscript 𝑠 0 subscript 𝑠 1 subscript 𝑑 0 subscript 𝑑 1 1.9E4 \overline{s}_{i,i}=s_{i,i}\times(s_{0}s_{1})/(d_{0}d_{1})=$1.9\text{\times}{10}^{4}$ .

Parameter	Value	Units
$s_0$	${10}^3$	$\mathrm{mRNA}\cdot h^{-1}$
$s_1$	$10$	$\mathrm{protein}\cdot h^{-1}\cdot {\mathrm{mRNA}}^{-1}$
$d_0$	$0.5$	$h^{-1}$
$d_1$	$0.1$	$h^{-1}$
$k_0$	$0.34$	$h^{-1}$
$k_1$	$2.15$	$h^{-1}$
$k_{\text{off}}$	$10$	$h^{-1}$
$m_{i,j}$	$2$ for $i\ne j$	–
$m_{i,i}$	$2$ (Fig 5, 6) or $3$ (Fig 8)	–
$s_{i,j}$	$0.01$ for $i\ne j$	proteins (normalized)
$s_{i,i}$	$0.095$ from eq. (31)	proteins (normalized)

Table 4. Table 4: Network parameters used in the examples.

Fig 5, 6	${\theta }_{1,1}$	${\theta }_{1,2}$	${\theta }_{2,1}$	${\theta }_{2,2}$
	$4$	$-8$	$-8$	$4$
Fig 8	${\theta }_{1,1}$	${\theta }_{1,2}$	${\theta }_{2,1}$	${\theta }_{2,2}$
1	$0$	$0$	$0$	$0$
2	$0$	$0$	$1$	$0$
3	$0$	$1$	$0$	$0$
4	$-0.1$	$1$	$1$	$-0.1$
5	$0$	$0$	$-1$	$0$
6	$0$	$-1$	$0$	$0$
7	$0$	$-1$	$-1$	$0$

Table 5. Table 5: Inferred network parameters used to generate Fig 8b. Each row refers to one of the ten datasets generated for testing. Colors indicate whether the parameters represent the correct topology (blue) or another one (orange) regarding the true networks (cf. Table  4 ).

Network 1
	${\theta }_{1,2}$	${\theta }_{2,1}$
1	$0$	$0$
2	$0$	$0$
3	$0$	$0$
4	$0$	$0$
5	$0.10$	$0$
6	$0$	$0$
7	$0$	$0$
8	$0$	$0$
9	$0$	$0$
10	$0$	$0$

Table 6. Table 6: Example of inferred networks in the presence of dropouts (30% of the whole dataset) generated by applying a Poisson noise and then a threshold to the “perfect” data. Such zeros where replaced using the procedure described in section  E.1 before inferring the networks.

Network 1
	${\theta }_{1,2}$	${\theta }_{2,1}$
1	$0$	$0$
2	$0$	$0.10$
3	$0$	$0$
4	$0$	$0$
5	$0.10$	$0$
6	$0$	$0$
7	$0.09$	$0.16$
8	$0$	$0$
9	$0$	$0$
10	$0$	$0$

Equations

$$\{\begin{array}{rl}E(t)& :0\stackrel{k_{\text{on}}}{\to }1,1\stackrel{k_{\text{off}}}{\to }0\\ M^{\prime }(t)& =s_{0}E(t)-d_{0}M(t)\\ P^{\prime }(t)& =s_{1}M(t)-d_{1}P(t)\end{array}$$

$$\{\begin{array}{rl}{E}_i(t)& :0\stackrel{k_{\text{on},i}}{\to }1,1\stackrel{k_{\text{off},i}}{\to }0\\ {M_i}^{\prime }(t)& =s_{0,i}{E}_i(t)-d_{0,i}{M}_i(t)\\ {P_i}^{\prime }(t)& =s_{1,i}{M}_i(t)-d_{1,i}{P}_i(t)\end{array}$$

$$G_i\stackrel{{\theta }_{i,0}}{\to }{G_i}^{\ast }\text{and}{G}_i+P_j\stackrel{{\theta }_{i,j}}{\to }{G_i}^{\ast }+P_j$$

$$\{\begin{array}{rl}{E}_i(t)& :0\stackrel{k_{\text{on},i}}{\to }1,1\stackrel{k_{\text{off},i}}{\to }0\\ {M_i}^{\prime }(t)& =d_{0,i}\left(E_i(t)-M_i(t)\right)\\ {P_i}^{\prime }(t)& =d_{1,i}\left(M_i(t)-P_i(t)\right)\end{array}$$

$${P_i}^{\prime }(t)=d_{1,i}\left(\frac{k_{\text{on},i}(\mathbf{P}(t))}{k_{\text{on},i}(\mathbf{P}(t))+k_{\text{off},i}(\mathbf{P}(t))}-P_i(t)\right)$$

$$\{\begin{array}{rl}{E_i}^{t+\delta t}& \sim \mathcal{B}\left({\pi }_i^t\right)\\ {M_i}^{t+\delta t}& =(1-d_{0,i}\delta t){M_i}^t+d_{0,i}\delta t{E_i}^t\\ {P_i}^{t+\delta t}& =(1-d_{1,i}\delta t){P_i}^t+d_{1,i}\delta t{M_i}^t\end{array}$$

$${\pi }_i^t=\frac{a_i^t}{a_i^t+b_i^t}+\left({E_i}^t-\frac{a_i^t}{a_i^t+b_i^t}\right)e^{-(a_i^t+b_i^t)\delta t}$$

$$\{\begin{array}{rl}{E}_i(t)& :0\stackrel{k_{\text{on},i}}{\to }1,1\stackrel{k_{\text{off},i}}{\to }0\\ {P_i}^{\prime }(t)& =d_{1,i}\left(E_i(t)-P_i(t)\right)\end{array}$$

$$M_i|\mathbf{P}\sim \mathrm{Beta}\left(\frac{k_{\text{on},i}(\mathbf{P})}{d_{0,i}},\frac{k_{\text{off},i}(\mathbf{P})}{d_{0,i}}\right)$$

$$u(y)=\prod _{i=1}^n\frac{{y_i}^{a_i(y)-1}{(1-y_i)}^{b_i(y)-1}}{B(a_i(y),b_i(y))}$$

$${\Phi }_i(y)=\exp ({\theta }_{i,i})\prod _{j\ne i}\frac{1+\exp ({\theta }_{i,j})(y_j/s_{i,j}{)}^{m_{i,j}}}{1+(y_j/s_{i,j}{)}^{m_{i,j}}}$$

$$k_{\text{on},i}(y)=\frac{k_{0,i}+k_{1,i}{\Phi }_i(y)(y_i/s_{i,i}{)}^{m_{i,i}}}{1+{\Phi }_i(y)(y_i/s_{i,i}{)}^{m_{i,i}}}.$$

$$u(y)=\prod _{i=1}^n\left(\sum _{r=0}^{c_i}{p}_{i,r}(y)\frac{{y_i}^{a_{i,r}-1}{(1-y_i)}^{b_i-1}}{B(a_{i,r},b_i)}\right)$$

$$p_{i,r}(y)=\frac{\left(\genfrac{}{}{0ex}{}{c_i}{r}\right)B(a_{i,r},b_i)({\Phi }_i(y)/s_{i,i}^{m_{i,i}}{)}^r}{{\sum }_{r^{\prime }=0}^{c_i}\left(\genfrac{}{}{0ex}{}{c_i}{r^{\prime }}\right)B(a_{i,r^{\prime }},b_i)({\Phi }_i(y)/s_{i,i}^{m_{i,i}}{)}^{r^{\prime }}}.$$

$$v(x,y)=\prod _{i=1}^n\frac{{x_i}^{{\tilde{a}}_i(y)-1}{(1-x_i)}^{{\tilde{b}}_i(y)-1}}{B({\tilde{a}}_i(y),{\tilde{b}}_i(y))}$$

$$\ell$$

$$x^{a-1}(1-x{)}^{b-1}$$

$$\approx x^{a-1}\exp (-bx)$$

$$\{\begin{array}{rl}{E}_i(t)& :0\stackrel{k_{\text{on},i}}{\to }1,1\stackrel{k_{\text{off},i}}{\to }0\\ {M_i}^{\prime }(t)& =s_{0,i}{E}_i(t)-d_{0,i}{M}_i(t)\\ {P_i}^{\prime }(t)& =s_{1,i}{M}_i(t)-d_{1,i}{P}_i(t)\end{array}$$

$$(E_1^0,\dots ,E_n^0)\in \{0,1{\}}^n,(M_1^0,\dots ,M_n^0)\in \prod _{i=1}^n\left[0,\frac{s_{0,i}}{d_{0,i}}\right],(P_1^0,\dots ,P_n^0)\in \prod _{i=1}^n\left[0,\frac{s_{0,i}{s}_{1,i}}{d_{0,i}{d}_{1,i}}\right],$$

$${\overline{M}}_i=\frac{d_{0,i}}{s_{0,i}}{M}_i\in [0,1]\text{and}{\overline{P}}_i=\frac{d_{0,i}{d}_{1,i}}{s_{0,i}{s}_{1,i}}{P}_i\in [0,1].$$

$${{\overline{M}}_i}^{\prime }(t)=\frac{d_{0,i}}{s_{0,i}}{M_i}^{\prime }(t)=d_{0,i}\left(E_i(t)-\frac{d_{0,i}}{s_{0,i}}{M}_i(t)\right)=d_{0,i}\left(E_i(t)-{\overline{M}}_i(t)\right)$$

$${{\overline{P}}_i}^{\prime }(t)=\frac{d_{0,i}{d}_{1,i}}{s_{0,i}{s}_{1,i}}{P_i}^{\prime }(t)=d_{1,i}\left(\frac{d_{0,i}}{s_{0,i}}{M}_i(t)-\frac{d_{0,i}{d}_{1,i}}{s_{0,i}{s}_{1,i}}{P}_i(t)\right)=d_{1,i}\left({\overline{M}}_i(t)-{\overline{P}}_i(t)\right).$$

$$\{\begin{array}{rl}{E}_i(t)& :0\stackrel{\overline{k_{\text{on},i}}}{\to }1,1\stackrel{\overline{k_{\text{off},i}}}{\to }0\\ {{\overline{M}}_i}^{\prime }(t)& =d_{0,i}\left(E_i(t)-{\overline{M}}_i(t)\right)\\ {{\overline{P}}_i}^{\prime }(t)& =d_{1,i}\left({\overline{M}}_i(t)-{\overline{P}}_i(t)\right)\end{array}$$

$$\overline{k_{\text{on},i}}({\overline{P}}_1,\dots ,{\overline{P}}_n)=k_{\text{on},i}\left(\frac{s_{0,1}{s}_{1,1}}{d_{0,1}{d}_{1,1}}{\overline{P}}_1,\dots ,\frac{s_{0,n}{s}_{1,n}}{d_{0,n}{d}_{1,n}}{\overline{P}}_n\right)$$

$$(E_i(t),{\overline{M}}_i(t),{\overline{P}}_i(t){)}_i$$

$${\left(E_i(t),\frac{s_{0,1}}{d_{0,1}}{\overline{M}}_i(t),\frac{s_{0,1}{s}_{1,1}}{d_{0,1}{d}_{1,1}}{\overline{P}}_i(t)\right)}_i.$$

$$\{\begin{array}{rl}{M}^{\prime }& =d_0(E-M)\\ P^{\prime }& =d_1(M-P)\end{array}$$

$$\{\begin{array}{rl}M(t)& =E+(M_0-E)e^{-d_0(t-t_0)}\\ P(t)& =E+(P_0-E)e^{-d_1(t-t_0)}+\frac{d_1}{d_0-d_1}(M_0-E)\left(e^{-d_1(t-t_0)}-e^{-d_0(t-t_0)}\right)\end{array}$$

$$P(t)\approx E+(P_0-E)e^{-d_1(t-t_0)}$$

$$\frac{\partial u}{\partial t}+\sum _{i=1}^n\frac{\partial \left(F_{i}u\right)}{\partial y_i}=\sum _{i=1}^n{K}_{i}u$$

Full text

Inferring gene regulatory networks from single-cell data: a mechanistic approach

Ulysse Herbach1,2,3111ulysse.herbach@ens-lyon.fr, Arnaud Bonnaffoux1,2,4, Thibault Espinasse3, Olivier Gandrillon1,2

1 Univ Lyon, ENS de Lyon, Univ Claude Bernard, CNRS UMR 5239, INSERM U1210, Laboratory of Biology and Modelling of the Cell, 46 allée d’Italie Site Jacques Monod, F-69007 Lyon, France

2 Inria Team Dracula, Inria Center Grenoble Rhône-Alpes, France

3 Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, 43 blvd. du 11 novembre 1918, F-69622 Villeurbanne cedex, France

4 The CoSMo company, 5 passage du Vercors, 69007 Lyon, France

1 Introduction

Inferring regulatory networks from gene expression data is a longstanding question in systems biology [1], with an active community developing many possible solutions. So far, almost all studies have been based on population-averaged data, which historically used to be the only possible way to observe gene expression. Technologies now allow us to measure mRNA levels in individual cells [2, 3, 4], a revolution in terms of precision. However, the network reconstruction task paradoxically remains more challenging than ever.

The main reason is that the variability in gene expression unexpectedly stands at a large distance from a trivial, limited perturbation around the population mean. It is now clear indeed that this variability can have functional significance [5, 6, 7] and should therefore not be ignored when dealing with gene network inference. In particular, as the mean is not sufficient to account for a population of cells, a deterministic model – e.g. ordinary differential equation (ODE) systems, often used in inference [8, 9] – is unlikely to faithfully inform about an underlying gene regulatory network. Whether such a deterministic approach could still be a valid approximation or not is a difficult question that may require some biological insight into the system under consideration [10]. Another key aspect when considering individual cells is that they generally have to be killed for measurements: from a statistical point of view, temporal single-cell data therefore should not be seen as a set of time series, but rather snapshots, i.e. independent samples from a time series of distributions.

On the other hand, single-cell data give the opportunity of moving one step further toward a more accurate physical description of gene expression. Molecular processes of gene expression are overall now well understood, in particular transcription, but precisely how stochasticity emerges is still somewhat of a conundrum. Harnessing variability in single-cell data is expected to allow for the identification of critical parameters and also to provide hints about the basic molecular processes involved [11, 12]. Moreover, the variability arising from perturbations in cell populations is often crucial for network reconstruction to succeed [13, 14] as the deterministic inference problem suffers from intrinsic limitations [15]. From this point of view, the same information is expected to be contained in the variability between cells in single-cell data. Some of the few existing single-cell inference methods follow this path, for example using asynchronous Boolean network models [16] or generating pseudo time series [17, 9]. In this article, we use a mechanistic approach in the sense that every part of our model has an explicit physical interpretation. Importantly, mRNA observations are not used as a proxy for proteins since both are explicitly modeled.

Besides, mechanistic models that are accurate enough to describe gene expression at the single-cell level usually do not consider interactions between genes. For example, the so-called “two-state” (aka random telegraph) model has been successfully used with single-cell RNA-seq data [18], but the joint distribution of a set of genes contains much more information than the marginal kinetics of individual genes: our aim is to exploit this information while keeping the mechanistic point of view.

Namely, we propose to view the inference as a fitting procedure for a mechanistic gene network model. Whereas the goal here is not to achieve global predictability performances (e.g. as in [19]), our framework makes it possible to explicitly implement many biological hypotheses, and to test them by going back and forth between simulations and experiments. The main point of this article is to show that a tractable statistical model for network inference from single-cell data can be derived through successive relevant approximations. Finally, we demonstrate that our approach is capable of extracting enough information out of in silico-simulated noisy single-cell data to correctly infer the structures of various two-gene networks.

2 Methods

In this part, we aim at deriving a tractable statistical model from a mechanistic one. We will use the two-state model for gene expression to build a “network of two-state models” by making the promoter switching rates depend on protein levels. Then, successive relevant simplifications will lead to an explicit approximation of a statistical likelihood.

2.1 A simple mechanistic model for gene regulatory networks

2.1.1 Basic block: stochastic expression of a single gene

Our starting point is the well-known two-state model of gene expression [20, 21, 22, 23], a refinement of the model introduced by [24] from pioneering single-cell experiments [25]. In this model, a gene is described by its promoter which can be either active (on) or inactive (off) – possibly representing a transcription complex being “bound” or “unbound” but it may be more complicated [26] – with mRNA being transcribed only during the active periods. Translation is added in a standard way, each mRNA molecule producing proteins at a constant rate. The resulting model (Fig. 1) can be entirely defined by the set of chemical reactions detailed in Table 1, where chemical species $G$, $G^{*}$, $M$ and $P$ respectively denote the inactive promoter, the active promoter, the amount of mRNA and proteins. The mathematical framework generally assumes stochastic mass-action kinetics [27] for all reactions, since they typically involve few molecules compared to Avogadro’s number. In this fully discrete setting, one can use the master equation to compute stationary distributions: for mRNA the exact distribution is a Beta-Poisson mixture [28], and an approximation is available for proteins when they degrade much more slowly than mRNA [29]. In addition, the time-dependent generating function of mRNA is known in closed form [30] and can be inverted in some cases to obtain the transient distribution [28].

In practice, the formulas involve hypergeometric series that are not straightforward to use in a statistical inference framework. Besides, these series essentially arise from the fact that such a discrete model has to enumerate all potential collisions between molecules (the stochastic mass-action assumption in the master equation). It is therefore natural to consider keeping only the most important source of noise, that is, keeping a molecular representation for rare species but describing abundant species at a higher level where molecular noise averages out to continuous quantities. A quick look at reactions in Table 1 indicates that the only rare species are $G$ and $G^{*}$, with quantities $[G]$ and $[G^{*}]$ being equal to 0 or 1 molecule and satisfying the conservation relation $[G]+[G^{*}]=1$. The other two, $M$ and $P$, are not conserved quantities in the model and reach a much wider range in biological situations [31], meaning that saturation constants $s_{0}/d_{0}$ and $s_{1}/d_{1}$ are much larger than 1 molecule.

Hence, letting $E(t)$, $M(t)$ and $P(t)$ denote the respective quantities of $G^{*}$, $M$ and $P$ at time $t$, we consider a hybrid version of the previous model, where $E$ has the same stochastic dynamics as before, but with $M$ and $P$ now following usual rate equations:

[TABLE]

This system simply switches between two ordinary differential equations, depending on the value of the two-state continuous-time Markov process $E(t)$, making it a Piecewise-Deterministic Markov Process (PDMP) [32]. From a mathematical perspective, model (1) rigorously approximates the original molecular model when $s_{0}/d_{0}$ and $s_{1}/d_{1}$ are large enough [33, 34] and interestingly, it has already been implicitly considered in the biological literature [22, 23]. Note also that the stationary distribution of mRNA is a scaled Beta distribution that is exactly the one of the Beta-Poisson mixture in the discrete model [28]. Similarly to a recent approach for a two-gene toggle switch [35], we will use (1) as a basic building block for gene networks.

When both $k_{\text{on}}\ll k_{\text{off}}$ and $d_{0}\ll k_{\text{off}}$, mRNA is transcribed by bursts, i.e. during short periods which make the mRNA quantity stay far from saturation. Hence, the amount transcribed within each burst is approximately proportional to the burst duration, whose mean is $1/k_{\text{off}}$ by definition: this justifies the quantity $s/k_{\text{off}}$ often being called “burst size” or “burst amplitude”. Furthermore, promoter active periods are much shorter than inactive ones so they can be seen as instantaneous, justifying the name “burst frequency” for the inverse of the mean inactive time $k_{\text{on}}$. We place ourselves in this situation as it often occurs in experiments [22, 23, 36, 37, 38]. Note however that these two notions are not clearly defined when relations $k_{\text{on}}\ll k_{\text{off}}$ and $d_{0}\ll k_{\text{off}}$ do not hold.

2.1.2 Adding interactions between genes: the network model

Now considering a given set of $n$ genes, a natural way of building a network is to assume that each gene $i$ produces specific mRNA $M_{i}$ and protein $P_{i}$, and to define a version of model (1) with its own parameters:

[TABLE]

Still, genes have static parameters and do not interact with each other. To get an actual network, we need to go one step further: reactions $G_{i}\to{G_{i}}^{*}$ and ${G_{i}}^{*}\to G_{i}$ are not assumed to be elementary anymore, but rather represent complex reactions involving proteins so that promoter parameters $k_{\text{on},i}$ and $k_{\text{off},i}$ now depend on proteins (Fig. 2a), and a fortiori on time. Our network model will correspond to the explicit definition, for all gene $i$, of functions $k_{\text{on},i}(P_{1},\dots,P_{n})$ and $k_{\text{off},i}(P_{1},\dots,P_{n})$. These functions shall also depend on network-specific parameters quantifying the interactions, thus making the link between “fitting a chemical model” and “inferring a network”. As a toy example, consider replacing $G_{i}\to{G_{i}}^{*}$ with two parallel elementary reactions

[TABLE]

for which applying the law of mass action directly gives $k_{\text{on},i}(P_{1},\dots,P_{n})=\theta_{i,0}+\theta_{i,j}P_{j}$. In a regulatory network (Fig. 2b), it would correspond to adding a directed edge from gene $j$ to gene $i$, with $\theta_{i,0}$ the basal parameter of gene $i$, and $\theta_{i,j}$ the strength of activation of gene $i$ by protein $P_{j}$. We emphasize that the action of $P_{j}$ on the promoter $G_{i}$ is not necessarily direct. For example, $P_{j}$ can instead indirectly modulate the amount/activity of a transcription factor: we suppose in this article that such hidden reactions are fast enough regarding gene expression dynamics so that protein $P_{j}$ is a relevant proxy for the transcription factor. Moreover, although we assume here that interactions can only happen at the level of $k_{\text{on},i}$ and $k_{\text{off},i}$, mainly for identifiability purposes, it is also possible to make $d_{1,i}$ and $s_{1,i}$ depend on proteins without fundamentally changing the mathematical approach (e.g. see [39, 40]).

In order to simplify notations, we normalize model (2) into a dimensionless equivalent model: we rewrite it in terms of new variables $\overline{M}_{i}=\frac{d_{0,i}}{s_{0,i}}M_{i}$ and $\overline{P}_{i}=\frac{d_{0,i}d_{1,i}}{s_{0,i}s_{1,i}}P_{i}$, which have values between [math] and $1$, and report this scale change in the definition of $k_{\text{on},i}$ and $k_{\text{off},i}$ (see Appendix A for details). In the remainder of this article, the new variables will still be denoted by $M_{i}$ and $P_{i}$ as no confusion arises. The resulting normalized model is:

[TABLE]

still omitting the dependence of $k_{\text{on},i}$ and $k_{\text{off},i}$ on $(P_{1}(t),\dots,P_{n}(t))$ for clarity. This form enlightens the fact that $s_{0,i}$ and $s_{1,i}$ are just scaling constants: given a path $(E_{i},M_{i},P_{i})_{i}$ of system (4), one can go back to the physical path by simply multiplying $M_{i}$ by $(s_{0,i}/d_{0,i})$ and $P_{i}$ by $(s_{0,i}/d_{0,i})\times(s_{1,i}/d_{1,i})$.

Therefore, we get a general network model where each link between two genes is directed and has an explicit biochemical interpretation in terms of molecular interactions. The previous example is very simplistic but one can use virtually any model of chromatin dynamics to derive a form for $k_{\text{on},i}$ and $k_{\text{off},i}$, involving hit-and-run reactions, sequential binding, etc. [41]. Such aspects are still far from being completely understood [42, 43, 44, 45] and this simple network model can hopefully be used to assess biological hypotheses. In the next part, we will introduce a more sophisticated interaction form based on an underlying probabilistic model, which is both “statistics-friendly” and interpretable as a non-equilibrium steady state of chromatin environment [43].

2.1.3 Some known mathematical results

Thanks to some recent theoretical results [40, 46], simple sufficient conditions on $k_{\text{on},i}$ and $k_{\text{off},i}$ ensure that the PDMP network model (4) is actually well-defined and that the overall joint distribution of $(E_{i},M_{i},P_{i})_{i}$ converges as $t\to+\infty$ to a unique stationary distribution, which will be the basis of our statistical approach. Namely, we assume in this article that $k_{\text{on},i}$ and $k_{\text{off},i}$ are continuous functions of $(P_{1},\dots,P_{n})$ and that they are greater than some positive constants. These conditions are satisfied in most interesting cases, including the above toy example (3) when $\theta_{i,0}>0$.

Contrary to creation rates $s_{0,i}$ and $s_{1,i}$, degradation rates $d_{0,i}$ and $d_{1,i}$ play a crucial role in the dynamics of the system. Intuitively, the ratios $(k_{\text{on},i}+k_{\text{off},i})/d_{0,i}$ and $d_{0,i}/d_{1,i}$ respectively control the buffering of promoter noise by mRNA and the buffering of mRNA noise by proteins. A common situation is when promoter and mRNA dynamics are fast compared to proteins, i.e. when $d_{0,i}\gg d_{1,i}$ with $(k_{\text{on},i}+k_{\text{off},i})/d_{0,i}$ fixed. At the limit, the promoter-mRNA noise is fully averaged by proteins and model (4) simplifies into a deterministic system [47]:

[TABLE]

where $\mathbf{P}(t)=(P_{1}(t),\dots,P_{n}(t))$. The diffusion limit, which keeps a residual noise, can also be rigorously derived [48]. Unsurprisingly, one recovers the traditional way of modelling gene regulatory networks with Hill-type interaction functions. Equation (5) is useful to get an insight into the behaviour of the system (4) for given $k_{\text{on},i}$ and $k_{\text{off},i}$, yet it should be used with caution. Indeed, the $d_{0,i}/d_{1,i}$ ratio has been shown to span a high range, averaging out to the value $d_{0,i}/d_{1,i}\approx 5$ in mammalian cells [31], for which taking the limit $d_{0,i}\gg d_{1,i}$ is not obvious. This is consistent with recent single-cell experiments showing a high variability of both mRNA and protein levels between cells [37]. In that sense, the PDMP model is much more robust than its deterministic/diffusion counterpart while keeping a similar level of mathematical complexity, which motivates our approach.

2.1.4 Simulation

We propose a simple algorithm to compute sample paths of our stochastic network model (4). It consists in a hybrid version of a basic ODE solver, making it efficient enough to perform massive simulations on large scale networks involving arbitrary numbers of molecules, which would be intractable with a classic molecule-based model (Fig. 3). The deterministic part of the algorithm is a standard explicit Euler scheme, while the stochastic part is based on the transient promoter distribution for single genes: this can be justified by the fact that during a small enough time interval, proteins remain almost constant so genes behave as if $k_{\text{on},i}$ and $k_{\text{off},i}$ were constant. We therefore use Bernoulli steps, in a similar way of a diffusion being simulated using gaussian steps.

After discretizing time with step $\delta t$, the numerical scheme is as follows. Starting from an initial state $({E_{i}}^{0},{M_{i}}^{0},{P_{i}}^{0})_{i}$, the update of the system from $t$ to $t+\delta t$ is given by:

[TABLE]

where the Bernoulli distribution parameter $\pi_{i}^{t}$ is derived by locally solving the master equation for the promoter [49], i.e.

[TABLE]

with the notation $a_{i}^{t}=k_{\text{on},i}({P_{1}}^{t},\dots,{P_{n}}^{t})$ and $b_{i}^{t}=k_{\text{off},i}({P_{1}}^{t},\dots,{P_{n}}^{t})$. Intuitively, the algorithm is valid when $\delta t\ll 1/\max_{i}\left\{{K_{\text{on},i}},{K_{\text{off},i}},{d_{0,i}},{d_{1,i}}\right\}$ where $K_{\text{on},i}$ and $K_{\text{off},i}$ denote the maximum values of functions $k_{\text{on},i}$ and $k_{\text{off},i}$.

2.2 Deriving a tractable statistical model

We will now adopt a statistical perspective in order to deal with gene network inference, considering a set of observed cells. If they are evolving in the same environment for a long enough time, we can reasonably assume that their mRNA and protein levels follow the stationary distribution of an underlying gene network: this distribution can be used as a statistical likelihood for the cells. Furthermore assuming no cell-cell interactions (which may of course depend on the experimental context), we obtain a standard statistical problem with independent samples. Since the stationary distribution of the stochastic network model (4) is well-defined but a priori not analytically tractable, we will derive an explicit approximation and then reduce our inference problem to a traditional likelihood-based estimation. We will do so in two cases: when there is no self-interaction, and for a specific form of auto-activation.

2.2.1 Separating mRNA and protein timescales

It is for the moment very rare to experimentally obtain the amount of proteins for many genes at the single-cell level. We will therefore assume here that only mRNAs are observed. To deal with this problem, we take the protein timescale as our reference by fixing $d_{1,i}$ and assume that promoter dynamics are faster than proteins, i.e. $(k_{\text{on},i}+k_{\text{off},i})\gg d_{1,i}$ in a biologically relevant way, say $(k_{\text{on},i}+k_{\text{off},i})/d_{1,i}>10$ (thus the deterministic limit (5) does not necessarily hold). Furthermore, in line with several recent experiments [50, 37], we assume that $d_{0,i}$ is sufficiently larger than $d_{1,i}$ so that the correlation between mRNAs and proteins produced by the gene is very small: model (4) then can be reduced by removing mRNA and making proteins directly depend on the promoters (see Appendix A). The result is

[TABLE]

which still admits the deterministic limit (5). Since mRNA dynamics are faster than proteins, one can also assume that, given protein levels $\mathbf{P}=(P_{1},\dots,P_{n})$, each mRNA level $M_{i}$ follows the quasi-steady state distribution

[TABLE]

corresponding to the single-gene model [39, 28] with constant parameters $k_{\text{on},i}(\mathbf{P})$ and $k_{\text{off},i}(\mathbf{P})$. Numerically, this approximation works well even for moderate values of $d_{0,i}$, such as $d_{0,i}=5\times d_{1,i}$ (see the Results section).

Biologically, equations (7) and (8) suggest that correlations between mRNA levels may not directly arise from correlations between promoters states (which in fact are weak because of $(k_{\text{on},i}+k_{\text{off},i})\gg d_{1,i}$), but rather originate from correlations between promoter parameters $k_{\text{on},i}$ and $k_{\text{off},i}$, which themselves depend on the protein joint distribution.

Table 2 sums up the successive modelling steps introduced so far. From now on, we will always assume the form (8) for the mRNA distribution, and thus our model is reduced to equation (7) which only involves proteins.

2.2.2 Hartree approximation

In this section, we present the Hartree approximation principle and provide an explicit formula in the particular case of no self-interaction. The simplified model (7) is still not analytically tractable, but it is now appropriate for employing the self-consistent proteomic field approximation introduced in [51, 52] and successfully applied in [53, 54]. More precisely, we will use its natural PDMP counterpart, which will be referred to as “Hartree approximation” since the main idea is similar to the Hartree approximation in physics [51]. It consists in assuming that genes behave as if they were independent from each other, but submitted to a common “proteomic field” created by all other genes. In other words, we transform the original problem of dimension $2^{n}$ into $n$ independent problems of dimension $2$ that are much easier to solve (see Appendix B for details).

When $k_{\text{on},i}$ and $k_{\text{off},i}$ do not depend on $P_{i}$ (i.e. no self-interaction), this approach results in approximating the protein stationary distribution of model (7) by the function

[TABLE]

where $y=(y_{1},\dots,y_{n})=(P_{1},\dots,P_{n})=\mathbf{P}$, $a_{i}(y)={k_{\text{on},i}(y)}/{d_{1,i}}$, $b_{i}(y)={k_{\text{off},i}(y)}/{d_{1,i}}$ and $\operatorname{B}$ is the standard Beta function. Note that promoter states have been integrated out since they are not required by equation (8).

The function $u$ is a heuristic approximation of a probability density function. It is only valid when interactions are not too strong, that is, when $k_{\text{on},i}$ and $k_{\text{off},i}$ are close enough to constants, and it becomes exact when they are true constants. Besides, it does not integrate to $1$ in general. However, this approximation turns out to be very robust in practice and it has the great advantage to be fully explicit (and significantly simpler than in the non-PDMP case), thus providing a promising base for a statistical model.

When $k_{\text{on},i}$ and $k_{\text{off},i}$ depend on $P_{i}$, one can still explicitly compute the Hartree approximation in many cases: we will give an example in the next section. Alternatively, it is always possible to use formula (9) even with self-interactions, giving a correct approximation when the feedback is not too strong, as for other proteins.

2.2.3 An explicit form for interactions

We now propose an explicit definition of functions $k_{\text{on},i}$ and $k_{\text{off},i}$. Recent work [36, 55, 56] showed that apparent increased transcription actually reflects an increase in burst frequency rather than amplitude. We therefore decided to model only $k_{\text{on},i}$ as an actual function and to keep $k_{\text{off},i}$ constant. In this view, the activation frequency of a gene can be influenced by ambiant proteins, whereas the active periods have a random duration that is dictated only by an intrinsic stability constant of the transcription machinery.

Our approach uses a description of the molecular activity around the promoter in a very similar way as Coulon et al. [42]. Accordingly, we make a quasi-steady state assumption to obtain $k_{\text{on},i}$. This idea based on thermodynamics was also used in the DREAM3 in-Silico Challenge [57] to simulate gene networks. However, only mean transcription rate was described (instead of promoter activity in our work), which is inappropriate to model bursty mRNA dynamics at the single-cell level.

We herein derive $k_{\text{on},i}$ from an underlying stochastic model for chromatin dynamics. We first introduce a set of abstract chromatin states, each state being associated with one of two possible rates of promoter activation, either a low rate $k_{0,i}$ or a high rate $k_{1,i}\gg k_{0,i}$. More specifically, such chromatin states may be envisioned as a coarse-grained description of the chromatin-associated parameters that are critical for transcription of gene $i$. Second, we assume a separation of timescales between the abstract chromatin model and the promoter activity, so that the promoter activation reaction depends only on the quasi-steady state of chromatin. In other words, the effective $k_{\text{on},i}$ is a combination of $k_{0,i}$ and $k_{1,i}$ which integrates all the chromatin states: its value depends on the probability of each state and a fortiori on the transitions between them. We propose a transition scheme which leads to an explicit form for $k_{\text{on},i}$, based on the idea that proteins can alter chromatin by hit-and-run reactions and potentially introduce a memory component. Some proteins thereby tend to stabilize it either in a “permissive” configuration (with rate $k_{1,i}$) or in a “non-permissive” configuration (with rate $k_{0,i}$), providing notions of activation and inhibition. A more precise definition and details of the derivation are provided in Appendix C.

The final form is the following. First, we define a function of every protein but $P_{i}$,

[TABLE]

which may represent the external input of gene $i$. Then, $k_{\text{on},i}$ is defined by

[TABLE]

Hence, when the input $\Phi_{i}(y)$ is fixed, $k_{\text{on},i}$ is a standard Hill function which describes how gene $i$ is self-activating, depending on the Hill coefficient $m_{i,i}$ (Fig. 4). The neutral value is set to $\Phi_{i}(y)=1$, so that for this particular value, $s_{i,i}$ is the usual dissociation constant. Moreover, if $\theta_{i,j}=0$ for all $j\neq i$, then $\Phi_{i}$ becomes the constant function $\Phi_{i}(y)=\exp(\theta_{i,i})$, and thus $\theta_{i,i}$ may be seen as a “basal” parameter, summing up all potential hidden inputs. On the contrary, if some $\theta_{i,j}>0$ (resp. $\theta_{i,j}<0)$, then $\Phi_{i}$ becomes itself an increasing (resp. decreasing) Hill-type function of protein $P_{j}$, where $m_{i,j}$ and $s_{i,j}$ again play their usual roles.

The $n\times n$ matrix $\theta=(\theta_{i,j})$ therefore plays the same role as the interaction matrix in traditional network inference frameworks [8]. For $i\neq j$, $\theta_{i,j}$ quantifies the regulation of gene $i$ by gene $j$ (activation if $\theta_{i,j}>0$, inhibition if $\theta_{i,j}<0$, no influence if $\theta_{i,j}=0$), and the diagonal term $\theta_{i,i}$ aggregates the “basal input” and the “self-activation strength” of gene $i$. Note that self-inhibition could be considered instead, but the choice has to be made before the inference since the self-interaction form is notoriously difficult to identify, especially in the stationary regime. In the remainder of this article, we assume that parameters $k_{0,i}$, $k_{1,i}$, $m_{i,j}$ and $s_{i,j}$ are known and we focus on inferring the matrix $\theta$.

A benefit of the interaction form (10) is to allow for a fully explicit Hartree approximation of the protein distribution (see Appendix C for details). In particular, if $m_{i,i}>0$ and $c_{i}=(k_{1,i}-k_{0,i})/(d_{1,i}m_{i,i})$ is a positive integer for all $i$, the approximation is given by

[TABLE]

with $a_{i,r}=((c_{i}-r)k_{0,i}+rk_{1,i})/(d_{1,i}c_{i})$, $b_{i}=k_{\text{off},i}/d_{1,i}$ and

[TABLE]

In other words, the Hartree approximation (11) is a product of gene-specific distributions which are themselves mixtures of Beta distributions: for gene $i$, the $a_{i,r}$ correspond to “frequency modes” ranging from $k_{0,i}$ to $k_{1,i}$, weighted by the probabilities $p_{i,r}(y)$. It is straightforward to check that inhibitors tend to select the low burst frequencies of their target ($a_{i,r}\approx k_{0,i}$) while activators select the high frequencies ($a_{i,r}\approx k_{1,i}$). If $m_{i,i}=0$ for some $i$, then $k_{\text{on},i}$ does not depend on $P_{i}$ so one just has to replace the $i$-th term in the product (11) with the single Beta form as in equation (9), which is equivalent to taking the limit $c_{i}\to+\infty$. Finally, when $m_{i,i}>0$ but $c_{i}$ is not an integer, using $\lceil c_{i}\rceil$ instead keeps a satisfying accuracy.

2.2.4 The statistical model in practice

Our statistical framework simply consists in combining the timescale separation (8) and the Hartree approximation (11) into a standard hidden Markov model. Indeed, conditionally to the proteins, mRNAs are independent and follow well-defined Beta distributions

[TABLE]

where $x=(x_{1},\dots,x_{n})=(M_{1},\dots,M_{n})=\mathbf{M}$, $\widetilde{a}_{i}(y)={k_{\text{on},i}(y)}/{d_{0,i}}$ and $\widetilde{b}_{i}(y)={k_{\text{off},i}(y)}/{d_{0,i}}$. Then one can use (11) to approximate the joint distribution of proteins. Hence, recalling the unknown interaction matrix $\theta$, the inference problem for $m$ cells with respective levels $(\mathbf{M}_{k},\mathbf{P}_{k})_{1\leqslant k\leqslant m}$ is based on the (approximate) complete log-likelihood:

[TABLE]

where we used conditional factorization and independence of the cells.

The basic statistical inference problem would be to maximize the marginal likelihood of mRNA with respect to $\theta$. Since this likelihood has no simple form, a typical way to perform inference is to use an Expectation-Maximization (EM) algorithm on the complete likelihood (13). However, the algorithm may be slow in practice because of the computation of expectations over proteins. A faster procedure consists in simplifying these expectations using the distribution modes: the resulting algorithm is often called “hard EM” or “classification EM” and is used in the Results section. Moreover, it is possible to encode some potential knowledge or constraints on the network by introducing a prior distribution $w(\theta)$. In this case, from Baye’s rule, one can perform maximum a posteriori (MAP) estimation of $\theta$ by using the same EM algorithm but adding the penalization term $\log(w(\theta))$ to $\ell$ during the Maximization step (see Appendix D.1 and the Results section). Alternatively, a full bayesian approach, i.e. sampling from the posterior distribution of $\theta$ conditionally to $(\mathbf{M}_{1},\dots,\mathbf{M}_{m})$, may also be considered using standard MCMC methods.

Taking advantage of the latent structure of proteins, we can also deal with missing data in a natural way: if the mRNA measurement of gene $i$ is invalid in a cell $k$ owing to technical problems, it is possible to ignore it by removing the $i$-th term in the conditional distribution of mRNAs (12). This only modifies the definition of $v$ for cell $k$ in equation (13), ensuring that all valid data is effectively used for each cell.

3 Results

In this part, we first compare the distribution of the mechanistic model (4) to the mRNA quasi-steady state combined with Hartree approximation for proteins, on a simple toggle-switch example. Then, we show that the single-gene model with auto-activation can fit marginal mRNA distributions from real data better than the constant-$k_{\text{on}}$ model. Finally, we successfully apply the inference procedure to various two-gene networks simulated from the mechanistic model.

3.1 Relevance of the approximate likelihood

Starting from the normalized mechanistic model (4), two approximations were used to derive the final statistical likelihood (13): the quasi-steady state assumption for mRNAs given protein levels, and the Hartree approximation for the joint distribution of proteins. Crucially, this approximate likelihood has to be close enough to the exact one in order to preserve the equivalence between inferring a network and fitting the mechanistic model. To get an idea of the accuracy, we considered a basic two-gene toggle switch defined by $k_{\text{on},i}$ following equation (10) with the interaction matrix given by $\theta_{1,1}=\theta_{2,2}=4$ and $\theta_{1,2}=\theta_{2,1}=-8$ (full parameter list in Appendix F). By computing sample paths (Fig. 5), we estimated the stationary distribution and compared it with our approximation, which appeared to be very satisfying, both for proteins and mRNAs (Fig. 6).

3.2 Fitting marginal mRNA distributions from real data

A particularity of single-cell data is to often exhibit bursty regimes for mRNA (meaning $k_{\text{on}}\ll k_{\text{off}}$ and $d_{0}\ll k_{\text{off}}$) and potentially also for proteins (adding $d_{1}\ll k_{\text{off}}$), which are well fitted by Gamma distributions [37]. At this stage, it is worth mentioning that the Gamma distribution can be seen as a limit case of the Beta distribution. Intuitively, when $b\gg 1$ and $b\gg a$ (typically $a=k_{\text{on}}/d_{0}$ and $b=k_{\text{off}}/d_{0}$), most of the mass of the distribution $\operatorname{Beta}(a,b)$ is located at $x\ll 1$ so we have the first order approximation

[TABLE]

and thus $\operatorname{Beta}(a,b)\approx\gamma(a,b)$. This way, formulas (11) and (12) can be easily transformed into Gamma-based distributions. Parameters $s_{0}$ and $k_{\text{off}}$ then aggregate in $k_{\text{off}}/s_{0}$ because of the scaling property of the Gamma distribution, so only this ratio has to be inferred: from an applied perspective, it simply represents a scale parameter for each gene. This remark leads to a possible preprocessing phase that can be used for estimating the crucial basal parameters of the network, without requiring the knowledge of such scale parameters (see Appendix E).

In addition, our network model is able to generate multiple modes while keeping such bursty regimes (Fig. 5), as noticeable in the stationary distribution (11). Interestingly, this feature has already been considered in the literature by empirically introducing mixture distributions [58, 59]. As a first step toward applications, we compared our model in the simplest case (independent genes with auto-activation) to marginal distributions of single-cell mRNA measurements from [38]. Our model was fitted and compared to the basic two-state model in the bursty regime, i.e. to a simple Gamma distribution: Fig. 7 shows the example of the LDHA gene. Although very close when viewed in raw molecule numbers, the distributions differ after applying the transformation $x\mapsto x^{\alpha}$ with $\alpha=1/3$, which tends to compress great values while preserving small values. The data becomes bimodal, suggesting the presence of two bursting regimes, a “normal” one and a very small “inhibited” one: the auto-activation model then performs better than the simple Gamma, which necessarily stays unimodal for $0<\alpha<1$. Note that the RTqPCR protocol used in [38] was shown to be far more sensitive than single-cell RNA-seq in the detection of low abundance transcripts [60]. Since the data also contains small nonzero values, this tends to support a true biological origin for the peak in zero. Besides, the case of distributions that are not bimodal until transformed also arises for proteins [61].

3.3 Application of the inference procedure

By construction of the mechanistic model, the interaction matrix $\theta$ can describe any oriented graph by explicitly defining causal quantitative links between genes, which is difficult to do within traditional statistical frameworks (e.g. bayesian networks or undirected Markov random fields). The logical downside is that identifiability issues seem inevitable. In a first attempt to assess this aspect, we implemented the inference method presented above and tested it on various two-gene networks, assuming auto-activation for each gene (i.e. $m_{i,i}>0$) with equation (10) to maximize variability without considering perturbations of the system (parameter list in Appendix F).

We decided to investigate the worst case scenario in terms of cell numbers. We are fully aware of the existence of technologies allowing to interrogate thousands of cells simultaneously, but most of the recent studies still rely upon a much smaller number of cells. For each network, we therefore simulated mRNA snapshot data for 100 cells using the full PDMP model (4). We then inferred the matrix $\theta$ using a “hard EM” algorithm based on the likelihood (13), that is, alternatively maximizing the likelihood with respect to $\theta$ and with respect to the (unknown) protein levels of each cell. A lasso-like penalization term, corresponding to a prior distribution, was added to the $\theta_{i,j}$ for $i\neq j$ to obtain true zeros – so that the inferred network topology is clear – and to prevent keeping both $\theta_{i,j}$ and $\theta_{j,i}$ when one is significantly weaker (see Appendix D.1 for details of the penalization and the whole procedure).

We obtained highly encouraging results since every structure was inferred with a high probability of success (Fig. 8), meaning that the non-diagonal (i.e. interaction) terms of $\theta$ had the right sign and were nonzero at the right places. A list of the inferred values is provided in Table 5 of Appendix F. It is very important at that stage to emphasize that we are not trying to infer $\theta$ exactly: we only assess whether it has a zero or nonzero value and its sign. Although the results tend to support the identifiability of the full matrix $\theta$ in this simple two-gene case, one has to be aware that the quantity we maximize (an approximate likelihood) is a priori nonconvex and can have several local maxima (i.e. networks that are relevant candidates to explain the data). The result of the inference thus can depend on the starting point: in this first approach we chose the null matrix to be the starting point for $\theta$, which corresponds to the – biologically relevant – expectation of “balanced” behaviors (e.g. we do not expect $\theta_{1,1}\ll\theta_{2,2}$). Alternatively, one can consider some probabilistic prior knowledge on $\theta$ to implement a (possibly rough) idea of parameter values from a Bayesian viewpoint: it is worth mentioning that any knockout information can be implemented this way in our model.

Finally, we assessed the inference behavior in the presence of dropouts, i.e. genes expressed at a low level in a cell that give rise to zeros after measurement [4]. Our first tests tend to indicate that our approach is robust regarding dropouts, in the sense that up to 30% of simulated dropouts does not drastically affect the estimation of $\theta$ once the other parameters have been estimated correctly (see Table 6 of Appendix F for an example).

4 Discussion

In this paper, we introduce a general stochastic model for gene regulatory networks, which can describe bursty gene expression as observed in individual cells. Instead of using ordinary differential equations, for which cells would structurally all behave the same way, we adopt a more detailed point of view including stochasticity as a fundamental component through the two-state promoter model. This model is but a simplification of the complexity of the real molecular processes [42]. Modifications have been proposed, from the existence of a refractory period [23] to its attenuation by nuclear buffering [62]. In bacteria, the two states originate from the accumulation of positive supercoiling on DNA which stops transcription [63]. In eukaryotes, although its molecular basis is not quite understood, the two-state model is a remarkable compromise between simplicity and the ability to capture real-life data [22, 36, 18, 37, 38]. Our PDMP framework appears to be conceptually very similar to the random dynamical system proposed in [64] but it has two major advantages: time does not have to be discretized, and the mathematical analysis is significantly easier. We also note that a similar framework appears in [65, 66] and that a closely related PDMP – which can be seen as the limit of our model for infinitely short bursts – has recently been described in [67].

We then derive an explicit approximation of the stationary distribution and propose to use it as a statistical likelihood to infer networks from single-cell data. The main ingredient is the separation of three physical timescales – chromatin, promoter/RNA, and proteins – and the core idea is to use the self consistent proteomic field approximation from [51, 52] in a slightly simpler mathematical framework, providing fully explicit formulas that make possible the massive computations usually needed for parameter inference. From this viewpoint, it is a rather simple approach and we hope it can be adapted or improved in more specific contexts, for example in the study of lineage commitment [68]. Besides, the main framework does not necessarily has to include an underlying chromatin model and thus it can in principle also be used to describe gene networks in procaryotes.

Mechanistic modelling and statistical inference

An important quality of the PDMP network model is that the simulation algorithm is comparable in speed with classic ODE and diffusion systems, while providing an effective approximation of the “perfect”, fully discrete, molecular counterpart [33, 35]. It is worth noticing that the PDMP – at least the promoter-mRNA system – naturally appears as an example of Poisson representation [69, 28], that is, not a simple approximation but rather the core component of the exact distribution of the discrete molecular model. Furthermore, such a simulation speed allowed us to compare our approximate likelihood with the true likelihood for a simple two-gene toggle switch, giving excellent results (Fig. 6). This obviously does not constitute a proof of robustness for every network: a proper quantitative (theoretical or numeric) comparison is beyond the scope of this article but would be extremely valuable. Intuitively, it should work for any number of genes, provided that interactions are not too strong.

Besides, some widely used ODE frameworks [57, 8, 17] can be seen as the fast-promoter limit of the PDMP model: this limit may not always hold in practice, especially in the bursty regime. In particular, Fig. 5 highlights the risk of using mRNA levels as a proxy for protein levels. It also explains why ordering single-cell mRNA measurements by pseudo-time may not always be relevant, as found in [38]. In [70], the authors use a hybrid model of gene expression to infer regulatory networks: it is very close to the diffusion limit of our reduced model (7) with the difference that the discrete component, called “promoter” by the authors, would correspond to the “frequency mode” in the present article, as visible for proteins in Fig. 5. From such a perspective, our approach adds a description of bursty mRNA dynamics that allows for fitting single-cell data such as in Fig. 7.

Finally, our method performed well for simple two-gene networks (Fig. 8), showing that part of the causal information remains present in the stationary distribution: this suggests that it is indeed possible to retrieve network structures with a mechanistic interpretation, even from bursty mRNA data.

Perspectives

We focused here on presenting the key ideas behind the general network model and the inference method: the logical next step is to apply it to real data and with a larger number of genes, which is the subject of work in progress in our group. In particular, we propose a functional preprocessing phase, detailed in Appendix E, that only requires the knowledge of the ratio $d_{0,i}/d_{1,i}$ to estimate all the relevant parameters before inferring $\theta$. The ratio between protein and mRNA degradation rates (or half-lives) hence appears to be the minimum required for such a mechanistic approach to be relevant. Depending upon the species, mRNA and protein half-lives values can be found in the literature (see e.g. [31] for human proteins half-lives), or should be estimated from ad hoc experiments.

From a computational point of view, the main challenge is the algorithmic complexity induced by the fact that proteins are not observed and have to be treated as latent variables. There is a priori no possibility of reducing this without loosing too much accuracy, and therefore some finely optimized algorithms may be required to make the method scalable. Furthermore, the identifiability properties of the interaction matrix $\theta$ seem difficult to derive theoretically. In this paper we focused on the stationary distribution for simplicity: importantly, several aspects such as time dependence (computing the Hartree approximation in transitory regime) or perturbations (changing the cell’s medium or performing knockouts [71], which can be naturally embedded in our framework) could greatly improve the practical identifiability.

From a biological point of view, our model does not really describe individual cells but rather a concatenation of trajectories obtained by following cells throughout divisions. Experiments suggest that it should be a relevant approximation, providing one considers mRNA and proteins levels in terms of concentrations instead of molecule numbers [72], which is made possible by the PDMP framework. In this view, the cell cycle results in increasing the apparent degradation rates – because of the increase in cell volume followed by division – and thus plays a crucial role for very stable proteins. However, at such a description level, many aspects of possible compensation mechanisms [73] and chromatin dynamics [74] remain to be elucidated. Regarding the latter aspect, our abstract chromatin states were not modeled from real-life data – chromatin composition for instance – but our approach is relevant in that partitioning into dual-type chromatin states as we did is now known as a pervasive feature of all eukaryotic genomes [75, 76, 77, 78].

Acknowledgments

We thank Geneviève Fourel (LBMC/ENSL) for critical reading of the manuscript and Anne-Laure Fougères (ICJ) for her constant support. We also thank all members of the SBDM and Dracula teams for enlightening discussions, and BioSyL Federation and Ecofect Labex for inspiring scientific events. This work was supported by funding from the Institut Rhônalpin des Systèmes Complexes (IXXI) and from French agency ANR (ICEBERG; ANR-IABI-3096).

Appendix A First simplifications

A.1 Normalizing the PDMP network model

In this section we detail the normalization of our network model. Recall that the original model is defined by

[TABLE]

where $k_{\text{on},i}=k_{\text{on},i}(P_{1},\dots,P_{n})$ and $k_{\text{off},i}=k_{\text{off},i}(P_{1},\dots,P_{n})$. First we observe that, given an initial condition

[TABLE]

the system stays in this set for all $t>0$, and we introduce the dimensionless variables:

[TABLE]

Then, since $s_{0,i}$, $s_{1,i}$, $d_{0,i}$ and $d_{1,i}$ are constants, we get

[TABLE]

and

[TABLE]

As a result, we obtain the normalized model:

[TABLE]

where the rescaled interaction function $\overline{k_{\text{on},i}}$ is defined by

[TABLE]

and $\overline{k_{\text{off},i}}$ is defined analogously. It is straightforward to see that, given a path

[TABLE]

of the normalized model (15), the corresponding path of the original model (14) is

[TABLE]

In this sense, both models are equivalent: in the main text and in the next sections, we always consider model (15) but forget the “bars” to keep the notations simple.

A.2 Separating mRNA and protein timescales

Here we justify the reduced network model involving only promoters and proteins, which is valid when $d_{1,i}\ll d_{0,i}$ for all gene $i$. A full proof is beyond the scope of this article but we provide a heuristic explanation. We temporarily drop the $i$ index for simplicity. Let $t_{1}\geqslant t_{0}\geqslant 0$ and $E\in\{0,1\}$, and suppose $E(t)=E$ for all $t\in[t_{0},t_{1}]$. Moreover, let $M_{0}=M(t_{0})\in[0,1]$ and $P_{0}=P(t_{0})\in[0,1]$. If $d_{1}<d_{0}$, the solution of the linear ODE system

[TABLE]

is given for $t\in[t_{0},t_{1}]$ by

[TABLE]

Hence, if $d_{1}\ll d_{0}$, we have

[TABLE]

using the fact that $|M_{0}-E|\leqslant 1$ and $|e^{-d_{1}(t-t_{0})}-e^{-d_{0}(t-t_{0})}|\leqslant 1$, and thus $P(t)$ approximates the solution of the differential equation ${P}^{\prime}=d_{1}\left({E}-{P}\right)$.

Appendix B Hartree approximation

B.1 Hartree approximation for the PDMP model

Before deriving the approximation, we introduce some notation. Let $n$ be the number of genes in the network, $\mathcal{E}=\{0,1\}^{n}$ and $\Omega={(0,1)}^{n}$. At time $t$, promoter and protein configurations are denoted by $E_{t}=(e_{1},\dots,e_{n})=e\in\mathcal{E}$ and $P_{t}=(y_{1},\dots,y_{n})=y\in\Omega$, respectively. The distribution of $(E_{t},P_{t})$ then evolves along time according to its Kolmogorov forward (aka master) equation, which is a linear partial differential equation (PDE) system in our case. This system is high dimensional ($|\mathcal{E}|=2^{n}$, the number of possible promoter configurations) but the associated linear operator contains lots of zeros. Using the tensor product notation $\otimes$, one can write down the equation in a compact form:

[TABLE]

where $u(t,y)=(u_{e}(t,y))_{e\in\mathcal{E}}\in\mathbb{R}^{2^{n}}\simeq(\mathbb{R}^{2})^{\otimes n}$ represents the probability density function (pdf) of $(E_{t},P_{t})$, and matrices $F_{i}(y_{i})$, $K_{i}(y)\in\mathcal{M}_{2^{n}}(\mathbb{R})\simeq\mathcal{M}_{2}(\mathbb{R})^{\otimes n}$ are defined by

[TABLE]

with

[TABLE]

The sum in the left side of equation (16) clearly corresponds to a deterministic transport term, while the right side corresponds to the stochastic transitions between promoter configurations.

Furthermore, the PDE system comes with the boundary condition

[TABLE]

and the probability condition

[TABLE]

The self-consistent “Hartree” approximation consists in splitting this $2^{n}$-dimensional problem into $n$ independent $2$-dimensional problems by “freezing” the $y_{j}$ for $j\neq i$ where $i$ is fixed, and then gathering the solutions by taking their tensor product to produce an approximation of the true pdf (see [52] for a heuristic explanation in the discrete protein setting). More precisely, one reduced problem is derived for each gene $i$ from (16)-(17)-(18):

[TABLE]

where $u^{i}(t,y)=(u_{0}^{i}(t,y),u_{1}^{i}(t,y))^{\top}\in\mathbb{R}_{+}^{2}$ satisfies the initial condition $u^{i}(0,y)=u^{i,0}(y)$, the boundary condition $F^{(i)}(y_{i})u^{i}(y)\to 0$ when $y_{i}\to 0$ or $1$, and the probability condition $\int_{0}^{1}[u_{0}^{i}(t,y)+u_{1}^{i}(t,y)]\operatorname{d}\!{y_{i}}=1$ for all $t\geqslant 0$ and $y_{1},\dots,y_{i-1},y_{i+1},\dots,y_{n}\in(0,1)$. Therefore, each $u^{i}$ is a pdf with respect to $(e_{i},y_{i})\in\{0,1\}\times(0,1)$ but not on $\mathcal{E}\times\Omega$. Finally, the Hartree approximation is given by

[TABLE]

where the equality holds if for all $i$, $k_{\text{on},i}$ and $k_{\text{off},i}$ only depend on $y_{i}$.

B.2 Solving the reduced problem

For the moment, the time-dependent closed-form solution of (19) is unavailable, but the unique stationary solution can be easily obtained if one knows a primitive of

[TABLE]

which is the nonzero eigenvalue of the matrix $M^{(i)}=K^{(i)}(F^{(i)})^{-1}$. Indeed, letting $v^{i}=F^{(i)}u^{i}$, the stationary equation for $v^{i}$ from (19) becomes

[TABLE]

and then, crucially using the fact that $M^{(i)}$ has a constant eigenvector $(-1,1)^{\top}$ associated with eigenvalue $\lambda_{i}$ (the other eigenvalue being 0), one can check that $v^{i}=e^{\varphi_{i}}(-1,1)^{\top}$ is a solution when $\frac{\partial\varphi_{i}}{\partial y_{i}}=\lambda_{i}$. If one has such a $\varphi_{i}$, the stationary solution of (19) is given by

[TABLE]

where ${Z_{i}}$ is the normalizing constant (which may still depend on $y_{j}$ for $j\neq i$). Note that the existence of a positive constant $\alpha$ such that $\min(k_{\text{on},i},k_{\text{off},i})\geqslant\alpha$ imposes the limit 0 for $\exp(\varphi_{i}(y))$ when $y_{i}\to 0$ or $1$, and thus the boundary condition is satisfied. We also obtain the promoter probabilities $p_{0,i}=p(e_{i}=0)=Z_{0,i}/{Z_{i}}$ and $p_{1,i}=p(e_{i}=1)=Z_{1,i}/Z_{i}$ where $Z_{0,i}=\int_{0}^{1}y_{i}^{-1}\exp(\varphi_{i}(y))\operatorname{d}\!{y_{i}}$, $Z_{1,i}=\int_{0}^{1}(1-y_{i})^{-1}\exp(\varphi_{i}(y))\operatorname{d}\!{y_{i}}$ and $Z_{i}=Z_{0,i}+Z_{1,i}$.

In particular, when $k_{\text{on},i}$ and $k_{\text{off},i}$ do not depend on $y_{i}$ (i.e. no self-interaction), we get

[TABLE]

which gives the classical solution

[TABLE]

with $a_{i}={k_{\text{on},i}(y)}/{d_{1,i}}$ and $b_{i}={k_{\text{off},i}(y)}/{d_{1,i}}$. This form makes clear the promoter probabilities $p_{0,i}$ and $p_{1,i}$ and the conditional distributions of protein $y_{i}$ given the promoter state $e_{i}=0$ or $1$, both being Beta distributions. Since the state is usually not observed, one usually considers the marginal pdf of $y_{i}$, which is also a Beta:

[TABLE]

Note that the conditional distribution of mRNA given proteins also has the form (23) since the PDMP equation is the same, although the argument is not the Hartree approximation but rather the more common quasi-steady state assumption.

B.3 Protein marginal distribution

Given the form of the solution (21), it is in fact always straightforward to integrate over promoters, even for the full (stationary) Hartree approximation (20), and we finally obtain

[TABLE]

where we recalled the possible dependence of $Z_{i}$ on some $y_{j}$. Hence, when $\varphi_{i}$ and $Z_{i}$ are known functions, one gets a fully explicit approximation of the joint protein distribution.

Appendix C Explicit interactions

Here we derive an explicit form for the interactions between genes, starting from a coarse-grained biochemical model. That is, for a given gene $i$, we focus on defining functions $k_{\text{on},i}(y_{1},\dots,y_{n})$ and $k_{\text{off},i}(y_{1},\dots,y_{n})$ where $y_{1},\dots,y_{n}$ denote the protein quantities. For simplicity, we drop the $i$ index in this section when there is no ambiguity.

C.1 Simple biochemical model

The basic idea is to slightly refine the two-state model of gene expression: in addition to the usual switching reactions (whose rates are $k_{\text{on}}$ and $k_{\text{off}}$), we consider a set of reversible transitions between some chromatin states (e.g. describing enhancer regions). Each chromatin state is then associated with a particular rate for the promoter activation reaction. For simplicity, we consider only two cases: a high rate $k_{1}$ (the chromatin will be said permissive) and a low rate $k_{0}\ll k_{1}$ (the chromatin will be said non-permissive). Once active, the promoter can switch off at a rate that is supposed to be independent from chromatin states. Finally, we assume that the chromatin transitions are due to fast interactions with ambiant proteins (binding, hit-and-run, etc.) so that the promoter-switching reactions always see chromatin in its quasi-stationary state. Effective rates $k_{\text{on}}$ and $k_{\text{off}}$ can therefore be obtained by averaging over chromatin states: this way, $k_{\text{off}}$ is still a constant and $k_{\text{on}}$ is now defined by

[TABLE]

where $p_{0}$ (resp. $p_{1}$) is the probability of the chromatin being non-permissive (resp. permissive).

We now define an explicit model for chromatin dynamics and compute its stationary distribution to derive $p_{0}$ and $p_{1}$ as functions of $y_{1},\dots,y_{n}$. We consider $2^{n}$ permissive configurations and $2^{n}$ non-permissive configurations as follows: for all $I\subset\mathcal{G}$ where $\mathcal{G}=\{1,\dots,n\}$, species $C_{I}$ (resp. $C^{*}_{I}$) stands for the chromatin being non-permissive (resp. permissive) and in state $I$. The underlying physics are the following: the chromatin has two “basal” configurations $C_{\emptyset}$ (non-permissive) and $C^{*}_{\emptyset}$ (permissive), which describe dynamics when no protein is present, according to the reactions

[TABLE]

Then, each protein $P_{j}$ is able to modify the chromatin state through a “hit-and-run” reaction, which is kept in memory by encoding the index $j$ in the list $I$, giving the state $C_{I}$ or $C^{*}_{I}$. Eventually, this memory can be lost by “emptying” $I$ step by step (going back to the basal configuration). That is, for all $I\subset\mathcal{G}$ and $j\in\mathcal{G}\setminus I$, we consider the reactions

[TABLE]

[TABLE]

The system then evolves with $[C_{I}],[C_{I}^{*}]\in\{0,1\}$ and $\sum_{I}[C_{I}]+[C_{I}^{*}]=1$, so that only one molecule is present at a time: its species therefore entirely describes the state of the system. Mathematically, we obtain a standard jump Markov process with $2^{n+1}$ states. For example, the case $n=2$ leads to the scheme of Figure 9, writing $\overline{a}_{j}=a_{j}[P_{j}]$ and $\overline{c}_{j}=c_{j}[P_{j}]$ for simplicity. The underlying idea is that, depending on $a_{j}$, $b_{j}$, $c_{j}$ and $d_{j}$, proteins will tend to stabilize the chromatin either in a permissive configuration or in a non-permissive one – providing notions of activation and inhibition. The basal reactions with rates $\alpha$ and $\beta$ sum up what we do not observe (i.e. what is likely to happen for the chromatin when none of the $P_{j}$ are present).

C.2 Stationary distribution

Letting $\mathcal{S}=\{0,1\}^{n+1}$, each state can be coded by a vector $s=(s_{0},s_{1},\dots,s_{n})\in\mathcal{S}$ where $s_{0}=1$ if the chromatin is permissive and 0 otherwise, and for $j\geqslant 1$, $s_{j}=1$ if it has been modified by protein $P_{j}$ and 0 otherwise. If all rates are positive, the system has a unique stationary distribution $\pi$ which can be exactly computed from the master equation. More precisely, the probability $\pi_{s}$ of the chromatin being in state $s\in\mathcal{S}$ is given by

[TABLE]

where $\lambda_{j}=a_{j}/b_{j}$, $\mu_{j}=c_{j}/d_{j}$ and $Z$ is a normalizing constant. Now going back to our initial intention of computing $k_{\text{on}}$, we are only interested in the probability for the chromatin to be permissive,

[TABLE]

and the probability for the chromatin to be non-permissive,

[TABLE]

Observing that each product term only depends on one $s_{j}$, these formulas collapse to

[TABLE]

and the distribution condition $p_{0}+p_{1}=1$ gives $Z=\alpha\prod_{j=1}^{n}(\lambda_{j}[P_{j}]+1)+\beta\prod_{j=1}^{n}(\mu_{j}[P_{j}]+1)$. We finally get

[TABLE]

From this formula, it is straightforward to see that $k_{\text{on}}$ will actually depend on a protein $P_{j}$ only if $\lambda_{j}\neq\mu_{j}$, that is, when reactions involving $P_{j}$ have unbalanced speeds and tend to favor either permissive configurations ($\lambda_{j}>\mu_{j}$) or non-permissive configurations ($\lambda_{j}<\mu_{j}$).

C.3 Higher order interactions

So far we only considered that the $P_{j}$ were interacting as monomers. If they in fact interact after forming dimers or other complexes, and if such complex-forming reactions are even faster than chromatin dynamics, one can take this into account by simply replacing $[P_{j}]$ in equation (25) with a function of $[P_{j}]$ corresponding to the quasi-stationary concentration of the complex. This approximation seems to be relevant to capture the overall dependence of $k_{\text{on}}$ on the proteins, the main point being to use a continuous description (e.g. rate equations) for proteins, which are abundant, while keeping a discrete (stochastic) description for chromatin. We chose to replace $[P_{j}]$ with $[P_{j}]^{m_{j}}$ where $m_{j}>0$, which gives our model a general Hill-type form. Note that $m_{j}=2$ (resp. $m_{j}=3$) may represent a correct approximation for $P_{j}$ interacting as a dimer (resp. a trimer) but in general $m_{j}$ does not necessarily have to be an integer.

C.4 The case of auto-activation

A this stage, it is possible to implement self-interaction for gene $i$ by taking $\lambda_{i}\neq\mu_{i}$ in (25) but this leads to obvious identifiability issues: in stationary state, one cannot really distinguish between auto-activation, auto-inhibition and basal level. To cope with these, we restrict ourselves to auto-activation by setting $c_{i}=d_{i}=0$ and keeping only the relevant chromatin states ($C_{I}^{*}$ for all $I$, and $C_{I}$ for $I$ such that $i\notin I$). The system still has a unique stationary distribution and the formula for $k_{\text{on}}$ corresponds to the case $\mu_{i}=0$ in (25). Then, starting from the fact that auto-activation is only relevant when the basal level is small enough (for a bistable behaviour to be possible), we take the limit $\alpha\ll 1$ while keeping $\alpha\lambda_{i}$ fixed: the formula becomes

[TABLE]

where $m_{i}>0$ if gene $i$ activates itself and $m_{i}=0$ otherwise.

C.5 Parameterization for inference

Parameters of equation (26) are still clearly not identifiable: in order to get a more minimal form, we introduce the following parameterization: $s_{j}={\mu_{j}}^{-1/m_{j}}$, $\theta_{j}=\log(\lambda_{j}/\mu_{j})$ for all $j\neq i$, and $s_{i}=(\beta/\alpha)^{1/m_{i}}$, $\theta_{i}=\log(\lambda_{i})$. After simplifying (26), we obtain

[TABLE]

where

[TABLE]

The new parameters have an intuitive meaning: $s_{j}$ can be seen as a threshold for the influence by protein $j$, and $\theta_{j}$ characterizes this influence via its sign and absolute value ($\theta_{j}=0$ implying that $k_{\text{on}}$ does not depend on protein $j$), with the exception that $s_{i}$ and $\theta_{i}$ aggregate a basal behaviour and an auto-activation strength.

Finally, we recall the notation $y_{j}=[P_{j}]$ and reintroduce the index $i$ of the gene of interest and add it to each parameter. Hence, for every gene $i$, the function $k_{\text{on},i}$ is defined by:

[TABLE]

with

[TABLE]

In our statistical framework, we assume that parameters $k_{0,i}$, $k_{1,i}$, $m_{i,j}$ and $s_{i,j}$ are known and we focus on inferring the matrix $\theta=(\theta_{i,j})\in\mathcal{M}_{n}(\mathbb{R})$, which is similar to the interaction matrix in usual gene network inference methods.

C.6 Explicit distribution for an auto-activation model

Here we derive the stationary distribution for a self-activating gene. For simplicity, we drop the $i$ index. In this model, $k_{\text{off}}$ is constant and we assume that there are some constants $\Phi\geqslant 0$, $m\geqslant 0$, $s>0$ and $k_{1}\gg k_{0}>0$ such that $k_{\text{on}}$ has the form

[TABLE]

so the stationary distribution can directly be used in the Hartree approximation of the network model (27), recalling that $\Phi$ has to be independent of the gene’s own protein but can depend on others. Letting $c=(k_{1}-k_{0})/(md_{1})>0$, we are in the case of the explicit solution (21) with

[TABLE]

so the protein distribution is

[TABLE]

To get a fully explicit result, i.e. to compute $Z$, we shall assume that $c$ is a positive integer. If it is not, one can get a satisfying approximation by taking $c=\lceil(k_{1}-k_{0})/(md_{1})\rceil$. Then, expanding (29) using the binomial theorem, we obtain

[TABLE]

where $a_{r}=((c-r)k_{0}+rk_{1})/(d_{1}c)$ and $b=k_{\text{off}}/d_{1}$, and a probabilistic representation of $\underline{u}$ in terms of a mixture of Beta distributions:

[TABLE]

where $f_{r}(y)={y}^{a_{r}-1}{(1-y)}^{b-1}/\operatorname{B}(a_{r},b)$ and $p_{r}=\binom{c}{r}\operatorname{B}(a_{r},b)(\Phi/s^{m})^{r}/Z$.

The dissociation constant $s$ is clearly redundant with the input $\Phi$. We fix the particular value

[TABLE]

for which the arbitrary neutral case $\Phi=1$ is “symmetric”, i.e. $p_{0}=p_{c}$. Note that $s$ actually only depends on the fundamental parameters $k_{0}$, $k_{1}$, $k_{\text{off}}$ and $d_{1}$ (and not on $c$ nor $m$). Figure 10 shows some examples of the resulting distribution, which can be bimodal or not, depending on the value of $c$ (or equivalently, $m$) when all other parameters are fixed.

Appendix D EM algorithm for network inference

D.1 EM algorithm for MAP estimation

Here we briefly recall the formulation of the Expectation-Maximization (EM) algorithm for maximum a posteriori (MAP) estimation. Consider the probabilistic hierarchical model defined by the distribution of proteins $p(y|\theta)$, the distribution of mRNA given proteins $p(x|y,\theta)$, and a prior distribution $p(\theta)$ on the parameters. Assuming we only observe $x$, we want to infer $\theta$ by MAP estimation, that is, find a mode – hopefully the highest – of the posterior distribution $p(\theta\,|\,x)$, which satisfies by Baye’s rule:

[TABLE]

As $p(\theta\,|\,x)$ has a too complex expression to be efficiently maximized, the EM algorithm rather uses $\ell_{\theta}(x,y)=\log(p(\theta,y\,|\,x))$ by iteratively computing $\theta^{t+1}=\arg\max_{\theta}\{Q(\theta,\theta^{t})\}$ given $\theta^{t}$, where

[TABLE]

A well-known result states that at each step we in fact maximize a lower bound of $p(\theta\,|\,x)$, which is the key point of the algorithm and makes it a particular case of “variational method” (see [79] for example). Now, since $p(x)$ (resp. $p(\theta)$) does not depend on $\theta$ (resp. $y$), it turns out that

[TABLE]

where $g(\theta)=-\log(p(\theta))$ and $\overline{Q}(\theta,\theta^{t})=\int[\log p(y\,|\,\theta)+\log p(x\,|\,y,\theta)]p(y\,|\,x,\theta^{t})\operatorname{d}\!{y}$ is the more standard quantity that appears in the “frequentist” EM algorithm for maximum likelihood estimation. Hence, considering a prior on $\theta$ simply results in adding a penalization term $g(\theta)$ during the M step in the algorithm.

For example, if we assume that $\theta_{i,j}$ for $i\neq j$ are independent and follow Laplace distributions, i.e. $p(\theta)=\prod_{i\neq j}\frac{\lambda}{2}\exp(-\lambda|\theta_{i,j}|)$, then $g(\theta)=\lambda\sum_{i\neq j}|\theta_{i,j}|\;+\;C$ where $C=n(n-1)\log(2/\lambda)$. Since $C$ does not depend on $\theta$, this is equivalent to the standard $L^{1}$ (lasso) penalization, which is well known to enforce the sparsity of the network.

D.2 Custom prior on the interactions

Here we consider a custom prior to deal with oriented interactions. Indeed, for every pair of nodes $\{i,j\}$ there are two possible interactions with respective parameters $\theta_{i,j}$ and $\theta_{j,i}$, but it is likely that only one is actually present in the true network. Hence, we want $\theta_{i,j}$ and $\theta_{j,i}$ to “compete” against each other so that only one is nonzero after MAP estimation, unless there is enough evidence in the data that both interactions are present. To this aim, we define the following prior:

[TABLE]

with $\lambda,\alpha\geqslant 0$. Thus $\alpha$ can be seen as a competition parameter, the case $\alpha=0$ leading to the standard lasso penalization parametrized by $\lambda$.

D.3 The algorithm in practice

As visible in (32), the true EM algorithm involves integration against the distribution $p(y\,|\,x,\theta)$, which does not allow for direct numerical integration because of the dimension ($y\in\mathbb{R}^{n}$). To overcome this problem, a first option is Monte Carlo integration – typically by MCMC – leading to a “stochastic EM” algorithm that is slow but accurate if samples are large enough. A faster option consists in approximating $p(y\,|\,x,\theta)$ by its highest mode, i.e. by the Dirac mass $\delta_{\widehat{y}}$ where $\widehat{y}=\arg\max_{y}\{p(y\,|\,x,\theta)\}$. Then it is worth noticing that since $p(y\,|\,x,\theta)\varpropto p(y\,|\,\theta)p(x\,|\,y,\theta)$, the whole procedure can be seen as performing a coordinate ascent on the function $(\theta,y)\mapsto p(\theta,y\,|\,x)$. We chose this option for the examples: it is sometimes called “hard” or “classification” EM, since a particular case leads to the well-known $k$-means clustering algorithm [80]. Unfortunately, theoretical foundations of the true EM algorithm are lost by the hard EM (we do not maximize a lower bound of $p(\theta\,|\,x)$ anymore), but it often gives satisfying results while requiring much less computational time.

In practice, the procedure is the following. Suppose we observe mRNA levels in $m$ independent cells, and let $\mathbf{x}_{k}\in\mathbb{R}^{n}$ (resp. $\mathbf{y}_{k}\in\mathbb{R}^{n}$) denote the mRNA (resp. protein) levels of cell $k$. In line with sections D.1-D.2 and letting $\mathbf{x}=(\mathbf{x}_{1},\dots,\mathbf{x}_{m})$ and $\mathbf{y}=(\mathbf{y}_{1},\dots,\mathbf{y}_{m})$ for simplicity, we define the objective function

[TABLE]

where the complete log-likelihood $\ell(\mathbf{x},\mathbf{y},\theta)$ and the penalization $g(\theta)$ are given by

[TABLE]

with $u(y,\theta)=p(y|\theta)$ and $v(x,y,\theta)=p(x|y,\theta)$.

The algorithm then simply consists in iterating the following two steps until convergence:

[TABLE]

The “approximate E step” (35) can be performed using a standard gradient method since $u$ and $v$ are smooth functions of $y$. The “penalized M step” (36) is a non-smooth maximization problem since $g$ is non-smooth, but it can be performed using a proximal gradient method detailed in the next section. The form of $\ell(\mathbf{x},\mathbf{y},\theta)$ is such that we just need to compute $\nabla\log u$ and $\nabla\log v$.

The formulas for $u$ and $v$ derived from the normalized model are given in the main text: they can be applied once the data has been normalized, i.e. after dividing each mRNA $i$ level by $s_{0,i}/d_{0,i}$. In the bursty regime, this scale parameter is neither identifiable nor necessary. Indeed, as explained in the main text, the “Beta-like” distributions collapse to “Gamma-like” ones which we provide below, and for which the scale parameter is identifiable.

D.3.1 Likelihood form in the basic case

In the basic case (no self-interaction), we have:

[TABLE]

and

[TABLE]

where $a_{i}=k_{\text{on},i}/d_{1,i}$, $b_{i}=(d_{0,i}/s_{1,i})\times(k_{\text{off},i}/s_{0,i})$, $\widetilde{a}_{i}=k_{\text{on},i}/d_{0,i}$ and $\widetilde{b}_{i}=k_{\text{off},i}/s_{0,i}$.

D.3.2 Likelihood form in the auto-activation case

In the auto-activation case (section C.6), the formula for $\log v$ is the same but $\log u$ is given by

[TABLE]

where $c_{i}=\lceil(k_{1,i}-k_{0,i})/(d_{1,i}m_{i,i})\rceil$, $a_{i,r}=((c_{i}-r)k_{0,i}+rk_{1,i})/(d_{1,i}c_{i})$ and $w_{i,r}=\binom{c_{i}}{r}(\Phi_{i}/s_{i,i}^{m_{i,i}})^{r}$.

D.3.3 Gradients

Explicit computation of the gradients is then straightforward (e.g. with $k_{\text{off},i}$ constant and $k_{\text{on},i}$, $\Phi_{i}$ defined by (27)-(28)) but leads to cumbersome formulas: we implemented them in Scilab and the code is available upon request.

D.4 Proximal gradient method

Here we recall a standard proximal gradient method [81] to solve the M step (36) and provide the proximal operator associated with $g(\theta)$. Note that the method seems to converge in practice, even if $g$ is not convex. It is based on the update

[TABLE]

where $\gamma>0$ is a step size (learing rate) and $\operatorname{prox}_{\gamma}$ is the proximal operator associated with $g(\theta)$, defined on $\Theta\simeq\mathbb{R}^{n^{2}-n}$ by

[TABLE]

In fact, for any $i,j\in\{1,\dots,n\}$ such that $i\neq j$, one can see that $\theta_{i,j}$ and $\theta_{j,i}$ appear in the minimized quantity as independent of all other $\theta$ components. Hence, one just has to compute

[TABLE]

and use it for any $(\tau_{1},\tau_{2})=(\tau_{i,j},\tau_{j,i})\in\mathbb{R}^{2}$ to obtain the corresponding components of $\operatorname{prox}_{\gamma}(\tau)$. Then, letting $\varepsilon=\lambda\gamma$ and assuming $\gamma$ small enough such that $\alpha\varepsilon<1$, we obtain

[TABLE]

with 9 cases for the value of $(h_{1},h_{2})$ depending on $(\tau_{1},\tau_{2})$, given by:

\left\{\begin{array}[]{l}\tau_{1}>\varepsilon\\ \tau_{1}>\varepsilon(1+\alpha(\tau_{2}-\varepsilon))\\ \tau_{2}>\varepsilon(1+\alpha(\tau_{1}-\varepsilon))\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=\tau_{1}-\varepsilon(1+\alpha(\tau_{2}-\varepsilon))\\ h_{2}=\tau_{2}-\varepsilon(1+\alpha(\tau_{1}-\varepsilon))\end{array}\right. 2. 2.

\left\{\begin{array}[]{l}\tau_{1}>\varepsilon\\ |\tau_{2}|\leqslant\varepsilon(1+\alpha(\tau_{1}-\varepsilon))\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=\tau_{1}-\varepsilon\\ h_{2}=0\end{array}\right. 3. 3.

\left\{\begin{array}[]{l}\tau_{1}>\varepsilon\\ \tau_{1}>\varepsilon(1+\alpha(-\tau_{2}-\varepsilon))\\ \tau_{2}<-\varepsilon(1+\alpha(\tau_{1}-\varepsilon))\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=\tau_{1}-\varepsilon(1+\alpha(-\tau_{2}-\varepsilon))\\ h_{2}=\tau_{2}+\varepsilon(1+\alpha(\tau_{1}-\varepsilon))\end{array}\right. 4. 4.

\left\{\begin{array}[]{l}|\tau_{1}|\leqslant\varepsilon(1+\alpha(-\tau_{2}-\varepsilon))\\ \tau_{2}<-\varepsilon\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=0\\ h_{2}=\tau_{2}+\varepsilon\end{array}\right. 5. 5.

\left\{\begin{array}[]{l}\tau_{1}<-\varepsilon\\ \tau_{1}<-\varepsilon(1+\alpha(-\tau_{2}-\varepsilon))\\ \tau_{2}<-\varepsilon(1+\alpha(-\tau_{1}-\varepsilon))\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=\tau_{1}+\varepsilon(1+\alpha(-\tau_{2}-\varepsilon))\\ h_{2}=\tau_{2}+\varepsilon(1+\alpha(-\tau_{1}-\varepsilon))\end{array}\right. 6. 6.

\left\{\begin{array}[]{l}\tau_{1}<-\varepsilon\\ |\tau_{2}|\leqslant\varepsilon(1+\alpha(-\tau_{1}-\varepsilon))\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=\tau_{1}+\varepsilon\\ h_{2}=0\end{array}\right. 7. 7.

\left\{\begin{array}[]{l}\tau_{1}<-\varepsilon\\ \tau_{1}<-\varepsilon(1+\alpha(\tau_{2}-\varepsilon))\\ \tau_{2}>\varepsilon(1+\alpha(-\tau_{1}-\varepsilon))\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=\tau_{1}+\varepsilon(1+\alpha(\tau_{2}-\varepsilon))\\ h_{2}=\tau_{2}-\varepsilon(1+\alpha(-\tau_{1}-\varepsilon))\end{array}\right. 8. 8.

\left\{\begin{array}[]{l}|\tau_{1}|\leqslant\varepsilon(1+\alpha(\tau_{2}-\varepsilon))\\ \tau_{2}>\varepsilon\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=0\\ h_{2}=\tau_{2}-\varepsilon\end{array}\right. 9. 9.

\left\{\begin{array}[]{l}|\tau_{1}|\leqslant\varepsilon\\ |\tau_{2}|\leqslant\varepsilon\end{array}\right.  $\Rightarrow$  \left\{\begin{array}[]{l}h_{1}=0\\ h_{2}=0\end{array}\right.

These 9 cases form a partition of $\mathbb{R}^{2}$ and are represented in Figure 11. One can check that the case $\alpha=0$ collapses to the usual proximal operator associated with lasso penalization.

To obtain the results of Fig 8, we used $\lambda=10$, $\alpha=5$ and $\gamma=10^{-4}$. In a broader context, one may use standard cross-validation to derive appropriate values for $\lambda$ and $\alpha$.

Appendix E Dealing with real data

In this section, we propose a pre-processing phase which would be required in order to apply our network inference method to real data. The first step ensures that the approximate likelihood is well-defined, while the second step consists in estimating the basal parameters appearing in functions $k_{\text{on},i}$, $k_{\text{off},i}$. Please note that inferring real networks from real data is beyond the scope of this paper and will be the subject of future papers.

E.1 Spreading zeros

The likelihood does not accept exact zeros (cf. section D.3). This is not a problem with continuous-type data (for instance, based on fluorescence measurements), but it becomes one when dealing with counts (e.g. RNA-Seq). We propose to replace such zeros with relevant positive values. Recall that the PDMP focuses on the promoter and neglects the local molecular noise at the mRNA level. It is therefore natural to consider that, given a value $M>0$ of mRNA level in the PDMP, the actual number $m$ of molecules in the cell is drawn from the Poisson distribution $\mathcal{P}(M)$. Then, a possible way to replace zeros is to go backwards, i.e. to draw a value $M$ from the PDMP distribution conditioned to $m=0$. Namely, we propose the following procedure to be applied independently for each gene:

Infer a gamma distribution $\gamma(a,b)$ (as a local approximation of the PDMP) from the whole data (possibly at a given time-point) using the standard method of moments; 2. 2.

Replace zeros with independent samples from the distribution $\gamma(a,b+1)$, conditioned to be smaller than the smallest positive value that was measured.

This procedure ensures that zeros are replaced with very small values and that no artificial correlation is introduced. The distribution $\gamma(a,b+1)$ comes from the fact that, if $\mathcal{L}(M)=\gamma(a,b)$ and $\mathcal{L}(m|M)=\mathcal{P}(M)$, then a simple computation gives $\mathcal{L}(M|m=0)=\gamma(a,b+1)$.

E.2 Estimating basal parameters

Here we describe a heuristic method to estimate the model-specific parameters (i.e. everything but the matrix $\theta$) when they cannot be measured through ad hoc experiments, in the case of the auto-activation form (27)-(28). Once again we refer to section D.3. Note that for the mechanistic approach to be relevant, one should know at least the ratio $d_{0,i}/d_{1,i}$, which can be obtained by measuring mRNA and protein half-lives. When even this is unavailable, we propose to use the default value $d_{0,i}/d_{1,i}=5$ (mean value derived from the literature, cf. main text).

The main idea consists in noticing that, when protein $i$ is described by the auto-activation model (27)-(28) (thus following the distribution (30)), mRNA $i$ in quasi-steady state happens to be well described by the same distribution class as (30), with the same $m_{i,i}$ and other parameters being divided by $d_{0,i}/d_{1,i}$. More precisely, we perform the following steps:

Estimate $\widetilde{a}_{0,i}=k_{0,i}/d_{0,i}$, $\widetilde{a}_{1,i}=k_{1,i}/d_{0,i}$, $\widetilde{b}_{i}=k_{\text{off},i}/s_{0,i}$, $\widetilde{c}_{i}$ and $\Phi$ from the likelihood

[TABLE]

This can be done for instance using an EM algorithm for each value of $\widetilde{c}_{i}$ in some range (e.g. $\widetilde{c}_{i}=1,2,\dots,10$), and then choosing the “$\arg\max$” tuple $(\widetilde{a}_{0,i},\widetilde{a}_{1,i},\widetilde{b}_{i},\widetilde{c}_{i},\Phi)$. Afterwards, $\widetilde{a}_{0,i}$, $\widetilde{a}_{1,i}$, $\widetilde{b}_{i}$ and $\widetilde{c}_{i}$ are stored ($\Phi$ only serves this step). 2. 2.

Consistently with the definition of the model, we set $m_{i,i}=(\widetilde{a}_{1,i}-\widetilde{a}_{0,i})/\widetilde{c}_{i}$,

[TABLE]

and we choose $b_{i}=\frac{d_{0,i}}{d_{1,i}}\cdot\widetilde{b}_{i}$. Then we define, as an approximation of (31) in the bursty regime,

[TABLE]

Note that such $b_{i}$ is not the “true” value regarding section D.3.2, as we would need to know $\frac{d_{1,i}}{s_{1,i}}$ to apply the formula $b_{i}=\frac{d_{1,i}}{s_{1,i}}\cdot\frac{d_{0,i}}{d_{1,i}}\cdot\widetilde{b}_{i}$. Fortunately, the network inference does not depend on this scale parameter since the Hill threshold $s_{i,i}$ is proportional to $1/b_{i}$. 3. 3.

Last step consists in extrapolating $m_{i,i}$ and $s_{i,i}$ to the remaining unknown parameters $m_{i,j}$ and $s_{i,j}$ (describing how gene $j$ influences gene $i$). Since the crucial point is their coherence with respect to the range of protein $j$, a relevant choice without additional knowledge is, for all $i\neq j$,

[TABLE]

Appendix F Parameter values

F.1 Models

F.2 Results