A framework for modelling gene regulation which accommodates non-equilibrium mechanisms

A graph-theoretic view of gene regulation

We offer in this section a non-technical account of the linear framework as applied to gene regulation. The technical details are provided, along with references, in the section on ‘Calculating microstate probabilities at steady state’.

The framework starts with a labelled, directed graph consisting of a collection of vertices with directed edges between pairs of vertices and labels on the edges (Figure 1, bottom). The graphs considered here have only finitely many vertices and the edges always go between distinct vertices, so that there are no self-loops. It is further assumed that each graph is connected, which means that, given any two vertices, there is always a path of edges between them, ignoring edge directions. A connected graph is not in disjoint pieces.

The vertices of the graph correspond to microstates, or snapshots of DNA and its accompanying proteins. Figure 1 (top) shows the range of features that can potentially be found in a microstate, including TFs, transcriptional co-regulators, RNA polymerase, nucleosomes, chromatin remodelling enzymes, DNA looping, various forms of post-translational modification and DNA methylation. The directed edges correspond to transitions between microstates arising from biochemical reactions taking place on chromatin, such as the binding and unbinding of TFs or co-regulators or post-translational modification or demodification of proteins bound to DNA. Directed graphs of this kind are often found in the literature as qualitative summaries of the behaviour of regulatory mechanisms. Such cartoons can be given a rigorous mathematical basis through the methods introduced here.

The labels on the edges supply quantitative information in the form of effective rate constants for the corresponding transitions. Each label has units of inverse time, as in per second. The rate of some transitions, such as binding events, can depend on the concentration of components in solution around DNA. The labels can therefore be compound expressions involving component concentrations as well as kinetic parameters. In this way biochemical non-linearity is accommodated in the labels. An important feature of the framework is that the numerical values of the parameters do not have to be known in advance. They can be treated as symbols and many properties of the system can be calculated in symbolic form. This permits analysis without having to measure or estimate the actual values of the parameters.

The level of granularity used for the microstates, and the corresponding transitions, is a matter of choice. It can range from coarse-grained descriptions of open and closed chromatin to fine-grained descriptions of DNA sequence, individual nucleosomes and specific histone modifications. The choice depends on the context, the available experimental methods and data and the biological questions being asked. The graph constitutes a mathematical model of the system being studied and is best thought of not as a description of reality but as a precise statement of the assumptions being made about that reality – a hypothesis – from which rigorous deductions can be made and experiments proposed [43].

Because there is only one molecule of DNA, the dynamical behaviour of microstates has to be understood in terms of probabilities. If we imagine watching DNA over time, the microstates will fluctuate as transitions take place due to random molecular events, such as binding or unbinding of components. Let us denote the probability of the system being in microstate i at time t by u _i (t). The following thought experiment may help to interpret this quantity. Imagine a large number of copies of the system being created in the identical starting condition at time 0, with the same initial microstate and the same protein components present in the surrounding solution at the same concentrations. As time progresses, the randomness of molecular events will cause the different copies of the system to diverge so that different microstates will be found in each system copy. The proportion of copies in which microstate i is found at time t is an approximation for u _i (t) and this approximation becomes more accurate as the number of copies is increased. In other words, u _i (t) measures how often microstate i will be found at time t, were it possible to repeatedly replay the system from its initial condition at time 0.

Probabilities can appear difficult to reason with but the graph-based framework offers a different way to think about them which may be more familiar. The vertices of the graph are regarded as chemical species with concentrations, the edges as chemical reactions and the labels as rate constants. Each reaction has only a single substrate and only a single product, like an isomerisation, so the graph describes a kind of one-dimensional chemistry. This macroscopic interpretation allows us to reason about concentrations and reactions but gives the same results as the microscopic interpretation in terms of probabilities and transitions. In other words, if we imagine placing concentrations of matter at each vertex and allowing the chemistry to work, then the change in concentrations over time is identical to the change in probabilities over time. The only thing we have to remember is that probabilities add up to 1 – the system must be in some microstate – so that the total concentration of matter at all vertices should be kept at 1. Because the reactions only move matter between vertices, and neither create nor destroy it, the total concentration remains the same over time (see Equation 2 below), so we only need to make it 1 to begin with.

It is easy to imagine that, no matter what initial concentrations of matter are distributed over the vertices, the one-dimensional chemistry will eventually reach a steady state, in which production and consumption of each species are in balance and the concentration of each species is unchanging. Such a steady state occurs no matter what the structure of the graph. In a general graph, the steady state can depend on the initial concentrations that were chosen at time 0, so that there is a memory of these initial conditions (see the section ‘Formation of an inherently bounded chromatin domain’). However, if the graph is strongly connected, such memory is lost and the steady state becomes independent of the initial conditions and depends only on the structure of the graph. A strongly connected graph is one in which any pair of vertices are connected, both ways, by a path of consecutive edges that all point in the same direction (Figure 2A). In effect, any two vertices can communicate with each other in both directions. Strong connectivity depends only on the edges and not on the labels.

A strongly connected graph can be arbitrarily large and complicated but its one-dimensional chemistry is particularly simple. The steady-state concentration of each species can be calculated in terms of the edge labels using certain sub-graphs called spanning trees (see Equation 7 below). Among other things, this shows that each microstate in a strongly connected graph has positive probability at steady state: if such a system is watched over time, each microstate will appear at steady state, even if that microstate had zero probability in the initial condition.

