Pandemic drugs at pandemic speed: infrastructure for accelerating COVID-19 drug discovery with hybrid machine learning- and physics-based simulations on high-performance computers.

The race to meet the challenges of the global pandemic has served as a reminder that the existing drug discovery process is expensive, inefficient and slow. There is a major bottleneck screening the vast number of potential small molecules to shortlist lead compounds for antiviral drug development. New opportunities to accelerate drug discovery lie at the interface between machine learning methods, in this case, developed for linear accelerators, and physics-based methods. The two in silico methods, each have their own advantages and limitations which, interestingly, complement each other. Here, we present an innovative infrastructural development that combines both approaches to accelerate drug discovery. The scale of the potential resulting workflow is such that it is dependent on supercomputing to achieve extremely high throughput. We have demonstrated the viability of this workflow for the study of inhibitors for four COVID-19 target proteins and our ability to perform the required large-scale calculations to identify lead antiviral compounds through repurposing on a variety of supercomputers.

Simulated Solute Tempering in Fully Polarizable Hybrid QM/MM Molecular Dynamics Simulations.

We successfully apply a solute tempering approach, which substantially reduces the large number of temperature rungs required in conventional tempering methods by solvent charge scaling, to hybrid molecular dynamics simulations combining quantum mechanics with molecular mechanics (QM/MM). Specifically, we integrate a combination of density functional theory (DFT) and polarizable MM (PMM) force fields into the simulated solute tempering (SST) concept. We show that the required DFT/PMM-SST weight parameters can be obtained from inexpensive calculations and that for alanine dipeptide (DFT) in PMM water three rungs suffice to cover the temperature range from 300 to 550 K.

Continuous Tempering Molecular Dynamics: A Deterministic Approach to Simulated Tempering.

Continuous tempering molecular dynamics (CTMD) generalizes simulated tempering (ST) to a continuous temperature space. Opposed to ST the CTMD equations of motion are fully deterministic and feature a conserved quantity that can be used to validate the simulation. Three variants of CTMD are discussed and compared by means of a simple test system. The implementation features of the most stable and simplest variant CTMD-U in the program package Iphigenie are described. Two applications--alanine dipeptide (Ac-Ala-NHMe) in explicit water and octa-alanine (Ac-(Ala)8-NHMe) simulated in a dielectric continuum--demonstrate the functionality of CTMD-U. Furthermore, they serve to evaluate its sampling efficiency. Here, CTMD-U outperforms ST by 35% and replica exchange even by 75%.

Linearly scaling and almost Hamiltonian dielectric continuum molecular dynamics simulations through fast multipole expansions.

Hamiltonian Dielectric Solvent (HADES) is a recent method [S. Bauer et al., J. Chem. Phys. 140, 104103 (2014)] which enables atomistic Hamiltonian molecular dynamics (MD) simulations of peptides and proteins in dielectric solvent continua. Such simulations become rapidly impractical for large proteins, because the computational effort of HADES scales quadratically with the number N of atoms. If one tries to achieve linear scaling by applying a fast multipole method (FMM) to the computation of the HADES electrostatics, the Hamiltonian character (conservation of total energy, linear, and angular momenta) may get lost. Here, we show that the Hamiltonian character of HADES can be almost completely preserved, if the structure-adapted fast multipole method (SAMM) as recently redesigned by Lorenzen et al. [J. Chem. Theory Comput. 10, 3244-3259 (2014)] is suitably extended and is chosen as the FMM module. By this extension, the HADES/SAMM forces become exact gradients of the HADES/SAMM energy. Their translational and rotational invariance then guarantees (within the limits of numerical accuracy) the exact conservation of the linear and angular momenta. Also, the total energy is essentially conserved-up to residual algorithmic noise, which is caused by the periodically repeated SAMM interaction list updates. These updates entail very small temporal discontinuities of the force description, because the employed SAMM approximations represent deliberately balanced compromises between accuracy and efficiency. The energy-gradient corrected version of SAMM can also be applied, of course, to MD simulations of all-atom solvent-solute systems enclosed by periodic boundary conditions. However, as we demonstrate in passing, this choice does not offer any serious advantages.

