Laura Pérez-Benito

We here introduce the Aquamarine (AQM) dataset, an extensive quantum-mechanical (QM) dataset that contains the structural and electronic information of 59,783 low-and high-energy conformers of 1,653 molecules with a total number of atoms ranging from 2 to 92 (mean:50.9), and containing up to 54 (mean:28.2) non-hydrogen atoms. To gain insights into the solvent effects as well as collective dispersion interactions for drug-like molecules, we have performed QM calculations supplemented with a treatment of many-body dispersion (MBD) interactions of structures and properties in the gas phase and implicit water. Thus, AQM contains over 40 global (molecular) and local (atom-in-a-molecule) physicochemical properties (including ground-state and response properties) per conformer computed at the tightly converged PBE0+MBD level of theory for gas-phase molecules, whereas PBE0+MBD with the modified Poisson-Boltzmann (MPB) model of water was used for solvated molecules. By addressing both molecule-solvent and dispersion interactions, AQM dataset can serve as a challenging benchmark for state-of-the-art machine learning methods for property modeling and \textit{de novo} generation of large (solvated) molecules with pharmaceutical and biological relevance.

Binding free energy calculations predict the potency of compounds to protein binding sites in a physically rigorous manner and see broad application in prioritizing the synthesis of novel drug candidates. Relative binding free energy calculations (RBFE) have emerged as an industry standard approach to achieve highly accurate rank-order predictions of the potency of related compounds; however, this approach requires that the ligands share a common scaffold and a common binding mode, restricting the methods' domain of applicability. This is a critical limitation, since complex modifications to the ligands, especially core hopping, are very common in drug design. Absolute Binding Free Energy calculations (ABFE) are an alternate method, which can be used for ligands that are not congeneric. However, ABFE suffer from a known problem of long convergence times, due to the need to sample additional degrees of freedom within each system, such as sampling rearrangements necessary to open and close the binding site. Here, we report on an alternative method for RBFE, called Separated Topologies (SepTop), which overcomes the issues in both of the aforementioned methods, by enabling large scaffold changes between ligands with a convergence time comparable to traditional RBFE. Instead of only mutating atoms that vary between two ligands, this approach performs two absolute free energy calculations at the same time in opposite directions, one for each ligand. Defining the two ligands independently allows the comparison of binding of diverse ligands without the artificial constraints of identical poses or a suitable atom-atom mapping. This approach also avoids the need to sample the unbound state of the protein, making it more efficient than absolute binding free energy calculations. Here, we introduce an implementation of SepTop. We developed a general and efficient protocol for running SepTop, and we demonstrated the method on four diverse, pharmaceutically relevant systems. Here we report the performance of the method, as well as our practical insights into strengths, weaknesses, and challenges of applying this method in an industrial drug design setting. We find that the accuracy of the approach is sufficiently high to rank order ligands with an accuracy comparable to traditional RBFE calculations, while maintaining the additional flexibility of SepTop.

Free energy calculations are rapidly becoming indispensable in structure-enabled drug discovery programs. As new methods, force fields, and implementations are developed, assessing their accuracy on real-world systems (benchmarking) becomes critical to provide users with an assessment of the accuracy expected when these methods are applied within their domain of applicability, and developers with a way to assess the expected impact of new methodologies. These assessments require construction of a benchmark—a set of well-prepared, high quality systems with corresponding experimental measurements designed to ensure the resulting calculations provide a realistic assessment of expected performance. To date, the community has not yet adopted a common standardized benchmark, and existing benchmark reports suffer from a myriad of issues, including poor data quality, limited statistical power, and statistically deficient analyses, all of which can conspire to produce benchmarks that are poorly predictive of real-world performance. Here, we address these issues by presenting guidelines for (1) curating experimental data to develop meaningful benchmark sets, (2) preparing benchmark inputs according to best practices to facilitate widespread adoption, and (3) analysis of the resulting predictions to enable statistically meaningful comparisons among methods and force fields. We highlight challenges and open questions that remain to be solved in these areas, as well as recommendations for the collection of new datasets that might optimally serve to measure progress as methods become systematically more reliable. Finally, we provide a curated, versioned, open, standardized benchmark set adherent to these standards (protein-ligand-benchmark) and an open source toolkit for implementing standardized best practices assessments (arsenic) for the community to use as a standardized assessment tool. While our main focus is benchmarking free energy methods based on molecular simulations, these guidelines should prove useful for assessment of the rapidly growing field of machine learning methods for affinity prediction as well.