A general graph, which is not strongly connected, breaks up naturally into maximal strongly connected sub-graphs, or strongly connected components (SCCs) (Figure 2B). Once matter has left a SCC under one-dimensional chemistry, it can never return to it, for otherwise the SCC would not be maximal. Hence, matter eventually accumulates on those SCCs from which there is no escape, which are the terminal SCCs. If a microstate is not in a terminal SCC, its steady-state probability is zero: if the system is watched over time, such microstates never appear in the steady state, even if they had positive probability in the initial condition. For the microstates that do lie in terminal SCCs, their steady-state probability may or may not be zero depending on the initial conditions. For instance, if matter is only placed on the vertices of one terminal SCC, it will remain there forever and cannot escape into any other SCC, whose vertices will have zero probability at all times.

A system that reaches thermodynamic equilibrium always has a strongly connected graph. The property of detailed balance, which must always hold at equilibrium, requires that each edge in the graph has a corresponding reverse edge, so that strong connectivity is guaranteed. If the labels on a pair of reversible edges are a and b, then the ratio a/b is a thermodynamic quantity that depends only on the free energy difference between the two microstates (see Equation 6 below). The steady-state probabilities depend only on these thermodynamic ratios and can be calculated as products of the ratios along paths in the graph, without the need for any spanning trees (see Equation 5 below). This gives the same result as equilibrium statistical mechanics. In this way, the framework provides a generalisation of equilibrium statistical mechanics for gene-regulation systems that are far from equilibrium.

Constructing graphs to describe gene regulation

Linear framework graphs are constructed from labelled edges, which arise from two kinds of transitions, as listed below. The main restrictive assumptions concern the interplay between mechanisms taking place in solution around chromatin and those taking place on chromatin itself. The basic approach is to assume that these can be uncoupled from each other. More relaxed assumptions can be made, using the methods of [35], but at the expense of considerably increased complexity.

Binding transitions

These represent the binding of a component L to a microstate (Figure 3A). The label is a=k[ L], where k is an on-rate and [ L] is the free concentration of L. We follow the thermodynamic formalism and assume, first, that components are neither synthesised nor degraded over the timescale of interest so that their total amounts are conserved quantities and, second, that the depletion of L can be ignored, so that the binding of a single molecule of L does not appreciably change its free concentration, [ L]. In other words, [ L]≈L _tot. Non-specific binding to DNA can significantly reduce the free concentration and if this is thought to jeopardise the no-depletion assumption, a more elaborate analysis is needed [36],[44].

Components can also engage in interactions such as oligomerisation. We again follow the thermodynamic formalism and assume that such reactions are fast compared to binding reactions on DNA, so that they have reached a rapid equilibrium. The label on the edge has the form a=k[ X], were k is an appropriate on-rate and X is the component form that binds to DNA (Figure 3B). [ X] can be calculated in terms of the concentrations of the underlying components using the rapid equilibrium assumption (Methods).

Calculating microstate probabilities at steady state

The mathematical details of the linear framework were developed in previous work [35]-[37], as reviewed in [38]. As this may not be familiar, and to keep this paper as self-contained as possible, the material is summarised here. Proofs of most of the assertions can be found in [37]. A graph of the kind constructed above, as in Figure 1, gives rise to a linear differential equation that describes how the probabilities of each microstate change in time. We first explain how this differential equation arises and then show how microstate probabilities can be calculated at steady state. The key formulas for the microstate probabilities are Equation 5 at equilibrium and Equation 7 away from equilibrium. We have italicised mathematical concepts that may be unfamiliar and have provided a glossary to explain these in the Methods.

Laplacian dynamics

Suppose we are given a graph G, as in Figure 4A, with vertices indexed 1,…,n. We typically use the index 1 for the reference microstate with no TFs bound and choose the order of the other microstates arbitrarily. The notation $i\stackrel{a}{\to }j$ signifies the edge with label a from source vertex i to target vertex j. A dynamics can be imposed on G in two equivalent ways. In the macroscopic interpretation, the vertices are chemical species and the edges are chemical reactions, which convert source species to target species. The edge labels are rate constants for the corresponding reactions, assuming mass-action kinetics. Since each reaction is uni-molecular, with only one substrate and one product, this one-dimensional chemistry yields a linear dynamics (Figure 4A),

$\frac{d}{\mathit{\text{dt}}}x\left(t\right)=\mathcal{ℒ}\left(G\right)·x\left(t\right),$

(1)

where x(t) is a column vector of species concentrations and $\mathcal{ℒ}\left(G\right)$ is an n×n matrix whose entries are labels, which is called the Laplacian matrix of G.

Since the dynamics inter-converts between species and neither creates matter nor destroys it, the total concentration does not change over time. The dynamics therefore satisfies the conservation law

$x_1\left(t\right)+\cdots +x_n\left(t\right)=u_{\text{tot}}.$

(2)

This corresponds to the columns of the Laplacian matrix adding up to 0 (Figure 4A), so that $1^t·\mathcal{ℒ}\left(G\right)=0$, where 1 signifies the all-ones column vector and ^t denotes the transpose operation, which turns a column vector into a row vector.

