Electron and Spin Delocalization in [Co6 Se8 (PEt3 )6 ]0/+1 Superatoms.

Molecular clusters can function as nanoscale atoms/superatoms, assembling into superatomic solids, a new class of solid-state materials with designable properties through modifications on superatoms. To explore possibilities on diversifying building blocks, here we thoroughly studied one representative superatom, Co6 Se8 (PEt3 )6 . We probed its structural, electronic, and magnetic properties and revealed its detailed electronic structure as valence electrons delocalize over inorganic [Co6 Se8 ] core while ligands function as an insulated shell. 59 Co SSNMR measurements on the core and 31 P, 13 C on the ligands show that the neutral Co6 Se8 (PEt3 )6 is diamagnetic and symmetric, with all ligands magnetically equivalent. Quantum computations cross-validate NMR results and reveal degenerate delocalized HOMO orbitals, indicating aromaticity. Ligand substitution keeps the inorganic core nearly intact. After losing one electron, the unpaired electron in [Co6 Se8 (PEt3 )6 ]+1 is delocalized, causing paramagnetism and a delocalized electron spin. Notably, this feature of electron/spin delocalization over a large cluster is attractive for special single-electron devices.

Highly efficient implementation of analytic nonadiabatic derivative couplings within the pseudospectral method.

A pseudospectral implementation of nonadiabatic derivative couplings in the Tamm-Dancoff approximation is reported, and the accuracy and efficiency of the pseudospectral nonadiabatic derivative couplings are studied. Our results demonstrate that the pseudospectral method provides mean absolute errors of 0.2%-1.9%, while providing a significant speedup. Benchmark calculations on fullerenes (Cn, n up to 100) using B3LYP achieved 10- to 15-fold, 8- to 17-fold, and 43- to 75-fold speedups for 6-31G**, 6-31++G**, and cc-pVTZ basis sets, respectively, when compared to the conventional spectral method.

Enhancing Hit Discovery in Virtual Screening through Absolute Protein-Ligand Binding Free-Energy Calculations.

In the hit identification stage of drug discovery, a diverse chemical space needs to be explored to identify initial hits. Contrary to empirical scoring functions, absolute protein-ligand binding free-energy perturbation (ABFEP) provides a theoretically more rigorous and accurate description of protein-ligand binding thermodynamics and could, in principle, greatly improve the hit rates in virtual screening. In this work, we describe an implementation of an accurate and reliable ABFEP method in FEP+. We validated the ABFEP method on eight congeneric compound series binding to eight protein receptors including both neutral and charged ligands. For ligands with net charges, the alchemical ion approach is adopted to avoid artifacts in electrostatic potential energy calculations. The calculated binding free energies correlate with experimental results with a weighted average of R2 = 0.55 for the entire dataset. We also observe an overall root-mean-square error (RMSE) of 1.1 kcal/mol after shifting the zero-point of the simulation data to match the average experimental values. Through ABFEP calculations using apo versus holo protein structures, we demonstrated that the protein conformational and protonation state changes between the apo and holo proteins are the main physical factors contributing to the protein reorganization free energy manifested by the overestimation of raw ABFEP calculated binding free energies using the holo structures of the proteins. Furthermore, we performed ABFEP calculations in three virtual screening applications for hit enrichment. ABFEP greatly improves the hit rates as compared to docking scores or other methods like metadynamics. The good performance of ABFEP in rank ordering compounds demonstrated in this work confirms it as a useful tool to improve the hit rates in virtual screening, thus facilitating hit discovery.

Accurate Quantum Chemical Reaction Energies for Lithium-Mediated Electrolyte Decomposition and Evaluation of Density Functional Approximations.

An important concern related to the performance of Li-ion batteries is the formation of a solid electrolyte interphase on the surface of the anode. This film is formed from the decomposition of electrolytes and can have important effects on the stability and performance. Here, we evaluate the decomposition pathway of ethylene carbonate and related organic electrolyte molecules using a series of density functional approximations and correlated wave function (WF) methods, including the coupled-cluster theory with single, double, and perturbative triple excitations [CCSD(T)] and auxiliary-field quantum Monte Carlo (AFQMC). We find that the transition state barrier associated with ring opening varies widely across different functionals, ranging from 3.01 to 17.15 kcal/mol, which can be compared to the value of 12.84 kcal/mol predicted by CCSD(T). This large variation underscores the importance of benchmarking against accurate WF methods. A performance comparison of all of the density functionals used in this study reveals that the M06-2X-D3 (a meta-hybrid GGA), CAM-B3LYP-D3 (a range-separated hybrid), and B2GP-PLYP-D3 (a double hybrid) perform the best, with average errors of about 1.50-1.60 kcal/mol compared to CCSD(T). We also compared the performance of the WF methods that are more scalable than CCSD(T), finding that DLPNO-CCSD(T) and phaseless AFQMC with a DFT trial wave function exhibit average errors of 1.38 and 1.74 kcal/mol, respectively.

