The Allostery Landscape: Quantifying Thermodynamic Couplings in Biomolecular Systems.

Allostery plays a fundamental role in most biological processes. However, little theory is available to describe it outside of two-state models. Here we use a statistical mechanical approach to show that the allosteric coupling between two collective variables is not a single number, but instead a two-dimensional thermodynamic coupling function that is directly related to the mutual information from information theory and the copula density function from probability theory. On this basis, we demonstrate how to quantify the contribution of specific energy terms to this thermodynamic coupling function, enabling an approximate decomposition that reveals the mechanism of allostery. We illustrate the thermodynamic coupling function and its use by showing how allosteric coupling in the alanine dipeptide molecule contributes to the overall shape of the Φ/Ψ free energy surface, and by identifying the interactions that are necessary for this coupling.

Introduction

Allostery plays a fundamental role in most biological processes and has been suggested to be present in nearly all proteins.[1] One of the best-studied allosteric phenomena is the activation of a receptor, which we will denote as R, by a ligand, denoted as L. The most common model for allostery in this system is the allosteric two-state model (ATSM).[2,3] We can construct a thermodynamic cycle for the process of ligand-induced activation of the receptor:We refer to this as an “allosteric cycle”. To describe the allostery in this system, the allosteric efficacy, α, can be calculated from the cycle aswhere (assuming the volume is constant), each equilibrium constant is a function of the difference in Helmholtz free energy, A, for the two states:Thus,For convenience, we will discuss the allosteric efficacy in terms of the quantity ΔΔA, which we will call the thermodynamic coupling.

ΔΔA is symmetric at equilibrium, that is, ΔΔA for receptor activation conditional on ligand binding is equivalent to the ΔΔA for ligand binding conditional on activation. Thus, the following two definitions of ΔΔA are equivalent:However, there is no reason to assume that the receptor activation is a two-state process. In fact, NMR experiments have revealed a multimodal activation process in the β2-adrenergic receptor (β2AR),[4] and quantitative mass spectroscopy experiments have revealed ligand-specific states in the same system.[5] These results, along with other evidence for additional states in β2AR and other receptors,[6−11] indicate that activation must be treated as either involving more than two discrete states, or even as involving a continuous conformational space.

Receptor activation involves not only multiple states but also multiple dimensions. The complex behavior of an allosteric receptor is thus unlikely to be well described by a single reaction coordinate; instead, the large number of potential conformational states may be best described by multiple collective variables (CVs; variables that are functions of the atomic coordinates) that are thermodynamically coupled in nontrivial ways. It is impossible to understand the molecular mechanism of ligand-induced receptor activation without explicitly considering the thermodynamic coupling between the ligand binding site and active site, so a minimal set of CVs should at least include one CV for each one of these sites. These CVs are intrinsic to the specific receptor, and their thermodynamic coupling arises from the complex network of molecular interactions that separates them spatially in the receptor’s structure. This intrinsic thermodynamic coupling is of great interest, as an understanding of the nature of this coupling can be used to inform both the design of ligands that modulate function in a highly specific manner, and the design of receptors with modified allosteric properties. We have previously represented this intrinsic thermodynamic coupling using the recently developed Allosteric Ising Model (AIM),[12] a two-state model of allostery that implicitly includes the potential energy of interaction between structural components. While the AIM and other statistical mechanical models of allostery, such as the ensemble allosteric model,[13] have allowed us to derive some analytical features of simple allosteric systems, a general method that does not rely on the two-state assumption is still needed to study the intrinsic thermodynamic coupling between structural components in real systems. Here we describe such a method and illustrate its capabilities by showing how allosteric coupling in the alanine dipeptide molecule contributes to the overall shape of the Φ/Ψ free energy surface, and by identifying the interactions that are necessary for this coupling and their contributions to the energetics.

Theoretical Developments

Derivation