In the microscopic interpretation, the vertices are microstates, the edges are transitions between microstates and the labels are infinitesimal transition rates for the corresponding edges. This means that, if $i\stackrel{a}{\to }j$ and Δ t is a time interval sufficiently small so that a Δ t<1, then the probability of taking the transition from state i to state j is approximately a Δ t and the approximation gets better as Δ t gets smaller (see Equation 15 in the glossary). This interpretation defines a continuous time, finite state Markov process. A Markov process gives rise to a master equation that describes how the microstate probabilities change over time. This master equation is identical to Equation 1, so that

$\frac{d}{\mathit{\text{dt}}}u\left(t\right)=\mathcal{ℒ}\left(G\right)·u\left(t\right),$

where u _i (t) is the probability of occurrence of microstate i at time t. The only difference with the macroscopic interpretation is that probabilities must always add up to 1, so that u _tot=1 in Equation 2. Matrices of Laplacian type often arise when master equations are used but the underlying graph, from which the Laplacian can always be derived, has not been exploited as we do here.

At thermodynamic equilibrium

If the graph represents a system that can reach thermodynamic equilibrium, then detailed balance must be satisfied [36]. This requires two conditions to hold. First, the graph must be reversible: if the graph has an edge $i\stackrel{a}{\to }j$, then it must also have a reverse edge, $j\stackrel{b}{\to }i$, corresponding to the same underlying biochemical reaction working in reverse. Note that reversible edges imply that the graph is strongly connected. Second, in any steady state, x ^∗, any such pair of reversible edges must be independently at equilibrium, with the forward flux in balance with the reverse flux, irrespective of any other edges involving i and j. Setting the two fluxes to be in balance, it follows that $x_j^{\ast }=(a/b)x_i^{\ast }$.

To determine ${\rho }_j^G$, choose any path of reversible edges from vertex 1 to vertex j,

$\begin{array}{cc}1& =i_1\underset{b_1}{\overset{a_1}{\rightleftharpoons }}{i}_2\underset{b_2}{\overset{a_2}{\rightleftharpoons }}\dots {\rightleftharpoons }_{b_{p-1}}^{a_{p-1}}{i}_p{\rightleftharpoons }_{b_p}^{a_p}{i}_{p+1}=j,\end{array}$

and let ${\rho }_j^G$ to be the corresponding product of label ratios,

${\rho }_j^G=\left(\frac{a_p}{b_p}\right)\left(\frac{a_{p-1}}{b_{p-1}}\right)\dots \left(\frac{a_2}{b_2}\right)\left(\frac{a_1}{b_1}\right).$

(5)

It follows from detailed balance that $x_j^{\ast }={\rho }_j^G{x}_1^{\ast }$, so that x ^∗=λ ρ ^G where $\lambda =x_1^{\ast }$. Hence, ρ ^G provides the required basis vector of $ker\mathcal{ℒ}\left(G\right)$, from which probabilities can be calculated using Equation 4. For this procedure to be consistent, ${\rho }_j^G$ must be independent of the chosen path from 1 to j. This is ensured by the cycle condition, which is a necessary consequence of detailed balance [36]. It is an important feature of being at thermodynamic equilibrium that history does not matter: any path to a microstate can be used to determine its equilibrium probability.

Equation 5 is equivalent to the thermodynamic formalism through van’t Hoff’s formula. If $i\stackrel{a}{\to }j$ and $j\stackrel{b}{\to }i$, then, at thermodynamic equilibrium,

$\left(\frac{x_j^{\ast }}{x_i^{\ast }}\right)=\left(\frac{a}{b}\right)=exp\left(\frac{-\mathrm{\Delta G}}{\mathit{\text{RT}}}\right),$

(6)

where Δ G is the free-energy difference between microstates j and i, R is the molar Boltzmann constant and T is the absolute temperature. The product of label ratios in Equation 5 is transformed, through the exponential function in Equation 6, into a sum of free energies, which determines the free energy of microstate j relative to that of the reference microstate 1. The denominator in Equation 4 is then the partition function of equilibrium statistical mechanics.

Thermodynamic equilibrium requires detailed balance but a graph can satisfy detailed balance without being at equilibrium. For instance, certain graph structures in which each edge is reversible, such as a sequence structure (Figure 5A) or, more generally, a tree structure (Figure 5B), always satisfy detailed balance (Methods). In such a graph the edges may involve dissipative mechanisms. However, although an edge $i\stackrel{a}{\to }j$ is accompanied by a reverse edge $i\stackrel{a}{\to }j$, these edges may not arise from an underlying biochemical reaction operating reversibly but from two separate dissipative reactions, such as phosphorylation and dephosphorylation, each acting irreversibly. The ratio a/b would no longer have a thermodynamic interpretation in terms of a free energy difference, as in Equation 6.

Non-strongly connected graphs

Graphs arising in gene regulation may not always be strongly connected (see the section ‘Formation of an inherently bounded chromatin domain’ and Figure 6C). Steady-state probabilities for non-strongly connected graphs can be calculated by considering the SCCs of G (Figures 2B and 4C). The SCCs inherit connections from the underlying graph but these connections can never form a cycle, for otherwise the SCCs would collapse into each other. It is therefore possible to identify terminal SCCs, from which there are no outgoing connections. The terminal SCCs yield steady states in the following way.

