Predicting Drug Residence Times from Unbinding Kinetics Simulations

The unbinding kinetics workflow script automatically performs the necessary tasks for setting up the system for MD, such as running the Desmond System Builder, setting up the appropriate equilibration protocol, and launching the job.

The parameters of the workflow generalize well across systems and detailed tweaking should not be necessary in most cases. Unless working with ligands with very long expected residence times, we recommend lowering the residence time cutoff from its default of 10¹² s to 10¹⁰ s to reduce the computational time spent on replicas which seem to be stuck.

The unbinding kinetics script workflow takes additional optional arguments that may be helpful in specific cases, such as adding a lipid bilayer membrane or tweaking the simulation parameters. See Running Unbinding Kinetics from the Command Line for an explanation of the most common arguments and the unbinding_kinetics command help for the full list.

As mentioned above, the workflow’s default parameters should work well in most cases. Here are a few heuristics from our experience for when it makes sense to deviate from the defaults:

• For large ligands or those with very long residence times, you should increase the number of replicas for the metadynamics stage (-total-mtd-replicas) and each replica’s maximum duration (-time, in picoseconds).
• For fragment-sized or weakly-binding ligands you may want to reduce the initial force constant for the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. simulations from 15 kcal/mol/A² to 5 kcal/mol/A². To do this, add -ramd-spring-fk 5.0 to the command above. See the publication for more details.
• For systems where you suspect there are multiple competing exit channels, increasing the number of replicas for the exploration stage (-total-ramd-replicas) can help improve the sampling of unbinding paths. The barrier wall applied to the ligand during the metadynamics stage is quite soft, so the workflow should be able to recover the correct statistics even if a suboptimal CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations is chosen in the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage (within reason).

When troubleshooting or exploring the system’s response to changing parameters, keep in mind that the metadynamics stage is much more computationally intensive than the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. exploration stage. You can run only the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage by appending -ramd-only to the launch command and inspect the results up to that stage as described in the next section to see how the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. parameters affect your system before re-running the job with the final parameters.

Please note that calculations may have been performed with a different version of the software than you are using and the results may not be exactly the same as those you produced in this tutorial. Especially for methods with a statistical component, such as MD simulations, an exact reproduction of a trajectory and its properties is only possible if the same random seed, software and hardware configuration are used. With sufficient sampling, all properties of the ensemble should be reproducible, while the properties of individual frames and the exact course of a trajectory may differ.

The unbinding kinetics workflow creates a directory for each of the steps of the protocol: initial structure (1), equilibration (2), RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. setup (3), RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. simulations (4), pathway clustering (5), metadynamics simulations (6), final predictions (7). The folder structure and each stage’s results are explained here. The results for many of the stages provide both a directory with the most commonly useful files and a .tgz archive with more in-depth data. For the sake of keeping the file size low, the tutorial zip archive only contains the full results for the later stages of the workflow discussed below.

For a successful run, the most practically relevant outputs are from stage 5 (unbinding pathway clustering results) and stage 7 (final predictions). The files from other stages can be useful when troubleshooting job failures or unsuccessful runs of the protocol.

The unbinding path summary shows the unbinding paths observed during the exploration stage in the various replicas as spheres tracing the center of mass of the ligand as it unbinds. Note that these paths are not pre-determined but result from applying force to displace the ligand in whichever direction it can be displaced. The paths are colored by which cluster they belong to and the clusters are listed in descending number of paths in the cluster.

In this case, there are multiple directions in which the ligand can exit the binding pocket. The most populated cluster (light blue spheres, 11 paths) clearly corresponds to the most dominant direction for unbinding, as it is also very similar to the dark blue, light green, and pink clusters (which together account for another 14 paths). The light blue cluster will therefore be used to construct the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations for the next part of the simulation. Note that the biasing potential during the metadynamics stage will not be very steep, initially more similar to giving the ligand a nudge in a particular direction than to keeping it on railroad tracks. So the ligand will still be quite free to move around while unbinding, making the method less reliant on the accuracy of the initial choice of CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations.

