A stochastic spectral analysis of transcriptional regulatory cascades

doi:10.1073/pnas.0811999106

A stochastic spectral analysis of transcriptional regulatory cascades

^aPrinceton Center for Theoretical Science, Princeton University, Princeton, NJ 08544; and
Departments of ^bPhysics and
^cApplied Physics and Applied Mathematics, Center for Computational Biology and Bioinformatics, Columbia University, New York, NY 10027

Edited by Peter G. Wolynes, University of California at San Diego, La Jolla, CA, and approved February 18, 2009 (received for review November 25, 2008)

Abstract

The past decade has seen great advances in our understanding of the role of noise in gene regulation and the physical limits to signaling in biological networks. Here, we introduce the spectral method for computation of the joint probability distribution over all species in a biological network. The spectral method exploits the natural eigenfunctions of the master equation of birth – death processes to solve for the joint distribution of modules within the network, which then inform each other and facilitate calculation of the entire joint distribution. We illustrate the method on a ubiquitous case in nature: linear regulatory cascades. The efficiency of the method makes possible numerical optimization of the input and regulatory parameters, revealing design properties of, e.g., the most informative cascades. We find, for threshold regulation, that a cascade of strong regulations converts a unimodal input to a bimodal output, that multimodal inputs are no more informative than bimodal inputs, and that a chain of up-regulations outperforms a chain of down-regulations. We anticipate that this numerical approach may be useful for modeling noise in a variety of small network topologies in biology.

Footnotes

¹To whom correspondence should be addressed. E-mail: awalczak{at}princeton.edu
Author contributions: A.M.W., A.M., and C.H.W. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/cgi/content/full/0811999106/DCSupplemental.
↵† The recursion can equivalently be performed in the reverse direction, with g_N = 0, g_N−1 = Np_N/p_N−1, and g_n−1 = (g_np_n + np_n − (n + 1)p_n+1)/p_n−1, where N is a cutoff in n. In several test cases we found that reverse recursion is more numerically stable than forward recursion at large N
↵‡ Setting x = e^ik₁ and y = e^ik₂ makes clear that the generating function is simply the Fourier transform.
↵§ G(x,y) can be recovered by projecting onto position space 〈x,y|, with 〈x|n〉 = xⁿ and 〈y|m〉 = y^m.
↵¶ The adjoint operations are 〈i|b̂_i⁺ = 〈i − 1|−〈i| for i ∈ {n,m}, 〈n|b̂_n⁻(n) = (n + 1)〈n + 1|−g_n〈n|, and 〈n,m|b̂_m⁻(n) = (m + 1)〈n,m + 1|−q_n〈n,m|.
↵∥ In position space the eigenfunctions are 〈x|j〉 = (x − 1)^je^(x−1) and 〈y|k〉 = (y − 1)^ke^(y−1). The operators b̂⁺ and $\text{[math]}$ ⁻ raise and lower in eigenspace: b̂_n⁺|j〉 = |j + 1〉 and $\text{[math]}$ |j〉 = j|j − 1〉 (or 〈j|b̂_n⁺ = 〈j − 1| and 〈j| $\text{[math]}$ = (j + 1)〈j + 1|), and similarly for n → m and j → k.
↵** The selection rules are derived by starting with 〈n|b̂_n⁺|j〉 or 〈j|b̂_n⁺|n〉 and allowing b̂_n⁺ to act both to the left and to the right. Alternatively, one may use $\text{[math]}$ , obtaining 〈n + 1|j〉 = (j〈n|j − 1〉 + $\text{[math]}$ 〈n|j〉)/(n + 1) and 〈j + 1|n〉 = (n〈j|n − 1〉− $\text{[math]}$ 〈j|n〉)/(j + 1), initialized with 〈n = 0|j〉 = (−1)^je⁻ and 〈j = 0|n〉 = 1. We find the latter relations yield smoother distributions p_nm for large cutoffs N and J.
↵†† Although the spectral method only involves the calculation of joint distributions between adjacent species in the cascade, the input-output distribution can be obtained by using p(n₁,n_ℓ) = ∑_{n_ℓ−1}p(n₁,n_ℓ−1)p(n_ℓ−1,n_ℓ)/p(n_ℓ−1), initialized with ℓ = 3 and run up to ℓ = L. This assumes p(n₁|n_ℓ−1,n_ℓ) = p(n₁|n_ℓ−1), which at worst (at ℓ = 3) is equivalent to the Markovian approximation, Eq. 2.
↵‡‡ There is a slight decrease in I* with 〈n〉 beginning near 〈n〉≈ 8 that is more pronounced at higher L. This is likely due to the decrease in accuracy of the Markovian approximation with increasing Δ (see SI Appendix), since large 〈n〉 requires large Δ. Calculations with the full joint distribution (via stochastic simulation) at L = 3 and 4 give qualitatively similar results, but with I* increasing monotonically with 〈n〉.
↵§§ The decrease of I* with L is consistent with, but not a direct consequence of, the data processing inequality (36), as each p*(n₁,…,n_L) results from a separate optimization for each subsequent choice of L.
↵* Although chemical kinetic rates can be nonlinear in the coordinate (e.g., through q_s or q_n in Eq. 1), the master equation itself is a linear equation in the unknown p.