Allosteric networks in thrombin distinguish procoagulant vs. anticoagulant activities

Results

We initially performed a set of six independent 20-ns CMD simulations for each of the three systems: thrombin, thrombin:TM56, and thrombin:TM456. For all three systems, a rather large conformational transition in the active site loops was observed during the equilibration procedure. Central to this conformational relaxation was a reorientation of both the extended 148_CT γ-loop (residues 178–192) and the 220_CT loop (residues 265–274), forming a more “open” active site pocket. Notably, the extended 148_CT γ-loop moves closer to the ABE1 70_CT loop (residues 93–109) (residue numbering details are provided in SI Text). This observation confirms the hypothesis that the presence of the inhibitor in the available X-ray crystal structures, which was removed before performing the CMD simulations, promotes a closed loop conformation. After this initial structural relaxation, the conformational dynamics of each of the three systems rapidly stabilized and the backbone rmsd to the average structure for the thrombin molecule was found to be 1.4 Å, 1.2 Å, and 1.2 Å for thrombin, thrombin:TM56, and thrombin:TM456, respectively.

The molecular ensembles generated from the CMD simulations were subjected to a community network analysis, which identifies communities (clusters of highly connected residues) based on residue-by-residue correlation and proximity (22). A measure of “betweenness” of the intercommunity edges (Methods) affords the identification and characterization of potential allosteric pathways. Community network analysis results for thrombin, thrombin:TM56, and thrombin:TM456 are shown in Fig. 2.

Several interesting and highly reproducible communities of residues were observed in the analyses in the majority of the CMD trajectories. In all cases, TM4, TM5, and TM6 are each represented by a single community. By contrast, the thrombin molecule in all systems consists of multiple interconnected communities divided between, most notably, the TM binding site (ABE1) and the proteolytic active site, including the catalytic triad (H57_CT, D102_CT, and S195_CT), the 60_CT insertion, and the active site loops [90_CT, 148_CT (γ), 170_CT, 186_CT, and 220_CT]. However, the number of thrombin communities differs from one system to the next. On average, isolated thrombin ostensibly comprised eight distinct large communities, whereas the thrombin molecule in thrombin:TM56 and thrombin:TM456 comprised only five to six large communities, on average. Similarly, the metric of cumulative betweenness of intercommunity edges, which is a direct measure of the strength of the potential allosteric interaction between communities, is greater in the two thrombin–TM complexes compared with thrombin alone. In the thrombin:TM56 system, a strong allosteric pathway is consistently observed linking the TM5, ABE1, and active site communities. The TM4 domain, which forms a strong allosteric connection to TM5, acts to enhance the strength of this allosteric pathway further. A second and very important difference between the community networks for thrombin:TM56 and thrombin:TM456 is the presence of a strong communication between TM5 and the 148_CT γ-loop of thrombin, which is only observed in the thrombin:TM456 construct.

The fast time-scale dynamics of thrombin in the isolated thrombin, thrombin:TM56, and thrombin:TM456 systems were analyzed by calculating backbone NH^N order parameters, which provide a quantitative measure of the extent of (internal) reorientation dynamics of the N-H bond vector per residue (Fig. 3A). Similar to the recent NMR/molecular dynamics simulation study on thrombin:PPACK, a strongly heterogeneous distribution of fast time-scale motions was observed in all three systems (21). In the isolated thrombin system, order parameters in the 30_CT loop (residues 53–61), which forms part of ABE1, and in the 60_CT insertion (residues 83–91) were as low as 0.66 and 0.64, respectively. By contrast, in both the thrombin–TM complexes, these order parameters were substantially higher, with the lowest value being 0.76. Reduced fast time-scale dynamics (and hence higher order parameters) in the 30_CT loop are readily explained because this region makes direct contact with the TM5 domain. However, the same cannot be said for the 60_CT insertion. Significantly, the community network analysis models for both thrombin:TM56 and thrombin:TM456 (Fig. 2) identified a strong allosteric pathway connecting TM5 to the active site of thrombin via the 60_CT insertion. Higher order parameters in these same loops were also observed experimentally (21). The dynamic communication observed between the active site, the 30_CT loop, and the 60_CT insertion is a prime example of entropic allostery. No other significant differences in the fast time-scale NH^N order parameters for the three systems were observed. Calculated backbone rms atomic fluctuations also revealed no other significant differences, and the average structures of thrombin for the three systems were remarkably similar (backbone rmsd < 0.5 Å).

