Our approach uses the linear framework for timescale separation (Gunawardena, 2012), details of which are provided in the 'Materials and methods' along with further references. We briefly outline the approach here.

In the linear framework, a suitable biochemical system is described by a finite directed graph with labelled edges. In our context, graph vertices represent microstates of the target molecule and graph edges represent transitions between microstates, for which the edge labels are the instantaneous transition rates. A linear framework graph specifies a finite-state, continuous-time Markov process, and any reasonable such Markov process can be described by such a graph. We will be concerned with the probabilities of microstates at steady state. These probabilities can be interpreted in two ways, which reflect the ensemble and single-molecule viewpoints of Figure 2. From the ensemble perspective, the probability is the proportion of target molecules which are in the specified microstate, once the molecular population has reached steady state, considered in the limit of an infinite population. From the single-molecule perspective, the probability is the proportion of time spent in the specified microstate, in the limit of infinite time. The equivalence of these definitions comes from the ergodic theorem for Markov processes (Stroock, 2014). These different interpretations may be helpful when dealing with different biological contexts: a population of haemoglobin molecules may be considered from the ensemble viewpoint, while an individual gene may be considered from the single-molecule viewpoint. As far as the determination of probabilities is concerned, the two viewpoints are equivalent.

The graph representation may also be seen as a discrete approximation of a continuous energy landscape, as in Figure 3, in which the target molecule is moving deterministically on a high-dimensional landscape in response to a potential, while being buffeted stochastically through interactions with the surrounding thermal bath (Frauenfelder et al., 1991). In mathematics, this approximation goes back to the work of Wentzell and Freidlin on large deviation theory for stochastic differential equations in the low noise limit (Ventsel' and Freidlin, 1970; Freidlin and Wentzell, 2012). It has been exploited more recently to sample energy landscapes in chemical physics (Wales, 2006) and in the form of Markov State Models arising from molecular dynamics simulations (Noé and Fischer, 2008; Sengupta and Strodel, 2018). In this approximation, the vertices correspond to the minima of the free energy up to some energy cut-off, the edges correspond to appropriate limiting barrier crossings and the labels correspond to transition rates over the barrier.

Figure 3

Download asset Open asset

The linear framework graph, or the accompanying Markov process, describes the time-dependent behaviour of the system. Our concern in the present paper is with systems which have reached a steady state of thermodynamic equilibrium, so that detailed balance, or microscopic reversibility, is satisfied. The assumption of thermodynamic equilibrium has been standard since allostery was introduced (Koshland et al., 1966; Monod et al., 1965) but has significant implications, as pointed out in the Introduction, and we will return to this issue in the Discussion. At thermodynamic equilibrium, we can dispense with dynamical information and work with what we call ‘equilibrium graphs’ (Figure 3). These are also directed graphs with labelled edges but the edge labels no longer contain dynamical information in the form of rates but rather ratios of forward to reverse rates. These ratios are determined by the minima of the free-energy landscape, with the equilibrium label on the edge from vertex $i$ to vertex $j$ being given by the formula in Figure 3 . Free energy is often expressed relative to a reference level, as we will do below, so it will be convenient to write the equilibrium label from $i$ to $j$ as

where ${\displaystyle \mathrm{\Delta }{\mathrm{\Phi }}_u}$ is the relative free-energy of vertex ${\displaystyle u}$, $k_B$ is Boltzmann’s constant and $T$ is the absolute temperature (Figure 3). Note that if the edge in question involves components from outside the graph itself, such as a ligand which binds to $i$ to yield $j$, then the chemical potential of the ligand will contribute to the free energy. This contribution will manifest itself in the presence of a ligand concentration term in the edge label, as seen in Figure 4. The equilibrium edge labels are the only parameters needed at thermodynamic equilibrium and the free energies of the vertices can be recovered from them, up to an additive constant. From now on, in the main text, when we say ‘graph’, we will mean ‘equilibrium graph’.

Figure 4

Download asset Open asset

We explain such graphs using our main example. Figure 4 shows the graph, $A$, for an allosteric ensemble, with multiple conformations $c_1,\mathrm{\cdots },c_N$ and multiple sites, $1,\mathrm{\cdots },n$, for binding of a single ligand ($n=3$ in the example). The graph vertices represent abstract conformations with patterns of ligand binding, denoted $(c_k,S)$, where the index $k$ designates the conformation with $1\le k\le N$, and $S\subseteq \{1,\mathrm{\cdots },n\}$ is the subset of bound sites. Directed edges represent transitions arising either from binding without change of conformation (‘vertical’ edges), $(c_k,S)\to (c_k,S\cup \{i\})$ where $i\notin S$, which occur for all conformations c_k, or from conformational change without binding (‘horizontal’ edges), $(c_k,S)\to (c_j,S)$ where $k\ne j$, which occur for all binding subsets $S$. Edges are shown in only one direction for clarity—when binding or unbinding is present, we use the direction of binding—but edges are always reversible, in accordance with thermodynamic equilibrium. Ignoring labels and thinking only in terms of vertices and edges, or ‘structure’, $A$ has a product form: the vertical subgraphs, $A^{c_k}$, consisting of those vertices with conformation c_k and all edges between them, all have the same structure and the horizontal subgraphs, $A_S$, consisting of those vertices with binding subset $S$ and all edges between them, also all have the same structure (Figure 4). Structurally speaking, we can think of $A$ as the graph product (Ahsendorf et al., 2014) of the vertical subgraph $A^{c_1}$ and the horizontal subgraph $A_{\mathrm{\varnothing }}$ (Figure 4).

In an allostery graph, ‘conformation’ is meant abstractly as any state for which binding association constants can be defined. It does not imply any particular atomic configuration of a target molecule nor make any commitments as to how the pattern of binding changes.

The product-form structure of the allostery graph reflects the ‘conformational selection’ viewpoint of MWC, in which conformations exist prior to ligand binding, rather than the ‘induced fit’ viewpoint of KNF, in which binding can induce new conformations. Considerable evidence now exists for conformational selection, in the form of transient, rarely populated conformations which exist prior to binding (Tzeng and Kalodimos, 2011). Induced fit may be incorporated within our graph-based approach by treating new conformations as always present but at extremely low probability. One of the original justifications for induced fit was that it enabled negative cooperativities, in contrast to conformational selection (Koshland and Hamadani, 2002), but we will show below that induced fit is not necessary for this and that negative HOCs arise naturally in our approach. Accordingly, the product-form structure of our allostery graphs is both convenient and powerful.

The edge labels are the non-dimensional ratios of the forward transition rate to the reverse transition rate; accordingly, the label for the reverse edge is the reciprocal of the label for the forward edge (Materials and methods). Labels may include the influence of components outside the graph, such as a binding ligand. For instance, the label for the binding edge $(c_k,S)\to (c_k,S\cup \{i\})$ is $x{K}_{c_k,i,S}$, where $x$ is the ligand concentration and $K_{c_k,i,S}$ is the association constant (Figure 1A), with dimensions of (concentration)⁻¹, as described in the Introduction. Horizontal edge labels are not individually annotated and need only be specified for the horizontal subgraph of empty conformations, $A_{\mathrm{\varnothing }}$, since all other labels are determined by detailed balance (Materials and methods).

The graph structure allows HOCs between binding events to be calculated, as suggested in the Introduction. We will define this first for the ‘intrinsic’ HOCs which arise in a given conformation and explain in the next section how ‘effective’ HOCs are defined for the ensemble. In conformation c_k, the intrinsic HOC for binding to site $i$, given that the sites in $S$ are already bound, denoted ${\omega }_{c_k,i,S}$, is defined by normalising the corresponding association constant to that for binding to site $i$ when nothing else is bound (Estrada et al., 2016),

HOCs are non-dimensional quantities. If $S$ has only a single site, say $S=\{j\}$, then the intrinsic HOC of order 1, ${\omega }_{c_k,i,\{j\}}$, is the classical pairwise cooperativity between sites $i$ and $j$. There is positive or negative intrinsic HOC if ${\omega }_{c_k,i,S}>1$ or ${\omega }_{c_k,i,S}<1$, respectively, and independence if ${\omega }_{c_k,i,S}=1$ (Figure 1A).

For any graph $G$, the steady-state probabilities of the vertices can be calculated from the edge labels. For each vertex, $v$, in $G$, the probability, ${\text{Pr}}_v(G)$, is proportional to the quantity, ${\mu }_v(G)$, obtained by multiplying the edge labels along any directed path of edges from a fixed reference vertex to $v$. It is a consequence of detailed balance that ${\mu }_v(G)$ does not depend on the choice of path in $G$. This implies algebraic relationships among the edge labels. These can be fully determined from $G$ and independent sets of parameters can be chosen (Materials and methods). For the allostery graph, a convenient choice vertically is those association constants $K_{c_k,i,S}$ with $i$ less than all the sites in $S$, denoted $i<S$; horizontal choices are discussed in the Materials and methods but are not needed for the main text.

Since probabilities must add up to 1, it follows that

Equation 3 yields the same result as equilibrium statistical mechanics, with the denominator being the partition function for the thermodynamic grand canonical ensemble. Equilibrium statistical mechanics typically focusses only on vertices and uses their free energies as the fundamental parameters. Directed graphs of the form considered here were previously used in Hill, 1966 and Schnakenberg, 1976 to study systems away from thermodynamic equilibrium, where the graph edges become essential to represent entropy production (Wong and Gunawardena, 2020). We find that the graph remains just as useful at thermodynamic equilibrium because binding and unbinding are the fundamental mechanisms through which information is integrated and these mechanisms must be represented by graph edges. Indeed, as the next section shows, graphs are invaluable for formulating higher-order concepts.

Our specification of an allostery graph allows for arbitrary conformational complexity and arbitrary interacting ligands (we consider only one ligand here for simplicity), with the independent association constants in each conformation being arbitrary and with arbitrary changes in these parameters between conformations. Moreover, the abstract nature of ‘conformation’, as described above, permits substantial generality. Allostery graphs can be formulated to encompass the two conformations of MWC (Marzen et al., 2013), nested models (Robert et al., 1987), the fluctuations of Cooper and Dryden, 1984 and more recent views of dynamical allostery (Tzeng and Kalodimos, 2011), the multiple domains of the Ensemble Allosteric Model developed by Hilser and colleagues (Hilser et al., 2012) and applied also to intrinsically disordered proteins (Motlagh et al., 2012), other ensemble models (LeVine and Weinstein, 2015; Tsai and Nussinov, 2014) and Markov State Models arising from molecular dynamics simulations (Noé and Fischer, 2008).

As mentioned in the Introduction, a systematic approach to higher-order effects using mutant-cycle analysis was developed in Horovitz and Fersht, 1990 and Horovitz and Fersht, 1992 and widely used subsequently (Carter, 2017). The HOCs presented above were introduced in our previous work (Estrada et al., 2016), and the present paper is concerned not with HOCs per se, but with effective HOCs that arise from an allosteric ensemble, as will be described below. Nevertheless, it may still be helpful to explain the relationship between our HOCs and the higher-order couplings arising from mutant-cycle analysis. We are grateful to an anonymous reviewer for making this point to us. The material which follows may be of particular interest to those familiar with the relevant literature but is not required for the main results of the paper.

Both HOCs and higher-order couplings can be seen as different ways of analysing the underlying free-energy landscape. Both approaches make essential use of directed graphs to organise this landscape. Figure 5A shows the labelled equilibrium graph for ligand binding to three sites in a single conformation, while Figure 5B shows a directed graph of the kind used in Horovitz and Fersht, 1990 for defining higher-order couplings for perturbations to three sites. The latter graphs are sometimes called ‘boxes’ (Horovitz and Fersht, 1990). We use ‘sites’ here for either individual residues or the modules described in Carter, 2017. Perturbations are typically mutations, such as replacement of an asparagine residue by alanine. The choice of replacement can make a difference to the results, but this is not usually depicted in graph representations like Figure 5B. The directed edges have rather different interpretations in the two examples in Figure 5: for the equilibrium graph in Figure 5A, a directed edge represents the biochemical process of ligand binding; for the coupling graph in Figure 5B, a directed edge represents an experimental perturbation. In both cases, the vertices have an associated free energy, denoted $\mathrm{\Delta }{\mathrm{\Phi }}_S$, where $S\subseteq \{1,\mathrm{\cdots },n\}$ is either the subset of bound sites in the equilibrium graph (Figure 5A) or the subset of perturbed sites in the coupling graph (Figure 5B). The $\mathrm{\Delta }$ notation is conventionally used in the literature to signify a free-energy difference (Equation 1) or free energy relative to a chosen zero level. A frequent choice of zero is the free energy of empty binding or of the unperturbed state, in which case $\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{\varnothing }}=0$, but we have not assumed this here. Note that the free energies of the equilibrium graph have a contribution from the ligand, which manifests itself in the dependence of the edge labels on the ligand concentration, $x$, while the free energies of the coupling graph do not. Despite this difference, the free energies provide in both cases the fundamental independent thermodynamic parameters, of which there are $2^n-1$ for $n$ sites, in terms of which both HOCs and higher-order couplings can be rigorously defined.

