Hamiltonian Mapping Revisited: Calibrating Minimalist Models to Capture Molecular Recognition by Intrinsically Disordered Proteins.

Molecular recognition by intrinsically disordered proteins (IDPs) plays a central role in many critical cellular processes. Toward achieving detailed mechanistic understanding of IDP-target interactions, here we employ the "Hamiltonian mapping" methodology, which is rooted in the weighted histogram analysis method (WHAM), for the fast and efficient calibration of structure-based models in studies of IDPs. By performing reference simulations on a given Hamiltonian, we illustrate for two model IDPs how this method can extrapolate thermodynamic behavior under a range of modified Hamiltonians, in this case representing changes in the binding affinity (Kd) of the system. Given sufficient conformational sampling in a single trajectory, Hamiltonian mapping accurately reproduces Kd values from direct simulation. This method may be generally applied to systems beyond IDPs in force field optimization and in describing changes in thermodynamic behavior as a function of external conditions for connection with experiment.

Molecular recognition involving intrinsically disordered proteins (IDPs) is a key component in signaling and regulation networks. Moreover, IDPs are increasingly implicated in human diseases, such as cancer and neurodegeneration.[1] Delineating the mechanisms by which IDPs interact with their biological partners thus has the potential to expand our knowledge of protein interactions in the cell and to treat human diseases.

Biomolecular simulation methods provide a means to characterize the atomic-level details of IDP recognition. However, the long biological time scale of binding events is beyond the reach of conventional and most enhanced sampling simulation approaches. Recently, this issue has been addressed through the application of structure-based models to study IDP interactions.[2−10] Since IDPs commonly fold into a stable structure upon interaction with a target protein, in these studies, an experimental structure of an IDP–target protein complex serves as a reference for defining a native residue–residue contact map. The success of structure-based models in characterizing IDP recognition requires a delicate balance between the strength of the interactions governing both the native intermolecular binding and intramolecular folding contacts.[4] A reasonable balance of forces is typically achieved by tuning the strength of the contacts such that the binding affinity (Kd) and secondary structure content from simulation match those values from experiment.[4,8] While tuning the model is essential for capturing the proper mechanism of molecular recognition, it often necessitates running many additional simulations until the model is properly calibrated, which can be both computationally expensive and time-consuming.

In 1998, Brooks and co-workers[11] developed a framework for exploring the space of Hamiltonians in the “vicinity” of a given reference system in the context of investigating the folding free-energy landscape of a coarse-grained model. This framework builds upon the weighted histogram analysis method (WHAM) of Swendsen[12] and has been applied in a number of related contexts to examine the influence of denaturants[13−15] and pH,[16] as well as other factors that affect the folding landscape of proteins and peptides.[17] In the present work, we utilize for the first time, to our knowledge, this formalism to explore a timely question of tuning a coarse-grained force field to recapitulate binding affinities between two protein partners. The approach is general and can also be applied to higher-order complexes.[8] Given sufficient conformational sampling, this approach extrapolates changes in the thermodynamic behavior of a system under modified Hamiltonians from just a single reference trajectory. We illustrate the utility of “Hamiltonian mapping” by comparing theoretically computed Kd values to direct simulations and to experiment for two model IDPs (Figure 1). The approach is general and can be used to study molecular recognition in other IDP systems and to help advance other biophysical modeling efforts.

Figure 1

KIX–c-Myb (left) and MLL–KIX (right) binary complexes (PDB ID: 2AGH).

