Biased Molecular Dynamics Simulations Research Articles

Potential acridinedione derivatives for the development of a heterobifunctional PROTAC for targeted degradation of PCSK9 protein

Proprotein convertase substilisin-like/kexin type 9 (PCSK9) interacts with the low-density lipoprotein (LDL) receptor by protein-protein interaction (PPI). The inhibition of this PPI lowers the LDL cholesterol levels, thereby reducing the risks of cardiovascular diseases. We identified in-house synthesized acridinedione scaffolds as binders of PCSK9 through microsecond timescale conventional and biased molecular dynamics (MD) simulations. The free energy of binding was calculated by the Molecular Mechanics Poisson-Boltzmann Surface Area and umbrella sampling methods for acridinedione molecules and compared with reference molecules (co-crystallized inhibitors). Three acridinedione molecules (DSPD-2, DSPD-4, and DSPD-6) were identified as potential binders of the PCSK9 allosteric binding site from preliminary docking and steered MD simulations. The stability of molecules inside the allosteric pocket was established by the low values of root mean square deviations of backbone C-α atoms. The acridinedione derivatives also showed a shorter distance of hydrogen bonds as compared to the standard molecules. The molecule DSPD-6 showed the lowest binding free energy as compared to other selected and standard molecules. The binding of the selected molecules at the allosteric binding pocket had no effect on structural parameters of LDL receptor-binding site of PCSK9 protein. Finally, a proteolysis targeting chimera (PROTAC) was developed by attaching a proteasome recruiting tag to DSPD-6 (lowest binding affinity) to achieve ligand-induced protein degradation. These results qualify the development of the DSPD-6 molecule as a heterobifunctional PROTAC for targeted degradation of PCSK9 protein. Moreover, our computational strategy provides a framework for the identification of molecules against conventionally undruggable PPI targets.

Molecular Mechanisms Underlying Medium-Chain Free Fatty Acid-Regulated Activity of the Phospholipase PlaF from Pseudomonas aeruginosa.

PlaF is a membrane-bound phospholipase A1 from Pseudomonas aeruginosa that is involved in remodeling membrane glycerophospholipids (GPLs) and modulating virulence-associated signaling and metabolic pathways. Previously, we identified the role of medium-chain free fatty acids (FFAs) in inhibiting PlaF activity and promoting homodimerization, yet the underlying molecular mechanism remained elusive. Here, we used unbiased and biased molecular dynamics simulations and free energy computations to assess how PlaF interacts with FFAs localized in the water milieu surrounding the bilayer or within the bilayer and how these interactions regulate PlaF activity. Medium-chain FFAs localized in the upper bilayer leaflet can stabilize inactive dimeric PlaF, likely through interactions with charged surface residues, as has been experimentally validated. Potential of mean force (PMF) computations indicate that membrane-bound FFAs may facilitate the activation of monomeric PlaF by lowering the activation barrier for changing into a tilted, active configuration. We estimated that the coupled equilibria of PlaF monomerization-dimerization and tilting at the physiological concentration of PlaF lead to the majority of PlaF forming inactive dimers when in a cell membrane loaded with decanoic acid (C10). This is in agreement with a suggested in vivo product feedback loop and gas chromatography-mass spectrometry profiling results, indicating that PlaF catalyzes the release of C10 from P. aeruginosa membranes. Additionally, we found that C10 in the water milieu can access the catalytic site of active monomeric PlaF, contributing to the competitive component of C10-mediated PlaF inhibition. Our study provides mechanistic insights into how medium-chain FFAs may regulate the activity of PlaF, a potential bacterial drug target.

How is Membrane Permeation of Small Ionizable Molecules Affected by Protonation Kinetics?