On the potentially transformative role of auxiliary-field quantum Monte Carlo in quantum chemistry: A highly accurate method for transition metals and beyond.

Approximate solutions to the ab initio electronic structure problem have been a focus of theoretical and computational chemistry research for much of the past century, with the goal of predicting relevant energy differences to within "chemical accuracy" (1 kcal/mol). For small organic molecules, or in general, for weakly correlated main group chemistry, a hierarchy of single-reference wave function methods has been rigorously established, spanning perturbation theory and the coupled cluster (CC) formalism. For these systems, CC with singles, doubles, and perturbative triples is known to achieve chemical accuracy, albeit at O(N7) computational cost. In addition, a hierarchy of density functional approximations of increasing formal sophistication, known as Jacob's ladder, has been shown to systematically reduce average errors over large datasets representing weakly correlated chemistry. However, the accuracy of such computational models is less clear in the increasingly important frontiers of chemical space including transition metals and f-block compounds, in which strong correlation can play an important role in reactivity. A stochastic method, phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC), has been shown to be capable of producing chemically accurate predictions even for challenging molecular systems beyond the main group, with relatively low O(N3 - N4) cost and near-perfect parallel efficiency. Herein, we present our perspectives on the past, present, and future of the ph-AFQMC method. We focus on its potential in transition metal quantum chemistry to be a highly accurate, systematically improvable method that can reliably probe strongly correlated systems in biology and chemical catalysis and provide reference thermochemical values (for future development of density functionals or interatomic potentials) when experiments are either noisy or absent. Finally, we discuss the present limitations of the method and where we expect near-term development to be most fruitful.

Pseudospectral implementations of long-range corrected density functional theory.

We have implemented pseudospectral density-functional theory (DFT) with long-range corrected DFT functionals (PS-LRC) in quantum mechanics package Jaguar, and applied it in the calculations of geometry optimizations, dimmer interaction energies, polarizabilities and first-order hyperpolarizabilities, harmonic vibrational frequencies, S1 and T1 excitation energies, singlet-triplet gaps, charge transfer numbers, oscillator strengths, reaction barrier heights, electron-transfer couplings, and charge-transfer excitation energies. From our accuracy benchmark analysis, PS grids, PS dealiasing functions, PS atomic corrections, PS multigrid strategy, PS length scales, and PS cutoff scheme perform well in PS DFT with LRC density functionals with very small and ignorable deviations when compared to the conventional spectral (CS) method. The timing benchmark study of S1 excitation energy calculations of fullerenes (Cn , n up to 540) demonstrates that PS-LRC achieves 1.4-8.4-fold speedups in SCF, 22-32-fold speedups in Tamm-Dancoff approximation, and 6-15-fold speedups in total wall clock time with an average error 0.004 eV of excitation energies compared to the CS method.

Highly efficient implementation of the analytical gradients of pseudospectral time-dependent density functional theory.

The accuracy and efficiency of time-dependent density functional theory (TDDFT) excited state gradient calculations using the pseudospectral method are presented. TDDFT excited state geometry optimizations of the G2 test set molecules, the organic fluorophores with large Stokes shifts, and the Pt-complexes show that the pseudospectral method gives average errors of 0.01-0.1 kcal/mol for the TDDFT excited state energy, 0.02-0.06 pm for the bond length, and 0.02-0.12° for the bond angle when compared to the results from conventional TDDFT. TDDFT gradient calculations of fullerenes (Cn, n up to 540) with the B3LYP functional and 6-31G** basis set show that the pseudospectral method provides 8- to 14-fold speedups in the total wall clock time over the conventional methods. The pseudospectral TDDFT gradient calculations with a diffuse basis set give higher speedups than the calculations for the same basis set without diffuse functions included.

Multiple Stable Isoprene-Ozone Complexes Reveal Complex Entrance Channel Dynamics in the Isoprene + Ozone Reaction.

