^{1,a)}and Hong Qian

^{1,b)}

### Abstract

A thorough kinetic analysis of the rate theory for stochastic self-regulating gene networks is presented. The chemical master equation kinetic model in terms of a coupled birth–death process is deconstructed into several simpler kinetic modules. We formulate and improve upon the rate theory of self-regulating genes in terms of perturbation theory. We propose a simple five-state scheme as a faithful caricature that elucidates the full kinetics including the “resonance phenomenon” discovered by Walczak *et al.* [Proc. Natl. Acad. Sci. U.S.A.102, 18926 (2005)]. The same analysis can be readily applied to other biochemical networks such as phosphorylation signaling with fluctuating kinase activity. Generalization of the present approach can be included in multiple time-scale numerical computations for large biochemical networks.

The authors are grateful to Jin Wang and Jianhua Xing for helpful discussions and Masaki Sasai and Patrick Warren for valuable comments.

I. INTRODUCTION

II. SELF-REGULATING GENE NETWORK AND WOW RESONANCE

III. ASYMPTOTIC SERIES SOLUTION

A. Nonadiabatic situation

B. Adiabatic situation

C. Escape rate at intermediate region and the condition for WOW resonance

IV. FURTHER INVESTIGATION OF WOW RESONANCE

A. Five-state model

B. Concatenation of five-state blocks

C. Use asymptotic series solution leading terms to fit five-state model parameters

V. DISCUSSION

### Key Topics

- Attractors
- 37.0
- Proteins
- 34.0
- Non adiabatic reactions
- 31.0
- Difference equations
- 11.0
- Perturbation theory
- 10.0

## Figures

Schematic diagrams showing kinetic isomorphism between a self-regulating gene network in (a) and a cellular signaling network based on phosphorylation–dephosphorylation cycle with feedback in (b). (a) Binding of a transcription factor (TF) monomer (χ = 1) or dimer (χ = 2) turns the gene state from 0 to 1; the gene states 0 and 1 determine the rate of TF synthesis with *g* _{0} and *g* _{1}, respectively. The TF is degraded with rate *k*. (b) A phosphorylated enzyme *E** activates a kinase from state *K* to *K*† by binding to the kinase; the *K* and *K*† in turn are the enzymes responsible for the phosphorylation of *E* with rate constants *k* _{0} and *k* _{1}. The dephosphorylation reaction has a rate of *k* _{3}[*P*] where *P* is a phosphatase. The solid lines represent biochemical reactions while the dashed lines represent information flow.

Schematic diagrams showing kinetic isomorphism between a self-regulating gene network in (a) and a cellular signaling network based on phosphorylation–dephosphorylation cycle with feedback in (b). (a) Binding of a transcription factor (TF) monomer (χ = 1) or dimer (χ = 2) turns the gene state from 0 to 1; the gene states 0 and 1 determine the rate of TF synthesis with *g* _{0} and *g* _{1}, respectively. The TF is degraded with rate *k*. (b) A phosphorylated enzyme *E** activates a kinase from state *K* to *K*† by binding to the kinase; the *K* and *K*† in turn are the enzymes responsible for the phosphorylation of *E* with rate constants *k* _{0} and *k* _{1}. The dephosphorylation reaction has a rate of *k* _{3}[*P*] where *P* is a phosphatase. The solid lines represent biochemical reactions while the dashed lines represent information flow.

The self-regulating gene model as a discrete state Markov process. In the theory of the chemical master equation (CME), this type of drawing is called a *chemical master equation graph* (Ref. 20).

The self-regulating gene model as a discrete state Markov process. In the theory of the chemical master equation (CME), this type of drawing is called a *chemical master equation graph* (Ref. 20).

First passage time distribution at nonadiabatic case (σ = 0.01), intermediate adiabaticity (σ = 10), and adiabatic case (σ = 10 000). Left panels: transition time distribution from d state to u state. Right panels: transition time distribution from u state to d state. The parameters used in this and all following computations, if not otherwise specified, are *g* _{ u } = 100, *g* _{ d } = 8, *k* = 1, *h* = 0.0002, *f* = 0.2852.

First passage time distribution at nonadiabatic case (σ = 0.01), intermediate adiabaticity (σ = 10), and adiabatic case (σ = 10 000). Left panels: transition time distribution from d state to u state. Right panels: transition time distribution from u state to d state. The parameters used in this and all following computations, if not otherwise specified, are *g* _{ u } = 100, *g* _{ d } = 8, *k* = 1, *h* = 0.0002, *f* = 0.2852.

Inverse of the real part of first eight nonzero eigenvalues ( from bottom to top) of the coupled birth–death process as a function of adiabaticity parameter σ. All the eigenvalues show a resonancelike behavior as adiabaticity increases, which means the resonance happens in all time scales. The tangling of eigenvalues in intermediate region shows the dynamics in nonadiabatic region and adiabatic region are essentially different.

Inverse of the real part of first eight nonzero eigenvalues ( from bottom to top) of the coupled birth–death process as a function of adiabaticity parameter σ. All the eigenvalues show a resonancelike behavior as adiabaticity increases, which means the resonance happens in all time scales. The tangling of eigenvalues in intermediate region shows the dynamics in nonadiabatic region and adiabatic region are essentially different.