Parameterization of the Hamiltonian Dielectric Solvent (HADES) Reaction-Field Method for the Solvation Free Energies of Amino Acid Side-Chain Analogs.

Optimization of the Hamiltonian dielectric solvent (HADES) method for biomolecular simulations in a dielectric continuum is presented with the goal of calculating accurate absolute solvation free energies while retaining the model's accuracy in predicting conformational free-energy differences. The solvation free energies of neutral and polar amino acid side-chain analogs calculated by using HADES, which may optionally include nonpolar contributions, were optimized against experimental data to reach a chemical accuracy of about 0.5 kcal mol(-1). The new parameters were evaluated for charged side-chain analogs. The HADES results were compared with explicit-solvent, generalized Born, Poisson-Boltzmann, and QM-based methods. The potentials of mean force (PMFs) between pairs of side-chain analogs obtained by using HADES and explicit-solvent simulations were used to evaluate the effects of the improved parameters optimized for solvation free energies on intermolecular potentials.

Utilizing fast multipole expansions for efficient and accurate quantum-classical molecular dynamics simulations.

Recently, a novel approach to hybrid quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations has been suggested [Schwörer et al., J. Chem. Phys. 138, 244103 (2013)]. Here, the forces acting on the atoms are calculated by grid-based density functional theory (DFT) for a solute molecule and by a polarizable molecular mechanics (PMM) force field for a large solvent environment composed of several 10(3)-10(5) molecules as negative gradients of a DFT/PMM hybrid Hamiltonian. The electrostatic interactions are efficiently described by a hierarchical fast multipole method (FMM). Adopting recent progress of this FMM technique [Lorenzen et al., J. Chem. Theory Comput. 10, 3244 (2014)], which particularly entails a strictly linear scaling of the computational effort with the system size, and adapting this revised FMM approach to the computation of the interactions between the DFT and PMM fragments of a simulation system, here, we show how one can further enhance the efficiency and accuracy of such DFT/PMM-MD simulations. The resulting gain of total performance, as measured for alanine dipeptide (DFT) embedded in water (PMM) by the product of the gains in efficiency and accuracy, amounts to about one order of magnitude. We also demonstrate that the jointly parallelized implementation of the DFT and PMM-MD parts of the computation enables the efficient use of high-performance computing systems. The associated software is available online.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description.

