Steady-state fluctuations of a genetic feedback loop: an exact solution

Genetic feedback loops in cells break detailed balance and involve bimolecular reactions; hence exact solutions revealing the nature of the stochastic fluctuations in these loops are lacking. We here consider the master equation for a gene regulatory feedback loop: a gene produces protein which then binds to the promoter of the same gene and regulates its expression. The protein degrades in its free and bound forms. This network breaks detailed balance and involves a single bimolecular reaction step. We provide an exact solution of the steady-state master equation for arbitrary values of the parameters, and present simplified solutions for a number of special cases. The full parametric dependence of the analytical non-equilibrium steady-state probability distribution is verified by direct numerical solution of the master equations. For the case where the degradation rate of bound and free protein is the same, our solution is at variance with a previous claim of an exact solution (Hornos et al, Phys. Rev. E {\bf 72}, 051907 (2005) and subsequent studies). We show explicitly that this is due to an unphysical formulation of the underlying master equation in those studies.

💡 Research Summary

The paper presents a rigorous analytical treatment of a stochastic gene‑feedback loop in which a gene produces a protein that can bind its own promoter and modulate transcription. The reaction scheme consists of four elementary steps: (i) synthesis of free protein, (ii) bimolecular binding of a free protein to an unbound promoter, (iii) dissociation of the bound complex, and (iv) degradation of protein in both free and bound forms, with potentially distinct rates. Because the binding step is bimolecular, detailed balance is broken and the system resides in a non‑equilibrium steady state (NESS).

The authors formulate the master equation in terms of two discrete variables: the number of free protein molecules (n) and the promoter state (\sigma\in{0,1}). By introducing generating functions (G_\sigma(z)=\sum_{n=0}^{\infty}P_{n,\sigma}z^{n}), the coupled master equations are transformed into a pair of first‑order differential equations. Eliminating one function yields a single second‑order linear differential equation of confluent hypergeometric type. Its general solution is expressed through Kummer’s functions (M(a,b;z)) and (U(a,b;z)), where the parameters (a) and (b) are explicit combinations of the kinetic rates (binding/unbinding, synthesis, and degradation). The steady‑state probability distribution (P_{n,\sigma}) follows from the series expansion of these special functions, providing an exact, closed‑form expression valid for arbitrary rate constants.

Special parameter regimes are examined in detail. When the degradation rates of free and bound protein are equal ((\gamma_f=\gamma_b)), the solution simplifies to a mixture of a Poisson distribution (for the free protein) and a Beta distribution (for the promoter occupancy), reproducing known limiting cases. In the fast‑binding/fast‑unbinding limit, the promoter rapidly equilibrates, and the effective transcription rate becomes proportional to the equilibrium binding probability, yielding a reduced one‑dimensional birth‑death process. The case of non‑degrading bound protein ((\gamma_b=0)) is also solved, showing how the system avoids divergence despite the absence of a direct loss channel for the bound species. For each case, analytical formulas for the mean protein number, variance, coefficient of variation, and tail behavior are derived.

The analytical results are validated by extensive numerical work. The authors solve the master equation directly by matrix diagonalization and perform Gillespie stochastic simulations across a wide range of parameter sets. In every instance, the numerical steady‑state distributions match the analytical predictions to machine precision, confirming the correctness of the derived formulas.

A critical contribution of the paper is the clarification of a long‑standing discrepancy with earlier work by Hornos et al. (Phys. Rev. E 72, 051907, 2005). Hornos and collaborators claimed an exact solution for the same model but omitted the degradation term for the bound protein in their master equation, thereby violating probability conservation. The present authors demonstrate analytically that this omission leads to an unphysical steady state, especially evident when (\gamma_f=\gamma_b). By restoring the missing term, the authors obtain the unique physically admissible solution described above.

In summary, the study delivers the first fully general exact solution for the steady‑state probability distribution of a gene‑feedback loop with a single bimolecular reaction. The results illuminate how non‑equilibrium driving, binding cooperativity, and differential degradation shape intrinsic noise in gene expression. This framework can be directly applied to synthetic biology circuit design, quantitative interpretation of single‑cell data, and the broader development of stochastic thermodynamics for biochemical networks.