Chemical master equation graph of the truncated self-regulating gene model which the perturbation equations [Eq. (???)] are based on. *n** is the number of proteins at the attractor at u state of the gene. *N* is the maximum allowed number of proteins in the system. In the computation of *T* _{on}, *n** at u line is the tination.

Chemical master equation graph of the truncated self-regulating gene model which the perturbation equations [Eq. (???)] are based on. *n** is the number of proteins at the attractor at u state of the gene. *N* is the maximum allowed number of proteins in the system. In the computation of *T* _{on}, *n** at u line is the tination.

(a) The nonadiabatic and adiabatic asymptotic series solution of the summation of transition rates *r* = *k* _{on} + *k* _{off}. The solid curve is the exact transition rate . The diamond and circle curves are the asymptotic series approximation at non-adiabatic and adiabatic cases, respectively. From top to bottom, the first one, two, and three leading terms are used in the asymptotic series. Note that the solution with one leading term (top graph) is the same as the approximation solution using fast equilibrium. (b) The relative errors of the asymptotic series solution. Error is defined as , where *k* ^{ptb} is the perturbation solution of transition rate, *k* ^{true} is the true transition rate.

(a) The nonadiabatic and adiabatic asymptotic series solution of the summation of transition rates *r* = *k* _{on} + *k* _{off}. The solid curve is the exact transition rate . The diamond and circle curves are the asymptotic series approximation at non-adiabatic and adiabatic cases, respectively. From top to bottom, the first one, two, and three leading terms are used in the asymptotic series. Note that the solution with one leading term (top graph) is the same as the approximation solution using fast equilibrium. (b) The relative errors of the asymptotic series solution. Error is defined as , where *k* ^{ptb} is the perturbation solution of transition rate, *k* ^{true} is the true transition rate.

(a) Kinetic schemes of five-state model. It exhibits the Walczak–Onuchic–Wolynes (WOW) resonance as well as all other features of the rate theory for self-regulating gene networks. (b) A six-state model that represents the full coupled birth–death process with two stable attractors. The transition from either attractor to the other is exactly a five-state model shown in (a).

(a) Kinetic schemes of five-state model. It exhibits the Walczak–Onuchic–Wolynes (WOW) resonance as well as all other features of the rate theory for self-regulating gene networks. (b) A six-state model that represents the full coupled birth–death process with two stable attractors. The transition from either attractor to the other is exactly a five-state model shown in (a).

Phase diagram of resonance region for reduced six-state scheme of self-regulating gene. The left panel is in parameter space of *h* _{3}/*f* _{3} vs *h* _{2}/*f* _{2}, and the right panel is in parameter space of *h* _{1}/*f* _{1} vs *h* _{2}/*f* _{2}. The vertical dash line in the left panel marks the maximum *h* _{2}/*f* _{2} value, which in turn bounds the region in the right panel. The solid lines are conditions of resonance [inequalities (36), (35)], and the dotted-dashed lines are conditions of potential [inequalities (34), (33)]. The other parameters are set as *g* _{1}/*k* _{1} = 3, *g* _{1} = *k* _{2}, *g* _{2} = *k* _{1}.

Phase diagram of resonance region for reduced six-state scheme of self-regulating gene. The left panel is in parameter space of *h* _{3}/*f* _{3} vs *h* _{2}/*f* _{2}, and the right panel is in parameter space of *h* _{1}/*f* _{1} vs *h* _{2}/*f* _{2}. The vertical dash line in the left panel marks the maximum *h* _{2}/*f* _{2} value, which in turn bounds the region in the right panel. The solid lines are conditions of resonance [inequalities (36), (35)], and the dotted-dashed lines are conditions of potential [inequalities (34), (33)]. The other parameters are set as *g* _{1}/*k* _{1} = 3, *g* _{1} = *k* _{2}, *g* _{2} = *k* _{1}.

The self-regulating gene model as a concatenation of a sequence of five-state model blocks. It is an approximation of the full self-regulating gene model, where the TF protein number fluctuation is unidirectional.

The self-regulating gene model as a concatenation of a sequence of five-state model blocks. It is an approximation of the full self-regulating gene model, where the TF protein number fluctuation is unidirectional.

The transition rate of approximate system approximates the true transition rate well. The solid and dash lines are numerical solutions of the transition rates from d state to u state and vise versa, respectively. The circles and diamonds symbols are transition rates between the same attractors in simplified system as shown in Fig. 9.

The transition rate of approximate system approximates the true transition rate well. The solid and dash lines are numerical solutions of the transition rates from d state to u state and vise versa, respectively. The circles and diamonds symbols are transition rates between the same attractors in simplified system as shown in Fig. 9.

Dotted line is exact numerical solution by solving the full MFPT backward equation. Diamonds are approximate solution by fitting the perturbation series to the five-state model, and the asterisks are the result of Walczak *et al.* in nonadiabatic case (Ref. 1).

Dotted line is exact numerical solution by solving the full MFPT backward equation. Diamonds are approximate solution by fitting the perturbation series to the five-state model, and the asterisks are the result of Walczak *et al.* in nonadiabatic case (Ref. 1).

Article metrics loading...

Full text loading...

Commenting has been disabled for this content