Stochastic Dynamics of Gene Switching and Energy Dissipation for Gene Expression

Stochastic dynamics of gene switching and energy dissipation for gene expression are largely unknown, mainly due to the complexity of non-equilibrium mechanisms. Here, based on an important double-deck loop model, the stochastic mechanisms of gene switching and energy dissipation are studied. First, the probability distributions of steady states are calculated theoretically. It is found that the signal can strengthen the choice of gene switching between the “off” and “on” states. Our analysis of energy consumption illustrates that, compared with the synthesis and degradation of proteins, the process of gene switching costs little energy. Our theoretical analysis reveals some interesting insights into the determination of cell state and energy dissipation for gene expression.

1. Introduction

Signal pathways play vital roles in life by cooperating to control more than one biochemical process while consuming free energy supplied by ATP or high-energy bonds to carry out different vital functions. Based on the core negative feedback control loop shared by various adaption biological systems, Lan et al. show that energy dissipation is indicated to stabilize the adapted state against noise (Lan et al., 2012). Further study explores the present analytic results on the non-equilibrium steady-state (NESS) of the model through mapping to a one-dimensional birth-death process, and the result suggests that the adaptation error can be reduced exponentially as the methylation range increases (Wang et al., 2015). In recent research, the number of phase coherent periods is found to be proportional to the free energy consumed per period (Cao et al., 2015). Increasingly numerous theoretical studies focus on the role of energy in biological information processes and biochemical signal transduction (Lan and Tu, 2013; Endres, 2017).

Biological information processes are complex. In the process of skeletal development, extracellular signals activate RhoA, and control the state of downstream genes mainly through two pathways: RhoA/SRF and RhoA/ROCK (Charrasse et al., 2002; Sordella et al., 2003; Tsai et al., 2013; Matsuoka and Yashiro, 2014). The marvelous phenomenon in this signal cascade is that those two pathways exert the opposite effects, as shown in Figure 1. Hence, such a signal cascade is composed of two competitive pathways and possesses the ability to accurately control the vital bio-processes (Wei et al., 1998, 2000, 2001; Meriane et al., 2000; Beqaj et al., 2002; Castellani et al., 2006; Charrasse et al., 2006). The study of stochastic dynamics and energy dissipation in this biological system represents an interesting topic. The competitive networks may present two kinetic characteristics: oscillation and bistable state. Ouyang et al. show a series of works on the non-equilibrium thermodynamics of oscillations within cells and their results revealed that the critical energy dissipation per period depends on both the frequency and strength of the exchange reaction which gave an optimal design for achieving maximum synchronization with a fixed energy budget (Cao et al., 2015; Fei et al., 2018; Zhang et al., 2020). Gene switch as another kind of competitive networks is a representative bistable state system and will be the focus of our study.

FIGURE 1

Figure 1. Signal pathways from RhoA to the target gene in the development of skeletal muscle. RhoA-GDP is activated by extracellular signals to RhoA-GTP (RhoA^*) while RhoA^* is deactivated back to RhoA-GDP as well. RhoA^* contributes to the development of skeletal muscles through two signal pathways: RhoA/ROCK regulates the target gene myoD negatively which is corresponding to the blue cascade, and RhoA/SRF regulates myoD positively which is corresponding to the red cascade.

In this paper, we propose a double-deck loop (DDL) model to describe signal cascades which are similar to those found in the development of skeletal muscle, as shown in Figure 2. By virtue of non-equilibrium statistical physics theory and stochastic dynamics, exact analytical solutions to the steady-state probability distribution are obtained and energy dissipation in the DDL is derived, which allows for many deeper discussions. Our aims are to reveal (i) the crucial factors that determine the state of gene switching in our model and (ii) the energy dissipation in biochemical reactions. We expect that these theoretical results could help us to understand the general principles of signaling selectivity and energy dissipation in gene regulation networks.

FIGURE 2

Figure 2. The detailed cascades in the double-deck loop (DDL) which is a simplified model based on the biochemical reaction networks in Figure 1. The orange stick is the gene in the “on” state and the blue one is the gene in the “off” state; ROCK is represented as the small green ball; Arrow 1 corresponds to the Rho/ROCK pathway in Figure 1 and Arrow 2 corresponds to the Rho/SRF pathway. The negative regulation from ROCK can push down the state of gene from the “on” state to the “off” state along Arrow 1 with rate k₂. The positive regulation from RhoA^* can also push up the state of gene along Arrow 2 with rate ak₁. The Arrow 3 and 4 represent the basic switching of gene's state with rates k_on and k_off respectively. The synthesis and degradation of ROCK will drive the system move along the horizontal directions (i.e., Arrow 5 and 6) with rates ak₃ and k₄, respectively. Since two parallel loops (i.e., the blue and red loops) can be found in our model, it is so called DDL.

2. The Model