To access microsecond-millisecond time-scale motions more typical of allosteric regulation, a series of AMD simulations was performed at two acceleration levels: a moderate acceleration level (level 1), which samples configurational dynamics up to time scales of several hundred nanoseconds, and an aggressive acceleration level (level 2), probing dynamics up to several tens of microseconds (details are provided in SI Text). After free energy weighting of the AMD trajectories, a representative ensemble of structures was obtained. Whereas no substantial domain reorientation was observed for the TM EGF domains in the CMD, the AMD trajectories showed some domain reorientation, although with limited extent of motion. From the free energy-weighted trajectories, averaged thrombin NH^N order parameters were calculated (Fig. 3A). In contrast to the fast time-scale dynamics, marked differences in the order parameters were observed for the three systems, particularly at the most aggressive acceleration level (Fig. 3B). Both TM-bound systems exhibit considerably less dynamics in the 30_CT loop (residues 53–61) and the 60_CT insertion (residues 83–91), consistent with what was observed for the fast time-scale order parameters. However, on slow time scales, thrombin:TM456 exhibits remarkably more re-orientation dynamics in the active site loop regions: the 90_CT loop (residues 122–132), the 148_CT loop (residues 184–192), the 170_CT loop (residues 213–221), the 186_CT loop (residues 228–234), and the 220_CT loop (residues 265–274). By contrast, thrombin:TM56 shows either no more or, in some regions, even less slow time-scale dynamics compared with the isolated thrombin system.

A residue-by-residue cross-correlation analysis (25) (Methods) applied to the free energy-weighted molecular ensembles obtained at the most aggressive acceleration level also revealed striking differences between the three systems (Fig. 4). In all three systems, we observed cross-correlated motion between the 30_CT loop and the 60_CT insertion, and the strength of this cross-correlation was markedly enhanced in the thrombin:TM456 system. Additionally, the thrombin:TM456 system exhibited a complex pattern of strongly correlated motions across the active site loops (Fig. 5) directly coinciding with the observation of enhanced slow time-scale dynamics in the level 2 AMD trajectories of the thrombin–TM456 complex (Fig. 3B). In contrast, correlations between the corresponding regions in both the thrombin:TM56 and the isolated thrombin systems were very weak or absent.

To confirm that the pattern of strongly correlated motions observed in the thrombin–TM456 complex was due to TM4, we extended our analysis to include the TM domains. To circumvent complications from the reorientation dynamics of the TM domains, multiple cross-correlation analyses were performed after superposing the molecular ensemble on different regions of the system, and the complete cross-correlated dynamics map (Fig. 6) was then reconstructed. The results clearly show that the complex network of correlated motions observed in thrombin:TM456 propagates from the TM4 domain through the TM5 domain and into the thrombin molecule (Fig. 6). Global reorientation dynamics of TM4 contribute significantly to the entire network of correlated motion by inducing global reorientation dynamics of TM5, and hence the specific orientation and interactions of TM5 with thrombin. Additionally, interdomain, residue-by-residue correlations are observed between regions of TM4 (residues 359–369) and TM5 (residues 404–414). Both global domain correlations and specific residue correlations in TM extend into ABE1 of thrombin, including the 30_CT loop (residues 53–61), the 60_CT insert (residues 83–91), and the 70_CT loop (residues 93–109). TM5 residues (404–414) are strongly correlated because they lie in close proximity to the 30_CT loop. In addition, Asp416 and Asp417 from TM5 make a significant contact with Thr105/Arg106 in thrombin. All these TM residues were identified previously as critical by alanine scanning mutagenesis (13). The network of correlated motion continues from ABE1 into the thrombin active site loops and remains especially strong in the 148_CT loop (residues 180–190) and includes the 90_CT loop (residues 122–132) and the 170_CT loop (residues 210–220).

