Kinetics for Drug Discovery: an industry-driven effort to target drug residence time

Imidazole[1,5-a]pyridine derivatives as EGFR tyrosine kinase inhibitors unraveled by umbrella sampling and steered molecular dynamics simulations

Although the use of the tyrosine kinase inhibitors (TKIs) has been proved that it can save live in a cancer treatment, the currently used drugs bring in many undesirable side-effects. Therefore, the search for new drugs and an evaluation of their efficiency are intensively carried out. Recently, a series of eighteen imidazole[1,5-a]pyridine derivatives were synthetized by us, and preliminary analyses pointed out their potential to be an important platform for pharmaceutical development owing to their promising actions as anticancer agents and enzyme (kinase, HIV-protease,…) inhibitors. In the present theoretical study, we further analyzed their efficiency in using a realistic scenario of computational drug design. Our protocol has been developed to not only observe the atomistic interaction between the EGFR protein and our 18 novel compounds using both umbrella sampling and steered molecular dynamics simulations, but also determine their absolute binding free energies. Calculated properties of the 18 novel compounds were in detail compared with those of two known drugs, erlotinib and osimertinib, currently used in cancer treatment. Inspiringly the simulation results promote three imidazole[1,5-a]pyridine derivatives as promising inhibitors into a further step of clinical trials.

Ligand Gaussian accelerated Molecular Dynamics 3 (LiGaMD3): Improved Calculations of Binding Thermodynamics and Kinetics of Both Small Molecules and Flexible Peptides

Binding thermodynamics and kinetics play critical roles in drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics of small molecules and flexible peptides using conventional Molecular Dynamics (cMD), due to limited simulation timescales. Based on our previously developed Ligand Gaussian accelerated Molecular Dynamics (LiGaMD) method, we present a new approach, termed "LiGaMD3", in which we introduce triple boosts into three individual energy terms that play important roles in small-molecule/peptide dissociation, rebinding and system conformational changes to improve the sampling efficiency of small-molecule/peptide interactions with target proteins. To validate the performance of LiGaMD3, MDM2 bound by a small molecule (Nutlin 3) and two highly flexible peptides (PMI and P53) were chosen as model systems. LiGaMD3 could efficiently capture repetitive small-molecule/peptide dissociation and binding events within 2 microsecond simulations. The predicted binding kinetic constant rates and free energies from LiGaMD3 agreed with available experimental values and previous simulation results. Therefore, LiGaMD3 provides a more general and efficient approach to capture dissociation and binding of both small-molecule ligand and flexible peptides, allowing for accurate prediction of their binding thermodynamics and kinetics.

Accurate Characterization of Binding Kinetics and Allosteric Mechanisms for the HSP90 Chaperone Inhibitors Using AI-Augmented Integrative Biophysical Studies

The binding kinetics of drugs to their targets are gradually being recognized as a crucial indicator of the efficacy of drugs in vivo, leading to the development of various computational methods for predicting the binding kinetics in recent years. However, compared with the prediction of binding affinity, the underlying structure and dynamic determinants of binding kinetics are more complicated. Efficient and accurate methods for predicting binding kinetics are still lacking. In this study, quantitative structure-kinetics relationship (QSKR) models were developed using 132 inhibitors targeting the ATP binding domain of heat shock protein 90α (HSP90α) to predict the dissociation rate constant (koff), enabling a direct assessment of the drug-target residence time. These models demonstrated good predictive performance, where hydrophobic and hydrogen bond interactions significantly influence the koff prediction. In subsequent applications, our models were used to assist in the discovery of new inhibitors for the N-terminal domain of HSP90α (N-HSP90α), demonstrating predictive capabilities on an experimental validation set with a new scaffold. In X-ray crystallography experiments, the loop-middle conformation of apo N-HSP90α was observed for the first time (previously, the loop-middle conformation had only been observed in holo-N-HSP90α structures). Interestingly, we observed different conformations of apo N-HSP90α simultaneously in an asymmetric unit, which was also observed in a holo-N-HSP90α structure, suggesting an equilibrium of conformations between different states in solution, which could be one of the determinants affecting the binding kinetics of the ligand. Different ligands can undergo conformational selection or alter the equilibrium of conformations, inducing conformational rearrangements and resulting in different effects on binding kinetics. We then used molecular dynamics simulations to describe conformational changes of apo N-HSP90α in different conformational states. In summary, the study of the binding kinetics and molecular mechanisms of N-HSP90α provides valuable information for the development of more targeted therapeutic approaches.

Accelerating Molecular Dynamics Simulations for Drug Discovery

Accurate prediction of ligand binding thermodynamics and kinetics is crucial in drug design. However, it remains challenging for conventional molecular dynamics (MD) simulations due to sampling issues. Gaussian accelerated MD (GaMD) is an enhanced sampling method that adds a harmonic boost to overcome energy barriers, which has demonstrated significant benefits in exploring protein-ligand interactions. Especially, the ligand GaMD (LiGaMD) applies a selective boost potential to the ligand nonbonded potential energy, significantly improving sampling for ligand binding and dissociation. Furthermore, a selective boost potential is applied to the potential of both ligand and protein residues around binding pocket in LiGaMD2 to further increase the sampling of protein-ligand interaction. LiGaMD and LiGaMD2 simulations could capture repetitive ligand binding and unbinding events within microsecond simulations, allowing to simultaneously characterize ligand binding thermodynamics and kinetics, which is expected to greatly facilitate drug design. In this chapter, we provide a brief review of the status of LiGaMD in drug discovery and outline its usage.

Attempting Well-Tempered Funnel Metadynamics Simulations for the Evaluation of the Binding Kinetics of Methionine Aminopeptidase-II Inhibitors

Understanding the unbinding kinetics of protein-ligand complexes is considered a significant approach for the design of ligands with desired specificity and safety. In recent years, enhanced sampling methods have emerged as effective tools for studying the unbinding kinetics of protein-ligand complexes at the atomistic level. MetAP-II is a target for the treatment of cancer for which not a single effective drug is available yet. The identification of the dissociation rate of ligands from the complexes often serves as a better predictor for in vivo efficacy than the ligands' binding affinity. Here, funnel-based restraint well-tempered metadynamics simulations were applied to predict the residence time of two ligands bound to MetAP-II, along with the ligand association and dissociation mechanism involving the identification of the binding hotspot during ligand egress. The ligand-egressing route revealed by metadynamics simulations also correlated with the identified pathways from the CAVER analysis and by the enhanced sampling simulation using PLUMED. Ligand 1 formed a strong H-bond interaction with GLU364 estimating a higher residence time of 28.22 ± 5.29 ns in contrast to ligand 2 with a residence time of 19.05 ± 3.58 ns, which easily dissociated from the binding pocket of MetAP-II. The results obtained from the simulations were consistent to reveal ligand 1 being superior to ligand 2; however, the experimental data related to residence time were close for both ligands, and no kinetic data were available for ligand 2. The current study could be considered the first attempt to apply an enhanced sampling method for the evaluation of the binding kinetics and thermodynamics of two different classes of ligands to a binuclear metalloprotein.