Following recent work,[4,8] we performed reference simulations of the KIX–c-Myb and MLL–KIX binary complexes[18] (Figure 1) using an intermolecular scaling factor, λ, equal to 1.1 and at three different time lengths (1.5, 15, and 450 μs). The intermolecular interactions from these reference simulations were then reweighted using the Hamiltonian mapping approach (see the Computational Methods section) by increasing λ. For each reweighted ensemble (i.e., representing a given λ value), we computed the binding affinity (Kd) and compared this value with the corresponding Kd value from 450 μs long direct simulations with the same scaling factor. Overall, Figure 2 shows that as λ increases, binding becomes tighter (lower log Kd values). Kd values calculated from Hamiltonian mapping using a λ = 1.1 reference ensemble matched the Kd values calculated from the direct simulations (Figure 2A–F). We also examined the scenario in which we increased λ to 1.5 in our reference simulation and decreased λ for the mapping (Figures 2G–I). With the exception of the 1.5 μs long simulation (Figure 2G), which failed Hamiltonian mapping due to the simulation being too short to sample unbound states, we also obtained accurate results in this case. In general, given sufficient sampling, the Hamiltonian mapping method was able to accurately reweight Kd values spanning several orders of magnitude and provide a reliable estimate of λ to quickly and efficiently match simulated Kd values to experiment.

Figure 2

Binding affinities for (A–C) KIX–c-Myb and (D–I) MLL–KIX extrapolated from Hamiltonian mapping (black boxes) compared with reference simulations at each λ value. The columns are organized by the time length of the reference simulations, which are indicated above the first row. In each panel, the dotted vertical line denotes the reference simulation to which Hamiltonian mapping was applied, and the gray horizontal bars correspond to the experimental Kd value (10 ± 2 μM for KIX–c-Myb and 2.8 ± 0.4 μM for MLL–KIX).[19] Note that the simulated Kd values in (D–F) are identical to those in (G–I), but the reference simulations for Hamiltonian mapping are different.

In the earliest work on the subject by Brooks and co-workers,[11] it was demonstrated that the thermodynamics of a large range of modified Hamiltonians could be extrapolated by simply reweighting one simulation using a revised expression of the WHAM equations. As noted previously and reiterated again here, it is important that the sampling in the reference simulation adequately covers the configuration space of interest. In other words, the sampling in the reference simulation must overlap significantly with the sampling of the modified Hamiltonian in order for the extrapolation method to be effective. This is generally not a problem with much longer coarse-grained simulations using modern computers. In all of the cases presented in the current work, the 15 μs long reference simulations, which took about 1 week to perform on a single CPU, provided a reliable estimate of Kd following reweighting. Moreover, it is important to point out that Hamiltonian mapping is not strictly limited to reweighting a single reference simulation as presented here. In fact, the sampling from multiple reference simulations with different λ scaling factors could be trivially combined and reweighted simultaneously using the modified WHAM equations. For example, we assume the scenario in which λ = 0.5 is the scaling factor that yields a Kd in agreement with the direct simulation and that two independent reference simulations produce ensembles that are either predominantly bound (λ = 1) or unbound (λ = 0). Then, instead of reweighting either reference simulation alone, which would produce a poor estimate of λ due to insufficient overlap in sampling, one could combine information from both reference simulations for the mapping to obtain an accurate estimate of the Kd.

To our knowledge, this is the first systematic study that demonstrates the true efficacy of the Hamiltonian mapping method through direct comparisons with reference simulations where the phase space was exhaustively sampled. We illustrate the utility of the method by examining the thermodynamics of binding in coarse-grained simulations of two model IDP systems. Given sufficient conformational sampling, Kd values extrapolated from Hamiltonian mapping are in excellent agreement with those achieved from direct simulation. Furthermore, we find that a moderate simulation time with the coarse-grained model and low computational cost are needed to obtain a reliable reference ensemble for Hamiltonian mapping, such that this approach could be practically applied to the rapid and efficient calibration of other IDP systems. We anticipate that the quantitative assessment of binding thermodynamics afforded by the method will play a valuable role in characterizing mechanisms of IDP recognition.

Beyond studies of IDPs, Hamiltonian mapping has the potential to facilitate other biophysical modeling efforts. For instance, the mapping procedure could significantly expedite the parameterization of force fields[4,8] (both coarse-grained and all-atom) or any study in which the Hamiltonian is being modified (e.g., Hamiltonian replica exchange simulations[20−23]). The approach also presents a practical tool for modeling changes in thermodynamics under varying external conditions (e.g., temperature, pH, denaturant concentration, etc.).[14−17] The unification of such efforts within the framework of Hamiltonian mapping in studies of protein interactions and folding would provide a robust link between simulation and experiment.