Discussion

The work presented in this paper provides unique and detailed insight into the functional dynamics of the thrombin:TM system, which not only rationalizes the available experimental data but explicitly identifies and defines the important role played by the TM4 domain in the activation of protein C. The results show two distinct but interrelated allosteric pathways that mediate the anticoagulant activity of the thrombin–TM complexes. One pathway connects the TM5 domain with ABE1, the 60_CT insertion, and the active site. This allosteric network, which exists in both the thrombin:TM56 and thrombin:TM456 systems, was previously identified using H/D exchange (16). Community network models further indicate that TM4 strengthens this allosteric pathway (Fig. 2), leading to a significant enhancement in the cross-correlated motion between the 30_CT loop and the 60_CT insertion in the thrombin:TM456 system compared with the isolated thrombin and thrombin:TM56 systems (Fig. 4). Indeed, binding of TM56 to thrombin stabilizes some thrombin active site loops on slow time scales, as shown by the backbone NH^N order parameters calculated from AMD simulations (Fig. 3B).

The addition of the fourth EGF-like domain (TM456) activates slow time-scale dynamics in the active site loops of thrombin via another allosteric pathway. Global reorientation and intradomain configurational dynamics of TM4 are strongly coupled to both the global and intradomain dynamics of TM5. The importance of several of the most highly cross-correlated TM residues in these regions was previously characterized in an NMR study, which found that backbone dynamics of Tyr358 and Gln359 are correlated with anticoagulant activity and that backbone dynamics of Tyr413 and Ile414 are inversely correlated with thrombin binding (24). This tight coupling between TM4 and TM5 was previously shown experimentally to depend on M388 in the TM4-TM5 linker (25, 26). The specific orientation and interactions of TM5 induced via coupling to TM4 mediate specific interactions with ABE1. The dynamics of ABE1 are, in turn, highly correlated with the dynamics of the thrombin active site loops only in the TM4-containing construct. Through this complex pattern of extended correlated motions, slow time-scale configurational dynamics of TM4 are directly linked to the thrombin active site, providing evidence for entropic allostery between TM4 and thrombin (Figs. 5 and 6). The community network analysis also reveals that when TM456 is bound, the information transfer between the loop regions of thrombin is significantly enhanced, resulting in the observation of fewer yet larger communities and increased potential for communication between communities. This enhanced communication is particularly strong for the path connecting the active site loops to TM4 via ABE1 and TM5 and, interestingly, includes a direct allosteric interaction between the 148_CT loop and TM5. The highly correlated motions between the active site loops on slower time scales observed in thrombin:TM456 overlap directly with the different communities identified in the community network model. In particular, the 90_CT loop coalesces with the active site into a single community consistent with H/D exchange experiments that identified subtle changes in the 90_CT loop only in the presence of TM456 (15). Our observation of different slow time-scale dynamic behavior of the active site loops of thrombin in thrombin:TM456 vs. thrombin:TM56 provides a mechanism for these heretofore inexplicable experimental results.

The slow conformational loop dynamics, mainly of the active site loops of thrombin, which only occur in thrombin:TM456, probably facilitate the association of small substrates and inhibitors (that may enter the active site without interacting directly with TM4), thereby increasing the association rate by several orders of magnitude. For the specific case of the significantly enhanced k_a of protein C, two complementary mechanisms exist. First, following the work of Fuentes-Prior et al. (10), the presence of TM4 forms an extended binding surface for protein C, providing optimal alignment for insertion into the active site and subsequent cleavage. This docking and optimal alignment mechanism probably works together with the allosteric mechanism identified in this study, in which the altered dynamics of the active site loops caused by TM binding increase the protein C k_a. The binding of protein C to TM4 may also affect dynamics within TM4 to enhance the allosteric network and further promote protein C cleavage.

Methods

AMD.