Elucidating protein-ligand binding kinetics based on returning probability theory

The returning probability (RP) theory, a rigorous diffusion-influenced reaction theory, enables us to analyze the binding process systematically in terms of thermodynamics and kinetics using molecular dynamics (MD) simulations. Recently, the theory was extended to atomistically describe binding processes by adopting the host-guest interaction energy as the reaction coordinate. The binding rate constants can be estimated by computing the thermodynamic and kinetic properties of the reactive state existing in the binding processes. Here, we propose a methodology based on the RP theory in conjunction with the energy representation theory of solution, applicable to complex binding phenomena, such as protein-ligand binding. The derived scheme of calculating the equilibrium constant between the reactive and dissociate states, required in the RP theory, can be used for arbitrary types of reactive states. We apply the present method to the bindings of small fragment molecules [4-hydroxy-2-butanone (BUT) and methyl methylthiomethyl sulphoxide (DSS)] to FK506 binding protein (FKBP) in an aqueous solution. Estimated binding rate constants are consistent with those obtained from long-timescale MD simulations. Furthermore, by decomposing the rate constants to the thermodynamic and kinetic contributions, we clarify that the higher thermodynamic stability of the reactive state for DSS causes the faster binding kinetics compared with BUT.

Non-equilibrium modalities of inhibition: Characterizing irreversible inhibition for the ErbB receptor family members

Most drug target interactions for clinically approved small-molecules are non-equilibrium slow-onset, tight-binding or irreversible in nature, with pronounced element of time-dependence of inhibition. Analysis of such modality of inhibition requires a continuous enzyme kinetic measurement that can yield complete progress curves and an automated high-throughput analysis pipeline. Given the increasing emphasis on designing non-equilibrium modes of inhibiting an enzyme target (especially irreversible), the above specified pipeline for data generation and analysis is essential for extracting parameters to guide decisions in early drug discovery. In this manuscript, the methodology and data analysis protocol from our irreversible inhibitor characterization campaigns for the ErbB receptor family members is presented. Guidance is provided on the appropriate design of assay to generate quality data, setting up the analysis and estimation of inactivation rate (kinact) and the pseudo-equilibrium binding affinity (KI) constant (or their ratio kinact/KI) in a high-throughput manner for the inhibitor interacting with the protein target of interest.

Structure-Kinetic Relationship for Drug Design Revealed by a PLS Model with Retrosynthesis-Based Pre-Trained Molecular Representation and Molecular Dynamics Simulation

Drug design based on kinetic properties is growing in application. Here, we applied retrosynthesis-based pre-trained molecular representation (RPM) in machine learning (ML) to train 501 inhibitors of 55 proteins and successfully predicted the dissociation rate constant (k_off) values of 38 inhibitors from an independent dataset for the N-terminal domain of heat shock protein 90α (N-HSP90). Our RPM molecular representation outperforms other pre-trained molecular representations such as GEM, MPG, and general molecular descriptors from RDKit. Furthermore, we optimized the accelerated molecular dynamics to calculate the relative retention time (RT) for the 128 inhibitors of N-HSP90 and obtained the protein-ligand interaction fingerprints (IFPs) on their dissociation pathways and their influencing weights on the k_off value. We observed a high correlation among the simulated, predicted, and experimental -log(k_off) values. Combining ML, molecular dynamics (MD) simulation, and IFPs derived from accelerated MD helps design a drug for specific kinetic properties and selectivity profiles to the target of interest. To further validate our k_off predictive ML model, we tested our model on two new N-HSP90 inhibitors, which have experimental k_off values and are not in our ML training dataset. The predicted k_off values are consistent with experimental data, and the mechanism of their kinetic properties can be explained by IFPs, which shed light on the nature of their selectivity against N-HSP90 protein. We believe that the ML model described here is transferable to predict k_off of other proteins and will enhance the kinetics-based drug design endeavor.

Predicting Protein-Ligand Binding and Unbinding Kinetics with Biased MD Simulations and Coarse-Graining of Dynamics: Current State and Challenges

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.

Structure-affinity and structure-residence time relationships of macrocyclic Gαq protein inhibitors

The macrocyclic depsipeptides YM-254890 (YM) and FR900359 (FR) are potent inhibitors of Gα_q/11 proteins. They are important pharmacological tools and have potential as therapeutic drugs. The hydrogenated, tritium-labeled YM and FR derivatives display largely different residence times despite similar structures. In the present study we established a competition-association binding assay to determine the dissociation kinetics of unlabeled G_q protein inhibitors. Structure-affinity and structure-residence time relationships were analyzed. Small structural modifications had a large impact on residence time. YM and FR exhibited 4- to 10-fold higher residence times than their hydrogenated derivatives. While FR showed pseudo-irreversible binding, YM displayed much faster dissociation from its target. The isopropyl anchor present in FR and some derivatives was essential for slow dissociation. These data provide a basis for future drug design toward modulating residence times of macrocyclic G_q protein inhibitors, which has been recognized as a crucial determinant for therapeutic outcome.

Predicting Biomolecular Binding Kinetics: A Review

Biomolecular binding kinetics including the association (k_on) and dissociation (k_off) rates are critical parameters for therapeutic design of small-molecule drugs, peptides, and antibodies. Notably, the drug molecule residence time or dissociation rate has been shown to correlate with their efficacies better than binding affinities. A wide range of modeling approaches including quantitative structure-kinetic relationship models, Molecular Dynamics simulations, enhanced sampling, and Machine Learning has been developed to explore biomolecular binding and dissociation mechanisms and predict binding kinetic rates. Here, we review recent advances in computational modeling of biomolecular binding kinetics, with an outlook for future improvements.

NanoB2 to monitor interactions of ligands with membrane proteins by combining nanobodies and NanoBRET

The therapeutic potential of ligands targeting disease-associated membrane proteins is predicted by ligand-receptor binding constants, which can be determined using NanoLuciferase (NanoLuc)-based bioluminescence resonance energy transfer (NanoBRET) methods. However, the broad applicability of these methods is hampered by the restricted availability of fluorescent probes. We describe the use of antibody fragments, like nanobodies, as universal building blocks for fluorescent probes for use in NanoBRET. Our nanobody-NanoBRET (NanoB2) workflow starts with the generation of NanoLuc-tagged receptors and fluorescent nanobodies, enabling homogeneous, real-time monitoring of nanobody-receptor binding. Moreover, NanoB2 facilitates the assessment of receptor binding of unlabeled ligands in competition binding experiments. The broad significance is illustrated by the successful application of NanoB2 to different drug targets (e.g., multiple G protein-coupled receptors [GPCRs] and a receptor tyrosine kinase [RTK]) at distinct therapeutically relevant binding sites (i.e., extracellular and intracellular).