To quantify the intrinsic thermodynamic coupling between CVs, we will derive expressions analogous to the allosteric efficacy for the coupled perturbation of discrete or continuous CVs away from their equilibrium distributions. Let represent the coordinates of the allosteric protein and its environment that define our system, which does not include any ligand that we consider here as external perturbations. The probability density of each microstate r⃗ is given by the Boltzmann distribution,where U(r⃗) is the potential energy function. The numerator is the Boltzmann factor denoted asThe free energy can be written as a functional of the Boltzmann factor function,We define a CV as a function X(r⃗) of the system’s coordinates that can be either continuous or discrete. For a continuous CV, the probability density function isFor a discrete CV, the probability mass function, p(x), is defined by an identical expression, but is bounded to be less than 1 everywhere. We can calculate the free energy conditional on a value of the CV asEquation can be rewritten in terms of either f(x) or p(x). Because we use the histogram method to estimate the probability mass function of the CVs in the application following this derivation, we will assume discrete CVs described using p(x) without loss of generality. The free energy becomesConsider a second CV, Y(r⃗), with analogous probability function and free energy definitions. A joint probability mass function for the two CVs can be written asso that the analogous free energy conditional on values of both CVs is

One can imagine X(r⃗) to describe the ligand binding site and Y(r⃗) to describe the active site of the protein; the binding of a ligand to the system then acts as an external perturbation to the distributions of these CVs. To quantify the intrinsic coupling between these CVs, we apply artificial perturbations to the equilibrium CV distributions such that one or both CVs become constrained to a given value. From the equilibrium state and these artificially perturbed states, we calculate the allosteric efficacy of the following thermodynamic cycle:We will refer to this class of thermodynamic cycles as “thermodynamic perturbation cycles”. The thermodynamic coupling of the perturbations at position (x,y) in the CV space, ΔΔA(x,y), can be calculated asEquation simplifies toThis is the mathematical definition we propose for the central quantity ΔΔA(x,y) that we call the thermodynamic coupling function for the CVs X(r⃗) and Y(r⃗). In two dimensions, eq defines what we call the allostery landscape (see Figure for an example).

Figure 1

Allostery in the alanine dipeptide. Our calculations reveal that the allosteric coupling between Φ and Ψ destabilizes the high free energy regions and stabilizes the αL and C′7ax states. (a) The alanine dipeptide molecule with the backbone dihedral angles Φ and Ψ indicated by arrows. The molecule is partitioned in three domains as indicated by the black lines (see text for details). (b) Free energy surface A(Φ,Ψ) calculated according to eq , with prominent states labeled in white. In all panels, Φ and Ψ are expressed in radians. (c) The allostery landscape representing the thermodynamic coupling between CVs Φ and Ψ, calculated according to eq . (d) The normalized allosteric coupling, calculated according to eq . In panels c and d, greyed-out regions represent data that are not surely different from zero, based on its 95%-confidence interval estimated by bootstrapping (see Methods).

It should be noted that the thermodynamic coupling function has a natural normalization when the CVs are discrete. If the two CVs are maximally coupled, constraining one CV will fully constrain the other. Thus, at maximum coupling,and thusWe can then normalize (16) to this upper bound to define the normalized allosteric coupling (AC),The AC ranges from 1 to −1 and matches the convention commonly used for positive and negative allostery; positive values indicate that constraining one CV reduces the free energy required to constrain the other, whereas negative values indicate that constraining one CV increases the free energy required to constrain the other. In essence, the magnitude of the AC describes what fraction of the maximal allostery is contributing to the free energy of the joint state, whereas the sign of the AC describes whether that allostery is positive or negative.

When applied to the biophysical information transmission process occurring in a receptor (i.e., the thermodynamic coupling between the ligand binding site and the active site), the definitions above indicate that the thermodynamic coupling function is negative if measuring the active site to be in the “active” state reduces the uncertainty associated with whether or not the ligand binding site is in the “bound” state. The two-state ligand-induced receptor activation model defined in eq can be just as easily described using the thermodynamic coupling function if the collective variables X(r⃗) and Y(r⃗) are defined to take only two discrete values (bound/unbound and on/off, respectively). In this context, the two-state allosteric efficacy in eq can be calculated from an allosteric cycle composed of four allosteric perturbation cycles:so thatThus, a large negative ΔΔAcycle for the allosteric coupling between activation and ligand binding indicates that when the receptor is in the active state, the uncertainty that a ligand is bound is greatly reduced, whereas when the receptor is in the inactive state, the uncertainty that a ligand is not bound is greatly reduced.