The details of the AMD protocol have been discussed previously (22, 23), and only a brief summary is provided here. In AMD, a reference or “boost energy,” E_b, is defined, which lies above the minimum of the potential energy surface (PES). At each step in the simulation, if the instantaneous potential energy, V(r), lies below the boost energy, a continuous nonnegative bias potential, ΔV(r), is added to the actual potential. If the potential energy is greater than the boost energy, it remains unaltered. The application of the bias potential results in a raising and flattening of the PES, decreasing the magnitude of the energy barriers and thereby accelerating the exchange between low-energy conformational states, although still maintaining the essential details of the potential energy landscape. Explicitly, the modified potential, V^∗ (r), on which the system evolves during an AMD simulation is given by (22):

and ΔV(r) is defined as:

The extent of acceleration is determined by the choice of the boost energy, E_b, and the acceleration parameter, α. Conformational space sampling can be enhanced by either increasing the boost energy or decreasing α. In the present work, we have implemented a “dual boost” AMD approach (29), in which two acceleration potentials are applied simultaneously to the system: The first acceleration potential is applied to the torsional terms only, and a second, weaker acceleration is applied across the entire potential. For each of the three systems, thrombin, thrombin:TM56, and thrombin:TM456 dual-boost AMD simulations were performed at two (torsional) acceleration levels. For the most aggressive AMD simulations (level 2), the specific acceleration parameters were defined as E_b(dih) − <V₀(dih)> = [4 kcal/mol * No. residues], and the acceleration parameter, α(dih), was set to one-fifth of this value. The specific choice of these AMD parameters was based on the recent NMR/AMD study of thrombin:PPACK (details are provided in SI Text) (21). For the second, moderate AMD simulations (level 1), the α(dih) parameter was kept the same and E_b(dih) − <V₀(dih)> was reduced to [2 kcal/mol * No. residues]. In all AMD simulations, the total background acceleration parameters were fixed at E_b(tot) − <V₀(tot)> = α(tot) = [0.16 kcal/mol * No. atoms in simulations cell]. For each of the three systems, six AMD simulations were performed at both acceleration levels for 10,000,000 steps (the equivalent of 20-ns standard molecular dynamics). The physical conditions and computational parameters used in the AMD simulations were identical to those defined above for the CMD simulations, and all molecular dynamics simulations were performed using an in-house modified version of the AMBER10 simulations suite (30). For each AMD trajectory, a corrected canonical ensemble was obtained by performing a Boltzmann free energy reweighting protocol using the bias potential block averaging method to remove statistical noise errors (details are provided in SI Text). In this way, six representative free energy-weighted molecular ensembles were generated at both acceleration levels, along with the six unbiased 20-ns CMD simulations for all three systems.

Trajectory Analysis.

Allosteric networks were characterized using a community network analysis approach previously applied to investigate allostery in tRNA–protein complexes and other protein systems (24, 31, 32). This approach constructs a dynamic contact map consisting of a network graph in which each residue is treated as a “node” connected by edges to other nodes when two residues are deemed to be “in contact.” The dynamic contact map is subsequently decomposed into communities (i.e., clusters of residues) of highly intraconnected but loosely interconnected nodes using the Girvan–Newman algorithm (33). Central to this method is calculation of edge betweenness, the number of shortest paths that cross an edge. The edge betweenness is calculated for all edges, and the edge with the greatest betweenness is removed. This process is repeated, and a modularity score is tracked to identify the division that results in the optimal community structure. Network graph calculations were performed using the python module NetworkX (33).

Residue-by-residue cross-correlations were calculated using the generalized cross-correlation approach applied to all backbone Cα atomic coordinates based on the mutual information method developed by Lange and Grubmüller (25) using the g_correlation module in GROMACS 3.3.3 (34).

The internal dynamics were monitored by calculating backbone NH^N order parameters (S²) from the different CMD and AMD simulations, which provide a quantitative measure of the extent of reorientational motion of the given bond vector (35). In all cases, molecular ensembles generated from the standard CMD simulations and the free energy-weighted AMD trajectories were superposed onto the heavy backbone atoms of all residues for the appropriate average structure. Order parameters were calculated as follows (36):

where μ_i are the Cartesian coordinates of the normalized internuclear vector of interest. The resulting order parameters were then averaged over all molecular dynamics/AMD trajectories.

Supplementary Material