Free Energy Perturbation Calculations Research Articles

Ionic hydration at ambient and higher pressures: Computed chemical potentials from simulations and finite-size effects

For the chemical potential of a hydrated mono-atomic ion, finite- size effects on simulation results obtained using molecular potential truncation are investigated. Free energy perturbation (FEP) calculations were carried out by scaling in two processes the Lennard-Jones (LJ) ion-water parameters and the ion-charge (q) for Br−, K+ and Ca2+ interacting with TIP4P water. Corrections which scale with q2 enable us to reduce finite-size effects. However, at ambient conditions, discrepancies which depend on q are shown by the corrected values when comparison is made with the experimental data of the Marcus compilation. Similar behavior was observed by extrapolating the original FEP results to an infinitely large system. Hence, these errors were assumed to depend on water density and corrected at high pressures. Consistency, within statistical uncertainties, is shown when comparing with results derived from computed volumetric quantities. Results are also compared with those derived from experimental values of excess volumes at ambient conditions.

Origin of Enantioselectivity in Engineered Cytochrome c-Catalyzed Carbon-Radical FePP Hydrolysis Revealed Using QM/MM (ABEEM Polarizable Force Field) and MD Simulations.

The origin of highly efficient asymmetric aminohydroxylation of styrene catalyzed by engineered cytochrome c is investigated by the developed Atom-Bond Electronegativity Equalization Method polarizable force field (ABEEM PFF), which is a combined outcome of electronic and steric effects. Model molecules were used to establish the charge parameters of the ABEEM PFF, for which the bond-stretching and angle-bending parameters were obtained by using a combination of modified Seminario and scan methods. The interactions between carbon-radical Fe-porphyrin (FePP) and waters are simulated by molecular dynamics, which shows a clear preference for the pre-R over the pre-S. This preference is attributed to the hydrogen-bond between the mutated 100S and 101P residues as well as van der Waals interactions, enforcing a specific conformation of the carbon-radical FePP complex within the binding pocket. Meanwhile, the hydrogen-bond between water and the nitrogen atom in the active intermediate dictates the stereochemical outcome. Quantum mechanics/molecular mechanics (QM/MM (ABEEM PFF)) and free-energy perturbation calculations elucidate that the 3RTS is characterized by sandwich-like structure among adjacent amino acid residues, which exhibits greater stability than crowed arrangement in 3STS and enables the R enantiomer to form more favorably. Thus, this study provides mechanistic insight into the catalytic reaction of hemoproteins.

Fastlomap: faster lead optimization mapper algorithm for large-scale relative free energy perturbation

AbstractIn recent years, free energy perturbation calculations have garnered increasing attention as tools to support drug discovery. The lead optimization mapper (Lomap) was proposed as an algorithm to calculate the relative free energy between ligands efficiently. However, Lomap requires checking whether each edge in the FEP graph is removable, which necessitates checking the constraints for all edges. Consequently, conventional Lomap requires significant computation time, at least several hours for cases involving hundreds of compounds, and is impractical for cases with more than tens of thousands of edges. In this study, we aimed to reduce the computational cost of Lomap to enable the construction of FEP graphs for hundreds of compounds. We can reduce the overall number of constraint checks required from an amount dependent on the number of edges to one dependent on the number of nodes by using the chunk check process to check the constraints for as many edges as possible simultaneously. Based on the analysis of the execution profiles, we also improved the speed of cycle constraint and diameter constraint checks. Moreover, the output graph is the same as that obtained using the conventional Lomap, enabling direct replacement of the original one with our method. With our improvement, the execution was hundreds of times faster than that of the original Lomap.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation.

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Adaptive lambda schemes for efficient relative binding free energy calculation.

The relative free energy perturbation (RFEP) calculation is one of the most theoretically sound computational chemistry approaches for the binding affinity prediction. However, its application is often hindered by the complexity of the calculation choices and the requirement of expertise in the field. Improper lambda scheme of RFEP may result in deviations from an accurate description of the perturbation process and is prone to erroneous affinity predictions. To address such challenges, an automated adaptive lambda method is proposed where the adaptive lambda schemes are obtained through a split-and-merge algorithm based on the pilot runs. The newly established workflow along with a series of improvements to the perturbation settings increases the consistency of the RFEP calculation results. Comparing the pilot and adaptive lambda schemes, the latter demonstrated improvements in convergence and reproducibility and lowered the mean unsigned error and the root-mean-square error. Overall, the adaptive lambda method is a reliable and robust choice to predict small molecule relative binding free energy and can be capitalized to benefit routine RFEP calculations for drug discovery projects.