Let T ₁,…,T _t denote the terminal SCCs. Each T _k is by definition strongly connected, so that it has a basis vector ${\rho }^{T_k}\in ker\mathcal{ℒ}\left(T_k\right)$, as given by Equation 7. We can now construct the vector ρ ^G,k that agrees with ${\rho }^{T_k}$ on those microstates that lie in T _k and which is zero on all other microstates (Figure 4C). The vectors ρ ^G,k provide a basis for the kernel of the Laplacian of G:

$ker\mathcal{ℒ}\left(G\right)=〈{\rho }^{G,1},\dots ,{\rho }^{G,t}〉.$

(8)

The dimension of the kernel is then t, the number of terminal SCCs. Note that, if i is any microstate that is not in a terminal SCC, then ${\rho }_i^{G,k}=0$ for each basis vector ρ ^G,k.

The t basis vectors in $ker\mathcal{ℒ}\left(G\right)$ are matched by t conservation laws. In contrast to Equation 2, which is the only conservation law when t=1, the additional conservation laws for t>1 depend on the structure of the graph. These additional laws can be algorithmically calculated from $\mathcal{ℒ}\left(G\right)$.

Any steady state x ^∗ can be expressed as a linear combination of the basis vectors in Equation 8. If these vectors are normalised to their respective totals, then, in the resulting expression for x ^∗,

$x^{\ast }=z_1\left(\frac{{\rho }^{G,1}}{1·{\rho }^{G,1}}\right)+\cdots +z_t\left(\frac{{\rho }^{G,t}}{1·{\rho }^{G,t}}\right),$

(9)

the coefficients z ₁,…,z _t are the values taken by the t conservation laws.

Regulation of steroid-hormone responsive genes

Ong et al. have put forward a theoretical framework for gene induction [46], motivated by studies of steroid-hormone receptors [51]. They use ad hoc methods, which are independent of previous work on gene regulation. We show here how their analysis can be generalised and simplified within the linear framework.

Recent work on steroid-hormone sensitive genes has revealed new co-regulators, such as the Ubiquitin conjugating enzyme, Ubc9, indicating the existence of multiple steps in addition to hormone-receptor binding to DNA [46]. Despite this additional complexity, gene-regulation functions [16], which describe how rates of gene expression depend on hormone concentration, are well fitted to Michaelis–Menten style functions, or first-order Hill dose–response curves (FHDCs) in the language of Ong et al., who use their theoretical framework to derive conditions under which such FHDCs arise.

They consider a sequence of reversible reactions (Figure 5A), representing the behaviour of the promoter of a hormone-sensitive gene. Such a sequence graph always satisfies detailed balance (Methods). We consider the more general case of an arbitrary graph G of reversible edges that satisfies detailed balance. This might be, for instance, a tree graph (Figure 5B), which also always satisfies detailed balance (Methods). If a general graph satisfies detailed balance it may not necessarily reach thermodynamic equilibrium and the edges of G may involve dissipative mechanisms.

which is a FHDC as a function of [ S].

The constants M _G and K _G have clear interpretations in terms of G. M _G is (evidently) the average rate of gene expression at saturation (i.e., when [ R S]=R _tot). Less obviously, K _G is K _R multiplied by the saturation probability of those microstates in which R is not bound. Additional file 1A gives the details of the proof and shows how the formulas in Ong et al. emerge from Equation 11. It also discusses how Ong et al. show, for the special case of a sequence, that g([ S]) remains a FHDC even if the no-depletion assumption is dropped at a concentration limiting step. Ong et al. also address other issues, such as inhibitory reactions, which are not discussed here.

The framework introduced here generalises and clarifies the work of Ong et al., showing how formulas like Equation 11 can be rigorously proved irrespective of the complexity of the underlying graph. The interpretation of the parameters in Equation 11 is new but emerges easily from our analysis (Additional file 1A). However, because detailed balance is assumed, the consequences of being away from equilibrium remain hidden, as we will see subsequently.

Formation of an inherently bounded chromatin domain

Our next application is to a model of chromatin organisation, with no explicit gene regulation. Hathaway et al. recently showed how a bounded chromatin domain could be nucleated in vivo and stably inherited as a form of epigenetic memory [47]. To explain the dynamics of such domains, they developed a mathematical model based on a linear array of 257 nucleosomes [47],[48]. This model is readily translated into our framework. We considered nucleosome arrays with varying numbers of sites n. We placed the nucleation site at the right-hand end of our array (Figure 6A). This is essentially similar to the left-hand half of the array of 2n−1 nucleosomes (for n=129) considered by Hathaway et al. The microstates correspond to array marking patterns, of which there are 2 ⁿ , while the edges correspond to mark nucleation, propagation and turnover (Figure 6A,B). Propagation and turnover were assumed uniform at all nucleosomes, at rates k+ and k_, respectively. However, nucleation was limited to the nucleation site at rate k+, so that some edges are not reversible. This irreversibility reflects the dissipative mechanism of histone marking and the non-equilibrium nature of the model. The graph does not satisfy detailed balance but is strongly connected.