Ligand Gaussian Accelerated Molecular Dynamics 2 (LiGaMD2): Improved Calculations of Ligand Binding Thermodynamics and Kinetics with Closed Protein Pocket

Ligand binding thermodynamics and kinetics are critical parameters for drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics from molecular simulations due to limited simulation timescales. Protein dynamics, especially in the ligand binding pocket, often plays an important role in ligand binding. Based on our previously developed Ligand Gaussian accelerated molecular dynamics (LiGaMD), here we present LiGaMD2 in which a selective boost potential was applied to both the ligand and protein residues in the binding pocket to improve sampling of ligand binding and dissociation. To validate the performance of LiGaMD2, the T4 lysozyme (T4L) mutants with open and closed pockets bound by different ligands were chosen as model systems. LiGaMD2 could efficiently capture repetitive ligand dissociation and binding within microsecond simulations of all T4L systems. The obtained ligand binding kinetic rates and free energies agreed well with available experimental values and previous modeling results. Therefore, LiGaMD2 provides an improved approach to sample opening of closed protein pockets for ligand dissociation and binding, thereby allowing for efficient calculations of ligand binding thermodynamics and kinetics.

Transferring chemical and energetic knowledge between molecular systems with machine learning

Predicting structural and energetic properties of a molecular system is one of the fundamental tasks in molecular simulations, and it has applications in chemistry, biology, and medicine. In the past decade, the advent of machine learning algorithms had an impact on molecular simulations for various tasks, including property prediction of atomistic systems. In this paper, we propose a novel methodology for transferring knowledge obtained from simple molecular systems to a more complex one, endowed with a significantly larger number of atoms and degrees of freedom. In particular, we focus on the classification of high and low free-energy conformations. Our approach relies on utilizing (i) a novel hypergraph representation of molecules, encoding all relevant information for characterizing multi-atom interactions for a given conformation, and (ii) novel message passing and pooling layers for processing and making free-energy predictions on such hypergraph-structured data. Despite the complexity of the problem, our results show a remarkable Area Under the Curve of 0.92 for transfer learning from tri-alanine to the deca-alanine system. Moreover, we show that the same transfer learning approach can also be used in an unsupervised way to group chemically related secondary structures of deca-alanine in clusters having similar free-energy values. Our study represents a proof of concept that reliable transfer learning models for molecular systems can be designed, paving the way to unexplored routes in prediction of structural and energetic properties of biologically relevant systems.

Advances in computational methods for ligand binding kinetics

Binding kinetic parameters can be correlated with drug efficacy, which in recent years led to the development of various computational methods for predicting binding kinetic rates and gaining insight into protein-drug binding paths and mechanisms. In this review, we introduce and compare computational methods recently developed and applied to two systems, trypsin-benzamidine and kinase-inhibitor complexes. Methods involving enhanced sampling in molecular dynamics simulations or machine learning can be used not only to predict kinetic rates, but also to reveal factors modulating the duration of residence times, selectivity, and drug resistance to mutations. Methods which require less computational time to make predictions are highlighted, and suggestions to reduce the error of computed kinetic rates are presented.

Heterobivalent Inhibitors of Acetyl-CoA Carboxylase: Drug Target Residence Time and Time-Dependent Antibacterial Activity

The relationship between drug-target residence time and the post-antibiotic effect (PAE) provides insights into target vulnerability. To probe the vulnerability of bacterial acetyl-CoA carboxylase (ACC), a series of heterobivalent inhibitors were synthesized based on pyridopyrimidine 1 and moiramide B (3) which bind to the biotin carboxylase and carboxyltransferase ACC active sites, respectively. The heterobivalent compound 17, which has a linker of 50 Å, was a tight binding inhibitor of Escherichia coli ACC (Kiapp 0.2 nM) and could be displaced from ACC by a combination of both 1 and 3 but not just by 1. In agreement with the prolonged occupancy of ACC resulting from forced proximity binding, the heterobivalent inhibitors produced a PAE in E. coli of 1-4 h in contrast to 1 and 3 in combination or alone, indicating that ACC is a vulnerable target and highlighting the utility of kinetic, time-dependent effects in the drug mechanism of action.

pH Regulates Ligand Binding to an Enzyme Active Site by Modulating Intermediate Populations

Understanding the mechanism of ligands binding to their protein targets and the influence of various factors governing the binding thermodynamics is essential for rational drug design. The solution pH is one of the critical factors that can influence ligand binding to a protein cavity, especially in enzymes whose function is sensitive to the pH. Using computer simulations, we studied the pH effect on the binding of a guanidinium ion (Gdm⁺) to the active site of hen egg-white lysozyme (HEWL). HEWL serves as a model system for enzymes with two acidic residues in the active site and ligands with Gdm⁺ moieties, which can bind to the active sites of such enzymes and are present in several approved drugs treating various disorders. The computed free energy surface (FES) shows that Gdm⁺ binds to the HEWL active site using two dominant binding pathways populating multiple intermediates. We show that the residues close to the active site that can anchor the ligand could play a critical role in ligand binding. Using a Markov state model, we quantified the lifetimes and kinetic pathways connecting the different states in the FES. The protonation and deprotonation of the acidic residues in the active site in response to the pH change strongly influence the Gdm⁺ binding. There is a sharp jump in the ligand-binding rate constant when the pH approaches the largest pK_a of the acidic residue present in the active site. The simulations reveal that, at most, three Gdm⁺ can bind at the active site, with the Gdm⁺ bound in the cavity of the active site acting as a scaffold for the other two Gdm⁺ ions binding. These results can aid in providing greater insights into designing novel molecules containing Gdm⁺ moieties that can have high binding affinities to inhibit the function of enzymes with acidic residues in their active site.

Qualitative Estimation of Protein-Ligand Complex Stability through Thermal Titration Molecular Dynamics Simulations

The prediction of ligand efficacy has long been linked to thermodynamic properties such as the equilibrium dissociation constant, which considers both the association and the dissociation rates of a defined protein-ligand complex. In the last 15 years, there has been a paradigm shift, with an increased interest in the determination of kinetic properties such as the drug-target residence time since they better correlate with ligand efficacy compared to other parameters. In this article, we present thermal titration molecular dynamics (TTMD), an alternative computational method that combines a series of molecular dynamics simulations performed at progressively increasing temperatures with a scoring function based on protein-ligand interaction fingerprints for the qualitative estimation of protein-ligand-binding stability. The protocol has been applied to four different pharmaceutically relevant test cases, including protein kinase CK1δ, protein kinase CK2, pyruvate dehydrogenase kinase 2, and SARS-CoV-2 main protease, on a variety of ligands with different sizes, structures, and experimentally determined affinity values. In all four cases, TTMD was successfully able to distinguish between high-affinity compounds (low nanomolar range) and low-affinity ones (micromolar), proving to be a useful screening tool for the prioritization of compounds in a drug discovery campaign.

Estimation of Drug-Target Residence Time by Targeted Molecular Dynamics Simulations

