thrombin is the central 
 protease in the cascade of blood coagulation 
 proteases . structurally , thrombin consists of a double -barrel 
 core surrounded by connecting loops and helices . genetic analysis 
 of the clotting factor genes demonstrates that the clotting proteases 
 of the chymotrypsinogen superfamily have evolved as a result of several 
 gene duplications , exon shuffling , and intron sliding events . prothrombin 
 has a unique exon organization and is thought to be the ancestral 
 gene of the clotting factor family . the 
 extended active site loops in thrombin are thought to have arisen 
 from insertions in the serine protease that evolved to impart greater 
 specificity . thrombin is produced in a low - activity 
 zymogen form that requires proteolytic cleavage to attain full activity . this cleavage event results in small overall changes to the molecular 
 architecture but results in a large change in dynamics wherein one 
 -barrel becomes more dynamic and the other becomes less dynamic . the result is a more perfectly formed active 
 site for rapid proteolytic cleavage activity . despite the highly specific 
 nature of thrombin activity , in association with allosteric modulators , in addition , 
 allostery is key to thrombin regulation , and misregulation 
 can lead to bleeding disorders or thrombosis . although , traditionally , 
 allostery was defined as occurring among subunits in a multisubunit 
 system such as hemoglobin , the phenomenon 
 of altered activity resulting from binding of a regulatory molecule 
 on the opposite side of a monomeric enzyme is now also recognized 
 as a form of allostery . several experimental and 
 computational approaches have hinted that 
 the solution structure of thrombin is a broad and malleable dynamic 
 ensemble . h / d exchange mass spectrometry showed that d - phe - pro - arg chloromethylketone ( ppack ) occupation of the active 
 site not only protected the active site loops but also propagated 
 to decreased exchange in several regions of the protein distant from 
 the active site . thrombin has two binding 
 sites distal to the active site ; exosite 1 is where thrombomodulin 
 binds , and exosite 2 is where heparin binds . isothermal titration 
 calorimetry ( itc ) experiments showed alteration of thermodynamic parameters 
 of ligands binding to thrombin exosites when the active site was occupied . binding of active site ligands altered the balance 
 of enthalpic and entropic contributions to binding of exosite 1 ligands 
 and vice versa . such thermodynamic compensation 
 phenomena are more likely if the allosteric mechanism is entropic 
 rather than enthalpic , suggesting that differences in the dynamic 
 properties of the system affect the ligand binding mechanism . indeed , x - ray crystallography shows no significant changes in the 
 thrombin structure upon ligand binding , providing further evidence 
 that it may exist as a malleable dynamic ensemble . nmr studies 
 and md simulations remain the most direct approaches 
 to investigating protein dynamics . a 
 recent nmr / md study on ppack - thrombin revealed a large degree of dynamic 
 motions , particularly in the active site loops , spanning time scales 
 from picoseconds to milliseconds . a computational 
 exploration revealed a strong dynamic component in the allosteric 
 regulation of thrombin by tm . in the 
 work presented here , we combine conventional md , nmr - calibrated accelerated 
 md ( amd ) , and analysis of the residual local frustration to further 
 explore the dynamics of thrombin , with particular interest in changes 
 that occur upon active site ligation . the combined approach allows 
 us to analyze a broad range of motional time scales . atomic coordinates 
 for thrombin were obtained from the protein data bank 1.9  x - ray 
 crystal structure [ pdb i d : 1ppb ] . both systems were placed 
 at the center of a periodically repeating box , and the simulation 
 cell size was defined such that the distance between the edge of the 
 simulation box and the surface of the solute was at least 12  . all simulations were performed in explicit solvent , and three cl counterions were introduced to obtain cell neutrality . six 20 ns conventional md simulations were performed using a different 
 random seed generator for the maxwellian distribution of atomic velocities 
 following standard energy minimization and equilibration procedures . electrostatic interactions were treated using the particle mesh ewald 
 ( pme ) method with a direct space sum 
 limit of 10  . the ff99sb force field was used for the solute residues , and the tip3p water force field was employed for the solvent molecules . in the 
 case of ppack - thrombin simulations , the conventional md simulations 
 were analyzed in the community analysis and also constituted the starting 
 point for the amd simulations . these simulations also provided the 
 average ( unbiased ) dihedral angle energy , v0(dih) and total energy v0(tot) values used to define the acceleration 
 parameters in the amd simulations described below . allosteric 
 networks were characterized using a community network analysis approach 
 previously applied to investigate allostery in trna / protein complexes 
 and other protein systems . this approach constructs a dynamic 
 contact map consisting of a network graph in which each residue is 
 treated as a two residues are deemed to be  in contact  throughout 
 the majority of the simulation . clusters of residues ) of highly intraconnected 
 but loosely interconnected nodes using the girvan central to this method is the 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 . amd 
 simulations were performed 
 as described previously using an in - house modified version of the 
 amber 10 code . amd approach is used , 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 . this dual boost 
 amd protocol represents a unified approach facilitating the efficient 
 sampling of both the torsional degrees of freedom and slow diffusive 
 motions in the solute . in total , six independent dual boost amd simulations 
 were performed for 10  000  000 steps ( the equivalent 
 of 20 ns of md ) for each system . the physical conditions , force fields , 
 and all other simulation parameters employed were identical to those 
 described for the conventional md simulations . the specific acceleration 
 parameters used in this study were [ eb(dih )  v0(dih) ] 
 = [ 4 kcal / mol  no . of solute residues ] and (dih ) = [ 0.8 
 kcal / mol  no . v0(tot) ] = (tot ) = [ 0.16 kcal / mol 
  no . of atoms in the simulation cell ] for the background total 
 acceleration . these acceleration parameters had been previously identified 
 as the optimal acceleration parameters for the reproduction of experimental 
 rdcs in ppack - thrombin , accessing configurational dynamics on time 
 scales up to 10s100s of microseconds . for each amd simulation , a corrected canonical ensemble was obtained 
 by reweighting each point in the configuration space on the modified 
 potential by the strength of the boltzmann factor of the bias energy , 
 exp[v(rt(i ) ) ] at that particular point , and the bias potential block averaging 
 method was employed to remove statistical noise errors . the internal dynamics present in the different 
 amd simulations of apo - thrombin and ppack - thrombin were assessed by 
 calculating order parameters , s , from 
 the free - energy weighted amd ensembles . members of each ensemble were 
 superimposed onto the backbone atoms ( n , c , c ) 
 of all heavy chain residues for the appropriate average structure , 
 and order parameters , s were calculated 
 aswhere i are the cartesian coordinates of the normalized 
 internuclear vector 
 of interest . others have shown that s values calculated from standard md simulations in this way were 
 in excellent agreement with experimental s values calculated using the lipari the order parameters 
 presented here are averaged over all six amd trajectories for each 
 system . residue - by - residue cross - correlations for the free - energy - weighted 
 amd ensembles were calculated using the generalized cross - correlation 
 approach applied to all backbone c atomic coordinates 
 based on the mutual information method developed by the grubmller 
 group using the g_correlation module 
 in gromacs 3.3.3 . an algorithm for 
 determining residual local frustration , that is , whether a contact 
 between amino acid residues is energetically optimized or not in the 
 folded state , was developed by the wolynes and komives groups some 
 time ago . residue 
 interactions by systematically perturbing the identity of individual 
 residues and evaluating the resulting total energy change . for the 
 work presented here index , in which the decoy set involves randomizing not just the identities 
 but also the distance and densities of the interacting amino acids i , j. this scheme effectively evaluates 
 the native pair with respect to a set of structural decoys that might 
 be encountered in the folding process . after constructing a histogram 
 of the energy of the decoys and comparing the distribution to the 
 native energy , cutoffs energetically favorable contacts between 
 residues are minimally frustrated , whereas highly frustrated contacts 
 are energetically unfavorable in the native state . depictions of the 
 contacts on structural models typically show minimally frustrated 
 contacts in green and highly frustrated contacts in red . the average 
 of the frustration scores over all of the contacts made by a particular 
 residue are also plotted in a per - residue format . a webserver is now 
 available for performing these computations . for the work presented here , the minimally and highly frustrated 
 contacts are depicted on the lowest - energy structure from the boltzmann - reweighted 
 ensemble of structures from the amd simulations . to compute the average 
 per - residue frustration , we clustered the boltzmann reweighted ensemble 
 and averaged the frustration scores of all contacts made by each residue 
 in the representative structure from the three most populated clusters . it is interesting to note that the residual local frustration varied 
 between members of the ensemble , and the error bars on the residual 
 local frustration plots represent one standard deviation . note : thrombin has several numbering schemes ; the chymotrypsin 
 numbering scheme is used in the text because it is used in pdb files . this numbering 
 scheme creates problems for the data presentation in this paper ; therefore , 
 sequential numbering ( of the light chain or a - chain followed by the 
 heavy chain ) is used in the plots , and sequential residue numbers 
 are given in parentheses throughout the paper for reference . we performed a set of six 
 independent 20 ns conventional md simulations for both apo - thrombin 
 and ppack - thrombin . during the equilibration procedure , a rather large 
 conformational transition in the active site loops was observed for 
 apo - thrombin that involved a reorientation of both the -loop 
 ( 178195 ) and the na binding loop ( 264271 ) , 
 forming a more  open  active site pocket . a community 
 network analysis approach was applied to identify 
 groups of residues undergoing correlated motions in ppack - thrombin 
 and apo - thrombin . representative community network analyses obtained 
 from conventional md simulations are shown in figure 1 . the flow of information in the physical network of the protein 
 was measured by the edge betweenness , defined as the number of shortest 
 paths that pass through the edge in the network , and is a direct measure 
 of the strength of intercommunity communication within the network 
 ( black lines in figure 1a and b ) . ppack ligation 
 causes consolidation of the community structure including the two 
 active site communities most proximal to the ppack binding site ( green 
 and brown , figure 1 ) . the community that includes 
 part of the na binding loop in apo - thrombin ( figure 1c , brown ) consolidates with the active site serine195ct ( 241 ) , the 70s loop ( 98113 ) , residues 191194ct of the -loop ( 237240 ) , and the n - terminal 
 residues 1719ct ( 3840 ) in the ppack - liganded 
 form ( brown , figure 1d ) . ppack ligation also 
 acts to consolidate the n - terminal -barrel , which is formed 
 by two separate communities in apo - thrombin and forms a large community 
 that also contains the 30s loop ( 5562 ) and part of the 60s 
 loop ( 8294 ) ( orange , figure 1d ) . in 
 summary , the community analysis revealed consolidation of the na binding site , the base of the -loop , and the n - terminus 
 of the heavy chain into one community and most of the active site 
 loops into a second community upon ppack binding . these two communities , 
 which are strongly connected in ppack - thrombin , unite the residues 
 required for proteolytic catalysis . the a - chain community also becomes 
 more strongly connected to the active site and the community containing 
 the 70s loop . on the basis of the substantial consolidation observed 
 in the community analysis upon ppack binding , we set out to examine 
 whether there were concomitant changes in dynamics . community analysis of 
 apo - thrombin ( a , c ) and ppack - thrombin ( b , d ) . the two - dimensional view of the communities in panels ( a ) and ( b ) 
 depicts the relative size of the communities ( based on the number 
 of residues ) as colored circles of varying sizes with the thickness 
 of the connecting lines representing the relative interconnectivity 
 among communities . panels ( c ) and ( d ) are structural representations 
 of communities superimposed on pdb 1ppb . residual dipolar couplings ( rdcs ) , 
 which report on an ensemble average over all orientations of the magnetic 
 dipole interaction vector up to the chemical shift coalescence limit , 
 provide useful experimental data for determining the ensemble of structures 
 that best represents the dynamic properties of a protein . when the experimentally derived rdcs are compared to rdcs that are 
 back - calculated from ensembles of structures generated from amd simulations , 
 the acceleration level that provides the most realistic representative 
 structural ensemble can be identified . we previously demonstrated 
 that the rdcs measured on ppack - thrombin did not agree well with the 
 available crystal structures ( r = 0.72 ) . 80 ) when the rdcs were back - calculated from an ensemble of 
 structures obtained by conventional md . however , remarkable agreement was obtained between the experimental 
 rdcs and the rdcs back - calculated from the ensemble of structures 
 obtained from an amd simulation at the optimal acceleration level 
 ( r = 0.92 ) . the ensembles obtained from such amd simulations of both ppack - thrombin 
 and apo - thrombin are shown in figure 2 , with 
 the loops colored according to the scheme shown below the structure . ensemble 
 of the 10 representative structures of the rmsd clusters 
 from the boltzmann reweighted structures from the rdc - calibrated amd 
 for apo - thrombin ( left ) and ppack - thrombin ( right ) . the loops are 
 colored according to the schematic under the structures . for ppack - thrombin , order parameters ( s ) from conventional md simulations agreed extremely well 
 with those 
 measured by nmr relaxation experiments that are limited to the picosecond however , the fact that ensembles obtained from 
 amd were required for good agreement with the rdc measurements suggested 
 that motions on longer time scales are contributing to the solution 
 structure . therefore , order parameters for the n  h bond vectors 
 were calculated from the amd ensembles ( samd2 ) . a comparison 
 of samd2 for 
 apo - thrombin to those obtained previously 
 for ppack - thrombin is shown in figure 3a . most of the active site loops in both forms are 
 highly flexible , yet a marked decrease in flexibility is observed 
 upon active site ligation with ppack ( figure 3b ) . as expected , ppack ligation caused significant ordering ( s > 0.1 ) of residues that directly contact 
 the 
 ppack arg side chain . the loops that surround the active site also 
 experience significant ordering , including the 60s loop ( 8294 ) , 
 the -loop , and the 180s loop ( 225239 ) . some regions 
 distal to the ppack also showed significant ordering upon ppack ligation , 
 including the a - chain , the 30s loop , residue 221ct ( 269 ) 
 of the sodium binding loop , and residues in the c - terminal helix ( figure 3c ) . ( a ) order parameters calculated from amd simulations of 
 apo - thrombin 
 ( red ) and ppack - thrombin ( black ) . ( b ) differences in order parameters 
 between the forms : s = sppack2  sapo2 from amd order 
 parameters . ( c ) those residues with a s > + 0.1 are marked with red spheres on the structure of 
 thrombin 
 ( pdb code 1ppb ) , indicating stabilization upon 
 active site occupation by ppack . the schematic of important surface 
 loops is provided above the graph . to identify residues undergoing 
 correlated motions on longer time scales , amd simulations were performed 
 that were optimized based on previous work comparing experimental 
 rdcs to those back calculated from amd simulations carried out at 
 different acceleration levels . the analysis 
 of apo - thrombin revealed correlated motions between the active site 
 loops , exosite 1 , and other distal sites . in particular , the entire 
 -loop appears to undergo strongly correlated motions with the 
 a - chain residues 1h1dct and 1214cct ( 15 and 2025 ) , the catalytic triad , h57ct , d102ct , s195ct ( 79 , 135 , 241 ) , the 60s loop , 
 the 70s loop , the 90s loop ( 127133 ) , the surface strand under 
 the 70s loop ( 145151 ) , the 170s loop ( 204219 ) , and 
 the 180s loop ( figure 4 , lower triangle ) . whereas 
 these correlated motions appeared to involve the entire -loop 
 in apo - thrombin , only the tip of the -loop ( residues 146149ect ( 182190 ) ) appears to be undergoing the same set 
 of motions in ppack - thrombin ( figure 4 , upper 
 triangle ) . in apo - thrombin , the 170s and , to a lesser extent , the 
 180s loop , which are strongly correlated to the -loop , are 
 also correlated with the a - chain and catalytic residues . upon ppack , 
 the 30s loop and the 60s loop are weakly correlated , but the 30s loop 
 is not correlated to the 70s and 90s loops . correlated motions between 
 the 30s and 60s loops are stronger in ppack - thrombin , and these extend 
 to the 70s and 90s loops ( figure 4 ) . analysis of 
 correlated motions performed on the amd trajectories 
 of ppack - thrombin ( top triangle ) and apo - thrombin ( bottom triangle ) . the motions range from 0.0 ( no correlation , white ) to 1.0 ( completely 
 correlated , black ) . the schematic diagram indicating the location 
 of surface loops is inserted above the correlated motions plot . the 
 black boxes indicate correlations of the 170s and 180s loops with 
 the a - chain and surface loops that are stronger in apo - thrombin than 
 in ppack - thrombin . the blue boxes indicate correlations of the 30s 
 loop with the 60s , 70s , and 90s loops that are stronger in ppack thrombin 
 than in apo - thrombin . we applied a previously derived algorithm to identify the residual 
 local frustration in representative structures from rmsd clusters 
 of the boltzmann reweighted amd simulation results . according to the principle of minimal frustration , contacts made in the folded native state should 
 be minimally frustrated , meaning that they are energetically favorable . we previously showed that while most contacts made in the native state 
 are , indeed , minimally frustrated , some 1015% of contacts 
 are energetically unfavorable ( i.e. , highly frustrated ) in the native 
 state . these highly frustrated contacts map to functional sites and 
 are thought to have been preserved in evolution . both apo- and ppack - thrombin 
 show regions of high frustration in many of the surface loops ( figure 5 ) . to discover whether regions of high frustration 
 also map to dynamic regions , we compared the average residual frustration 
 across representative structures from the three most populated rmsd 
 clusters derived from the amd simulations to the order parameters . the order parameters derived from conventional md simulations ( sns2 ) agree very well with order parameters derived from nmr relaxation 
 experiments on thrombin and reveal the 
 disorder resulting from motions in the nanosecond time regime . the 
 order parameters derived from the rdc - calibrated amd simulations ( samd2 ) reveal the disorder resulting from motions on longer time scales . the sns2 did not correspond well to the regions of high residual frustration ; 
 however , the correspondence with the samd2 is remarkable 
 ( figure 6 ) . these results highlight that regions 
 of high residual frustration may have evolved to allow slow time scale 
 motions to occur with relative energetic ease , as previously suggested 
 by wolynes and colleagues . analysis of 
 residual local frustration in the lowest - energy 
 structure from the rdc - calibrated amd ensemble 
 of apo - thrombin ( left ) and ppack - thrombin ( right ) . the contacts that 
 are minimally frustrated are shown in green , and the contacts that 
 are highly frustrated are shown in red . thin lines represent water - mediated 
 contacts . the active site catalytic residues are shown in magenta . comparison of the order parameters reflecting 
 nanosecond time scale 
 motions versus longer time scale motions ( samd2 ) , with the average 
 per residue fraction of highly frustrated contacts for the three lowest - energy 
 structures from the amd simulation . ( a ) the sns ( gray ) and samd ( red ) for apo - thrombin are 
 compared to the average fraction of highly frustrated contacts ( cyan ) . ( b ) the sns2 ( gray ) and samd2 ( black ) for 
 ppack - thrombin are compared to the average fraction of highly frustrated 
 contacts ( blue ) . we 
 used a combination of community network analysis , rdc - calibrated 
 amd simulations , and analysis of residual frustration to explore the 
 dynamic ensemble of thrombin in solution from the nanosecond to the 
 microsecond time regime . comparative analysis of apo - thrombin and 
 the active - site - ligated thrombin systems identified differences in 
 the dynamic fluctuations and changes in correlated motions upon active 
 site ligation . analysis of thrombin complexed with the relatively 
 small substrate 
 analogue ppack showed a substantial rearrangement of the community 
 structure . ppack binding consolidated the two active site communities 
 most proximal to the ppack binding site as well as the na binding loop with the active site serine , the 70s loop , part of 
 the -loop , and part of the a - chain . interestingly , 
 this consolidation resulted in more residual local frustration near 
 the active site ( figure 5 ) . the generalized 
 cross - correlation analysis predicted that many 
 of the thrombin surface loops are undergoing correlated motions ( figure 4 ) . the analysis of residual local frustration revealed 
 that all of these loops are highly frustrated in the native structure 
 of thrombin ( figure 5 ) . the striking correlation 
 between the long timescale dynamics and the residual local frustration 
 suggests that evolution has selected for energetically unfavorable 
 contacts within the surface loops in order to facilitate larger amplitude 
 slower motions . upon ppack binding , many of the loops retain dynamics 
 and also remain highly frustrated . the surface loops exhibiting 
 correlated motions are not necessarily 
 structurally proximal , hinting at how allosteric regulation by binding 
 of ligands at sites distal to the active site might occur . upon ppack 
 binding , the active site community consolidates with the 70s loop 
 where fibrinogen binds to stimulate its own cleavage and where thrombomodulin 
 binds to alter substrate specificity . indeed , previous experiments 
 have demonstrated a thermodynamic coupling between the 70s loop and 
 the active site . another loop that is 
 distal to the active site but that differs between apo- and ppack - thrombin 
 is the -loop . residues leading up to the -loop are highly 
 frustrated and highly dynamic in apo - thrombin but become much less 
 dynamic in the ppack - bound form , with only the tip of the -loop 
 remaining dynamic . although the strongly correlated motions in apo - thrombin 
 between the -loop and the 170s loop almost completely disappear 
 upon ppack ligation , correlations are strengthened between the tip 
 of the loop and nearly every other loop in thrombin , including the 
 30s loop , the 60s loop , the 70s loop , the -strand connecting 
 the 60s and 70s loops , and the 90s loop . thus , we can speculate that 
 the reduced dynamics result from the binding of ppack , causing these 
 regions to sample a more restricted conformational space , and that 
 this reduced sampling leads to stronger correlated motions extending 
 between the active site and the surface loops involved in allosteric 
 regulation . the order parameter analysis revealed , as expected , that 
 active site occupation decreases the microsecond motions in the surface 
 loops ( figure 3 ) , consistent with amide h / d 
 exchange experiments that showed that binding of ppack dampened exchange 
 in all of the surface loops including those distal to the active site . millisecond 
 time scale could explain how active site occupation lowers the entropic 
 penalty for binding an exosite 1 allosteric modulator . thermodynamic measurements of binding show that 
 the allosteric pathway between these two sites operates in both directions and changes in order parameters in the thrombin 
 surface loops upon ppack binding to the active site are similar to 
 those seen upon tm56 binding to exosite 1 . the serine protease architecture is that of a double -barrel 
 fold in which the -strands are connected by surface loops . the active site is between the two -barrels , each of which 
 has several -strands connected by the surface loops . such a 
 -strand architecture provides a direct connection between distal 
 surface loops , allowing transmission of allosteric information by 
 subtle redistribution of the dynamic ensemble . such local unfolding or cracking was predicted some time 
 ago to be the fundamental mechanism of allostery , and residual local 
 frustration is indeed found at such sites . it will be interesting to see if dynamic allostery is present in 
 other proteins with mostly -strand structures similar to thrombin .