In the process of skeletal development, RhoA is activated by extracellular signals from RhoA-GDP to RhoA-GTP (RhoA^⋆), and RhoA^⋆ is deactivated back to RhoA-GDP as well. RhoA^⋆ plays a key role in process of skeletal development and contributes to the regulation of the expression of muscle-specific genes both through RhoA/SRF and RhoA/ROCK signal pathways. RhoA/ROCK pathway triggers a negative control function on the target gene while RhoA/SRF pathway induces a positive regulation. The relevant biochemical processes are shown in Figure 1. The red signal cascade represents RhoA/SRF pathway and the blue one represents RhoA/Rock pathway.

To derive our mathematical model, some assumptions are put forward. (i) The switching dynamics between GTP and GDP have been discussed thoroughly (Lan and Tu, 2013) which has the same dynamic mechanism with the switching between RhoA^⋆ and RhoA-GDP. This aspect is neglected in our model, since we focus on the selectivity between different gene's state modes and energy dissipation of RhoA/ROCK and RhoA/SRF signal pathways. These two pathways are adopted to regulate the development of skeletal muscle. (ii) It is well-known that a large number of genes are involved in this biological process (Matsui et al., 1996). If all genes are considered, it will complicate the modeling and theoretical analysis. We hypothesize that the development state of skeletal muscle can be represented with the state of gene myoD. Muscle-specific genes begin to be expressed when the state of gene myoD is “on”, otherwise these genes are closed when the state of gene myoD is “off”. (iii) Since the role of ROCK in biological activities is vital (Leung et al., 1995; Aelst and D'Souza-Schorey, 1997; Kaibuchi et al., 1999; Cloutier et al., 2010), the detailed biochemical process of RhoA/SRF pathway is neglected and the regulation of RhoA/SRF pathway is simplified to the direct regulation of RhoA^* on gene myoD as shown in Figure 1. Therefore, we mainly aim to discuss the biochemical reactions including the switch of myoD state and synthesis/degradation of ROCK. Based on the above assumptions, the detailed biochemical equations are as follows:

$\begin{array}{l}{\text{RhoA}}^{*}\stackrel{k_3}{\to }{\text{RhoA}}^{*}+\text{ROCK},\text{ROCK}\stackrel{k_4}{\to }\varnothing \\ G_{\text{off}}+{\text{RhoA}}^{*}\stackrel{k_1}{\to }{G}_{\text{on}}+{\text{RhoA}}^{*},G_{\text{off}}\stackrel{k_{\text{on}}}{\to }{G}_{\text{on}}\\ G_{\text{on}}+\text{ROCK}\stackrel{k_2}{\to }{G}_{\text{off}}+\text{ROCK},G_{\text{on}}\stackrel{k_{\text{off}}}{\to }{G}_{\text{off}}.& (1)\end{array}$

The “on” state of myoD is indicated as G_on in Equation (1), while the “off” one is G_off. k₁ is the transition rate of gene state under the positive control and k₂ is the transition rate of gene under the negative control. The basic switching rates between G_on and G_off are k_off and k_on. k₃ is the synthesis rate of ROCK and k₄ is the degradation rate of it.

As shown in Equation (1), every state of gene can be achieved through two ways: G_on can be achieved both through the promotion of RhoA^* and the basic switching; G_off can be achieved both through the repressive control of ROCK and the basic switching. Considering the synthesis and degradation of ROCK which correspond to the increasing and decreasing of the small green balls in Figure 2, two parallel loops can be found. One is clockwise [i.e., $G_{\text{on}}(m)\stackrel{a{k}_3}{\to }{G}_{\text{on}}(m+1)\stackrel{k_2}{\to }{G}_{\text{off}}(m+1)\stackrel{k_4}{\to }{G}_{\text{off}}(m)\stackrel{a{k}_1}{\to }{G}_{\text{on}}(m)$] which is represented as a blue loop in Figure 2. The other one is anticlockwise [i.e., $G_{\text{on}}(m)\stackrel{k_{\text{off}}}{\to }{G}_{\text{off}}(m)\stackrel{a{k}_3}{\to }{G}_{\text{off}}(m+1)\stackrel{k_{\text{on}}}{\to }{G}_{\text{on}}(m+1)\stackrel{k_4}{\to }{G}_{\text{on}}(m)$] which is represented as red loop in Figure 2. Since these two loops are two parallel decks between G_on and G_off, this theoretical model is called “double-deck loop (DDL)” in this paper. Moreover, m is the number of ROCK, and a is the concentration of RhoA^*. Based on above biochemical equations, the chemical master equations (i.e., CME) can be presented as