Relationships to the Mutual Information and the Copula

The mutual information, which is often used to quantify allostery,[14−17] is defined asInterestingly, ΔΔA(x,y) in eq 16 is proportional to the logarithmic term in 22, which is known as the pointwise mutual information (PMI),[18]Like the allosteric efficacy, the PMI is symmetric (i.e., the order of variables x and y does not matter). To understand the PMI from the perspective of information theory, one can consider the information gained due to the reduction in uncertainty associated with measuring a variable. This information gain by measuring X(r⃗) to be equal to x isHowever, if two variables are measured, and those variables are dependent on each other, the amount of information gained by measuring the second variable will be different from the amount that would be gained if it were measured alone. For example, if Y(r⃗) was measured to be y, the probability distribution of X(r⃗) is now conditioned on Y(r⃗) = y, and thus the information gained by measuring X(r⃗) becomesThe PMI is the difference in the information gain,The mutual information is the PMI weighted by the joint probability density function. Consequently, the mutual information gives a high weight to the thermodynamic coupling of perturbations of high equilibrium probability states and low weight to those of low equilibrium probability. This is important for the mechanistic interpretation of allosteric couplings that are quantified only by their mutual information, as functionally significant perturbations do not necessarily drive the protein toward a region of its intrinsic CV space that is already high probability prior to perturbation. In fact, perturbations such as ligands generally drive the system away from the unbound equilibrium (e.g., where the inactive state is preferred to the active state), so the mutual information would end up giving larger weight to less functionally relevant states. In such cases, when considering only the protein’s degrees of freedom, the mutual information is not a good quantification of the intrinsic thermodynamic couplings that mediate the system’s response to ligand binding. Even in the simplest case of the allosteric coupling between ligand binding and activation as described in the ATSM, the mutual information between ligand binding and activation will depend on the affinity of the ligand, and will go to 0 as the affinity goes to either 0 or ∞, independent of the allosteric efficacy of the ligand. Therefore, we argue that it becomes preferable instead to analyze the entire 2-dimensional thermodynamic coupling surface, ΔΔA(x,y), which we call the “allostery landscape”, as it contains information regarding the allosteric efficacy for all possible perturbations to the distribution of those CVs.

The thermodynamic coupling function is also related to the copula density function from probability theory.[19] The copula density function of a bivariate probability distribution iswhere P(x) and P(y) are the cumulative distribution functions of X(r⃗) and Y(r⃗).

Any multivariate distribution can be expressed as a set of marginal probability distributions and a copula that defines the dependency between them,[19] and the entropy of the copula distribution is equivalent to the mutual information.[20] Thus, a multivariate thermodynamic coupling function behaves like a copula, defining the information transmission properties of the allosteric system. The relationship between the thermodynamic coupling function and fundamental concepts in information theory and probability theory suggest that past work in these fields may be able to be adapted for biophysical applications and provide new insights into allostery.

Contribution of Specific Energy Terms

Having introduced a quantification of the allosteric coupling between two CVs with the allostery landscape, a major mechanistic question still remains. What features of the structure and energetics of a given system define the thermodynamic coupling function? To answer this question, we derive the change in thermodynamic coupling function when a generic biasing potential energy term Ubias(r⃗) is added to the system’s total potential energy function, U(r⃗). The change in thermodynamic coupling, eq 16, at any point in the CV space can be estimated using a free energy perturbation approach, which we will refer to here as “biasing” to avoid confusion with “perturbing” that refers here to constraining the system at (x,y) in the CV space in eq . The change in free energy of the system when a biasing potential is added isThe change in free energy of the perturbed states can be similarly written asThus, the biased thermodynamic coupling function isWe wish to understand the contribution to the thermodynamic coupling function of a structural feature of interest, or of a specific interaction between structural elements of the system. Assuming this feature of interest is described by a specific energy term Uint(r⃗) of the total potential energy function U(r⃗), we can use eq with a biasing potential that is equal and opposite to that energy term, Ubias(r⃗) = −Uint(r⃗). In the next section, we use this approach to quantify the contribution of specific interactions in the alanine dipeptide system by mapping the corresponding change in thermodynamic coupling, ΔΔΔA(x,y) = ΔΔA(x,y) – ΔΔAbias(x,y).

