Frequency adaptive metadynamics for the calculation of rareevent kinetics
Abstract
The ability to predict accurate thermodynamic and kinetic properties in biomolecular systems is of both scientific and practical utility. While both remain very difficult, predictions of kinetics are particularly difficult because rates, in contrast to free energies, depend on the route taken and are thus not amenable to all enhanced sampling methods. It has recently been demonstrated that it is possible to recover kinetics through so called ‘infrequent metadynamics’ simulations, where the simulations are biased in a way that minimally corrupts the dynamics of moving between metastable states. This method, however, requires the bias to be added slowly, thus hampering applications to processes with only modest separations of timescales. Here we present a frequencyadaptive strategy which bridges normal and infrequent metadynamics. We show that this strategy can improve the precision and accuracy of rate calculations at fixed computational cost, and should be able to extend rate calculations for much slower kinetic processes.
I Introduction
Accessing long timescales in molecular dynamics (MD) simulations remains a longstanding challenge. Limited by the short time steps of a few femtoseconds, and despite recent progress in methods and hardware LindorffLarsen et al. (2011); Buch, Giorgino, and De Fabritiis (2011); Tiwary et al. (2015a); Paul et al. (2017), it remains difficult to access the timescales of milliseconds and beyond. This leaves many applications outside the realm of MD, including many structural rearrangements in biomolecules and binding and unbinding events of ligands and drug molecules. Many of these reactions are so called ‘rare events’ where the time scale of the reaction is orders of magnitude longer than the time it takes to cross the barrier between the states. In such cases, most of the simulation time ends up being used in simulating the local fluctuations inside free energy basins.
To study phenomena that involve basintobasin transitions that occur on long timescales, one typically has to rely on some combination of machine parallelism Dror et al. (2012); Lane et al. (2013), advanced strategies for analysing simulations Prinz et al. (2011); TrendelkampSchroer and Noé (2016); Chong and Zuckerman (2016) and enhanced sampling methods Bernardi, Melo, and Schulten (2015); Valsson, Tiwary, and Parrinello (2016); De Vivo et al. (2016); Bruce et al. (2018) that allow for efficient exploration of phase space. Metadynamics is one such popular and now commonly used enhanced sampling method that involves the periodic application of a historydependent biasing potential to a few selected degrees of freedom, typically also called collective variables (CVs) Laio and Parrinello (2002). Through this bias, the system is discouraged from getting trapped in low energy basins, and one can observe processes that would be far beyond the timescales accessible in normal MD, while still maintaining complete atomic resolution. Originally designed to explore and reconstruct the equilibrium free energy surface Bussi and Branduardi (2015), it has recently been shown that metadynamics can also be used to calculate kinetic properties Tiwary and Parrinello (2013); McCarty et al. (2015); Tiwary et al. (2015b, a); Mondal, Tiwary, and Berne (2016); Fleming, Tiwary, and Pfaendtner (2016); Tung and Pfaendtner (2016); Sun et al. (2017); Fu, L. Oliveira, and Pfaendtner (2017); Wang, Martins, and LindorffLarsen (2017); Tiwary (2017); Casasnovas et al. (2017); Tiwary, Mondal, and Berne (2017). In particular, inspired by the pioneering work of Grubmüller Grubmüller (1995) and Voter Voter (1997), it has recently been shown that the unbiased kinetics can be correctly recovered from metadynamics simulations using a method called infrequent metadynamics (InMetaD)Tiwary and Parrinello (2013).
The basic idea in InMetaD is to add a bias to the free energy landscape sufficiently infrequently so that only the free energy basins, but not the free energy barriers, experience the biasing potential, , where is the simulation time and is one or more CV’s chosen to distinguish the different free energy basins. By adding bias infrequently enough compared to the time spent in barrier regions, the landscape can be filled up to speed up the transitions without perturbing the sequence of statetostate transitions. The effect of the bias on time can then be reweighted through an acceleration factor:
(1) 
Here is the inverse temperature and the angular brackets denote an average over a metadynamics run up until the simulation time .
In the application of InMetaD to recover the correct rates, there are two key assumptions. First, the statetostate transitions are of the rare event type: namely, the system is trapped in a basin for a duration long enough that memory of the previous history is lost, and when the system does translocate into another basin, it does so rather rapidly. Second, it is important that no substantial bias is added to the transition state (TS) region during the simulation. This requirement can be achieved by adjusting the bias deposition frequency, determined by the time between two instances of bias deposition. If the deposition frequency is kept low enough, it becomes possible to keep the TS region unbiased. This second requirement, however, also means that it takes longer to fill up the basin, revealing one of the practical limitations of applying InMetaD. It is this second requirement that we address and improve in this work.
In particular, we asked ourselves the question, can we design a bias deposition scheme so that the frequency is high near the bottoms of the free energy basins, but decreases gradually so as to lower the risk of biasing the TS regions? Our scheme, which we term FrequencyAdaptive Metadynamics (FaMetaD), is illustrated using a ligand (benzene) unbinding from the T4 lysozyme (T4L) L99A mutant as an example (Fig. 1). In particular, we designed a strategy that uses a high frequency at the beginning of the simulation (to fill up the basin quickly) and then slows progressively down (to minimize the risk of perturbing the TS). In this way, we aim to improve the reliability, accuracy and robustness of the calculations without additional computational cost.
Ii Methods
Our adaptive frequency scheme reads:
(2) 
where is the initial deposition time between adding a bias, similar to the relatively short time in normal metadynamics, and is a cutoff value for the frequency equal to or larger than the deposition time originally proposed in InMetaD. These two parameters bridge the normal metadynamics and InMetaD, by controlling the minimal and maximal bias deposition frequency. is the instantaneous acceleration factor at simulation time in a metadynamics simulation.
Through the starting deposition time , we can modulate the enhancement in our ability to fill free energy wells, relative to an InMetaD run performed with a constant stride of . For practical considerations the choice of is dependent on the computational resources. The key free parameter here is , the threshold value for the acceleration factor to trigger the gradual change from normal (frequent) metadynamics to InMetaD. The choice of requires some considerations. On one hand, should be as large as possible to delay the switch from normal metadynamics to InMetaD so as to get significant enhancement in basin filling. On the other hand, too large could lead to problems that a transition might occur when the frequency is still less than the typical duration of a reactive trajectory crossing over from one basin to another.
We now estimate a value for in terms of the expected transition time available from either experimental measurements or an estimate, the simulation time, (determined also by computational resources), and a ‘safety coefficient’ . We seek the following to hold true:
(3) 
If we use a proteindrug system as an example, and we (i) can estimate the residence time to be roughly on the order of a second (), (ii) aim to observe the transition within a 100 ns metadynamics run (), and (iii) use to counteract the risk that a transition time falls in the long tail of Poisson distribution, we thus obtain . As we demonstrate by an example below, if one chooses too high a value of this leads to perturbed (and erroneous) kinetics; an error that may be detected by testing whether the observed distribution of transition times follows a Poisson distribution Salvalaglio, Tiwary, and Parrinello (2014).
As per Eq. (2), the frequency is also changed as a function of . In practice, we observed significant fluctuation of , that might cause problems if is small at the transition time point. Therefore we finally modified Eq. (2) to become a monotonously increasing function:
(4) 
where is the instantaneous deposition time at step N in a metadynamics simulation with a MD time step . The frequencyadaptive scheme was implemented in a development version of the PLUMED2.2 codeTribello et al. (2014).
Parameters: (ps), (ps), h (kJ/mol)  ()  Pvalue  Cost (s)  Set  