$\begin{array}{l}\frac{d{P}_0(m,t)}{dt}=k_{\text{off}}{P}_1(m,t)-k_{\text{on}}{P}_0(m,t)+m{k}_2{P}_1(m,t)\\ -a{k}_1{P}_0(m,t)+a{k}_3{P}_0(m-1,t)-a{k}_3{P}_0(m,t)\\ +(m+1)k_4{P}_0(m+1,t)-m{k}_4{P}_0(m,t);\\ \frac{d{P}_1(m,t)}{dt}=-k_{\text{off}}{P}_1(m,t)+k_{\text{on}}{P}_0(m,t)-m{k}_2{P}_1(m,t)\\ +a{k}_1{P}_0(m,t)+a{k}_3{P}_1(m-1,t)-a{k}_3{P}_1(m,t)\\ +(m+1)k_4{P}_1(m+1,t)-m{k}_4{P}_1(m,t)& (2)\end{array}$

where P₀(m, t) is the probability of the gene “off” state and P₁(m, t) are the probability of the gene “on” states.

3. Results

3.1. The Analytical Solutions of the DDL Model

Regulatory networks generally consist of interactional signal pathways. Different signal pathways may dominate cell fate in different circumstances. Based on the biology processes in the development of skeletal (Wei et al., 1998, 2000, 2001; Meriane et al., 2000; Beqaj et al., 2002; Castellani et al., 2006; Charrasse et al., 2006), the DDL model has both activation and inhibitory signal pathways originating from the same input. Therefore, it is interesting to identify the crucial factors in the selection of different signaling pathways. To address this problem, the probability distributions under different a are derived. In response to a persistent input, the synthesis and degradation of ROCK induce the evolution of ROCK, resulting in a steady state of the number of ROCK m. Using the method of probability-generating functions (Qian, 2007; Huang et al., 2015), the analytical expressions of steady-state probability distributions P₀(m) and P₁(m) are obtained

$\begin{array}{l}{P}_0(m)=\frac{1}{k_{\text{on}}+a{k}_1}\left(\frac{a{k}_2{k}_3}{k_2+1}+k_{\text{off}}\right)A_0\frac{1}{m!}{\left(\frac{a{k}_3}{k_2+1}\right)}^m\\ \times {\displaystyle \sum _{l=0}^m}{k}_2^l{C}_m^l\frac{{\left(\alpha \right)}_l}{{\left(\beta \right)}_l}{}_{1}F_1\left(\alpha +l,\beta +l;{\omega }_2\right)\\ -\frac{1}{k_{\text{on}}+a{k}_1}\frac{a{k}_2{k}_3}{k_2+1}\frac{\alpha }{\beta }{A}_0\frac{1}{m!}{\left(\frac{a{k}_3}{k_2+1}\right)}^m\\ \times {\displaystyle \sum _{l=0}^m}{k}_2^l{C}_m^l\frac{{\left(\alpha +1\right)}_l}{{\left(\beta +1\right)}_l}{}_{1}F_1\left(\alpha +1+l,\beta +1+l;{\omega }_2\right)\\ +\frac{a{k}_2{k}_3}{k_{\text{on}}+a{k}_1}\frac{\alpha }{\beta }{A}_0\frac{1}{\left(m-1\right)!}{\left(\frac{a{k}_3}{k_2+1}\right)}^{m-1}\\ \times {\displaystyle \sum _{l=0}^{m-1}}{k}_2^l{C}_{m-1}^l\frac{{\left(\alpha +1\right)}_l}{{\left(\beta +1\right)}_l}{}_{1}F_1\left(\alpha +1+l,\beta +1+l;{\omega }_2\right);\\ P_1\left(m\right)=A_0\frac{1}{m!}{\left(\frac{a{k}_3}{k_2+1}\right)}^m\\ \times {\displaystyle \sum _{l=0}^m}{k}_2^l{C}_m^l\frac{{\left(\alpha \right)}_l}{{\left(\beta \right)}_l}{}_{1}F_1\left(\alpha +l,\beta +l;{\omega }_2\right)& (3)\end{array}$

where $\alpha =a{k}_2{k}_3/{(k_2+1)}^2+(k_{\text{on}}+k_{\text{off}}+a{k}_1)/(k_2+1)$,β = α + 1, and ${\omega }_2=-a{k}_2{k}_3/{(k_2+1)}^2$. Here, all the parameters are normalized by k₄, i.e., k_i = k_i/k₄ which is different with Equation (2). $C_m^l$ is the binomial coefficient in Equation (3), (γ)_l is the Pochhammer symbol defined as (γ)_l = Γ(γ + l)/Γ(γ) with Γ(γ) being the Gamma function, and ₁F₁(α, β; ω₂) is a confluent hypergeometric function (Huang et al., 2015). A₀ is the normalization constant as follows

$\begin{array}{l}{A}_0^{-1}=e^{\frac{a{k}_3}{k_2+1}}[\frac{1}{k_{\text{on}}+a{k}_1}\left(\frac{a{k}_2{k}_3}{k_2+1}+k_{\text{off}}\right){}_{1}F_1(\alpha ,\beta ;{\omega }_1)\\ -\frac{1}{k_{\text{on}}+a{k}_1}\frac{a{k}_2{k}_3}{k_2+1}\frac{\alpha }{\beta }{}_{1}F_1(\alpha +1,\beta +1;{\omega }_1)\\ +\frac{a{k}_2{k}_3}{k_{\text{on}}+a{k}_1}\frac{\alpha }{\beta }{}_{1}F_1(\alpha +1,\beta +1;{\omega }_1)]+e^{\frac{a{k}_3}{k_2+1}}{}_{1}F_1(\alpha ,\beta ;{\omega }_1)& (4)\end{array}$