Accurately characterizing isoprene ozonolysis continues to challenge atmospheric chemists. The reaction is believed to be a spontaneous, concerted cycloaddition. However, little information is available about the entrance channel and isoprene-ozone complexes thought to define the long-range portion of the reaction coordinate. Our coupled cluster and auxiliary field quantum Monte Carlo calculations predict multiple stable isoprene-ozone van der Waals complexes for trans-isoprene in the gas phase with moderate association energies. These results indicate that long-range dynamics in the isoprene-ozone entrance channel can impact the overall reaction in the troposphere and provide the spectroscopic information necessary to extend the microwave characterization of isoprene ozonolysis to prereactive complexes. At the air-water interface, Born-Oppenheimer molecular dynamics simulations indicate that the cycloaddition reaction between ozone and trans-isoprene follows a stepwise mechanism, which is quite distinct from our proposed gas-phase mechanism and occurs on a femtosecond time scale. The stepwise nature of isoprene ozonolysis on the aqueous surface is more consistent with the DeMore mechanism than with the Criegee mechanism suggested by the gas-phase calculations, suggesting that the reaction media may play an important role in the reaction. Overall, these predictions aim to provide a missing fundamental piece of molecular insight into isoprene ozonolysis, which has broad tropospheric implications due to its critical role as a nighttime source of hydroxyl radicals.

Advancing Drug Discovery through Enhanced Free Energy Calculations.

A principal goal of drug discovery project is to design molecules that can tightly and selectively bind to the target protein receptor. Accurate prediction of protein-ligand binding free energies is therefore of central importance in computational chemistry and computer aided drug design. Multiple recent improvements in computing power, classical force field accuracy, enhanced sampling methods, and simulation setup have enabled accurate and reliable calculations of protein-ligands binding free energies, and position free energy calculations to play a guiding role in small molecule drug discovery. In this Account, we outline the relevant methodological advances, including the REST2 (Replica Exchange with Solute Temperting) enhanced sampling, the incorporation of REST2 sampling with convential FEP (Free Energy Perturbation) through FEP/REST, the OPLS3 force field, and the advanced simulation setup that constitute our FEP+ approach, followed by the presentation of extensive comparisons with experiment, demonstrating sufficient accuracy in potency prediction (better than 1 kcal/mol) to substantially impact lead optimization campaigns. The limitations of the current FEP+ implementation and best practices in drug discovery applications are also discussed followed by the future methodology development plans to address those limitations. We then report results from a recent drug discovery project, in which several thousand FEP+ calculations were successfully deployed to simultaneously optimize potency, selectivity, and solubility, illustrating the power of the approach to solve challenging drug design problems. The capabilities of free energy calculations to accurately predict potency and selectivity have led to the advance of ongoing drug discovery projects, in challenging situations where alternative approaches would have great difficulties. The ability to effectively carry out projects evaluating tens of thousands, or hundreds of thousands, of proposed drug candidates, is potentially transformative in enabling hard to drug targets to be attacked, and in facilitating the development of superior compounds, in various dimensions, for a wide range of targets. More effective integration of FEP+ calculations into the drug discovery process will ensure that the results are deployed in an optimal fashion for yielding the best possible compounds entering the clinic; this is where the greatest payoff is in the exploitation of computer driven design capabilities. A key conclusion from the work described is the surprisingly robust and accurate results that are attainable within the conventional classical simulation, fixed charge paradigm. No doubt there are individual cases that would benefit from a more sophisticated energy model or dynamical treatment, and properties other than protein-ligand binding energies may be more sensitive to these approximations. We conclude that an inflection point in the ability of MD simulations to impact drug discovery has now been attained, due to the confluence of hardware and software development along with the formulation of "good enough" theoretical methods and models.

Highly efficient implementation of pseudospectral time-dependent density-functional theory for the calculation of excitation energies of large molecules.

