SEEKR: Simulation Enabled Estimation of Kinetic Rates, A Computational Tool to Estimate Molecular Kinetics and its Application to Trypsin-Benzamidine Binding
We present the Simulation Enabled Estimation of Kinetic Rates (SEEKR) package, a suite of open-source scripts and tools designed to enable researchers to perform multi-scale computation of the kinetics of molecular binding, unbinding, and transport using a combination of molecular dynamics, Brownian dynamics, and milestoning theory. To demonstrate its utility, we compute the kon, koff, and \xce\x94Gbind for the protein trypsin with its noncovalent binder, benzamidine, and examine the kinetics and other results generated in the context of the new software, and compare our findings to previous studies performed on the same system. We compute a kon estimate of 2.1\xc2\xb10.3\xe2\x80\xa2107 M\xe2\x88\x921s\xe2\x88\x921, a koff estimate of 83\xc2\xb114 s\xe2\x88\x921, and a \xce\x94Gbind of \xe2\x88\x927.4\xc2\xb10.2 kcal\xe2\x80\xa2mol\xe2\x88\x921, all of which compare closely to the experimentally measured values of 2.9\xe2\x80\xa2107 M\xe2\x88\x921s\xe2\x88\x921, 600\xc2\xb1300 s\xe2\x88\x921, and \xe2\x88\x926.7 kcal\xe2\x80\xa2mol\xe2\x88\x921, respectively.
Introduction<!>Theory<!>Description of the SEEKR package<!><!>Description of the SEEKR package<!>Trypsin structure preparation and SEEKR creation of milestoning structures<!>MD simulations<!>BD simulations<!>Milestoning calculations<!>Results<!>Discussion<!>Conclusions
<p>Elucidating the kinetics and thermodynamics of binding and unbinding processes between a biomolecule and a substrate remains an important challenge in the field of molecular biophysics. Countless processes within the cell involve the association of a biomolecule with a metabolite, signaling molecule, toxin, drug, or another biomolecule.1–2 Many of these interactions have important kinetic considerations: for instance, the speed of reactions or the residence time of an intermolecular encounter.2–3</p><p>Significant effort has been expended to accurately estimate the thermodynamics of binding using a variety of methods, particularly in the field of drug discovery, where the identification of a tight binder is an integral step towards obtaining a potential drug molecule that would accomplish a desired medical result.4–8 While the thermodynamics of binding, encapsulated in the quantity of the free energy ΔGbind of receptor-ligand complex formation, is an important factor in the binding process, a comprehensive understanding of the binding process requires consideration of binding kinetics and reaction rates.</p><p>Many theoretical approaches and simulation methods have been used to estimate both the thermodynamics and kinetics of binding. For instance, specialized machinery and long molecular dynamics (MD) simulations can be used in a 'brute force' approach, although it is relatively costly compared to other methods.9–12 Markov models13–19 can also be used to investigate the kinetics of binding,20–22 as can milestoning.7 Additional clever methodologies can be used to speed the computation using MD.6, 8, 23–25 Brownian dynamics (BD) can also be used to approach the problem of binding kinetics,5, 26–29 as can Smoluchowski equation solvers.30</p><p>Our past work31–33 has focused on using a multi-scale combination of MD and BD, unified through the theoretical framework of milestoning. In our previous study, we presented a hybrid MD/BD/milestoning methodology to conduct our investigations into the kinetics of binding between superoxide dismutase and its natural substrate, the superoxide anion, and between troponin C and its natural substrate, the calcium ion.31 Here, we make available a software package, SEEKR, that implements this method with significant improvements in automation, usability, and analysis. We demonstrate the utility of SEEKR by applying it to estimate the kon, the koff, and ΔGbind between the serine protease trypsin and its ligand, benzamidine. In addition to the SEEKR software to perform milestoning calculations on any receptor-ligand system, we also make available a user guide, tutorial, and workflow to allow users to repeat our simulations and analysis for the trypsin-benzamidine system, and compute kinetics and thermodynamics for additional receptor-ligand systems.</p><!><p>The rationale and methodology behind our usage of milestoning to estimate kinetics using both MD and BD has been described recently in detail31 for multiple applications. Our implementation to the trypsin-benzamidine receptor-ligand system in this study was adapted with few changes, the majority consisting of improvements in software efficiency.</p><p>In the case of bimolecular association, the kinetics of binding and unbinding can be represented respectively by two quantities, kon and koff, which are frequently depicted according to the following equation: 
Eq. 1A+B⇄koffkonAB</p><p>Which is shorthand for specifying that the values kon and koff function as parameters within the following differential equations: 
Eq. 2d[AB]dt=kon[A][B]-koff[AB]
Eq. 3d[A]dt=koff[AB]-kon[A][B]
Eq. 4d[B]dt=koff[AB]-kon[A][B]</p><p>Where [A], [B], and [AB] represent the concentrations of chemical species A, B, and their complex AB. The kon and koff relate to the dissociation constant KD, and by extension, a free energy of association ΔGbind:6</p><p>Where R is the gas constant, T is temperature, and K⊝ is a factor equal to one, in units of concentration.</p><p>The theory of milestoning has been formulated to compute kinetic and thermodynamic details of a process if the states of that process are represented as carefully chosen surfaces in phase space. These surfaces are known as "milestones".34 In this study, the milestones are represented as concentric spherical shells (figure 1) that encapsulate the binding site of the receptor. These spherical milestones are used for the computation of kon, koff, and ΔGbind. Milestoning theory allows us to approach the problem of kinetics by utilizing a multi-scale strategy. We use highly-detailed, but computationally expensive MD simulations to observe transitions between milestones closer to the binding site so that molecular flexibility will be a component of the transitions between milestones. We then use BD for the larger and more widely-spaced milestones far from the binding site, where fast sampling of long trajectories is required and rigid body dynamics and implicit solvent are adequate21, 28, 35 to model transition times and probabilities. In this way, we take advantage of fully flexible MD where molecular flexibility is required, and also take advantage of the computation efficiency of BD where molecular flexibility is less important. Milestoning is the theory that combines the MD and BD components, by allowing statistics to be obtained in each regime independently, and then unifying the statistics through a rigorous theory that is agnostic to the method that was used to obtain them. Since the statistics of each milestone are obtained independently from the others, and since milestoning theory is a robust framework that can utilize information obtained by either Brownian or Newtonian dynamics,36 we can choose whichever simulation method is most appropriate and convenient for that milestone.</p><p>By sampling transition statistics and times between the milestones using numerous short simulations, one can construct a transition kernel K that represents the transition probabilities and an incubation time vector 〈t〉 that represents the average times of a system traversing the milestones.37–38</p><p>The transition kernel K is a square matrix whose elements are constructed according to the following formula: 
Eq. 6Kij=ni→j∑kni→k</p><p>Where ni
→
j is the number of trajectories that begin at a given milestone i and end at an adjacent milestone j. And the incubation time vector 〈t〉 has elements that are constructed according to the following formula: 
Eq. 7〈t〉i=∑ltl∑kni→k</p><p>Where tl is the time of the l'th successful forward trajectory starting at milestone i, and ni
→
k, as before, is the number of trajectories beginning at milestone i and ending at milestone k. Therefore, 〈t〉i represents the average time spent by the system after crossing i and before crossing any other milestone.</p><p>In order to compute a free energy profile along the milestones, we must first obtain the stationary flux vector qstat along the milestones by computing the principle eigenvector of K.</p><p>Then qstat must be multiplied elementwise by 〈t〉 to find the stationary probability vector pstat.</p><p>Finally, pstat,i relates to the relative free energy ΔGi at milestone i according to the following: 
Eq. 10ΔGi=-RTln(pstat,i/pstat,ref)</p><p>Where the index of pstat,ref is any reference state, such as the lowest energy, bound state. The value of pstat,ref is found by applying Eq. 9 to the chosen reference state.</p><p>To compute the kon, we utilize the formula that is also used in BD theory26: 
Eq. 11kon=k(b)β</p><p>Where k(b) is computed using the following formula: 
Eq. 12k(b)=[∫b∞eW(r)kBT4πr2D(r)dr]-1</p><p>The value k(b) represents the rate constant at which the ligand particles are crossing the b-surface, W(r) and D(r) are the potential of mean force and diffusion coefficient, respectively, that the ligand experiences at a distance r from the center of the receptor beyond the b-surface.26
D(r) is computed by generating a Rotne-Prager diffusion tensor to approximate the hydrodynamics of a two body interaction in a viscous medium.39 The value k(b) is computed automatically in BrownDye.</p><p>To find β, which represents the proportion of ligands crossing the b-surface that continue on to bind to the receptor, a starting probability vector q0 must be obtained in BD simulations by running a large number of conventional BD simulations where ligand molecules are started on a b-surface surrounding a receptor molecule. As the simulations run, and the proportion of trajectories that touch the outermost milestone(s) that encompasses a binding site on the biomolecule, rather than escaping to an infinite distance, are counted. In this case, q0 becomes: 
Eq. 13q0=[0,…,0,q0,i,0,…,0,q0,j,0,…,0,q0,∞,0,…,0]T</p><p>Where i and j are the indices of one or more of these outermost site-encompassing milestones, q0,i, q0,j, are the probabilities that a BD trajectory started on the b-surface descend and touch these milestones, and q0,∞ is the probability that a trajectory diffuses away to an infinite distance. All the entries in q0 must be normalized such that their sum equals a value of one. An "infinity" state in both vector q0 and in matrix K, represents the condition in which the ligand has escaped to an infinite distance from the receptor.</p><p>Next the transition matrix K must be modified to a new matrix K̂ such that the milestones representing the bound and "infinity" states are sink states. That is, they all must have a probability of one that they transition only to themselves, and a zero probability to transition to anything else.</p><p>Once K̂ and q0 are properly defined, we compute the static flux vector21
q∞.</p><p>Finally, we obtain β: 
Eq. 16β=∑iq∞,i</p><p>Where i is the index of one of the bound states.</p><p>To compute the koff, we must return to the initial definition of matrix K as specified in Eq. 6. But it must be modified by introducing a "draining" state i by changing K into a draining matrix K̃ according to the following: 
Eq. 17K∼ij=0,∀j</p><p>That is, once we have decided that i is the draining state, we set that entire column of the matrix K̃ to zeros, while all other columns are kept the same as they were in K. In the SEEKR implementation, the outermost non-infinite milestone is considered to be the draining state. Then, we compute a mean first passage time (MFPT) τ: 
Eq. 18τ=p0(I-K∼T)-1〈t〉</p><p>Where p0 is a starting distribution of probabilities along each milestone, and K̃T is the transpose of matrix K̃. We set p0,i to be 1 if i was a bound state, and set p0,i to be equal to 0 otherwise. The MFTP τ is equivalent to a residence time of the ligand within the binding site, and can be related to the koff according to the following relation: 
Eq. 19koff=1τ</p><!><p>SEEKR is a collection of scripts and files designed to automate the preparation and analysis of ligand-receptor kinetic calculations that use a multi-scale MD/BD/milestoning framework.</p><p>SEEKR does not run the simulations themselves, but instead relies on the well-established NAMD40 and BrownDye41 programs. In this case, SEEKR is more of a specialist interface or tool that automates the cumbersome process of preparing, running, and analyzing a particular type of multi-scale milestoning calculation so that researchers will be able to run them more easily than if the process were done manually.</p><p>SEEKR programs are classified into three general categories:</p><!><p>Preparation: These scripts and modules accept input from the user in order to construct all the necessary files needed by both NAMD and BrownDye to run their respective simulations. The files are organized into a file tree whose branches represent the various independent milestones, which simulation method is being used (MD or BD), and the various stages of the calculations. When run, the user will have all the required files arranged and poised for simulation and milestoning calculations.</p><p>Running: Other scripts aid the user in running the MD and BD simulations locally and on supercomputers. For instance, SEEKR contains a script to prepare the submission of the computationally-intensive MD simulation jobs to a SLURM supercomputer queue, and when the allotted time runs out, the script prepares all the necessary resubmission files for one, some, or all of the milestones with a single command. Other scripts use previous BD trajectory output to prepare and run ensembles of BD simulations from first hitting point distributions (FHPD).</p><p>Analysis: When all the simulations are complete, the user can run an analysis script that descends into the file tree, gathering all the simulation output. It then combines this information to construct the milestoning model, and performs all the milestoning and error calculations, providing the user with kinetic and thermodynamic information, including kon, koff, and the free energy profile. It also has the option to perform convergence analysis on these values. Additional analysis scripts can be utilized to generate a single file containing the ligand equilibrium distribution or FHPD of each milestone for easy visualization.</p><!><p>The Python scripts have been tested using Python 2.7 and can be safely run in any version of Python 2 at version 2.7 or later. The remainder of the scripts are written in TCL, particularly those interfacing with NAMD, which has a TCL-based interface. SEEKR also uses the Numpy, Scipy, and MDAnalysis python libraries. The Adaptive Poisson-Boltzmann Solver (APBS)42 is used to generate the electrostatic potential maps for input to BrownDye, and the AmberTools program LEaP43 is also used to prepare structures for MD simulation.</p><!><p>Atomic coordinates of the trypsin-benzamidine system were obtained from the high resolution crystal structure Protein Databank (PDB) ID: 3PTB.44 Hydrogens were added using Molprobity with ring flips allowed.45–46 The system was then further prepared using LEaP with the Amber forcefield, ff14SB.47 Disulfide bonds were added manually. The appropriate protonation states of ASP, GLU, and HIS residues at a pH of 7.7 were determined using PROPKA.48–49 This pH was selected to align with the experimental conditions of Guillian and Thusius50. The structure was then solvated in a truncated octahedron of TIP4Pew51–52 waters and eight Cl− ions were added to neutralize the overall charge. The benzamidine ligand was parameterized using Antechamber with the GAFF force field.52–53 The total size of the system was approximately 23,000 atoms. To allow for relaxation from the crystallographic starting structure, the benzamidine ligand was removed and a 20 ns simulation of the apo structure was performed at a constant temperature of 298 K using the Langevin thermostat and a constant pressure of 1 atm using the Langevin piston with a damping coefficient of 5 ps−1. A representative structure from this simulation was then used as the SEEKR input structure to generate all the necessary inputs for the MD simulations to be run using NAMD, and the BD simulations using Browndye.</p><p>The benzamidine bound-state coordinates were defined from the center of mass of the alpha carbons of residues 190, 191, 192, 195, 213, 215, 216, 219, 220, 224, 228 of PDB: 3PTB because these residues form the binding pocket in the bound-state crystal structure by manual inspection. Spherical milestones were defined with radii of 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 14 Å, with the origin being the bound state coordinates defined above. This spacing of the milestones was chosen to facilitate the simulation of transitions between milestones while still ensuring the Markov assumptions required by formal milestoning theory. Ten copies of the apo structure were generated, each with the benzamidine ligand inserted on one of the ten spherical milestones (figure 2A). Water molecules that clashed with the ligand structure were removed. The first nine milestones correspond to the MD simulation regime, with the innermost milestone (1 Å) representing the bound state, as the center of mass of the bound benzamidine ligand falls well within the 1 Å sphere that defines this milestone (figure 2B). Furthermore, in a ~170 ns unrestrained MD simulation with the ligand in the bound pose, the 1 Å sphere contained the center of mass of the ligand over 71% of the simulation. The tenth and outermost milestone (14 Å) corresponds to the BD simulation regime. The distribution along any milestone where BD was started was constructed by first running conventional BD simulations and obtaining the distribution of hitting points along that milestone.</p><p>The b-surface is a relatively large spherical shell that encloses the entire receptor molecule, with a radius of sufficient size that the entire surface sits well out into the bulk solvent where forces between the ligand and receptor would be largely unaffected by molecular orientation, and are therefore centrosymmetric.</p><!><p>A modified version of NAMD 2.11 was used for all MD calculations. The numerous MD inputs, including input files, integrator parameters, boundary conditions, temperature and pressure controls, etc. are either defined by the user or set by SEEKR to default values. Relevant settings and procedures implemented for each milestone in the MD regime are described here.</p><p>For each milestone system generated by SEEKR as described above, the solvent molecules were allowed to relax around the newly placed benzamidine ligand by minimizing for 5000 steps with both the ligand and receptor restrained. The solvent was then further relaxed through a series of 2 ps heating simulations, where the temperature was increased from 298 K to 350 K and then cooled back to 298 K in 10 K increments, keeping the atoms of the ligand and receptor restrained. Following this relaxation of the solvent, an equilibrium distribution of the ligand on the milestone surface was obtained from 1 μs of constant volume simulation at a temperature of 298 K where a harmonic spring force of 90 kcal•mol−1•Å−2 was imposed to restrain the ligand at the appropriate radius from the binding site center for each milestone to generate an equilibrium distribution (figure 3A). This is also known as the umbrella sampling stage. From this equilibrium distribution, a FHPD (figure 3B) was obtained by selecting 4700 position and velocity configurations from times 60 ns – 1 μs of the equilibrium trajectory and allowing them to propagate backwards in time by reversing their velocities at constant energy and volume (reverse stage). Any trajectories that struck another milestone before re-crossing the milestone from which they originated were counted as part of the FHPD. All members of the FHPD were then brought back to their original positions and velocities and subsequently allowed to propagate forward in time at constant energy and volume (forward stage). When a simulation crossed its starting milestone again, it was then monitored for transitions to adjacent milestones and the incubation time for these transitions was also recorded. Once a trajectory crossed an adjacent milestone, the simulation was terminated. Any trajectories in this forward stage that crossed adjacent milestones before re-crossing their starting milestone were rejected. The 1, 1.5, 2, 4, and 10 Å milestones produced results with significantly fewer transitions than the other milestones. Therefore, to improve the robustness of our statistics, we performed additional reverse and forward simulations where 10 more trajectories were initiated at random Maxwell-Boltzmann velocities from each equilibrium distribution point, in addition to the one described above (a total of 470,000 reversals for each of these milestones), increasing the number of transitions observed. For each milestone, successful forward stage statistics were inserted into the transition kernel K and incubation time vector 〈t〉.</p><!><p>All BD calculations were conducted with BrownDye, a software package specializing in the rigid-body diffusion of two biological molecules in an implicit solvent.41 The electric potential map used as input for the BD simulation was calculated with the APBS version 1.4. All BD inputs, as well as the necessary APBS inputs for creation of the electrostatics map, are user defined in the SEEKR input file or generated as SEEKR default values.</p><p>In an attempt to recreate the ionic conditions used in the experiment,50 a nonlinear APBS calculation was run at 298 K, with a solvent dielectric of 78 and a solute dielectric of 2, with the following ions: Ca2+ at a concentration of 0.02 mM with a charge of +2.0 e and a radius of 1.14 Å, Cl− at a concentration of 0.10 mM with a charge of −1.0 e and a radius of 1.67 Å, and tris at a concentration of 0.06 mM with a charge of +1.0 e and a radius of 4.0 Å.54 At the specified concentrations, these ions generate a Debye length of 8 Å, which is used as input to BrownDye. Both the b-surface BD simulations and BD trajectories starting from a milestone ran with a solvent dielectric 78 and a solute dielectric of 2, at 298 K. We ran three additional sets of BD simulations at different ionic concentrations to examine the effect of ionic strength in the BD simulations on the kon. Therefore, three additional simulations were run: one with an ion concentration of zero, another with half of the ion concentrations of the experimental procedure, and another with double the ion concentration of the experimental procedure. Although an electrolyte solution technically has a Debye length equal to infinity, we approximated the Debye length with a value of 99 Å in the BrownDye program.</p><p>For each kon calculation, we performed 106 BD simulations initiated at random points distributed on the b-surface, which were used to construct the vector q0 in Eq. 13. Once these simulations completed, the trajectories that successfully reached the outermost milestone were used as that milestone's FHPD. From that FHPD, an additional 106 BD trajectories were run until reaching the second-outermost milestone or escaping to the q-surface. These statistics were also included in the transition kernel K and incubation time vector 〈t〉.</p><!><p>Using the statistics obtained from all the milestones in both the MD and BD regimes, the SEEKR software was used to construct the milestoning model and compute the kon, koff, ΔGbind, and other quantities of interest. Additional scripts used to generate some of the figures and data are also included in the SEEKR package. Error estimates were computed according to our previously defined procedure31.</p><p>The vast majority of the procedure outlined in the Materials and Methods section is automated within the SEEKR software package.</p><!><p>Using the MD/BD/milestoning methodology through the SEEKR interface yielded a kon of 2.1±0.3•107 M−1s−1 for the trypsin-benzamidine system. This value deviates from the experimentally measured kon for the same system at 2.9•107 M−1s−1 by a factor of ~1.5 (no experimental error margins were reported). We also estimate a koff of 83±14 s−1, which is within an order of magnitude of the experimentally determined value of 600±300 s−1 though our value is slower than expected. similar phenomenon is observed in other computational koff estimations of this system. An examination of the effect of ionic concentration on the kon convergence of the rate constants as a function of the length of umbrella sampling performed is provided in the SI. Using Eq. 5, we obtain a ΔGbind estimate of −7.3±0.2 kcal•mol−1 from a Kd of 4.3±1.2•10−6 M compared to the experimental ΔGbind of −6.71±0.05 kcal•mol−1, computed at 298 K using Eq. 5 and an experimental Kd of 1.2±0.1•10−5 M.50</p><p>In addition, we obtained a relative free energy at each of the milestones along the binding pathway using the vector pstat in combination with Eq. 10. This free energy profile is displayed in figure 4.</p><p>Aside from the predicted thermodynamic and kinetic quantities, we used the trajectories generated during the SEEKR run to make other observations about the system during the binding and unbinding process.</p><p>By removing the benzamidine molecule and the solvent, we used POVME255 to provide pocket volume measurement and characterization during the course of the MD runs. The same origin and radius of the inclusion region that defined the binding pocket were used for all umbrella sampling trajectories. The pocket itself remains relatively rigid when the benzamidine is deep in the binding site during the umbrella sampling stage, however, more variation in volume was observed when the benzamidine was constrained to a milestone nearer to the entrance of the opening of the binding site (figure 5).</p><p>Closer analysis of the umbrella sampling trajectories for the 6, 10, and 12 Å milestones in conjunction with the POVME data indicates sampling of multiple conformations of the trypsin S1 binding pocket (figure 6A, 6B, and 6C). The binding pocket conformation is primarily dependent on the motion of two loops; the loop containing TRP215 and the loop containing ASP189, a critical residue for benzamidine recognition. The opening and closing of the S1 pocket is greatly influenced by the orientation of TRP215 When oriented downward as in figure 6A, the S1 pocket is open. This is the conformation observed in the crystal structure 3PTB with benzamidine bound. When TRP215 rotates upwards as in figure 6B, the binding pocket is closed, and pocket volume significantly decreases. The dramatic change in pocket volume for the 10 Å milestone also occurs when TRP215 moves to close the S1 binding site.</p><p>We also observe the formation of an S1* pocket, that results from the motion of these two loops (figure 6C). This pocket provides an alternate binding pathway, in which benzamidine can approach ASP189 from a different orientation. These observations are in agreement with the study of Plattner and Noé22 where these results were observed through several hundred independent MD trajectories totaling over 100 μs of aggregate simulation time.</p><p>We also observed significant positional and rotational sampling by the benzamidine along most of the milestones during the umbrella sampling stages. This information can provide an idea for the likelihood of pathways that benzamidine follows on its route to binding. Figure 3A shows the equilibrium distribution along each of the milestones, and figure 3B shows the FHPD for each of the milestones. Figure 7 shows the angle between a vector pointing along the amidine group and a vector pointing out from the opening of the binding site as a function of time during the equilibrium simulations. Several flips are observed in all but the lowest milestones, where benzamidine rotation was restricted because these milestones are located deep within the binding pocket. The 10 Å also experiences a decrease in rotational sampling because benzamidine is interacting extensively with TRP215 and thus adopts an orientation that favors stacking of the aromatic rings.</p><p>The crystal structure of the trypsin/benzamidine complex shows the amidine group pointing downward toward the binding site (figure 2B). This structural feature is confirmed by our own simulations, and a relatively narrow arrangement of ligand orientations are observed along the lowest milestone.</p><p>The entire calculation cost approximately 1.4 million CPU hours on the Stampede supercomputer and local machines, with a total MD cost of approximately 19 μs of simulation.</p><!><p>Compared to the experimental kon, our estimated kon is slower by about a factor of 1.3, but falls well within an order of magnitude. We attempted to closely recreate the experimental ionic conditions within our simulations, which has a pronounced effect on the kon (details in the SI). Our kon of 2.2±0.3•107 M−1s−1 is much closer to the experimental value of 2.9•107 M−1s−1 than the kons obtained by Buch et. al.20 (15±2•107 M−1s−1) and comparable to what was obtained by Plattner et. al.22 (6.4±1.6•107 M−1s−1), although ours was obtained with significantly less computational resources, smaller by an order of magnitude. Our result is also very close to what was obtained by Tiwary et. al.25 (1±1•107 M−1s−1). Our estimated koff of 83±14 s−1 is within an order of magnitude of the experimental koff, far closer than the value obtained by Buch et. al. (9.5±3.3•104s−1), and comparable to the values obtained by Plattner et. al. (131±109 s−1), Teo et. al.24 (260±240 s−1), and Tiwary et. al (9.1±2.5 s−1). To our knowledge, this is the first successful estimate of koff using a hybrid MD/BD/milestoning model.</p><p>An advantage of our approach is that both koff and kon can be determined from the same calculation. We can use our calculated koff and kon values in Eq. 5 to obtain an entirely computationally-determined dissociation constant KD of 3.8±0.8•10−6 M, and by extension a free energy of binding ΔGbind estimate of −7.4±0.2 kcal•mol−1. This is in good agreement with the experimental KD of 1.2•10−5 M, which when put through eq. 5 at a temperature of 298 K, yields a free energy of −6.7 kcal•mol−1.</p><p>The accurate determination of kinetics using milestoning requires the proper generation of equilibrium and FHPD distributions. It is important to ensure adequate sampling in the generation of equilibrium distributions. Figure 3A shows the equilibrium distribution of benzamidine center-of-mass along the 1 Å to 12 Å milestones in the MD regime. The benzamidine appears to have explored all solvent-accessible regions along the milestones. Along with positional sampling, the observed diversity of benzamidine orientation in figure 7 indicates that the ligand orientational degree of freedom is well-sampled in all but the lowest milestones. In addition to the ligand, it is important that receptor conformations that may affect ligand binding are also well sampled. By using POVME2, we observed conformational states that have been observed in other studies such as the S1* pocket (figure 6).22 We do not however observe any complete binding events via the S1* pocket, presumably as a result of our simplified spherical milestoning model. This may provide some explanation as to why our calculated rates are somewhat slower than experiment, as we do not capture this alternate pathway. However, we may reasonably assume that we are capturing most of the effects of slower receptor conformational changes and subsequently, that our kinetics predictions are reasonable.</p><p>While, of course, verification of SEEKR as a computational kinetics and thermodynamics estimator will need to be performed on additional systems, this similarity between experimental and theoretical free energies and rate constants in our accessible and highly parallel framework is encouraging.</p><!><p>In this work, we use our multi-scale MD/BD/milestoning methodology to examine ligand-protein binding events with a larger, more complex, and more drug-like ligand than in our previous work. Furthermore, we present the first successful koff calculation to within one order of magnitude of experiment using this approach. Using the obtained values of kon and koff, and entirely computational estimate of KD and ΔGbind in good agreement with experiment were obtained. These results are further evidence that the MD/BD/milestoning methodology can be successfully applied to the investigation of binding and unbinding kinetics in receptor-ligand systems. We also present the SEEKR software package, which automates much of the preparation, submission, and analysis of these types of calculations. We have made SEEKR freely available and open-source on Github, and hope that it will be used and improved by the community to run predictive multi-scale milestoning calculations. SEEKR downloads, tutorials, and the user guide may be found at http://amarolab.ucsd.edu/seekr.</p>