The Master Equation Approach and the Langevin Approach

7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 Courtesy of Alexander Van Oudenaarden and Juan Pedraza. Used with permission. 1. THE MASTER EQUATION APPROACH The master equation corresponds to the statement that the probability of being in a given state changes depending on the probabilities of transition to and from any other state in the system. It provides the full probability distribution when it can be directly solved. Unfortunately, this is not often the case, so we must settle for some of the moments of the distribution. These are easily obtained from the generating function, so we will work with the master equation in a form in which it depends on the generating function rather than the distribution. The genetic network is defined by N state variables n1 .. nN and M rate constants k1 .. kM. The variables denote the number of copies of a certain chemical species such as mRNAs or proteins. Before applying the master equation approach to determine the noise properties of a genetic network we will start by obtaining the master equation in the generating function form for some elementary chemical equations: Synthesis from a template In numerous genetic reactions, such as transcription and translation, mRNAs and proteins are synthesized from a template (DNA and mRNA, respectively). After synthesis the number of templates is not changed. The corresponding reaction is therefore: k A → A+ B Molecule A produces molecule B at a rate k (in units of (concentration× time) ). The master equation describes how the probability to be in state [n1, n2] (n1 A molecules, n2 B molecules) at time t changes in time. For the reaction above: −1 p& (n1 , n2 ,t) = − kn1 p (n1 , n2 ,t) + kn1 p (n1 , n2 − 1, t) The first term reflects a transition from state [n1, n2] to state [n1, n2+1] and therefore leads to a decrease in p (n1 , n2 ,t) . The second term denotes the transition [n1, n2-1] → [n1, n2] and leads to an increased p (n1 , n2 ,t) . The master equation above is linear and can be solved for the moments by constructing the moment generating function. In general for N system variables the moment generating function is given by: F ( z1 , z 2 ,..., z N , t) = ∑z n1 1 n1 ,n2 ,...,n N z 2n2 ...z NnN p (n1 , n 2 ,..., n N , t) , 1 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 where the sum runs over all possible states for each ni (in this case, from 0 to ∞ ). This function has the following useful properties: F 1 = 1, where 1 ∂F = ni , ∂z i 1 ∂2F = ni (ni −1) , ∂z i2 1 ∂2F ∂z i ∂z j = ni n j , 1 means that the function is evaluated at zj = 1 for all j. These expressions justify the name “moment generating”: we can obtain the moments of the probability distribution by evaluating the partial derivatives of the function. Multiplying the master equation above by z1n1 z2n2 on both sides gives: ∑z n1 1 n1 ,n2 z 2n2 p& (n1 , n 2 , t) = − k ∑ z1n1 z 2n2 n1 p(n1 , n 2 , t) + k ∑ z1n1 z 2n2 n1 p(n1 , n 2 − 1, t) n1 ,n2 n1 ,n2 This equation can be simplified significantly by realizing that: ( ) ∂F ∂F = ∑ n1 z1n1 −1 p (n1 , n2 , t) ⇒ z1 = ∑ n1 z1n1 p (n1 , n2 ,t) ∂z1 n1 ∂z1 n1 ∑n z z p (n1 , n2 − 1, t) = z1 n1 n2 1 1 2 n1 = 0 , n2 =0 = z1 ∑ (n z )z n1 −1 1 1 n1 = 0 , n'2 =0 ∑ (n z )z n1 −1 1 1 n1 = 0 , n2 =−1 n'2 +1 2 n2 2 p (n1 , n2 − 1, t) p(n1 , n'2 , t) = z1 z 2 ∂F ∂z1 where the change in the lower limit of the sum for n2 is allowed because p(n1 ,−1, t) = 0. This leads to: ∂F F& (z1 , z2 ,t) = kz1 (z 2 −1) ∂z1 In the special case of synthesis from a fixed number of templates (n1 = n), the equation for the moment generating function reduces to: F& (z 2 , t) = kn(z 2 − 1)F This equation can be explicitly solved, but in itself it does not represent the full process. We therefore will obtain the expressions for the other terms before combining them to model a real situation. Degradation Now consider the degradation reaction: 2 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 γ B → 0. This reaction can represent two different processes: degradation, where molecule B is converted into a species which is not part of the subset of interest, and dilution, where it is physically separated from the volume of interest. In latter context γ is the degradation rate and ln(2)/γ the half-life of the molecule. The master equation for this reaction is: p& (n1 ,t) = −γn1 p(n1 , t) + γ (n1 + 1) p(n1 + 1, t) Using the same strategy as above the time evolution of the moment generation function yields: ∂F F& ( z1 , t ) = −γ ( z1 − 1) ∂z1 Forward reaction, conservation of total number of molecules Now consider the reaction: k A → B where no + n1 = n = const. Since the total number n is conserved, the system is defined by only one variable. We will use n2 as the single variable to define this system. For the reaction above: p& (n 2 , t) = − k (n − n 2 ) p (n 2 , t) + k (n − n 2 + 1) p (n2 − 1, t) This leads1 to:   ∂F ∂F ∂F F& ( z 2 , t) = − knF + kz2 + k (n + 1) z2 F − kz2  z 2 + F  = kn( z2 − 1)F − kz2 (z2 − 1) ∂z1 ∂z2   ∂z2 Based on these elementary reactions larger chemical networks can be built up. The results above are summarized in Table 1. In this case, the sums only go up to n, instead of ∞ . However, the extra terms that appear when applying the change of variables cancel with each other. 1 3 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 Reaction Type I k A → A+ B II γ B → 0 III .IV k A → B (n1 + n2= n = const., in terms of n2) k1 → A B ← k −1 F& = ∂F ∂z1 ∂F − γ ( z1 − 1) ∂z1 kz1 ( z 2 − 1) kn( z 2 − 1) F − kz 2 ( z 2 − 1) k1 ( z 2 − z1 ) ∂F ∂z 2 ∂F ∂F + k −1 ( z1 − z 2 ) ∂z1 ∂z 2 Table 1. Moment generating function equations for elementary reactions. The master equation for each reaction type produces different terms which can be combined to model more complex processes. Noise properties of a constitutively expressed gene Based on the results for these elementary reactions the equation for the moment generating functions of more complex networks can be easily deduced. First let us consider a constitutive expressed gene in a single copy in the chromosome of a bacterium. In this case the state of this system at any time is defined by the number of mRNA molecules r and number of proteins p for that gene. mRNA molecules are synthesized off the template DNA strand at a rate kR and are translated at a rate kP. The mRNA and protein degradation are described by the destruction rates γ R and γ P respectively (Fig. 1). Figure 1. Basic model for constitutive expression of a single gene. Only four individual reactions are considered: creation of mRNA from a DNA template, creation of proteins from individual mRNA molecules, and the degradation/dilution of both species. 4 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 Based on the results in Table 1 the moment generating function can be deduced directly: ∂F ∂F ∂F − γ R ( z1 − 1) F& ( z1 , z2 , t ) = k R ( z1 − 1) F + k P z1 ( z2 − 1) − γ P ( z2 − 1) ∂z1 ∂z1 ∂z 2 The first two terms are the transcription and translation reactions (Table 1, type I) and the last two terms model degradation of mRNA and proteins respectively (Table 1, type II). Below the equation will be solved for the moments in the steady state ( F& = 0 ). In this case: k R (1 − z1 ) F = k P z1 ( z 2 − 1) ∂F ∂F ∂F − γ R ( z1 − 1) − γ P ( z 2 − 1) ∂z1 ∂z1 ∂z2 The mean mRNA level r and protein level p are found by taking the derivative with respect to z1 and z2 respectively: k R (1 − z1 ) ∂F ∂2F ∂F ∂2F − k R F = [k P z1 ( z2 − 1) − γ R ( z1 − 1)] 2 + [k P ( z 2 − 1) − γ R ] − γ P ( z 2 − 1) ∂z1 ∂z1 ∂z1 ∂z2∂z1 evaluating both expressions at z1 = z2 =1 gives: r = p = kR γR kR kP γ Rγ P The results are consistent with the equivalent deterministic system: r& = k R − γ R r p& = k P r − γ P p The fluctuations in mRNA and proteins level are found by differentiating the above equations again with respect to z1 and z2 and evaluating at z1 = z2 = 1: r2 − r 2 = r rp − r p = p2 − p 2 p 1+ γ R /γ P   k /γ = p  P R + 1  1+ γ P /γ R 5 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 Further moments can be obtained sequentially in this manner. Note that for a random variable with Poissonian distribution, all moments are equal, so the variance over the mean equals one, as is the case for the mRNA in this model. The protein number fluctuates with a higher than poissonian noise, the correction determined primarily by the term k P / γ R (the ”burst size”), which corresponds to the average number of proteins produced per mRNA [9]. In simple cases like this, the moments can also be obtained as a function of time. For a single gene, the noise out of equilibrium can be 40% larger than its steady state value in the limit of short mRNA lifetimes [9]. A more detailed modeling of this process could include more intermediate detailed processes, such as the random steps that a ribosome takes along an mRNA, but most turn out to have little effect when compared in simulations. However, when a repressor or activator is present, its binding and unbinding might have to be included in the model, for this can be a major source of noise. It is in this context that the terms shown in Table 1, type III, are needed. Furthermore, the repressor concentration itself might be fluctuating, in which case we have to consider the entire system of genes. Linearized matrix formulation The method above can also be used for interacting systems of genes, but solving it is not straightforward unless the connections are linear. Alternatively, if the system is at a stable point in steady state, the interaction can be linearized around the steady state value. A practical way of writing this out is in matrix form. The transition probabilities for species xi are given by fi(x1,x2,...,xn) for creation and γi for destruction, and A and Γ are the ∂f matrices defined by Aij = i and Γij = γ iδ ij . Letting x be the vector of chemical ∂x j x1 , x2 ... species, the linearized macroscopic equations are then given by x& = ( A − Γ ) x . Note that since in many cases the macroscopic equations include constant creation terms. If the system is linear it might be necessary to include an additional variable, which is not fluctuating and allows the inclusion of the constant terms in the compact matrix form. As an illustration of this, the matrices for the single gene case are 0 A = k R  0 0 0 kP 0 0 0   0  , Γ = 0 γ R 0 0 0 0 0  . γ P  where the state vector is x T = (d , r , p ) where d is the gene copy number. This constant state coordinate needs not to represent an actual chemical; for a system where many species have a constant creation rate, these rates can all be placed in the first column of A (setting d=1 and Γ1 j = Γ j1 = 0 ). An example of this is the matrix for the case of two 6 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 interacting genes, linearized around steady state, with fixed gene copy numbers d1 and d2 respectively, and where the first gene (r1,p1) represses the second (r2,p2) with transfer function f(p1):  0 d k  1 R1  0 A= d 2 k R 2    0 0 0 k P1 where k R 2 = f ( p1 ) − 0 0 0 ∂f d2 ∂p1 0 0 0 ∂f ∂p1 0 0 0 0 p1 kP2 0 0 0 , 0   0 0 0 0  0 0 0 γ 0 0 0  R1  Γ = 0 0 γ P1 0 0 ,   0 γ R2 0  0 0 0 0 0 0 γ P 2  p1 . p1 Written in terms of these matrices, the master equation in generating function form would be  ∂F ∂F  − ∑ Aij z j F& = ∑ (1 − z i ) Γii .  ∂z  ∂ z i j i j   At steady state, F& = 0 , and taking the derivative with respect to zl we obtain:  ∂2F ∂2F ∂F   ∂F ∂F  − ∑ Aij z j − Ail − Γll − ∑ Alj z j 0 = ∑ (1 − z i ) Γii  ∂z ∂z ∂z j ∂z l ∂z l   ∂z l ∂z j  j j i i l  Setting all zi=1, we have for each i 0 = Γii ∂F ∂F − ∑ Aij z j ⇒ 0 = (Γ − A)∇ F 1 = (Γ − A) x1 , ∂z i 1 j ∂z j 1 corresponding to the macroscopic result. Similarly, differentiating again and evaluating at zi=1,  ∂F ∂2F ∂2F  0 = Γii − Ail − ∑ Aij z j  ∂z i ∂z l ∂zl ∂z j ∂z l j 1 1  ( ) ( 2 2    +  Γ ∂ F − A z ∂ F − A ∂F ∑j lj j ∂z ∂z li   ll ∂z l ∂zi ∂z i j i 1 1  1 = (Γ − A)∇∇ T F 1 − AΘ F 1 + (Γ − A)∇∇ T F 1 − AΘ F 1 )    1, T 7 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 where Θij = δ ij ∂ . These linear equations can be solved for the means, variances and ∂zi correlations. This approach is very general and the resulting matrix equations can be solved directly2. However, even for the case of just two interacting genes this requires a 5x5 matrix system as shown, so it gets cumbersome for larger systems even though most entries are zero. Using symbolic matrix manipulation software it is straightforward to obtain the desired expressions, so for known parameters this is a good method for obtaining values without further approximations. 2. THE LANGEVIN APPROACH An alternate approach that allows for a more straightforward interpretation and scales easily to different levels of detail is the use of a Langevin equation. The Langevin approach consists essentially of adding a noise term to the deterministic equations. This noise term can represent the effect of the intrinsic fluctuations [20] or the external inputs of the system [21]. For x, the concentration of some chemical species, x& = f ( x) → x& = f ( x) + q ( x)ε (t ) , where the random variable ε (t ) is determined by its statistical properties. Formally, this can be any random process, but in practice we assume white-noise statistics, which will give approximate values for the first two moments. The conditions for white noise are: ε (t ) = 0, where ε (t )ε (t + τ ) = δ (τ ), denotes an ensemble average. Since we are interested in the steady state fluctuations, we will assume the coefficient of the noise term to be constant3, i.e. evaluated at x ss . For the case of our basic model of the single gene, we have two macroscopic equations representing mRNA and protein creation, respectively: r& = k R − γ R r → r& = k R − γ R r + qrε r p& = k P r − γ P p → p& = k P r − γ P p + q pε p , 2 This can be summarized in a very practical way [13] in terms of the logarithmic gains to obtain an equation which reflects the resulting components of the noise. 3 For the case where q(x) is not constant, the stochastic differential equation will be understood to follow the Stratonovich interpretation [19,22]. This allows a general Fokker-Planck equation to be written in this form, but will not be necessary in the cases of interest. 8 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 where the coefficients of the noise terms are to be determined. Clearly, r = k R / γ R and p = k P r / γ P from the condition of zero mean for the noise term. The difference with the steady state δr = r − r follows the equation δr& + γ Rδr = qr ε r Fourier transforming, we obtain qr εˆr , iω + γ R δr̂ (ω ) = so after multiplying by the complex conjugate and taking the average, δr̂(ω ) 2 2 = qr 2 ω + γ R2 The steady state fluctuations are given by the inverse Fourier transform with t = 04: 1 δr = 2π 2 ∞ 2 2 q q ∫−∞ω 2 +r γ R2 dω = 2πγr R ∞ 2 q dx ∫−∞ x 2 +1 = 2γr R . But since the production of mRNA is in this model a single step, independent random process, it has a Poisson distribution, so the variance equals the mean, which implies 2 qr k 2 = R ⇒ qr = 2k R . 2γ R γ R For the number of proteins, we have δp& + γ Pδp = δr + q pε p , δrˆ( w) + q p εˆ p , δp̂(ω ) = iω + γ P but in this case we also need to notice that δrˆ(w)εˆ*p = δrˆ(w) εˆ*p = 0 , since these are two independent random processes with zero mean. So in this case, δp̂ (ω ) 4 2 = δr̂ (ω ) 2 + qp ω 2 + γ P2 2 = (ω qr 2 2 )( + γ R2 ω 2 + γ P2 )+ ω qp 2 2 + γ P2 . From the Wiener-Khintchine theorem; see [23]. 9 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 2 2 Using qr = 2k R and q p = 2k p kR γR , (since this represents the internal noise and for a fixed number of mRNAs the production of proteins is also a Poissonian process). Performing the inverse transform5, 1 δp = 2π 2 ∞ ∫ (ω qr 2 −∞ 2 )( + γ R2 ω 2 + γ P2 ) + qp 2 ω 2 + γ P2 dω = 2 kP kR  1 1 k k − + P R 2 2  γ R − γ P  γ P γ R  γ Pγ R ( for comparison with the previous result, note that p = ) kPkR γ Pγ R , so this can be rewritten as   k /γ δp 2 = p  P R + 1  1+ γ P /γ R This is identical to the result obtained by the master equation. This method can be readily generalized for many interacting genes when the system is fluctuating around a steady state. As an example, we will analyze the case where one gene represses a second gene. Let y0, y1 be the protein numbers of each gene, and let f(y0) be the rate of creation of the second protein as a function of the first. This means that the equations describing this system are y& 0 = k0 − γ 0 y0 , y&1 = f ( y0 ) − γ 1 y1 . Note that the equations include the entire process of producing a protein, so mRNA levels are no longer explicitly calculated. Including the Langevin noise term and looking at the fluctuations from steady state, δy& 0 = −γ 0δy0 + q0ε 0 , δy&1 = f ( y0 ) − f ( y0 ) − γ 1δy1 + q1ε1 ≈ c0δy0 − γ 1δy1 + q1ε1 , where c0 = df dy0 and each noise term has the same conditions as before. This y0 linearization is valid at each stable point, but not for transitions between different stable points or for limit cycles. For very small numbers n of chemicals this also breaks down, because since this processes are mostly Poissonian, the fluctuations are of order n so a Taylor expansion might not be valid. Fourier transforming and taking the square and the average as before, we get 5 1 2π ∫ (ω dω 2 +γ ) 2 n = 1 2 πγ 2n−1 Γ(n −1/ 2) , where Γ(n ) = Γ(n −1)(n −1) , Γ(1) = 1 and Γ(1/ 2) = π . Γ(n ) 10 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 δŷ0 (ω ) 2 δŷ1 (ω ) 2 2 = = q0 2 ω + γ 02 c02 δŷ 0 (ω ) 2 + q1 2 = ω 2 + γ 12 (ω c02 q0 2 2 2 )( + γ 12 ω 2 + γ 02 ) q + 21 2 ω +γ1 The correlation between the genes can also be calculated, from δŷ1δŷ * 0 c0 δŷ0 (ω )  c δŷ (w) + q1εˆ1  * δŷ0 (w) = =  0 0 iω + γ 1 iω + γ 1   2 = c0 q 0 2 (iω + γ 1 )(ω 2 + γ 02 ) , where the term εˆ1δyˆ 0* vanishes because the fluctuations in the first gene are independent from the internal fluctuations in the second gene. In many cases, the decay time will be determined primarily by the dilution time, so it will be the same for all genes. This assumption simplifies the expressions that are obtained upon transforming back: δy 2 0 1 = 2π δy12 = 1 2π δy1δy0 = ∞ 2 ∞ c02 q0 2 q0 q0 ∫−∞ω 2 + γ 2 dω = 2γ ∫ (ω −∞ 1 2π 2 ∞ 2 +γ 2 2 2 ) 2 + 2 c02 q0 q q1 d ω = + 1 3 2 2 2γ 4γ ω +γ c0 q0 2 ∫ (iω + γ )(ω 2 −∞ +γ 2 ) dω = c0 q0 2π (iω − γ ) c q ∫ (ω + γ ) dω = 4γ 2 ∞ 0 2 −∞ 2 2 2 0 2 , where the irrational part of the integral vanishes because of parity. From our previous results we know that for a single gene, δy 2 0 2   k /γ q0 2 = = y 0  Po Ro + 1 ≈ y 0 (b0 + 1) ⇒ q 0 ≈ 2γ y 0 (b0 + 1) , 2γ   1 + γ / γ Ro where bi is the burst size for gene i. For basic Hill-type repression, k1 − k1 f ( y0 ) = + B1 ⇒ c0 = h h 1 + ( y 0 / Y1/ 2 ) 1 + ( y0 / Y1/ 2 ) (  y0  2   Y1/ 2 )     h −1 h = y1 Y1/ 2 2  y0  Y  1/ 2     h −1 − hγ 2 k1Y1/ 2 where k1+B1 is the maximum creation rate, Y1/2 is the half induction point, h is the Hill coefficient and B1 is the basal transcription level. Assuming that the internal noise for the second gene alone has the same form, 11 7.81/8.591/9.531 Systems Biology – A. van Oudenaarden and Mukund Thattai – MIT– October 2004 q1 ≈ 2γ y1 (b1 + 1) , 2 the variance and correlation can be explicitly written as δy 2 1  y =  0  Y1 / 2 δy1δy0 = y1     2h 2  hγ    2 k  1  y0  Y  1/ 2 2 y1 y0 2 2(b0 + 1) + y1 (b1 + 1) h  − hγ   2k (b0 + 1) . 1  Note that we need the parameters of the macroscopic equations plus an “internal” parameter for each gene, bi = k Pi / γ Ri which depends on the parameters of the macroscopic equations for each gene. 12