Hathaway et al. used a Monte Carlo simulation to generate stochastically a succession of microstates, from which steady-state probabilities were estimated as the frequencies with which microstates appear. They found that, if k+/k_≤1.5, marking persisted in a stochastically fluctuating but inherently bounded domain near the nucleation site, reflecting what was found experimentally.

Monte Carlo simulation is an efficient method for studying very large graphs: an array of 257 nucleosome has a graph with approximately 10⁷⁷ microstates. However, the linear framework provides mathematical access to the steady-state probabilities for any array size and this yields insights that are not easily found by simulation. For instance, the ratio k+/k_ appears as a convenience in the simulations [48]. However, for a nucleosome array of n sites, the spanning trees in the corresponding graph (Figure 6A) have 2 ⁿ −1 edges, each of which is labelled k+ or k_. Dividing Equation 7 by ${\left(k\text{\_}\right)}^{2^n-1}$, it is evident that the steady-state probabilities in Equation 4 depend only on the ratio k+/k_ and not on the individual rates. The importance of the ratio becomes readily apparent within our framework.

More significantly, Hathaway et al. proposed a modification to their model to explain the inherited stability of the domain after the nucleating stimulus was removed. They imposed a stabilisation of the nucleosome mark through a transition to a hypothetical new marked state, whose turnover was inhibited (Figure 6C, left). Each nucleosome can now be in one of three states and the graph has 3 ⁿ microstates (Figure 6C, right, for n=2). Because turnover is prevented by the stabilised mark, the graph is no longer strongly connected. If nucleation is stopped, as was done in the simulation, then the resulting graph has two terminal SCCs, each consisting of a single extreme microstate, one in which the entire nucleosome array is unmarked and the other in which the entire array is stably marked. According to Equation 9, all other microstates have zero steady-state probability.

Which of the two extreme microstates is reached in a simulated trajectory depends on the microstate in which nucleation is stopped. If some nucleosome has become stably marked in that microstate, then it cannot become unmarked, so the trajectory can only reach the completely stably marked microstate. This is likely to happen once the inherently bounded domain is established, unless the stabilisation rate, k ^∗, is so low that no stable mark has appeared. In their simulation, Hathaway et al. chose k ^∗ to be low compared to propagation and turnover but not so low that stable marks had not appeared by the time nucleation was stopped. They concluded that the inherently bounded domain was stably maintained in the absence of the initial nucleating stimulus. Our analysis shows that this conclusion is incorrect. Once nucleation is stopped, the bounded domain becomes a transient phenomenon, which eventually expands to fill the whole array. It is conceivable that a bound on the domain size is maintained for sufficiently long to still be biologically relevant. But this places the stabilising rate k ^∗ in a double bind: it must be sufficiently high so as to stabilise the domain, yet sufficiently low so as not to destroy its boundedness too quickly. Such fine-tuning of rate constants is inherently fragile and we think it is more likely that other mechanisms are at work to ensure stable inheritance of the inherently bounded domain.

Our framework allows these conclusions to be reached by elementary mathematical deductions, without the need for the numerical simulations undertaken by Hathaway et al.

Regulation of yeast PHO5

We now turn back to gene regulation and to one of the very few models in which a non-equilibrium mechanism has been rigorously analysed without assuming detailed balance. Pho5 is an acid phosphatase in Saccharomyces cerevisiae that is expressed under phosphate-starvation conditions. Kim and O’Shea undertook a quantitative analysis of PHO5 regulation by the transcription factor Pho4, using a construct detached from the phosphate-response pathway [52] (Figure 7A).

To calculate the PHO5 gene-regulation function, Kim and O’Shea constructed a stochastic master equation based on a graph of transitions between DNA states. They pointed out that the nucleosomal transitions were dissipative and in some cases irreversible under their assumptions, so that detailed balance could not be assumed. Accordingly, they determined steady-state probabilities using the Symbolic Math Toolbox in MATLAB.

Kim and O’Shea’s graph of transitions is readily translated into our linear framework (Figure 7B). They assumed that the binding of Pho4 saturates according to a Hill function, which can be accommodated in a similar way to Figure 3B. The non-binding reactions correspond to unbinding of Pho4 (Figure 3C), or to nucleosomal assembly or disassembly (Figure 3F). The graph is strongly connected, a point not mentioned by Kim and O’Shea, but as noted above for Equation 7, this ensures that the steadystate probability of each microstate is positive. They assumed that PHO5 is transcribed when there is no nucleosome occluding the TATA box, so that, in the average in Equation 10, g _i =1 for the microstates 2, 3, 7, 8, 9 and 12 on the right in Figure 7B and g _i =0 for those on the left. We used our own software written in the programming language Python to enumerate the spanning trees by a fast algorithm and then used the polynomial algebra capabilities of Mathematica to calculate the microstate probabilities and the gene-regulation function (Methods). This gave an identical result to Kim and O’Shea’s MATLAB calculation (H Kim, personal communication, January 2013). This strongly suggests that what can be done for the yeast PHO5 gene can be systematically undertaken for other genes with non-equilibrium features, with the solution now being understood explicitly through Equation 7, without recourse to MATLAB.