Consider the special case of a particular energy term of interest that is a function of a CV, Z(r⃗). This corresponds to a biasing potential of the form Ubias(r⃗) = −Uint(Z(r⃗)). We have found (see Appendix) that this results in the following biased thermodynamic coupling functions:The second result is interesting because in this case ΔΔΔA(x,y) becomes a function of y only. A corresponding result is found if X(r⃗) and Y(r⃗) are permuted.

If Z(r⃗) is independent of either X(r⃗) or Y(r⃗),

If Z(r⃗) is conditionally independent of Y(r⃗) given X(r⃗), or if Z(r⃗) = X(r⃗),

Importantly, these findings indicate that the influence of any specific energy term on the thermodynamic coupling between two CVs can be clearly defined from the difference between the unbiased and biased thermodynamic coupling functions. The two conditions described above indicate that unbiasing by a potential energy term that mediates an allosteric coupling will have a two-dimensional effect on the thermodynamic coupling function (i.e., the second term in eq 30 is dependent on both X(r⃗) and Y(r⃗)). We know that, owing to nonadditive effects, the free energy contributions of such coupled energy terms cannot be rigorously deconvoluted.[21] Therefore, the contributions obtained by biasing the allosteric coupling function (shown in Figure ), cannot be taken as a sensu stricto decomposition of ΔΔΔA(x,y). This method nonetheless allows for the detailed analysis of the mechanism of allosteric coupling and can be applied generally across any system whose conformational ensemble can be sampled using methods such as molecular dynamics (MD).

Figure 2

Contributions of specific interactions to the allosteric coupling function of the alanine dipeptide. The contribution of the specific interactions identified for each panel was calculated according to eq , revealing the prominent role of bonded interactions between the termini and channel. (a) Contribution of the dihedral energy term corresponding to the definition of the angle Φ (C–N–Cα–C). The corresponding interaction energy as a function of Φ and Ψ, biased free energy surface, the biased allosteric coupling function, and the associated AC are shown in Figure S2. (b) Contribution of the nonbonded interaction energy between the termini. Additional plots are in Figure S3. (c) Contribution of the nonbonded interactions between the termini and the channel. Additional plots are in Figure S4. (d) Contribution of the bonded interactions involving atoms from both the termini and the channel. Additional plots are in Figure S5. In all panels, greyed-out regions represent data that are not surely different from zero, based on its 95%-confidence interval estimated by bootstrapping (see Methods section).

As we show below with the example of the alanine dipeptide in solution, MD can provide the data necessary for our formalism to be implemented for specific allosteric proteins of interest. The need for this type of data in order to describe allostery is foreshadowed by the discussion accompanying Weber’s introduction[22] of a multistate model of cooperativity in oligomers in 1972. He noted that a thermodynamic description of allostery based on the protein’s intrinsic conformational equilibrium (referred to as “tautomerization”) was “far more restricted in scope and less useful”[22] than a model based on its oligomerization equilibrium. Weber argued that “characterization of tautomerizations in molecules involving as large a number of potential degrees of freedom as proteins is confined to some structural detail that happens to be observable”[22] and that it appeared “unlikely than any method can furnish us with a precise determination of the tautomerization constants necessary to characterize completely the system”.[22] Indeed, in order to utilize the formalism derived above, it is necessary to observe a large number of degrees of freedom simultaneously with high accuracy and precision so that a free energy landscape can be estimated. Although such detail had not yet been described for proteins in 1972, the first molecular dynamics (MD) simulation of a protein was published[23] just five years later in 1977. Using MD simulations, we are now able to estimate the equilibrium distribution and corresponding free energy landscape of one or more CVs in many proteins, while simultaneously observing both the coordinates of all atoms in the system and their interaction energies, which is essential to the application of the formalism derived above.