where ω₁ is a constant with the expression ${\omega }_1=a{k}_2{k}_3/(k_2+1)-a{k}_2{k}_3/{(k_2+1)}^2$. The details of analysis for chemical master equations are presented in the Supplementary Material of this paper. The above results will be checked through the structure of our model in the following part.

The character of signaling cascades in the development of skeletal muscle can provide us with information to verify our analytical solutions. The changes of ROCK, as an upstream component of the signaling cascades, follow a basic process of birth and death. The statistical law of a birth and death process is that the probability distribution about m is a standard Poisson distribution. In order to test this, we calculate the total probability P(m) which is provided by P₀(m) + P₁(m). This represents the statistical law of ROCK and can be simplified from Equation (3) as:

$\begin{array}{ll}P(m)=\frac{1}{m!}{(a{k}_3)}^m{e}^{-a{k}_3}.& (5)\end{array}$

This is a standard Poisson distribution. Furthermore, ROCK should be in its steady state m_s = ak₃ most of the time. This means that the peaks of P₀(m) and P₁(m) focus on m_s. According to Equation (3), the values of P₀(m) and P₁(m) with different parameters are computed through “Mathematica” which are shown in Figure 3. It is obvious that the peaks of distributions occur at the steady state value of m_s. It's worth mentioning that the roughness of the curves shown in Figures 3a,c,d is caused by the computational accuracy of “Mathematica” when it is used to calculate the confluent hypergeometric function rather than biological or physical factors. We can also get these curves through “Matlab” which appear very smooth. However, compared with “Mathematica”, “Matlab” fails to calculate confluent hypergeometric function when m is too large. In order to verify these curves, the corresponding curves with the same parameters obtained by Monte Carlo simulation are shown in Figures S1a,b. It is obviously that the curves in Figures S1a,b closely resemble the ones in Figures 3a,b, respectively.

FIGURE 3

Figure 3. Distributions of probability as a function of m. The synthetic rate of ROCK k₃ = 5 in (a,b) and k₃ = 1.5 in (c,d). The value of other parameters can be found in Table S1. P₀(m) is the probability of the gene “off” state indicated with orange curve, P₁(m) is the probability of gene “on” state indicated with blue curve.

The above discussion confirms the reliability of our theoretical results in Equation (3). Based on those results, we will try to explore the selectivity of different pathways and energy consumption in the following section.

3.2. Stochastic Dynamics of Gene Switching in DDL

The steady-state probability distributions with different stimulation strength are displayed in Figure 3. ${\sum }_m{P}_0(m)$ and ${\sum }_m{P}_1(m)$ which correspond to the areas under the curves of P₀(m) and P1(m) are respectively the probabilities of the gene's “off” state and “on” state. We use them to define two gene modes: Mode I and Mode II. Mode I denotes that the probability of the gene's “off” state is larger than the one of the “on” state, and Mode II denotes that the probability of the gene 's “on” state is larger than the one of the “off” state. As shown in Figures 3a,b, when k₃ = 5, the gene's state is Mode I. Conversely, if k₃ = 1.5, the gene's state becomes Mode II shown in Figures 3c,d. Compared with k₃, even if the parameter a varies widely, the gene's state mode is not changed. Figure 4A shows the areas of these two modes. It is obviously that the boundary between them is almost a horizontal line where k₃ ≈ 2.1. This means that the mode of the gene's state is determined primarily by k₃ and the parameter a has little effect on the selection of gene's state modes. Since k₃ and a represent the synthesis rate of the negative controller ROCK and the strength of external stimulations, respectively, the gene's state mode depends almost exclusively on the synthesis rate of the negative controller ROCK rather than the external stimulations.

FIGURE 4

Figure 4. The influences of k₃ and a on gene's “off” and “on” states. (A) The areas of Mode I and Mode II in a-k₃ plane; (B) The heat map of gene's state dominance factor δ in a-k₃ plane. The black line is the boundary between Mode I and Mode II; (C) The curve between δ and a in Mode I; (D) The curve between δ and a in Mode II. The values of other parameters are listed in Table S1.