Having calculated the gene-regulation function using our framework, we sought to compare it to the experimental data acquired by Kim and O’Shea [52]. They used their synthetic construct (Figure 7A, with details in the caption) to measure the PHO5 gene-regulation function. In response to doxycycline, individual cells expressed Pho4-YFP, which was treated as the input to the gene-regulation function, and this induced the expression of CFP from the Pho4-responsive promoter in the construct. CFP was treated as the output as a proxy for Pho5. By using different doses of doxycycline to cover a range of Pho4-YFP expression levels, the gene-regulation function was assembled from single-cell measurements. Kim and O’Shea also measured the gene-regulation function of five other variant promoters, in which the low-affinity and high-affinity sites for Pho4 binding were either exchanged or removed.

Kim and O’Shea estimated the threshold and the maximum expression level of each variant by fitting their experimental data to a Hill function, whose Hill coefficient was found to be nearly 2 for all variants. They then fitted the estimated threshold and maximum values to the calculated gene-regulation function for each variant and found good agreement ([52], Figure 5). We were curious as to how well the gene-regulation function itself would fit the data. This is a more challenging question because the data are noisy and the gene-regulation function is very complicated (see below). To address this, we first smoothed the data. We then used numerical optimisation to find excellent quantitative fits to each variant individually (Figure 8, red curves) but could only undertake a manual fit to all variants collectively, which yielded the parameter values in Equation 16 (Methods). The collective fit was considerably poorer (Figure 8, black curves). While this broadly confirms Kim and O’Shea’s more coarse-grained analysis, it also suggests that the individual variants may exhibit more nuanced behaviours, which are better described by distinct parameter values.

History-dependent complexity away from equilibrium

The simplicity of these highest-order terms is deceptive, however. The numerator of Equation 12 has 261 distinct monomials while the denominator has 500 distinct monomials. Indeed, the graph in Figure 7B has 53,376 spanning trees in total. We see that the calculated PHO5 gene-regulation function is very complicated – the full details shown in Additional file 1C cover six pages – despite the model having only two binding sites and two nucleosomes. Because Kim and O’Shea did not provide the gene-regulation function in their original paper, these features are revealed here for the first time.

The linear framework allows us to understand this surprising explosion in complexity. At equilibrium, Equation 5 shows that any single path to a microstate can be used to calculate its steady-state probability. As a physicist would say, free energy at equilibrium is a function of the microstate, not of the route through which that microstate is reached. In marked contrast, away from equilibrium, Equation 7 shows that every spanning tree rooted at that microstate is required. In this case, all routes to the microstate become relevant and microstate probabilities depend in a more intricate way on the structure of the graph. Equation 7 takes care of the bookkeeping. The number of spanning trees increases very rapidly with the size of a graph: the complete undirected graph on n vertices (i.e., the graph in which there is an undirected edge between each pair of distinct vertices) has n ⁿ⁻² spanning trees in total. This worse than exponential increase manifests itself in the complexity of the PHO5 gene-regulation function.

It is important to appreciate, however, that it is not the complexity or the size of a graph that is the dominant factor in explaining the complexity found here. If we imposed additional edges on the graph in Figure 7B so as to make all the edges reversible, this would only make the graph more complex. If we then imposed detailed balance, which restricts the values of the parameters, the equilibrium probabilities would be given by Equation 5 rather than Equation 7 and the gene-regulation function could be written down in a few lines. The complexity uncovered here depends crucially on being far from thermodynamic equilibrium.

Additional study of PHO5 has shown that nucleosomes decouple the threshold for PHO5 expression from its dynamic range [53]. However, this kind of behaviour can be recapitulated within the thermodynamic formalism [54]. This suggests that the full implications of non-equilibrium behaviour, as revealed by the complexity of the PHO5 gene-regulation function, have not yet been uncovered experimentally. To suggest experimental options, we need ways to decompose the complexity found in Additional file 1C and to attribute aspects of it to specific biochemical mechanisms. Approximation methods may help in particular cases [55] but new ideas are needed for addressing the complexity barrier systematically, to which we now turn.

Graph independence leads to reduced complexity

Gene regulation often takes a modular form, with repeated binding sites, reiterated motifs and multiple enhancers [56],[57]. The microstate probabilities and the resulting gene-regulation function could become extremely complicated, especially if the modules are operating far from equilibrium. There is, however, one context in which simplification may be expected. This occurs when modules operate independently of each other, so that whatever takes place within one module does not affect what takes place in any other module. For instance, developmental genes are often regulated by multiple enhancers, which sometimes appear to act independently of each other [58].

where a ₁,…,a _m can be arbitrary numbers. Evidently, the partition function in Equation 13 is considerably less complex and easier to understand. In the light of this result for equilibrium systems, we wanted to find a generalisation in which the modules are no longer individual binding sites but are represented by potentially complex graphs, which may not be at thermodynamic equilibrium. Such modules might correspond, for instance, to independent enhancers.