According to the pH-partition hypothesis, the aqueous solution adjacent to a membrane is a mixture of the ionization states of the permeating molecule at fixed Henderson-Hasselbalch concentrations, such that each state passes through the membrane in parallel with its own specific permeability. An alternative view, based on the assumption that the rate of switching ionization states is instantaneous, represents the permeation of ionizable molecules via an effective Boltzmann-weighted average potential (BWAP). Such an assumption is used in constant-pH molecular dynamics simulations. The inhomogeneous solubility-diffusion framework can be used to compute the pH-dependent membrane permeability for each of these two limiting treatments. With biased WTM-eABF molecular dynamics simulations, we computed the potential of mean force and diffusivity of each ionization state of two weakly basic small molecules: nicotine, an addictive drug, and varenicline, a therapeutic for treating nicotine addiction. At pH = 7, the BWAP effective permeability is greater than that determined by pH-partitioning by a factor of 2.5 for nicotine and 5 for varenicline. To assess the importance of ionization kinetics, we present a Smoluchowski master equation that includes explicitly the protonation and deprotonation processes coupled with the diffusive motion across the membrane. At pH = 7, the increase in permeability due to the explicit ionization kinetics is negligible for both nicotine and varenicline. This finding is reaffirmed by combined Brownian dynamics and Markov state model simulations for estimating the permeability of nicotine while allowing changes in its ionization state. We conclude that for these molecules the pH-partition hypothesis correctly captures the physics of the permeation process. The small free energy barriers for the permeation of nicotine and varenicline in their deprotonated neutral forms play a crucial role in establishing the validity of the pH-partitioning mechanism. Essentially, BWAP fails because ionization kinetics are too slow on the time scale of membrane crossing to affect the permeation of small ionizable molecules such as nicotine and varenicline. For the singly protonated state of nicotine, the computational results agree well with experimental measurements (P1 = 1.29 × 10-7 cm/s), but the agreement for neutral (P0 = 6.12 cm/s) and doubly protonated nicotine (P2 = 3.70 × 10-13 cm/s) is slightly worse, likely due to factors associated with the aqueous boundary layer (neutral form) or leaks through paracellular pathways (doubly protonated form).

Estimating Free Energy Barriers for Heterogeneous Catalytic Reactions with Machine Learning Potentials and Umbrella Integration.

Predicting the rate constants of elementary reaction steps is key for the computational modeling of catalytic processes. Within transition state theory (TST), this requires an accurate estimation of the corresponding free energy barriers. While sophisticated methods for estimating free energy differences exist, these typically require extensive (biased) molecular dynamics simulations that are computationally prohibitive with the first-principles electronic structure methods that are typically used in catalysis research. In this contribution, we show that machine-learning (ML) interatomic potentials can be trained in an automated iterative workflow to perform such free energy calculations at a much reduced computational cost as compared to a direct density functional theory (DFT) based evaluation. For the decomposition of CHO on Rh(111), we find that thermal effects are substantial and lead to a decrease in the free energy barrier, which can be vanishingly small, depending on the DFT functional used. This is in stark contrast to previously reported estimates based on a harmonic TST approximation, which predicted an increase in the barrier at elevated temperatures. Since CHO is the reactant of the putative rate limiting reaction step in syngas conversion on Rh(111) and essential for the selectivity toward oxygenates containing multiple carbon atoms (C2+ oxygenates), our results call into question the reported mechanism established by microkinetic models.

Modeling the Cellular Uptake of Functionalized Carbon Nanohorns Loaded with Cisplatin through a Breast Cancer Cell Membrane.

The cisplatin encapsulation into carbon nanohorns (CNH) is a promising nanoformulation to circumvent the drug dissipation and to specifically accumulate it in tumor sites. Herein, biased molecular dynamics simulations were used to analyze the transmembrane transport of the CNH loaded with cisplatin through a breast cancer cell membrane prototype. The simulations revealed a four-stage mechanism: approach, insertion, permeation, and internalization. Despite the lowest structural disturbance of the membrane provided by the nanocarrier, the average free energy barrier for the translocation was 55.2 kcal mol-1, suggesting that the passive process is kinetically unfavorable. In contrast, the free energy profiles revealed potential wells of -6.8 kcal mol-1 along the insertion stage in the polar heads region of the membrane, which might enhance the retention of the drug in tumor sites; therefore, the most likely cisplatin delivery mechanism should involve the adsorption and retention of CNH on the surface of cancer cells, allowing the loaded cisplatin be slowly released and passively transported through the cell membrane.

Autoencoders for dimensionality reduction in molecular dynamics: Collective variable dimension, biasing, and transition states.