OPLS/2020 Force Field for Unsaturated Hydrocarbons, Alcohols, and Ethers.

The OPLS all-atom force field was updated and applied to modeling unsaturated hydrocarbons, alcohols, and ethers. Testing has included gas-phase conformational energetics, properties of pure liquids, and free energies of hydration. Monte Carlo statistical mechanics (MC) calculations were used to model 60 liquids. In addition, a robust, automated procedure was devised to compute the free energies of hydration with high precision via free-energy perturbation (FEP) calculations using double annihilation. Testing has included larger molecules than in the past, and parameters are reported for the first time for some less common groups including alkynes, allenes, dienes, and acetals. The average errors in comparison with experimental data for the computed properties of the pure liquids were improved with the modified force field (OPLS/2020). For liquid densities and heats of vaporization, the average unsigned errors are 0.01 g/cm3 and 0.2 kcal/mol. The average error and signed error for free energies of hydration are both 1.2 kcal/mol. As noted before, this reflects a systematic overestimate of the hydrophobicity of organic molecules when the parametrization is done to minimize the errors for properties of pure liquids. Implications for the modeling of biomolecular systems with standard force fields are considered.

Accurate Computation of Thermodynamic Activation Parameters in the Chorismate Mutase Reaction from Empirical Valence Bond Simulations.

Chorismate mutase (CM) enzymes have long served as model systems for benchmarking new methods and tools in computational chemistry. Despite the enzymes' prominence in the literature, the extent of the roles that activation enthalpy and entropy play in catalyzing the conversion of chorismate to prephenate is still subject to debate. Knowledge of these parameters is a key piece in fully understanding the mechanism of chorismate mutases. Within this study, we utilize EVB/MD free energy perturbation calculations at a range of temperatures, allowing us to extract activation enthalpies and entropies from an Arrhenius plot of activation free energies of the reaction catalyzed by a monofunctional Bacillus subtilis CM and the promiscuous enzyme isochorismate pyruvate lyase of Pseudomonas aeruginosa. In comparison to the uncatalyzed reaction, our results show that both enzyme-catalyzed reactions exhibit a substantial reduction in activation enthalpy, while the effect on activation entropy is relatively minor, demonstrating that enzyme-catalyzed CM reactions are enthalpically driven. Furthermore, we observe that the monofunctional CM from B. subtilis more efficiently catalyzes this reaction than its promiscuous counterpart. This is supported by a structural analysis of the reaction pathway at the transition state, from which we identified key residues explaining the enthalpically driven nature of the reactions and also the difference in efficiencies between the two enzymes.

Discovery of (R)-4-(8-methoxy-2-methyl-1-(1-phenylethy)-1H-imidazo[4,5-c]quinnolin-7-yl)-3,5-dimethylisoxazole as a potent and selective BET inhibitor for treatment of acute myeloid leukemia (AML) guided by FEP calculation

The functions of the bromodomain and extra terminal (BET) family of proteins have been proved to be involved in various diseases, particularly the acute myeloid leukemia (AML). In this work, guided by free energy perturbation (FEP) calculation, a methyl group was selected to be attached to the 1H-imidazo[4,5-c]quinoline skeleton, and a series of congeneric compounds were synthesized. Among them, compound 10 demonstrated outstanding activity against BRD4 BD1 with an IC50 value of 1.9 nM and exhibited remarkable antiproliferative effects against MV4-11 cells. The X-ray cocrystal structure proved that 10 occupied the acetylated lysine (KAc) binding cavity and the WPF shelf of BRD4 BD1. Additionally, 10 displayed high selectivity towards BET family members, effectively inhibiting the growth of AML cells, promoting apoptosis, and arresting the cell cycle at the G0/G1 phase. Further mechanistic studies demonstrated that compound 10 could suppress the expression of c-Myc and CDK6 while enhancing the expression of P21, PARP, and cleaved PARP. Moreover, 10 exhibited remarkable pharmacokinetic properties and significant antitumor efficacy in vivo. Therefore, compound 10 may represent a new, potent and selective BET bromodomain inhibitor for the development of therapeutics to treat AML.

