Atoms in molecules in real space: a fertile field for chemical bonding.

In this perspective, we review some recent advances in the concept of atoms-in-molecules from a real space perspective. We first introduce the general formalism of atomic weight factors that allows unifying the treatment of fuzzy and non-fuzzy decompositions under a common algebraic umbrella. We then show how the use of reduced density matrices and their cumulants allows partitioning any quantum mechanical observable into atomic or group contributions. This circumstance provides access to electron counting as well as energy partitioning, on the same footing. We focus on how the fluctuations of atomic populations, as measured by the statistical cumulants of the electron distribution functions, are related to general multi-center bonding descriptors. Then we turn our attention to the interacting quantum atom energy partitioning, which is briefly reviewed since several general accounts on it have already appeared in the literature. More attention is paid to recent applications to large systems. Finally, we consider how a common formalism to extract electron counts and energies can be used to establish an algebraic justification for the extensively used bond order-bond energy relationships. We also briefly review a path to recover one-electron functions from real space partitions. Although most of the applications considered will be restricted to real space atoms taken from the quantum theory of atoms in molecules, arguably the most successful of all the atomic partitions devised so far, all the take-home messages from this perspective are generalizable to any real space decompositions.

Toward Reliable and Insightful Entropy Calculations on Flexible Molecules.

The absolute entropy of a flexible molecule can be approximated by the sum of a rigid-rotor-harmonic-oscillator (RRHO) entropy and a Gibbs-Shannon entropy associated to the Boltzmann distribution for the occupation of the conformational energy levels. Herein, we show that such partitioning, which has received renewed interest, leads to accurate entropies of single molecules of increasing size provided that the conformational part is estimated by means of a set of discretization and expansion techniques that are able to capture the significant correlation effects among the torsional motions. To ensure a reliable entropy estimation, we rely on extensive sampling as that produced by classical molecular dynamics simulations on the microsecond time scale, which is currently affordable for small- and medium-sized molecules. According to test calculations, the gas-phase entropy of simple organic molecules is predicted with a mean unsigned error of 0.9 cal/(mol K) when the RRHO entropies are computed at the B3LYP-D3/cc-pVTZ level. Remarkably, the same protocol gives small errors [<1 cal/(mol K)] for the extremely flexible linear alkane molecules (CnH2n+2, n = 14, 16, and 18). Similarly, we obtain well-converged entropies for a more challenging test of drug molecules, which exhibit more pronounced correlation effects. We also perform equivalent entropy calculations on a 76 amino acid protein, ubiquitin, by taking advantage of the cutoff-dependent formulation of an expansion technique (correlation-consistent multibody local approximation, CC-MLA), which incorporates genuine correlation effects among the neighboring dihedral angles. Moreover, we show that insightful descriptors of the coupled torsional motions can be obtained with the CC-MLA approach.

QM/MM Energy Decomposition Using the Interacting Quantum Atoms Approach.

The interacting quantum atoms (IQA) method decomposes the quantum mechanical (QM) energy of a molecular system in terms of one- and two-center (atomic) contributions within the context of the quantum theory of atoms in molecules. Here, we demonstrate that IQA, enhanced with molecular mechanics (MM) and Poisson-Boltzmann surface-area (PBSA) solvation methods, is naturally extended to the realm of hybrid QM/MM methodologies, yielding intra- and inter-residue energy terms that characterize all kinds of covalent and noncovalent bonding interactions. To test the robustness of this approach, both metal-water interactions and QM/MM boundary artifacts are characterized in terms of the IQA descriptors derived from QM regions of varying size in Zn(II)- and Mg(II)-water clusters. In addition, we analyze a homologous series of inhibitors in complex with a matrix metalloproteinase (MMP-12) by carrying out QM/MM-PBSA calculations on their crystallographic structures followed by IQA energy decomposition. Overall, these applications not only show the advantages of the IQA QM/MM approach but also address some of the challenges lying ahead for expanding the QM/MM methodology.

A Quantum Chemical Topology Picture of Intermolecular Electrostatic Interactions and Charge Penetration Energy.

Based on the Interacting Quantum Atoms approach, we present herein a conceptual and theoretical framework of short-range electrostatic interactions, whose accurate description is still a challenging problem in molecular modeling. For all the noncovalent complexes in the S66 database, the fragment-based and atomic decomposition of the electrostatic binding energies is performed using both the charge density of the dimers and the unrelaxed densities of the monomers. This energy decomposition together with dispersion corrections gives rise to a pairwise approximation to the total binding energy. It also provides energetic descriptors at varying distance that directly address the atomic and molecular electrostatic interactions as described by point-charge or multipole-based potentials. Additionally, we propose a consistent definition of the charge penetration energy within quantum chemical topology, which is mainly characterized in terms of the intramolecular electrostatic energy. Finally, we discuss some practical implications of our results for the design and validation of electrostatic potentials.