Drug-target residence time has emerged as a key selection factor in drug discovery since the binding duration of a drug molecule to its protein target can significantly impact its in vivo efficacy. The challenge in studying the residence time, in early drug discovery stages, lies in how to cost-effectively determine the residence time for the systematic assessment of compounds. Currently, there is still a lack of computational protocols to quickly estimate such a measure, particularly for large and flexible protein targets and drugs. Here, we report an efficient computational protocol, based on targeted molecular dynamics, to rank drug candidates by their residence time and to obtain insights into ligand-target dissociation mechanisms. The method was assessed on a dataset of 10 arylpyrazole inhibitors of CDK8, a large, flexible, and clinically important target, for which the experimental residence time of the inhibitors ranges from minutes to hours. The compounds were correctly ranked according to their estimated residence time scores compared to their experimental values. The analysis of protein-ligand interactions along the dissociation trajectories highlighted the favorable contribution of hydrophobic contacts to residence time and revealed key residues that strongly affect compound residence time.

Design, Synthesis, and Pharmacological Evaluation of Second-Generation Soluble Adenylyl Cyclase (sAC, ADCY10) Inhibitors with Slow Dissociation Rates

Soluble adenylyl cyclase (sAC: ADCY10) is an enzyme involved in intracellular signaling. Inhibition of sAC has potential therapeutic utility in a number of areas. For example, sAC is integral to successful male fertility: sAC activation is required for sperm motility and ability to undergo the acrosome reaction, two processes central to oocyte fertilization. Pharmacologic evaluation of existing sAC inhibitors for utility as on-demand, nonhormonal male contraceptives suggested that both high intrinsic potency, fast on and slow dissociation rates are essential design elements for successful male contraceptive applications. During the course of the medicinal chemistry campaign described here, we identified sAC inhibitors that fulfill these criteria and are suitable for in vivo evaluation of diverse sAC pharmacology.

Ligand Unbinding Pathway and Mechanism Analysis Assisted by Machine Learning and Graph Methods

We present two methods to reveal protein-ligand unbinding mechanisms in biased unbinding simulations by clustering trajectories into ensembles representing unbinding paths. The first approach is based on a contact principal component analysis for reducing the dimensionality of the input data, followed by identification of unbinding paths and training a machine learning model for trajectory clustering. The second approach clusters trajectories according to their pairwise mean Euclidean distance employing the neighbor-net algorithm, which takes into account input data bias in the distances set and is superior to dendrogram construction. Finally, we describe a more complex case where the reaction coordinate relevant for path identification is a single intraligand hydrogen bond, highlighting the challenges involved in unbinding path reaction coordinate detection.

Assessing potency and binding kinetics of soluble adenylyl cyclase (sAC) inhibitors to maximize therapeutic potential

In mammalian cells, 10 different adenylyl cyclases produce the ubiquitous second messenger, cyclic adenosine monophosphate (cAMP). Amongst these cAMP-generating enzymes, bicarbonate (HCO3 -)-regulated soluble adenylyl cyclase (sAC; ADCY10) is uniquely essential in sperm for reproduction. For this reason, sAC has been proposed as a potential therapeutic target for non-hormonal contraceptives for men. Here, we describe key sAC-focused in vitro assays to identify and characterize sAC inhibitors for therapeutic use. The affinity and binding kinetics of an inhibitor can greatly influence in vivo efficacy, therefore, we developed improved assays for assessing these efficacy defining features.

Insight into Thermodynamic and Kinetic Profiles in Small-Molecule Optimization

Structure-activity relationships (SARs) and structure-property relationships (SPRs) have been considered the most important factors during the drug optimization process. For medicinal chemists, improvements in the potencies and druglike properties of small molecules are regarded as their major goals. Among them, the binding affinity and selectivity of small molecules on their targets are the most important indicators. In recent years, there has been growing interest in using thermodynamic and kinetic profiles to analyze ligand-receptor interactions, which could provide not only binding affinities but also detailed binding parameters for small-molecule optimization. In this perspective, we are trying to provide an insight into thermodynamic and kinetic profiles in small-molecule optimization. Through a highlight of strategies on the small-molecule optimization with specific cases, we aim to put forward the importance of structure-thermodynamic relationships (STRs) and structure-kinetic relationships (SKRs), which could provide more guidance to find safe and effective small-molecule drugs.

SEEKR2: Versatile Multiscale Milestoning Utilizing the OpenMM Molecular Dynamics Engine

We present SEEKR2 (simulation-enabled estimation of kinetic rates version 2)—the latest iteration in the family of SEEKR programs for using multiscale simulation methods to computationally estimate the kinetics and thermodynamics of molecular processes, in particular, ligand-receptor binding. SEEKR2 generates equivalent, or improved, results compared to the earlier versions of SEEKR but with significant increases in speed and capabilities. SEEKR2 has also been built with greater ease of usability and with extensible features to enable future expansions of the method. Now, in addition to supporting simulations using NAMD, calculations may be run with the fast and extensible OpenMM simulation engine. The Brownian dynamics portion of the calculation has also been upgraded to Browndye 2. Furthermore, this version of SEEKR supports hydrogen mass repartitioning, which significantly reduces computational cost, while showing little, if any, loss of accuracy in the predicted kinetics.

Fast approximative methods for study of ligand transport and rational design of improved enzymes for biotechnologies

Acceleration of chemical reactions by the enzymes optimized using protein engineering represents one of the key pillars of the contribution of biotechnology towards sustainability. Tunnels and channels of enzymes with buried active sites enable the exchange of ligands, ions, and water molecules between the outer environment and active site pockets. The efficient exchange of ligands is a fundamental process of biocatalysis. Therefore, enzymes have evolved a wide range of mechanisms for repetitive conformational changes that enable periodic opening and closing. Protein-ligand interactions are traditionally studied by molecular docking, whereas molecular dynamics is the method of choice for studying conformational changes and ligand transport. However, computational demands make molecular dynamics impractical for screening purposes. Thus, several approximative methods have been recently developed to study interactions between a protein and ligand during the ligand transport process. Apart from identifying the best binding modes, these methods also provide information on the energetics of the transport and identify problematic regions limiting the ligand passage. These methods use approximations to simulate binding or unbinding events rapidly (calculation times from minutes to hours) and provide energy profiles that can be used to rank ligands or pathways. Here we provide a critical comparison of available methods, showcase their results on sample systems, discuss their practical applications in molecular biotechnologies and outline possible future developments.

Public Data Set of Protein-Ligand Dissociation Kinetic Constants for Quantitative Structure-Kinetics Relationship Studies