Due to the inherently stochastic nature of MD simulations and the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. approach, the exact number of paths per cluster will fluctuate between runs of the workflow. If you observe very even splits between clusters representing very different unbinding paths, you may want to extend the sampling in the exploration phase by running the above command with the -ramd-only and -total-ramd-replicas <number> flags to increase the number of replicas for the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage from the default of 30. Compared to the metadynamics stage, the exploration stage is fairly computationally cheap. Once you’re confident that the sampling is robust, you can re-run the job with the original command, adding any optional flags you’ve decided to use.

In case of a tie, the workflow uses each of the most populated clusters to determine a separate CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations and runs metadynamics on both unbinding paths. After the metadynamics simulations complete, the results and quality metrics from both paths are compared to determine the final result.

The most important metric to check here are the simulation times (given both individually for the replicas and as a summary). The force applied to the ligand is increased over the course of the simulation to force it out of the pocket. Once the distance between the ligand and the receptor C_α atoms reaches 8.0 Å, it is declared as unbound and the simulation is automatically terminated. Optimally, this should occur within 1-2 ns of RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. simulation. Our Ligand 1 has a very short residence time of only a few seconds, so the simulation time average of 300 ps across the replicas is reasonable and expected.

Shorter unbinding times (<100 ps) are likely for fragment-like or weakly-bound ligands, in which case you may want to re-run the simulation with a smaller force constant to avoid immediate unbinding and improve sampling (see above). If the simulation is taking very long, you can increase the force constant, but you should also have a more detailed look at what might prevent the ligand from unbinding (e.g. an incorrectly reconstructed part of the structure or an occluded binding pocket).

This file reports the predicted residence time in seconds, as well as quality metrics to gauge the reliability of the results. In our case, the residence time prediction of 0.664 ± 0.284 s is within one order of magnitude of the experimental value of 7.14 s, i.e. within the error bars of most kinetics measurements.