Loops Mediate Agonist-Induced Activation of the Stimulator of Interferon Genes Protein.

The stimulator of interferon genes (STING) is an important therapeutic target for cancer diseases. The activated STING recruits downstream tank-binding kinase 1 (TBK1) to trigger several important immune responses. However, the molecular mechanism of how agonist molecules mediate the STING-TBK1 interactions remains elusive. Here, we performed molecular dynamics simulations to capture the conformational changes of STING and TBK1 upon agonist binding. Our simulations revealed that multiple helices (α5-α7) and especially three loops (loop 6, loop 8, and C-terminal tail) of STING participated in the allosteric mediation of the STING-TBK1 interactions. Consistent results were also observed in the simulations of the constitutive activating mutant of STING (R284S). We further identified α5 as a key region in this agonist-induced activation mechanism of STING. Free-energy perturbation calculations of multiple STING agonists demonstrated that an alkynyl group targeting α5 is a determinant for agonist activities. These results not only offer deeper insights into the agonist-induced allosteric mediation of STING-TKB1 interactions but also provide a guidance for future drug development of this important therapeutic target.

Dynamic extracellular vestibule of human SERT: Unveiling druggable potential with high-affinity allosteric inhibitors

The serotonin transporter (SERT) tightly regulates synaptic serotonin levels and has been the primary target of antidepressants. Binding of inhibitors to the allosteric site of human SERT (hSERT) impedes the dissociation of antidepressants bound at the central site and may enhance the efficacy of such antidepressants to potentially reduce their dosage and side effects. Here, we report the identification of a series of high-affinity allosteric inhibitors of hSERT in a unique scaffold, with the lead compound, Lu AF88273 (3-(1-(2-(1H-indol-3-yl)ethyl)piperidin-4-yl)-6-chloro-1H-indole), having 2.1 nM allosteric potency in inhibiting imipramine dissociation. In addition, we find that Lu AF88273 also inhibits serotonin transport in a noncompetitive manner. The binding pose of Lu AF88273 in the allosteric site of hSERT is determined with extensive molecular dynamics simulations and rigorous absolute binding free energy perturbation (FEP) calculations, which show that a part of the compound occupies a dynamically formed small cavity. The predicted binding location and pose are validated by site-directed mutagenesis and can explain much of the structure-activity relationship of these inhibitors using the relative binding FEP calculations. Together, our findings provide a promising lead compound and the structural basis for the development of allosteric drugs targeting hSERT. Further, they demonstrate that the divergent allosteric sites of neurotransmitter transporters can be selectively targeted.

Structure-Based Ligand Discovery Targeting the Transmembrane Domain of Frizzled Receptor FZD7.

Despite the essential roles of Frizzled receptors (FZDs) in mediating Wnt signaling in embryonic development and tissue homeostasis, ligands targeting FZDs are rare. A few antibodies and peptide modulators have been developed that mainly bind to the family-conserved extracellular cysteine-rich domain of FZDs, while the canonical binding sites in the transmembrane domain (TMD) are far from sufficiently addressed. Based on the recent structures of FZDs, we explored small-molecule ligand discovery by targeting TMD. From the ChemDiv library with ∼1.6 million compounds, we identified compound F7H as an antagonist of FZD7 with an IC50 at 1.25 ± 0.38 μM. Focusing on this hit, the structural dissection study, together with computing studies such as molecular docking, molecular dynamics simulation, and free energy perturbation calculations, defined the binding pocket with key residue recognition. Our results revealed the structural basis of ligand recognition and demonstrated the feasibility of structure-guided ligand discovery for FZD7-TMD.