Protein-ligand binding affinity reflects the equilibrium thermodynamics of the protein-ligand binding process. Binding/unbinding kinetics is the other side of the coin. Computational models for interpreting the quantitative structure-kinetics relationship (QSKR) aim at predicting protein-ligand binding/unbinding kinetics based on protein structure, ligand structure, or their complex structure, which in principle can provide a more rational basis for structure-based drug design. Thus far, most of the public data sets used for deriving such QSKR models are rather limited in sample size and structural diversity. To tackle this problem, we have compiled a set of 680 protein-ligand complexes with experimental dissociation rate constants (k off), which were mainly curated from the references accumulated for updating our PDBbind database. Three-dimensional structure of each protein-ligand complex in this data set was either retrieved from the Protein Data Bank or carefully modeled based on a proper template. The entire data set covers 155 types of protein, with their dissociation kinetic constants (k off) spanning nearly 10 orders of magnitude. To the best of our knowledge, this data set is the largest of its kind reported publicly. Utilizing this data set, we derived a random forest (RF) model based on protein-ligand atom pair descriptors for predicting k off values. We also demonstrated that utilizing modeled structures as additional training samples will benefit the model performance. The RF model with mixed structures can serve as a baseline for testifying other more sophisticated QSKR models. The whole data set, namely, PDBbind-koff-2020, is available for free download at our PDBbind-CN web site (http://www.pdbbind.org.cn/download.php).

Design, Synthesis, and Biological Characterization of Inhaled p38α/β MAPK Inhibitors for the Treatment of Lung Inflammatory Diseases

The identification of novel inhaled p38α/β mitogen-activated protein kinases (MAPK) (MAPK14/11) inhibitors suitable for the treatment of pulmonary inflammatory conditions has been described. A rational drug design approach started from the identification of a novel tetrahydronaphthalene series, characterized by nanomolar inhibition of p38α with selectivity over p38γ and p38δ isoforms. SAR optimization of 1c is outlined, where improvements in potency against p38α and ligand-enzyme dissociation kinetics led to several compounds showing pronounced anti-inflammatory effects in vitro (inhibition of TNFα release). Targeting of the defined physicochemical properties allowed the identification of compounds 3h, 4e, and 4f, which showed, upon intratracheal instillation, low plasma levels, prolonged lung retention, and anti-inflammatory effects in a rat acute model of a bacterial endotoxin-induced pulmonary inflammation. Compound 4e, in particular, displayed remarkable efficacy and duration of action and was selected for progression in disease models of asthma and chronic obstructive pulmonary disease (COPD).

Determination of Multidirectional Pathways for Ligand Release from the Receptor: A New Approach Based on Differential Evolution

Steered molecular dynamics (SMD) simulation is a powerful method in computer-aided drug design as it can be used to access the relative binding affinity with high precision but with low computational cost. The success of SMD depends on the choice of the direction along which the ligand is pulled from the receptor-binding site. In most simulations, the unidirectional pathway was used, but in some cases, this choice resulted in the ligand colliding with the complex surface of the exit tunnel. To overcome this difficulty, several variants of SMD with multidirectional pulling have been proposed, but they are not completely devoid of disadvantages. Here, we have proposed to determine the direction of pulling with a simple scoring function that minimizes the receptor–ligand interaction, and an optimization algorithm called differential evolution is used for energy minimization. The effectiveness of our protocol was demonstrated by finding expulsion pathways of Huperzine A and camphor from the binding site of Torpedo California acetylcholinesterase and P450cam proteins, respectively, and comparing them with the previous results obtained using memetic sampling and random acceleration molecular dynamics. In addition, by applying this protocol to a set of ligands bound with LSD1 (lysine specific demethylase 1), we obtained a much higher correlation between the work of pulling force and experimental data on the inhibition constant IC50 compared to that obtained using the unidirectional approach based on minimal steric hindrance.

Kinetic Modeling of Time-Dependent Enzyme Inhibition by Pre-Steady-State Analysis of Progress Curves: The Case Study of the Anti-Alzheimer's Drug Galantamine

The Michaelis–Menten model of enzyme kinetic assumes the free ligand approximation, the steady-state approximation and the rapid equilibrium approximation. Analytical methods to model slow-binding inhibitors by the analysis of initial velocities have been developed but, due to their inherent complexity, they are seldom employed. In order to circumvent the complications that arise from the violation of the rapid equilibrium assumption, inhibition is commonly evaluated by pre-incubating the enzyme and the inhibitors so that, even for slow inhibitors, the binding equilibrium is established before the reaction is started. Here, we show that for long drug-target residence time inhibitors, the conventional analysis of initial velocities by the linear regression of double-reciprocal plots fails to provide a correct description of the inhibition mechanism. As a case study, the inhibition of acetylcholinesterase by galantamine, a drug approved for the symptomatic treatment of Alzheimer’s disease, is reported. For over 50 years, analysis based on the conventional steady-state model has overlooked the time-dependent nature of galantamine inhibition, leading to an erroneous assessment of the drug potency and, hence, to discrepancies between biochemical data and the pharmacological evidence. Re-examination of acetylcholinesterase inhibition by pre-steady state analysis of the reaction progress curves showed that the potency of galantamine has indeed been underestimated by a factor of ~100.

RNA-ligand molecular docking: advances and challenges

With rapid advances in computer algorithms and hardware, fast and accurate virtual screening has led to a drastic acceleration in selecting potent small molecules as drug candidates. Computational modeling of RNA-small molecule interactions has become an indispensable tool for RNA-targeted drug discovery. The current models for RNA-ligand binding have mainly focused on the docking-and-scoring method. Accurate docking and scoring should tackle four crucial problems: (1) conformational flexibility of ligand, (2) conformational flexibility of RNA, (3) efficient sampling of binding sites and binding poses, and (4) accurate scoring of different binding modes. Moreover, compared with the problem of protein-ligand docking, predicting ligand binding to RNA, a negatively charged polymer, is further complicated by additional effects such as metal ion effects. Thermodynamic models based on physics-based and knowledge-based scoring functions have shown highly encouraging success in predicting ligand binding poses and binding affinities. Recently, kinetic models for ligand binding have further suggested that including dissociation kinetics (residence time) in ligand docking would result in improved performance in estimating in vivo drug efficacy. More recently, the rise of deep-learning approaches has led to new tools for predicting RNA-small molecule binding. In this review, we present an overview of the recently developed computational methods for RNA-ligand docking and their advantages and disadvantages.

Practical Protocols for Efficient Sampling of Kinase-Inhibitor Binding Pathways Using Two-Dimensional Replica-Exchange Molecular Dynamics

Molecular dynamics (MD) simulations are increasingly used to study various biological processes such as protein folding, conformational changes, and ligand binding. These processes generally involve slow dynamics that occur on the millisecond or longer timescale, which are difficult to simulate by conventional atomistic MD. Recently, we applied a two-dimensional (2D) replica-exchange MD (REMD) method, which combines the generalized replica exchange with solute tempering (gREST) with the replica-exchange umbrella sampling (REUS) in kinase-inhibitor binding simulations, and successfully observed multiple ligand binding/unbinding events. To efficiently apply the gREST/REUS method to other kinase-inhibitor systems, we establish modified, practical protocols with non-trivial simulation parameter tuning. The current gREST/REUS simulation protocols are tested for three kinase-inhibitor systems: c-Src kinase with PP1, c-Src kinase with Dasatinib, and c-Abl kinase with Imatinib. We optimized the definition of kinase-ligand distance as a collective variable (CV), the solute temperatures in gREST, and replica distributions and umbrella forces in the REUS simulations. Also, the initial structures of each replica in the 2D replica space were prepared carefully by pulling each ligand from and toward the protein binding sites for keeping stable kinase conformations. These optimizations were carried out individually in multiple short MD simulations. The current gREST/REUS simulation protocol ensures good random walks in 2D replica spaces, which are required for enhanced sampling of inhibitor dynamics around a target kinase.

A Scintillation Proximity Assay for Real-Time Kinetic Analysis of Chemokine–Chemokine Receptor Interactions

Chemokine receptors are extensively involved in a broad range of physiological and pathological processes, making them attractive drug targets. However, despite considerable efforts, there are very few approved drugs targeting this class of seven transmembrane domain receptors to date. In recent years, the importance of including binding kinetics in drug discovery campaigns was emphasized. Therefore, kinetic insight into chemokine–chemokine receptor interactions could help to address this issue. Moreover, it could additionally deepen our understanding of the selectivity and promiscuity of the chemokine–chemokine receptor network. Here, we describe the application, optimization and validation of a homogenous Scintillation Proximity Assay (SPA) for real-time kinetic profiling of chemokine–chemokine receptor interactions on the example of ACKR3 and CXCL12. The principle of the SPA is the detection of radioligand binding to receptors reconstituted into nanodiscs by scintillation light. No receptor modifications are required. The nanodiscs provide a native-like environment for receptors and allow for full control over bilayer composition and size. The continuous assay format enables the monitoring of binding reactions in real-time, and directly accounts for non-specific binding and potential artefacts. Minor adaptations additionally facilitate the determination of equilibrium binding metrics, making the assay a versatile tool for the study of receptor–ligand interactions.

Estimation of binding rates and affinities from multiensemble Markov models and ligand decoupling

Accurate and efficient simulation of the thermodynamics and kinetics of protein-ligand interactions is crucial for computational drug discovery. Multiensemble Markov Model (MEMM) estimators can provide estimates of both binding rates and affinities from collections of short trajectories but have not been systematically explored for situations when a ligand is decoupled through scaling of non-bonded interactions. In this work, we compare the performance of two MEMM approaches for estimating ligand binding affinities and rates: (1) the transition-based reweighting analysis method (TRAM) and (2) a Maximum Caliber (MaxCal) based method. As a test system, we construct a small host-guest system where the ligand is a single uncharged Lennard-Jones (LJ) particle, and the receptor is an 11-particle icosahedral pocket made from the same atom type. To realistically mimic a protein-ligand binding system, the LJ ϵ parameter was tuned, and the system was placed in a periodic box with 860 TIP3P water molecules. A benchmark was performed using over 80 µs of unbiased simulation, and an 18-state Markov state model was used to estimate reference binding affinities and rates. We then tested the performance of TRAM and MaxCal when challenged with limited data. Both TRAM and MaxCal approaches perform better than conventional Markov state models, with TRAM showing better convergence and accuracy. We find that subsampling of trajectories to remove time correlation improves the accuracy of both TRAM and MaxCal and that in most cases, only a single biased ensemble to enhance sampled transitions is required to make accurate estimates.

Molecular Simulations and Drug Discovery of Adenosine Receptors

G protein-coupled receptors (GPCRs) represent the largest family of human membrane proteins. Four subtypes of adenosine receptors (ARs), the A1AR, A2AAR, A2BAR and A3AR, each with a unique pharmacological profile and distribution within the tissues in the human body, mediate many physiological functions and serve as critical drug targets for treating numerous human diseases including cancer, neuropathic pain, cardiac ischemia, stroke and diabetes. The A1AR and A3AR preferentially couple to the Gi/o proteins, while the A2AAR and A2BAR prefer coupling to the Gs proteins. Adenosine receptors were the first subclass of GPCRs that had experimental structures determined in complex with distinct G proteins. Here, we will review recent studies in molecular simulations and computer-aided drug discovery of the adenosine receptors and also highlight their future research opportunities.

Ligands with different dimeric configurations potently activate the EphA2 receptor and reveal its potential for biased signaling

The EphA2 receptor tyrosine kinase activates signaling pathways with different, and sometimes opposite, effects in cancer and other pathologies. Thus, highly specific and potent biased ligands that differentially control EphA2 signaling responses could be therapeutically valuable. Here, we use EphA2-specific monomeric peptides to engineer dimeric ligands with three different geometric configurations to combine a potential ability to differentially modulate EphA2 signaling responses with the high potency and prolonged receptor residence time characteristic of dimeric ligands. The different dimeric peptides readily induce EphA2 clustering, autophosphorylation and signaling, the best with sub-nanomolar potency. Yet, there are differences in two EphA2 signaling responses induced by peptides with different configurations, which exhibit distinct potency and efficacy. The peptides bias signaling when compared with the ephrinA1-Fc ligand and do so via different mechanisms. These findings provide insights into Eph receptor signaling, and proof-of-principle that different Eph signaling responses can be distinctly modulated.

Implications of ligand-receptor binding kinetics on GLP-1R signalling

G protein-coupled receptors (GPCRs) are the largest class of membrane proteins and in recent years there has been a growing appreciation of the importance in understanding temporal aspects of GPCR behaviour, including the kinetics of ligand binding and downstream receptor mediated signalling. Class B1 GPCRs are activated by peptide agonists and are validated therapeutic targets for numerous diseases. However, the kinetics of ligand binding and how this is linked to downstream activation of signalling cascades is not routinely assessed in development of peptide agonists for this receptor class. The glucagon-like peptide-1 receptor (GLP-1R) is a prototypical class B1 GPCR and a validated target for treatment of global health burdens, including type 2 diabetes and obesity. In this study we examined the kinetics of different steps in GLP-1R activation and subsequent cAMP production mediated by a series of GLP-1R peptide agonists, including the ligand-receptor interaction, ligand-receptor-mediated G protein engagement and conformational change and cAMP production. Our results revealed GLP-1R peptide agonist dissociation kinetics (K_off), but not association kinetics (K_on), were positively correlated with the onset of receptor-G protein coupling/conformational change, onset of cAMP production and duration of cAMP signalling. Thus, this study advances the understanding of molecular events that couple GLP-1R ligand binding to intracellular signaling, with the findings likely to have implications for mechanistic understanding of agonist action at other related class B1 GPCRs.

In the footsteps of an inhibitor unbinding from the active site of human carbonic anhydrase II

The crystal structure of human carbonic anhydrase (HCA) II bound to an inhibitor molecule, 6-hydroxy-2-thioxocoumarin (FC5), shows FC5 to be located in a hydrophobic pocket at the active site. The present work employs classical molecular dynamics (MD) simulation to follow the FC5 molecule for 1 μs as it unbinds from its binding location, adopts the path of substrate/product diffusion (path 1) to leave the active site at around 75 ns. It is then found to undergo repeated binding and unbinding at different locations on the surface of the enzyme in water. Several transient excursions through different regions of the enzyme are also observed prior to its exit from the active site. These transient paths are combined with functionally relevant cavities/channels to enlist five additional pathways (path 2-6). Pathways 1-6 are subsequently explored using steered MD and umbrella sampling simulations. A free energy barrier of 0.969 kcal mol^-1 is encountered along path 1, while barriers in the range of 0.57-2.84 kcal mol^-1 are obtained along paths 2, 4 and 5. We also analyze in detail the interaction between FC5 and the enzyme along each path as the former leaves the active site of HCA II. Our results indicate path 1 to be the major exit pathway for FC5, although competing contributions may also come from the paths 2, 4 and 5.Communicated by Ramaswamy H. Sarma.

A guide to enzyme kinetics in early drug discovery

Drugs interact with their target of interest to bring about the desired phenotypic outcome that results in disease alleviation. Traditionally, most lead optimization exercises were driven by affinity measures (like IC₅₀ ) to inform structure-activity relationship (SAR)-guided medicinal chemistry. However, an IC₅₀ value is a thermodynamic estimate measured under equilibrium conditions that can vary as a function of substrate concentration and/or time (the latter especially for nonequilibrium modalities). Further, like other thermodynamic estimates, it is a state-function that is indifferent to the path traversed from the initial state to the final state. This can be a cause for concern in drug discovery given the predominance of nonequilibrium interactions and the open thermodynamic nature of the human system. Under such situations, employing rates along with equilibrium constants (or IC₅₀ values) would be far more relevant to capture the time evolution of the small molecule's interaction with the target of interest. These rates are generally typified by the rate of association, rate of dissociation and the residence time of the small molecule on the target (target occupancy). These parameters, when combined with the concept of target vulnerability, therapeutic window, pharmacokinetic profile of the small molecule, estimates of endogenous ligand and target turnover, will shed critical insights into the kinetics and dynamics of a small molecule's interaction with the protein, and allow realistic modelling of the system to enable optimizations and dosing decisions. With that aim, this guide will attempt to introduce the traditional role of mechanistic enzymology within drug discovery and emphasize the importance of kinetics in guiding SAR-based optimizations. It will also present initial ideas on how kinetic investigation should be positioned relative to the temporal span of a drug-discovery pipeline to leverage maximal utility from the investment in time and effort.

Computing Rates and Understanding Unbinding Mechanisms in Host-Guest Systems

The long time scale associated with ligand residence times renders their computation challenging. Therefore, the influence of factors like solvation and steric hindrance on residence times is not fully understood. Here, we demonstrate in a set of model host-guest systems that the recently developed Gaussian mixture based enhanced sampling allows residence times to be computed and enables an understanding of their unbinding mechanism. We observe that guest unbinding often proceeds via a series of intermediate states that can be labeled by the number of water molecules present in the binding cavity. In several cases the residence time is correlated to the water trapping times in the cavity.

Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased “downhill” trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand–protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

In Vitro and In Silico Analysis of the Residence Time of Serotonin 5-HT7 Receptor Ligands with Arylpiperazine Structure: A Structure-Kinetics Relationship Study

During the last decade, the kinetics of drug-target interaction has received increasing attention as an important pharmacological parameter in the drug development process. Several studies have suggested that the lipophilicity of a molecule can play an important role. To date, this aspect has been studied for several G protein-coupled receptors (GPCRs) ligands but not for the 5-HT₇ receptor (5-HT₇R), a GPCR proposed as a valid therapeutic target in neurodevelopmental and neuropsychiatric disorders associated with abnormal neuronal connectivity. In this study, we report on structure-kinetics relationships of a set of arylpiperazine-based 5-HT₇R ligands. We found that it is not the overall lipophilicity of the molecule that influences drug-target interaction kinetics but rather the position of polar groups within the molecule. Next, we performed a combination of molecular docking studies and molecular dynamics simulations to gain insights into structure-kinetics relationships. These studies did not suggest specific contact patterns between the ligands and the receptor-binding site as determinants for compounds kinetics. Finally, we compared the abilities of two 5-HT₇R agonists with similar receptor-binding affinities and different residence times to stimulate the 5-HT₇R-mediated neurite outgrowth in mouse neuronal primary cultures and found that the compounds induced the effect with different timing. This study provides the first insights into the binding kinetics of arylpiperazine-based 5-HT₇R ligands that can be helpful to design new 5-HT₇R ligands with fine-tuning of the kinetic profile.

Integrating Incompatible Assay Data Sets with Deep Preference Learning

A large amount of bioactivity assay data is already accumulated in public databases, but the integration of these data sets for quantitative structure-activity relationship (QSAR) studies is not straightforward due to differences in experimental methods and settings. We present an efficient deep-learning-based approach called Deep Preference Data Integration (DPDI). For integrating outcome variables of different assay types, a surrogate variable is introduced, and a neural network is trained such that the total order induced by the surrogate variable is maximally consistent with given data sets. In a task of predicting efficacy of factor Xa inhibitors, DPDI successfully integrated 2959 molecules distributed in 129 assay data sets. In most of our experiments, data integration improved prediction accuracy strongly in interpolation and extrapolation tasks, indicating that DPDI is an effective tool for QSAR studies.

Predicting Residence Time of GPCR Ligands with Machine Learning

Drug-target residence time, the duration of binding at a given protein target, has been shown in some protein families to be more significant for conferring efficacy than binding affinity. To carry out efficient optimization of residence time in drug discovery, machine learning models that can predict that value need to be developed. One of the main challenges with predicting residence time is the paucity of data. This chapter outlines all of the currently available ligand kinetic data, providing a repository that contains the largest publicly available source of GPCR-ligand kinetic data to date. To help decipher the features of kinetic data that might be beneficial to include in computational models for the prediction of residence time, the experimental evidence for properties that influence residence time are summarized. Finally, two different workflows for predicting residence time with machine learning are outlined. The first is a single-target model trained on ligand features; the second is a multi-target model trained on features generated from molecular dynamics simulations.

Current Approaches of Building Mechanistic Pharmacodynamic Drug-Target Binding Models

Mechanistic pharmacodynamic models that incorporate the binding kinetics of drug-target interactions have several advantages in understanding target engagement and the efficacy of a drug dose. However, guidelines on how to build and interpret mechanistic pharmacodynamic drug-target binding models considering both biological and computational factors are still missing in the literature. In this chapter, current approaches of building mechanistic PD models and their advantages are discussed. We also present a methodology on how to select a suitable model considering both biological and computational perspectives, as well as summarize the challenges of current mechanistic PD models.

Challenges and frontiers of computational modelling of biomolecular recognition

Biomolecular recognition including binding of small molecules, peptides and proteins to their target receptors plays a key role in cellular function and has been targeted for therapeutic drug design. However, the high flexibility of biomolecules and slow binding and dissociation processes have presented challenges for computational modelling. Here, we review the challenges and computational approaches developed to characterise biomolecular binding, including molecular docking, molecular dynamics simulations (especially enhanced sampling) and machine learning. Further improvements are still needed in order to accurately and efficiently characterise binding structures, mechanisms, thermodynamics and kinetics of biomolecules in the future.

Using mathematical modeling to estimate time-independent cancer chemotherapy efficacy parameters

One of the primary cancer treatment modalities is chemotherapy. Unfortunately, traditional anti-cancer drugs are often not selective and cause damage to healthy cells, leading to serious side effects for patients. For this reason more targeted therapeutics and drug delivery methods are being developed. The effectiveness of new treatments is initially determined via in vitro cell viability assays, which determine the IC 50  of the drug. However, these assays are known to result in estimates of IC 50  that depend on the measurement time, possibly resulting in over- or under-estimation of the IC 50 . Here, we test the possibility of using cell growth curves and fitting of mathematical models to determine the IC 50  as well as the maximum efficacy of a drug ( ε max ). We measured cell growth of MCF-7 and HeLa cells in the presence of different concentrations of doxorubicin and fit the data with a logistic growth model that incorporates the effect of the drug. This method leads to measurement time-independent estimates of IC 50  and ε max , but we find that ε max  is not identifiable. Further refinement of this methodology is needed to produce uniquely identifiable parameter estimates.

Enhanced Sampling Approach to the Induced-Fit Docking Problem in Protein-Ligand Binding: The Case of Mono-ADP-Ribosylation Hydrolase Inhibitors

Enhanced sampling methods can predict free-energy landscapes associated with protein/ligand binding, characterizing the involved intermolecular interactions in a precise way. However, these in silico approaches can be challenged by induced-fit effects. Here, we present a variant of volume-based metadynamics tailored to tackle this problem in a general and efficient way. The validity of the approach is established by applying it to substrate/enzyme complexes of pharmacological relevance: mono-ADP-ribose (ADPr) in complex with mono-ADP-ribosylation hydrolases (MacroD1 and MacroD2), where induced-fit phenomena are known to be significant. The calculated binding free energies are consistent with experiments, with an absolute error smaller than 0.5 kcal/mol. Our simulations reveal that in all circumstances, the active loops, delimiting the boundaries of the binding site, undergo significant conformation rearrangements upon ligand binding. The calculations further provide, for the first time, the molecular basis of ADPr specificity and the relative changes in its experimental binding affinity on passing from MacroD1 to MacroD2 and all its mutants. Our study paves the way to the quantitative description of induced-fit events in molecular recognition.

Characterisation of bis(4-aminoquinoline)s as α1A adrenoceptor allosteric modulators

The development of sub-type selective α₁ adrenoceptor ligands has been hampered by the high sequence similarity of the amino acids forming the orthosteric binding pocket of the three α₁ adrenoceptor subtypes, along with other biogenic amine receptors. One possible approach to overcome this issue is to target allosteric sites on the α₁ adrenoceptors. Previous docking studies suggested that one of the quinoline moieties of a bis(4-aminoquinoline), comprising a 9-carbon methylene linker attached via the amine groups, could interact with residues outside of the orthosteric binding site while, simultaneously, the other quinoline moiety bound within the orthosteric site. We therefore hypothesized that this compound could act in a bitopic manner, displaying both orthosteric and allosteric binding properties. To test this proposition, we investigated the allosteric activity of a series of bis(4-aminoquinoline)s with linker lengths ranging from 2 to 12 methylene units (designated C2-C12). A linear trend of increasing [³H]prazosin dissociation rate with increasing linker length between C7 and C11 was observed, confirming their action as allosteric modulators. These data suggest that the optimal linker length for the bis(4-aminoquinoline)s to occupy the allosteric site of the α_1A adrenoceptor is between 7 and 11 methylene units. In addition, the ability of C9 bis(4-aminoquinoline) to modulate the activation of the α_1A adrenoceptor by norepinephrine was subsequently examined, showing that C9 acted as a non-competitive antagonist. Our findings indicate that the bis(4-aminoquinolines) are acting as allosteric modulators of orthosteric ligand binding, but not efficacy, in a bitopic manner.

Atomistic description of molecular binding processes based on returning probability theory

The efficiency of molecular binding such as host-guest binding is commonly evaluated in terms of kinetics, such as rate coefficients. In general, to compute the coefficient of the overall binding process, we need to consider both the diffusion of reactants and barrier crossing to reach the bound state. Here, we develop a methodology of quantifying the rate coefficient of binding based on molecular dynamics simulation and returning probability (RP) theory proposed by Kim and Lee [J. Chem. Phys. 131, 014503 (2009)]. RP theory provides a tractable formula of the rate coefficient in terms of the thermodynamic stability and kinetics of the intermediate state on a predefined reaction coordinate. In this study, the interaction energy between reactants is utilized as the reaction coordinate, enabling us to effectively describe the reactants' relative position and orientation on one-dimensional space. Application of this method to the host-guest binding systems, which consist of β-cyclodextrin and small guest molecules, yields the rate coefficients consistent with the experimental results.

Multiscale Approach for Computing Gated Ligand Binding from Molecular Dynamics and Brownian Dynamics Simulations

We develop an approach to characterize the effects of gating by a multiconformation protein consisting of macrostate conformations that are either accessible or inaccessible to ligand binding. We first construct a Markov state model of the apo-protein from atomistic molecular dynamics simulations from which we identify macrostates and their conformations, compute their relative macrostate populations and interchange kinetics, and structurally characterize them in terms of ligand accessibility. We insert the calculated first-order rate constants for conformational transitions into a multistate gating theory from which we derive a gating factor γ that quantifies the degree of conformational gating. Applied to HIV-1 protease, our approach yields a kinetic network of three accessible (semi-open, open, and wide-open) and two inaccessible (closed and a newly identified, "parted") macrostate conformations. The parted conformation sterically partitions the active site, suggesting a possible role in product release. We find that the binding kinetics of drugs and drug-like inhibitors to HIV-1 protease falls in the slow gating regime. However, because γ = 0.75, conformational gating only modestly slows ligand binding. Brownian dynamics simulations of the diffusional association of eight inhibitors to the protease─having a wide range of experimental association constants (∼10⁴-10¹⁰ M^-1 s^-1)─yields gated rate constants in the range of ∼0.5-5.7 × 10⁸ M^-1 s^-1. This indicates that, whereas the association rate of some inhibitors could be described by the model, for many inhibitors either subsequent conformational transitions or alternate binding mechanisms may be rate-limiting. For systems known to be modulated by conformational gating, the approach could be scaled computationally efficiently to screen association kinetics for a large number of ligands.