Stochastic model reduction using a modified Hill-type kinetic rate law

Abstract

In the present work, we address a major challenge facing the modeling of biochemical reaction networks: when using stochastic simulations, the computational load and number of unknown parameters may dramatically increase with system size and complexity. A proposed solution to this challenge is the reduction of models by utilizing nonlinear reaction rate laws in place of a complex multi-reaction mechanism. This type of model reduction in stochastic systems often fails when applied outside of the context in which it was initially conceived. We hypothesize that the use of nonlinear rate laws fails because a single reaction is inherently Poisson distributed and cannot match higher order statistics. In this study we explore the use of Hill-type rate laws as an approximation for gene regulation, specifically transcription repression. We matched output data for several simple gene networks to determine Hill-type parameters. We show that the models exhibit inaccuracies when placed into a simple feedback repression model. By adding an additional abstract reaction to the models we account for second-order statistics. This split Hill rate law matches higher order statistics and demonstrates that the new model is able to more accurately describe the mean protein output. Finally, the modified Hill model is shown to be modular and models retain accuracy when placed into a larger multi-gene network. The work as presented may be used in gene regulatory or cell-signaling networks, where multiple binding events can be captured by Hill kinetics. The added benefit of the proposed split-Hill kinetics is the improved accuracy in modeling stochastic effects. We demonstrate these benefits with a few specific reaction network examples.

Figure 1

The Hill parameter fits for models 1M, 1D, 2M, and 2D. The parameters K_M and n are given in the plot along with R² values for quantitative comparison. The steady-state of the product from elementary stochastic simulations (blue dots), and the linear fit from Eq. 3 (red dashed line).In the figure, X refers to X_a (monomer).

Figure 2

Models 1M, 1D, 2M, and 2D, mean product count versus time. Elementary (blue solid line) and reduced (red circle line) models. Each line represents an average over 10 000 trajectories. Elementary simulations used Hy3S, reduced simulation used a basic Gillespie next-reaction SSA. Relative root mean squared differences (rRMS) values are provided for comparison.

Figure 3

Models 1M, 1D, 2M, and 2D, mean product count versus time. Elementary (blue solid line) and split-Hill (green circle line) models. For the split-Hill model k⁺ = 0.0086 s⁻¹ (1M), 0.0075 s⁻¹ (1D), 0.000185 s⁻¹ (2M), and 0.00028 s⁻¹ (2D). Each line represents an average over 10 000 trajectories. Relative root mean squared differences (rRMS) values are provided for comparison.

Figure 4

Elementary (blue solid line), Hill (red circle-line), and split-Hill (green dashed line) models, mean product count versus time for gene A. The inset figure describes the reaction network. Each line represents an average over 10 000 trajectories.