By the definition of Mode I, the gene is more likely to be in the “off” state than the “on” state. In other words, the ‘off” state of gene is dominant in Mode I. Similarly, the “on” state of gene is dominant in Mode II. To quantify the dominance of the gene's state, we define a gene's state dominance factor δ = |P_0,max − P_1,max|/P_1,max. The larger peak in P_i(m) (i = 1, 2) means more dominance as shown in Figure 3. Therefore, δ can be used to measure the dominance of the gene's state. Next, we will discuss the influences of the synthesis rate of the negative controller ROCK and the external stimulations on gene's state dominance. The values of δ in a-k₃ plane correspond to the scale values of the color bar in Figure 4B. As shown in Figure 4B, with increasing k₃ (i.e., the synthesis rate of the negative controller ROCK), δ (i.e., the dominance of the gene's “off” state) increases in Mode I. Conversely, as k₃ increases, δ decreases in Mode II. That suggests that when the synthesis rate of the negative controller ROCK increases, the dominance of the gene's “off” state is gradually weakening in Mode II, until the mode of the gene's state changes to Mode I and the dominance of the gene's “on” state increases gradually. Compared with k₃, the influence of the strength of external stimulations RhoA^* (i.e., a) on the dominance of the gene's state (i.e., δ) is not obvious enough in Figure 4B. Therefore, the relation between a and δ with different parameters k₃ which correspond respectively to Mode I and Mode II is shown in Figures 4C,D (Note that we also use Monte Carlo simulation to verify the trend of δ with a in Mode I which is shown in Figure S2). Similar to the case of k₃, with increasing the strength of external stimulations RhoA^* (i.e., a), the dominance of the gene's state (i.e., δ) increases in Mode I and decreases in Mode II. That means that although external stimulations RhoA^* has little effect on the selection of gene's state modes, it can fine-tune the dominance of the gene's state in its respective modes.

In summary, the core factor for the stochastic dynamics of gene switching in DDL-type biochemical networks is identified in our work. Reaction rates are responsible for the selectivity between different gene's states, while the external signal stimulation fine-tunes the choice in its respective modes. The cooperation between signals and network maintains the vital process in an orderly manner.

3.3. Energy Dissipation in DDL

It is intuitively obvious that living biochemical systems need free energy (Gui et al., 2016, 2018). From the viewpoint of thermodynamics, gene expression is essentially a non-equilibrium process due to feedback or feedforward regulation that breaks detailed balance and thus necessarily consumes energy (Lu et al., 2017). But how is energy actually utilized during the regulation of gene expression in the development of skeletal muscle? To our knowledge, few works have touched upon this point. The composition of total energy consumption may help us grasp the selection mechanism between different biochemical processes.

From the definition of entropy $S(P_i)=k_{B}T{\sum }_i{P}_{i}ln{P}_i$, the entropy production rate ε_p(t) is given as follows (Ge and Qian, 2010):

$\begin{array}{ll}{\epsilon }_P(t)={\displaystyle \sum _{i,j}}(P_i(t)q_{ij}-P_j(t)q_{ji})ln\left(\frac{q_{ij}}{q_{ji}}\right).& (6)\end{array}$

P_i(t) is the probability of the system in state i at time t, while q_ji is the transport rate from state j to i. k_BT is set to be 1 for convenience in our work. ε_p(t) is the sum of energy dissipated in the biochemical network. Furthermore, the entropy production rate for a non-equilibrium steady-state system can be calculated as

$\begin{array}{ll}EP={\displaystyle \sum _{(\sigma ,{\sigma }^{\prime })}}P(\sigma )k(\sigma ,{\sigma }^{\prime })log\left(\frac{k(\sigma ,{\sigma }^{\prime })}{k({\sigma }^{\prime },\sigma )}\right)& (7)\end{array}$

where k(σ, σ′) is the transition probability from state σ to σ′.

Considering the detailed biochemical reactions in our model, we derive the EP of the DDL network as

$\begin{array}{l}EP={\sum }_m[P_0(m)a{k}_3\ln \left(\frac{a{k}_3}{m+1}\right)+P_0(m)m\ln \left(\frac{m}{a{k}_3}\right)\\ +P_0(m)(a{k}_1+k_{\text{on}})\ln \left(\frac{a{k}_1+k_{\text{on}}}{k_{\text{off}}+m{k}_2}\right)\\ +P_1(m)a{k}_3\ln \left(\frac{a{k}_3}{m+1}\right)+P_1(m)m\ln \left(\frac{m}{a{k}_3}\right)\\ +P_1(m)(k_{\text{off}}+m{k}_2)\ln \left(\frac{k_{\text{off}}+m{k}_2}{a{k}_1+k_{\text{on}}}\right)].& (8)\end{array}$

Next, we discuss the influences of the strength of external stimulation on the total energy dissipation EP in Mode I and Mode II respectively. As shown in Figure 5A, EP increases and its rate of increase decreases with increasing a in Mode I. This suggests that when the strength of external stimulation increases, the system consumes more and more energy to response it. Compared with the small strengthen of external stimulation, there is less growth of energy for the system corresponding to the large one in Mode I. In contrast to Mode I, both EP and its reduction rate decrease in Mode II as a increases (Figure 5B). This means although the system requires less and less energy with the increase of the strengthen of external stimulation, a small amount of energy is still needed to sustain it in Mode II.

FIGURE 5

Figure 5. The influences of the strength of external stimulation on energy dissipation in Mode I and II. The values of all parameters are listed in Table S1. (A,B) The relation curves between total energy dissipation EP and a in Mode I and II; (C,D) The relation curves between the percentages of energy dissipation in the synthesis-degradation process of ROCK (i.e., EP_m) in total energy dissipation (i.e., EP) and a in Mode I and II.

Since the control of MyoD gene expression can be divided into two parts: one is the synthesis and degradation of ROCK and the second is the switching of gene state, the process of energy dissipation can be decomposed into three state transitions: (m − 1, off) ⇌ (m, off) ⇌ (m + 1, off), (m − 1, on) ⇌ (m, on) ⇌ (m + 1, on) and (m, on) ⇌ (m, off). Here, (m, on) represents the state when the gene is “on” and the level of ROCK is m. Note that the first two formulas are related to the synthesis and degradation of ROCK. According to these three state transitions, the total energy dissipation EP is decomposed into three terms as follows:

$\begin{array}{l}\begin{array}{l}E{P}_1={\displaystyle \sum _m}[P_0(m)a{k}_3\log \left(\frac{a{k}_3}{m+1}\right)+P_0(m)m\log \left(\frac{m}{a{k}_3}\right)]\hfill \end{array}\\ \begin{array}{l}E{P}_2={\displaystyle \sum _m}[P_1(m)a{k}_3\log \left(\frac{a{k}_3}{m+1}\right)+P_1(m)m\log \left(\frac{m}{a{k}_3}\right)]\hfill \end{array}\\ \begin{array}{l}\begin{array}{l}E{P}_3={\displaystyle \sum _m}[P_0(m)a{k}_1+k_{on}\log \left(\frac{a{k}_3+k_{on}}{k_{\text{off}}+m{k}_2}\right)\\ +P_1(m)(k_{\text{off}}+m{k}_2)\log \left(\frac{k_{\text{off}}+m{k}_2}{a{k}_1+k_{on}}\right)].\end{array}\hfill \end{array}& (9)\end{array}$

Since the synthetic rate of ROCK k₃ which has been normalized by its degradation rate k₄ just appears in the formulas of EP₁ and EP₂, we define EP_m = EP₁ + EP₂ and use it to represent the energy dissipation in the synthesis-degradation process of ROCK. In the following part, we will study the allocation of total energy dissipation corresponding to different strengthen of external stimulations by the comparison between EP and EP_m. The Figure S3 shows the trends of the total energy dissipation (i.e., EP) and the energy dissipation in the synthesis-degradation process of ROCK (i.e., EP_m) with the increase of the strengthen of external stimulations (i.e., a) in Mode I and II. It is obviously that their trends are consistent as a increases in its respective modes. Specifically, EP and EP_m increase simultaneously in Mode I and decrease simultaneously in Mode II with the enhancing of the strengthen of external stimulations (i.e., a). Moreover, the difference between EP and EP_m diminishes both in Mode I and Mode II when a increases. The percentages of energy dissipation in the synthesis-degradation process of ROCK (i.e., EP_m) in total energy dissipation (i.e., EP) shown in Figures 5C,D are more than 60% both in Mode I and Mode II with different strengthen of external stimulations. In other words, the synthesis-degradation process of ROCK consume more energy than the third state transitions (i.e., EP₃) with different strengthen of external stimulations. Furthermore, Figures 5C,D also show that the percentage of energy which is consumed by the synthesis-degradation process of ROCK increases as the strengthen of external stimulations increases until almost no energy is consumed in process of the third state transitions both in Mode I and Mode II.

4. Discussion and Conclusion

In our work, a double-deck loop model is constructed. Due to the stochastic nature of bio-processes (Wang et al., 2017; Yao et al., 2018a,b), we have calculated the steady-state probability distributions of ROCK protein through the method of probability-generating functions for chemical master equations. The crucial factors in the stochastic dynamics of gene switching are identified. It is found that the weights between different pathways (i.e., the internal reaction rates) in DDL are the key point governing the state of gene switching, while an external stimulus fine-tunes this choice preference. Furthermore, the energy consumption in DDL is also discussed. Our results show that most of the energy is required for synthesis and degradation of ROCK, however, a very small amount of energy consumption is required for the basic transition processes of downstream genes between “on” and “off” states. This is because the ROCK processes are not in equilibrium and do not follow detailed balance. But the inter-conversion between “on” state and “off” state is indeed in equilibrium and follows detailed balance. In other words, the two terms in EP₃ defined in Equation (8) cancel each other out because of detailed balance. The theoretical findings about selectivity between different gene states and energy dissipation will be advantageous for our understanding of cell fate determination. Our next steps are to conduct closely related experiments about the development of skeletal muscle and to combine our theoretical study with experimental observations and data.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

QL, MY, and JS conceived and designed the work. QL, JS, FY, LY, YG, and RG carried out computer implementation and theoretical analysis. QL, JS, FY, LY, YG, and RG interpreted the simulation results. MY supervised the project. QL, FY, LY, YG, JS, and RG wrote the original manuscript. QL, MY, and JS contributed to the writing of final manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Grant Nos. 11675060, 91730301, 11804111, and 11805091).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00676/full#supplementary-material