We have developed and implemented pseudospectral time-dependent density-functional theory (TDDFT) in the quantum mechanics package Jaguar to calculate restricted singlet and restricted triplet, as well as unrestricted excitation energies with either full linear response (FLR) or the Tamm-Dancoff approximation (TDA) with the pseudospectral length scales, pseudospectral atomic corrections, and pseudospectral multigrid strategy included in the implementations to improve the chemical accuracy and to speed the pseudospectral calculations. The calculations based on pseudospectral time-dependent density-functional theory with full linear response (PS-FLR-TDDFT) and within the Tamm-Dancoff approximation (PS-TDA-TDDFT) for G2 set molecules using B3LYP/6-31G*(*) show mean and maximum absolute deviations of 0.0015 eV and 0.0081 eV, 0.0007 eV and 0.0064 eV, 0.0004 eV and 0.0022 eV for restricted singlet excitation energies, restricted triplet excitation energies, and unrestricted excitation energies, respectively; compared with the results calculated from the conventional spectral method. The application of PS-FLR-TDDFT to OLED molecules and organic dyes, as well as the comparisons for results calculated from PS-FLR-TDDFT and best estimations demonstrate that the accuracy of both PS-FLR-TDDFT and PS-TDA-TDDFT. Calculations for a set of medium-sized molecules, including Cn fullerenes and nanotubes, using the B3LYP functional and 6-31G(**) basis set show PS-TDA-TDDFT provides 19- to 34-fold speedups for Cn fullerenes with 450-1470 basis functions, 11- to 32-fold speedups for nanotubes with 660-3180 basis functions, and 9- to 16-fold speedups for organic molecules with 540-1340 basis functions compared to fully analytic calculations without sacrificing chemical accuracy. The calculations on a set of larger molecules, including the antibiotic drug Ramoplanin, the 46-residue crambin protein, fullerenes up to C540 and nanotubes up to 14×(6,6), using the B3LYP functional and 6-31G(**) basis set with up to 8100 basis functions show that PS-FLR-TDDFT CPU time scales as N(2.05) with the number of basis functions. © 2016 Wiley Periodicals, Inc.

Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field.

Designing tight-binding ligands is a primary objective of small-molecule drug discovery. Over the past few decades, free-energy calculations have benefited from improved force fields and sampling algorithms, as well as the advent of low-cost parallel computing. However, it has proven to be challenging to reliably achieve the level of accuracy that would be needed to guide lead optimization (∼5× in binding affinity) for a wide range of ligands and protein targets. Not surprisingly, widespread commercial application of free-energy simulations has been limited due to the lack of large-scale validation coupled with the technical challenges traditionally associated with running these types of calculations. Here, we report an approach that achieves an unprecedented level of accuracy across a broad range of target classes and ligands, with retrospective results encompassing 200 ligands and a wide variety of chemical perturbations, many of which involve significant changes in ligand chemical structures. In addition, we have applied the method in prospective drug discovery projects and found a significant improvement in the quality of the compounds synthesized that have been predicted to be potent. Compounds predicted to be potent by this approach have a substantial reduction in false positives relative to compounds synthesized on the basis of other computational or medicinal chemistry approaches. Furthermore, the results are consistent with those obtained from our retrospective studies, demonstrating the robustness and broad range of applicability of this approach, which can be used to drive decisions in lead optimization.

On achieving high accuracy and reliability in the calculation of relative protein-ligand binding affinities.

We apply a free energy perturbation simulation method, free energy perturbation/replica exchange with solute tempering, to two modifications of protein-ligand complexes that lead to significant conformational changes, the first in the protein and the second in the ligand. The approach is shown to facilitate sampling in these challenging cases where high free energy barriers separate the initial and final conformations and leads to superior convergence of the free energy as demonstrated both by consistency of the results (independence from the starting conformation) and agreement with experimental binding affinity data. The second case, consisting of two neutral thrombin ligands that are taken from a recent medicinal chemistry program for this interesting pharmaceutical target, is of particular significance in that it demonstrates that good results can be obtained for large, complex ligands, as opposed to relatively simple model systems. To achieve quantitative agreement with experiment in the thrombin case, a next generation force field, Optimized Potentials for Liquid Simulations 2.0, is required, which provides superior charges and torsional parameters as compared to earlier alternatives.

Covalent O-H bonds as electron traps in proton-rich rutile TiO2 nanoparticles.

The cation in the electrolyte of the dye-sensitized solar cell (DSSC) has a profound effect on electron trapping and transport behavior in TiO2 nanocrystalline film; this is one of the important factors that determines the overall efficiency of DSSCs. Here, we present a quantum mechanical investigation on the structures and energetics of proton-induced electron trap states and the thermodynamical barrier heights for the ambipolar diffusion of proton/electron pair using a large cluster model for the computations. Our calculations indicate that protons react with TiO2 to form covalent O-H bonds. This is in contrast to the reaction of Li(+) with TiO2, in which case the alkali metal is more accurately described as a simple coordinating cation. The covalent O-H bonding leads both to deeper electron trap states and to significantly higher barriers for the diffusion of carriers. These results are qualitatively consistent with experimental observations, and they extend our understanding of the cation effect in DSSCs at an atomic level of detail.