Discovery of antiviral SARS-CoV-2 main protease inhibitors by structure-guided hit-to-lead optimization of carmofur

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) main protease (Mpro) has been targeted for the development of anti-SARS-CoV-2 agents against COVID-19 infection because Mpro processes essential viral polyproteins and plays a key role in SARS-CoV-2 replication. In this study, we report the development of novel SARS-CoV-2 Mpro inhibitors derived from carmofur, a previously identified compound that has shown moderate potency as a covalent inhibitor of SARS-CoV-2 Mpro. To employ a structure-guided drug design strategy, a putative intact binding mode of carmofur at catalytic active site of Mpro was initially predicted by docking simulation. Based on the predicted binding mode, a series of carmofur derivatives aiming to occupy the Mpro substrate binding regions were investigated for structure-activity relationship analysis. As a result, an indole-based derivative, speculated to interact with the S4 binding pocket, 21b (IC50 = 1.5 ± 0.1 μM) was discovered. Its structure was further modified and evaluated in silico by combining docking simulation, free energy perturbation calculation and subpocket interaction analysis to optimize the interactions at the S2 and S4 binding pockets. Among the newly designed novel derivatives, 21h and 21i showed the best inhibitory potencies against Mpro with IC50 values of 0.35 and 0.37 μM, respectively. Moreover, their antiviral activities were confirmed with EC50 values of 20–30 μM in the SARS-CoV-2-infected cell-based assay, suggesting that these novel Mpro inhibitors could be applied as potential lead compounds for the development of substantial anti-SARS-CoV-2 agents.

FEP Protocol Builder: Optimization of Free Energy Perturbation Protocols Using Active Learning.

Significant improvements have been made in the past decade to methods that rapidly and accurately predict binding affinity through free energy perturbation (FEP) calculations. This has been driven by recent advances in small-molecule force fields and sampling algorithms combined with the availability of low-cost parallel computing. Predictive accuracies of ∼1 kcal mol-1 have been regularly achieved, which are sufficient to drive potency optimization in modern drug discovery campaigns. Despite the robustness of these FEP approaches across multiple target classes, there are invariably target systems that do not display expected performance with default FEP settings. Traditionally, these systems required labor-intensive manual protocol development to arrive at parameter settings that produce a predictive FEP model. Due to the (a) relatively large parameter space to be explored, (b) significant compute requirements, and (c) limited understanding of how combinations of parameters can affect FEP performance, manual FEP protocol optimization can take weeks to months to complete, and often does not involve rigorous train-test set splits, resulting in potential overfitting. These manual FEP protocol development timelines do not coincide with tight drug discovery project timelines, essentially preventing the use of FEP calculations for these target systems. Here, we describe an automated workflow termed FEP Protocol Builder (FEP-PB) to rapidly generate accurate FEP protocols for systems that do not perform well with default settings. FEP-PB uses an active-learning workflow to iteratively search the protocol parameter space to develop accurate FEP protocols. To validate this approach, we applied it to pharmaceutically relevant systems where default FEP settings could not produce predictive models. We demonstrate that FEP-PB can rapidly generate accurate FEP protocols for the previously challenging MCL1 system with limited human intervention. We also apply FEP-PB in a real-world drug discovery setting to generate an accurate FEP protocol for the p97 system. FEP-PB is able to generate a more accurate protocol than the expert user, rapidly validating p97 as amenable to free energy calculations. Additionally, through the active-learning workflow, we are able to gain insight into which parameters are most important for a given system. These results suggest that FEP-PB is a robust tool that can aid in rapidly developing accurate FEP protocols and increasing the number of targets that are amenable to the technology.

Free Energy Perturbation Calculations of Mutation Effects on SARS-CoV-2 RBD::ACE2 Binding Affinity

The strength of binding between human angiotensin converting enzyme 2 (ACE2) and the receptor binding domain (RBD) of viral spike protein plays a role in the transmissibility of the SARS-CoV-2 virus. In this study we focus on a subset of RBD mutations that have been frequently observed in infected individuals and probe binding affinity changes to ACE2 using surface plasmon resonance (SPR) measurements and free energy perturbation (FEP) calculations. Our SPR results are largely in accord with previous studies but discrepancies do arise due to differences in experimental methods and to protocol differences even when a single method is used. Overall, we find that FEP performance is superior to that of other computational approaches examined as determined by agreement with experiment and, in particular, by its ability to identify stabilizing mutations. Moreover, the calculations successfully predict the observed cooperative stabilization of binding by the Q498R N501Y double mutant present in Omicron variants and offer a physical explanation for the underlying mechanism. Overall, our results suggest that despite the significant computational cost, FEP calculations may offer an effective strategy to understand the effects of interfacial mutations on protein–protein binding affinities and, hence, in a variety of practical applications such as the optimization of neutralizing antibodies.