The heat shock protein 90 (Hsp90) is a molecular chaperone that controls the folding and activation of client proteins using the free energy of ATP hydrolysis. The Hsp90 active site is in its N-terminal domain (NTD). Our goal is to characterize the dynamics of NTD using an autoencoder-learned collective variable (CV) in conjunction with adaptive biasing force Langevin dynamics. Using dihedral analysis, we cluster all available experimental Hsp90 NTD structures into distinct native states. We then perform unbiased molecular dynamics (MD) simulations to construct a dataset that represents each state and use this dataset to train an autoencoder. Two autoencoder architectures are considered, with one and two hidden layers, respectively, and bottlenecks of dimension k ranging from 1 to 10. We demonstrate that the addition of an extra hidden layer does not significantly improve the performance, while it leads to complicated CVs that increase the computational cost of biased MD calculations. In addition, a two-dimensional (2D) bottleneck can provide enough information of the different states, while the optimal bottleneck dimension is five. For the 2D bottleneck, the 2D CV is directly used in biased MD simulations. For the five-dimensional (5D) bottleneck, we perform an analysis of the latent CV space and identify the pair of CV coordinates that best separates the states of Hsp90. Interestingly, selecting a 2D CV out of the 5D CV space leads to better results than directly learning a 2D CV and allows observation of transitions between native states when running free energy biased dynamics.

CpSRP43 Is Both Highly Flexible and Stable: Structural Insights Using a Combined Experimental and Computational Approach.

The novel multidomain protein, cpSRP43, is a unique subunit of the post-translational chloroplast signal recognition particle (cpSRP) targeting pathway in higher plants. The cpSRP pathway is responsible for targeting and insertion of light-harvesting chlorophyll a/b binding proteins (LHCPs) to the thylakoid membrane. Upon emergence into the stroma, LHCPs form a soluble transit complex with the cpSRP heterodimer, which is composed of cpSRP43 and cpSRP54. cpSRP43 is irreplaceable as a chaperone to LHCPs in their translocation to the thylakoid membrane and remarkable in its ability to dissolve aggregates of LHCPs without the need for external energy input. In previous studies, cpSRP43 has demonstrated significant flexibility and interdomain dynamics. In this study, we explore the structural stability and flexibility of cpSRP43 using a combination of computational and experimental techniques and find that this protein is concurrently highly stable and flexible. In addition to microsecond-level unbiased molecular dynamics (MD), biased MD simulations based on system-specific collective variables are used along with biophysical experimentation to explain the basis of the flexibility and stability of cpSRP43, showing that the free and cpSRP54-bound cpSRP43 has substantially different conformations and conformational dynamics.

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.

Interaction at a distance: Xenon migration in Mb

The transport of ligands, such as NO or O2, through internal cavities is essential for the function of globular proteins, including hemoglobin, myoglobin (Mb), neuroglobin, truncated hemoglobins, or cytoglobin. For Mb, several internal cavities (Xe1 through Xe4) were observed experimentally and they were linked to ligand storage. The present work determines barriers for xenon diffusion and relative stabilization energies for the ligand in the initial and final pocket, linking a transition depending on the occupancy state of the remaining pockets from both biased and unbiased molecular dynamics simulations. It is found that the energetics of a particular ligand migration pathway may depend on the direction in which the transition is followed and the occupancy state of the other cavities. Furthermore, the barrier height for a particular transition can depend in a non-additive fashion on the occupancy of either cavity A or B or simultaneous population of both cavities, A and B. Multiple repeats for the Xe1 → Xe2 transition reveal that the activation barrier is a distribution of barrier heights rather than one single value, which is confirmed by a distribution of transition times for the same transition from unbiased simulations. Dynamic cross correlation maps demonstrate that correlated motions occur between adjacent residues or through space, residue Phe138 is found to be a gate for the Xe1 → Xe2 transition, and the volumes of the internal cavities vary along the diffusion pathway, indicating that there is dynamic communication between the ligand and the protein. These findings suggest that Mb is an allosteric protein.

Path separation of dissipation-corrected targeted molecular dynamics simulations of protein-ligand unbinding.

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.