References

Aelst, L. V., and D'Souza-Schorey, C. (1997). Rho GTPases and signaling networks. Genes Dev. 11, 2295–2322. doi: 10.1101/gad.11.18.2295

PubMed Abstract | CrossRef Full Text | Google Scholar

Beqaj, S., Jakkaraju, S., Mattingly, R. R., Pan, D., and Schuger, L. (2002). High RhoA activity maintains the undifferentiated mesenchymal cell phenotype, whereas RhoA down-regulation by laminin-2 induces smooth muscle myogenesis. J. Cell Biol. 156, 893–903. doi: 10.1083/jcb.200107049

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Wang, H., Ouyang, Q., and Tu, Y. (2015). The free-energy cost of accurate biochemical oscillations. Nat. Phys. 11, 772–778. doi: 10.1038/nphys3412

PubMed Abstract | CrossRef Full Text | Google Scholar

Castellani, L., Salvati, E., Alema, S., and Falcone, G. (2006). Fine regulation of RhoA and rock is required for skeletal muscle differentiation. J. Biol. Chem. 281, 15249–15257. doi: 10.1074/jbc.M601390200

PubMed Abstract | CrossRef Full Text | Google Scholar

Charrasse, S., Comunale, F., Grumbach, Y., Poulat, F., Blangy, A., and Gauthier-Rouviere, C. (2006). RhoA GTPase regulates M-cadherin activity and myoblast fusion. Mol. Biol. Cell 17, 749–759. doi: 10.1091/mbc.e05-04-0284