Figure 5

Download asset Open asset

The definition is easiest for HOCs. Equation 1 tells us that the edge label, $x{K}_{i,S}$, is given by

We omit the single conformation from subscripts for clarity. It follows from Equation 2 that HOCs can be written in terms of free energies as follows:

HOCs are non-dimensional quantities associated to graph edges. As noted above, there are algebraic relationships among them arising from detailed balance at thermodynamic equilibrium. An independent set of parameters is formed by restricting to those for which $i<S$, of which there are $2^n-n-1$. Taken together with the $n$ ‘bare’ association constants for initial ligand binding, $K_{i,\mathrm{\varnothing }}$, they form a complete set of $2^n-1$ independent parameters for the free-energy landscape. It follows from Equations 4 and 5 that these parameters can be used to recover the fundamental free energies, so that the two sets of parameters are mathematically equivalent.

Mutant-cycle studies often refer to both Horovitz and Fersht, 1990 and Horovitz and Fersht, 1992, which present apparently different measures of higher-order coupling. The second of these papers introduces what we will refer to as the ‘residual free energy’ of a vertex and denote $\mathrm{\Delta }{\varphi }_S$. This is the free energy remaining at vertex $S$ after accounting for the contributions from all proper subsets of $S$. The residual free energy may be concisely defined recursively, starting from $\mathrm{\Delta }{\varphi }_{\mathrm{\varnothing }}=\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{\varnothing }}$, by

