Modeling stochastic dynamics in biochemical systems with feedback using maximum caliber

Abstract

Complex feedback systems are ubiquitous in biology. Modeling such systems with mass action laws or master equations requires information rarely measured directly. Thus rates and reaction topologies are often treated as adjustable parameters. Here we present a general stochastic modeling method for small chemical and biochemical systems with emphasis on feedback systems. The method, Maximum Caliber (MaxCal), is more parsimonious than others in constructing dynamical models requiring fewer model assumptions and parameters to capture the effects of feedback. MaxCal is the dynamical analogue of Maximum Entropy. It uses average rate quantities and correlations obtained from short experimental trajectories to construct dynamical models. We illustrate the method on the bistable genetic toggle switch. To test our method, we generate synthetic data from an underlying stochastic model. MaxCal reliably infers the statistics of the stochastic bistability and other full dynamical distributions of the simulated data, without having to invoke complex reaction schemes. The method should be broadly applicable to other systems.

Fig. 1. The genetic toggle switch

DNA plasmid is shown with promoters p_A and p_B of genes g_A and g_B that, when transcribed, produce proteins A and B, respectively. The gene-promoter complex for species A(B) is denoted α(β) in Eq.(1). In this gene circuit, production of A inhibits B or production of B inhibits A. In Eq.(2), this is shown by having A bind to the gene-promoter complex of B, and vice-versa.

Fig.2. Effect of negative feedback strength on production fluctuations of A

Here we look at how F, production rate fluctuations normalized by its mean squared, varies with N and K using h_α = h_β = −4.605, h_A = h_B = 7.60.

Fig.3. Levels of promoter unbound protein A(red) and levels of unbound protein B(green)

(a) Toggle-switch trajectory generated by Gillespie model. Computer-generated trajectories using the parameters d = 0.005, k = 0.1, f = 2 and r = 0.01. The time trace was broken up into time intervals δt = 0.1. Applying MaxCal to the trajectories from (a), we extracted the Lagrange multiplier K = −0.380 and use h_α = h_β = −4.605 and h_A = h_B = 7.6. (b) MaxCal model with parameters extracted from the simulations. Representative trajectories generated using those values in the MaxCal model.

Fig.4. Correlation strength changes hopping dynamics

Given the same Gillespie model as in Fig. (3a), if we now use a value of: (a) K 50% smaller(−0.19) or (b) 50% larger (−0.57) than those extracted by MaxCal, the predicted model behavior is very different.

Fig.5. MaxCal predicts static and dynamical distribution functions

(a) Distribution of particle numbers, (b) Distribution of dwell times. (Green) Gillespie ‘raw data’. (Red) Maximum Caliber model. Parameters are those used in Fig. (3).

Fig.6. Gillespie and MaxCal time traces for a different parameter regime

(a) Gillespie time trace. Gillespie model, using d = 0.005, k = 0.1, f = 100 and r = 2. (b) MaxCal model time trace. Corresponding MaxCal model with the extracted parameter K = −0.294 using h_α = h_β = −4.605 and h_A = h_B = 7.6.

Fig. 7. Prediction of dynamical and static distrbutions for a different parameter regime

(a) Distribution of particle numbers. (b) Distribution of dwell times. (Green) Gillespie simulations. (Red) MaxCal model. Parameters are those used in Fig. (6).