Multicolor live-cell chemical imaging by isotopically edited alkyne vibrational palette.

Vibrational imaging such as Raman microscopy is a powerful technique for visualizing a variety of molecules in live cells and tissues with chemical contrast. Going beyond the conventional label-free modality, recent advance of coupling alkyne vibrational tags with stimulated Raman scattering microscopy paves the way for imaging a wide spectrum of alkyne-labeled small biomolecules with superb sensitivity, specificity, resolution, biocompatibility, and minimal perturbation. Unfortunately, the currently available alkyne tag only processes a single vibrational "color", which prohibits multiplex chemical imaging of small molecules in a way that is being routinely practiced in fluorescence microscopy. Herein we develop a three-color vibrational palette of alkyne tags using a (13)C-based isotopic editing strategy. We first synthesized (13)C isotopologues of EdU, a DNA metabolic reporter, by using the newly developed alkyne cross-metathesis reaction. Consistent with theoretical predictions, the mono-(13)C ((13)C≡(12)C) and bis-(13)C ((13)C≡(13)C) labeled alkyne isotopologues display Raman peaks that are red-shifted and spectrally resolved from the originally unlabeled ((12)C≡(12)C) alkynyl probe. We further demonstrated three-color chemical imaging of nascent DNA, RNA, and newly uptaken fatty-acid in live mammalian cells with a simultaneous treatment of three different isotopically edited alkynyl metabolic reporters. The alkyne vibrational palette presented here thus opens up multicolor imaging of small biomolecules, enlightening a new dimension of chemical imaging.

A first-principles polarized Raman method for determining whether a uniform region of a sample is crystalline or isotropic.

Previous methods for determining whether a uniform region of a sample is crystalline or isotropic-what we call the "state of internal orientation" S-require a priori knowledge of properties of the purely crystalline and purely isotropic states. In addition, these methods can be ambiguous in their determination of state S for particular materials and, for a given material, the spectral methods can be ambiguous when using particular peaks. Using first-principles Raman theory, we have discovered a simple, non-resonance, polarized Raman method for determining the state S that requires no information a priori and will work unambiguously for any material using any vibrational mode. Similar to the concept behind "magic angle spinning" in NMR, we have found that for a special set of incident/analyzed polarizations and scattering angle, the dependence of the Raman modulation depth M on the sample composition-and, for crystalline regions, the unit cell orientation-falls out completely, leaving dependence on only whether the region is crystalline (M = 1) or isotropic (M = 0). Further, upon scanning between homogeneous regions or domains within a heterogeneous sample, our signal M is a clear detector of the region boundaries, so that when combined with methods for determining the orientations of the crystalline domains, our method can be used to completely characterize the molecular structure of an entire heterogeneous sample to a very high certainty. Interestingly, our method can also be used to determine when a given mode is vibrationally degenerate. While simulations on realistic terthiophene systems are included to illustrate our findings, our method should apply to any type of material, including thin films, molecular crystals, and semiconductors. Finally, our discovery of these relationships required derivations of Raman intensity formulas that are at least as general as any we have found, and herein we present our comprehensive formulas for both the crystalline and isotropic states.

Successful prediction of the intra- and extracellular loops of four G-protein-coupled receptors.

We present results of the restoration of all crystallographically available intra- and extracellular loops of four G-protein-coupled receptors (GPCRs): bovine rhodopsin (bRh), the turkey ß-1 adrenergic receptor (ß1Ar), and the human ß-2 adrenergic (ß2Ar) and A2A adenosine (A2Ar) receptors. We use our Protein Local Optimization Program (PLOP), which samples conformational space from first principles to build sets of loop candidates and then discriminates between them using our physics-based, all-atom energy function with implicit solvent. We also discuss a new kind of explicit membrane calculation developed for GPCR loops that interact, either in the native structure or in low-energy false-positive structures, with the membrane, and thus exist in a multiphase environment not previously incorporated in PLOP. Our results demonstrate a significant advance over previous work reported in the literature, and of particular note we are able to accurately restore the extremely long second extracellular loop (ECL2), which is also key for GPCR ligand binding. In the case of ß2Ar, accurate ECL2 restoration required seeding a small helix into the loop in the appropriate region, based on alignment with the ß1Ar ECL2 loop, and then running loop reconstruction simulations with and without the seeded helix present; simulations containing the helix attain significantly lower total energies than those without the helix, and have rmsds close to the native structure. For ß1Ar, the same protocol was used, except the alignment was done to ß2Ar. These results represent an encouraging start for the more difficult problem of accurate loop refinement for GPCR homology modeling.