Background: Three amino acid differences between rodent and human APP affect medically important features, including β-secretase cleavage of APP and Aβ peptide aggregation (De Strooper et al., EMBO J 14:4932-38, 1995; Ueno et al., Biochemistry 53:7523-30, 2014; Bush, 2003, Trends Neurosci 26:207–14). Most rodent models for Alzheimer’s disease (AD) are, therefore, based on the human APP sequence, expressed from artificial mini-genes randomly inserted in the rodent genome. While these models mimic rather well various biochemical aspects of the disease, such as Aβ-aggregation, they are also prone to overexpression artifacts and to complex phenotypical alterations, due to genes affected in or close to the insertion site(s) of the mini-genes (Sasaguri et al., EMBO J 36:2473-87, 2017; Goodwin et al., Genome Res 29:494-505, 2019). Knock-in strategies which introduce clinical mutants in a humanized endogenous rodent APP sequence (Saito et al., Nat Neurosci 17:661-3, 2014) represent useful improvements, but need to be compared with appropriate humanized wildtype (WT) mice. Methods: Computational modelling of the human β-CTF bound to BACE1 was used to study the differential processing of rodent and human APP. We humanized the three pivotal residues we identified G676R, F681Y and R684H (labeled according to the human APP770 isoform) in the mouse and rat genomes using a CRISPR-Cas9 approach. These new models, termed mouse and rat App^hu/hu, express APP from the endogenous promotor. We also introduced the early-onset familial Alzheimer’s disease (FAD) mutation M139T into the endogenous Rat Psen1 gene. Results: We show that introducing these three amino acid substitutions into the rodent sequence lowers the affinity of the APP substrate for BACE1 cleavage. The effect on β-secretase processing was confirmed as both humanized rodent models produce three times more (human) Aβ compared to the original WT strain. These models represent suitable controls, or starting points, for studying the effect of transgenes or knock-in mutations on APP processing (Saito et al., Nat Neurosci 17:661-3, 2014). We introduced the early-onset familial Alzheimer’s disease (FAD) mutation M139T into the endogenous Rat Psen1 gene and provide an initial characterization of Aβ processing in this novel rat AD model. Conclusion: The different humanized APP models (rat and mouse) expressing human Aβ and PSEN1 M139T are valuable controls to study APP processing in vivo allowing the use of a human Aβ ELISA which is more sensitive than the equivalent system for rodents. These animals will be made available to the research community.

Single chemical entities with potential to simultaneously interact with two binding sites are emerging strategies in medicinal chemistry. We have designed, synthesized and functionally characterized the first bitopic ligands for the CB₂ receptor. These compounds selectively target CB₂ versus CB₁ receptors. Their binding mode was studied by molecular dynamic simulations and site‐directed mutagenesis.