Integrating microsecond timescale classical and biased molecular dynamics simulations to screen potential molecules for BRD4-BD1

Contemporary pre-clinical drug discovery is complemented by the combinatorial use of numerous computational techniques linked to computer-aided drug design and discovery. We combined inherently different computational analyses derived from multi-scale conventional, steered and umbrella sampling simulations to identify potential drug candidates for the first bromodomain of the human bromodomain containing protein 4 (BRD4-BD1). The BRD4-BD1 protein is a recognized therapeutic target for various diseases including cancer, neurological disorders, inflammation, and obesity. The in-house synthesized benzosuberene-sulfone (BSS) molecules were assessed for their potential to act as inhibitors against the BRD4-BD1 protein. We compared the molecular interactions and docking scores of BSS analogues with six different co-crystal inhibitors of the BRD4-BD1 protein. Further, the thermodynamic free energy of binding for each complex was estimated by the end-state Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) approach. Moreover, umbrella sampling simulations were utilized to validate the results obtained from MM/PBSA analyses and to demonstrate the unbinding/binding process of the most potent BSS analogue from the binding pocket of BRD4-BD1 protein. Our results show the potential of organosulfur molecules to be developed as inhibitors of the BRD4-BD1 protein and endorse their evaluation by suitable in-vitro and in-vivo studies.

Binding affinity estimation from restrained umbrella sampling simulations

The protein–ligand binding affinity quantifies the binding strength between a protein and its ligand. Computer modeling and simulations can be used to estimate the binding affinity or binding free energy using data- or physics-driven methods or a combination thereof. Here we discuss a purely physics-based sampling approach based on biased molecular dynamics simulations. Our proposed method generalizes and simplifies previously suggested stratification strategies that use umbrella sampling or other enhanced sampling simulations with additional collective-variable-based restraints. The approach presented here uses a flexible scheme that can be easily tailored for any system of interest. We estimate the binding affinity of human fibroblast growth factor 1 to heparin hexasaccharide based on the available crystal structure of the complex as the initial model and four different variations of the proposed method to compare against the experimentally determined binding affinity obtained from isothermal titration calorimetry experiments.

A look at the face of the molten globule: Structural model of the Helicobacter pylori apoflavodoxin ensemble at acidic pH.

Molten globule (MG) is the name given to a compact, non‐native conformation of proteins that has stimulated the imagination and work in the protein folding field for more than 40 years. The MG has been proposed to play a central role in the folding reaction and in important cell functions, and to be related to the onset of misfolding diseases. Due to its inherent intractability to high‐resolution studies, atomistic structural models have not yet been obtained. We present here an integrative atomistic model of the MG formed at acidic pH by the apoflavodoxin from the human pathogen Helicobacter pylori. This MG has been previously shown to exhibit the archetypical expansion, spectroscopic and thermodynamic features of a molten conformation. To obtain the model, we have analyzed the stability of wild‐type and 55 apoflavodoxin mutants to derive experimental equilibrium Φ values that have been used in biased molecular dynamics simulations to convert the native conformation into an MG ensemble. The ensemble has been refined to reproduce the experimental hydrodynamic radius and circular dichroism (CD) spectrum. The refined ensemble, deposited in PDB‐Dev, successfully explains the characteristic 1H‐nuclear magnetic resonance (NMR) and near‐UV CD spectral features of the MG as well as its solvent‐accessible surface area (SASA) change upon unfolding. This integrative model of an MG will help to understand the energetics and roles of these elusive conformations in protein folding and misfolding. Interestingly, the apoflavodoxin MG is structurally unrelated to previously described partly unfolded conformations of this protein, exemplifying that equilibrium MGs need not to reflect the properties of kinetic intermediates.

Molecular Dynamics Simulations of Interactions between Human Telomeric i-Motif Deoxyribonucleic Acid and Functionalized Graphene.

