Markovian Promoter Models: A Mechanistic Alternative to Hill Functions in Gene Regulatory Networks

TL;DR

This paper introduces Markovian Promoter Models as a mechanistic alternative to phenomenological Hill functions in gene regulatory networks, achieving CME compatibility by coarse-graining promoter dynamics into discrete states while coupling to deterministic or stochastic protein dynamics. A key innovation is parameterizing promoter kinetics from ChEC-seq dwell times, enabling in vivo inference of kinetic rates and facilitating data-driven, scalable modeling for whole-cell contexts. The framework is validated across seven diverse systems (GAL, repressilator, Goodwin, toggle switch, I1-FFL, p53-Mdm2, NF-κB), demonstrating comparable stochastic fidelity to full SSA while delivering 10–100× speedups. These results establish a mechanistic, scalable, and data-parameterizable approach for incorporating promoter-level stochasticity into large-scale models, with significant implications for whole-cell simulations and systems biology insights into ultrasensitivity, pulses, and oscillations.

Abstract

Gene regulatory networks are typically modeled using ordinary differential equations (ODEs) with phenomenological Hill functions to represent transcriptional regulation. While computationally efficient, Hill functions lack mechanistic grounding and cannot capture stochastic promoter dynamics. We present a hybrid Markovian-ODE framework that explicitly models discrete promoter states while maintaining computational tractability. Uniquely, we parameterize this model using fractional dwell times derived from ChEC-seq data, enabling the inference of in vivo kinetic rates from steady-state chromatin profiling. Our approach tracks individual transcription factor binding events as a continuous-time Markov chain, linked to deterministic molecular dynamics. We validate this framework on seven gene regulatory systems spanning basic to advanced complexity: the GAL system, repressilator, Goodwin oscillator, toggle switch, incoherent feed-forward loop, p53-Mdm2 oscillator, and NF-$κ$B pathway. Comparison with stochastic simulation algorithm (SSA) ground truth demonstrates that Markovian promoter models achieve similar accuracy to full stochastic simulations while being 10-100$\times$ faster. Our framework provides a mechanistic foundation for gene regulation modeling and enables investigation of promoter-level stochasticity in complex regulatory networks.

Figures (11)

• Figure 1: Dose-response curves for GAL system across four modeling approaches. (A-E) Steady-state protein levels for GAL1, GAL2, GAL3, GAL80, and reporter gene as a function of external galactose concentration. Blue: ODE model with Hill functions. Red: Markovian hybrid model. (Comparative statistics with SSA and CME are provided in Table \ref{['tab:gal_correlations']}). All models exhibit ultrasensitive dose-response curves with similar EC$_{50}$ values. The Markovian model accurately captures both the mean behavior (matching ODE) and stochastic fluctuations (matching SSA variance). Note the gene-specific responses: GAL1 and GAL2 (4-5 binding sites) show sharper transitions than GAL3 and GAL80 (1 binding site), consistent with cooperative binding.
• Figure 2: Quantitative benchmarking of GAL system dynamics. Time series for Reporter gene following step increase in external galactose (0 → 2.0 mM) at t = 10 min. Black lines: SSA (n=3 replicates). Red dashed line: Markovian Model. The Markovian model accurately captures the induction kinetics and the stochastic variability observed in the full SSA simulation (which uses the original CME-ODE formulation).
• Figure 3: Promoter state dynamics for GAL1 gene. (A) Heatmap showing probability distribution over all 15 promoter states (rows) as a function of time (columns). State (k,m) indicates k Gal4p bound and m Gal80p bound. Color intensity represents probability. (B) Single stochastically sampled promoter state trajectory showing discrete jumps between states. (C) Gal4p dimer (blue) and Gal80p dimer (orange) concentrations driving promoter transitions. Following galactose induction at t=10 min, the promoter rapidly transitions from empty state (0,0) to active states (k>m), then partially to repressed states (k=m) as Gal80p accumulates. The equilibrium distribution at t>100 min shows predominance of partially active states (2-3 Gal4p bound, 0-1 Gal80p bound), consistent with experimental ChIP data venturelli2012synergistic.
• Figure 4: Quantitative benchmarking of the Markovian hybrid model against SSA ground truth. (Left) Repressilator: Stationary distribution of cI protein. (Middle) Goodwin Oscillator: Stationary distribution of mRNA X. (Right) Toggle Switch: Bimodal distribution of Protein U. Gray filled areas represent SSA (Ground Truth) distributions. Colored lines show Markovian model results with varying time steps $\Delta t$. The hybrid model accurately reproduces the stochastic distributions, including the bimodal switching regime of the Toggle switch, provided $\Delta t$ is sufficiently small.
• Figure 5: Promoter state validation for Repressilator. (A) Promoter state agreement between Markovian model and SSA ground truth, quantified by Cohen's kappa statistic. All three genes show $\kappa > 0.75$ (substantial agreement). (B) Confusion matrix showing state-by-state agreement. (C) Switching comparisons: Markovian model (red) vs. SSA (green). The Markovian model accurately captures the switching dynamics, with mean switching rates of 8.3 ± 1.2 switches per period, matching SSA (8.1 ± 1.4 switches per period).