The computational prediction of relative binding free energies is a crucial goal for drug discovery, and G protein-coupled receptors (GPCRs) are arguably the most important drug target class. However, they present increased complexity to model compared to soluble globular proteins. Despite breakthroughs, experimental X-ray crystal and cryo-EM structures are challenging to attain, meaning computational models of the receptor and ligand binding mode are sometimes necessary. This leads to uncertainty in understanding ligand–protein binding induced changes such as, water positioning and displacement, side chain positioning, hydrogen bond networks, and the overall structure of the hydration shell around the ligand and protein. In other words, the very elements that define structure activity relationships (SARs) and are crucial for accurate binding free energy calculations are typically more uncertain for GPCRs. In this work we use free energy perturbation (FEP) to predict the relative binding free energies for ligands of two different GPCRs. We pinpoint the key aspects for success such as the important role of key water molecules, amino acid ionization states, and the benefit of equilibration with specific ligands. Initial calculations following typical FEP setup and execution protocols delivered no correlation with experiment, but we show how results are improved in a logical and systematic way. This approach gave, in the best cases, a coefficient of determination (R²) compared with experiment in the range of 0.6–0.9 and mean unsigned errors compared to experiment of 0.6–0.7 kcal/mol. We anticipate that our findings will be applicable to other difficult-to-model protein ligand data sets and be of wide interest for the community to continue improving FE binding energy predictions.

In silico binding site location and pose prediction for a molecule targeted at a large protein surface is a challenging task. We report a blind test with two peptidomimetic molecules that bind the flu virus hemagglutinin (HA) surface antigen, JNJ7918 and JNJ4796 (recently disclosed in van Dongen et al., Science, 2019, 363). Tests with a series of conventional approaches such as rigid (receptor) docking against available X-ray crystal structures or against an ensemble of structures generated by quick methodologies (NMA, homology modeling) gave mixed results, due to the shallowness and flexibility of the binding site and the sheer size of the target. However, tests with our Monte Carlo platform PELE in two protocols involving either exploration of the whole protein surface (global exploration), or the latter followed by refinement of best solutions (local exploration) yielded remarkably good results by locating the actual binding site and generating binding modes that recovered all native contacts found in the X-ray structures. Thus, the Monte Carlo scheme of PELE seems promising as a quick methodology to overcome the challenge of identifying entirely unknown binding sites and modes for protein–protein disruptors.

Ligand binding affinity calculations based on molecular dynamics (MD) simulations and non-physical (alchemical) thermodynamic cycles have shown great promise for structure-based drug design. However, their broad uptake and impact is held back by the notoriously complex setup of the calculations. Only a few tools other than the free energy perturbation approach by Schrödinger Inc. (referred to as FEP+) currently enable end-to-end application. Here, we present for the first time an approach based on the open-source software pmx that allows to easily set up and run alchemical calculations for diverse sets of small molecules using the GROMACS MD engine. The method relies on theoretically rigorous non-equilibrium thermodynamic integration (TI) foundations, and its flexibility allows calculations with multiple force fields. In this study, results from the Amber and Charmm force fields were combined to yield a consensus outcome performing on par with the commercial FEP+ approach. A large dataset of 482 perturbations from 13 different protein–ligand datasets led to an average unsigned error (AUE) of 3.64 ± 0.14 kJ mol⁻¹, equivalent to Schrödinger's FEP+ AUE of 3.66 ± 0.14 kJ mol⁻¹. For the first time, a setup is presented for overall high precision and high accuracy relative protein–ligand alchemical free energy calculations based on open-source software.

A systematic and statistically robust protocol is applied for the evaluation of free energy calculations with and without replica‐exchange. The protocol is based on ensemble averaging to generate accurate assessments of the uncertainties in the predictions. Comparison is made between FEP+ and TIES—free energy perturbation and thermodynamic integration with enhanced sampling—the latter with and without the so‐called “enhanced sampling” based on replica‐exchange protocols. Standard TIES performs best for a reference set of targets and compounds; no benefits accrue from replica‐exchange methods. Evaluation of FEP+ and TIES with REST—replica‐exchange with solute tempering—reveals a systematic and significant underestimation of free energy differences in FEP+, which becomes increasingly large for long duration simulations, is confirmed by extensive analysis of previous publications, and raises a number of questions pertaining to the accuracy of the predictions with the REST technique not hitherto discussed.

Machine learning approach tailored for ranking congeneric series based on 3D-convolutional neural networks tested it on over 3246 ligands and 13 targets.