We see from Equation 6 that $\mathrm{\Delta }{\varphi }_{\{i\}}=\mathrm{\Delta }{\mathrm{\Phi }}_{\{i\}}-\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{\varnothing }}$ and that ${\displaystyle \mathrm{\Delta }{\varphi }_{\{i,j\}}=\mathrm{\Delta }{\mathrm{\Phi }}_{\{i,j\}}-(\mathrm{\Delta }{\mathrm{\Phi }}_{\{i\}}+\mathrm{\Delta }{\mathrm{\Phi }}_{\{j\}})+\mathrm{\Delta }{\mathrm{\Phi }}_{\mathrm{\varnothing }}}$. ${\displaystyle \mathrm{\Delta }{\varphi }_S}$ may be calculated directly from $\mathrm{\Delta }{\mathrm{\Phi }}_X$ but, as the previous example suggests, overlapping contributions of the actual free energies must be cancelled out (Horovitz and Fersht, 1992, Equation 4),

To see why Equation 7 is a consequence of Equation 6, note first that Equation 7 gives the correct result for $S=\mathrm{\varnothing }$. It may then be recursively checked by assuming it holds for $X\subset S$ and substituting into Equation 6 to check that it holds for $S$. Each subset $Y\subset S$ contributes a term $\pm \mathrm{\Delta }{\mathrm{\Phi }}_Y$ arising from $\mathrm{\Delta }{\varphi }_X$ for each $X$ that satisfies $Y\subseteq X\subset S$. The sign of $\mathrm{\Delta }{\mathrm{\Phi }}_Y$ coming from Equation 7 is ${(-1)}^{\mathrm{\#}(X)-\mathrm{\#}(Y)}$. These terms almost completely cancel each other out because, letting $p=\mathrm{\#}(S)-\mathrm{\#}(Y)$,

Taking into account the additional sign coming from Equation 6, we recover Equation 7 for $S$. This proves recursively that Equation 7 is the solution of Equation 6 in terms of free energies.

We can go further to show how $\mathrm{\Delta }{\varphi }_S$ is expressed in terms of HOCs. For this, we must assume that $q=\mathrm{\#}(S)>1$. When $q=1$, ligand binding contributes to $\mathrm{\Delta }{\varphi }_S$, but when $q>1$ that is no longer the case, as we will see. Choose any site $i\in S$. The summation in Equation 7 involves $2^q$ terms $\mathrm{\Delta }{\mathrm{\Phi }}_Y$. It can be reorganised into a sum of $2^{q-1}$ terms of the form $\pm (\mathrm{\Delta }{\mathrm{\Phi }}_{Z\cup \{i\}}-\mathrm{\Delta }{\mathrm{\Phi }}_Z)$, where $Z\subseteq S\backslash \{i\}$. The sign of these terms is given by the sign of $\mathrm{\Delta }{\mathrm{\Phi }}_{Z\cup \{i\}}$ coming from Equation 7 and is therefore ${(-1)}^{\mathrm{\#}(S)-\mathrm{\#}(Z)-1}$. It is easy to see that, because $q>1$, there must be equal numbers of +1 and −1 signs. It follows from Equation 4 that

where the double exponent just means that the right-hand side is a ratio in which those terms for which $\mathrm{\#}(S)-\mathrm{\#}(Z)$ is odd go in the numerator and those terms for which $\mathrm{\#}(S)-\mathrm{\#}(Z)$ is even go in the denominator. Using Equation 2, we can rewrite $K_{i,Z}$ as ${\displaystyle K_{i,\mathrm{\varnothing }}{\omega }_{i,Z}}$. Since there are equal numbers of each sign, we can cancel each occurrence of $x{K}_{i,\mathrm{\varnothing }}$ between numerator and denominator to yield a formula for residual free energies in terms of HOCs when $\mathrm{\#}(S)>1$:

The choice of $i\in S$ in Equation 8 is arbitrary. As an illustration of Equation 8, recalling from Equation 5 that ${\omega }_{i,\mathrm{\varnothing }}=1$, we see that

Equations 8 and 9 show how the residual free energy is built up from binding at any given site to the hierarchy of subsets of the remaining sites.

Residual free energies can be thought of as a measure of collective synergy between sites (Horovitz and Fersht, 1992). They are associated to graph vertices and constitute $2^n-1$ independent parameters, with no algebraic relationships between them. It follows from Equations 6 and 7 that they are mathematically equivalent to the fundamental free energies. Residual free energies have also been independently described for other purposes in Equation 4 of Martini, 2017.

The higher-order couplings introduced in Horovitz and Fersht, 1990 appear at first sight to be quite different from the residual free energies introduced in Horovitz and Fersht, 1992. The couplings are described by examples for low orders, as are typically encountered in practice (Horovitz and Fersht, 1990). We provide a general definition here by introducing a slightly more complex version. A coupling is associated to a pair, consisting of, first, a vertex, $Z\subseteq \{1,\mathrm{\cdots },n\}$, and, second, an ordered sequence of distinct sites, $(i_1,\mathrm{\cdots },i_k)$, none of which are in $Z$, so that $Z\cap \{i_1,\mathrm{\cdots },i_k\}=\mathrm{\varnothing }$. The vertex $Z$ should be thought of as an ‘offset’ within the coupling graph and the sites, $i_1,\mathrm{\cdots },i_k$ as specifying an ordered sequence of perturbations undertaken around $Z$. Higher-order couplings are conventionally used in the literature only for $Z=\mathrm{\varnothing }$, but this more complex version is needed for the definition in Equation 11 below. Associated to such a pair $Z,(i_1,\mathrm{\cdots },i_k)$ is a $k$th order coupling, which we will denote by ${\mathrm{\Delta }}^k{\gamma }_{Z,(i_1,\mathrm{\cdots },i_k)}$. We start by defining the first-order coupling, ${\mathrm{\Delta }}^1{\gamma }_{Z,(i_1)}$, for any $Z$ satisfying the restriction above, in terms of the free energy,

With that in hand, we can define for $k\ge 2$, again for any $Z$ satisfying the restriction

where it is clear that $Z\cup \{i_k\}$ must be disjoint from $\{i_1,\mathrm{\cdots },i_{k-1}\}$, so that the right-hand side of Equation 11 is recursively well defined. Unravelling Equations 11 and 10, we see that

which corresponds when $Z=\mathrm{\varnothing }$ to Equation 1 of Horovitz and Fersht, 1990. With some more work, it can be seen that Equation 11 reproduces the $k=3$ and $k=4$ examples in Horovitz and Fersht, 1990. Equation 12 expresses the intuition behind higher-order coupling, that it measures the effect of a perturbation relative to the unperturbed state, hierarchically for a sequence of perturbations.

It can be seen quite easily from Equations 5 and 12 that

We note from Equation 13 that ‘order’ is counted differently between HOCs and conventional higher-order couplings: when $Z=\mathrm{\varnothing }$, Equation 13 relates a higher-order coupling with $k=2$ to a HOC of order 1. Substituting Equation 13 into Equation 11 and continuing the recursion, we find that

at which point the similarity with Equation 9 becomes evident and the pattern emerges. It can be shown by direct substitution in Equation 11 that the following general formula holds, which expresses higher-order couplings in terms of HOCs for any $k\ge 2$: 

Comparing Equation 14 with Equation 8 we see that, despite their very different definitions in Equations 11 and 6, conventional higher-order couplings are the same as residual free energies. Indeed, for $k\ge 1$,

Equation 15 may seem strange because a higher-order coupling is defined in terms of an ordered sequence of perturbations, $(i_1,\mathrm{\cdots },i_k)$, while a residual free energy depends only on the subset of sites, $\{i_1,\mathrm{\cdots },i_k\}$, without respect to the order of sites. It is a consequence of detailed balance at thermodynamic equilibrium that the order in which the perturbations are undertaken does not matter. For example, it is clear from Equation 12 that ${\mathrm{\Delta }}^2{\gamma }_{\mathrm{\varnothing },(i_1,i_2)}={\mathrm{\Delta }}^2{\gamma }_{\mathrm{\varnothing },(i_2,i_1)}$. More generally, if ρ is any permutation of the perturbed sites, so that ρ is a bijective function, $\rho :\{i_1,\mathrm{\cdots },i_k\}\to \{i_1,\mathrm{\cdots },i_k\}$, then it can be shown that

Note that Equation 16 follows from Equation 15 when $Z=\mathrm{\varnothing }$. This property of invariance under permutation is referred to as ‘symmetry’ in Horovitz and Fersht, 1990 and is similar to the algebraic relations which give rise to the independent HOCs, ${\omega }_{i,S}$ with $i<S$, as described previously.

The equality between the higher-order couplings introduced in Horovitz and Fersht, 1990 and the residual free energies introduced in Horovitz and Fersht, 1992, as described in Equation 15, is presumably well known to those in the field. It seems to be implicitly assumed in Horovitz and Fersht, 1992, but we have not found a clear statement of it in the literature. It would be difficult to formulate one in the absence of a general definition of higher-order coupling, as we have given in Equation 11. The formulas above may therefore be of some value in offering a rigorous treatment.

Each of the measures we have discussed, HOCs, residual free energies and higher-order couplings, offers a different way of analysing the free-energy landscape using the graphs in Figure 5. HOCs are associated to graph edges; residual free energies are associated to graph vertices; and higher-order couplings are associated to sequences of sites, at least when symmetries are ignored. As we have seen above, the three measures are mathematically equivalent. However, they are useful for different purposes. HOCs tell us about the integration of binding information; residual free energies capture the collective synergy between sets of sites; and higher-order couplings show how these same synergies can be extracted from a sequence of experimental perturbations. One advantage of HOCs is that they are non-dimensional quantities in terms of which it is straightforward to calculate the other measures. By doing so, we were able to show rigorously that higher-order couplings are also residual free energies (Equation 15).

Having explained how various higher-order measures are related to each other, we return to the question of how effective cooperativity arises from allosteric ensembles with multiple conformations. For this problem, HOCs are much easier to use than either residual free energies or higher-order couplings. With Equations 8 and 14 now available, effective residual free energies or effective higher-order couplings may be calculated from the effective HOCs that we construct below, but we will not exploit this capability in the present paper.

As MWC showed, even if there is no intrinsic cooperativity in any conformation, an effective cooperativity can arise from the ensemble. This is usually detected in the shape of the binding function (Figure 2A). Here, we introduce a method of coarse graining through which effective cooperativities can be rigorously defined. We illustrate this for the allostery graph, $A$, and explain the general coarse-graining method in the Materials and methods. For allostery, the idea is to treat the horizontal subgraphs, $A_S$, as the vertices of a new coarse-grained graph, $A^{\varphi }$, (Figure 4, bottom right). There is an edge between two vertices in $A^{\varphi }$, if, and only if, there is an edge in $A$ between the corresponding horizontal subgraphs. It is not hard to see that $A^{\varphi }$ is identical in structure to any of the vertical subgraphs $A^{c_k}$. We can think of $A^{\varphi }$ as if it represents a single effective conformation to which ligand is binding, and we can index each vertex of $A^{\varphi }$ by the corresponding subset of bound sites, $S$. The key point, as explained in detail in the Materials and methods, is that it is possible to assign labels to the edges in $A^{\varphi }$ so that

with $A^{\varphi }$ being at thermodynamic equilibrium under these label assignments. According to Equation 17, the probability of being in a coarse-grained vertex of $A^{\varphi }$ is identical to the overall probability of being in any of the corresponding vertices of $A$. This is exactly the property a coarse graining should satisfy at steady state. It is not difficult to see why a procedure like this should work. The coarse-graining formula in Equation 17 tells us the expected probability distribution on the coarse-grained graph, $A^{\varphi }$. Equation 3 can then be used to back out the equilibrium labels on the edges of $A^{\varphi }$ which give rise to this probability distribution. We provide a more direct way of achieving the same result in Equation 40. This assignment of labels to $A^{\varphi }$ is the only way to ensure Equation 17 at equilibrium, so that the coarse graining is both systematic and unique. The Materials and methods gives a more careful treatment for coarse graining any linear framework graph, which may not itself be at thermodynamic equilibrium.

Our coarse-graining procedure offers a general method for calculating how effective behaviour emerges, at thermodynamic equilibrium, from a more detailed underlying mechanism. This procedure is likely to be broadly useful for other studies. We note that it applies only to the steady state. It does not provide a coarse graining of the underlying dynamics, which is a much harder problem.

Because $A^{\varphi }$ resembles the graph for ligand binding at a single conformation, we can calculate HOCs for $A^{\varphi }$—equivalently, effective HOCs for $A$—just as we did above, by normalising the effective association constants. Once the dust of calculation has settled (Materials and methods), we find that $A$ has effective association constants and effective HOCs:

The quantity ${\mu }_S(A^{c_k})$ is calculated by multiplying labels over paths, as above, within the vertical subgraph $A^{c_k}$. The terms within angle brackets, of the form $⟨X(c_k)⟩$, where $X(c_k)$ is some function over conformations c_k, denote averages over the steady-state probability distribution of the horizontal subgraph: $⟨X(c_k)⟩={\sum }_{1\le k\le N}X(c_k){\text{Pr}}_{c_k}(A_{\mathrm{\varnothing }})$. The right-hand formula in Equation 18 for the effective HOCs has a suggestive structure: it is an average of a product divided by the product of the averages. The effective parameters in Equation 18 provide a biophysical language in which the integrative capabilities of any ensemble can be rigorously specified.

The functional viewpoint is readily recovered from the ensemble. A generalised MWC formula can be given in terms of effective HOCs, from which the classical two-conformation MWC formula is easily derived (Materials and methods). Some expected properties of effective HOCs are also easily checked (Materials and methods). First, ${\omega }_{i,S}^{\varphi }$ is independent of ligand concentration, $x$. Second, there is no effective HOC for binding to an empty conformation, so that ${\omega }_{i,\mathrm{\varnothing }}^{\varphi }=1$. Third, if there is only one conformation c₁, then the effective HOC reduces to the intrinsic HOC, so that ${\omega }_{i,S}^{\varphi }={\omega }_{c_1,i,S}$.

More illuminating are the effective HOCs for the MWC model. We consider any conformational ensemble which is MWC-like: there is no intrinsic HOC in any conformation, so that ${\omega }_{c_k,i,S}=1$ and $K_{c_k,i,S}=K_{c_k,i,\mathrm{\varnothing }}$; and the bare association constants are identical at all sites, so that we can set $K_{c_k,i,\mathrm{\varnothing }}=K_{c_k}$. There may, however, be any number of conformations, not just the two conformations of the classical MWC model. It then follows that ${\omega }_{i,S}^{\varphi }$ depends only on the size of $S$, so that we can write ${\omega }_{i,S}^{\varphi }$ as ${\omega }_s^{\varphi }$, where $s=\mathrm{\#}(S)$ is the order of cooperativity. Equation 18 then simplifies to (Materials and methods)

We see that, although there is no intrinsic HOC in any conformation, effective HOC of each order arises from the moments of $K_{c_k}$ over the probability distribution on $A_{\mathrm{\varnothing }}$. In particular, Equation 19 shows that the effective pairwise cooperativity is ${\omega }_1^{\varphi }=⟨{(K_{c_k})}^2⟩/{⟨K_{c_k}⟩}^2$.

In studies of G-protein coupled receptor (GPCR) allostery, Ehlert relates ‘empirical’ to ‘ultimate’ levels of explanation by a procedure similar to our coarse graining, but with only two conformations, and calculates a ‘cooperativity constant’ which is the same as ${\omega }_1^{\varphi }$ (Ehlert, 2016). Gruber and Horovitz calculate ‘successive ligand binding constants’ for the two-conformation MWC model which are the same as effective association constants, $K_s^{\varphi }$, (Gruber and Horovitz, 2018) (Materials and methods). To our knowledge, these are the only other calculations of effective allosteric quantities. We note that Equation 19 applies to all HOCs, not just pairwise, and to any MWC-like ensemble, not just those with two conformations.

The classical MWC model yields only positive cooperativity (Koshland and Hamadani, 2002; Monod et al., 1965), as measured in the functional perspective (Figure 2A). We find that MWC-like ensembles yield positive effective HOCs of all orders. Strikingly, these effective HOCs increase with increasing order of cooperativity: provided $K_{c_k}$ is not constant over conformations (Materials and methods),

This shows that ensembles with independent and identical sites, including the two-conformation MWC model, can effectively implement high orders and high levels of positive cooperativity. Equation 20 is very informative, and we return to it in the Discussion.

It is often suggested that negative cooperativity requires a different kind of ensemble to those considered here, such as one allowing KNF-style induced fit (Koshland and Hamadani, 2002). However, if two sites are independent but not identical, so that $K_{c_k,1,\mathrm{\varnothing }}\ne K_{c_k,2,\mathrm{\varnothing }}$, then, with just two conformations, the effective pairwise cooperativity can become negative. Indeed, ${\omega }_{1,\{2\}}^{\varphi }<1$, if, and only if, the values of the association constants are not in the same relative order in the two conformations (Materials and methods). Negative effective cooperativity can arise from non-identical sites and does not need a special kind of ensemble.

Equation 18 shows that effective HOCs of any order can arise for a conformational ensemble but does not reveal what values they can attain. Can they vary arbitrarily? The question can be rigorously posed as follows. Suppose that we are considering $n$ binding sites and that numbers ${\beta }_i>0$, for $1\le i\le n$, and ${\alpha }_{i,S}>0$, for $i<S$, are chosen at will. Does there exist a conformational ensemble such that the bare effective association constants satisfy $K_{i,\mathrm{\varnothing }}^{\varphi }={\beta }_i$, and the independent effective HOCs satisfy ${\omega }_{i,S}^{\varphi }={\alpha }_{i,S}$?

To address this question, we assume that there is no intrinsic HOC, so as not to introduce cryptically what we want to generate. It follows that the sites cannot be identical, for otherwise Equation 20 shows that integrative flexibility is impossible. Accordingly, the bare association constants, $K_{c_k,i,\mathrm{\varnothing }}$ for $1\le i\le n$, can be treated as $n$ free parameters in each conformation c_k. If there are $N$ conformations in the ensemble, then there are $N-1$ free parameters coming from the horizontal edges (Materials and methods). Dimensional considerations imply that the effective HOCs cannot take arbitrary values if $n(N-1)<2^n-1$. Conversely, we prove the following flexibility theorem: any pattern of values can be realised by an allosteric ensemble with no intrinsic cooperativity, to any required degree of accuracy, provided there are enough conformations with the right probability distribution and the right patterns of bare association constants.

To see why this is possible, we outline the argument here and give rigorous details in Theorem 1 in the Materials and methods. Other arguments may of course be possible and the details presented here should not be thought of as the only way for the results to hold. We will use an allostery graph $A$ whose conformations are indexed by subsets $T\subseteq \{1,\mathrm{\cdots },n\}$ and denoted $c_T$. Both binding subsets and conformations will then be indexed by subsets of $\{1,\mathrm{\cdots },n\}$. To avoid confusion, we will use $S$ to label binding subsets and $T$ to label conformations, so that a vertex of $A$ will be $(c_T,S)$. The allostery graph for the case $n=2$ is shown in Figure 6. We will focus on the horizontal subgraph of empty conformations, $A_{\mathrm{\varnothing }}$, because that is what is needed for calculating effective HOCs using Equation 18. We will take the reference vertex of $A_{\mathrm{\varnothing }}$ to be $c_{\mathrm{\varnothing }}$. Recall from what was explained previously that the product of the equilibrium labels along any path in $A_{\mathrm{\varnothing }}$ from the reference vertex to the vertex $c_T$ is the quantity ${\mu }_{c_T}(A_{\mathrm{\varnothing }})$, in terms of which the steady-state probabilities of $A_{\mathrm{\varnothing }}$ are given by Equation 3. Let ${\lambda }_T={\mu }_{c_T}(A_{\mathrm{\varnothing }})$. These quantities are $2^n-1$ free parameters whose values we are going to assign. They are more convenient for our purposes than an independent set of equilibrium labels for $A_{\mathrm{\varnothing }}$. By Equation 3,

Figure 6

Download asset Open asset

The other free parameters that we need are $n$ quantities, ${\kappa }_1,\mathrm{\cdots },{\kappa }_n>0$, to which we will subsequently assign values, in terms of which we will define the intrinsic association constants. We will assume that the sites are independent in each conformation, so that all intrinsic HOCs of $A$ are 1. It follows that $K_{c_T,i,S}=K_{c_T,i,\mathrm{\varnothing }}$. We then set $K_{c_T,i,\mathrm{\varnothing }}={\kappa }_i$ if $i\in T$, and $K_{c_T,i,\mathrm{\varnothing }}=\epsilon {\kappa }_i$ if $i\notin T$. Here, ε is a small positive quantity which can be chosen to determine the degree of accuracy to which the ${\beta }_i$ and ${\alpha }_{i,S}$ are approximated. In the calculations which follow, we will only be interested in terms which do not involve ε as a factor. Because the sites are independent in each conformation, it follows that, in the vertical subgraph, $A^{c_T}$, at any conformation $c_T$, ${\mu }_S(A^{c_T})=({\prod }_{i\in S}{\kappa }_i)x^{\mathrm{\#}(S)}$, whenever $S\subseteq T$. However, if $S⊈T$, then ${\mu }_S(A^{c_T})$ acquires factors of ε and so ${\mu }_S(A^{c_T})\approx 0$, where ≈ means simply that the related quantities become equal as ε becomes very small. In this case, for our purposes, ${\mu }_S(A^{c_T})$ is negligible whenever $S⊈T$. Figure 6 illustrates how this plays out in the allostery graph for $n=2$.

To calculate the effective association constants, the left-hand formula in Equation 18 shows that we must evaluate the averages $⟨K_{c_T,i,S}.{\mu }_S(A^{c_T})⟩$ and $⟨{\mu }_S(A^{c_T})⟩$. Using Equation 21,

The only terms in the sum which do not involve ε as a factor are those $T$ for which $S\subseteq T$. Furthermore, the definition of ${\mu }_S(A^{c_T})$ given above shows that these terms do not depend on $T$. Similarly, using Equation 21 again,

and the only terms in the sum which do not involve ε as a factor are those for which $S\cup \{i\}\subseteq T$. These terms also do not depend on $T$. It follows from Equation 18 that

where we have ignored all terms involving ε as a factor.

Equation 22 tells us two things. First, that the effective association constants are approximately proportional to the corresponding κ’s. Hence, if the proportionality constants, which depend only on the ${\lambda }_T$, are determined, we can choose the ${\kappa }_i$ so as to make the bare effective association constants $K_{i,\mathrm{\varnothing }}^{\varphi }$ approximately equal to ${\beta }_i$. Second, Equation 22 tells us that the effective HOCs, ${\omega }_{i,S}^{\varphi }$, are independent of the ${\kappa }_i$ and depend only on the ${\lambda }_T$,

It remains for us to assign values to the ${\lambda }_T$ so that the effective HOCs become approximately equal to the α’s.

To do this, assume that, for the conformation $c_T$, the subset $T$ is written as $T=\{i_1,\mathrm{\cdots }{i}_k\}$, where the indices are in increasing order, ${\displaystyle i_1<i_2<\cdots <i_k}$. Because of this ordering, the quantities ${\alpha }_{i_j,\{i_{j+1},\mathrm{\cdots },i_k\}}$ are given to us by hypothesis. Hence, we can define

Here, δ is another small positive quantity, similar to ε, which can be chosen to set the degree of accuracy to which the β’s and α’s are approximated. As with ε, we will treat as negligible terms in which δ is a factor. Figure 6 illustrates Equation 24 for the case $n=2$.

It can be seen from Equation 24 that ${\sum }_{X\subseteq T}{\lambda }_T={\lambda }_X(1+U)$, where $U$ has a factor of δ and is therefore negligible as δ becomes very small, $U\approx 0$. It then follows from Equation 23 that

where we have used $U$ as a generic symbol for quantities which are negligible as δ becomes very small. By Equation 24, ${\lambda }_{S\cup \{i\}}={\alpha }_{i,S}\delta {\lambda }_S$, so that

This establishes part of what is required. For the other part, we can return to Equation 22 and set

from which it follows from Equation 22 that

Equations 26 and 27 show that the effective association constants and effective HOCs can take arbitrary positive values to any desired degree of accuracy, as determined by ε and δ. This establishes the flexibility theorem. The Materials and methods provides a more careful treatment in Theorem 1, which rigorously establishes the approximation as ε and δ become very small.

Figures 7 and 8 together illustrate the flexibility theorem. Figure 7A shows three arbitrarily chosen patterns of effective parameters for a target molecule with four ligand binding sites. Figure 7B shows the corresponding overall binding functions (black curves) together with the individual site-specific binding functions (coloured curves). As a matter of thermodynamics, the overall binding function is always an increasing function of ligand concentration. In contrast, the site-specific binding functions may increase or decrease depending on the combinations of positive and negative effective HOCs in Figure 7A, and thereby show more clearly the complexity arising from those different combinations. The implementation of the effective parameters by an allosteric ensemble, as specified by the flexibility theorem, is illustrated in Figure 8. Figure 8A shows the allosteric ensemble for $n=4$ sites as a product graph with 16 binding patterns and 16 conformations. Figure 8B shows the intrinsic association constants in each conformation coming from the proof of the flexibility theorem, to an accuracy of 0.01. Figure 8C confirms that this allosteric ensemble exactly reproduces the overall binding functions in Figure 7B.

Figure 7

Download asset Open asset

Figure 8

Download asset Open asset

In respect of the dimensional argument made previously, the allostery graph used in the proof above has $2^n-1$ free parameters for $A_{\mathrm{\varnothing }}$ and the ${\kappa }_1,\mathrm{\cdots },{\kappa }_n$ are a further $n$ free parameters, giving $2^n-1+n$ free parameters in total. This is more than the minimal required number of $2^n-1$ but not by much. It remains an interesting open question whether a conformational ensemble can be constructed, perhaps with more free parameters, which gives the effective HOCs exactly, rather than only approximately. One consequence of the definitions of $K_{c_T,i,\mathrm{\varnothing }}$ and of ${\lambda }_T$ in Equation 24 is that the parameters of the allosteric ensemble become exponentially small, as is evident for the examples in Figure 8B. Another interesting question is whether alternative constructions can be found which do not exhibit such a broad range of parameter values. Irrespective of these questions, the proof given above confirms that there is no fundamental biophysical limitation to achieving any pattern of values to any desired degree of accuracy. Accordingly, a central result of the present paper is that sufficiently complex allosteric ensembles can implement any form of information integration that is achievable without energy expenditure.

As mentioned in the Introduction, the starting point for the present paper was to account for experimental data on gene expression. Studies in Drosophila have shown that the Hunchback gene, in response to the maternal TF Bicoid, is sharply expressed in a way that is well fitted, after appropriate normalisation, to a Hill function, ${\mathscr{H}}_h(x)=x^h/(1+x^h)$. This sharp expression underlies the initial patterning of anterior-posterior stripes in the early Drosophila embryo. Estimated values for the Hill coefficient, $h$, vary depending on the experimental construct and time of measurement but are typically in the range $4\le h\le 8$ during early nuclear cycle 14 (Tran et al., 2018). The relevant promoter is believed to have $n=6$ Bicoid binding sites, and the mechanistic basis for the sharpness is the subject of considerable interest. We showed in previous work that, if the promoter was assumed to have six Bicoid binding sites and to be operating at thermodynamic equilibrium, then the highest Hill coefficient that could be achieved of $h=6$, at the so-called Hopfield barrier, required HOCs for Bicoid binding of order up to 5 (Estrada et al., 2016). In particular, pairwise cooperativities, which had previously been invoked to account for the sharpness (Gregor et al., 2007), are not sufficient to explain the data. Left open by this previous work was a molecular mechanism which could create the high-order HOCs required for Hill functions. We have seen above that allosteric ensembles can create any pattern of HOCs, so it is natural to ask if there are allosteric ensembles which yield good approximations to Hill functions.

We implemented a numerical optimisation algorithm to find binding functions which approximated Hill functions (Materials and methods). Hill functions are naturally normalised so that ${\mathscr{H}}_h(1)=0.5$, so we followed the procedure introduced previously (Estrada et al., 2016) of normalising concentration to its value at half-maximum: if the normalised binding function is denoted $f(x)$, then $f(1)=0.5$. Figure 9 shows results for an allosteric ensemble with four conformations for ligand binding to six sites. The ensemble has no intrinsic cooperativity in any conformation, so that $K_{c_k,i,S}=K_{c_k,i,\mathrm{\varnothing }}$ for any binding subset $S\subseteq \{1,\mathrm{\cdots },6\}$, while the bare association constants, $K_{c_k,i,\mathrm{\varnothing }}$, differ between the conformations (Figure 9B). This gives $4\times 6=24$ free parameters together with an additional three free parameters for the independent equilibrium labels on the horizontal subgraph $A_{\mathrm{\varnothing }}$ (Figure 9A). We limited the parameter ranges so that the $K_{c_k,i,\mathrm{\varnothing }}$ were in the range $[{10}^{-4},{10}^4]$ and the equilibrium labels of $A_{\mathrm{\varnothing }}$ were in the range $[{10}^{-6},{10}^6]$. With these settings, it was not difficult to find normalised binding functions which are very well fitted by the Hill function, ${\mathscr{H}}_h(x)$, for Hill coefficients $h=4$, 5 and 6 (Figure 9D).

Figure 9

Download asset Open asset

We were able to find multiple sets of parameters which yielded excellent fits; Figure 9 shows two representative examples for each Hill coefficient. It is evident that very different numerical ensembles (Figure 9B) can give almost identical binding functions (Figure 9D). This reinforces the point made in the Introduction that the binding function, or some associated measure of its shape, such as a Hill coefficient, are aggregate measures which give little insight into how binding information is being integrated. For this, the patterns of effective parameters provide more detailed information. As can be seen from Figure 9C, effective HOCs of all orders up to 5 are needed for each Hill function, as suggested previously (Estrada et al., 2016), with predominantly positive effective HOCs, ${\omega }_{i,S}^{\varphi }>1$, and varying amounts of independence, ${\omega }_{i,S}^{\varphi }=1$.

It is interesting to ask what role the size of the ensemble plays in approximating Hill functions. We cannot give a definitive answer but can make some observations. We were able to approximate ${\mathscr{H}}_6$ with a two-conformation ensemble with six sites but only with much wider parametric ranges. It was also more difficult in terms of optimisation time to find a good fit, and we did not find multiple fits. This suggests that the larger the ensemble the easier it is to approximate Hill functions with limited parameter ranges. It is also conceivable that the size of the ensemble may have to increase with the number of binding sites to retain control over the parametric ranges. We must leave such issues to subsequent work. While our results are numerical, and therefore limited to the ensemble we have analysed, it seems clear that allosteric ensembles provide a molecular mechanism that can closely approximate Hill functions with the required high orders of effective cooperativity, thereby providing a solution to our original question. Since Hill functions are widely used to fit data, the potential for an underlying allosteric mechanism may be broadly useful.

Graphs will always be connected, so that they cannot be separated into sub-graphs between which there are no edges. The set of vertices of a graph $G$ will be denoted by $\nu (G)$. For a general graph, the vertices will be indexed by numbers $1,\mathrm{\cdots },N\in \nu (G)$ and vertex 1 will be taken to be the reference vertex. Particular kinds of graphs, such as the allostery graphs discussed in the paper, may use a different indexing. An edge from vertex $i$ to vertex $j$ will be denoted $i\to j$ and the label on that edge by $\mathrm{\ell }(i\to j)$. A subscript, as in $i{\to }_{G}j$, may be used to specify which graph is under discussion. When discussing graphs, we used the word ‘structure’ to refer to properties that depend on vertices and edges only, ignoring the labels.

A graph gives rise to a dynamical system by assuming that each edge is a chemical reaction under mass-action kinetics with the label as the rate constant. Since each edge has only a single source vertex, the corresponding dynamics is linear and can be represented by a linear differential equation in matrix form:

Here, $G$ is the graph, $u$ is a vector of component concentrations and $ℒ(G)$ is the Laplacian matrix of $G$. Since material is only moved between vertices, there is a conservation law, ${\sum }_i{u}_i(t)=u_{tot}$. By setting $u_{tot}=1$, $u$ can be treated as a vector of probabilities. In such a stochastic setting, Equation 28 is the master equation (Kolmogorov forward equation) of the underlying Markov process. This is a general representation: given any well-behaved Markov process on a finite state space, there is a graph, whose vertices are the states, for which Equation 28 is the master equation.

The linear dynamics in Equation 28 gives the linear framework its name and is common to all applications. The treatment of the external components, which appear in the edge labels and which introduce nonlinearities, depends on the application. For the case of allostery treated here, we make the same assumptions as in thermodynamics for the grand canonical ensemble, with each ligand being present in a reservoir from which binding and unbinding to graph vertices does not change its free concentration. In this case, the edge labels are effectively constant. The same assumptions are implicitly used in other studies of allostery.

The dynamics in Equation 28 always tends to a steady state, at which $du/dt=0$, and, under the fundamental timescale separation, it is assumed to have reached a steady state. If the graph is strongly connected, it has a unique steady state up to a scalar multiple, so that $dim\ker ℒ(G)=1$. Strong connectivity means that, given any two distinct vertices, $i$ and $j$, there is a path of directed edges from $i$ to $j$, $i=i_1\to i_2\to \mathrm{\cdots }\to i_{k-1}\to i_k=j$. Under strong connectivity, a representative steady state for the dynamics, $\rho (G)\in \ker ℒ(G)$, may be calculated in terms of the edge labels by the Matrix Tree Theorem. We omit the corresponding expression as it is not needed here, but it can be found in any of the references given above. This expression holds whether or not the steady state is one of thermodynamic equilibrium. However, at thermodynamic equilibrium, the description of the steady state simplifies considerably because detailed balance holds. This means that the graph is reversible, so that, if $i\to j$, then also $j\to i$, and each pair of such edges is independently in flux balance, so that

This ‘microscopic reversibility’ is a fundamental property of thermodynamic equilibrium. Note that a reversible, connected graph is necessarily strongly connected.

Take any path of reversible edges from the reference vertex 1 to some vertex $i$, $1=i_1\rightleftharpoons i_2\rightleftharpoons \mathrm{\cdots }\rightleftharpoons i_{k-1}\to i_k=i$, and let ${\mu }_i(G)$ be the product of the label ratios along the path:

It is straightforward to see from Equation 29 that ${\mu }_i(G)$ does not depend on the chosen path and that ${\rho }_i(G)={\mu }_i(G){\rho }_1(G)$. The vector $\mu (G)$ is therefore a scalar multiple of $\rho (G)$ and so also a steady state for the dynamics. The detailed balance formula in Equation 29 also holds for μ in place of ρ. At thermodynamic equilibrium, the only parameters needed to describe steady states are label ratios.

This observation about label ratios leads to the concept of an equilibrium graph. Suppose that $G$ is a linear framework graph which can reach thermodynamic equilibrium and is therefore reversible (above). $G$ gives rise to an equilibrium graph, $ℰ(G)$, as follows. The vertices and edges of $ℰ(G)$ are the same as those of $G$, but the edge labels in $ℰ(G)$, which we will refer to as ‘equilibrium edge labels’ and denote ${\displaystyle {\ell }_{\mathrm{e}\mathrm{q}}(i\to j)}$, are the label ratios in $G$. In other words,

Scheme 1 illustrates the relationship between the linear framework graph and the corresponding equilibrium graph. Note that the equilibrium edge labels of $ℰ(G)$ are non-dimensional and that ${\displaystyle {\ell }_{\mathrm{e}\mathrm{q}}(j\to i)={\ell }_{\mathrm{e}\mathrm{q}}(i\to j{)}^{-1}}$. The equilibrium edge labels are the essential parameters for describing a state of thermodynamic equilibrium.

Scheme 1

Download asset Open asset

These parameters are not independent because Equation 29 implies algebraic relationships among them. Indeed, Equation 29 is equivalent to the following ‘cycle condition’, which we formulate for $ℰ(G)$: given any cycle of edges, $i_1\to i_2\to \mathrm{\cdots }\to i_{k-1}\to i_1$, the product of the equilibrium edge labels along the cycle is always 1:

This cycle condition is equivalent to the detailed balance condition in Equation 29 and either condition is equivalent to $G$ being at thermodynamic equilibrium.

There is a systematic procedure for choosing a set of equilibrium edge label parameters which are both independent, so that there are no algebraic relationships among them, and also complete, so that all other equilibrium edge labels can be algebraically calculated from them. Recall that a spanning tree of $G$ is a connected subgraph, $T$, which contains each vertex of $G$ (spanning) and which has no cycles when edge directions are ignored (tree). Any strongly connected graph has a spanning tree and the number of edges in such a tree is one less than the number of vertices in the graph. Since $G$ and $ℰ(G)$ have the same vertices and edges, they have identical spanning trees. The equilibrium edge labels ${\displaystyle {\ell }_{\mathrm{e}\mathrm{q}}(i{\to }_{T}j)}$, taken over all edges $i\to j$ of $T$, form a complete and independent set of parameters at thermodynamic equilibrium. In particular, if $G$ has $N$ vertices, there are $N-1$ independent parameters at thermodynamic equilibrium.

In the main text, we defined an equilibrium allostery graph, $A$ (Figure 4), without specifying a corresponding linear framework graph, $G$, for which $ℰ(G)=A$. Because label ratios are used in an equilibrium graph, there is no unique linear framework graph corresponding to it. However, some choice of transition rates, $\mathrm{\ell }(i{\to }_{G}j)$ and $\mathrm{\ell }(j{\to }_{G}i)$, can always be made such that their ratio is ${\displaystyle {\ell }_{\mathrm{e}\mathrm{q}}(i{\to }_{\mathcal{ℰ}\mathcal{(}\mathcal{𝒢}\mathcal{)}}j)}$. Hence, some linear framework graph $G$ can always be defined such that $ℰ(G)=A$. In some of the constructions below, we will work with the linear framework graph, $G$, rather than with the equilibrium graph $A$ and will then show that the construction does not depend on the choice of $G$.

The steady-state probability of vertex $i$, ${\text{Pr}}_i(G)$, can be calculated from the steady state of the dynamics by normalising, so that

where the first formula holds for any strongly connected graph and the second formula also holds if the graph is at thermodynamic equilibrium. In the latter case, Equation 29 holds and $\mu (G)$ can be defined by Equation 30. The second formula in Equation 33 corresponds to Equation 3. If the graph is at thermodynamic equilibrium, the equilibrium edge labels may be interpreted thermodynamically, as illustrated in Figure 3 and discussed in the main text (Equation 1):

If Equation 34 is used to expand the second formula in Equation 33, it gives the specification of equilibrium statistical mechanics for the grand canonical ensemble, with the denominator being the partition function.

It will be helpful to let $\mathrm{\Pi }(G)$ and $\mathrm{\Psi }(G)$ denote the corresponding denominators in Equation 33, so that $\mathrm{\Pi }(G)={\rho }_1(G)+\mathrm{\cdots }+{\rho }_N(G)$ for any strongly connected graph and $\mathrm{\Psi }(G)={\mu }_1(G)+\mathrm{\cdots }+{\mu }_N(G)$ for a graph which is at thermodynamic equilibrium. We will refer to $\mathrm{\Pi }(G)$ and $\mathrm{\Psi }(G)$ as partition functions. It follows from Equation 33 that

depending on the context.

We can choose any spanning tree in the horizontal subgraph of empty conformations, $A_{\mathrm{\varnothing }}$. As explained above, the equilibrium labels on the edges of this tree define a complete set of $N-1$ independent parameters for $A_{\mathrm{\varnothing }}$. As for the vertical subgraphs, $A^{c_k}$, which all have the same structure, consider the subgraph of $A^{c_k}$ consisting of all edges, together with the corresponding source and target vertices, of the form, $(c_k,S)\to (c_k,S\cup \{i\})$, where $\mathrm{\varnothing }\subseteq S\subset \{1,\mathrm{\cdots },n\}$ and $i$ is less than all the sites in $S$ ($i<S$). It is not difficult to see that this subgraph is a spanning tree of $A^{c_k}$ (Estrada et al., 2016, SI, §3.2). Accordingly, the association constants, $K_{c_k,i,S}$ from Equation 36, with $i<S$, form a complete set of independent parameters for $A^{c_k}$. Because of the product structure of $A$, adjoining the spanning trees in $A^{c_k}$, for each conformation c_k with $1\le k\le N$, to the spanning tree in $A_{\mathrm{\varnothing }}$, yields a spanning tree in $A$. Hence, the independent parameters for $A^{c_k}$ together with the $N-1$ independent parameters for $A_{\mathrm{\varnothing }}$ are also collectively independent as parameters for $A$. It follows from the description of labels above that these parameters are also complete for $A$, so that any equilibrium label in $A$ can be expressed in terms of them.