Because of the observation that kinetic rates can be better correlated with in vivo drug efficacy than simple affinities, significant efforts have been made to develop experimental and computational methods to predict drug-target residence times over the past few years. Here, we summarize the discussions and reflections from an international Centre Européen de Calcul Atomique et Moléculaire (CECAM) Workshop, which took place in March 2025, with participants from both academia and industry, including experimentalists as well as computational scientists, and was dedicated to computational methods for protein-ligand binding kinetics prediction. As currently standing challenges, we identify the need for standardized benchmark datasets, the need to move from simple to more complex and biomedically relevant targets and the question of how academia and industry can work together to move the field of drug-target binding kinetics forward. We also discuss what level of accuracy can be expected from computational methods and highlight that the field would benefit from blind challenges to enable a fair comparison of different computational methods to predict kinetic rates.

Lipocalins are a family of proteins found in mammals that are essential for the binding and transport of various molecules, but the mechanisms underlying their target recognition are still unclear. To answer this question, we studied odorant-binding proteins (OBPs), a specific type of lipocalin involved in chemical communication and olfaction. Using an integrative approach combining numerical modelling and experimental validation, we identified key structural regions that regulate the entry of molecules into the binding pocket. Modification of these regions disrupts molecular recognition, highlighting their importance for function. In addition, we found that changes in distant parts of the protein influence binding, shedding light on allosteric mechanisms. These results advance our understanding of lipocalin function and open up avenues for the design of proteins with targeted binding properties.

Predicting the molecular friction and energy landscapes under nonequilibrium conditions is key to coarse-graining the dynamics of selective solute transport through complex, fluctuating, and responsive media, e.g., polymeric materials such as hydrogels, cellular membranes, or ion channels. The analysis of equilibrium ensembles already allows such a coarse-graining for very mild nonequilibrium conditions. However, in the presence of stronger external driving and/or inhomogeneous setups, the transport process is governed apart from a potential of mean force also by a nontrivial position- and velocity-dependent friction. It is therefore important to find suitable and efficient methods to estimate the mean force and the friction landscape, which can then be used in a low-dimensional, coarse-grained Langevin framework to predict the system's transport properties and timescales. In this work, we evaluate different coarse-graining approaches based on constant-velocity constraint simulations for generating such estimates using two model systems, which are a 1D responsive barrier as a minimalistic model and a single tracer driven through a 3D bead-spring polymer membrane as a more sophisticated problem. Finally, we demonstrate that the estimates from 3D constant-velocity simulations yield the correct velocity-dependent friction, which can be directly utilized for coarse-grained (1D) Langevin simulations with constant external driving forces.

We analyze the coarse-grained equations of motion of molecular systems subject to external driving. As exemplary processes, we study by means of targeted and steered molecular dynamics simulations the dissociation of a sodium-chloride ion pair in water and ligand-protein unbinding of trypsin-benzamidine. We derive an exact generalization of Mori's Langevin equation that contains the memory kernel of the stationary process, an additive driving force, and a non-equilibrium response force describing the effects of the perturbed environment. We show that both the fluctuating force in the stationary case and the non-equilibrium response force in the driven cases exhibit spatial structure in their first and second moments. The latter depends sensitively on the employed driving protocols. For sodium chloride, we find that the first moment of the non-equilibrium response force matches the mean force for slow constrained pulling. In contrast, for all tested restrained pulling protocols, significant differences arise between the two properties in both systems. We conclude that the non-equilibrium response of the solvent needs to be taken into account carefully when analyzing data from pulling simulations.

Finding process pathways in molecular simulations such as the unbinding paths of small molecule ligands from their binding sites at protein targets in a set of trajectories via unsupervised learning approaches requires the definition of a suitable similarity measure between trajectories. Here, we evaluate the performance of four such measures with varying degree of sophistication, i.e., Euclidean and Wasserstein distances, Procrustes analysis, and dynamic time warping, when analyzing trajectory data from two different biased simulation driving protocols in the form of constant velocity constraint targeted MD and steered MD. In a streptavidin-biotin benchmark system with known ground truth clusters, Wasserstein distances yielded the best clustering performance, closely followed by Euclidean distances, both being the most computationally efficient similarity measures. In a more complex A2a receptor-inhibitor system, however, the simplest measure, i.e., Euclidean distances, was sufficient to reveal meaningful and interpretable clusters.