Application: The Alanine Dipeptide

Estimating the Thermodynamic Coupling Function

To illustrate the use and utility of the thermodynamic coupling function, we analyzed the allostery landscape of the alanine dipeptide in solution. The alanine dipeptide is a popular model system for testing enhanced sampling and free energy methods as the entire system can be described well by only two CVs, the Φ and Ψ torsion angles along the bonds connecting the alanine Cα atom to the capped N- and C-terminus, respectively (see Figure a). Here, in analogy to larger allosteric systems, we consider that Φ captures the state of the N-terminal domain and Ψ the state of the C-terminal domain, and we ask the question of how the N-terminal and C-terminal domains of the protein are allosterically coupled. Despite the small size of the system, the irregular features of the Φ/Ψ free energy surface (see Figure b) indicate that these CVs are thermodynamically coupled in a nontrivial way, and thus the alanine dipeptide is an ideal model system for illustrating the power of the thermodynamic coupling function.

We constructed the 2-dimensional Φ/Ψ probability density function of alanine dipeptide in water from five independent 50 ns trajectories produced with driven adiabatic free energy dynamics[24,25] (see Methods). Following a protocol that we previously demonstrated to yield well-converged free energy surfaces up to 40 kJ/mol above the global minimum,[26,27] we reconstructed the free energy landscape shown in Figure b using the reweighted histogram estimator. To investigate which features of the alanine dipeptide free energy landscape are due to thermodynamic coupling between the two angles, we used eq to calculate the thermodynamic coupling function shown in Figure c. Significant allosteric couplings are evident in the regions of the left-handed α-helix (known as αL, which should not be confused with the symbol for allosteric efficacy) and the region labeled C7ax′ (shown on Figure b, slightly lower than the C7ax conformation stabilized in the gas phase[27]). This indicates that if Φ is driven to the (0 to 2 rad/0° to 120°) region, the transition of Ψ to the (0 to 2 rad/0° to 120°) and (−1 to −2 rad/–60° to −180°) regions becomes more favorable.

In the normalized AC landscape of the alanine dipeptide, calculated according to eq 19 and shown in Figure d, the αL and C7ax′, regions have couplings of around 0.4, indicating that a substantial amount of the maximal theoretically possible Φ/Ψ allostery contributes to these state’s stabilities. Thus, while these regions have a relatively low probability, our analysis suggests that the allosteric coupling accounts for the small but significant populations of αL and C7ax′ conformations that appear at equilibrium. In addition, there appears to be significant coupling present at the transition region between αL and PPII, and to a lesser extent between αR and C7ax′, indicating that these transitions may also be facilitated by allostery. We also see significant negative allosteric coupling in the high free energy regions, which indicates that an unfavorable thermodynamic coupling between the CVs contributes to the high free energy of these regions.

The mutual information, eq , between Φ and Ψ is 1.11 ± 0.01 nats (95% confidence interval from bootstrapping), or for better comparison to the thermodynamic coupling, 0.29 ± 0.03 kJ/mol at 300 K. It should be noted that these values are quite low compared to the numerous regions of high thermodynamic coupling and normalized AC (|ΔΔA(x,y)| > 6 kJ/mol, |AC(x,y)| > 0.3, see Figure ), and thus utilizing the mutual information alone understates the thermodynamic coupling between Φ and Ψ. Mapping the quantity summed over in eq (see Figure S1) shows that the major contributions to the mutual information come from very localized regions of the CV space. Thus, using a single number to quantify the coupling between Φ and Ψ misses the fact that the allostery landscape has significant regions of both negative and positive coupling. This can be important if, for example, one seeks to design a ligand that allosterically stabilizes a lower probability state.

The Influence of Specific Interactions on the Allosteric Coupling between Termini