Fragment Linker Prediction Using the Deep Encoder-Decoder Network for PROTACs Drug Design.

A drug discovery and development pipeline is a prolonged and complex process that remains challenging for both computational methods and medicinal chemists and has not been able to be resolved using computational methods. Deep learning has been utilized in various fields and achieved tremendous success in designing novel molecules in the pharmaceutical industry. Herein, we use state-of-the-art techniques to propose a deep neural network, AIMLinker, to rapidly design and generate meaningful drug-like proteolysis targeting chimeras (PROTACs) analogs. The model extracts the structural information from the input fragments and generates linkers to incorporate them. We integrate filters in the model to exclude nondruggable structures guided via protein-protein complexes while retaining molecules with potent chemical properties. The novel PROTACs subsequently pass through molecular docking, taking root-mean-square deviation (RMSD), relative Gibbs free energy (ΔΔGbinding), molecular dynamics (MD) simulation, and free energy perturbation (FEP) calculations as the measurement criteria for testing the robustness and feasibility of the model. The generated novel PROTACs molecules possess similar structural information with superior binding affinity to the binding pockets compared to the existing CRBN-dBET6-BRD4 ternary complexes. We demonstrate the effectiveness of the methodology of leveraging AIMLinker to design novel compounds for PROTACs molecules exhibiting better chemical properties compared to the dBET6 crystal pose.

Unbiased molecular dynamics simulation of a first-in-class small molecule inhibitor binds to oncostatin M

Small molecule inhibitors (SMIs) targeting oncostatin M (OSM) signaling pathway represent new therapeutics to combat cancer, inflammatory bowel disease (IBD) and CNS disease. Recently, the first-in-class SMI named SMI-10B that target OSM and block its interaction with receptor (OSMR) were reported. However, the binding pocket and interaction mode of the compound on OSM remain poorly understood, which hampering the rational design of SMIs that target OSM. Here, using SMI-10B as a probe, the multiple pockets on OSM for small molecules binding were extensively explored by unbiased molecular dynamics (MD) simulations. Then, the near-native structure of the complex was identified by molecular mechanics generalized Born surface area (MM/GBSA) binding energy funnel. Moreover, the binding stabilities of the protein-ligand complexes in near- and non-native conformations were verified by additional independent MD runs and absolute free energy perturbation (FEP) calculation. In summary, the unique feature of SMI-10B spontaneously binds to OSM characterized here not only provide detailed information for understanding the molecular mechanism of SMI-10B binding to OSM, but also will facilitate the rational design of novel and more potent SMIs to block OSM signaling.

Recent Advances in Alchemical Binding Free Energy Calculations for Drug Discovery.

Rigorous physics-based methods to calculate binding free energies of protein-ligand complexes have become a valued component of structure-based drug design. Relative and absolute binding free energy calculations have been deployed prospectively in support of solving diverse drug discovery challenges. Here we review recent applications of binding free energy calculations to fragment growing and linking, scaffold hopping, binding pose validation, virtual screening, covalent enzyme inhibition, and positional analogue scanning. Furthermore, we discuss the merits of using protein models and highlight recent efforts to replace costly binding free energy calculations with predictions from machine learning models trained on a limited number of free energy perturbation or thermodynamic integration calculations thereby allowing for extended chemical space exploration.

Improving the CHARMM36m force field for arginine-phosphate interactions.