Protein–ligand (un)binding kinetics have been demonstrated to correlate better with drug efficacy than protein–ligand affinity. Consequently, the prediction of such kinetics via molecular dynamics simulations is of recent interest for pharmaceutical research. In this chapter, the theoretical basis for the calculation of such kinetics as well as biased molecular dynamics methods aiming for such predictions are introduced. The challenges involved in such predictions are highlighted and compared to binding affinity calculations. State-of-the-art of kinetics-predicting simulation methods are reviewed concerning both their capabilities and shortcomings.

Understanding the dynamics of biomolecular complexes, e.g., of protein-ligand (un)binding, requires the comprehension of paths such systems take between metastable states. In MD simulations, paths are usually not observable per se, but they need to be inferred from simulation trajectories. Here, we present a novel approach to cluster trajectories based on a community detection algorithm that necessitates only the definition of a single parameter. The unbinding of the streptavidin-biotin complex is used as a benchmark system and the A2a adenosine receptor in complex with the inhibitor ZM241385 as an elaborate application. We demonstrate how such clusters of trajectories correspond to pathways and how the approach helps in the identification of reaction coordinates for a considered (un)binding process.

To sample rare events, dissipation-corrected targeted molecular dynamics (dcTMD) applies a constant velocity constraint along a one-dimensional reaction coordinate s, which drives an atomistic system from an initial state into a target state. Employing a cumulant approximation of Jarzynski's identity, the free energy ΔG(s) is calculated from the mean external work and dissipated work of the process. By calculating the friction coefficient Γ(s) from the dissipated work, in a second step, the equilibrium dynamics of the process can be studied by propagating a Langevin equation. While so far dcTMD has been mostly applied to study the unbinding of protein-ligand complexes, here its applicability to rare conformational transitions within a protein and the prediction of their kinetics are investigated. As this typically requires the introduction of multiple collective variables {xj} = x, a theoretical framework is outlined to calculate the associated free energy ΔG(x) and friction Γ(x) from dcTMD simulations along coordinate s. Adopting the α-β transition of alanine dipeptide as well as the open-closed transition of T4 lysozyme as representative examples, the virtues and shortcomings of dcTMD to predict protein conformational transitions and the related kinetics are studied.

The prediction of drug-target binding and unbinding kinetics that occur on time scales between milliseconds and several hours is a prime challenge for biased molecular dynamics simulation approaches. This Perspective gives a concise summary of the theory and the current state-of-the-art of such predictions via biased simulations, of insights into the molecular mechanisms defining binding and unbinding kinetics as well as of the extraordinary challenges predictions of ligand kinetics pose in comparison to binding free energy predictions.

Protein dynamics have been investigated on a wide range of time scales. Nano- and picosecond dynamics have been assigned to local fluctuations, while slower dynamics have been attributed to larger conformational changes. However, it is largely unknown how fast (local) fluctuations can lead to slow global (allosteric) changes. Here, fast molecule-spanning dynamics on the 100 to 200 ns time scale in the heat shock protein 90 (Hsp90) are shown. Global real-space movements are assigned to dynamic modes on this time scale, which is possible by a combination of single-molecule fluorescence, quasi-elastic neutron scattering and all-atom molecular dynamics (MD) simulations. The time scale of these dynamic modes depends on the conformational state of the Hsp90 dimer. In addition, the dynamic modes are affected to various degrees by Sba1, a co-chaperone of Hsp90, depending on the location within Hsp90, which is in very good agreement with MD simulations. Altogether, this data is best described by fast molecule-spanning dynamics, which precede larger conformational changes in Hsp90 and might be the molecular basis for allostery. This integrative approach provides comprehensive insights into molecule-spanning dynamics on the nanosecond time scale for a multi-domain protein.

The effect of an externally applied directional force on molecular friction is so far poorly understood. Here, we study the force-driven dissociation of the ligand-protein complex biotin-streptavidin and identify anisotropic friction as a not yet described type of molecular friction. Using AFM-based stereographic single molecule force spectroscopy and targeted molecular dynamics simulations, we find that the rupture force and friction for biotin-streptavidin vary with the pulling angle. This observation holds true for friction extracted from Kramers’ rate expression and by dissipation-corrected targeted molecular dynamics simulations based on Jarzynski’s identity. We rule out ligand solvation and protein-internal friction as sources of the angle-dependent friction. Instead, we observe a heterogeneity in free energy barriers along an experimentally uncontrolled orientation parameter, which increases the rupture force variance and therefore the overall friction. We anticipate that anisotropic friction needs to be accounted for in a complete understanding of friction in biomolecular dynamics and anisotropic mechanical environments.