We used the product graph construction to capture the concept of independence. Let G and H be any two graphs which represent two modules within a gene regulation system. We make no assumptions about the graphs, which do not have to be at equilibrium and do not have to be strongly connected. The product graph G×H is constructed as follows (Figure 9). It has vertices (i,j), where i is a vertex in G and j is a vertex in H. The vertices are enumerated lexicographically, so that (i,j)<(i ^′,j ^′) if either i<i ^′ or i=i ^′ and j<j ^′. For each labelled edge $i_1\stackrel{a}{\to }{i}_2$ in G and for every vertex j in H, the labelled edge $(i_1,j)\stackrel{a}{\to }(i_2,j)$ is created in G×H. The retention of the same label a on these edges ensures that the transition from (i ₁,j) to (i ₂,j) occurs independently of j and always at the same rate, which captures the independence assumption. Similarly, for each labelled edge $j_1\stackrel{a}{\to }{j}_2$ in H and for every vertex i in G, the labelled edge $(i,j_1)\stackrel{b}{\to }(i,j_2)$ is created in G×H. These are the only edges in G×H.

This uses the Kronecker product of two vectors, x⊗y, defined by (x⊗y)_(i,j)=x _i y _j (Figure 9). If either G or H are not strongly connected then G×H will not be strongly connected. A basis for the Laplacian kernel of G×H is then given by the Kronecker products ρ ^G,i⊗ρ ^H,j between each pair of basis vectors from each respective kernel. The precise product theorem is stated and proved in Additional file 1B.

which can be divided out to give the much simpler expressions in Figure 9. The basis vectors from the product theorem are substantially less complicated, both in degree and in the numbers of monomials, than those from Equation 7.

This product theorem is important because it shows that a system that is far from equilibrium may still have simple expressions for its microstate probabilities. What is required is that the system has independent modules within it. This suggests a starting point for addressing the complexity challenge identified above, as reviewed further in the Discussion below.

Calculating labelling functions

which gives the formula for Φ([ L]) shown in Figure 3B. Rapid equilibrium amounts to a timescale separation, which uncouples the dynamics of the interactions in solution from those on DNA. The rapid equilibrium equations for more complicated interactions can often be formulated in terms of the linear framework, which can then be used to calculate [ X].

Glossary of mathematical concepts

Markov process. A time-varying probability distribution over a set of states in which the probability of reaching a given state in the next time step depends only on the current state. If time varies continuously then the next time step is interpreted infinitesimally, by taking a small unit of time, Δ t, and letting this tend to zero. The Markov property says that history does not matter in making the choice of which state comes next in time. However, history may be essential for determining the steady-state probabilities, as happens when the system is far from thermodynamic equilibrium.

With this notation, the probability of occurrence of microstate i at time t, which was denoted u _i (t) in the main text, is given by u _i (t)=Pr(X(t)=i).

Master equation. The probability of being in microstate i at time t+Δ t, u _i (t+Δ t), can be calculated in terms of u _j (t) and the infinitesimal transition rate from j to i, taking into account all microstates j that have an edge to i. The resulting differential equation, obtained by letting Δ t→0, which describes the forward evolution of probabilities over time, is the master equation, or Kolmogorov forward equation, of the Markov process [68]. The equivalence between the master equation of X(t) and Laplacian dynamics is proved in ([37], Corollary 2).

Kernel. If M is an n×n matrix acting on column vectors of size n, then the kernel of M, kerM, is the subspace of column vectors that become zero when multiplied by M: kerM={v | M·v=0}.

Vertex i is said to be strongly connected to j if i⇝j and j⇝i. Strong connectivity is an equivalence relation on the vertices and the equivalence classes are called the SCCs of G. A graph is strongly connected if it has only one SCC. The graph in Figure 4B is strongly connected.

Cycle condition. If a graph describes a system that can reach thermodynamic equilibrium then it must satisfy detailed balance, as described in the main text. If detailed balance holds, then, in any cycle of reversible edges, the product of the labels going clockwise around the cycle must equal the product of the labels going counterclockwise around the cycle. Conversely, if a graph has reversible edges and the cycle conditions holds, then detailed balance is satisfied for any steady state of the graph. This is proved in ([36], Supporting Information).

Sequence/tree of reversible edges. A graph consisting of reversible edges, which are arranged in a sequence (Figure 5A) or, more generally, in a tree structure (Figure 5B), automatically satisfies detailed balance, irrespective of the edge labels. The argument for a sequence was presented in [69] but is easily generalised to a tree. Given a reversible edge, $i\stackrel{a}{\to }j$ and $j\stackrel{b}{\to }i$, and a steady state x ^∗, the net flux through the reversible edge is $a{x}_i^{\ast }-b{x}_j^{\ast }$. If the reversible edge is a leaf of the tree structure then there can be no net flux leaving the tree from that edge. Hence, $x_i^{\ast }=(b/a)x_j^{\ast }$. This reversible edge is therefore at equilibrium. This holds irrespective of the labels a and b. Arguing in this way by induction from the leaves, each reversible edge in the tree is independently at equilibrium, so that detailed balance holds.

Rooted spanning trees. A spanning tree of a graph G is a sub-graph that contains each vertex of G (spanning) and that has no cycles when edge directions are ignored (tree). A spanning tree is rooted at vertex j in G if j is the only vertex with no outgoing edges. A graph is strongly connected if, and only if, it has at least one rooted spanning tree at each vertex ([37], Lemma 1). Figure 4B shows a strongly connected graph, together with the spanning trees rooted at each vertex.