The work deals with molecular dynamics (MD) simulations of protonated, human telomeric i-motif deoxyribonucleic acid (DNA) with functionalized graphene. We studied three different graphene sheets: unmodified graphene with hydrogen atoms attached to their edges and two functionalized ones. The functionalization of graphene edge consists in attaching partially protonated or dissociated amine and carboxyl groups. We found that in all cases the protonated i-motif adsorbs strongly on the graphene surface. The biased MD simulations showed that the work necessary to drag the i-motif out from amine-doped graphene is about twice larger than that in other cases. In general, the system i-motif/amine-doped graphene stands out from the rest, e.g., in this case, the i-motif adsorbs its side with 3′ and 5′ ends oriented in the opposite to surface direction. In other cases, the DNA fragment is adsorbed to graphene by 3′ and 5′ ends. In all cases, the adsorption on graphene influences the i-motif internal structure by changing the distances between i-motif strands as well as stretching or shortening the DNA chain, but only in the case of amine-doped graphene the adsorption affects internal H-bonds formed between nucleotides inside the i-motif structure.

Molecular simulations of fluoxetine in hydrated lipid bilayers, as well as in aqueous solutions containing β-cyclodextrin

Fluoxetine, which is a well-known antidepressant drug, is studied in hydrated cholesterol-free and cholesterol-containing lipid bilayers through unbiased and biased atomistic molecular dynamics simulations. The latter are conducted for the calculation of the potential of mean force (PMF) of fluoxetine along an axis perpendicular to the two leaflets of the bilayer. The PMF indicates that the drug prefers to reside inside the lipid phase and allows us to calculate important thermodynamic properties, such as the Gibbs energy difference of partitioning from the water to the lipid phase and the Gibbs energy barrier for hopping events between the two leaflets of the bilayer. The results from the biased simulations are in accord with the mass density profiles calculated from the unbiased simulations. Moreover, we estimate the effect of fluoxetine mole fraction on the order parameters of the lipid alkyl chains and on the area per lipid. It is also found that fluoxetine forms a hydrogen bond network with lipids and water molecules penetrating into the lipid phase. In addition, fluoxoetine is studied in detail in aqueous solutions containing β-cyclodextrin. It is observed from unbiased molecular dynamics simulations that the two aforementioned molecules form a noncovalent complex spontaneously and the calculated binding free energy is in agreement with the literature.

Counteractive Effects of Choline Geranate (CAGE) ILs and Ethanol on Insulin’s Stability—A Leap Forward towards Oral Insulin Formulation

Choline geranate (CAGE) ionic liquids (ILs) stabilize insulin, thereby aiding its oral delivery, whereas ethanol (EtOH) affects its stability by disrupting the hydrophobic interactions. In this study, cognizance of the stabilization mechanism of insulin dimer in the presence of both CAGE ILs and EtOH mixtures is achieved through biased and unbiased molecular dynamics (MD) simulations. Here, two order parameters are employed to study the insulin dimer dissociation using well-tempered metadynamics (WT-MetaD). The stability of insulin is found to be strongly maintained until a 0.20 mole fraction of EtOH. Besides, higher concentrations of EtOH marginally affect the insulin stability. Moreover, geranate anions form a higher number of H-bonding interactions with water molecules, which aids insulin stabilization. Conversely, the addition of EtOH minimizes the water-mediated H-bonding interactions of geranate. Additionally, geranate traps the EtOH molecules, thereby preventing the interactions between insulin and EtOH. Furthermore, the free energy landscape (FEL) reveals the absence of dimer dissociation along with noticeable deviations in the distances R and the number of contacts Q. The dimerization free energy of insulin was calculated to be −16.1 kcal/mol at a 0.20 mole fraction of EtOH. Moreover, increments in mole fractions of EtOH effectuate a decrease in the insulin stability. Thus, the present study represents CAGE ILs as efficient insulin dimer stabilizes at low concentrations of EtOH.

Free Energy Landscape of H2A-H2B Displacement From Nucleosome