PubMed Abstract | CrossRef Full Text | Google Scholar

Charrasse, S., Meriane, M., Comunale, F., Blangy, A., and Gauthier-Rouviere, C. (2002). N-cadherin-dependent cell-cell contact regulates rho GTPases and β-catenin localization in mouse C2C12 myoblasts. J. Cell Biol. 158, 953–965. doi: 10.1083/jcb.200202034

PubMed Abstract | CrossRef Full Text | Google Scholar

Cloutier, M., Tremblay, M., and Piedboeuf, B. (2010). ROCK2 is involved in accelerated fetal lung development induced by in vivo lung distension. Pediatr. Pulmonol. 45, 966–976. doi: 10.1002/ppul.21266

PubMed Abstract | CrossRef Full Text | Google Scholar

Endres, R. G. (2017). Entropy production selects nonequilibrium states in multistable systems. Sci. Rep. 7:14437. doi: 10.1038/s41598-017-14485-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Fei, C., Cao, Y., Ouyang, Q., and Tu, Y. (2018). Design principles for enhancing phase sensitivity and suppressing phase fluctuations simultaneously in biochemical oscillatory systems. Nat. Commun. 9:1434. doi: 10.1038/s41467-018-03826-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, H., and Qian, H. (2010). Physical origins of entropy production, free energy dissipation, and their mathematical representations. Phys. Rev. E 81:051133. doi: 10.1103/PhysRevE.81.051133

PubMed Abstract | CrossRef Full Text | Google Scholar

Gui, R., Li, Z., Hu, L., Cheng, G., Liu, Q., Xiong, J., et al. (2018). Noise decomposition algorithm and propagation mechanism in feed-forward gene transcriptional regulatory loop. Chin. Phys. B 27, 92–103. doi: 10.1088/1674-1056/27/2/028706

PubMed Abstract | CrossRef Full Text | Google Scholar

Gui, R., Liu, Q., Yao, Y., Deng, H., Ma, C., Jia, Y., et al. (2016). Noise decomposition principle in a coherent feed-forward transcriptional regulatory loop. Front. Physiol. 30:600. doi: 10.3389/fphys.2016.00600

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, L., Yuan, Z., Yu, J., and Zhou, T. (2015). Fundamental principles of energy consumption for gene expression. Chaos 25:123101. doi: 10.1063/1.4936670

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaibuchi, K., Kuroda, S., and Amano1, M. (1999). Regulation of the cytoskeleton and cell adhesion by the rho family GTPases in mammalian cells. Annu. Rev. Biochem. 68, 459–486. doi: 10.1146/annurev.biochem.68.1.459

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, G., Sartori, P., Neumann, S., Sourjik, V., and Tu, Y. (2012). The energy–speed–accuracy trade-off in sensory adaptation. Nat. Phys. 8, 422–428. doi: 10.1038/nphys2276

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, G., and Tu, Y. (2013). The cost of sensitive response and accurate adaptation in networks with an incoherent type-1 feed-forward loop. J. R. Soc. Interface 10:20130489. doi: 10.1098/rsif.2013.0489