Noncovalent binding interactions between arginine and phosphate, such as the anchoring of transmembrane helices to phospholipid membrane interfaces or the binding of SH2 domain proteins with phosphorylated peptides, are essential for structural integrity or signal recognition in many biological processes. Using 120 ns long MD simulations, we find that arginine and phosphate have multiple binding conformations. The binding conformations are clustered into five groups and we compute binding energies for these clusters. Free energy perturbation (FEP) calculations are used to estimate the binding free energy using the current CHARMM36m and its CUFIX corrected version. Umbrella sampling (US) is used to generate the potential of mean force (PMF) along the inter-group distance reaction coordinate. Our results were compared with data from isothermal calorimetry (ITC) and suggest that neither the CHARMM36m nor the CUFIX-correct version match experiment and the ITC value is in between these force field estimates. Our current work focuses on the generation of improved force field parameters for arginine-phosphate interactions and experimental validation of these parameters for peptide-liposome binding.

Regioselective metabolism of cannabinoids by cytochrome P450 2J2.

Cannabichromene (CBC) is a major non-psychoactive phytocannabinoid extracted from the cannabis plant. CBC has been shown to have anti-inflammatory effects; hence it can serve as a non-narcotic alternative to tetrahydrocannabinol (THC). Nevertheless, the metabolic pathways and metabolites of CBC have not been fully characterized. Metabolism of phytocannabinoids typically occurs via a superfamily of hemoproteins called cytochrome P450s (CYP450s). Our liquid chromatography-mass spectrometry (LC/MS) experiments show that CBC is metabolized to three different products by various human CYP450s. We identified that CYP450 2J2 (CYP2J2) predominantly produces 8'-hydroxy-CBC, an allylic hydroxide, out of the three metabolites. To investigate the CBC metabolic regioselectivity in CYP2J2, we leveraged atomistic molecular dynamics (MD) and modeling of the system. Using 200 unique conformations of CYP2J2 obtained from MD, we performed ensemble molecular docking of CBC to reveal three putative binding modes. We then assessed the stability of the obtained CBC binding modes by simulating them individually for 1 μs. Our simulations captured an average distance of 4.8 Å between the CYP2J2 heme and the 8’ carbon atom on CBC for one of the predominant binding modes, at which the allylic hydroxide metabolite occurs experimentally. At the regioselective binding mode, we observe conserved hydrogen bonding networks, mainly between the CBC hydroxyl group and the backbone oxygen of I487, as well as N379, stabilized by a water molecule. We further corroborated the stability of the regioselective binding modes by calculating relative binding free energies using alchemical free energy perturbation calculations. In conjunction, our computational and experimental methodologies elucidate a likely mode for regioselective binding and metabolism of CBC by CYP2J2, furthering the way towards accessible plant-based non-psychoactive anti-inflammatory therapeutics.

MemCross: Accelerated Weight Histogram method to assess membrane permeability

Passive permeation events across biological membranes are determining steps in the pharmacokinetics of xenobiotics. To reach an accurate and rapid prediction of membrane permeation coefficients of drugs is a complex challenge, which can efficiently support drug discovery. Such predictions are indeed highly valuable as they may guide the selection of potential leads with optimum bioavailabilities prior to synthesis. Theoretical models exist to predict these coefficients. Many of them are based on molecular dynamics (MD) simulations, which allow calculation of permeation coefficients through the evaluation of both the potential of mean force (PMF) and the diffusivity profiles. However, these simulations still require intensive computational efforts, and novel methodologies should be developed and benchmarked. Free energy perturbation (FEP) method was recently shown to estimate PMF with a significantly reduced computational cost compared to the adaptive biasing force method. This benchmarking was achieved with small molecules, namely short-chain alcohols. Here, we show that to estimate the PMF of bulkier, drug-like xenobiotics, conformational sampling is a critical issue. To reach a sufficient sampling with FEP calculations requires a relatively long time-scale, which can lower the benefits related to the computational gain. In the present work, the Accelerated Weight Histogram (AWH) method was employed for the first time in all-atom membrane models. The AWH-based protocol, named MemCross, appears affordable to estimate PMF profiles of a series of drug-like xenobiotics, compared to other enhanced sampling methods. The continuous exploration of the crossing pathway by MemCross also allows modeling subdiffusion by computing fractional diffusivity profiles. The method is also versatile as its input parameters are largely insensitive to the molecule properties. It also ensures a detailed description of the molecule orientations along the permeation pathway, picturing all intermolecular interactions at an atomic resolution. Here, MemCross was applied on a series of 12 xenobiotics, including four weak acids, and a coherent structure-activity relationship was established.