Nucleosome reconstitution plays an important role in many cellular functions. As an initial step, H2A-H2B dimer displacement, which is accompanied by disruption of many of the interactions within the nucleosome, should occur. To understand how H2A-H2B dimer displacement occurs, an adaptively biased molecular dynamics (ABMD) simulation was carried out to generate a variety of displacements of the H2A-H2B dimer from the fully wrapped to partially unwrapped nucleosome structures. With regards to these structures, the free energy landscape of the dimer displacement was investigated using umbrella sampling simulations. We found that the main contributors to the free energy were the docking domain of H2A and the C-terminal of H4. There were various paths for the dimer displacement which were dependent on the extent of nucleosomal DNA wrapping, suggesting that modulation of the intra-nucleosomal interaction by external factors such as histone chaperones could control the path for the H2A-H2B dimer displacement. Key residues which contributed to the free energy have also been reported to be involved in the mutations and posttranslational modifications (PTMs) which are important for assembling and/or reassembling the nucleosome at the molecular level and are found in cancer cells at the phenotypic level. Our results give insight into how the H2A-H2B dimer displacement proceeds along various paths according to different interactions within the nucleosome.

Unveiling the Releasing Processes of Pt(II)-Based Anticancer Drugs from Oxidized Carbon Nanohorn: An In Silico Study.

About half of all cancer chemotherapies currently applied involve medication with the three worldwide approved Pt(II)-based drugs, cisplatin (cddp), carboplatin (cpx), and oxaliplatin (oxa), due to their notable antitumor activity for several cancers. However, this wide application is accompanied by severe side effects, such as nephrotoxicity, myelosuppression, and neurotoxicity, as a result of their low bioavailability and selectivity for cancer cells. To mitigate these drawbacks, the use of chemically functionalized carbon nanohorns (CNH) as nanocarriers represents a potential formulation since CNH has been noted for their biodegradability, biocompatibility, low toxicity, and cavities dimensionally compatible with small drugs. This work reports energetic and dynamic analyses of complexes formed by oxidized CNH (CNHox) and the cddp, cpx, and oxa drugs. Using unbiased molecular dynamics (MD) simulations, we show that the encapsulated formulations (cddp@CNHox, cpx@CNHox, and oxa@CNHox) were more stable by ∼11.0 kcal mol-1 than the adsorbed ones (cddp > CNHox, cpx > CNHox, and oxa > CNHox). This high stability, mainly governed by van der Waals interactions, was responsible for the drug confinement during the entire simulation time (200 ns). The biased MD simulations of the inclusion complexes confirmed the nonspontaneity of the drug release since the potentials of mean force (PMF) indicated the endergonic character of this process. Additionally, the releasing energy profiles pointed out that the free energy barrier (ΔΔG≠) for the escape from CNHox cavity follows the order oxa > cpx ∼ cddp, with the value for the oxa complex (21-26 kcal mol-1) found to be about 36 and 30% larger than those for cpx and cddp, respectively. While the approximate residence time (tres) of the oxa drug inside the CNHox cavity was 5.45 × 108 s, the same measure for the cddp and cpx drugs was 5.3 × 105 and 1.60 × 103 s. Simulations also revealed that the escape of oxa with the oxalate group facing the nanowindow was the most unfavorable process, giving tres = 1.09 × 109 s. Besides reinforcing and extending the nanovectorization of cddp, cpx, and oxa in CNHox for cancer chemotherapies, all features considered may provide interpretations for experimental data and encourage new investigations aiming to propose less aggressive treatments for oncological diseases.

Molecular mechanism of allosteric modulation for the cannabinoid receptor CB1.

Given the promising clinical value of allosteric modulators of G protein-coupled-receptors (GPCRs), mechanistic understanding of how these modulators alter GPCR function is of significance. Here, we report the crystallographic and cryo-electron microscopy structures of the cannabinoid receptor CB1 bound to the positive allosteric modulator (PAM) ZCZ011. These structures show that ZCZ011 binds to an extrahelical site in the transmembrane 2 (TM2)-TM3-TM4 surface. Through (un)biased molecular dynamics simulations and mutagenesis experiments, we show that TM2 rearrangement is critical for the propagation of allosteric signals. ZCZ011 exerts a PAM effect by promoting TM2 rearrangement in favor of receptor activation and increasing the population of receptors that adopt an active conformation. In contrast, ORG27569, a negative allosteric modulator (NAM) of CB1, also binds to the TM2-TM3-TM4 surface and exerts a NAM effect by impeding the TM2 rearrangement. Our findings fill a gap in the understanding of CB1 allosteric regulation and could guide the rational design of CB1 allosteric modulators.

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.