Hill functions for stochastic gene regulatory networks from master equations with split nodes and time-scale separation

The deterministic Hill function depends only on the average values of molecule numbers. To account for the fluctuations in the molecule numbers, the argument of the Hill function needs to contain the means, the standard deviations, and the correlations. Here we present a method that allows for stochastic Hill functions to be constructed from the dynamical evolution of stochastic biocircuits with specific topologies. These stochastic Hill functions are presented in a closed analytical form so that they can be easily incorporated in models for large genetic regulatory networks. Using a repressive biocircuit as an example, we show by Monte Carlo simulations that the traditional deterministic Hill function inaccurately predicts time of repression by an order of two magnitudes. However, the stochastic Hill function was able to capture the fluctuations and thus accurately predicted the time of repression.

Figure 1

Flowchart for obtaining stochastic Hill functions.

Figure 2

Splitting the product node.

Figure 3

The factorial moments of the complex, $C$, have a polynomial dependence on the input molecules' (1 and 2) factorial moments: $F_1\left(t\right)$, $F_2\left(t\right)$, $F_{11}\left(t\right)$, $F_{12}\left(t\right)$, and $F_{22}\left(t\right)$. The output, $F_A\left(t\right)$ and $F_{AA}\left(t\right)$, are rational functions of 1 and 2.

Figure 4

(a) Molecules 1 and 2 combine to form the complex $C$, which is the output of the cascade. Each node represents an intermediate complex formation, denoted by 3 and 4, Ref. [23]. Consecutive additions of 2 increase the power of the polynomial that represent the state of $C$

Figure 5

Molecule $C$ is given by the cascade. Its levels stay unchanged while it interacts with molecule $\beta$ to produce molecule $\alpha$. In (a), we see molecule $\beta$ drives molecule $A$. Thus, $C$ acts as a repressor for $A$. If we move $A$ to be driven by $\alpha$, then $C$ works as an activator for $A$. Notice that creation of $A$ is not accompanied by an annihilation of $\alpha$ or $\beta$, as is commonly seen in the Michaelis-Menten process. Adding an additional action line to include the missing annihilation would only create, in the mathematical formulas, sums of the form $k_{+}+b$. Since $b\to \infty$ and $k_{+}$ is finite, the effect of an additional action line is irrelevant. Similar approaches are taken for the deterministic case in Refs. [11, 28, 29, 30].

Figure 6

The majority of the paths for molecule $C$ terminate quickly when the molecule level falls to zero. The deterministic model does not capture this quick degradation. The negative and positive autoregulation loops on $C$ are $k_n=1.0\times {10}^{-2}$ and $k_p=9.9\times {10}^{-3}$, respectively. The transition probabilities in the RTF are $a=1.0$ and $b=0.1$. Initial values are $q_C=1$, $q_{\alpha }=1$, and $q_{\beta }=19$.

Figure 7

The deterministic Hill function ($\lambda =1$) does not match the Monte Carlo simulations, which are closer to the case of $\lambda =0$. The minimum average distance between the LC-method stochastic Hill function and the Monte Carlo simulations corresponds to $\lambda =0.41$.

Figure 8

The IFFL biocircuit. $X$ activates $W$ whereas $Y$ represses it. The activation and repression are describe by Hill functions with three different thresholds, ${\xi }_i,i=1,2,3$. $g\left(t\right)$ is a deterministic input signal. The degradation coefficients for $X$, $Y$, and $W$ are $b_1$, $b_2$, and $b_3$, respectively.

Figure 9

The continuous line represents the dynamics produced by the stochastic Hill function, whereas the traditional Hill function case is plotted as a dotted line. The top lines and the bottom lines depict $F_V+{\sigma }_V$ and $F_V-{\sigma }_V$, respectively, with $V=X$ or $V=Y$.

Figure 10

The parameters for the stochastic Hill functions from (43) are ${\xi }_1=5,{\xi }_2=5,{\xi }_3=20$ for the RTF thresholds, $\mathrm{\Gamma }=5$ for the cascade strength, and $\lambda =0.5$. The specific chosen numerical parameters are not crucial for the overall conclusion that there is a distinction between the use of stochastic versus traditional Hill functions.

Figure 11

The effect of $H_r(F_Y,F_{YY},{\xi }_3)$ on $W$ is more powerful than the effect on $H_a(F_X,F_{XX},{\xi }_2)$. The output of I1-FFL is small for the traditional Hill function because $H_r(F_Y,F_{YY},{\xi }_3)$ is almost zero during the duration of the pulse (dotted line).

Figure 12

$F_{XY}$ together with $H_a$ and $H_r$ generates the production of $W$, (43).