Influence of charge configuration on substrate binding to SARS-CoV-2 main protease.

While the state-of-the-art computational simulations support the neutral state for the catalytic dyad of the SARS-CoV-2 main protease, the recently-reported neutron structure exhibits a zwitterionic form. To better compare the structural and dynamical features of the two charge configurations, we perform a Molecular Dynamics study of the dimeric enzyme in complex with a peptide substrate. The simulations show that the enzyme charge configuration from the neutron structure is not compatible with a catalytically-competent binding mode for peptide substrates.

Understanding the Conformational Properties of Fluorinated Polypeptides: Molecular Modelling of Unguisin A.

In this work, we investigate the conformational properties of unguisin A, a natural macrocyclic heptapeptide that incorporates a γ-aminobutyric acid (Gaba), and four of its difluorinated stereoisomers at the Gaba residue. According to nuclear magnetic resonance (NMR) experiments, their secondary structure depends dramatically on the stereochemistry of the fluorinated carbon atoms. However, many molecular details of the structure and flexibility of these systems remain unknown, so that a rationale of the conformational changes induced by the fluorine atoms in the macrocycle is still missing. To fill this gap, we apply enhanced molecular dynamics (MD) techniques to explore the peptide conformational space in dimethyl sulfoxide solution followed by 4-8 µs of conventional MD simulations that provide extensive equilibrium sampling. The simulations, which compare reasonably well with the NMR-based observations, show that the secondary structure of the macrocycle is altered substantially upon fluorination, except for the (S,S) diastereomer. It also turns out that the conformations of the fluorinated peptides are visited during the enhanced MD simulation of natural unguisin A, suggesting thus that conformations accessible to the unsubstituted macrocyclic peptide may be selected by fluorination. Therefore, computational characterization of the macrocyclic peptides could be helpful in the rational design of stereoselective fluorinated peptides with fine-tuned conformation and activity.

SARS-CoV-2 Main Protease: A Molecular Dynamics Study.

Herein, we investigate the structure and flexibility of the hydrated SARS-CoV-2 main protease by means of 2.0 µs molecular dynamics (MD) simulations in explicit solvent. After having performed electrostatic pKa calculations on several X-ray structures, we consider both the native (unbound) configuration of the enzyme and its noncovalent complex with a model peptide, Ace-Ala-Val-Leu-Gln∼Ser-Nme, which mimics the polyprotein sequence recognized at the active site. For each configuration, we also study their monomeric and homodimeric forms. The simulations of the unbound systems show that the relative orientation of domain III is not stable in the monomeric form and provide further details about interdomain motions, protomer-protomer interactions, inter-residue contacts, accessibility at the catalytic site, etc. In the presence of the peptide substrate, the monomeric protease exhibits a stable interdomain arrangement, but the relative orientation between the scissile peptide bond and the catalytic dyad is not favorable for catalysis. By means of comparative analysis, we further assess the catalytic impact of the enzyme dimerization, the actual flexibility of the active site region, and other structural effects induced by substrate binding. Overall, our computational results complement previous crystallographic studies on the SARS-CoV-2 enzyme and, together with other simulation studies, should contribute to outline useful structure-activity relationships.

Aptamers targeting protein-specific glycosylation in tumor biomarkers: general selection, characterization and structural modeling.

Detecting specific protein glycoforms is attracting particular attention due to its potential to improve the performance of current cancer biomarkers. Although natural receptors such as lectins and antibodies have served as powerful tools for the detection of protein-bound glycans, the development of effective receptors able to integrate in the recognition both the glycan and peptide moieties is still challenging. Here we report a method for selecting aptamers toward the glycosylation site of a protein. It allows identification of an aptamer that binds with nM affinity to prostate-specific antigen, discriminating it from proteins with a similar glycosylation pattern. We also computationally predict the structure of the selected aptamer and characterize its complex with the glycoprotein by docking and molecular dynamics calculations, further supporting the binary recognition event. This study opens a new route for the identification of aptamers for the binary recognition of glycoproteins, useful for diagnostic and therapeutic applications.

Alkali and Alkaline-Earth Cations in Complexes with Small Bioorganic Ligands: Ab†Initio Benchmark Calculations and Bond Energy Decomposition.