Terminal strongly connected components. Let [ j] denote the SCC of G containing vertex j. In other words, [ j] is the equivalence class of vertex j under the relation of strong connectivity, as defined above. The SCC [ i] is said to precede [ j], denoted [ i]≼ [ j], if either [ i]= [ j] or some vertex in [ i] ultimately reaches some vertex in [ j]: i ^′⇝j ^′ where i ^′∈ [ i] and j ^′∈ [ j]. Precedence defines a partial order on the SCCs of the graph G. We can therefore speak of the terminal SCCs, which are those that do not precede any other SCC. The graph in Figure 4C has three SCCs of which two are terminal (asterisks), while the graph in Figure 6C has five SCCs of which two are terminal (asterisks).

Calculating the PHO5gene-regulation function

The gene-regulation function of the PHO5 example was calculated using the matrix-tree formula in Equation 7 and is shown in full in Additional file 1C. Software for enumerating spanning trees is available in packages like MATLAB, Mathematica and Maple, but we found these to be incapable of dealing with the large number of trees that arise. We therefore implemented in Python the fast algorithm developed by Takeaki Uno [70]. The resulting program reads a text file containing a description of a graph as a collection of labelled edges and, for each vertex in the graph, writes a text file listing the spanning trees rooted at that vertex. We also implemented an accompanying Mathematica notebook, which reads the graph description and the spanning tree files and assembles each ${\rho }_i^G$ as a polynomial function of the edge labels. The gene-regulation function can then be calculated using standard Mathematica functions for manipulating polynomial expressions. The Python program and the Mathematica notebook are freely available from our web site [71].

Fitting to the experimental data of Kim and O’Shea

Kim and O’Shea constructed 12 promoter variants ([52], Figure 3a). Six of these variants place a high affinity (H), low affinity (L) or deleted (X) Pho4-binding site in the positions corresponding to UASp1 and UASp2 in Figure 7A. The remaining six variants use sites occluded by nucleosome -3, which is not modelled in Figure 7, and we did not analyse these variants. The wild-type promoter in Figure 7 corresponds to variant LH.

We obtained the experimental data in the form of an Excel spreadsheet [72]. This gives the raw fluorescence values for YFP, CFP and RFP (yellow, cyan and red fluorescent proteins, respectively) for about 400 to 500 cells for each variant under different doxycycline concentrations. The RFP was attached to a chromatin protein to mark the nucleus and the RFP value was used to normalise the YFP and CFP values on a per-cell basis to control against imaging variations. We used a ±7 moving average to smooth the data and scaled each variant to its maximum expression level for the plots shown in Figure 8.

Each of the six variants gives rise to a graph, which uses the same labels as the wild type (Figure 7B). The labels b and c are the rates of Pho4 dissociation from the low-affinity and high-affinity sites, respectively. Kim and O’Shea assumed that the Pho4 association rate, a, is the same for both sites. If the Pho4 binding sites are changed in a variant, the labels b and c occur on different edges of the wild-type graph, while if a Pho4 binding site is deleted, some vertices become inaccessible and the graph changes from the 12-vertex wild-type graph to a graph with eight vertices. We used the wild-type 12-vertex gene-regulation function and a new eight-vertex gene-regulation function calculated using Equation 7. We then changed the labels b and c in these two gene-regulation functions, as required, to generate the gene-regulation function for each of the six variants (details in the accompanying Mathematica notebook).

The Mathematica notebook in which these calculations were undertaken is freely available from our web site [71]. It provides the normalised experimental data, the smoothed experimental data and the individual and collective fits of the variant gene-regulation functions to the corresponding data.

Imposing equilibrium on the Hodges–Crabtree model

As explained in the main text, to impose equilibrium is to require that detailed balance holds. This means, first, that all edges in the graph must be reversible and, second, that the cycle condition (described in the glossary above) is satisfied. The graph of microstates for an array of three nucleosomes is shown in Figure 6B and we follow the notation introduced there in which microstates are denoted by bit strings, indicating whether (bit = 1) or not (bit = 0) a nucleosome is marked. Edges only occur between microstates that differ by a single bit, corresponding to nucleation or mark propagation, when the number of bits increases by 1 and the edge has label k+, or to mark turnover, when the number of bits decreases by 1 and the edge has label k_ (Figure 6A). Irreversibility only arises for some of the latter edges, when an isolated site, whose immediate neighbours are unmarked, loses its mark (for instance, 5→1, 3→1 and 6→2 in Figure 6B).

In traversing this path, if an edge increases the number of bits in the microstate by 1, then the label encountered must be k+, while if an edge decreases the number of bits by 1, then the label must be k_. Since the path is a cycle, the number of edges with label k+ must equal the number of edges with label k_. Furthermore, for each edge with label k+, respectively, k_, the reverse edge has label k_, respectively, k+. But then the product of the labels going clockwise around the cycle must equal the product of the labels going counterclockwise around the cycle and the cycle condition is satisfied. The graph therefore satisfies detailed balance in any steady state.

We see from this that the probability of a microstate depends only on the number of bits that are marked, rather than which bits are marked and, consequently, there can be no inherent bound on the size of the marked domain.