Electron density learning of non-covalent systems
Chemists continuously harvest the power of non-covalent interactions to control phenomena in both the micro-and macroscopic worlds. From the quantum chemical perspective, the strategies essentially rely upon an in-depth understanding of the physical origin of these interactions, the quantification of their magnitude and their visualization in real-space. The total electron density r(r) represents the simplest yet most comprehensive piece of information available for fully characterizing bonding patterns and noncovalent interactions. The charge density of a molecule can be computed by solving the Schr ödinger equation, but this approach becomes rapidly demanding if the electron density has to be evaluated for thousands of different molecules or very large chemical systems, such as peptides and proteins. Here we present a transferable and scalable machine-learning model capable of predicting the total electron density directly from the atomic coordinates. The regression model is used to access qualitative and quantitative insights beyond the underlying r(r) in a diverse ensemble of sidechain-sidechain dimers extracted from the BioFragment database (BFDb). The transferability of the model to more complex chemical systems is demonstrated by predicting and analyzing the electron density of a collection of 8 polypeptides.
Introduction<!>Methods<!>Computational details<!>Results and discussion<!>Prediction on polypeptides<!>Conclusion<!>Conflicts of interest
<p>Non-covalent interactions (NCIs) govern a multitude of chemical phenomena and are key components for constructing molecular architectures. 1 Their importance fostered an intense research effort to accurately quantify their magnitude and develop an intuitive characterization of their physical nature using quantum chemistry. [2][3][4][5][6] Among the different approaches to characterize non-covalent interactions, one of the simplest and most generally applicable takes as a starting point the electron density r(r) that encodes, in principle, all the information needed to fully characterize a chemical system. 7 Despite the fact that the universal functional relationship between total energy and r(r) remains unknown, existing approximations within the framework of Kohn-Sham DFT (KS-DFT) 8 do permit access to all molecular properties within a reasonable degree of accuracy. [9][10][11] Properties that can be derived exactly from the electron density distribution include molecular and atomic electrostatic moments (e.g., charges, dipole, quadrupoles), electrostatic potentials and electrostatic interaction energies. Knowledge of these quantities is fundamental in diverse chemical applications, including the computation of the IR intensities, 12 the identication of binding sites in host-guest compounds, [13][14][15] and the exact treatment of electrostatics within molecular simulations. 16 Moreover, analyzing the deformation of r(r) in the presence of an external eld provides access to another set of fundamental properties, namely molecular static (hyper) polarizabilities and, thus, to the computation of Raman spectra 17 and non-linear optical properties. [18][19][20][21] The natural representation of the electron density in real space makes it especially suitable for accessing spatial information about structural and electronic molecular properties, including X-ray structure renement [22][23][24][25][26][27] and representations using scalar elds. 6 Routinely used examples include the quantum theory of atoms in molecules (QTAIM), 28,29 the density overlap region indicator (DORI), 30 and the non-covalent interaction (NCI) index. 31,32 r(r) is generally obtained by solving the electronic structure problem through ab initio computations. The main advantage of this approach is that it returns the variationally optimized electronic density for a given Hamiltonian. Yet, ab initio computations can become increasingly burdensome if r(r) has to be evaluated for thousands of different molecules or very large chemical systems, such as peptides and proteins. These large scale problems are typically tackled using a more scalable approach that consists of either using linear scaling techniques such as Mezey's molecular electron density LEGO assembler (MEDLA) 33,34 and adjustable density matrix assembler (ADMA), [35][36][37] as well as approaches based on localized molecular orbitals, such as ELMO. [38][39][40][41] Another methodology belonging to this second category involves the use of experimental techniques, such as X-ray diffraction, to probe the electron density and subsequently reconstructing r(r) through multipolar models [42][43][44] and pseudo-atomic libraries, such as ELMAM, [45][46][47][48] ELMAM2, 49,50 UBDB, 51,52 Invarioms 53 and SBFA. 54 While successful, these two methodologies have intrinsic limits: the rst is unable to capture the deformations of the charge density due to intermolecular interactions unless a suitable fragment is generated ad hoc, while the second relies on experimental data and is difficult to extend to thousands of different chemical systems at once. Recently, the development of several machinelearning models targeting the electron density has effectively established a third promising methodology, with the potential to overcome the limitations of the more traditional approaches.</p><p>The rst machine-learning model of r(r) was developed on the basis of the Hohenberg-Kohn mapping between the nuclear potential and the electron density. 55,56 Although successful, the choice of the nuclear potential as a representation of the different molecular conformations and the expansion of the electron density in an orthogonal plane-wave basis effectively constrained this landmark model to relatively small and rigid molecules with limited transferability to larger systems. Recently, we proposed an atom-centered, symmetry-adapted Gaussian process regression 57 (SA-GPR) framework explicitly targeting the learning of the electron density. 58 Using an optimized non-orthogonal basis set, pseudo-valence electron densities could be predicted in a linear-scaling and transferable manner, meaning that the model is able to tackle much larger chemical systems than those used to train the regression model. A third approach, that can also achieve transferability between different systems, uses a direct grid-based representation of the atomic environment to learn and predict the electronic density in each point of the molecular space. [59][60][61] Representing the density eld on a large set of grids points rather than on a basis set effectively avoids the introduction of a basis set error, but also dramatically increases the computational effort.</p><p>One should also consider that machine learning, being a data-driven approach, requires high-quality, diverse reference data. Fortunately, several specialized benchmark databases that target NCIs have appeared over the past decade. From the original S22 (ref. 62) to NCIE53, 63 S66, 64 NBC10/NBC10ext, [65][66][67] and S12L, 68,69 the evolution of these datasets has, generally, followed a prescription of increasing the number of entries, principally by including subtler interactions and/or larger systems. In this respect, the databases of Friesner, 70 Head-Gordon, 71 Shaw, 72 and the recent BFDb of Sherrill, 73 constitute a special category because of their exceptional size (reaching thousands of entries) which are now sufficiently large to be compatible with machine-learning applications. Beyond their conceptual differences, each of these benchmark sets aims at improving the capability of electronic structure methods to describe the energetic aspects of non-covalent interactions.</p><p>In this work, we introduce a dramatic improvement of our previous density-learning approach by making the regression machinery of r(r) compatible with density-tting auxiliary basis sets. These specialized basis sets are routinely used in quantum chemistry to approximate two-center one-electron densities.</p><p>Here, the auxiliary basis sets are used directly to represent the electron densities that enter our machine-learning model, with the additional advantage of avoiding the arbitrary basis set optimization procedures on the machine-learning side. This enhanced framework leverages the transferability of our symmetry-adapted regression method and is capable of learning the all-electron density across a vast spectrum of 2291 chemically diverse dimers formed by sidechain-sidechain interactions extracted from the BioFragment Database (BFDb). 73 The performance of the method is demonstrated through the reproduction of r(r) between and within each monomer forming the dimers. The accuracy of the predicted densities is assessed by computing density-based scalar elds and electrostatic potentials, while the errors made with respect to the reference densities are computed by direct integration on threedimensional grids. As a major breakthrough, the model is used to predict the charge density of a set of 8 polypeptides ($100 atoms) at DFT accuracy in few minutes.</p><!><p>Gaussian process regression (GPR) can be extended to encode all the fundamental symmetries of the O(3) group, effectively allowing machine-learning of all the molecular properties that transform as spherical tensors under rotation and inversion operations. 57,74 In the specic case of the electron density, the scheme relies upon the decomposition of the eld into additive, atom-centered contributions and the subsequent prediction of the corresponding expansion coefficients. 58 In SA-GPR, each molecule is represented as a collection of atom-centered environments, whose relationships and similarities are measured by symmetry adapted kernels. An in-depth discussion about how a symmetry adapted regression model of the electron density can be constructed is reported in the ESI. †</p><p>The decomposition of the electron density in continuous atom-centered basis functions is the cornerstone of the scalability and transferability of our SA-GPR model. Besides being generally desirable, these properties are actually crucial to accurately describe the chemical diversity present in the Bio-Fragment database within a reasonable computational cost. On the other hand, the projection of the density eld onto a basis set leads to an additional error on top of that which can be ascribed to machine learning. In practice, all the efforts placed into achieving a negligible machine-learning error are futile if the overall accuracy of the model is dictated by a large basis set decomposition error.</p><p>Standard quantum chemical basis sets are generally optimized to closely reproduce the behavior of atomic orbitals 75 and results in unacceptable errors if used to decompose the electronic density (Fig. 1). In contrast, specialized basis sets used in the density tting approximation (also known as resolution-ofthe-identity (RI) approximation) [76][77][78][79][80][81][82] are specically optimized to represent a linear expansion of one-electron charge densities obtained from the product of atomic orbitals. Using the RIauxiliary basis sets {f RI k }, the total electron density eld can be expressed as:</p><p>where, D ab is the one-electron reduced density matrix and d ab k are the RI-expansion coefficients. Given a molecular geometry, the value of the basis functions can be readily computed at each point of space, leaving the c k expansion coefficients as the only ingredient needed by the machine-learning model to fully determine r(r) (more details in the ESI †).</p><p>As shown in Fig. 1, the use of the RI-auxiliary basis sets results in nearly two orders of magnitude increase in the overall accuracy with respect to the corresponding standard basis set. The addition of diffuse functions marginally improves the performance of the decomposition, but leads to instabilities of the overlap matrix (high condition number) and increases dramatically the number of basis functions per atom.</p><p>In practice, Weigend's cc-pVQZ/JKFIT 81 basis set (henceforth: cc-pVQZ-RI) offers the best trade-off between accuracy and computational demand and therefore represents the best choice for the density decomposition.</p><!><p>The dataset of molecular dimers has been selected from the side-chain side-chain interaction (SSI) subset of the BioFragment database (BFDb). 73 The original set is made of 3380 dimers formed by amino-acids side-chain fragments taken from 47 different protein structures. Dimers with more than 25 atoms as well as those containing sulfur atoms were not considered. While the total number of sulfur-containing structures is too small to enable the machine-learning model to accurately capture its rich chemistry, the inclusion of the larger systems does not increase dramatically the chemical diversity of the dataset. The nal dataset contains a total of 2291 dimers.</p><p>As shown in Fig. 2, the complete set of 2291 dimers spans a large variety of dominant interaction types, ranging from purely dispersion dominated complexes (in blue) to mixed-inuence (green and yellow) to hydrogen-bonded and charged systems (red). We retain the same classication criteria as in the original database to attribute the nature of the dominant interaction.</p><p>For each dimer, the reference full-electron density has been computed at the uB97X-D/cc-pVQZ level using the resolution of identity approximation for the Coulomb and exchange potential (RI-JK). This implies that RI-auxiliary functions up to l ¼ 5 are included for carbon, nitrogen and oxygen atoms while auxiliary functions up to l ¼ 4 are used for hydrogen atoms.</p><!><p>The training set for the density-learning model was chosen by randomly picking 2000 dimers out of a total of 2291 possibilities. The remaining 291 were used to test the accuracy of the predictions. Given the tremendous number of possible atomic environments ($40 000) associated with such a chemically diverse database, a subset of M reference environments was selected to reduce the dimensionality of the regression problem   (see ESI †). To assess the consequences of this dimensionality reduction, the learning exercise was performed on three different sizes M ¼ {100, 500, 1000} for the reference atomic environments. Fig. 3 summarizes the performance of the machine learning algorithm, expressed in terms of the mean absolute difference between the predicted and the reference densities (QM). Here, only the machine-learning error is shown as the reference densities derive from the RI-expansion of the computed ab initio densities. Since the test set contains molecules of different sizes, the contribution of each dimer has been weighted considering the ratio between its number of electrons and the total number of electrons in the test set.</p><p>where the sum is performed over the 291 dimers of the test set, N e is the total number of electrons, N i e is the number of electrons in a dimer, r i QM (r) and r i ML (r) are, respectively, the ab initio and the predicted density amplitudes at a point. Both integrals of eqn (2) are evaluated in real-space over a cubic grid with step size of 0.1 bohr in all direction and at least 6 Å between any atom and the cube border.</p><p>As shown in the rst panel of Fig. 3, 100 training dimers were sufficient to reach saturation of the density error around 0.5% for M ¼ 100. This result already outperforms the level of accuracy reached in our previous work, which is remarkable given the large chemical diversity of the dataset and the consideration of all-electron densities. Learning curves obtained with M ¼ 500 and M ¼ 1000 show steeper slopes, approaching saturation at about 2000 training dimers with errors that were reduced to $0.2-0.3%. The predicted full-electron densities are ve times more accurate than the previous predictions of valence-only densities (approximately 1%). 58 A more detailed analysis of the M ¼ 1000 learning curve reveals a strong dependence on the nature of the dominant interaction (Fig. 3). Specically, a stronger non-local character in the interaction yields a larger error. This is especially prevalent for dimers dominated by electrostatic interactions (i.e., hydrogen bonds, charged systems), which are characterized by errors that are twice as large as those found in other regimes.</p><p>The origin of this slow convergence arises from two factors. First, only about 20% of the dimers are dominantly bound by electrostatics. 73 The priority of the regression model is thus to minimize the error on the other classes. Second, there is a fundamental dichotomy between the local nature of our symmetry-adapted learning scheme and the long-range nature of the interactions. In fact, the electron density encodes information about the whole chemical system at once, while the machine-learning model represents molecules as a collection of 4 Å wide atom-centered environments. This difference in the spatial reach of the information encoded in the target and in the representation is a limitation. In this respect, a global molecular representation, which includes the whole chemical system, would be more suitable, but this would imply renouncing to the scalability and transferability of the model. Given a large enough training set, however, our SA-GPR model is able to capture the density deformations due to the eld generated by the neighboring molecule. The reason is rooted in the intrinsic locality of density deformations and in the concept of "nearsightedness" 83,84 of all local electronic properties, which constitutes a theoretical justication for a local decomposition of such quantities.</p><p>The fundamental advantage of setting the electron density as the machine-learning target is the broad spectrum of chemical Fig. 4 DORI maps of representative dimers for each type of dominant interaction (DORI isovalue: 0.9). Isosurfaces are color-coded 31 with sgn(l 2 ) r(r) in the range from attractive À0.02 a.u. (red) to repulsive 0.02 a.u. (blue). In particular, sgn(l 2 )r(r) < 0 characterizes covalent bonds or strongly attractive NCIs (e.g. H-bonds); sgn(l 2 )r(r) $ 0 indicates weak attractive interactions (van der Waals); sgn(l 2 )r(r) > 0 repulsive NCIs (e.g. steric clashes).</p><p>properties that are directly derivable from r(r). For instance, the predicted charge densities are the key ingredient in densitydependent scalar elds aimed at visualizing and characterizing interactions between atoms and molecules in real space. Examples of the density overlap region indicator (DORI) 30 are given in Fig. 4 for representative dimers. Compared to the rather featureless r(r), DORI reveals ne details of the electronic structure, which constitute a more sensitive probe for the quality of the machine-learning predictions. In particular, it reveals density overlaps (or clashes) associated with bonding and non-covalent regions on equal footing through the behavior of the local wave-vector (Vr(r)/r(r)). [85][86][87] As shown in Fig. 4, the intra-and intermolecular DORI domains obtained with the SA-GPR densities are indistinguishable from those in the ab initio maps. This performance is especially impressive for the density clashes associated with low-density values, as is typical for the non-covalent domains. All the features are well captured by the predicted densities ranging from large and delocalized basins typical of the van der Waals complexes (in green) to the compact and directional domains typical of electrostatic interactions to intramolecular steric clashes (e.g. phenol, mixed regime). A quantitative measure of the DORI accuracy for the most characteristic basin of each type of interaction is reported in the ESI. † Overall, these results illustrate that the residual 0.2% mean absolute percentage error does not signicantly affect the density amplitude in the valence and intermolecular regions that are accurately described by the SA-GPR model. The highest amplitude errors are concentrated near the nuclei in the region dominated by the core-density uctuations.</p><p>The versatility of the machine-learning prediction is further illustrated by using the predicted densities to compute the molecular electrostatic potential (ESP) for the same representative dimers (Fig. 5). ESP maps based on predicted densities agree quantitatively with the ab initio reference and correctly attribute the sign and magnitude of the electrostatic potential in all regions of space. Importantly, the accuracy of the ESP magnitude remains largely independent of the dominant interaction type. This is especially relevant for charged dimers (electrostatics) as it demonstrates that despite slower convergence of the learning curve for this category, the achieved accuracy of the model is sufficient to describe the key features of the electrostatic potential.  The most widespread applications of ESP maps exploit qualitative information (e.g., identication of the molecular regions most prone to electrophilic/nucleophilic attack) but the electrostatic potentials can be related to quantitative properties such as the degree of acidity of hydrogen bonds and the magnitude of binding energies. [88][89][90][91][92] As a concrete example related to structure-based drug design, we used a recent model that estimates the strength of the stacking interactions between heterocycles and aromatic amino acid side-chains directly from the ESP maps. 88,91,93 This model derives the stacking energies of drug-like heterocycles from the maximum and mean value of their ESP within a surface delimited by molecular van der Waals volume (at 3.25 Å above the molecular plane). 88 Following this procedure, we used the ESP derived from the ML predicted densities to compute the binding energies between a representative heterocycle included in our dataset, the tryptophan sidechain, and the three aromatic amino acid side-chains (Fig. 6).</p><p>Comparison between ab initio and ML predicted stacking interaction energies shows that the deviations in the ESP maps lead to minor errors on the order of 0.05 kcal mol À1 . The largest deviations in the ESP would appear further away from the molecule, beyond the region exploited for the computation of the energy descriptors (i.e., the sum of the atomic van der Waals radii). The predicted ESP shows larger relative deviations far from the nuclei owing to the error propagation of the density predictions r(r) to the electrostatic potential f(r). This can be best understood in the reciprocal space, where the deviations of the potential at a given wave-vector k are related to the density error by d f(k) ¼ 4pdr(k)/k 2 . Because of the k À2 scaling, the error on f(k) increases as k / 0, implying that larger relative errors of the electrostatic potential are expected in regions of space where f(r) is slowly varying (i.e., thus determined by the longwavelength components).</p><!><p>The tremendous advantage of the atom-centered density decomposition is to deliver a machine-learning model that depends only on the different atomic environments and not on the identity of the molecules included in the training set. Thanks to its transferability, the model provides access to density information of large macromolecules, at the sole price of including sufficient diversity, that can capture the chemical complexity of a larger system. The predictive power of this extrapolation procedure is demonstrated by using the machinelearning model exclusively trained on the 2291 BFDb dimers to predict the electron density of 8 polypeptides taken from the Protein DataBank (PDB). 94 The performance of the ML model for each macromolecule, labeled by their PBD ID, is reported in Fig. 8.  Overall, the predictions lead to a low average error of only 1.5% for the 8 polypeptides, which is in line with the highest density errors obtained on the BFDb test set. Relevantly, the largest discrepancies are obtained for 3WNE, which is the only cyclopeptide of the set. The origin of these differences can be understood by performing a more detailed analysis of a representative polypeptide, the leu-enkephalin (4OLR). The errors in this percentage range do not affect the density-based properties, such as the spatial analysis of the non-covalent interactions with scalar elds (Fig. 7 top right panel). Yet, the density differences indicate that the highest absolute errors occur along the amino acid backbone (Fig. 7 lower panels). In addition, the analysis of the relative error with the Walker-Mezey L(a,a 0 ) index 34 shows the highest similarity at the core (99.3%), slowly decreasing while approaching the non-covalent domain (96.3%) (Fig. 7 top le panel). The L(a,a 0 ) index complements the density difference information by showing that the actual density amplitudes and the prediction error do not decrease at the same rate. Nevertheless, the loss of relative accuracy remains modest and the quality of the density is mainly governed by the predictions along the peptide backbone, which are especially sensitive for the more strained 3WNE cyclopeptide. Although similar chemical environments were included in the training set, the error is mainly determined by the lack of an explicit peptide bond motif and cyclopeptides in the training set. While this limitation could be addressed by ad hoc modication of the training set, the overall performance of the machine-learning model is rather exceptional as it provides in only a few minutes, instead of almost a day (about 500 times faster for e.g. enkephalin with the same functional and basis set), electron densities of DFT quality for large and complex molecular systems. For comparison, the superposition of atomic densities (i.e., the promolecular approach), which has been used to qualitatively analyze non-covalent interactions in peptides and proteins (e.g. ref. 32) lead to much larger mean absolute percentage errors (17 times higher, see Fig. S1 in the ESI †).</p><!><p>Given its central role in electronic structure methods, the total electron density is a very promising target for machine learning, since accurate predictions of r(r) give access to all the information needed to characterize a chemical system. Among the many possible properties that can be computed from the electron density, the patterns arising from non-covalent interactions constitute a particular challenge for machine learning models owing to their long-range nature and subtle physical origin. An effective ML model should be transferable across different systems, efficient in learning from relatively small training sets, and accurate in predicting r(r) both in the quicklyvarying region around the atomic nuclei, in the tail andcrucially for the study of non-covalent interactionsin those regions that are characterized by low densities and low density gradients. In this work, we have presented a model that fullls all of these requirements, based on an atom-centered decomposition of the density with a quadruple-zeta resolution-ofidentity basis set, a symmetry-adapted Gaussian Process regression ML scheme, and training on a diverse database of 2000 sidechain-sidechain dimers extracted from the BioFragment database.</p><p>The model reaches a 0.3% accuracy on a validation set, that is sufficient to investigate density-based ngerprints of NCIs, and to evaluate the electrostatic potential with sufficient accuracy to quantitatively estimate residue-residue interactions. The transferability of the model is demonstrated by predicting, at a cost that is orders of magnitude smaller than by explicit electronic structure calculations, the electron density for a demonstrative set of oligopeptides, with an accuracy sufficient to reliably visualize bonding patterns and non-covalent domains using the DORI scalar eld. Even though the model reaches an impressive accuracy (0.5% mean absolute percentage error) for dimers that are predominantly bound by electrostatic interactions, the comparatively larger error suggests that future work should focus on resolving the dichotomy between the local machine learning framework and the long-range nature of the intermolecular interactions.</p><!><p>There are no conicts to declare.</p>