Computational Methods

Hamiltonian Mapping Formalism. Hamiltonian mapping[11] demonstrates how the standard WHAM formalism[12] can be modified and used to extrapolate information about a biased Hamiltonian from the unmodified reference Hamiltonian, H0(x), where x represents the molecular coordinates. WHAM is typically used along with umbrella sampling,[24] in which a set of L additional restraining potentials, {V} = V1(x), ..., V(x), are added to H0(x) in order to enhance conformational sampling in rarely visited regions of phase space. In this case, the modified Hamiltonian takes the formwhere {λ} is the set of coupling parameters used for scaling individual biasing potentials, λ0 = 1, and H0(x) = V0(x). Then, for R independent simulations, each performed at a temperature T = 1/kBβ (where kB is the Boltzmann constant) and with differing {λ}, the bias is removed by solving the following pair of WHAM equations self-consistently[12]andT may be constant or vary across the R simulations. (In the current work, T = 300 K for all simulations.) ξ is a particular value along a progress variable of interest, n is the total number of configurations in the mth simulation, and N({V},ξ) is the histogram count of configurations with {V} and ξ. Finally, f is related to the Helmholtz free energy for the mth simulation.[11,12]

In contrast to standard WHAM, in the Hamiltonian mapping framework, we proceed in the opposite direction; we perform one or more simulations of the system under the original Hamiltonian H0(x) and then add a bias to the statistical averages computed from H0(x) to examine thermodynamic behavior under a modified Hamiltonian. For this scenario, eq 2 is recast asEssentially, the histogram counts N now come from the sampling of H0(x), and the sum of the scaled restraining potentials in the denominator reduces to H0(x). Equations 4 and 3 are then iterated as above for standard WHAM.

Coarse-Grained Potential. To model binding and folding events of the KIX–c-Myb and MLL–KIX binary complexes, we employ the sequence-flavored Go̅-like model of Karanicolas and Brooks.[25] Briefly, in this model, each residue is represented by a single bead centered at the Cα position and with the mass of the corresponding amino acid. An additive potential describes the bonded and non-bonded interactions between the beads. Virtual bonds and angles are defined by harmonic potentials with reference values determined by the Cα coordinates in the experimental structure. A potential for virtual dihedral angles is based upon backbone dihedral angle probability distributions from the Protein Data Bank (PDB) for the 400 possible amino acid pairs and thus is independent of the specific topology of the system. Residue pairs separated by three or more bonds interact through the following potentialwhere r is the distance between residues i and j, σ corresponds to the distance between the two residues at which V is minimum, and ε is the interaction strength at σ. Residues located in close proximity in the experimental structure (i.e., forming a “native contact”) interact favorably with ε values based on the statistical contact energies of Miyazawa and Jernigan,[26] while residue pairs not in close proximity experience a slight repulsive interaction that takes the form of a typical 12–6 Lennard-Jones potential. Further description of the Go̅-like model and details of the simulation setup can be found in refs (25) and (8), respectively.

Tuning Intermolecular Binding Forces. We further consider the non-bonded interaction potential (eq 5) for native contacts between the IDP and the receptor in tuning the model to match the experimental Kd. For this set of intermolecular contacts, we scale the interaction strength, ε, by a factor λ so as to control the proportion of bound and unbound states throughout simulation. From the fraction of unbound states, pu, the Kd is computed aswhere [protein] is the concentration of protein in moles per liter (10–3 mol/L in the current study).[4,8] A configuration is considered unbound if zero intermolecular native contacts are formed and if the distance between the centers of mass of the IDP and of the target protein is greater than the cutoff for calculating the non-bonded interactions (25 Å). We also note that the strength of intramolecular contacts within the IDP are tuned in a similar manner to reproduce the experimental estimate of residual helical structure in the unbound state.[8]