Unbiased MD  T=300K  112^{a}  300  
InMetaD  ,  165  0.40.3  2.8  
,  123  0.30.2  1.5  A  
,  196  0.30.2  0.8  B  
,  165  0.10.2  0.5  C  
,  3215  0.10.1  0.3  D  
,  9238  0.020.05  0.2  E  
,  9861  0.010.01  0.1  F  
FaMetaD  ,,,  153  0.40.3  1.5  A 
,,,  143  0.50.3  0.6  B  
,,,  196  0.20.2  0.4  C  
,,,  247  0.10.1  0.3  D  
,,,  2910  0.10.2  0.2  E  
,,,  2411  0.020.03  0.1  F 

From ref. Wang, Martins, and LindorffLarsen (2017).
Iii Results and Discussion
iii.1 Test on a fourstate model system
To benchmark our results, we first consider a fiveresidue peptide (NmeAla3Ace) as a model system with a nontrivial free energy landscape involving multiple conformational states Wang, Martins, and LindorffLarsen (2017). We consider the slowest statetostate transition time, , which has been estimated to be from unbiased MD Wang, Martins, and LindorffLarsen (2017), as the benchmark target. We used the same computational setup and CVs as previously described Wang, Martins, and LindorffLarsen (2017).
As a baseline, we used InMetaD with a fixed height of Gaussian bias potential ( kJ/mol) but with different times of bias deposition ranging from 2 ps to 200 ps (40 runs for each parameter set) (Table 1 and Fig. S1). The reliability of the calculated transition times were verified a posteriori using a KolmogorovSmirnov test Salvalaglio, Tiwary, and Parrinello (2014) to examine whether their cumulative distribution function is Poissonian. If the value is low (e.g. less than 0.05) then it is likely that the distribution has been perturbed because of too aggressive application of the biasing potential. As expected, we find that simulations that used a low frequency ( ps) gave rise to a consistent estimation of , with values close to those observed in unbiased MD, while the simulations with high frequency ( ps) resulted in substantially longer and less reliable times (Table 1 and Fig. S1). This correlation between the value and the deposition time in the InMetaD simulations can be explained by the fact that the higher bias frequency results in higher risk of biasing the TS regions and more nonPoissonian distribution. We also note that the transition times appear to be overestimated in these cases.
To compare the FaMetaD simulations with the InMetaD baseline, we designed six sets (A to F) of FaMetaD simulations to have comparable computational costs (Table 1 and Fig. 2). Overall, the results of FaMetaD show the same trends as InMetaD and support the same conclusions. In other words, both InMetaD and FaMetaD reveal the tradeoff between reliability, accuracy and computational cost. There are, however, some important differences that highlight the improvement obtained through FaMetaD. First, in each case there is a small but notable improvement in the reliability of the calculations, as judged by the greater values that suggest that FaMetaD perturbs the TS less than InMetaD. Most important, however, is the accuracy of the calculations. While both InMetaD and FaMetaD achieve accurate results when using the most conservative parameters (sets A–C), there are dramatic differences when more aggressive parameters are used (D–F). Importantly, we observe a much more ‘gracious’ decline in accuracy with more aggressive parameters in FaMetaD comparing to InMetaD. Thus together these results suggest that the frequencyadaptive scheme can further improve not only the reliability but also the accuracy of the calculation, without the need of increasing computational burden, and also that the reliability is more robust to the choice of the simulation parameters.
iii.2 Application on proteinligand binding
Methods  Parameters: (ps), (ps), (kJ/mol)  Time (ms)  pvalue  Cost 