PubMed Abstract | CrossRef Full Text | Google Scholar

Leung, T., Manser, E., Tan, L., and Lim, L. (1995). A novel serine/threonine kinase binding the Ras-related RhoA GTPase which translocates the kinase to peripheral membranes. J. Biol. Chem. 270, 29051–29054. doi: 10.1074/jbc.270.49.29051

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, L., Jia, Y., Liu, W., and Yang, L. (2017). Mixed stimulus-induced mode selection in neural activity driven by high and low frequency current under electromagnetic radiation. Complexity. 2017:7628537. doi: 10.1155/2017/7628537

CrossRef Full Text | Google Scholar

Matsui, T., Amano, M., Yamamoto, T., Chihara, K., Nakafuku, M., Ito, M., et al. (1996). Rho-associated kinase, a novel serine/threonine kinase, as a putative target for small GTP binding protein rho. EMBO J. 15, 2208–2216. doi: 10.1002/j.1460-2075.1996.tb00574.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Matsuoka, T., and Yashiro, M. (2014). Rho/rock signaling in motility and metastasis of gastric cancer. World J. Gastroenterol. 20, 13756–13766. doi: 10.3748/wjg.v20.i38.13756

PubMed Abstract | CrossRef Full Text | Google Scholar

Meriane, M., Roux, P., Primig, M., Fort, P., and Gauthier-Rouviere, C. (2000). Critical activities of rac1 and Cdc42Hs in skeletal myogenesis: antagonistic effects of JNK and p38 pathways. Mol. Biol. Cell 11, 2513–2831. doi: 10.1091/mbc.11.8.2513

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, H. (2007). Phosphorylation energy hypothesis: open chemical systems and their biological functions. Annu. Rev. Phys. Chem. 58, 113–142. doi: 10.1146/annurev.physchem.58.032806.104550

PubMed Abstract | CrossRef Full Text | Google Scholar

Sordella, R., Jiang, W., Chen, G.-C., Curto, M., and Settleman, J. (2003). Modulation of rho gtpase signaling regulates a switch between adipogenesis and myogenesis. Cell 113, 147–158. doi: 10.1016/S0092-8674(03)00271-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, S. H., Huang, P. H., Peng, Y. J., Chang, W. C., Tsai, H. Y., Leu, H. B., et al. (2013). Zoledronate attenuates angiotensin ii-induced abdominal aortic aneurysm through inactivation of Rho/ROCK-dependent JNK and NF-kB pathway. Cardiovasc. Res. 100, 501–510. doi: 10.1093/cvr/cvt230

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L.-F., Qiu, K., and Jia, Y. (2017). Effects of time delays in a mathematical bone model. Chin. Phys. B 26, 251–257. doi: 10.1088/1674-1056/26/3/030503

CrossRef Full Text | Google Scholar

Wang, S., Lan, Y., and Tang, L. (2015). Energy dissipation in an adaptive molecular circuit. J. Stat. Mech. 2015:P07025. doi: 10.1088/1742-5468/2015/00/P07025

CrossRef Full Text | Google Scholar

Wei, L., Roberts, W., Wang, L., Yamada, M., Zhang, S., Zhao, Z., et al. (2001). Rho kinases play an obligatory role in vertebrate embryonic organogenesis. Development 128, 2953–2962.

PubMed Abstract | Google Scholar

Wei, L., Zhou, W., Croissant, J. D., Johansen, F.-E., Prywes, R., Balasubramanyam, A., et al. (1998). Rhoa signaling via serum response factor plays an obligatory role in myogenic differentiation. J. Biol. Chem. 273, 30287–30294. doi: 10.1074/jbc.273.46.30287

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, L., Zhou, W., Wang, L., and Schwartz, R. J. (2000). β₁-integrin and pi 3-kinase regulate rhoa-dependent activation of skeletal α-actin promoter in myoblasts. J. Physiol. Heart Circ. Physiol. 278, H1736–H1743. doi: 10.1152/ajpheart.2000.278.6.H1736

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, Y., Cao, W., Pei, Q., Ma, C., and Yi, M. (2018a). Break up of spiral wave and order-disorder spatial pattern transition induced by spatially uniform cross-correlated sine-wiener noises in a regular network of hodgkin-huxley neurons. Complexity 2018:8793928. doi: 10.1155/2018/8793298

CrossRef Full Text | Google Scholar

Yao, Y., Yang, L., Wang, C., Liu, Q., Gui, R., Xiong, J., et al. (2018b). Subthreshold periodic signal detection by bounded noise-induced resonance in the Fitzhugh–nagumo neuron. Complexity 2018:5632650. doi: 10.1155/2018/5632650

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, D., Cao, Y., and Ouyang, Q. Y. T. (2020). The energy cost and optimal design for synchronization of coupled molecular oscillators. Nat. Phys. 16, 95–100. doi: 10.1038/s41567-019-0701-7

CrossRef Full Text | Google Scholar