We present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys. 118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys. 1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic anti-polarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian anti-polarization densities are calculated by a self-consistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another self-consistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ(i) of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys. 2, 351 (1934)] and with computationally expensive grid-based numerical solutions for simple model systems in dielectric continua including a di-peptide (Ac-Ala-NHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ(i). A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with so-called generalized Born methods. A follow-up paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics.

In Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys. 140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call "Hamiltonian dielectric solvent" (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be self-consistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the self-consistency iterations are accelerated by the use of the so-called direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADES-MD, and compare in a sample application on Ac-Ala-NHMe the HADES-MD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

Coupling density functional theory to polarizable force fields for efficient and accurate Hamiltonian molecular dynamics simulations.

Hybrid molecular dynamics (MD) simulations, in which the forces acting on the atoms are calculated by grid-based density functional theory (DFT) for a solute molecule and by a polarizable molecular mechanics (PMM) force field for a large solvent environment composed of several 10(3)-10(5) molecules, pose a challenge. A corresponding computational approach should guarantee energy conservation, exclude artificial distortions of the electron density at the interface between the DFT and PMM fragments, and should treat the long-range electrostatic interactions within the hybrid simulation system in a linearly scaling fashion. Here we describe a corresponding Hamiltonian DFT/(P)MM implementation, which accounts for inducible atomic dipoles of a PMM environment in a joint DFT/PMM self-consistency iteration. The long-range parts of the electrostatics are treated by hierarchically nested fast multipole expansions up to a maximum distance dictated by the minimum image convention of toroidal boundary conditions and, beyond that distance, by a reaction field approach such that the computation scales linearly with the number of PMM atoms. Short-range over-polarization artifacts are excluded by using Gaussian inducible dipoles throughout the system and Gaussian partial charges in the PMM region close to the DFT fragment. The Hamiltonian character, the stability, and efficiency of the implementation are investigated by hybrid DFT/PMM-MD simulations treating one molecule of the water dimer and of bulk water by DFT and the respective remainder by PMM.

Mechanistic insights into the hydrolysis of a nucleoside triphosphate model in neutral and acidic solution.

Nucleoside triphosphate hydrolysis is an essential component of all living systems. Despite extensive research, the exact modus and mechanism of this ubiquitous reaction still remain elusive. In this work, we examined the detailed hydrolysis mechanisms of a model nucleoside triphosphate in acidic and neutral solution by means of ab initio simulations. The timescale of the reaction was accessed through use of an accelerated sampling method, metadynamics. Both hydrolyses were found to proceed via different mechanisms; the acidic system reacted by means of concerted general acid catalysis (found to be a so-called D(N)A(N)A(H)D(xh) mechanism), whereas the neutral system reacted by way of a different mechanism (namely, D(N)*A(N)D(xh)A(H)). A neighboring water molecule took on the role of a general base in both systems, which has not been seen before but is a highly plausible reaction path, meaning that substrate-assisted catalysis was not observed in the bulk water environment.

Vibrational spectra of phosphate ions in aqueous solution probed by first-principles molecular dynamics.

We have carried out "first-principles" Born-Oppenheimer molecular dynamics (BOMD) simulations of the phosphate ions H2PO4⁻ and HPO4²â» in liquid water and have calculated their IR spectra by Fourier transform techniques from the trajectories. IR bands were assigned by a so-called "generalized normal coordinate analysis". The effects of including Hartree-Fock (HF) exchange into the density functional theory (DFT) computation of forces were studied by comparing results obtained with the well-known BP, BLYP, and B3LYP functionals. The neglect of dispersion in the functionals was empirically corrected. The inclusion of HF exchange turned out to yield dramatically improved and, thus, quite accurate descriptions of the IR spectra observed for H2PO4⁻ and HPO4²â» in aqueous solution. An analysis of earlier computational results (Klähn, M. et al. J. Phys. Chem. A 2004, 108, 6186-6194) on these vibrational spectra, which had been obtained in a hybrid setting combining a BP description of the respective phosphate with a simple molecular mechanics (MM) model of its aqueous environment, revealed three different sources of error, (i) the BP force field of the phosphates is much too soft and would have required a substantial scaling of frequencies, (ii) the oversimplified water force field entailed incorrect solvation structures and, thus, qualitatively wrong patterns of solvatochromic band shifts, and (iii) quantitative frequency computations additionally required the inclusion of HF exchange. Thus, the results of the B3LYP BOMD simulations do not only characterize physical properties like the IR spectra or the solvation structures of the phosphate systems but also provide clues for the future design of simplified but nevertheless reasonably accurate DFT/MM methods applicable to phosphates.

Understanding the Origins of Dipolar Couplings and Correlated Motion in the Vibrational Spectrum of Water.

The combination of vibrational spectroscopy and molecular dynamics simulations provides a powerful tool to obtain insights into the molecular details of water structure and dynamics in the bulk and in aqueous solutions. Applying newly developed approaches to analyze correlations of charge currents, molecular dipole fluctuations, and vibrational motion in real and k-space, we compare results from nonpolarizable water models, widely used in biomolecular modeling, to ab initio molecular dynamics. For the first time, we unfold the infrared response of bulk water into contributions from correlated fluctuations in the three-dimensional, anisotropic environment of an average water molecule, from the OH-stretching region down to the THz regime. Our findings show that the absence of electronic polarizability in the force field model not only results in differences in dipolar couplings and infrared absorption but also induces artifacts into the correlated vibrational motion between hydrogen-bonded water molecules, specifically at the intramolecular bending frequency. Consequently, vibrational motion is partially ill-described with implications for the accuracy of non-self-consistent, a posteriori methods to add polarizability.

DFT/MM description of flavin IR spectra in BLUF domains.

A class of photoreceptors occurring in various organisms consists of domains that are blue light sensing using flavin (BLUF). The vibrational spectra of the flavin chromophore are spectroscopically well characterized for the dark-adapted resting states and for the light-adapted signaling states of BLUF domains in solution. Here we present a theoretical analysis of such spectra by applying density functional theory (DFT) to the flavin embedded in molecular mechanics (MM) models of its protein and solvent environment. By DFT/MM we calculate flavin spectra for seven different X-ray and NMR structures of BLUF domains occurring in the transcriptional antirepressor AppA and in the blue light receptor B (BlrB) of the purple bacterium Rb. Sphaeroides as well as in the phototaxis photoreceptor Slr1694 of the cyanobacterium Synechocystis. By considering the dynamical stabilities of associated all-atom simulation models and by comparing calculated with observed vibrational spectra, we show that two of the considered structures (both AppA) are obviously erroneous and that specific features of two further crystal structures (BlrB and Slr1694) cannot represent the states of the respective BLUF domains in solution. Thereby, the conformational transitions elicited by solvation are identified. In this context we demonstrate how hydrogen bonds of varying strengths can tune in BLUF domains the C═O stretching frequencies of the flavin chromophore. Furthermore we show that the DFT/MM spectra of the flavin calculated for two different AppA BLUF conformations, which are called Trp(in) and Met(in), fit very well to the spectroscopic data observed for the dark and light states, respectively, if (i) polarized MM force fields, which are calculated by an iterative DFT/MM procedure, are employed for the flavin binding pockets and (ii) the calculated frequencies are properly scaled. Although the associated analysis indicates that the Trp(in) conformation belongs to the dark state, no clear light vs dark distinction emerges for the Met(in) conformation. In this connection, a number of methodological issues relevant for such complex computations are thoroughly discussed showing, in particular, why our current descriptions could not decide the light vs dark question for Met(in).

Assigning predissociation infrared spectra of microsolvated hydronium cations H3O(+)⋠(H2)n (n=0, 1, 2, 3) by ab initio molecular dynamics.

Messenger predissociation spectroscopy is an important experimental method to obtain vibrational spectra of molecular ions or complexes such as protonated water clusters H(+)⋅(H(2)O)(n) in the gas phase. However, the molecular properties and thus the linear infrared spectra may be modified upon microsolvation with typical messengers such as H(2) molecules or noble gas atoms. Employing ab initio molecular dynamics for the H(2)-microsolvated hydronium ion, we investigate these effects explicitly as a function of an increasing number of messengers up to filling the first microsolvation shell, that is, for H(3)O(+)⋅(H(2))(n) (n=0, 1, 2, 3). We find that microsolvation with H(2) lowers the inversion barrier of the hydronium core, which governs the inversion tunnel splitting due to umbrella motion, and thus accelerates the inversion dynamics. By comparison to experiment a comprehensive band assignment for the O-H stretching region is given, and thereby the observed blueshift of stretching bands with increasing n is explained. Furthermore, detailed analyses reveal intricate intra- and intermolecular anharmonic mode couplings induced by the messengers, which yield a rich vibrational dynamics in these, at first glance, simple systems. Finally, the virtues but also the shortcomings of the ab initio molecular dynamics approach to vibrational spectroscopy are discussed.

IR spectra of flavins in solution: DFT/MM description of redox effects.

The functional reactions in blue light photoreceptors generally involve transiently reduced flavins exhibiting characteristic infrared (IR) spectra. To approach a theoretical understanding, here we apply density functional theory (DFT) to flavin radicals embedded in a molecular mechanics (MM) model of an aqueous solution. Combining a DFT/MM approach with instantaneous normal-mode analyses (INMA), we compute the IR solution spectra of anionic and neutral flavin radicals. For a set of mid-IR marker bands, we identify those changes of spectral locations, intensities, and widths, which are caused by sequentially adding an electron and a proton to the oxidized flavin. Comparisons with experimental IR solution spectra of flavin radicals show the accuracy of our DFT/MM-INMA approach and allow us to assign the observed bands. The room temperature ensembles of solvent cages required for the INMA calculations of the IR spectra are generated in an MM setting from molecular dynamics (MD) simulations. For the solvated flavin radicals, these MD simulations employ MM force fields derived from DFT/MM calculations.

Generalized Normal Coordinates for the Vibrational Analysis of Molecular Dynamics Simulations.

The computation of vibrational spectra via molecular dynamics (MD) simulations has made lively progress in recent years. In particular, infrared spectra are accessible employing ab initio MD, for which only the total dipole moment has to be computed "on the fly" from the electronic structure along the trajectory. The analysis of such spectra in terms of the normal modes of intramolecular motion, however, still poses a challenge to theory. Here, we present an algorithm to extract such normal modes from MD trajectories by combining several ideas available in the literature. The algorithm allows one to compute both the normal modes and their vibrational bands without having to rely on an equipartition assumption, which hampered previous methods. Our analysis is based on a tensorial definition of the vibrational density of states, which spans both the frequency resolved cross- and auto-correlations of the molecular degrees of freedom. Generalized normal coordinates are introduced as orthonormal transforms of mass-weighted coordinates, which minimize their mutual cross-correlations. The generalized normal coordinates and their associated normal modes are iteratively constructed by a minimization scheme based on the Jacobi diagonalization. Furthermore, the analysis furnishes mode local temperatures, which provide not only a measure for the convergence of the computed intensities but also permits one to correct these intensities a posteriori toward the ensemble limit. As a first non-trivial test application we analyze the infrared spectrum of isoprene based on ab initio MD, which is an important building block of various dye molecules in molecular biology.

Density functional theory combined with molecular mechanics: the infrared spectra of flavin in solution.

The photophysics and photochemistry of flavin dyes determine the functional dynamics of a series of blue light photoreceptors that include the so-called BLUF (blue light sensors using flavin) domains. To enable molecular dynamics (MD) simulation studies of such signaling processes, we derived molecular mechanics (MM) models of flavin chromophores from density functional theory (DFT). Two 300 K ensembles of lumiflavin (LF) in aqueous solution were generated by extended MM-MD simulations using different MM potentials for the water. In a DFT/MM hybrid setting, in which LF was treated by DFT and the polarizing environment at atomistic resolution by MM, we applied instantaneous normal mode analyses (INMA) to these ensembles. From these data we determined the inhomogeneously broadened solution spectra as mixtures of Gaussian bands using a novel automated procedure for mode classification. Comparisons with vibrational spectra available in the literature on native and isotopically labeled flavins in aqueous solution serve us to determine suitable frequency scaling factors and to analyze the accuracy of our scaled DFT/MM-INMA approach. We show that our approach not only agrees with established computational descriptions but also extends such methods substantially by giving access to inhomogeneous line widths and band shapes.

Quantum-induced symmetry breaking explains infrared spectra of CH(5)(+) isotopologues.

For decades, protonated methane, CH(5)(+), has provided new surprises and challenges for both experimentalists and theoreticians. This is because of the correlated large-amplitude motion of its five protons around the carbon nucleus, which leads to so-called hydrogen scrambling and causes a fluxional molecular structure. Here, the infrared spectra of all its H/D isotopologues have been measured using the 'Laser Induced Reactions' technique. Their shapes are found to be extremely dissimilar and depend strongly on the level of deuteration (only CD(5)(+) is similar to CH(5)(+)). All the spectra can be reproduced and assigned based on ab initio quantum simulations. The occupation of the topologically different sites by protons and deuterons is found to be strongly non-combinatorial and thus non-classical. This purely quantum-statistical effect implies a breaking of the classical symmetry of the site occupations induced by zero-point fluctuations, and this phenomenon is key to understanding the spectral changes studied here.