Set 1: Benzene Binding ()  
InMetaD^{b}  ,  95  0.10.1  4.4s 
FaMetaD  ,,,  147  0.20.1  3.0s 
Set 2: Benzene Unbinding ()  
InMetaD^{b}  ,  16859  0.40.3  6.7s 
FaMetaD  ,,,  17668  0.30.2  5.5s 
Set 3: Indole Unbinding ()  
InMetaD  ,  10287  0.20.2  4.5s 
FaMetaD  ,,,  16895  0.20.1  2.1s 
26s 

In all simulations, the ligand concentration is mM.

The results are from our recent work Wang, Martins, and LindorffLarsen (2017).
We consider now the more complex case of a ligand binding to or escaping from a buried internal cavity in the L99A mutant of T4 lysozyme, processes which occur on timescales of milliseconds or more Feher, Baldwin, and Dahlquist (1996); Bouvignies et al. (2011); Wang, Papaleo, and LindorffLarsen (2016). To benchmark the results, we performed three sets of metadynamics simulations, up to 26s in total (including 12 s InMetaD simulations of T4L L99A with benzene from our recent work Wang, Martins, and LindorffLarsen (2017)). In each set, we compared the results of InMetaD with that of FaMetaD. In set 1 and 2 we calculated the time constants for binding () and unbinding () of benzene, while we in set 3 calculated the time for indole to escape the pocked (). We performed 20 independent runs for each set of simulations (using the CHARMM22* force field Piana, LindorffLarsen, and Shaw (2011) for protein and CGenFF force field Vanommeslaeghe et al. (2010) for the ligands) to collect the transition times, from which we obtained and (Table 2).
Again, we find consistency between the results of InMetaD and FaMetaD, with e.g. and to be ms and ms, respectively. As previously described Wang, Martins, and LindorffLarsen (2017), we can use these values to estimate the binding free energy to be kJ/mol, a value that agrees well with calorimetric Morton, Baase, and Matthews (1995) and NMR Feher, Baldwin, and Dahlquist (1996) measurements.
In set 1, the InMetaD simulation was performed with ps and kJ/mol. The frequencyadaptive scheme allows us to perform FaMetaD with weaker bias ( kJ/mol) and ending with longer and more conservative deposition time ( ps). This parameter set used less simulation time but resulted in a more reliable estimation. In set 2, the FaMetaD simulation was performed using the same ps and kJ/mol as that used in the InMetaD simulation, but with . Given that ms, ns and , according to Eq. (3), this allows us to judge that is a fairly conservative parameter. The estimated values of in this set are indistinguishable between InMetaD and FaMetaD within the error bars but at slightly less computational cost. In set 3, we used similar parameters as in set 2 except a bit more aggressive . Again, this parameter set resulted in of ms from FaMetaD that is reasonably close to the estimation from InMetaD. Remarkably, FaMetaD allowed us to reduce more than half the computational cost, without loss of accuracy. Overall, the application in the case of L99A binding with two ligands suggests that the frequencyadaptive scheme can improve the reliabilityaccuracyefficency balance of metadynamics on kinetics calculation.
Iv Conclusions
Many biological processes occur far from equilibrium, and kinetic properties can play an important role in biology and biochemistry. For example, there has been continued interest in determining ligand residence times in the context of drug optimization Copeland, Pompliano, and Meek (2006), and millisecond conformational dynamics in an enzyme has been shown to underlie an intriguing phenomenon of kinetic cooperativity of relevance to disease Larion et al. (2015).
The principle of microscopic reversibility, however, sets limits to how the kinetic properties can be varied independently of thermodynamics. Thus, for example, increasing the residence time of a ligand will also increase its thermodynamic affinity, unless there is also a simultaneous drop in the rate of binding. Thus, in practice it may be in many cases be difficult to disentangle the effects of kinetics and thermodynamics.
Taken together, the above considerations suggest the need for improved methods for understanding and ultimately predicting the rate constants of biological processes. Further, the ability to calculate the kinetics of conformational exchange or ligand binding and unbinding provides an alternative approach to determine equilibrium properties Wang, Martins, and LindorffLarsen (2017). For these and other reasons, several methods have been developed to enable the estimation of kinetics from simulations Bruce et al. (2018).
We have here proposed a modification to the already very powerful InMetaD algorithm, which leads to a further improvement in the reliability and accuracy in recovering unbiased transition rates from biased metadynamics simulations. The basic idea in the resulting FaMetaD approach is that, by filling up the basin more rapidly in the beginning and only using an infrequent bias near the barrier, we may spend the computational time where it is needed. We anticipate that our scheme will prove particularly useful in two different aspects. First, by enabling a decreased bias in the transition state region, we obtain more accurate kinetics at fixed computational cost, leading also to the observed robustness of our approach. Second, in the case of large barriers, it would be prohibitively slow to use a very infrequent bias through the entire duration of the simulation. Thus, the FaMetaD provides a practical approach to study rare events that involve escaping deep free energy minima, and which occur on timescales well beyond those accessible to current simulation methods. Here we have opted to examine processes that can be studied using both InMetaD and FaMetaD, but in the future we aim to apply this approach to barrier crossing events that occur on even longer timescales. With more examples in hand, we also expect that in the future it may be possible to design improved biasing schemes to increase the computational efficiency further, and at the same time retain the robustness of the FaMetaD approach.
V Acknowledgement
The authors thank Tristan Bereau and Claudio Perego for a critical reading of the manuscript. K.L.L. acknowledges funding by a HallasMøller Stipend from the Novo Nordisk Foundation and the BRAINSTRUC initiative from the Lundbeck Foundation. M.P. acknowledges funding from the National Centre for Computational Design and Discovery of Novel Materials MARVEL and European Union grant ERC2014AdG670227/VARMET.
References
 LindorffLarsen et al. (2011) K. LindorffLarsen, S. Piana, R. O. Dror, and D. E. Shaw, “How fastfolding proteins fold,” Science 334, 517–520 (2011).
 Buch, Giorgino, and De Fabritiis (2011) I. Buch, T. Giorgino, and G. De Fabritiis, ‘‘Complete reconstruction of an enzymeinhibitor binding process by molecular dynamics simulations,” Proc Natl Acad Sci 108, 10184–10189 (2011).
 Tiwary et al. (2015a) P. Tiwary, V. Limongelli, M. Salvalaglio, and M. Parrinello, “Kinetics of protein–ligand unbinding: Predicting pathways, rates, and ratelimiting steps,” Proc Natl Acad Sci 112, E386–E391 (2015a).
 Paul et al. (2017) F. Paul, C. Wehmeyer, E. T. Abualrous, H. Wu, M. D. Crabtree, J. Schöneberg, J. Clarke, C. Freund, T. R. Weikl, and F. Noé, “Proteinpeptide association kinetics beyond the seconds timescale from atomistic simulations,” Nature Communications 8, 1095 (2017).
 Dror et al. (2012) R. O. Dror, R. M. Dirks, J. Grossman, H. Xu, and D. E. Shaw, “Biomolecular simulation: a computational microscope for molecular biology,” Annu Rev Biophys 41, 429–452 (2012).
 Lane et al. (2013) T. J. Lane, D. Shukla, K. A. Beauchamp, and V. S. Pande, “To milliseconds and beyond: challenges in the simulation of protein folding,” Curr Opin Struct Biol 23, 58–65 (2013).
 Prinz et al. (2011) J.H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, C. Schütte, and F. Noé, “Markov models of molecular kinetics: Generation and validation,” J Chem Phys 134, 174105 (2011).
 TrendelkampSchroer and Noé (2016) B. TrendelkampSchroer and F. Noé, “Efficient estimation of rareevent kinetics,” Physical Review X 6, 011009 (2016).
 Chong and Zuckerman (2016) L. T. Chong and D. M. Zuckerman, “Weighted ensemble simulation methods and applications,” Annu Rev Biophys 46 (2016).
 Bernardi, Melo, and Schulten (2015) R. C. Bernardi, M. C. Melo, and K. Schulten, “Enhanced sampling techniques in molecular dynamics simulations of biological systems,” Biochim Biophys Acta, Gen Subj 1850, 872–877 (2015).
 Valsson, Tiwary, and Parrinello (2016) O. Valsson, P. Tiwary, and M. Parrinello, “Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint,” Annu Rev Phys Chem 67, 159–184 (2016).
 De Vivo et al. (2016) M. De Vivo, M. Masetti, G. Bottegoni, and A. Cavalli, “Role of molecular dynamics and related methods in drug discovery,” J Med Chem (2016).
 Bruce et al. (2018) N. J. Bruce, G. K. Ganotra, D. B. Kokh, S. K. Sadiq, and R. C. Wade, “New approaches for computing ligand–receptor binding kinetics,” Curr Opin Struct Biol 49, 1–10 (2018).
 Laio and Parrinello (2002) A. Laio and M. Parrinello, “Escaping freeenergy minima,” Proc Natl Acad Sci 99, 12562–12566 (2002).
 Bussi and Branduardi (2015) G. Bussi and D. Branduardi, “Freeenergy calculations with metadynamics: theory and practice,” Reviews in Computational Chemistry Volume 28 , 1–49 (2015).
 Tiwary and Parrinello (2013) P. Tiwary and M. Parrinello, “From metadynamics to dynamics,” Phys Rev Lett 111, 230602 (2013).
 McCarty et al. (2015) J. McCarty, O. Valsson, P. Tiwary, and M. Parrinello, “Variationally optimized freeenergy flooding for rate calculation,” Phys Rev Lett 115, 070601 (2015).
 Tiwary et al. (2015b) P. Tiwary, J. Mondal, J. A. Morrone, and B. Berne, “Role of water and steric constraints in the kinetics of cavity–ligand unbinding,” Proc Natl Acad Sci 112, 12015–12019 (2015b).
 Mondal, Tiwary, and Berne (2016) J. Mondal, P. Tiwary, and B. Berne, “How a kinase inhibitor withstands gatekeeper residue mutations,” J Am Chem Soc 138, 4608–4615 (2016).
 Fleming, Tiwary, and Pfaendtner (2016) K. L. Fleming, P. Tiwary, and J. Pfaendtner, “New approach for investigating reaction dynamics and rates with ab initio calculations,” J Phys Chem A 120, 299–305 (2016).
 Tung and Pfaendtner (2016) H.J. Tung and J. Pfaendtner, “Kinetics and mechanism of ionicliquid induced protein unfolding: application to the model protein hp35,” Molecular Systems Design & Engineering 1, 382–390 (2016).
 Sun et al. (2017) H. Sun, Y. Li, M. Shen, D. Li, Y. Kang, and T. Hou, “Characterizing drug–target residence time with metadynamics: how to achieve dissociation rate efficiently without losing accuracy against timeconsuming approaches,” J Chem Inf Model 57, 1895–1906 (2017).
 Fu, L. Oliveira, and Pfaendtner (2017) C. D. Fu, L. F. L. Oliveira, and J. Pfaendtner, “Determining energy barriers and selectivities of a multipathway system with infrequent metadynamics,” J Chem Phys 146, 014108 (2017).
 Wang, Martins, and LindorffLarsen (2017) Y. Wang, J. M. Martins, and K. LindorffLarsen, “Biomolecular conformational changes and ligand binding: from kinetics to thermodynamics,” Chemical Science 8, 6466–6473 (2017).
 Tiwary (2017) P. Tiwary, “Molecular determinants and bottlenecks in the dissociation dynamics of biotinstreptavidin,” J Phys Chem B (2017).
 Casasnovas et al. (2017) R. Casasnovas, V. Limongelli, P. Tiwary, P. Carloni, and M. Parrinello, “Unbinding kinetics of a p38 map kinase type ii inhibitor from metadynamics simulations,” J Am Chem Soc 139, 4780–4788 (2017).
 Tiwary, Mondal, and Berne (2017) P. Tiwary, J. Mondal, and B. Berne, “How and when does an anticancer drug leave its binding site?” Sci Adv 3, e1700014 (2017).
 Grubmüller (1995) H. Grubmüller, “Predicting slow structural transitions in macromolecular systems: Conformational flooding,” Phys Rev E 52, 2893 (1995).
 Voter (1997) A. F. Voter, “Hyperdynamics: Accelerated molecular dynamics of infrequent events,” Phys Rev Lett 78, 3908 (1997).
 Salvalaglio, Tiwary, and Parrinello (2014) M. Salvalaglio, P. Tiwary, and M. Parrinello, “Assessing the reliability of the dynamics reconstructed from metadynamics,” J Chem Theo Comput 10, 1420–1425 (2014).
 Tribello et al. (2014) G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi, “Plumed 2: New feathers for an old bird.” Comput Phys Commun 185, 604–613 (2014).
 Feher, Baldwin, and Dahlquist (1996) V. A. Feher, E. P. Baldwin, and F. W. Dahlquist, “Access of ligands to cavities within the core of a protein is rapid,” Nat Struct Biol 3 (1996).
 Bouvignies et al. (2011) G. Bouvignies, P. Vallurupalli, D. F. Hansen, B. E. Correia, O. Lange, A. Bah, R. M. Vernon, F. W. Dahlquist, D. Baker, and L. E. Kay, “Solution structure of a minor and transiently formed state of a t4 lysozyme mutant,” Nature 477, 111–114 (2011).
 Wang, Papaleo, and LindorffLarsen (2016) Y. Wang, E. Papaleo, and K. LindorffLarsen, “Mapping transiently formed and sparsely populated conformations on a complex energy landscape,” eLife 5, e17505 (2016).
 Piana, LindorffLarsen, and Shaw (2011) S. Piana, K. LindorffLarsen, and D. E. Shaw, “How robust are protein folding simulations with respect to force field parameterization?” Biophys J 100, L47–L49 (2011).
 Vanommeslaeghe et al. (2010) K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, and I. Vorobyov, “Charmm general force field: A force field for druglike molecules compatible with the charmm allatom additive biological force fields,” J Comput Chem 31, 671–690 (2010).
 Morton, Baase, and Matthews (1995) A. Morton, W. A. Baase, and B. W. Matthews, “Energetic origins of specificity of ligand binding in an interior nonpolar cavity of t4 lysozyme,” Biochemistry 34, 8564–8575 (1995).
 Copeland, Pompliano, and Meek (2006) R. A. Copeland, D. L. Pompliano, and T. D. Meek, “Drug–target residence time and its implications for lead optimization,” Nature reviews Drug discovery 5, 730 (2006).
 Larion et al. (2015) M. Larion, A. L. Hansen, F. Zhang, L. BruschweilerLi, V. Tugarinov, B. G. Miller, and R. Brüschweiler, “Kinetic cooperativity in human pancreatic glucokinase originates from millisecond dynamics of the small domain,” Angewandte Chemie 127, 8247–8250 (2015).