To understand which structural features contribute to the thermodynamic coupling of Φ and Ψ, we decomposed the alanine dipeptide into three structural components, as shown in Figure a: (i) the N-terminus, which includes all atoms on the N-terminal side of the Cα carbon, (ii) the C-terminus, which includes all atoms on the C-terminal side of Cα, and (iii) the “channel”, which includes Cα as well as the hydrogen and methyl side chain bound to it. These three structural components can mediate the Φ/Ψ thermodynamic coupling through three different mechanisms: (i) direct nonbonded interaction of the termini (estimated with a dielectric constant ε = 60), (ii) indirect interaction of the termini through nonbonded interactions with the channel, and (iii) indirect interaction of the termini through bonded interaction with the channel. Thus, we estimated the potential energy Uint(r⃗) contributed by each of these groups of energy terms for each frame along the trajectories and reweighted the free energy landscape with an equal and opposite biasing potential Ubias(r⃗) according to eq , see Figures S2–S4. We then calculated the contribution of Uint(r⃗) to the thermodynamic coupling landscape ΔΔA(x,y) using eqs and 30, as represented on Figure in the form of ΔΔΔA(x,y) = ΔΔA(x,y) – ΔΔAbias(x,y). As a control, we also biased using the negative of the potential energy term corresponding exactly to the definition of the Φ(r⃗) CV, that is, the C–N–Cα−C dihedral angle. This reweighting resulted in one-dimensional variations of the thermodynamic coupling function along the Ψ axis, as expected from eq (see Figure a).

Interestingly, while chemical intuition may suggest that the direct interaction of the termini is the major mediator of the thermodynamic coupling, we find that the direct nonbonded interaction only contributes to the negative thermodynamic coupling surrounding the central forbidden region, as shown in Figure b. The only other significant change to the thermodynamic coupling function is one-dimensional bands at Ψ ∼ 1 and Ψ ∼ −1 rad, which indicates that the interaction itself may to some extent be indirectly coupled to Φ through its direct dependency on Ψ.

Figure d shows that the bonded interactions between the termini and the channel are the most significant contributors to both the positive and negative thermodynamic coupling between the termini, while the nonbonded interactions between the termini and the channel (Figure c) do not significantly contribute to the thermodynamic coupling. These results suggest that Φ and Ψ become thermodynamically coupled due to the energetics of the bonds, angles, and dihedrals composed of atoms shared between each terminus and the channel. For example, the Φ and Ψ dihedrals each share the angle formed by three central atoms. Different combinations of Φ and Ψ frustrate this central angle to different extents, leading to a thermodynamic coupling between the two. We however note that in the alanine dipeptide system, the energy terms described above are tightly coupled with each other, as well as with other energy terms (such as the internal bonded energy of the channel). Therefore, the contributions represented in Figure do not represent an exact decomposition of ΔΔA(x,y) and must be regarded as useful cues for the qualitative understanding of how allosteric coupling can be established between two domains of a molecular system.

Conclusions

We have derived a thermodynamic coupling function based on the allosteric efficacy that quantifies the allosteric coupling between two continuous or discrete CVs. We find that the thermodynamic coupling function is related to both the pointwise mutual information and the copula, and is best represented in the form of an allostery landscape, in units of free energy. Such a representation reveals the allosteric response to all possible perturbations of the CVs. We showed that the allostery landscape of the Φ and Ψ dihedral angles of the alanine dipeptide contains positive allosteric couplings that appear to stabilize the αL and C7ax′ conformations, and negative allosteric couplings that coincide with the high regions of the Φ/Ψ free energy landscape. On the basis of the formalism we developed, we were able to attribute features of this thermodynamic coupling function to specific interaction energy terms, thus allowing interpretation of the allosteric landscape. It is important to note that the criterion introduced here for determining whether a specific interaction mediates an allosteric coupling is more rigorous than our previous n-body information-based criterion.[14] While the 3-body information between three CVs is in fact a function of biased thermodynamic coupling functions (see Appendix, eq ), if Z(r) is conditionally independent of X(r) or Y(r) given the other CV, the 3-body information will be maximal. Consequently, the 3-body information does not permit to determine definitively whether Z(r) mediates a thermodynamic coupling between X(r) and Y(r), or if one of the CVs mediates a thermodynamic coupling between Z(r) and the other. Specifically, the 3-body information criterion will include some number of false positives (as we have previously described[14]), whereas all structural features that correspond to a potential energy term and have a two-dimensional influence on the thermodynamic coupling function can be considered to be effective mediators of the thermodynamic coupling.