Predicted and Experimental NMR Chemical Shifts at Variable Temperatures: The Effect of Protein Conformational Dynamics.

NMR chemical shifts provide a sensitive probe of protein structure and dynamics but remain challenging to predict and interpret. We examine the effect of protein conformational distributions on 15N chemical shifts for dihydrofolate reductase (DHFR), comparing QM/MM predicted shifts with experimental shifts in solution as well as frozen distributions. Representative snapshots from MD trajectories exhibit variation in predicted 15N chemical shifts of up to 25 ppm. The average over the fluctuations is in significantly better agreement with room temperature solution experimental values than the prediction for any single optimal conformations. Meanwhile, solid-state NMR (SSNMR) measurements of frozen solutions at 105 K exhibit broad lines whose widths agree well with the widths of distributions of predicted shifts for samples from the trajectory. The backbone torsion angle ψi-1 varies over 60° on the picosecond time scale, compensated by φi. These fluctuations can explain much of the shift variation.

Parsing Dynamics of Protein Backbone NH and Side-Chain Methyl Groups using Molecular Dynamics Simulations.

Experimental NMR spectroscopy and theoretical molecular dynamics (MD) simulations provide complementary insights into protein conformational dynamics and hence into biological function. The present work describes an extensive set of backbone NH and side-chain methyl group generalized order parameters for the Escherichia coli ribonuclease HI (RNH) enzyme derived from 2-µs microsecond MD simulations using the OPLS4 and AMBER-FF19SB force fields. The simulated generalized order parameters are compared with values derived from NMR 15N and 13CH2D spin relaxation measurements. The squares of the generalized order parameters, S2 for the N-H bond vector and Saxis2 for the methyl group symmetry axis, characterize the equilibrium distribution of vector orientations in a molecular frame of reference. Optimal agreement between simulated and experimental results was obtained by averaging S2 or Saxis2 calculated by dividing the simulated trajectories into 50 ns blocks (∼five times the rotational diffusion correlation time for RNH). With this procedure, the median absolute deviations (MAD) between experimental and simulated values of S2 and Saxis2 are 0.030 (NH) and 0.061 (CH3) for OPLS4 and 0.041 (NH) and 0.078 (CH3) for AMBER-FF19SB. The MAD between OPLS4 and AMBER-FF19SB are 0.021 (NH) and 0.072 (CH3). The generalized order parameters for the methyl group symmetry axis can be decomposed into contributions from backbone fluctuations, between-rotamer dihedral angle transitions, and within-rotamer dihedral angle fluctuations. Analysis of the simulation trajectories shows that (i) backbone and side chain conformational fluctuations exhibit little correlation and that (ii) fluctuations within rotamers are limited and highly uniform with values that depend on the number of dihedral angles considered. Low values of Saxis2, indicative of enhanced side-chain flexibility, result from between-rotamer transitions that can be enhanced by increased local backbone flexibility.

Using AlphaFold and Experimental Structures for the Prediction of the Structure and Binding Affinities of GPCR Complexes via Induced Fit Docking and Free Energy Perturbation.

Free energy perturbation (FEP) remains an indispensable method for computationally assaying prospective compounds in advance of synthesis. However, before FEP can be deployed prospectively, it must demonstrate retrospective recapitulation of known experimental data where the subtle details of the atomic ligand-receptor model are consequential. An open question is whether AlphaFold models can serve as useful initial models for FEP in the regime where there exists a congeneric series of known chemical matter but where no experimental structures are available either of the target or of close homologues. As AlphaFold structures are provided without a bound ligand, we employ induced fit docking to refine the AlphaFold models in the presence of one or more congeneric ligands. In this work, we first validate the performance of our latest induced fit docking technology, IFD-MD, on a retrospective set of public experimental GPCR structures with 95% of cross-docks producing a pose with a ligand RMSD ≤ 2.5 Å in the top two predictions. We then apply IFD-MD and FEP on AlphaFold models of the somatostatin receptor family of GPCRs. We use AlphaFold models produced prior to the availability of any experimental structure from this family. We arrive at FEP-validated models for SSTR2, SSTR4, and SSTR5, with RMSE around 1 kcal/mol, and explore the challenges of model validation under scenarios of limited ligand data, ample ligand data, and categorical data.

Robust prediction of relative binding energies for protein-protein complex mutations using free energy perturbation calculations.

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.