Herein, we report a computational database for the complexes of alkali [Li(I), Na(I), K(I)] and alkaline-earth [Be(II), Mg(II) and Ca(II)] cations with 25 small ligands with varying charge and donor atoms ("O", "N", and "S") that provides geometries and accurate bond energies useful to analyze metal-ligand interactions in proteins and nucleic acids. The role of the ligand→metal charge transfer, the equilibrium bond distance, the electronegativity of the donor atom, the ligand polarizability, and the relative stability of the complexes are discussed in detail. The interacting quantum atoms (IQA) method is used to decompose the binding energy into electrostatic and quantum mechanical contributions. In addition, bond energies are also estimated by means of multipolar electrostatic calculations. No simple correlation exists between bond energies and structural/electronic descriptors unless the data are segregated by the type of ligand or metal. The electrostatic attraction of some molecules (H2 O, NH3 , CH3 OH) towards the metal cations is well reproduced using their (unrelaxed) atomic multipoles, but the same comparison is much less satisfactory for other ligands (e. g. benzene, thiol/thiolate groups, etc.). Besides providing reference structures and bond energies, the database can contribute to validate molecular mechanics potentials capable of yielding a balanced description of alkali and alkaline-earth metals binding to biomolecules.

Fluorine conformational effects characterized by energy decomposition analysis.

Electrostatic and stereoelectronic effects associated with fluorine atoms can be exploited as conformational tools for the design of shape-controlled functional molecules. To gain further insight into the nature and strength of these effects, we use the Interacting Quantum Atoms (IQA) method augmented with the semiclassical pairwise dispersion potential to decompose the conformational energies of fluoro-substituted molecules into fragment-based energy contributions, which include deformation/distortion terms and the electrostatic, exchange-correlation and dispersion interactions. The studied molecules comprise various F-CH2-CH2-X and F-CH2-CO-X systems, as well as selected conformers of an α,ß-difluoro-γ-amino-acid derivative that is potentially useful for the design of shape-controlled bioactive amino acids and peptides. We identify the most relevant exchange-correlation and/or electrostatic interaction terms contributing to the stability of the various conformers, and we show that IQA can be used to assess the gauche/anti or trans/cis preferences in molecules with two or more rotatable bonds as well as to study the roles played by other concomitant effects (e.g., CH/OH/NHF contacts). For the α,ß-difluoro-γ-amino acid derivatives, our theoretical analysis indicates that the gauche/anti and trans/cis effects associated with fluorine bonds can be significantly attenuated by other specific intra-molecular contacts.

Affinity Calculations of Cyclodextrin Host-Guest Complexes: Assessment of Strengths and Weaknesses of End-Point Free Energy Methods.

The end-point methods like MM/PBSA or MM/GBSA estimate the free energy of a biomolecule by combining its molecular mechanics energy with solvation free energy and entropy terms. On the one hand, their performance largely depends on the particular system of interest, and despite numerous attempts to improve their reliability that have resulted in many variants, there is still no clear alternative to improve their accuracy. On the other hand, the relatively small cyclodextrin host-guest complexes, for which high-quality binding calorimetric data are usually available, are becoming reference models for testing the accuracy of free energy methods. In this work, we further assess the performance of various MM/PBSA-like approaches as applied to cyclodextrin complexes. To this end, we select a set of complexes between ß-cyclodextrin and 57 small organic molecules that has been previously studied with the binding energy distribution analysis method in combination with an implicit solvent model ( Wickstrom, L.; He, P.; Gallicchio, E.; Levy, R. M. J. Chem. Theory Comput. 2013 , 9 , 3136 - 3150 ). For each complex, a conventional 1.0 µs molecular dynamics simulation in explicit solvent is performed. Then we employ semiempirical quantum chemical calculations, several variants of the MM-PB(GB)SA methods, entropy estimations, etc., to assess the reliability of the end-point affinity calculations. The best end-point protocol in this study, which combines DFTB3 energies with entropy corrections, yields estimations of the binding free energies that still have substantial errors (RMSE = 2.2 kcal/mol), but it exhibits a good prediction capacity in terms of ligand ranking ( R2 = 0.66) that is close to or even better than that of rigorous free energy methodologies. Our results can be helpful to discriminate between the intrinsic limitations of the end-point methods and other sources of error, such as the underlying energy and continuum solvation methods.

Interacting Quantum Atoms Approach and Electrostatic Solvation Energy: Assessing Atomic and Group Solvation Contributions.