The concepts developed here are very general and are applicable to larger molecular systems, provided enough sampling is available and the functionally relevant CVs are known. This second condition is especially noteworthy for cases in which a complete functional description involves multimolecular considerations. For example, to apply this method to a specific GPCR, one must first identify the CVs that best describe the agonist binding process and the G protein activation process. However, in order to understand agonist-induced activation of a G protein by a GPCR, one may need to also understand the thermodynamic coupling between the agonist binding CV and additional CVs that describe the process of G protein binding to the GPCR as well as the process of G protein activation. Our new theoretical formalism and its computational implementation remains applicable despite such complications, and can serve as a powerful tool in understanding the molecular mechanisms of the many proteins in which allostery is essential to biological function. It has the potential to identify novel allosteric sites that modulate functionally important reaction coordinates, and such capabilities can help achieve a large variety of end points. Examples include the design of novel therapeutic agents that allosterically modulate their specific targets in new ways, as well as the elucidation of allosteric mechanisms to guide the design of novel, synthetic allosteric proteins.

Methods

The alanine dipeptide (N-acetyl-alanine-N′-methyl amide) was modeled with the all-atom charmm22* force field[28] and solvated in explicit TIP3P water molecules.[29] charmm22* was chosen as it is able to reproduce an accurate alanine dipeptide free energy landscape without utilizing the CMAP[30] term used by other force fields. We chose to avoid force fields using the CMAP term as it induces a trivial thermodynamic coupling through a direct interaction between Φ and Ψ, rather than allowing it to emerge from separate terms of the traditional potential energy function. Molecular dynamics simulation were performed using the Charmm port[31] in the Gromacs 4.5 program[32] with particle-mesh Ewald[33] treatment of electrostatics and Lennard-Jones interactions switched off between 10 and 12 Å.

The systems were maintained at temperature T = 300 K with Nosé–Hoover chain thermostats.[34] Similarly to our previous study on dipeptides,[27] enhanced sampling was achieved with driven adiabatic free energy dynamics[24,25] (dAFED), also known as temperature accelerated molecular dynamics[35] (TAMD), implemented in the PLUMED plugin.[36] Two collective variables (CVs), defined as the backbone dihedral angles Φ and Ψ, were coupled (harmonic constant 1000 kJ/mol/rad2) to heavy fictitious particles (pseudomass 50 amu·nm2/rad2) held at temperature Ts = 600 K by generalized Gaussian moment thermostats (order 2).[37] Simulations were conducted in five independent replicates of 50 ns each after a standard equilibration phase starting with independent initial velocities. Free energy surfaces (FESs) in the (Φ,Ψ) plane were reconstructed[26] using the reweighted histogram smoothed with multivariate Gaussian kernel regression in Matlab (release 2015b, The MathWorks, Inc., Natick, Massachusetts, United States). A cutoff of 40 kJ/mol was used for the FESs, above which sampling was too poor for reliable surface estimation.

In principle, estimating an observable from a dAFED/TAMD simulation requires binning the observable values in the CV space, and reweighting each bin by a function of the FES at this point.[38] However, ΔΔA(x,y) in eq depends only on the probability mass at 300 K in the CV space, p(Φ,Ψ). This can be derived directly from the density obtained from the dAFED/TAMD simulation, padb(Φ,Ψ), by rescaling and renormalizing,Due to the surface smoothing steps, propagation of uncertainties is not practical for estimating confidence intervals on the allostery landscape. Instead, we use the bootstrapping approach.[39] Specifically, because observations from MD time series are notoriously not independent, we use block bootstrapping;[40] that is, we generate artificial samples by drawing at random (with replacement) segments of trajectory of 1 ns in length. Then, for each bin in the (Φ,Ψ) plane, we estimate a 95% confidence interval for the allosteric coupling function and for the AC based on the standard deviation among the bootstrapped samples. If in a given bin this confidence interval includes the value zero, the existence of an allosteric effect cannot be assessed with certainty in this bin and we represent it in a greyed-out color in panels c and d of Figure and in Figure .