Of the quality metrics, the most important is the ratio of simulation replicas in which the ligand successfully unbinds (#Unbound.Sims) to the total number of replicas (#Sims), i.e. the success rate. The p-value can be a helpful secondary quality metric.

Unbinding success rate:
Once the calculated residence time reaches 10¹² s, the simulation is terminated. Failure to unbind is usually a result of emerging degrees of freedom that are orthogonal to the path collective variable defined in the exploration stage. Example causes are large receptor rearrangements, changing hydration patterns, or highly flexible ligands. The ratio of unsuccessful simulations for a given ligand should be <10%, and is a robust metric for how confident you can be of the overall results and residence time predictions.

If you see unsuccessful simulations, you should check the <jobname>_7-out/job.log (after extracting the <jobname>_7-out.tgz archive) for the Status of each replica to understand why the simulation was unsuccessful. The possible statuses are:

• “OK” for a successful simulation,
• “Too slow” for simulations where the maximum residence time cutoff (by default 10¹² s) is reached, or
• “Bound” for simulations where the maximum simulation time (by default 200 ns) is reached.

If you observe unsuccessful runs due to exceeding the maximum simulation time, this does not mean that the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations is bad, as this cutoff is motivated by conserving computational resources rather than any scientific reasons. You can recover the correct statistics by increasing this cutoff (-time) – see the next section of the tutorial for an example where this is necessary. If replicas fail due to the residence time cutoff, this is an indication that you should investigate the system in more detail.

In this case, all 10 replicas result in a successful unbinding, which means we can be quite confident in the results.

p-value:
In a nutshell, this value shows how closely the observed unbinding resembles a well-sampled first-order kinetics process. For an intuitive explanation of what this means, consider the following plot (automatically generated by the workflow as <jobname>_7/<jobname>.png):

The blue curve shows the sampling we’ve performed with the successful replicas, the orange curve corresponds to an ideal first-order kinetics process. The p-value describes how similar the blue curve is to the orange curve.

The p-value should always be interpreted together with the unbinding success rate, resulting in three possible cases:

• High success rate (>90%) and high p-value (>>0.05) indicate the ideal scenario of good sampling of a singular transition state.
• High success rate and low p-value can occur when the unbinding event does not actually follow first-order kinetics, when egress paths mix, or other factors skew the statistics. This is arguably the most frequently observed scenario, and your model will likely still be predictive. Comparing the model’s predictions to experimental data (if available) is always recommended.
  In this case, you can have a closer look at the distribution of residence times between the replicas by checking the <jobname>_7-out/job.log (after extracting the <jobname>_7-out.tgz archive) for any outliers.
• If the success rate is low, the sampling is poor in any case and the p-value is irrelevant.
  In this case, you need to investigate your system in more detail, in particular the exit channel. If the channel corresponding to the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations is blocked but you suspect another channel exists, you could start over from the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage, increasing the sampling as described above to find new exit paths or make sure the statistics are stable. If the channel is truly occluded without significant receptor reorganization, your system may be beyond the domain of applicability of the workflow.

See the publications by Smith et al and Salvalaglio et al for more details on the p-value as a metric of metadynamics reliability.

In this example, our results show a p-value of 0.059, which together with our high success rate indicates that the simulation results are reasonable.

As mentioned above, in the case that multiple path clusters from the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage are tied in population, the metadynamics stage is run for all corresponding collective variables. All sets of statistics are reported in the <jobname>.dat file. The final result is chosen by sorting the CVs into  confidence categories based on the quality metrics (success rate ≶ 70% and/or p ≶ 0.05) and choosing the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations from the highest-confidence category which has the shortest predicted residence time.

The MetaD: p38-Lig1 group contains the ligand-receptor structure, the path CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations (path_cv_0) shown as spheres tracking the ligand center of mass along the path found in the exploration stage, and a sub-group (path_cv_0_mtd) tracking the actual center of mass of the ligand in each replica where it successfully unbinds. Each Entry from the path_cv_0_mtd group is named with the index of the corresponding replica.

Note that the predicted residence time and p-value are also visible in the Maestro Project Table as the residence time sec and KS p value properties of the path_cv_0 Entry.

You can see that in some replica (e.g. mtd_0-05), the ligand actually takes an egress path that is closer to one of the less-populated clusters from the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage than to the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations used in the metadynamics simulation. This is intentional and results from the softness of the biasing potential applied to the ligand during metadynamics to keep it near the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations. The ligand is thus free to explore the available space and even exit through another pore if it “wants to”. It is normal for the ligand to go off the rails like this in a few replicas, but if most of the trajectories show the ligand passing through a different pore than the one defined by the CVcollective variable, a term for the abstract coordinate used to bias the system during metadynamics simulations, this is a warning sign that the sampling in the RAMDrandom acceleration molecular dynamics; a biased molecular dynamics simulation in which an accelerating force is periodically applied to part of the system to displace it. stage may have been insufficient.

You can investigate the unbinding trajectory of any of the replicas in detail (e.g. to identify transient ligand-receptor interactions). You can find the full trajectories of all replicas by extracting the <jobname>_6-out.tgz archive in the job directory.

In summary, you can see that ligands such as BIRB796 require more thorough sampling in both the RAMD and metadynamics stages, resulting in higher computational costs to obtain an accurate prediction. For reference, here is the full command to launch the unbinding kinetics workflow with the recommendations listed above for this trickier case:

$SCHRODINGER/unbinding_kinetics <input>.mae -total-ramd-replicas 50 -total-mtd-replicas 20 -time 600000 -HOST <your_cpu_host> -SUBHOST <your_gpu_host> -OPLSDIR <your_opls_dir> -JOBNAME <jobname>

When troubleshooting your own system, keep in mind that the metadynamics stage is much more computationally costly than the RAMD stage, so testing the effect of changing parameters or improving the sampling is much more efficient if done for the RAMD stage.