The interacting quantum atoms (IQA) method decomposes the total energy of a molecular system in terms of one- and two-center (atomic) contributions within the context of the quantum theory of atoms in molecules. Here we incorporate electrostatic continuum solvent effects into the IQA energy decomposition. To this end, the interaction between the solute electrostatic potential and the solvent screening charges as defined within the COSMO solvation model is now included in a new version of the PROMOLDEN code, allowing thus to apply IQA in combination with COSMO-quantum chemical methods as well as to partition the electrostatic solvation energy into effective atomic and group contributions. To test the robustness of this approach, we carry out COSMO-HF/aug-cc-pVTZ calculations followed by IQA calculations on more than 400 neutral and ionic solutes extracted from the MNSol database. The computational results reveal a detailed atomic mapping of the electrostatic solvation energy that is useful to assess to what extent the solvation energy can be decomposed into atomic and group contributions of various parts of a solute molecule, as generally assumed by empirical methodologies that estimate solvation energy and/or logP values.

Application of the Interacting Quantum Atoms Approach to the S66 and Ionic-Hydrogen-Bond Datasets for Noncovalent Interactions.

The interacting quantum atoms (IQA) method can assess, systematically and in great detail, the strength and physics of both covalent and noncovalent interactions. The lack of a pair density in density functional theory (DFT), which precludes the direct IQA decomposition of the characteristic exchange-correlation energy, has been recently overcome by means of a scaling technique, which can largely expand the applicability of the method. To better assess the utility of the augmented IQA methodology to derive quantum chemical decompositions at the atomic and molecular levels, we report the results of Hartree-Fock (HF) and DFT calculations on the complexes included in the S66 and the ionic H-bond databases of benchmark geometry and binding energies. For all structures, we perform single-point and geometry optimizations using HF and selected DFT methods with triple-ζ basis sets followed by full IQA calculations. Pairwise dispersion energies are accounted for by the D3 method. We analyze the goodness of the HF-D3 and DFT-D3 binding energies, the magnitude of numerical errors, the fragment and atomic distribution of formation energies, etc. It is shown that fragment-based IQA decomposes the formation energies in comparable terms to those of perturbative approaches and that the atomic IQA energies hold the promise of rigorously quantifying atomic and group energy contributions in larger biomolecular systems.

Ligand Strain and Entropic Effects on the Binding of Macrocyclic and Linear Inhibitors: Molecular Modeling of Penicillopepsin Complexes.

Using extensive molecular dynamics simulations, we investigate the structure and dynamics of the complexes formed between penicillopepsin and two peptidomimetic inhibitors: a linear compound, isovaleryl(P4)-valine(P3)-asparagine(P2)-leucine(P1)-phosphonate-phenylalanine(P1'), and its macrocylic analog that includes a methylene bridge between the Asn(P2) and Phe(P1') side chains. The macrocyclic inhibitor, which has a 420-fold stronger affinity than that of the acyclic one, has been considered to lower the entropic penalty for binding. To better understand this binding preference, the solution structure of the inhibitors is studied by molecular dynamics simulations. Subsequently, we assess the influence of the enzyme/inhibitor contacts, the enzyme-induced inhibitor strain, the variation of the ligand configurational entropy and the enzyme reorganization by combining molecular-mechanics Poisson-Boltzmann surface area and normal mode calculations with conformational entropy calculations. We find that there is no relevant entropic stabilization on the binding of the cyclic inhibitor with respect to the acyclic analog because the methylene bridge does not reduce appreciably the conformational flexibility of the free inhibitor. The most important factors explaining the stronger affinity of the macrocyclic inhibitor are the conformational filtering and the lower ligand strain induced by the methylene bridge.

Molecular Dynamics Studies of Matrix Metalloproteases.

Matrix metalloproteases are multidomain enzymes with a remarkable proteolytic activity located in the extracellular environment. Their catalytic activity and structural properties have been intensively studied during the last few decades using both experimental and theoretical approaches, but many open questions still remain. Extensive molecular dynamics simulations enable the sampling of the configurational space of a molecular system, thus contributing to the characterization of the structure, dynamics, and ligand binding properties of a particular MMP. Based on previous computational experience, we provide in this chapter technical and methodological guidelines that may be useful to and stimulate other researchers to perform molecular dynamics simulations to help address unresolved questions concerning the molecular mode of action of MMPs.

Conformational and entropy analyses of extended molecular dynamics simulations of α-, ß- and γ-cyclodextrins and of the ß-cyclodextrin/nabumetone complex.

