Abstract
In previous work, we have introduced a “linear framework” for time-scale separation in biochemical systems, which is based on a labelled, directed graph, G, and an associated linear differential equation, \(dx/dt = \mathcal{L}(G)\cdot x\), where \(\mathcal{L}(G)\) is the Laplacian matrix of G. Biochemical nonlinearity is encoded in the graph labels. Many central results in molecular biology can be systematically derived within this framework, including those for enzyme kinetics, allosteric proteins, G-protein coupled receptors, ion channels, gene regulation at thermodynamic equilibrium, and protein post-translational modification. In the present paper, in response to new applications, which accommodate nonequilibrium mechanisms in eukaryotic gene regulation, we lay out the mathematical foundations of the framework. We show that, for any graph and any initial condition, the dynamics always reaches a steady state, which can be algorithmically calculated. If the graph is not strongly connected, which may occur in gene regulation, we show that the dynamics can exhibit flexible behavior that resembles multistability. We further reveal an unexpected equivalence between deterministic Laplacian dynamics and the master equations of continuous-time Markov processes, which allows rigorous treatment within the framework of stochastic, single-molecule mechanisms.
Similar content being viewed by others
References
Ackers, G. K., Johnson, A. D., & Shea, M. A. (1982). Quantitative model for gene regulation by lambda phage repressor. Proc. Natl. Acad. Sci. USA, 79, 1129–1133.
Agaev, R. P., & Chebotarev, P. Y. (2000). The matrix of maximum out forests of a digraph and its applications. Autom. Remote Control, 61, 1424–1450.
Ahsendorf, T., Wong, F., Eils, R., & Gunawardena, J. (2013, in preparation). A framework for modelling eukaryotic gene regulation that accommodates non-equilibrium mechanisms.
Bintu, L., Buchler, N. E., Garcia, G. G., Gerland, U., Hwa, T., Kondev, J., Kuhlman, T., & Phillips, R. (2005a). Transcriptional regulation by the numbers: applications. Curr. Opin. Genet. Dev., 15, 125–135.
Bintu, L., Buchler, N. E., Garcia, G. G., Gerland, U., Hwa, T., Kondev, J., & Phillips, R. (2005b). Transcriptional regulation by the numbers: models. Curr. Opin. Genet. Dev., 15, 116–124.
Bott, R., & Mayberry, J. P. (1954). Matrices and trees. In O. Morgenstern (Ed.), Economic activity analysis (pp. 391–400). New York: Wiley
Cairns, B. R. (2009). The logic of chromatin architecture and remodelling at promoters. Nature, 461, 193–198.
Chebotarev, P. (2010). Comment on ‘Consensus and cooperation in networked multi-agent systems’. Proc. IEEE, 98, 1353–1354.
Chebotarev, P., & Agaev, R. (2002). Forest matrices around the Laplacian matrix. Linear Algebra Appl., 356, 253–274.
Chebotarev, P. Y., & Agaev, R. P. (2009). Coordination in multiagent systems and Laplacian spectra of digraphs. Autom. Remote Control, 70, 469–483.
Chen, W. K. (1971). Applied graph theory. In Applied mathematics and mechanics, Amsterdam: North-Holland.
Chung, F. R. K. (1997). Spectral graph theory. Regional conference series in mathematics: Vol. 92. Providence: Am. Math. Soc.
Colquhoun, D. (2006). The quantitative analysis of drug-receptor interactions: a short history. Trends Pharmacol. Sci., 27, 149–157.
Cornish-Bowden, A. (1995). Fundamentals of enzyme kinetics (2nd ed.). London: Portland Press.
Dasgupta, T., Croll, D. H., Owen, J. A., Vander Heiden, M. G., Locasale, J. W., Alon, U., Cantley, L. C., & Gunawardena, J. (2013). A fundamental trade off in covalent switching and its circumvention in glucose homeostasis. Submitted.
Feinberg, M., & Horn, F. (1977). Chemical mechanism structure and the coincidence of the stoichiometric and kinetic subspace. Arch. Ration. Mech. Anal., 66, 83–97.
Gertz, J., Siggia, E. D., & Cohen, B. A. (2009). Analysis of combinatorial cis-regulation in synthetic and genomic promoters. Nature, 457, 215–218.
Gunawardena, J. (2012). A linear framework for time-scale separation in nonlinear biochemical systems. PLoS ONE, 7, e36321.
He, X., Samee, M. A. H., Blatti, C., & Sinha, S. (2010). Thermodynamics-based models of transcriptional regulation by enhancers: the roles of synergistic activation, cooperative binding and short-range repression. PLoS Comput. Biol., 6, e1000935.
Hill, T. L. (1966). Studies in irreversible thermodynamics IV. Diagrammatic representation of steady state fluxes for unimolecular systems. J. Theor. Biol., 10, 442–459.
Hill, T. L. (1985). Cooperativity theory in biochemistry: steady-state and equilibrium systems. Springer series in molecular biology. New York: Springer.
Hirsch, M. W., & Smale, S. (1974). Differential equations, dynamical systems and linear algebra. Pure and applied mathematics. San Diego: Academic Press.
Hopfield, J. J. (1974). Kinetic proofreading: a new mechanism for reducing errors in biosynthetic processes requiring high specificity. Proc. Natl. Acad. Sci. USA, 71, 4135–4139.
Horn, R. A., & Johnson, C. A. (1985). Matrix analysis. Cambridge: Cambridge University Press.
Jaffe, A. (1965). Divergence of perturbation theory for bosons. Commun. Math. Phys., 1, 127–149.
Janssens, H., Hou, S., Jaeger, J., Kim, A. R., Myasnikova, E., Sharp, D., & Reinitz, J. (2006). Quantitative and predictive model of transcriptional control of the drosophila melanogaster even skipped gene. Nat. Genet., 38, 1159–1165.
van Kampen, N. G. (1992). Stochastic processes in physics and chemistry. Amsterdam: Elsevier.
Kelly, F. P. (2011). Reversibility and stochastic networks. Cambridge: Cambridge University Press.
Kenakin, T. (2005). New concepts in drug discovery: collateral efficacy and permissive antagonism. Nat. Rev. Drug Discov., 4, 919–927.
Kim, H. D., & O’Shea, E. K. (2008). A quantitative model of transcription factor-activated gene expression. Nat. Struct. Mol. Biol., 15, 1192–1198.
King, E. L., & Altman, C. (1956). A schematic method of deriving the rate laws for enzyme-catalyzed reactions. J. Phys. Chem., 60, 1375–1378.
Kirchhoff, G. (1847). Über die Auflösung der Gleichungen, auf welche man bei der Untersuchung der linearen Verteilung galvanischer Ströme geführt wird. Ann. Phys. Chem., 72, 497–508.
Koshland, D. E., Némethy, G., & Filmer, D. (1966). Comparison of experimental binding data and theoretical models in proteins containing subunits. Biochemistry, 5, 365–385.
Kuhlman, T., Zhang, Z., Saier, M. H. Jr., & Hwa, T. (2007). Combinatorial transcriptional control of the lactose operon of Escherichia coli. Proc. Natl. Acad. Sci. USA, 104, 6043–6048.
Lam, F. H., Steger, D. J., & O’Shea, E. K. (2008). Chromatin decouples promoter threshold from dynamic range. Nature, 453, 246–250.
Lean, A. D., Stadel, J. M., & Lefkowitz, R. J. (1980). A ternary complex model explains the agonist-specific binding properties of the adenylate cyclase-coupled β-adrenergic receptor. J. Biol. Chem., 255, 7108–7117.
Lewis, G. N. (1925). A new principle of equilibrium. Proc. Natl. Acad. Sci. USA, 11, 179–183.
Magnus, J. R., & Neudecker, H. (1988). Matrix differential calculus with applications in statistics and econometrics. Chichester: Wiley
Merris, R. (1994). Laplacian matrices of graphs: a survey. Linear Algebra Appl., 198, 143–176.
Michaelis, L., & Menten, M. (1913). Die kinetik der Invertinwirkung. Biochem. Z., 49, 333–369.
Mirny, L. (2010). Nucleosome-mediated cooperativity between transcription factors. Proc. Natl. Acad. Sci. USA, 107(22), 534–539.
Monod, J., Wyman, J., & Changeux, J. P. (1965). On the nature of allosteric transitions: a plausible model. J. Mol. Biol., 12, 88–118.
Moon, J. W. (1970). Counting Labelled Trees. Canadian mathematical monographs: Vol. 1. Ottawa: Canadian Mathematical Society.
Nishikawa, T., & Motter, A. E. (2010). Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions. Proc. Natl. Acad. Sci. USA, 107(10), 342–347.
Olfati-Saber, R., Fax, J. A., & Murray, R. M. (2007). Consensus and cooperation in networked multi-agent systems. Proc. IEEE, 95, 215–233.
Olfati-Saber, R., & Murray, R. M. (2004). Consensus problems in networks of agents with switching topology and time-delays. IEEE Trans. Autom. Control, 49, 1520–1533.
Pecora, L. M., & Carroll, T. L. (1998). Master stability functions for synchronized coupled systems. Phys. Rev. Lett., 80, 2109–2112.
Raveh-Sadka, T., Levo, M., & Segal, E. (2009). Incorporating nucleosomes into thermodynamic models of transcription regulation. Genome Res., 19, 1480–1496.
Schnakenberg, J. (1976). Network theory of microscopic and macroscopic behaviour of master equation systems. Rev. Mod. Phys., 48, 571–586.
Segal, E., Raveh-Sadka, T., Schroeder, M., Unnerstall, U., & Gaul, U. (2008). Predicting expression patters from regulatory sequence in Drosophila segmentation. Nature, 451, 535–540.
Segal, E., & Widom, J. (2009). From DNA sequence to transcriptional behaviour: a quantitative approach. Nat. Rev. Genet., 10, 443–456.
Setty, Y., Mayo, A. E., Surette, M. G., & Alon, U. (2003). Detailed map of a cis-regulatory input function. Proc. Natl. Acad. Sci. USA, 100, 7702–7707.
Sherman, M. S., & Cohen, B. A. (2012). Thermodynamic state ensemble models of cis-regulation. PLoS Comput. Biol., 8, e1002407.
Stamatoyannopoulos, J. (2012). What does our genome encode? Genome Res., 22, 1602–1611.
Thomson, M., & Gunawardena, J. (2009a). The rational parameterisation theorem for multisite post-translational modification systems. J. Theor. Biol., 261, 626–636.
Thomson, M., & Gunawardena, J. (2009b). Unlimited multistability in multisite phosphorylation systems. Nature, 460, 274–277.
Tirosh, I., & Barkai, N. (2008). Two strategies for gene regulation by promoter nucleosomes. Genome Res., 18, 1084–1091.
Tolman, R. C. (1938). The principles of statistical mechanics. Oxford: Clarendon Press.
Tutte, W. T. (1948). The dissection of equilateral triangles into equilateral triangles. Proc. Camb. Philol. Soc., 44, 463–482.
Tutte, W. T. (2001). Graph theory. Encyclopedia of mathematics and its applications: Vol. 21. Cambridge: Cambridge University Press.
Whitin, T. M. (1954). An economic application of ‘Matrices and trees’. In O. Morgenstern (Ed.), Economic activity analysis (pp. 401–418). New York: Wiley
Xu, Y., & Gunawardena, J. (2012). Realistic enzymology for post-translational modification: zero-order ultrasensitivity revisited. J. Theor. Biol., 311, 139–152.
Zaher, H. S., & Green, R. (2009). Fidelity at the molecular level: lessons from protein synthesis. Cell, 136, 746–762.
Zinzen, R. P., Senger, K., Levine, M., & Papatsenko, D. (2006). Computational models for neurogenic gene expression in the Drosophila embryo. Curr. Biol., 16, 1358–1365.
Acknowledgements
Author information
Authors and Affiliations
Corresponding author
Appendix
Appendix
We include for completeness a proof of the Matrix-Tree theorem as stated in Theorem 1. There are many versions of this and the one given here is adapted from David Perkinson’s handout,Footnote 1 which is itself adapted from the proof in Tutte’s book (Tutte 2001). This has the merits of being elementary, concise, and transparent. Let G be a graph on the vertices 1,…,n and let the symbols e ij denote the labels on the edges: j→i, if, and only, if e ij ≠0. The entries of \(\mathcal{L}(G)\) are then given by Eq. (2).
The proof of Theorem 1 is in two steps, first for the special case of a principal minor, which is the essential part of the argument. Some additional notation will be helpful. If A is any (not necessarily square) matrix let A [u,v] denote the matrix obtained from A by removing row u and column v. It is easier to work with the negative Laplacian, \(L = -\mathcal{L}(G)\), to keep control of the signs. Since the minor is the determinant of a (n−1)×(n−1) matrix, this change incurs a sign of (−1)n−1, so we can rewrite Eq. (16) for the jth principal minor as
which is what we want to prove.
Proof
The argument does not depend on j, so take j=1 to simplify the notation. The determinant can be written in the standard form as
where π:{2,…,n}→{2,…,n} is a permutation on {2,…,n}, S n−1 is the symmetric group of all such permutations and \(\operatorname{sign}(\pi) = \pm1\) is the sign of the permutation. A permutation can be uniquely expressed as a product of disjoint cycles, where the cycle (i 1 i 2⋯i k ), with k>1, is the permutation that sends i v to i v+1 for v<k and cycles around to send i k to i 1. A permutation may also fix an element by sending j to j and this is denoted by a degenerate cycle, (j). For example,
A nonzero summand in Eq. (43), corresponding to the permutation π, can be written in terms of nonzero symbols e ij ≠0 as follows. As mentioned above, π is a product of disjoint cycles. Each cycle gives rise to a product of nonzero symbols,
where the signs arise because L is the negative of the Laplacian matrix defined by Eq. (2). Each fixed point of π, π(j)=j, introduces into the overall product the sum of all the nondiagonal symbols in column j of L,
of which at least one symbol e vj must be nonzero in order for the summand itself to be nonzero. Each cycle of odd length contributes a sign of −1 through Eq. (44), while each cycle of even length contributes a sign of −1 to \(\operatorname{sign}(\pi)\). The resulting overall product, when expanded out, is therefore a sum of monomials of the form
where each symbol is nonzero, \(e_{k_{i}i} \neq0\), and t is the number of cycles in π.
Monomials of the form in Eq. (46) have a natural interpretation in terms of the graph G. They correspond to subgraphs of G in which there is exactly one outgoing edge from the vertices 2,…,n and no outgoing edge from vertex 1. Lets call these special subgraphs. Given a special subgraph γ, it has an associated monomial e γ which is the product of the symbols on all the edges and it is evident that this correspondence defines a bijection between special subgraphs and the monomials in Eq. (46).
Special subgraphs have a distinctive structure: If we start at any vertex other than 1, then there is a unique path leading away from this vertex. Because G has only finitely many vertices, this path must either intersect itself to form a cycle of edges or reach vertex 1 and stop. It follows that each connected component of a special subgraph either contains an unique cycle or is a tree rooted at 1. (Recall that the connected components of a directed graph are the disjoint pieces into which the graph falls if edge directions are ignored.) Of course, the cycle components may also contain tree-like spurs that lead to its unique cycle. Equation (44) shows that each cycle of edges in G may be uniquely identified with a permutation cycle.
Suppose given a permutation π∈S n−1 and a special subgraph γ. If all the cycles of π appear as the cycles of cycle components of γ, then the symbols of any other edges in γ, including those of any other cycle components of γ, can be constructed in π by using the elements fixed by π, since these introduce, through the terms in Eq. (45), all the symbols needed for these remaining edges. It follows that e γ is one of the monomials that appear when π is expanded into a sum of monomials through Eq. (46). Conversely, if e γ appears in the monomial expansion of π then the cycles of π must occur among the cycles of cycle components of γ.
Given a special subgraph γ, we can now determine the multiplicity of e γ in the determinant in Eq. (43). Suppose that γ has s cycle components, where s>0. Any permutation π∈S n−1 which uses exactly t of these cycles gives rise to e γ with multiplicity (−1)t according to Eq. (46). There are \(\binom{s}{t}\) ways of choosing t cycle components out of s. Considering all possibilities for t from t=0 to t=s, the total multiplicity of e γ is
This amazing cancellation is the heart of the Matrix-Tree theorem.
The only remaining possibility is that γ has no cycle components, so that s=0. But then γ is a spanning tree rooted at 1. Moreover, e γ only appears as a monomial in the expansion of the identity permutation, whose monomial expansion is given by
Evidently, e γ appears in Eq. (48) with multiplicity 1. Of course, Eq. (48) gives rise to monomials other than those coming from spanning trees rooted at 1 but these other monomials disappear through cancellation because of Eq. (47). We are left with exactly the sum of monomials over rooted spanning trees in Eq. (42). □
To complete the proof of Theorem 1, it remains to deal with a non-principal minor, detL [i,j], which contributes the sign (−1)i+j to Eq. (16). This has nothing to do with Laplacian matrices and depends only on the fact that \(1\cdot\mathcal{L}(G) = 0\). Some further notation will be helpful. Let A be a n×(n−1) matrix and let A i∗ denote its i-th row. We want to consider various (n−1)×(n−1) matrices made from these rows and we will write these matrices by listing their rows in order. For instance, the (n−1)×(n−1) submatrix of A in which the first row is omitted, is written [A 2∗|A 3∗|⋯|A n∗]. We will also use the convention that, in an ordered sequence, a bar over a term signifies its absence, so we can also write this same submatrix as \([\overline{A_{1*}}|A_{2*}| \cdots| A_{n*} ]\).
Lemma 4
If A is any n×(n−1) matrix for which 1⋅A=0, then
Proof
Let \(B = [ A_{1*}+A_{i*} |A_{2*}| \cdots|\overline{A_{i*}}| \cdots| A_{n*} ]\). By construction, 1⋅B=0, so that detB=0. Because the determinant is multilinear,
The determinant is also antisymmetric for interchange of rows, so that
It takes i−2 successive row interchanges to bring A i∗ back to where it originally was in the list, from which the result follows. □
Now let A be the n×(n−1) matrix obtained from \(L = -\mathcal{L} (G)\) by removing the jth column. Then, using Lemma 4 twice, we see that
Combining this with Eq. (42) yields Eq. (16) and proves Theorem 1.
Rights and permissions
About this article
Cite this article
Mirzaev, I., Gunawardena, J. Laplacian Dynamics on General Graphs. Bull Math Biol 75, 2118–2149 (2013). https://doi.org/10.1007/s11538-013-9884-8
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11538-013-9884-8