Protein-ligand (un)binding simulations are a recent focus of biased molecular dynamics simulations. Such binding and unbinding can occur via different pathways in and out of a binding site. Here, we present a theoretical framework on how to compute kinetics along separate paths and on how to combine the path-specific rates into global binding and unbinding rates for comparison with experimental results. Using dissipation-corrected targeted molecular dynamics in combination with temperature-boosted Langevin equation simulations [S. Wolf et al., Nat. Commun. 11, 2918 (2020)] applied to a two-dimensional model and the trypsin-benzamidine complex as test systems, we assess the robustness of the procedure and discuss the aspects of its practical applicability to predict multisecond kinetics of complex biomolecular systems.

The friction coefficient of fluids may become a function of the velocity at increased external driving. This non-Newtonian behavior is of general theoretical interest and of great practical importance, for example, for the design of lubricants. Although the effect has been observed in large-scale atomistic simulations of bulk liquids, its theoretical formulation and microscopic origin are not well understood. Here, we use dissipation-corrected targeted molecular dynamics, which pulls apart two tagged liquid molecules in the presence of surrounding molecules, and analyze this nonequilibrium process via a generalized Langevin equation. The approach is based on a second-order cumulant expansion of Jarzynski's identity, which is shown to be valid for fluids and therefore allows for an exact computation of the friction profile as well of the underlying memory kernel. We show that velocity-dependent friction in fluids results from an intricate interplay of near-order structural effects and the non-Markovian behavior of the friction memory kernel. For complex fluids such as the model lubricant C40H82, the memory kernel exhibits a stretched-exponential long-time decay, which reflects the multitude of timescales of the system.

Allosteric communication between distant protein sites represents a key mechanism of biomolecular regulation and signal transduction. Compared to other processes such as protein folding, however, the dynamical evolution of allosteric transitions is still not well understood. As example of allosteric coupling between distant protein regions, we consider the global open-closed motion of the two domains of T4 lysozyme, which is triggered by local motion in the hinge region. Combining extensive molecular dynamics simulations with a correlation analysis of interresidue contacts, we identify a network of interresidue distances that move in a concerted manner. The cooperative process originates from a cogwheel-like motion of the hydrophobic core in the hinge region, which constitutes a flexible transmission network. Through rigid contacts and the protein backbone, the small local changes of the hydrophobic core are passed on to the distant terminal domains and lead to the emergence of a rare global conformational transition. As in an Ising-type model, the cooperativity of the allosteric transition can be explained via the interaction of local fluctuations.

Photoproteins such as bacteriorhodopsin (bR) and rhodopsin (Rho) need to effectively dissipate photoinduced excess energy to prevent themselves from damage. Another well-studied seven transmembrane (TM) helices protein is the β2 adrenergic receptor (β2AR), a G protein-coupled receptor for which energy dissipation paths have been linked with allosteric communication. To study the vibrational energy transport in the active and inactive states of these proteins, a master equation approach [J. Chem. Phys.2020, 152, 045103] is employed, which uses scaling rules that allow us to calculate energy transport rates solely based on the protein structure. Despite their overall structural similarity, the three 7TM proteins reveal quite different strategies to redistribute excess energy. While bR quickly removes the energy using the TM7 helix as a "lightning rod", Rho exhibits a rather poor energy dissipation, which might eventually require the hydrolysis of the Schiff base between the protein and the retinal chromophore to prevent overheating. Heating the ligand adrenaline of β2AR, the resulting energy transport network of the protein is found to change significantly upon switching from the active state to the inactive state. While the energy flow may highlight aspects of the inter-residue couplings of β2AR, it seems not particularly suited to explain allosteric phenomena.

While allostery is of paramount importance for protein signaling and regulation, the underlying dynamical process of allosteric communication is not well understood. The PDZ3 domain represents a prime example of an allosteric single-domain protein, as it features a well-established long-range coupling between the C-terminal α3-helix and ligand binding. In an intriguing experiment, Hamm and co-workers employed photoswitching of the α3-helix to initiate a conformational change of PDZ3 that propagates from the C-terminus to the bound ligand within 200 ns. Performing extensive nonequilibrium molecular dynamics simulations, the modeling of the experiment reproduces the measured time scales and reveals a detailed picture of the allosteric communication in PDZ3. In particular, a correlation analysis identifies a network of contacts connecting the α3-helix and the core of the protein, which move in a concerted manner. Representing a one-step process and involving direct α3-ligand contacts, this cooperative transition is considered as the elementary step in the propagation of conformational change.