Herein, we report the results of 5.0 µs molecular dynamics simulations of native α-, ß- and γ-cyclodextrins (CDs) in explicit water solvent that are useful to describe, in a comparative manner, the distorted geometry of the CD molecules in aqueous solution, the width and fluctuations of their cavities, and the number of cavity waters. By discretizing the time evolution of the dihedral angles, the rate of conformational change of the torsional motions and the conformational entropy are calculated for the three CDs, thus allowing the analysis of the extent of the MD sampling and the entropic significance of the CD flexibility. To obtain a first estimation of the conformational and entropy changes in the host molecule upon ligand binding, the inclusion complex formed between ß-CD and nabumetone is also studied. Overall, the simulations complement previous experimental results on the structure and dynamics of native CDs, and together with the results obtained for the inclusion complex, provide insight into the entropic effects at work on the binding equilibria between CDs and guest ligands.

Role of the Protonation State on the Structure and Dynamics of Albumin.

Human serum albumin undergoes reversible conformational transitions associated with ligand binding or pH changes. Among them, the neutral to base (N → B) transition occurring between pH 7 and pH 9 seems to be relevant for its function as a carrier. Unfortunately, a detailed atomic model for the B-form is still lacking, and several open questions remain concerning the charge distribution of the N-form. In this work, we report comparable molecular models for the N and B conformations that are built using continuum electrostatic calculations of pKa values and extended molecular dynamics (MD) simulations. Our computational models, which are critically assessed in terms of the available experimental observations relative to the N → B transition, reveal interesting similarities and differences between the N- and B-forms of HSA and highlight the importance of setting proper charge configurations in MD simulations.

Unraveling the distinctive features of hemorrhagic and non-hemorrhagic snake venom metalloproteinases using molecular simulations.

Snake venom metalloproteinases are important toxins that play fundamental roles during envenomation. They share a structurally similar catalytic domain, but with diverse hemorrhagic capabilities. To understand the structural basis for this difference, we build and compare two dynamical models, one for the hemorrhagic atroxlysin-I from Bothrops atrox and the other for the non-hemorraghic leucurolysin-a from Bothrops leucurus. The analysis of the extended molecular dynamics simulations shows some changes in the local structure, flexibility and surface determinants that can contribute to explain the different hemorrhagic activity of the two enzymes. In agreement with previous results, the long Ω-loop (from residue 149 to 177) has a larger mobility in the hemorrhagic protein. In addition, we find some potentially-relevant differences at the base of the S1' pocket, what may be interesting for the structure-based design of new anti-venom agents. However, the sharpest differences in the computational models of atroxlysin-I and leucurolysin-a are observed in the surface electrostatic potential around the active site region, suggesting thus that the hemorrhagic versus non-hemorrhagic activity is probably determined by protein surface determinants.

Molecular modeling of bioorganometallic compounds: thermodynamic properties of molybdocene-glutathione complexes and mechanism of Peptide hydrolysis.

The computational study of bioinorganic complexes between transition metals and flexible ligands is still challenging, given that, besides requiring extensive conformational searches, the treatment of metal-ligand bonds demands the application of quantum chemical methods. Herein, the adducts formed between molybdocene, which exhibits antitumor activity and reacts with thiol groups to give stable water-soluble complexes, and the tripeptide glutathione, which is a major source of biological thiols, are studied. Conformational searches are performed using the semiempirical PM6 method followed by geometry optimizations and single-point calculations using density functional theory methods. In addition, molecular dynamics simulations of the molybdocene-glutathione complex involved in the regioselective hydrolysis of the Cys-Gly linkage are performed in explicit solvent. The reactive process is also studied theoretically on cluster models of both the molybdocene-bound and the free peptide.

Extensive simulations of the full-length matrix metalloproteinase-2 enzyme in a prereactive complex with a collagen triple-helical peptide.

Collagen hydrolysis catalyzed by matrix metalloproteinases is an important and complex process involved in a variety of physiological and pathological conditions. To contribute to its characterization at the molecular level, herein we analyze three different models for the complex formed between the full-length matrix metalloproteinase-2 (MMP-2) enzyme and a synthetic triple-helical peptide (fTHP-5). The considered MMP-2/fTHP-5 complexes mainly differ in the location of the C-terminal hemopexin-like domain, but in all of them, the middle α-chain of the substrate (B-chain) is placed within the active site groove. We performed extended molecular dynamics (MD) simulations to determine the most likely rearrangements of the MMP-2 domains in response to the presence of the triple helix. The relative stability of the MD models is assessed in terms of molecular mechanics Poisson-Boltzmann calculations and approximate estimations of configurational entropy. In addition, the most significant MMP-2···fTHP-5 interactions at the catalytic and noncatalytic domains are also analyzed to gather some clues about the role of the different domains during collagenolysis.