www.ietdl.org Published in IET Systems Biology Received on 30th October 2013 Revised on 23rd December 2013 Accepted on 22nd January 2014 doi: 10.1049/iet-syb.2013.0050
Special Issue on Selected papers from The 7th IEEE International Conference on Systems Biology (ISB 2013) ISSN 1751-8849
Overshoot in biological systems modelled by Markov chains: a non-equilibrium dynamic phenomenon Chen Jia1,2, Minping Qian1, Daquan Jiang1,3 1
LMAM, School of Mathematical Sciences, Peking University, Beijing 100871, People’s Republic of China Beijing International Center for Mathematical Research, Beijing 100871, People’s Republic of China 3 Center for Statistical Science, Peking University, Beijing 100871, People’s Republic of China E-mail: [email protected] 2
Abstract: A number of biological systems can be modelled by Markov chains. Recently, there has been an increasing concern about when biological systems modelled by Markov chains will perform a dynamic phenomenon called overshoot. In this study, the authors found that the steady-state behaviour of the system will have a great effect on the occurrence of overshoot. They showed that overshoot in general cannot occur in systems that will finally approach an equilibrium steady state. They further classified overshoot into two types, named as simple overshoot and oscillating overshoot. They showed that except for extreme cases, oscillating overshoot will occur if the system is far from equilibrium. All these results clearly show that overshoot is a non-equilibrium dynamic phenomenon with energy consumption. In addition, the main result in this study is validated with real experimental data.
1
Introduction
Recent advances in the single-cell and single-molecule experiments have shown that biological systems in living cells are inherently stochastic [1–8]. It is widely observed that a number of biological systems can transition stochastically among multiple states. These systems are often mathematically characterised by the stochastic model of Markov chains, or in the language of physics, master equations. Typical examples of Markov chain systems include enzyme kinetics with the Michaelis-Menten mechanism [9] or the general modifier mechanism of Botts and Morales [10], phosphorylation–dephosphorylation kinetics of proteins [11], conformational changes of receptors [12] and stochastic state transitions of cell populations [13]. In recent years, scientists have become increasingly concerned about a kind of dynamic phenomenon in biological systems which is known as overshoot or biochemical adaptation [14–19]. Overshoot of a biological system refers to the dynamic phenomenon that the output, which is usually the average of an observable, of the system varies with time in a non-monotonic way and exceeds both its initial and steady-state values over a certain period of time. An intuitive description of overshoot is depicted in Fig. 1, where the output of the system first rises to a peak and then declines to a lower plateau. Overshoot is widely observed in numerous biological systems. Typical examples of overshoot include chemotaxis of bacteria [20, 21], osmotic sensing in yeast [22], calcium dose response of the inositol trisphosphate receptors [23, 24] and the ryanodine receptors [25–27], and hormone dose response of the hormone receptors [28]. Overshoot is an 138 & The Institution of Engineering and Technology 2014
important biological function possessed by many living systems. It allows the system to detect environmental changes more accurately, enables the system to respond to environmental fluctuations more rapidly and protects the system from irreversible damages caused by unfavourable conditions. To better understand how overshoot is achieved in biochemical feedback networks, several research groups studied the relationship between overshoot and the network topology [16–18]. Tang and co-workers [18] searched all possible three-node network topologies and found that overshoot is most likely to occur in two types of networks: the negative feedback loop and the incoherent feedforward loop. Tu and co-workers [19] studied the stochastic dynamics of the negative feedback loop in detail and found that the negative feedback mechanism breaks detailed balance, and thus always operates out of equilibrium with energy dissipation. These two works give us a hint that overshoot may tend to occur in biological systems which are far from equilibrium. To avoid misinterpretation, we point out here that the word ‘equilibrium’ appearing in this paper is referred to as the concept of equilibrium steady state in statistical physics, where the steady state of a system is called equilibrium (non-equilibrium) if the detailed balance condition is satisfied (broken). In fact, the concept of overshoot has long been suggested and studied in control theory [29, 30], electronics [31] and signal processing. In these disciplines, overshoot was studied in various deterministic systems composed of several ordinary differential equations and some overshoot conditions were provided. However, the traditional study of overshoot under the framework of deterministic systems does not capture the physical essence of overshoot. Many IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org
Fig. 1 Overshoot in biological systems Overshoot refers to the dynamic phenomenon that the output of the system varies with time in a non-monotonic way and exceeds both its initial and steady-state values over a certain period of time
important topics about overshoot, such as its relation to the breakdown of detailed balance and to energy dissipation, can only be clearly seen in stochastic systems. Up till now, there is still a lack of a general analysis of overshoot, which captures its physical essence, in biological systems which are inherently stochastic. In this paper, we presented a general analysis of overshoot in biological systems which can be modelled by Markov chains with two, three or multiple states. We found that the steady-state behaviour of the system will have a great effect on the occurrence of overshoot. We made it clear that overshoot in general cannot occur in systems which will finally approach an equilibrium steady state. This explains why overshoot can only be observed in systems with three or more states and cannot be observed in systems with only two states. We further classified overshoot into two types, named as simple overshoot and oscillating overshoot. We found that except for extreme cases, oscillating overshoot will occur in systems far from equilibrium. All the above results clearly show that overshoot is a non-equilibrium dynamic phenomenon and thus sustained energy consumption is required for the system to perform this important biological function. In addition, we used the experimental data of SUM159 human breast cancer cell line to validate the main theoretical result in this paper.
2
where pi(t) is the probability of the system being in state i at time t. We denote by πi the initial probability of state i, that is, the probability of the system being in state i at time t = 0 and denote by μi the steady-state probability of state i, that is, the probability of the system being in state i when time t is sufficiently large. Then f(0) = Ni=1 pi fi is the initial output of the system and f(1) = Ni=1 mi fi is the steady-state output of the system. According to Fig. 1, the system performs overshoot if and only if there exists a time t such that f (t) is larger than both f(0) and f(∞). Next, we shall use three specific examples to help the readers understand the Markov chain model discussed above. Example 1: We consider the well-known Michaelis-Menten enzyme kinetics k−1
k2
E + S O ES E + P k1
(2)
where E is an enzyme involved in converting the substrate S into the product P. If there is only one enzyme molecule, it
Model
In this paper, we consider biological systems that can be mathematically modelled as continuous-time Markov chains with multiple states. The characteristic property of the Markov chain is that it retains no memory of where it has been in the past. This means that where the system will go next only depends on its current state, but not depends on its prior states. We assume that the Markov chain system can be found in N states, 1, 2, …, N, and thus the system can transition stochastically among these states. The specific meaning of these states varies from systems to systems. Each state can represent a binding state of an enzyme molecule [10], a conformational state of a receptor molecule [14, 15, 25–27], a cellular state of a cell population [13] and so on. We further assume that the system has an observable f. If the system is in state i, the observation of the system will be fi. The output f (t) of the system at time t is then the weighted average of the observations of all states
f(t) =
N
pi (t)fi
i=1
IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
(1)
Fig. 2 Markov chain systems with two, three or multiple states a Markov chain systems with two states b Markov chain systems with three states c–e Examples of biological systems that can be modelled by Markov chains c The catalytic cycle of the Michaelis-Menten enzyme mechanism, where E is the enzyme, S is the substrate, ES is the enzyme-substrate complex and [S] is the concentration of the substrate d Cell-state dynamics of human breast cancer cells An individual breast cancer cell can transition stochastically among three differentiation states: stem-like (S), basal (B) and luminal (L) states The transition rates between states are estimated from the data set of the SUM 159 human breast cancer cell line [13] The transition rates are shown per cell division e The MWC allosteric model which describes the conformational changes of receptors with n − 1 identical subunits Each receptor can bind to an attractant (small solid circle) If the attractant binding site is occupied, all subunits will switch together from the inactivated configuration (solid square) to the activated configuration (large solid circle) In addition, each subunit can bind to a ligand (solid ellipse) According to whether all subunits are activated or inactivated and the number of subunits which have bound to the ligand, each receptor may convert among 2n possible states The upper n states are activated states and the lower n states are inactivated states 139
& The Institution of Engineering and Technology 2014
www.ietdl.org may convert between two states: the free enzyme E and the enzyme-substrate complex ES. Then from the perspective of a single enzyme molecule, the Michaelis-Menten kinetics can be represented by the catalytic cycle illustrated in Fig. 2c. We denote by E0 = [E] + [ES] the total enzyme concentration. Then pE(t) = [E]/E0 and pES(t) = [ES]/E0 represent the probability of a single enzyme molecule being in state E and state ES at time t, respectively. Mathematically, the catalytic cycle illustrated in Fig. 2c is nothing but a Markov chain with two states. According to the law of mass action, the probability flux from state E to state ES at time t is k1[S]pE(t), and the probability flux from state ES to state E at time t is (k−1 + k2)pES(t), where we add k−1 and k2 since there are two ways of transition from state ES to state E. Based on the expressions of the probability fluxes, we easily see that the transition rate from state E to state ES is k1[S] and the transition rate from state ES to state E is k−1 + k2. In the Michaelis-Menten enzyme system, the output f (t) is often chosen as the instantaneous rate of the product formation
f(t) =
d[P] = k2 [ES] = k2 E0 pES (t) dt
(3)
where the second equality is because of the law of mass action. Example 2: Recent research shows that human breast cancer cells within individual tumours can exist in three differentiation states that differ in functional attributes: stem-like (S), basal (B) and luminal (L) states. An individual breast cancer cell can transition stochastically among these three states. Mathematically, Lander and co-workers [13] modelled the cell-state transitions and dynamics in the breast cancer cell population as a three-state Markov chain depicted in Fig. 2d. In the breast cancer cell system, the output f (t) is often chosen as the proportion of stem-like, basal or luminal cells
f(t) = pi (t),
i = S, B, L
(4)
Example 3: Structural studies show that receptors in living cells are often protein complexes with multiple subunits. Mathematically, Monod, Wyman and Changeux modelled the allosteric protein interactions in receptors as a Markov chain with multiple states. The Monod–Wyman–Changeux (MWC) allosteric model assumes that each receptor consists of n − 1 identical subunits, each of which can switch between two configurations, an activated one and an inactivated one. The MWC model further assumes that each receptor has an attractant binding site. Once an attractant binds to the receptor, all subunits will switch together from being inactivated to being activated. In addition, the MWC model assumes that each subunit has a ligand binding site. A ligand can bind to a subunit in either configuration, but the dissociation constants are different. According to whether all subunits are activated or inactivated and the number of subunits which have bound to the ligand, the conformational changes of each receptor can be modelled as a Markov chain with 2n states. The transition diagram of the MWC model when n = 5 is depicted in Fig. 2e, where each solid square represents an inactivated subunit and each large (solid) circle represents an activated subunit. Thus the upper n states in Fig. 2e are activated states and the lower n states are inactivated states. In the MWC model, the output 140 & The Institution of Engineering and Technology 2014
f(t) is often chosen as the average activity, that is, the sum of probabilities of those activated states f(t) =
n
pi (t)
(5)
i=1
3
Results
3.1 General analysis and overshoot in two-state systems We now use the theory of Markov chains to present a general analysis of overshoot. We know that the dynamics of the probability distribution p(t) = ( p1(t), …, pN(t)) of the Markov chain system is governed by the master equation of the matrix form dp(t) = p(t)Q dt
(6)
where Q = (qij )N×N is the transition rate matrix of the Markov chain system and qij is the transition rate from state i to state j. We assume without loss of generality that the transition rate matrix Q has N linear independent eigenvectors since any matrix can be approximated by a matrix satisfying this condition with arbitrarily high accuracy. Under this assumption, the output f (t) of the system is nothing but a linear combination of exponential functions
f(t) = p(t)f T = p(0)etQ f T =
N
ci eli t
(7)
i=1
where f = ( f1,…, fN) is a vector whose components are the observations of all states, λ1, …, λN are all eigenvalues of the transition rate matrix Q, and c1, …, cN are N constants. The well-known Perron–Frobenius theorem [32] in matrix theory claims that one of the eigenvalues of the transition rate matrix Q must be 0 and the real parts of other eigenvalues are all negative. In the following discussion, we always assume that λ1 = 0. Thus the output f (t) of the system can be rewritten that
f(t) = c1 +
N
ci eli t
(8)
i=2
We easily see that ci is real if λi is real, and ci and cj are conjugate complex numbers if λi and λj are conjugate complex numbers. An interesting phenomenon widely observed in experiments is that overshoot cannot be observed in Markov chain systems with only two states and can be observed in Markov chain systems with three or more states [33]. We now use the previous analysis to explain this phenomenon. According to (8), the output of a two-state system depicted in Fig. 2a is given by
f(t) = c1 + c2 el2 t
(9)
which is obviously a monotonic function since λ2 must be a real number. This suggests that a Markov chain system with only two states will never perform overshoot. IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org 3.2
Equilibrium and non-equilibrium systems
We have seen from previous discussions that Markov chain systems with only two states will never perform overshoot and overshoot can only occur in systems with three or more states. This raises a natural question of what is the essential difference between systems with two states and systems with three or more states. We noted that there is an apparent difference between systems with two states and systems with three or more states. A two-state system is always an equilibrium system which will finally approach an equilibrium steady state, whereas a system with three or more states may be a non-equilibrium system which will finally approach a non-equilibrium steady state [34–36]. Here, we have used the concepts of equilibrium and non-equilibrium steady states in non-equilibrium statistical physics, where the steady state of a system is called equilibrium if each pair of states i and j of the system satisfies the detailed balance condition
mi qij = mj q ji
(10)
(otherwise, the steady state of the system is called nonequilibrium) where μi is the steady-state probability of state i and qij is the transition rate from state i to state j. In the following discussion, a system which will finally approach an equilibrium (non-equilibrium) steady state is referred to as an equilibrium (non-equilibrium) system. From the point of view of statistical mechanics, an equilibrium system in the steady state is microscopic reversible and does not consume energy, whereas a non-equilibrium system in the steady state is microscopic irreversible and always consumes energy. Mathematically, whether a system is a non-equilibrium system or not is linked to the net flux the system [34]. To make the readers understand the concept of the net flux, we now limit our discussion to three-state systems depicted in Fig. 2b. In a three-state system, we denote by c the cycle 1 → 2 → 3 → 1 and denote by c− the reverse cycle 1 → 3 → 2 → 1. Let wc(t) denote the number of cycle c formed by the system up to time t. Then the flux wc of cycle c is defined to be wc (t) t1 t
wc = lim
(11)
We clearly see that the flux wc of cycle c represents the number of the forming of cycle c per unit time. Similarly, we can define the flux wc− for the reverse cycle c−. The net flux J of a three-state system is then defined to be J = wc − wc−, which represents the net number of the forming of cycle c per unit time. Interestingly, the net flux J of a three-state system can be represented by the steady-state probabilities and the transition rates as J = m1 q12 − m2 q21 = m2 q23 − m3 q32 = m3 q31 − m1 q13
(12)
This relation is a special case of the well-known circulation
a1 = −
decomposition theorem [34] in the Markov chain theory. From (10) and (12), we easily see that a three-state system is a non-equilibrium system if and only if the net flux J fails to be zero. This shows that the net flux J characterises how far the system is away from equilibrium. Generally speaking, equilibrium systems in the steady state do not consume energy, whereas non-equilibrium systems always consume energy to maintain non-zero net fluxes. Finally, we state a simple but important fact about equilibrium systems which will be frequently used in following discussions. In an equilibrium system, one eigenvalue of the transition rate matrix Q must be zero and all other eigenvalues must be negative real numbers. In other words, the eigenvalues λ1, …, λN of the transition rate matrix Q of an equilibrium system must satisfy
l1 = 0, l2 , . . . , lN , 0 3.3
Simple overshoot
We have seen from previous discussions that systems with only two states will always approach an equilibrium steady state and will never perform overshoot. However, the situation is totally different in systems with three states. According to (8), the output f (t) of a three-state system has the general form of
f(t) = c1 + c2 el2 t + c3 el3 t
(14)
If λ2 and λ3 are negative real numbers and the signs of c2 and c3 are not the same, then f(t) may not be a monotonic function and the system may perform overshoot. This type of overshoot is named as simple overshoot (Fig. 3a). To gain a deeper insight into simple overshoot, we consider three-state systems depicted in Fig. 2b starting from state 1 whose transition rate matrix has three real eigenvalues satisfying λ1 = 0 and λ2, λ3 < 0. To simplify notations, we make an additional assumption that the transition rates of the system satisfy q21 = q31. We impose this condition onto the system for two reasons. Firstly, this assumption ensures that the transition rate matrix has three real eigenvalues. Secondly, this assumption reduces the complexities of calculations and formulations to a remarkable extent. Under this assumption, the three eigenvalues of the transition rate matrix have explicit expressions of λ1 = 0, λ2 = −(q12 + q13 + q21) and λ3 = −(q23 + q32 + q21). Furthermore, the output f(t) of the system can be explicitly calculated as
f(t) = a1 + a2 el2 t + J a3 el3 t
(15)
where J is the net flux of the system, and α1, α2, α3 are three constants with the following expressions (see (16))
a2 =
(q12 + q13 )[(q12 − q32 )(f1 − f2 ) + (q13 − q23 )(f1 − f3 )] l2 (l2 − l3 ) (17)
a3 =
l2 (f2 − f3 ) q21 (l2 − l3 )
(18)
We are now in the position to consider when overshoot will
q21 l3 f1 + (q12 q31 + q12 q32 + q32 q13 )f2 + (q12 q23 + q13 q23 + q21 q31 )f3 l2 l3
IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
(13)
(16)
141
& The Institution of Engineering and Technology 2014
www.ietdl.org
Fig. 3 Overshoot in Markov chain systems with three states a Simple overshoot and oscillating overshoot in three-state systems If λ2 and λ3 are negative real numbers, then the system may perform simple overshoot If λ2 and λ3 are conjugate complex numbers, then the system will perform oscillating overshoot b Overshoot in the SUM 159 human breast cancer cell line The proportion of luminal cells experiences a transient increase after isolation of stem-like cells followed by a slow decrease to its steady-state value
occur in an equilibrium three-state system. We recall that a three-state system is an equilibrium system if and only if the net flux J is zero. According to (15), the output of an equilibrium three-state system can be simplified as
f(t) = c1 + c2 el2 t
(19)
which is obviously a monotonic function since λ2 is a real number. This shows that under mild conditions, an equilibrium system with three states will never perform overshoot. 3.4
conjugate complex eigenvalues when the net flux J is sufficiently large. To this end, let j1 = μ1q12, j2 = μ2q23 and j3 = μ3q31 be the clockwise probability fluxes of the system. According to (12), the transition rate matrix Q can be represented by μ1, μ2, μ3, j1, j2, j3 and J as ⎛ ⎜ ⎜ ⎜ Q=⎜ ⎜ ⎜ ⎝
Oscillating overshoot
We have seen from previous discussions that a three-state system may perform simple overshoot if λ2 and λ3 are negative real numbers. However, more interesting is the case when λ2 and λ3 are conjugate complex numbers. If we denote them by λ2 = λ + ωi and λ3 = λ − ωi, then the output f(t) of the system can be rewritten as
f(t) = c1 + 2|c2 |elt cos (vt + w)
(20)
where λ is a negative real number and j is the argument of c2. This expression clearly shows that after the output rises to a peak and returns toward its initial value, it may rise and decline again to form a damped oscillation. This type of overshoot is named as oscillating overshoot (Fig. 3a). Next, we shall present a detailed discussion of oscillating overshoot in equilibrium and non-equilibrium three-state systems. According to (13), the eigenvalues of the transition rate matrix of an equilibrium system are all real numbers. This shows that λ2 and λ3 cannot be conjugate complex numbers and thus oscillating overshoot will never occur in equilibrium three-state systems. The situation is totally different if we consider three-state systems which will finally approach a steady state far from equilibrium. In other words, we consider a three-state system with a sufficiently large net flux J, since the net flux characterises how far the system is away from equilibrium. We have seen that a three-state system will perform oscillating overshoot if the transition rate matrix has a pair of conjugate complex eigenvalues. What we need to do next is to study when the transition rate matrix has a pair of 142 & The Institution of Engineering and Technology 2014
−
j3 + j1 + J m1 j1 + J m2 j3 m3
j1 m1 j + j2 + J − 1 m2 j2 + J m3
⎞ j3 + J ⎟ m1 ⎟ ⎟ j2 ⎟ ⎟ m2 ⎟ j2 + j3 + J ⎠ − m3 (21)
Let pQ(λ) = det(λI − Q) be the characteristic polynomial of the transition rate matrix Q. Straightforward calculations show that p Q (l ) =
1 l(m1 m2 m3 l2 + Al + B) m1 m2 m3
(22)
where A = (m1 m2 + m2 m3 + m3 m1 )J + (m1 m2 (j2 + j3 ) + m2 m3 (j3 + j1 ) + m3 m1 (j1 + j2 ))
(23)
B = J 2 + (j1 + j2 + j3 )J + (j1 j2 + j2 j3 + j3 j1 )
(24)
and
We know from elementary algebra that the transition rate matrix Q has a pair of conjugate complex eigenvalues if and only if the discriminant Δ = A 2 − 4μ1μ2μ3B < 0. We easily calculate that D = dJ 2 + aJ + b
(25)
d = (m1 m2 + m2 m3 + m3 m1 )2 − 4m1 m2 m3
(26)
where
From (25), we make a crucial observation that when the net flux J is sufficiently large, the transition rate matrix Q has a IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org pair of conjugate complex eigenvalues if and only if δ is negative. The remaining question is to answer when δ is negative. We can prove that δ is a quantity satisfying −
1 1 ≤d≤ 27 16
(27)
We can further prove that δ attains its minimum of δmin = −1/ 27 when μ1 = μ2 = μ3 = 1/3, and attains its maximum of δmax = 1/16 when μ1 = μ2 = 1/2 and μ3 = 0. This clearly shows that the more uniform the three steady-state probabilities, μ1, μ2 and μ3, the smaller the value of δ. Thus except for the extreme case that the three steady-state probabilities are extremely scattered, δ is always negative and the system will perform oscillating overshoot when the net flux J is sufficiently large. This suggests that oscillating overshoot in general will occur in systems far from equilibrium. 3.5
Overshoot in systems with multiple states
In previous discussions, we are mainly concerned about when overshoot will occur in systems with two or three states. We have seen that a two-state system and an equilibrium three-state system in general cannot perform overshoot. This shows that overshoot is a non-equilibrium dynamic phenomenon in systems with two or three states. This raises a natural question of whether we can obtain the same conclusion in systems with multiple states. We point out that if a Markov chain system has more than three states, the mathematical complexity will become so high that it is almost impossible to present a complete discussion of overshoot. Thus in this section, we only focus on an important class of multiple-state systems, namely, the MWC allosteric model depicted in Fig. 2e, which is widely used in modelling the conformational changes of receptors in living cells. We see from Fig. 2e that some transition rates of the MWC model are regulated by the attractant concentration I, which is the input of the system. Let Ii and If be two input levels with 0 ≤ Ii < If. In experiments, we are always concerned about whether the system will perform overshoot in response to a step increase of the input from Ii to If. Interestingly, we can prove that if the MWC model is an equilibrium system, then it will never perform overshoot no matter what Ii and If are chosen. For a proof of this result, please see Section 5. This result clearly shows that an equilibrium system with multiple states in general cannot perform overshoot. All the results in the above sections indicate that overshoot, whether in systems with two, three or multiple states, is a non-equilibrium dynamic phenomenon, and thus sustained energy consumption is required for living systems to perform this important biological function. 3.6 Validation of main results with experimental data As mentioned in previous discussions, human breast cancer cells within individual tumours can transition stochastically among three differentiation states: stem-like (S), basal (B) and luminal (L) states. Mathematically, the cell-state transitions and dynamics in the breast cancer cell population can be modelled as a three-state Markov chain depicted in Fig. 2d. To validate our theoretical results with real experimental data, we shall use the data set of the SUM159 human breast cancer cell line from recent published work [13] to study overshoot in the breast cancer cell system depicted in Fig. 2d. IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
From the experimental data, we estimate the transition rates between any pair of differentiation states of the breast cancer cell system, as illustrated in Fig. 2d. For more details on the data usage and calculation, please see Section 5. The output f(t) of the breast cancer cell system is chosen to be the proportion of stem-like, basal, or luminal cells: f(t) = pi(t), i = S, B, L. Lander and co-workers found that if the stem-like cells are isolated at a particular time, then the breast cancer cell system will perform overshoot [13]. To see this, we depict the time course of the proportions of stem-like, basal and luminal cells in Fig. 3b, from which we see that the proportion of the luminal cells experiences a transient increase after isolation of the stem-like cells followed by a slow decrease to its steady-state value. In order to see whether overshoot in the breast cancer cell system is a non-equilibrium dynamic phenomenon, we need to calculate the net flux of the system. From the experimental data, the net flux of the breast cancer cell system is estimated as J = −0.0028 (see Section 5), which fails to be zero. This shows that overshoot in the SUM159 human breast cancer cell line is indeed a non-equilibrium dynamic phenomenon, which is consistent with the main theoretical result in this paper.
4
Discussion
Living systems are highly dissipative, exchanging materials and energy with their environments and consuming energy to carry out various important biological functions. If the exchange with the environment is sustained, then the living system always approaches a steady state far from equilibrium. Recent studies show that many important biological phenomena, such as coherence resonance in excitable systems [35], unidirectional movement of molecular motors [37] and switching behaviour of the general modifier mechanism of Botts and Morales [10], fail to occur in systems which will finally approach an equilibrium steady state. This suggests that sustained energy consumption is required for living systems to perform these important biological functions. In this paper, we presented a general discussion of the dynamic phenomenon of overshoot in biological systems modelled by Markov chains (master equations). We found that the steady-state behaviour of the system will have a great effect on the occurrence of overshoot. We made it clear that overshoot in general cannot occur in systems which will finally approach an equilibrium steady state. We validated this result by showing that a two-state system, an equilibrium three-state system, and an equilibrium MWC model will never perform overshoot. We further classified overshoot into two types, named as simple overshoot and oscillating overshoot. We found that oscillating overshoot in general will occur in systems far from equilibrium. All the above results clearly show that overshoot is a non-equilibrium dynamic phenomenon and thus sustained energy consumption is required for the system to perform this important biological function. Further analysis is expected for deeper insights into the physical mechanisms of overshoot in biological systems.
5 5.1
Methods Overshoot in the equilibrium MWC model
In this section, we shall prove that an equilibrium MWC model cannot perform overshoot. Readers who are 143
& The Institution of Engineering and Technology 2014
www.ietdl.org not interested in the mathematical derivation can skip this part. We note that when the input level is elevated from Ii to If, the system is driven from one steady state to another steady state. Let Qi and Qf = (qij )2n×2n be the transition rate matrices of the MWC model under input levels Ii and If, respectively. Let π = (π1, …, π2n) denote the initial distribution of the system and let μ = (μ1, …, μ2n) denote the final distribution of the system. Then π and μ are, respectively, the steady-state distributions of the system under input levels Ii and If. We know from (13) that the eigenvalues λ1, … λ2n of the transition rate matrix Qf of an equilibrium MWC model satisfy λ1 = 0 and λ2,⋯ , λ2n < 0. Let M = diag(μ1, …, μ2n) be a diagonal matrix whose diagonal elements are μ1, …, μ2n, respectively. The detailed balance condition μiqij = μjqji indicates that the matrix √
1 1 mi qij − (28) S = M 2 Qf M 2 = √
mj
This relation and the orthogonality relation (29) suggest that for m = 2, …, 2n,
n √
am = c mi rmi
n 2n Ii √
√
mk rmk + mk rmk I f k=1 i=1 k=n+1 2 n Ii √
−1 mi rmi ≤ 0 =c If i=1
(34) Since λ1 = 0 and αm ≤ 0 for m = 2, …, 2n, the output f (t) of the system is a weighted sum of 2n − 1 exponential functions with negative powers and non-positive weights. This shows that f (t) must be a monotonic function. Thus an equilibrium MWC model can never perform overshoot. 5.2 Estimation of the transition rates of the breast cancer cell system
2n×2n
is a symmetric matrix whose eigenvalues are λ1, …, λ2n. According to the theory of linear algebra, there exists an orthogonal matrix R = (rij )2n×2n, such that RSR T is a diagonal matrix D = diag(λ1, …, λ2n). This shows that the 2n rows of the orthogonal matrix R are, respectively, the 2n unit eigenvectors of the symmetric matrix S. We easily see that the unit eigenvector of S corresponding to the √
√
eigenvalue λ1 = 0 is m1 , . . . , m2n . Thus the orthogonality of the matrix R leads to 2n √
mk rmk = 0,
m = 2, . . . , 2n
(29)
k=1
According to previous discussions, the probability pi(t) of state i at time t can be calculated as 1 1 pi (t) = petQf i = pM −2 R−1 etD RM 2
i
2n 2n √
pk −l t = mi rmi √
rmk e m m k m=1 k=1
n
pi (t) =
i=1
2n
am e−lm t
(31)
m=1
n √
am = mi rmi i=1
2n pk √
rmk mk k=1
(32)
Note that π and μ are, respectively, the steady-state distributions of the system under input levels Ii and If. According to the detailed balance condition, we easily see that there exists a constant c > 0, such that I pk = c i mk , If
pk = cmk ,
k = 1, . . . , n (33) k = n + 1, . . . , 2n
144 & The Institution of Engineering and Technology 2014
1 Qn n=0
where
P = exp (Q) =
(30)
Thus the output f (t) of the equilibrium MWC model is
f(t) =
In this section, we shall use the data set of the SUM159 human breast cancer cell line from recent published work [13] to estimate the transition rates between any pair of differentiation states of the breast cancer cell system. The data set includes the time-course measurements of the differentiation states of all cells in the breast cancer cell system per cell division. According to the time-course data collected after in vitro culture of 6 days, the transition probabilities per cell division between any pair of differentiation states can be estimated as p SS = 0.58, pSB = 0.35, pSL = 0.07, pBS = 0.01, pBB = 0.99, pBL = 0.00, pLS = 0.04, pLB = 0.49 and pLL = 0.47. We denote by P the transition probability matrix per cell division whose component pij represents the transition probability from state i to state j per cell division. Mathematically, the transition rate matrix Q of the breast cancer cell system and the transition probability matrix P per cell division is related by
n!
(35)
The remaining question is to estimate Q from P, which is a complicated mathematical problem called the embedding problem of Markov chains. Readers who are interested in how to estimate Q from P may refer to Metzner et al. [38]. Using the method provided in [38], we can estimate the transition rates between any pair of differentiation states of the breast cancer cell system, as illustrated in Fig. 2d. According to the transition rates, the steady-state probabilities of the three differentiation states can be estimated as μS = 0.023, μB = 0.971 and μL = 0.005. Combining the steady-state probabilities and the transition rates, the net flux of the breast cancer cell system can be estimated as J = μSqSB − μBqBS = μBqBL − μLqLB = μLqLS − μSqSL = −0.0028.
6
Acknowledgment
The authors gratefully acknowledge financial support from the NSFC 11271029 and the NSFC 11171024. The first author also acknowledges financial support from the Academic Award for Young PhD Researchers granted by the Ministry of Education of China. IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org 7
References
1 McAdams, H.H., Arkin, A.: ‘Stochastic mechanisms in gene expression’, Proc. Natl. Acad. Sci., 1997, 94, (3), pp. 814–819 2 Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S.: ‘Stochastic gene expression in a single cell science signalling’, Science, 2002, 297, (5584), pp. 1183–1186 3 Ozbudak, E.M., Thattai, M., Kurtser, I., Grossman, A.D., van Oudenaarden, A.: ‘Regulation of noise in the expression of a single gene’, Nature Genet., 2002, 31, (1), pp. 69–73 4 Paulsson, J.: ‘Summing up the noise in gene networks’, Nature, 2004, 427, (6973), pp. 415–418 5 Raser, J.M., O’Shea, E.K.: ‘Noise in gene expression: origins, consequences and control’, Science, 2005, 309, (5743), pp. 2010–2013 6 Kærn, M., Elston, T.C., Blake, W.J., Collins, J.J.: ‘Stochasticity in gene expression: from theories to phenotypes’, Nature Rev. Genet., 2005, 6, (6), pp. 451–464 7 Cai, L., Friedman, N., Xie, X.S.: ‘Stochastic protein expression in individual cells at the single molecule level’, Nature, 2006, 440, (7082), pp. 358–362 8 Xie, X.S., Choi, P.J., Li, G.-W., Lee, N.K., Lia, G.: ‘Single-molecule approach to molecular biology in living bacterial cells’, Annu. Rev. Biophys., 2008, 37, pp. 417–444 9 Cornish-Bowden, A.: ‘Fundamentals of Enzyme Kinetics’ (Wiley-Blackwell, 2012, 4th edn.) 10 Jia, C., Liu, X.-F., Qian, M.-P., Jiang, D.-Q., Zhang, Y.-P.: ‘Kinetic behavior of the general modifier mechanism of Botts and Morales with non-equilibrium binding’, J. Theory Biol., 2012, 296, pp. 13–20 11 Beard, D.A., Qian, H.: ‘Chemical Biophysics: Quantitative Analysis of Cellular Systems’ (Cambridge University Press, 2008, 1st edn.) 12 Keener, J.P., Sneyd, J.: ‘Mathematical Physiology’ (Springer, 1998, 1st edn.) 13 Gupta, P.B., Fillmore, C.M., Jiang, G., et al.: ‘Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells’, Cell, 2011, 146, (4), pp. 633–644 14 Segel, L.A., Goldbeter, A., Devreotes, P.N., Knox, B.E.: ‘A mechanism for exact sensory adaptation based on receptor modification’, J. Theory Biol., 1986, 120, (2), pp. 151–179 15 Knox, B.E., Devreotes, P.N., Goldbeter, A., Segel, L.A.: ‘A molecular mechanism for sensory adaptation based on ligand-induced receptor modification’, Proc. Natl. Acad. Sci., 1986, 83, (8), pp. 2345–2349 16 Behar, M., Hao, N., Dohlman, H.G., Elston, T.C.: ‘Mathematical and computational analysis of adaptation via feedback inhibition in signal transduction pathways’, Biophys. J., 2007, 93, (3), pp. 806–821 17 François, P., Siggia, E.D.: ‘A case study of evolutionary computation of biochemical adaptation’, Phys. Biol., 2008, 5, (2), pp. 026009 18 Ma, W., Trusina, A., El-Samad, H., Lim, W.A., Tang, C.: ‘Defining network topologies that can achieve biochemical adaptation’, Cell, 2009, 138, (4), pp. 760–773 19 Lan, G., Sartori, P., Neumann, S., Sourjik, V., Tu, Y.: ‘The energy-speed-accuracy trade-off in sensory adaptation’, Nature Phys., 2012, 8, (5), pp. 422–428
IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
20 Tu, Y., Shimizu, T.S., Berg, H.C.: ‘Modeling the chemotactic response of Escherichia coli to time-varying stimuli’, Proc. Natl. Acad. Sci., 2008, 105, (39), pp. 14855–14860 21 Tu, Y.: ‘Quantitative modeling of bacterial chemotaxis: signal amplification and accurate adaptation’, Annu. Rev. Biophys., 2013, 42, pp. 337–359 22 Mettetal, J.T., Muzzey, D., Gómez-Uribe, C., van Oudenaarden, A.: ‘The frequency dependence of osmo-adaptation in Saccharomyces cerevisiae’, Science, 2008, 319, (5862), pp. 482–484 23 Marchant, J.S., Taylor, C.W.: ‘Rapid activation and partial inactivation of inositol trisphosphate receptors by inositol trisphosphate’, Biochemistry, 1998, 37, (33), pp. 11524–11533 24 Adkins, C.E., Taylor, C.W.: ‘Lateral inhibition of inositol 1,4,5-trisphosphate receptors by cytosolic Ca2+’, Curr. Biol., 1999, 9, (19), pp. 1115–1118 25 Gyorke, S., Fill, M.: ‘Ryanodine receptor adaptation: control mechanism of Ca2+-induced Ca2+ release in heart’, Science, 1993, 260, (5109), pp. 807–809 26 Keizer, J., Levine, L.: ‘Ryanodine receptor adaptation and Ca2+-induced Ca2+ release-dependent Ca2+ oscillations’, Biophys. J., 1996, 71, (6), pp. 3477–3487 27 Fill, M., Zahradníková, A., Villalba-Galea, C.A., Zahradník, I., Escobar, A.L., Györke, S.: ‘Ryanodine receptor adaptation’, J. General Physiol., 2000, 116, (6), pp. 873–882 28 Li, Y., Goldbeter, A.: ‘Frequency specificity in intercellular communication. Influence of patterns of periodic signaling on target cell responsiveness’, Biophys. J., 1989, 55, (1), pp. 125–145 29 Ogata, K.: ‘Discrete-Time Control Systems’ (Prentice-Hall, 1995, 2nd edn.) 30 Golnaraghi, F., Kuo, B.C.: ‘Automatic Control Systems’ (Wiley, 2009, 9th edn.) 31 Allen, P.E., Holberg, D.R.: ‘Cmos Analog Circuit Design’ (Oxford University Press, USA, 2011, 3rd edn.) 32 Berman, A., Plemmons, R.J.: ‘Nonnegative Matrices in the Mathematical Sciences’ (Society for Industrial and Applied Mathematics, 1987, 1st edn.) 33 Zhou, D., Wu, D., Li, Z., Qian, M., Zhang, M.Q.: ‘Population dynamics of cancer cells with cell-state conversions between cancer stem cells and non-stem cancer cells’, Quant. Biol., 2013, 1, (3), pp. 201–208 34 Jiang, D.-Q., Qian, M., Qian, M.-P.: ‘Mathematical Theory of Nonequilibrium Steady States: On the Frontier of Probability and Dynamical Systems’ (Springer, 2004, 1st edn.) 35 Zhang, X.-J., Qian, H., Qian, M.: ‘Stochastic theory of nonequilibrium steady states and its applications. Part I’, Phys. Rep., 2012, 510, (1), pp. 1–86 36 Ge, H., Qian, M., Qian, H.: ‘Stochastic theory of nonequilibrium steady states. Part II: applications in chemical biophysics’, Phys. Rep., 2012, 510, (3), pp. 87–118 37 Qian, H.: ‘The mathematical theory of molecular motor movement and chemomechanical energy transduction’, J. Math. Chem., 2000, 27, (3), pp. 219–234 38 Metzner, P., Dittmer, E., Jahnke, T., Schütte, C.: ‘Generator estimation of Markov jump processes’, J. Comput. Phys., 2007, 227, (1), pp. 353–375
145
& The Institution of Engineering and Technology 2014
Special Issue on Selected papers from The 7th IEEE International Conference on Systems Biology (ISB 2013) ISSN 1751-8849
Overshoot in biological systems modelled by Markov chains: a non-equilibrium dynamic phenomenon Chen Jia1,2, Minping Qian1, Daquan Jiang1,3 1
LMAM, School of Mathematical Sciences, Peking University, Beijing 100871, People’s Republic of China Beijing International Center for Mathematical Research, Beijing 100871, People’s Republic of China 3 Center for Statistical Science, Peking University, Beijing 100871, People’s Republic of China E-mail: [email protected] 2
Abstract: A number of biological systems can be modelled by Markov chains. Recently, there has been an increasing concern about when biological systems modelled by Markov chains will perform a dynamic phenomenon called overshoot. In this study, the authors found that the steady-state behaviour of the system will have a great effect on the occurrence of overshoot. They showed that overshoot in general cannot occur in systems that will finally approach an equilibrium steady state. They further classified overshoot into two types, named as simple overshoot and oscillating overshoot. They showed that except for extreme cases, oscillating overshoot will occur if the system is far from equilibrium. All these results clearly show that overshoot is a non-equilibrium dynamic phenomenon with energy consumption. In addition, the main result in this study is validated with real experimental data.
1
Introduction
Recent advances in the single-cell and single-molecule experiments have shown that biological systems in living cells are inherently stochastic [1–8]. It is widely observed that a number of biological systems can transition stochastically among multiple states. These systems are often mathematically characterised by the stochastic model of Markov chains, or in the language of physics, master equations. Typical examples of Markov chain systems include enzyme kinetics with the Michaelis-Menten mechanism [9] or the general modifier mechanism of Botts and Morales [10], phosphorylation–dephosphorylation kinetics of proteins [11], conformational changes of receptors [12] and stochastic state transitions of cell populations [13]. In recent years, scientists have become increasingly concerned about a kind of dynamic phenomenon in biological systems which is known as overshoot or biochemical adaptation [14–19]. Overshoot of a biological system refers to the dynamic phenomenon that the output, which is usually the average of an observable, of the system varies with time in a non-monotonic way and exceeds both its initial and steady-state values over a certain period of time. An intuitive description of overshoot is depicted in Fig. 1, where the output of the system first rises to a peak and then declines to a lower plateau. Overshoot is widely observed in numerous biological systems. Typical examples of overshoot include chemotaxis of bacteria [20, 21], osmotic sensing in yeast [22], calcium dose response of the inositol trisphosphate receptors [23, 24] and the ryanodine receptors [25–27], and hormone dose response of the hormone receptors [28]. Overshoot is an 138 & The Institution of Engineering and Technology 2014
important biological function possessed by many living systems. It allows the system to detect environmental changes more accurately, enables the system to respond to environmental fluctuations more rapidly and protects the system from irreversible damages caused by unfavourable conditions. To better understand how overshoot is achieved in biochemical feedback networks, several research groups studied the relationship between overshoot and the network topology [16–18]. Tang and co-workers [18] searched all possible three-node network topologies and found that overshoot is most likely to occur in two types of networks: the negative feedback loop and the incoherent feedforward loop. Tu and co-workers [19] studied the stochastic dynamics of the negative feedback loop in detail and found that the negative feedback mechanism breaks detailed balance, and thus always operates out of equilibrium with energy dissipation. These two works give us a hint that overshoot may tend to occur in biological systems which are far from equilibrium. To avoid misinterpretation, we point out here that the word ‘equilibrium’ appearing in this paper is referred to as the concept of equilibrium steady state in statistical physics, where the steady state of a system is called equilibrium (non-equilibrium) if the detailed balance condition is satisfied (broken). In fact, the concept of overshoot has long been suggested and studied in control theory [29, 30], electronics [31] and signal processing. In these disciplines, overshoot was studied in various deterministic systems composed of several ordinary differential equations and some overshoot conditions were provided. However, the traditional study of overshoot under the framework of deterministic systems does not capture the physical essence of overshoot. Many IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org
Fig. 1 Overshoot in biological systems Overshoot refers to the dynamic phenomenon that the output of the system varies with time in a non-monotonic way and exceeds both its initial and steady-state values over a certain period of time
important topics about overshoot, such as its relation to the breakdown of detailed balance and to energy dissipation, can only be clearly seen in stochastic systems. Up till now, there is still a lack of a general analysis of overshoot, which captures its physical essence, in biological systems which are inherently stochastic. In this paper, we presented a general analysis of overshoot in biological systems which can be modelled by Markov chains with two, three or multiple states. We found that the steady-state behaviour of the system will have a great effect on the occurrence of overshoot. We made it clear that overshoot in general cannot occur in systems which will finally approach an equilibrium steady state. This explains why overshoot can only be observed in systems with three or more states and cannot be observed in systems with only two states. We further classified overshoot into two types, named as simple overshoot and oscillating overshoot. We found that except for extreme cases, oscillating overshoot will occur in systems far from equilibrium. All the above results clearly show that overshoot is a non-equilibrium dynamic phenomenon and thus sustained energy consumption is required for the system to perform this important biological function. In addition, we used the experimental data of SUM159 human breast cancer cell line to validate the main theoretical result in this paper.
2
where pi(t) is the probability of the system being in state i at time t. We denote by πi the initial probability of state i, that is, the probability of the system being in state i at time t = 0 and denote by μi the steady-state probability of state i, that is, the probability of the system being in state i when time t is sufficiently large. Then f(0) = Ni=1 pi fi is the initial output of the system and f(1) = Ni=1 mi fi is the steady-state output of the system. According to Fig. 1, the system performs overshoot if and only if there exists a time t such that f (t) is larger than both f(0) and f(∞). Next, we shall use three specific examples to help the readers understand the Markov chain model discussed above. Example 1: We consider the well-known Michaelis-Menten enzyme kinetics k−1
k2
E + S O ES E + P k1
(2)
where E is an enzyme involved in converting the substrate S into the product P. If there is only one enzyme molecule, it
Model
In this paper, we consider biological systems that can be mathematically modelled as continuous-time Markov chains with multiple states. The characteristic property of the Markov chain is that it retains no memory of where it has been in the past. This means that where the system will go next only depends on its current state, but not depends on its prior states. We assume that the Markov chain system can be found in N states, 1, 2, …, N, and thus the system can transition stochastically among these states. The specific meaning of these states varies from systems to systems. Each state can represent a binding state of an enzyme molecule [10], a conformational state of a receptor molecule [14, 15, 25–27], a cellular state of a cell population [13] and so on. We further assume that the system has an observable f. If the system is in state i, the observation of the system will be fi. The output f (t) of the system at time t is then the weighted average of the observations of all states
f(t) =
N
pi (t)fi
i=1
IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
(1)
Fig. 2 Markov chain systems with two, three or multiple states a Markov chain systems with two states b Markov chain systems with three states c–e Examples of biological systems that can be modelled by Markov chains c The catalytic cycle of the Michaelis-Menten enzyme mechanism, where E is the enzyme, S is the substrate, ES is the enzyme-substrate complex and [S] is the concentration of the substrate d Cell-state dynamics of human breast cancer cells An individual breast cancer cell can transition stochastically among three differentiation states: stem-like (S), basal (B) and luminal (L) states The transition rates between states are estimated from the data set of the SUM 159 human breast cancer cell line [13] The transition rates are shown per cell division e The MWC allosteric model which describes the conformational changes of receptors with n − 1 identical subunits Each receptor can bind to an attractant (small solid circle) If the attractant binding site is occupied, all subunits will switch together from the inactivated configuration (solid square) to the activated configuration (large solid circle) In addition, each subunit can bind to a ligand (solid ellipse) According to whether all subunits are activated or inactivated and the number of subunits which have bound to the ligand, each receptor may convert among 2n possible states The upper n states are activated states and the lower n states are inactivated states 139
& The Institution of Engineering and Technology 2014
www.ietdl.org may convert between two states: the free enzyme E and the enzyme-substrate complex ES. Then from the perspective of a single enzyme molecule, the Michaelis-Menten kinetics can be represented by the catalytic cycle illustrated in Fig. 2c. We denote by E0 = [E] + [ES] the total enzyme concentration. Then pE(t) = [E]/E0 and pES(t) = [ES]/E0 represent the probability of a single enzyme molecule being in state E and state ES at time t, respectively. Mathematically, the catalytic cycle illustrated in Fig. 2c is nothing but a Markov chain with two states. According to the law of mass action, the probability flux from state E to state ES at time t is k1[S]pE(t), and the probability flux from state ES to state E at time t is (k−1 + k2)pES(t), where we add k−1 and k2 since there are two ways of transition from state ES to state E. Based on the expressions of the probability fluxes, we easily see that the transition rate from state E to state ES is k1[S] and the transition rate from state ES to state E is k−1 + k2. In the Michaelis-Menten enzyme system, the output f (t) is often chosen as the instantaneous rate of the product formation
f(t) =
d[P] = k2 [ES] = k2 E0 pES (t) dt
(3)
where the second equality is because of the law of mass action. Example 2: Recent research shows that human breast cancer cells within individual tumours can exist in three differentiation states that differ in functional attributes: stem-like (S), basal (B) and luminal (L) states. An individual breast cancer cell can transition stochastically among these three states. Mathematically, Lander and co-workers [13] modelled the cell-state transitions and dynamics in the breast cancer cell population as a three-state Markov chain depicted in Fig. 2d. In the breast cancer cell system, the output f (t) is often chosen as the proportion of stem-like, basal or luminal cells
f(t) = pi (t),
i = S, B, L
(4)
Example 3: Structural studies show that receptors in living cells are often protein complexes with multiple subunits. Mathematically, Monod, Wyman and Changeux modelled the allosteric protein interactions in receptors as a Markov chain with multiple states. The Monod–Wyman–Changeux (MWC) allosteric model assumes that each receptor consists of n − 1 identical subunits, each of which can switch between two configurations, an activated one and an inactivated one. The MWC model further assumes that each receptor has an attractant binding site. Once an attractant binds to the receptor, all subunits will switch together from being inactivated to being activated. In addition, the MWC model assumes that each subunit has a ligand binding site. A ligand can bind to a subunit in either configuration, but the dissociation constants are different. According to whether all subunits are activated or inactivated and the number of subunits which have bound to the ligand, the conformational changes of each receptor can be modelled as a Markov chain with 2n states. The transition diagram of the MWC model when n = 5 is depicted in Fig. 2e, where each solid square represents an inactivated subunit and each large (solid) circle represents an activated subunit. Thus the upper n states in Fig. 2e are activated states and the lower n states are inactivated states. In the MWC model, the output 140 & The Institution of Engineering and Technology 2014
f(t) is often chosen as the average activity, that is, the sum of probabilities of those activated states f(t) =
n
pi (t)
(5)
i=1
3
Results
3.1 General analysis and overshoot in two-state systems We now use the theory of Markov chains to present a general analysis of overshoot. We know that the dynamics of the probability distribution p(t) = ( p1(t), …, pN(t)) of the Markov chain system is governed by the master equation of the matrix form dp(t) = p(t)Q dt
(6)
where Q = (qij )N×N is the transition rate matrix of the Markov chain system and qij is the transition rate from state i to state j. We assume without loss of generality that the transition rate matrix Q has N linear independent eigenvectors since any matrix can be approximated by a matrix satisfying this condition with arbitrarily high accuracy. Under this assumption, the output f (t) of the system is nothing but a linear combination of exponential functions
f(t) = p(t)f T = p(0)etQ f T =
N
ci eli t
(7)
i=1
where f = ( f1,…, fN) is a vector whose components are the observations of all states, λ1, …, λN are all eigenvalues of the transition rate matrix Q, and c1, …, cN are N constants. The well-known Perron–Frobenius theorem [32] in matrix theory claims that one of the eigenvalues of the transition rate matrix Q must be 0 and the real parts of other eigenvalues are all negative. In the following discussion, we always assume that λ1 = 0. Thus the output f (t) of the system can be rewritten that
f(t) = c1 +
N
ci eli t
(8)
i=2
We easily see that ci is real if λi is real, and ci and cj are conjugate complex numbers if λi and λj are conjugate complex numbers. An interesting phenomenon widely observed in experiments is that overshoot cannot be observed in Markov chain systems with only two states and can be observed in Markov chain systems with three or more states [33]. We now use the previous analysis to explain this phenomenon. According to (8), the output of a two-state system depicted in Fig. 2a is given by
f(t) = c1 + c2 el2 t
(9)
which is obviously a monotonic function since λ2 must be a real number. This suggests that a Markov chain system with only two states will never perform overshoot. IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org 3.2
Equilibrium and non-equilibrium systems
We have seen from previous discussions that Markov chain systems with only two states will never perform overshoot and overshoot can only occur in systems with three or more states. This raises a natural question of what is the essential difference between systems with two states and systems with three or more states. We noted that there is an apparent difference between systems with two states and systems with three or more states. A two-state system is always an equilibrium system which will finally approach an equilibrium steady state, whereas a system with three or more states may be a non-equilibrium system which will finally approach a non-equilibrium steady state [34–36]. Here, we have used the concepts of equilibrium and non-equilibrium steady states in non-equilibrium statistical physics, where the steady state of a system is called equilibrium if each pair of states i and j of the system satisfies the detailed balance condition
mi qij = mj q ji
(10)
(otherwise, the steady state of the system is called nonequilibrium) where μi is the steady-state probability of state i and qij is the transition rate from state i to state j. In the following discussion, a system which will finally approach an equilibrium (non-equilibrium) steady state is referred to as an equilibrium (non-equilibrium) system. From the point of view of statistical mechanics, an equilibrium system in the steady state is microscopic reversible and does not consume energy, whereas a non-equilibrium system in the steady state is microscopic irreversible and always consumes energy. Mathematically, whether a system is a non-equilibrium system or not is linked to the net flux the system [34]. To make the readers understand the concept of the net flux, we now limit our discussion to three-state systems depicted in Fig. 2b. In a three-state system, we denote by c the cycle 1 → 2 → 3 → 1 and denote by c− the reverse cycle 1 → 3 → 2 → 1. Let wc(t) denote the number of cycle c formed by the system up to time t. Then the flux wc of cycle c is defined to be wc (t) t1 t
wc = lim
(11)
We clearly see that the flux wc of cycle c represents the number of the forming of cycle c per unit time. Similarly, we can define the flux wc− for the reverse cycle c−. The net flux J of a three-state system is then defined to be J = wc − wc−, which represents the net number of the forming of cycle c per unit time. Interestingly, the net flux J of a three-state system can be represented by the steady-state probabilities and the transition rates as J = m1 q12 − m2 q21 = m2 q23 − m3 q32 = m3 q31 − m1 q13
(12)
This relation is a special case of the well-known circulation
a1 = −
decomposition theorem [34] in the Markov chain theory. From (10) and (12), we easily see that a three-state system is a non-equilibrium system if and only if the net flux J fails to be zero. This shows that the net flux J characterises how far the system is away from equilibrium. Generally speaking, equilibrium systems in the steady state do not consume energy, whereas non-equilibrium systems always consume energy to maintain non-zero net fluxes. Finally, we state a simple but important fact about equilibrium systems which will be frequently used in following discussions. In an equilibrium system, one eigenvalue of the transition rate matrix Q must be zero and all other eigenvalues must be negative real numbers. In other words, the eigenvalues λ1, …, λN of the transition rate matrix Q of an equilibrium system must satisfy
l1 = 0, l2 , . . . , lN , 0 3.3
Simple overshoot
We have seen from previous discussions that systems with only two states will always approach an equilibrium steady state and will never perform overshoot. However, the situation is totally different in systems with three states. According to (8), the output f (t) of a three-state system has the general form of
f(t) = c1 + c2 el2 t + c3 el3 t
(14)
If λ2 and λ3 are negative real numbers and the signs of c2 and c3 are not the same, then f(t) may not be a monotonic function and the system may perform overshoot. This type of overshoot is named as simple overshoot (Fig. 3a). To gain a deeper insight into simple overshoot, we consider three-state systems depicted in Fig. 2b starting from state 1 whose transition rate matrix has three real eigenvalues satisfying λ1 = 0 and λ2, λ3 < 0. To simplify notations, we make an additional assumption that the transition rates of the system satisfy q21 = q31. We impose this condition onto the system for two reasons. Firstly, this assumption ensures that the transition rate matrix has three real eigenvalues. Secondly, this assumption reduces the complexities of calculations and formulations to a remarkable extent. Under this assumption, the three eigenvalues of the transition rate matrix have explicit expressions of λ1 = 0, λ2 = −(q12 + q13 + q21) and λ3 = −(q23 + q32 + q21). Furthermore, the output f(t) of the system can be explicitly calculated as
f(t) = a1 + a2 el2 t + J a3 el3 t
(15)
where J is the net flux of the system, and α1, α2, α3 are three constants with the following expressions (see (16))
a2 =
(q12 + q13 )[(q12 − q32 )(f1 − f2 ) + (q13 − q23 )(f1 − f3 )] l2 (l2 − l3 ) (17)
a3 =
l2 (f2 − f3 ) q21 (l2 − l3 )
(18)
We are now in the position to consider when overshoot will
q21 l3 f1 + (q12 q31 + q12 q32 + q32 q13 )f2 + (q12 q23 + q13 q23 + q21 q31 )f3 l2 l3
IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
(13)
(16)
141
& The Institution of Engineering and Technology 2014
www.ietdl.org
Fig. 3 Overshoot in Markov chain systems with three states a Simple overshoot and oscillating overshoot in three-state systems If λ2 and λ3 are negative real numbers, then the system may perform simple overshoot If λ2 and λ3 are conjugate complex numbers, then the system will perform oscillating overshoot b Overshoot in the SUM 159 human breast cancer cell line The proportion of luminal cells experiences a transient increase after isolation of stem-like cells followed by a slow decrease to its steady-state value
occur in an equilibrium three-state system. We recall that a three-state system is an equilibrium system if and only if the net flux J is zero. According to (15), the output of an equilibrium three-state system can be simplified as
f(t) = c1 + c2 el2 t
(19)
which is obviously a monotonic function since λ2 is a real number. This shows that under mild conditions, an equilibrium system with three states will never perform overshoot. 3.4
conjugate complex eigenvalues when the net flux J is sufficiently large. To this end, let j1 = μ1q12, j2 = μ2q23 and j3 = μ3q31 be the clockwise probability fluxes of the system. According to (12), the transition rate matrix Q can be represented by μ1, μ2, μ3, j1, j2, j3 and J as ⎛ ⎜ ⎜ ⎜ Q=⎜ ⎜ ⎜ ⎝
Oscillating overshoot
We have seen from previous discussions that a three-state system may perform simple overshoot if λ2 and λ3 are negative real numbers. However, more interesting is the case when λ2 and λ3 are conjugate complex numbers. If we denote them by λ2 = λ + ωi and λ3 = λ − ωi, then the output f(t) of the system can be rewritten as
f(t) = c1 + 2|c2 |elt cos (vt + w)
(20)
where λ is a negative real number and j is the argument of c2. This expression clearly shows that after the output rises to a peak and returns toward its initial value, it may rise and decline again to form a damped oscillation. This type of overshoot is named as oscillating overshoot (Fig. 3a). Next, we shall present a detailed discussion of oscillating overshoot in equilibrium and non-equilibrium three-state systems. According to (13), the eigenvalues of the transition rate matrix of an equilibrium system are all real numbers. This shows that λ2 and λ3 cannot be conjugate complex numbers and thus oscillating overshoot will never occur in equilibrium three-state systems. The situation is totally different if we consider three-state systems which will finally approach a steady state far from equilibrium. In other words, we consider a three-state system with a sufficiently large net flux J, since the net flux characterises how far the system is away from equilibrium. We have seen that a three-state system will perform oscillating overshoot if the transition rate matrix has a pair of conjugate complex eigenvalues. What we need to do next is to study when the transition rate matrix has a pair of 142 & The Institution of Engineering and Technology 2014
−
j3 + j1 + J m1 j1 + J m2 j3 m3
j1 m1 j + j2 + J − 1 m2 j2 + J m3
⎞ j3 + J ⎟ m1 ⎟ ⎟ j2 ⎟ ⎟ m2 ⎟ j2 + j3 + J ⎠ − m3 (21)
Let pQ(λ) = det(λI − Q) be the characteristic polynomial of the transition rate matrix Q. Straightforward calculations show that p Q (l ) =
1 l(m1 m2 m3 l2 + Al + B) m1 m2 m3
(22)
where A = (m1 m2 + m2 m3 + m3 m1 )J + (m1 m2 (j2 + j3 ) + m2 m3 (j3 + j1 ) + m3 m1 (j1 + j2 ))
(23)
B = J 2 + (j1 + j2 + j3 )J + (j1 j2 + j2 j3 + j3 j1 )
(24)
and
We know from elementary algebra that the transition rate matrix Q has a pair of conjugate complex eigenvalues if and only if the discriminant Δ = A 2 − 4μ1μ2μ3B < 0. We easily calculate that D = dJ 2 + aJ + b
(25)
d = (m1 m2 + m2 m3 + m3 m1 )2 − 4m1 m2 m3
(26)
where
From (25), we make a crucial observation that when the net flux J is sufficiently large, the transition rate matrix Q has a IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org pair of conjugate complex eigenvalues if and only if δ is negative. The remaining question is to answer when δ is negative. We can prove that δ is a quantity satisfying −
1 1 ≤d≤ 27 16
(27)
We can further prove that δ attains its minimum of δmin = −1/ 27 when μ1 = μ2 = μ3 = 1/3, and attains its maximum of δmax = 1/16 when μ1 = μ2 = 1/2 and μ3 = 0. This clearly shows that the more uniform the three steady-state probabilities, μ1, μ2 and μ3, the smaller the value of δ. Thus except for the extreme case that the three steady-state probabilities are extremely scattered, δ is always negative and the system will perform oscillating overshoot when the net flux J is sufficiently large. This suggests that oscillating overshoot in general will occur in systems far from equilibrium. 3.5
Overshoot in systems with multiple states
In previous discussions, we are mainly concerned about when overshoot will occur in systems with two or three states. We have seen that a two-state system and an equilibrium three-state system in general cannot perform overshoot. This shows that overshoot is a non-equilibrium dynamic phenomenon in systems with two or three states. This raises a natural question of whether we can obtain the same conclusion in systems with multiple states. We point out that if a Markov chain system has more than three states, the mathematical complexity will become so high that it is almost impossible to present a complete discussion of overshoot. Thus in this section, we only focus on an important class of multiple-state systems, namely, the MWC allosteric model depicted in Fig. 2e, which is widely used in modelling the conformational changes of receptors in living cells. We see from Fig. 2e that some transition rates of the MWC model are regulated by the attractant concentration I, which is the input of the system. Let Ii and If be two input levels with 0 ≤ Ii < If. In experiments, we are always concerned about whether the system will perform overshoot in response to a step increase of the input from Ii to If. Interestingly, we can prove that if the MWC model is an equilibrium system, then it will never perform overshoot no matter what Ii and If are chosen. For a proof of this result, please see Section 5. This result clearly shows that an equilibrium system with multiple states in general cannot perform overshoot. All the results in the above sections indicate that overshoot, whether in systems with two, three or multiple states, is a non-equilibrium dynamic phenomenon, and thus sustained energy consumption is required for living systems to perform this important biological function. 3.6 Validation of main results with experimental data As mentioned in previous discussions, human breast cancer cells within individual tumours can transition stochastically among three differentiation states: stem-like (S), basal (B) and luminal (L) states. Mathematically, the cell-state transitions and dynamics in the breast cancer cell population can be modelled as a three-state Markov chain depicted in Fig. 2d. To validate our theoretical results with real experimental data, we shall use the data set of the SUM159 human breast cancer cell line from recent published work [13] to study overshoot in the breast cancer cell system depicted in Fig. 2d. IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
From the experimental data, we estimate the transition rates between any pair of differentiation states of the breast cancer cell system, as illustrated in Fig. 2d. For more details on the data usage and calculation, please see Section 5. The output f(t) of the breast cancer cell system is chosen to be the proportion of stem-like, basal, or luminal cells: f(t) = pi(t), i = S, B, L. Lander and co-workers found that if the stem-like cells are isolated at a particular time, then the breast cancer cell system will perform overshoot [13]. To see this, we depict the time course of the proportions of stem-like, basal and luminal cells in Fig. 3b, from which we see that the proportion of the luminal cells experiences a transient increase after isolation of the stem-like cells followed by a slow decrease to its steady-state value. In order to see whether overshoot in the breast cancer cell system is a non-equilibrium dynamic phenomenon, we need to calculate the net flux of the system. From the experimental data, the net flux of the breast cancer cell system is estimated as J = −0.0028 (see Section 5), which fails to be zero. This shows that overshoot in the SUM159 human breast cancer cell line is indeed a non-equilibrium dynamic phenomenon, which is consistent with the main theoretical result in this paper.
4
Discussion
Living systems are highly dissipative, exchanging materials and energy with their environments and consuming energy to carry out various important biological functions. If the exchange with the environment is sustained, then the living system always approaches a steady state far from equilibrium. Recent studies show that many important biological phenomena, such as coherence resonance in excitable systems [35], unidirectional movement of molecular motors [37] and switching behaviour of the general modifier mechanism of Botts and Morales [10], fail to occur in systems which will finally approach an equilibrium steady state. This suggests that sustained energy consumption is required for living systems to perform these important biological functions. In this paper, we presented a general discussion of the dynamic phenomenon of overshoot in biological systems modelled by Markov chains (master equations). We found that the steady-state behaviour of the system will have a great effect on the occurrence of overshoot. We made it clear that overshoot in general cannot occur in systems which will finally approach an equilibrium steady state. We validated this result by showing that a two-state system, an equilibrium three-state system, and an equilibrium MWC model will never perform overshoot. We further classified overshoot into two types, named as simple overshoot and oscillating overshoot. We found that oscillating overshoot in general will occur in systems far from equilibrium. All the above results clearly show that overshoot is a non-equilibrium dynamic phenomenon and thus sustained energy consumption is required for the system to perform this important biological function. Further analysis is expected for deeper insights into the physical mechanisms of overshoot in biological systems.
5 5.1
Methods Overshoot in the equilibrium MWC model
In this section, we shall prove that an equilibrium MWC model cannot perform overshoot. Readers who are 143
& The Institution of Engineering and Technology 2014
www.ietdl.org not interested in the mathematical derivation can skip this part. We note that when the input level is elevated from Ii to If, the system is driven from one steady state to another steady state. Let Qi and Qf = (qij )2n×2n be the transition rate matrices of the MWC model under input levels Ii and If, respectively. Let π = (π1, …, π2n) denote the initial distribution of the system and let μ = (μ1, …, μ2n) denote the final distribution of the system. Then π and μ are, respectively, the steady-state distributions of the system under input levels Ii and If. We know from (13) that the eigenvalues λ1, … λ2n of the transition rate matrix Qf of an equilibrium MWC model satisfy λ1 = 0 and λ2,⋯ , λ2n < 0. Let M = diag(μ1, …, μ2n) be a diagonal matrix whose diagonal elements are μ1, …, μ2n, respectively. The detailed balance condition μiqij = μjqji indicates that the matrix √
1 1 mi qij − (28) S = M 2 Qf M 2 = √
mj
This relation and the orthogonality relation (29) suggest that for m = 2, …, 2n,
n √
am = c mi rmi
n 2n Ii √
√
mk rmk + mk rmk I f k=1 i=1 k=n+1 2 n Ii √
−1 mi rmi ≤ 0 =c If i=1
(34) Since λ1 = 0 and αm ≤ 0 for m = 2, …, 2n, the output f (t) of the system is a weighted sum of 2n − 1 exponential functions with negative powers and non-positive weights. This shows that f (t) must be a monotonic function. Thus an equilibrium MWC model can never perform overshoot. 5.2 Estimation of the transition rates of the breast cancer cell system
2n×2n
is a symmetric matrix whose eigenvalues are λ1, …, λ2n. According to the theory of linear algebra, there exists an orthogonal matrix R = (rij )2n×2n, such that RSR T is a diagonal matrix D = diag(λ1, …, λ2n). This shows that the 2n rows of the orthogonal matrix R are, respectively, the 2n unit eigenvectors of the symmetric matrix S. We easily see that the unit eigenvector of S corresponding to the √
√
eigenvalue λ1 = 0 is m1 , . . . , m2n . Thus the orthogonality of the matrix R leads to 2n √
mk rmk = 0,
m = 2, . . . , 2n
(29)
k=1
According to previous discussions, the probability pi(t) of state i at time t can be calculated as 1 1 pi (t) = petQf i = pM −2 R−1 etD RM 2
i
2n 2n √
pk −l t = mi rmi √
rmk e m m k m=1 k=1
n
pi (t) =
i=1
2n
am e−lm t
(31)
m=1
n √
am = mi rmi i=1
2n pk √
rmk mk k=1
(32)
Note that π and μ are, respectively, the steady-state distributions of the system under input levels Ii and If. According to the detailed balance condition, we easily see that there exists a constant c > 0, such that I pk = c i mk , If
pk = cmk ,
k = 1, . . . , n (33) k = n + 1, . . . , 2n
144 & The Institution of Engineering and Technology 2014
1 Qn n=0
where
P = exp (Q) =
(30)
Thus the output f (t) of the equilibrium MWC model is
f(t) =
In this section, we shall use the data set of the SUM159 human breast cancer cell line from recent published work [13] to estimate the transition rates between any pair of differentiation states of the breast cancer cell system. The data set includes the time-course measurements of the differentiation states of all cells in the breast cancer cell system per cell division. According to the time-course data collected after in vitro culture of 6 days, the transition probabilities per cell division between any pair of differentiation states can be estimated as p SS = 0.58, pSB = 0.35, pSL = 0.07, pBS = 0.01, pBB = 0.99, pBL = 0.00, pLS = 0.04, pLB = 0.49 and pLL = 0.47. We denote by P the transition probability matrix per cell division whose component pij represents the transition probability from state i to state j per cell division. Mathematically, the transition rate matrix Q of the breast cancer cell system and the transition probability matrix P per cell division is related by
n!
(35)
The remaining question is to estimate Q from P, which is a complicated mathematical problem called the embedding problem of Markov chains. Readers who are interested in how to estimate Q from P may refer to Metzner et al. [38]. Using the method provided in [38], we can estimate the transition rates between any pair of differentiation states of the breast cancer cell system, as illustrated in Fig. 2d. According to the transition rates, the steady-state probabilities of the three differentiation states can be estimated as μS = 0.023, μB = 0.971 and μL = 0.005. Combining the steady-state probabilities and the transition rates, the net flux of the breast cancer cell system can be estimated as J = μSqSB − μBqBS = μBqBL − μLqLB = μLqLS − μSqSL = −0.0028.
6
Acknowledgment
The authors gratefully acknowledge financial support from the NSFC 11271029 and the NSFC 11171024. The first author also acknowledges financial support from the Academic Award for Young PhD Researchers granted by the Ministry of Education of China. IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
www.ietdl.org 7
References
1 McAdams, H.H., Arkin, A.: ‘Stochastic mechanisms in gene expression’, Proc. Natl. Acad. Sci., 1997, 94, (3), pp. 814–819 2 Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S.: ‘Stochastic gene expression in a single cell science signalling’, Science, 2002, 297, (5584), pp. 1183–1186 3 Ozbudak, E.M., Thattai, M., Kurtser, I., Grossman, A.D., van Oudenaarden, A.: ‘Regulation of noise in the expression of a single gene’, Nature Genet., 2002, 31, (1), pp. 69–73 4 Paulsson, J.: ‘Summing up the noise in gene networks’, Nature, 2004, 427, (6973), pp. 415–418 5 Raser, J.M., O’Shea, E.K.: ‘Noise in gene expression: origins, consequences and control’, Science, 2005, 309, (5743), pp. 2010–2013 6 Kærn, M., Elston, T.C., Blake, W.J., Collins, J.J.: ‘Stochasticity in gene expression: from theories to phenotypes’, Nature Rev. Genet., 2005, 6, (6), pp. 451–464 7 Cai, L., Friedman, N., Xie, X.S.: ‘Stochastic protein expression in individual cells at the single molecule level’, Nature, 2006, 440, (7082), pp. 358–362 8 Xie, X.S., Choi, P.J., Li, G.-W., Lee, N.K., Lia, G.: ‘Single-molecule approach to molecular biology in living bacterial cells’, Annu. Rev. Biophys., 2008, 37, pp. 417–444 9 Cornish-Bowden, A.: ‘Fundamentals of Enzyme Kinetics’ (Wiley-Blackwell, 2012, 4th edn.) 10 Jia, C., Liu, X.-F., Qian, M.-P., Jiang, D.-Q., Zhang, Y.-P.: ‘Kinetic behavior of the general modifier mechanism of Botts and Morales with non-equilibrium binding’, J. Theory Biol., 2012, 296, pp. 13–20 11 Beard, D.A., Qian, H.: ‘Chemical Biophysics: Quantitative Analysis of Cellular Systems’ (Cambridge University Press, 2008, 1st edn.) 12 Keener, J.P., Sneyd, J.: ‘Mathematical Physiology’ (Springer, 1998, 1st edn.) 13 Gupta, P.B., Fillmore, C.M., Jiang, G., et al.: ‘Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells’, Cell, 2011, 146, (4), pp. 633–644 14 Segel, L.A., Goldbeter, A., Devreotes, P.N., Knox, B.E.: ‘A mechanism for exact sensory adaptation based on receptor modification’, J. Theory Biol., 1986, 120, (2), pp. 151–179 15 Knox, B.E., Devreotes, P.N., Goldbeter, A., Segel, L.A.: ‘A molecular mechanism for sensory adaptation based on ligand-induced receptor modification’, Proc. Natl. Acad. Sci., 1986, 83, (8), pp. 2345–2349 16 Behar, M., Hao, N., Dohlman, H.G., Elston, T.C.: ‘Mathematical and computational analysis of adaptation via feedback inhibition in signal transduction pathways’, Biophys. J., 2007, 93, (3), pp. 806–821 17 François, P., Siggia, E.D.: ‘A case study of evolutionary computation of biochemical adaptation’, Phys. Biol., 2008, 5, (2), pp. 026009 18 Ma, W., Trusina, A., El-Samad, H., Lim, W.A., Tang, C.: ‘Defining network topologies that can achieve biochemical adaptation’, Cell, 2009, 138, (4), pp. 760–773 19 Lan, G., Sartori, P., Neumann, S., Sourjik, V., Tu, Y.: ‘The energy-speed-accuracy trade-off in sensory adaptation’, Nature Phys., 2012, 8, (5), pp. 422–428
IET Syst. Biol., 2014, Vol. 8, Iss. 4, pp. 138–145 doi: 10.1049/iet-syb.2013.0050
20 Tu, Y., Shimizu, T.S., Berg, H.C.: ‘Modeling the chemotactic response of Escherichia coli to time-varying stimuli’, Proc. Natl. Acad. Sci., 2008, 105, (39), pp. 14855–14860 21 Tu, Y.: ‘Quantitative modeling of bacterial chemotaxis: signal amplification and accurate adaptation’, Annu. Rev. Biophys., 2013, 42, pp. 337–359 22 Mettetal, J.T., Muzzey, D., Gómez-Uribe, C., van Oudenaarden, A.: ‘The frequency dependence of osmo-adaptation in Saccharomyces cerevisiae’, Science, 2008, 319, (5862), pp. 482–484 23 Marchant, J.S., Taylor, C.W.: ‘Rapid activation and partial inactivation of inositol trisphosphate receptors by inositol trisphosphate’, Biochemistry, 1998, 37, (33), pp. 11524–11533 24 Adkins, C.E., Taylor, C.W.: ‘Lateral inhibition of inositol 1,4,5-trisphosphate receptors by cytosolic Ca2+’, Curr. Biol., 1999, 9, (19), pp. 1115–1118 25 Gyorke, S., Fill, M.: ‘Ryanodine receptor adaptation: control mechanism of Ca2+-induced Ca2+ release in heart’, Science, 1993, 260, (5109), pp. 807–809 26 Keizer, J., Levine, L.: ‘Ryanodine receptor adaptation and Ca2+-induced Ca2+ release-dependent Ca2+ oscillations’, Biophys. J., 1996, 71, (6), pp. 3477–3487 27 Fill, M., Zahradníková, A., Villalba-Galea, C.A., Zahradník, I., Escobar, A.L., Györke, S.: ‘Ryanodine receptor adaptation’, J. General Physiol., 2000, 116, (6), pp. 873–882 28 Li, Y., Goldbeter, A.: ‘Frequency specificity in intercellular communication. Influence of patterns of periodic signaling on target cell responsiveness’, Biophys. J., 1989, 55, (1), pp. 125–145 29 Ogata, K.: ‘Discrete-Time Control Systems’ (Prentice-Hall, 1995, 2nd edn.) 30 Golnaraghi, F., Kuo, B.C.: ‘Automatic Control Systems’ (Wiley, 2009, 9th edn.) 31 Allen, P.E., Holberg, D.R.: ‘Cmos Analog Circuit Design’ (Oxford University Press, USA, 2011, 3rd edn.) 32 Berman, A., Plemmons, R.J.: ‘Nonnegative Matrices in the Mathematical Sciences’ (Society for Industrial and Applied Mathematics, 1987, 1st edn.) 33 Zhou, D., Wu, D., Li, Z., Qian, M., Zhang, M.Q.: ‘Population dynamics of cancer cells with cell-state conversions between cancer stem cells and non-stem cancer cells’, Quant. Biol., 2013, 1, (3), pp. 201–208 34 Jiang, D.-Q., Qian, M., Qian, M.-P.: ‘Mathematical Theory of Nonequilibrium Steady States: On the Frontier of Probability and Dynamical Systems’ (Springer, 2004, 1st edn.) 35 Zhang, X.-J., Qian, H., Qian, M.: ‘Stochastic theory of nonequilibrium steady states and its applications. Part I’, Phys. Rep., 2012, 510, (1), pp. 1–86 36 Ge, H., Qian, M., Qian, H.: ‘Stochastic theory of nonequilibrium steady states. Part II: applications in chemical biophysics’, Phys. Rep., 2012, 510, (3), pp. 87–118 37 Qian, H.: ‘The mathematical theory of molecular motor movement and chemomechanical energy transduction’, J. Math. Chem., 2000, 27, (3), pp. 219–234 38 Metzner, P., Dittmer, E., Jahnke, T., Schütte, C.: ‘Generator estimation of Markov jump processes’, J. Comput. Phys., 2007, 227, (1), pp. 353–375
145
& The Institution of Engineering and Technology 2014
