Microsolvation of Histidine—A Theoretical Study of Intermolecular Interactions Based on AIM and SAPT Approaches

: Histidine is unique among amino acids because of its rich tautomeric properties. It participates in essential enzymatic centers, such as catalytic triads. The main aim of the study is the modeling of the change of molecular properties between the gas phase and solution using microsolvation models. We investigate histidine in its three protonation states, microsolvated with 1:6 water molecules. These clusters are studied computationally, in the gas phase and with water as a solvent (Polarizable Continuum Model, PCM) within the Density Functional Theory (DFT) framework. The structural analysis reveals the presence of intra- and intermolecular hydrogen bonds. The Atoms-in-Molecules (AIM) theory is employed to determine the impact of solvation on the charge ﬂow within the histidine, with emphasis on the similarity of the two imidazole nitrogen atoms—topologically not equivalent, they are revealed as electronically similar due to the heterocyclic ring aromaticity. Finally, the Symmetry-Adapted Perturbation Theory (SAPT) is used to examine the stability of the microsolvation clusters. While electrostatic and exchange terms dominate in magnitude over polarization and dispersion, the sum of electrostatic and exchange term is close to zero. This makes polarization the factor governing the actual interaction energy. The most important ﬁnding of this study is that even with microsolvation, the polarization induced by the presence of implicit solvent is still signiﬁcant. Therefore, we recommend combined approaches, mixing explicit water molecules with implicit models.


Introduction
Intermolecular interactions play an important role in many processes at the molecular level, being involved in, for example, the formation of organometallic structures and supramolecular complexes [1,2]. Other important issues are the protein-ligand interactions, where weak but numerous non-covalent interactions stabilize the conformation, the tautomeric equilibria, and the distribution of atomic charges [3]. Another aspect that should be discussed is the solute-solvent mutual relation, which can be decisive in many processes in nature [4,5]. In general, the intermolecular forces are divided into several classes, depending on the strength of the interaction: ion-ion (electrostatic), hydrogen bonding, dipole-dipole (dipole-induced dipole), van der Waals and dispersion forces [6]. The theoretical description of such interactions has been developed, but there are still open questions and problems, which should be studied [7,8]. Moreover, new intermolecular interactions have been discovered and introduced, e.g., halogen bonds [9,10] and charge inverted hydrogen bonds [11,12]. Here, we discuss the intermolecular forces within the solute-solvent interplay.
The theoretical solvation models are divided into two groups: implicit and explicit [13]. The implicit models are usually used to analyze the structural changes and estimate the free energy profiles of studies are concerned mostly with structures, tautomeric transitions, vibrational signatures and-to some extent-the AIM analysis (in [36]). Our intent is different: taking a small subset of the rich conformational space of microsolvated clusters, we would like to observe the effect and nature of the changes induced in the electronic structure of histidine by the presence of a water micro-environment.

Materials and Methods
The model of hisitidine HIP, HIE and HID forms (protonated respectively at both imidazole sites for HIP, at Nδ for HID, at Nε for HIE, see Figure 1) was constructed on the basis of its derivative deposited in the Cambridge Crystallographic Data Centre (CCDC) [38], denoted as ACHIST21 [39] with the deposition number 1840707. Next, the Density Functional Theory (DFT) [31,32] was applied for the quantum-mechanical simulations. The energy minimizations were performed, on the basis of the DFT functional ωB97XD [40] with the triple-zeta split valence basis set 6-311+G(d,p), by Pople [41]. The energy minimizations were performed in the gas phase and with the solvent reaction field for the neutral (protonated at either Nδ or Nε) and protonated (cationic) forms of histidine. The standard structural convergence criteria (maximum force component less than 0.00045 a.u., root-mean-square force less than 0.0003 a.u., maximum step component less than 0.0018 a.u., and root-mean-square step less than 0.0012 a.u.) were employed. In order to reproduce the solvent's influence on the molecular structure of histidine, the Polarizable Continuum Model (PCM) [14], in its IEF-PCM formulation, and water as a solvent were used. Beside the geometry optimization, the harmonic frequencies were calculated for both phases to confirm that the obtained structures correspond to the minimum on the Potential Energy Surface (PES).  The second part of the study is devoted to the investigation of the intermolecular interactions. The three series of histidine-water (one to six water molecules) complexes were constructed for this purpose, assuming the cationic (protonated) form of histidine or two neutral forms of histidine protonated at either Nδ or Nε. The models were prepared with the Molden program's assistance [42] by addition of water molecules to the histidine in the places indicated by the presence of hydrogen bond donors/acceptors, in a manner similar to the procedure described in [36]. The energy minimizations of the obtained models were performed according to the procedure described above. The simulations were carried out in two phases as well. The Atoms-in-Molecules (AIM) theory [33] was applied to study the topology of the complexes, as well as to confirm the presence of the intra-and intermolecular short contacts and hydrogen bonds. The wavefunctions for the Atoms-in-Molecules (AIM) investigations were computed using ωB97XD [40] functional and 6-311+G(d,p) basis set [41].
The Supplementary Material contains also the validation of the computational model; the results of optimizations with the 6-311+G(d,p) basis set and a rich set of DFT functionals; the B3LYP [43,44] and PBE [45] for the gas phase, followed by PCM studies with dispersion-corrected B3LYP-D3 [46] and the Minnesota M06-2X functional [47]. These data serve only as a test proving that the choice of the DFT functional induces quantitative, not qualitative, changes in the solvation shell.
The interaction energy decomposition in the investigated complexes was performed using the Symmetry-Adapted Perturbation Theory (SAPT) [34] at the SAPT2 level [48]. The interaction energy at the SAPT2 level is defined in the following manner-see also Equations (1) to (8) in [48]: (1) where perturbational notation is used in the superscript, e.g., the E 12 elst,r term is an electrostatic term arising from the first-order intermolecular and second-order intramolecular perturbations. The types of terms are electrostatics, exchange (Pauli repulsion), induction (polarization) and dispersion, denoted in the subscripts as elst, exch (ex), ind and disp respectively. The δE HF term contains mostly higher-order corrections (especially for the interactions of induced multipoles), and together with the preceding four contributions, yields the Hartree-Fock interaction energy. The subscript r indicates that this particular term is calculated taking into account orbital relaxation (response) effects-see Equation (118) in [34].
The SAPT2 calculations were preceded by reoptimizations of the histidine-water complexes at the MP2/aug-cc-pVTZ level [49,50], using the density fitting acceleration of the MP2 scheme (DF-MP2) [51]. The optimizations were concluded with the interaction energy estimation at the DF-MP2/aug-cc-pVTZ level using the supermolecular approach, and the Boys-Bernardi counterpoise correction [52] for the basis set superposition error, taking the histidine as one "monomer" and a set of water molecules as the second "monomer" of the "dimer". Unfortunately, the computational demands of the SAPT2 scheme did not allow us to compute the interaction energies for the complexes with 5 or 6 water molecules. The structures optimized at the DF-MP2/aug-cc-pVTZ level served then to calculate supermolecular DF-MP2 and SAPT2 interaction energies, using a smaller aug-cc-pVDZ basis set. The calculations with the smaller basis set were performed for all complexes, so that we could estimate the basis set effect on the SAPT2 results.
The quantum-mechanical simulations were performed using the Gaussian 16 (G16) suite of programs [53], while the AIM calculations were done on the basis of AimAll [54] program. The Psi4. 1.2.1 [55] program was applied to carry out the SAPT calculations and preceding DF-MP2 optimizations. The visualizations of the structures of histidine and its complexes were prepared using the VMD 1.9.3. [56] and Gimp [57] programs. Remark: the Nδ and Nε imidazole nitrogen atoms of histidine are further denoted as ND and NE for simplicity, in a manner similar to the standard PDB atom naming scheme. The hydrogen atoms connected to them are denoted as HD and HE in the Figures 1 and 2 presented in the study. Only atoms of interest are indicated in the Figures 1 and 2, and their numbering scheme was prepared only for the study. Depending on the protonation and tautomeric state, the histidine moiety is labeled as HIP (cationic protonated state), HID (Nδ-H protonation) and HIE (Nε-H protonation), using the convention of the Amber force fields.

Metric Parameters Analysis
The structures of the investigated histidine forms (HIP, HIE and HID) are presented in Figure 1. The DFT simulations for a single molecule of the analyzed three forms were performed in the gas phase, and using the PCM solvation model. The ωB97XD/6-311+G(d,p) level of theory used in this study is very similar to the level used by Lee et al. [37], who chose the 6-311++G(d,p) basis set. The triple-zeta valence basis set, with diffuse and polarization functions, is, on the one hand, already sufficient for the DFT (a more flexible basis set would be chosen for post-Hartree-Fock approaches), and on the other hand, small enough to allow massive calculations in the future studies, for example when analyzing hundreds of snapshots from molecular dynamics trajectories. Further, the ωB97XD functional was chosen due to its performance in the description of non-covalent interactions, especially in the field of the microsolvation of amino acids [37].
The summary of the short contacts found between the imidazole ring and the backbone nitrogen (N1) atom forming a quasi-ring is presented in Table 1, while Table S1I presents an additional dataset obtained with diverse DFT functionals.
In the case of histidine HIP and HID forms, an intramolecular hydrogen bond is present (which was also confirmed by the AIM theory; see text below for more details). However, because of the steric effects, the hydrogen bond is not a linear one. The HIE form does not contain an intramolecular hydrogen bond due to a lack of the corresponding hydrogen atom (HD atom) in the imidazole ring. There is a possible interaction between the nitrogen (ND) atom from the imidazole ring with the NH 2 group. The ND . . . N1 interatomic distance is elongated upon the presence of the polar environment in the case of the histidine HIP form, while for the neutral HID and HIE tautomers, it is shortened upon the transition from the gas phase to PCM. There is no doubt that the presence of the intramolecular short contacts stabilizes the conformation of the histidine residue.
The second part of the analysis is devoted to the histidine (HIP, HIE and HID forms) surrounded by water molecules. The graphical representation of the microsolvation models is presented in the following figures: Figure 2 of the text body, and Figures S2I-S5I of the Supplementary Material. The metric parameters of histidine interacting with water molecules are listed in Table 2 and Tables S2I-S4I. The interactions between water molecules were not taken into account in these tables. Table 1. Metric parameters (interatomic distances are given in [Å]) of intramolecular short contacts obtained as a result of DFT gas phase and with solvent reaction field (PCM with water as a solvent) simulations with ωB97XD functional and 6-311+G(d,p) basis set. For geometric details see Figure 1.   Let us start the discussion of the microsolvation from the histidine HIP form. The selected metric parameters obtained as a result of the DFT simulations are presented in Table 2 and Table S2I. As is visible, the number of water molecules, and the inclusion of the polar environment due to the application of the continuum solvation model, influences the values of the metric parameters. The intramolecular hydrogen bond is not present when we introduced five and six water molecules into the vicinity of the ND-HD bond and NH 2 group. This critical number of water molecules corresponds well with the observation in [37] that six water molecules are necessary to stabilize the zwitterionic form. The study of Rai et al. [36] only includes up to two water molecules, but already significant polarization effects are visible in the histidine residue. The introduction of water molecules to the gas phase models shortened the ND . . . N1 interatomic distance. The intramolecular HD . . . N1 hydrogen bond was shortened when the number of water molecules was added to the model (one to four water molecules). Concerning the results of the simulations with the application of the PCM model, the intramolecular ND-HD . . . N1 interaction was not observed with five or six water molecules in the models. The ND . . . N1 interatomic distance and the HD . . . N1 hydrogen bond were both elongated.

Metric Parameters
The data in Table 2 reveal also some peculiarities in the network of hydrogen bonds connecting the protonated HIP residue with water molecules. A monovalent ion dissolved in water acts as a factor governing the orientation and location of the surrounding water shell. For cations, we expect that the water oxygen atoms are directed towards the cation, while the hydrogen atoms should face the dissolved anionic moiety. A simple (spherical) cation, such as Na + , would form contacts of similar lengths with diverse water molecules. The irregularly-shaped imidazole ring does not have to follow this line of reasoning, and indeed, the ND . . . O and NE . . . O distances differ quite significantly in the gas phase HIP-5 and HIP-6 complexes. This behavior changes sharply when the PCM model is used in addition to microsolvation-the two above-mentioned distances become quite similar. This symmetrization of the microsolvation shell is connected with the symmetrization of the imidazole ring's electronic properties, but it also signifies that even six water molecules is not enough to saturate the polarization of the HIP residue. Table S3I summarizes selected metric parameters of the histidine HIE tautomer. As we could expect, there is a difference in the interatomic distances and angle values of the atoms involved in the intermolecular hydrogen bonds' formation. The PCM's implicit solvation influence is noticeable when comparing the gas phase and solvent reaction field simulations. For example, as is shown in the Table S3I, the O-H . . . O2 intermolecular hydrogen bond was noted in the gas phase simulations, but it was not present for the optimized model with the assistance of the PCM. A similar conclusion could be drawn for the histidine HID tautomer (Table S4I). Comparing the values of the metric parameters obtained after two-phase simulations with various numbers of water molecules, one can observe changes in the interatomic distances and the valence angle value. Moreover, as is shown in Tables S2I-S4I, some interactions (intermolecular hydrogen bonds) are present in the gas phase models, but they were not recognized after the geometry optimization with the PCM solvation model. This signifies that even six water molecules do not saturate the polarization of the histidine residue, and either more explicit water molecules or some form of implicit model (e.g., PCM) is required. Thus, our computational results and metric parameters analysis provided a deeper insight into the histidine residue-water molecules interactions.

Atoms-In-Molecules (AIM) Analysis
The analysis of the electron density via the framework of the AIM theory was performed in this study for several reasons. First, the topology of the electron density field reveals the presence of intraand intermolecular interactions, especially hydrogen bonds. Second, the electron density flow can be evaluated by comparing the AIM net atomic charges in the gas phase and solvent reaction field (PCM model) simulations. Third, this also allows an estimation of the histidine charge transfer from or towards the solvent. These three issues will be discussed below.
The topology of the electron density field is shown schematically in the Figures S8I-S10I in the Supplementary Material. The depicted elements of the molecular graphs are bond paths (solid or dashed), bond critical points (BCPs-green spheres), and less frequent ring critical points (RCPs) or cage critical points (CCPs). The bond paths connect the interacting atoms, and it is important to note that the use of the solid vs. the dashed line is based on the simple criterion of the electron density at the BCP. Thus, one can appreciate that some of the hydrogen bonds paths are marked as a solid bond path, which suggests that these bonds are stronger and important, and exemplified in molecular conformations. This especially concerns the intramolecular ND...N1 contacts, but also some of the carboxylic O1-H1...O (water) interactions. On the other hand, some of the bond paths are significantly non-linear (bent), and they connect rather surprising atom pairs (e.g., two oxygen atoms) (see especially the interactions of the carboxylic O2 atom in the first two panels of Figure S8I). Moreover, in these cases the corresponding BCPs and RCPs are located in close proximity. Even a small structural change could potentially lead to their coalescence (annihilation) and the topological change of the molecular graph. A canonical study of Koch and Popelier on the issue of the AIM-based detection of hydrogen bonds [58] includes correct topology of the electron density as the first criterion to be met by a hydrogen-bonded atom pair. This is an indication that in our case these interactions are, from the point of view of the AIM theory, rather weak, and not of the hydrogen bond type. Their presence stems rather from the purely mathematical necessity of fulfilling the Poincaré-Hopf formula, which is also the case for sterically crowded H . . . H interactions [59]. Now, we will cover the results of the numerical integration of electronic density within atomic basins, which yields AIM net atomic charges. First, let us note that in all the covered cases (18 structures in both gas phase and PCM with water as a solvent), the sum of all the atomic charges deviates from the ideal value (either +1 or 0) by no more than 0.003 e. Such accuracy of the total integration assures us that the charge flow of magnitude larger than 0.01 e is significant, and not due to numerical integration errors. With this in mind, we will turn to the results gathered in Tables 3-5. The last column of these tables, which is the sum of the AIM net charges for the histidine residue, is an indicator of the charge transfer between the histidine and the microsolvation shell of water molecules. In the case of the cationic histidine HIP (Table 3), the cation acts as an electrophile, and the electron density is transferred from the water molecules via hydrogen bonds to the amino acid. The magnitude of this charge transfer is between 0.04 and 0.11 e, above the assumed "significance threshold" of 0.01 e, and it changes with no obvious correlation with the microsolvated cluster size or the number of hydrogen bonds. Table 3. AIM-derived net atomic charges (in the elementary charge units, e) for the selected interaction centers of microsolvated histidine in its protonated form (HIP). Values in bold typeface indicate atoms engaged in a hydrogen bond. Last column: the sum of AIM net charges for the histidine residue. The erratic behavior observed for the cationic HIP is also present for both tautomeric forms of neutral histidine (Tables 4 and 5), but in these cases the amino acid acts most frequently as a weak donor of the electron density (nucleophile), and transfers between 0.02 and 0.09 e to the water shell. Only for the cases where one or two water molecules are in contact with the HID tautomer does the histidine act as an electrophile. This diversity and the ease of modification of the electron-donating or -accepting properties is in good agreement with the role of histidine as a versatile amino acid, capable of acting as, for example, a proton relay in catalytic triads of many enzymes.

HIP-(H 2 O) n Simulation
A detailed understanding of the interplay between the solvent and the histidine requires an analysis of the charge flow around the selected interaction centers of the histidine. We have chosen to investigate more deeply three interaction centers: the vicinity of the imidazole nitrogen atoms, and the carboxylic oxygen atom O2 (retaining its C=O double bond, thus more similar to the situation in a protein backbone). By this method, we also investigate a possible similarity between the properties of the two imidazole nitrogen atoms of histidine. These atoms, Nδ and Nε (denoted ND and NE, respectively), are not symmetrically equivalent from the point of view of the molecular topology, but the presence of a strong electron delocalization in the heterocyclic aromatic ring allows us to expect the electronic behavior of the ND and NE atoms to be very similar. This similarity should be best seen in the protonated histidine (HIP , Table 3), where ND and NE are in the same protonated state. However, the microsolvation with one, two or three water molecules results in the formation of an intramolecular hydrogen bond between the ND-HD moiety and the backbone amine nitrogen atom, while the NE-HE atoms are not hydrogen-bonded. We can see that the presence of the hydrogen bonding results in a significant polarization; for example, in the gas phase HIP-1, the net charges for ND and NE are −1.246 e and −1.206 e, respectively, while the net charges for HD and HE are 0.537 e and 0.476 e. An inclusion of the PCM solvent reaction field reduces this polarization, and the charges of the corresponding atoms are more similar to each other. In agreement with the notion that explicit solvation saturates the capabilities of the donor/acceptor moiety, we see that when a particular moiety (ND-HD or NE-HE) is engaged in a hydrogen bond, it is less affected by the PCM implicit solvent. The most external atoms, in this case HD and HE, are most affected by the PCM solvation, but-roughly summarizing-the effect of PCM implicit water is a 0.02 e net charge decrease of the unsaturated (not-hydrogen-bonded) hydrogen atom, while this decrease is only ca. 0.01 e when such an atom is between explicit donor and acceptor atoms. This is also visible for the carboxyl O2 atom, which, when not engaged in explicit hydrogen bonding, is made more negative (by ca. 0.02 e) under the PCM solvation, but this polarization is not greater than 0.01 e in HIP-5 and HIP-6, where this atom participates in an explicit hydrogen bonding. Let us conclude with the possibility of similarity between the electronic properties of ND-HD and NE-HE; indeed, for the cluster sizes of four or more water molecules, when both moieties are engaged in explicit hydrogen bonding, the atomic charges of the nitrogen atoms are very similar (differing by no more than 0.009 e), indicating an equalization of their electronic properties. The HD and HE hydrogen atoms are not as similar; the net charge differences between them are on average 0.015 e, and can be as large as 0.035 e, but the screening effect of the explicit solvation is still visible. These results agree well with the structural equalization of the N-H . . . O distances between the histidine imidazole ring and the water shell, induced by the use of the PCM model, as discussed in Section 3.1.
The two neutral tautomers, protonated at either the NE (HIE tautomer, Table 4) or ND (HID tautomer, Table 5) nitrogen atom, mostly follow the rules of the behavior outlined above for the cationic HIP. In these two cases, however, it is possible to examine the polarizability of the nitrogen atom with its lone pair, a potential hydrogen bond acceptor. Due to a very flexible electron distribution, enhanced by the lone pair, the non-protonated nitrogen atoms (ND for HIE, NE for HID) under the PCM water solvent model become more negatively charged: by ca. 0.035 e for ND in HIE, and by ca. 0.06 e for NE in HID. This almost twofold difference between the properties of ND and NE shows that, despite the aromatic delocalization of the electron density in the imidazole ring, these nitrogen atoms are not fully equivalent in the electrostatic sense. Finally, it is evident that the protonation of one of the imidazole nitrogen atoms makes this atom polarized; for example, for the HIE-1 structure ( Table 4, no hydrogen bonds to the considered atoms), the protonated NE atom has a net charge more negative by 0.08 e than the non-protonated ND atom. The corresponding difference between ND(-HD) and NE in HID-1 (Table 5) is even larger (0.11 e). Summarizing the AIM part of the study, we would like to stress that the PCM model does not provide the same large degree of polarization, which is obtained with the explicit microsolvation and hydrogen bonding. On the other hand, even with six water molecules, the histidine residue is still visibly polarized by the PCM environment. Thus, it seems advisable to use in further studies a combination of explicit microsolvation and the PCM implicit solvent.

SAPT Analysis of Interaction Energies in the Microsolvated His-(H2O) n Complexes
The reliable assessment of the interaction energies, E int , within the histidine-water complexes has been carried out using the Symmetry-Adapted Perturbation Theory treatment. This approach, although it cannot account for the PCM solvation, is one of the important tools for studying the nature of intermolecular forces. The SAPT level employed in this study is SAPT2, which is roughly equivalent to the MP2 calculus, but it allows for the calculation of separate contributions to the interaction energy. First, we have reoptimized the DFT structures using the MP2/aug-cc-pVTZ level of the theory, and then we used the MP2 structures for the SAPT2 calculations. However, we have found that the calculations at the SAPT2/aug-cc-pVTZ level were not possible for the systems with five or six water molecules. Thus, we have retained the MP2/aug-cc-pVTZ geometries, but we have estimated the interaction energies with the MP2 and SAPT2 techniques with smaller aug-cc-pVDZ basis sets. All these results are summarized in Table S3I. In view of the need to resort to a smaller basis set, it is necessary to verify that the results with smaller and larger basis sets are compatible. This is true in the case of the microsolvated histidine-the MP2 interaction energies obtained with the double-zeta and triple-zeta basis sets are correlated with an r 2 correlation coefficient larger than 0.9998, and the larger basis set yields interactions ca. 8% stronger than the smaller basis. The same degree of correlation is found for the SAPT2 results (r 2 > 0.9995), and the triple-zeta basis set yields 10% stronger interactions than does the double-zeta SAPT2 level. Moreover, the double-zeta MP2 and SAPT2 interaction energies are very strongly correlated (r 2 > 0.9999), and the corresponding MP2 and SAPT2 results differ by no more than 2.3% (0.35% on average). These results confirm that the total SAPT2 energy corresponds closely to the MP2 result, and that the smaller basis set yields results strongly compatible with the aug-cc-pVTZ basis. It should thus be safe to assume that the SAPT2/aug-cc-pVDZ results are also of good quality in case of the His-(H 2 O) 5 and His-(H 2 O) 6 clusters, where we do not have the triple-zeta results.
Establishing the consistency of our procedure, we will now focus on the effect of enlarging the microsolvated His-(H 2 O) n cluster. For this purpose, Table 6 contains the values of the interaction energy divided by the number of water molecules in a given system. In all the three protonation states of histidine, the first H 2 O molecule contributes the most to the stability of the complex, and further water molecules bring smaller energetic gains. However, this trend is not monotonic, and deviations are most prominent for the HIP and HIE at n = 5 water molecules. At n = 6 we already observe the saturation of the hydrogen bond donors/acceptors in the histidine residue, and we do not expect the "per water molecule" contributions to grow significantly. However, due to the inability of the SAPT scheme to deal with the PCM formalism, we cannot establish the degree of polarization incurred by the continuum solvent model. Table 6. The Basis Set Superposition Error-corrected MP2 interaction energies and SAPT2 interaction energies (total per complex, and per one water molecule) calculated for the structures of histidine-(H 2 O) n (n = 1-6) complexes optimized at the MP2/aug-cc-pVTZ level of theory. Nomenclature: HIP-n-protonated histidine, HID-n-the Nδ-H tautomer, HIE-n-the Nε-H tautomer, n-number of water molecules. All values in kcal/mol.

MP2/aug-cc-pVDZ
MP2/aug-cc-pVTZ SAPT2/aug-cc-pVDZ SAPT2/aug-cc-pVTZ Let us note that, as one could expect from the usual magnitudes of the ion-dipole and dipole-dipole forces, the interactions of the protonated histidine (HIP-1 to HIP-6) are stronger (by 10% to 60%) than their counterparts with neutral tautomers HID and HIE; for example, with one water molecule we obtain an SAPT2 Eint for HIP-1, HIE-1 and HID-1 equal respectively to −12.57, −7.98 and −10.31 kcal/mol. This charge-induced strengthening is, curiously, not correlated with the system size.
Unfortunately, the previous studies on the microsolvation of histidine [36,37] are mostly concerned with the tautomeric equilibria, not with the nature of the interactions. We cannot directly compare our calculated interaction energy values with the results of these studies. However, [36] contains data on the estimated strengths of particular hydrogen bonds, evaluated from the electron density using the Espinosa formula [60]. The resulting interactions seem to be stronger than those of our study; for example, HID and HIE hydrated with one water molecule should have a stabilization of ca. −8 to −10.3 kcal/mol at the SAPT2 level, but the estimate given in [36] is close to −16 kcal/mol. The same happens with two water molecules: our SAPT2 results predict E int equal to −13 to −15 kcal/mol, while the data in [36] provide values ranging from −13 to −25 kcal/mol. This comparison cannot be taken too far; the Espinosa formula yields a local parameter, which is an estimate of the energy for a particular bond, with no guarantee that some other region of the interface between molecules is not strongly repulsive.
As noted in the Materials and Methods section, the SAPT approach allows us to identify the physical phenomena behind the interaction energy. A detailed partitioning of E int is given in Table 7 and Figure 3 (for a quick visual examination of the importance of diverse terms and their relative changes), while Figure S11I shows the percentages of the contributions of four principal terms to the E int , conserving their sign (thus the dimensionless E int is always −1). These terms are: ELST (Coulombic interaction of nuclei and frozen electron densities of the monomers), EXCH (Pauli repulsion of the frozen monomer densities), IND (induction/mutual polarization, taking into account a relaxation of monomer electron densities and associated changes to the Pauli repulsion) and DISP (dispersion contributions)-see Equations (2)-(5) in [48] for a detailed list of the contributions included in these four general interaction energy components. The last term, DISP, is the least sensitive to the histidine protonation state. The permanent and frozen electric moments contribute to the ELST term, but in the studied case of the relatively "soft", polarizable residue of histidine, this term is mostly cancelled out by the Pauli repulsion EXCH. In absolute values, ELST and EXCH dominate over the effects of polarization and dispersion. On the other hand, a rather good balance between ELST and EXCH makes polarization the most important factor for the actual interaction magnitude. This is especially true for the protonated HIP residue, wherein IND can be two times larger than the dispersion contribution. On average, for the HIP-n complexes, the magnitude of the dispersion term is ca. 60% of the induction, while for the neutral systems DISP is ca. 80% of the IND term. Surprisingly, the ELST term is not dramatically different between the charged and neutral complexes (although HIP-5 and HIP-6 exhibit very large ELST values), and we conclude that the increased stability of the charged complexes results from the ability of the HIP cation to polarize the surrounding network of the water molecules.

Conclusions
A detailed study of a histidine molecule microsolvated explicitly by up to six water molecules is here presented, and the effect of implicit PCM solvation with water as a solvent is also included. All three tautomers (neutral HID, HIE and protonated HIP) are considered in our study. The structural investigations show that the explicit solvation with a low number of water molecules leads to the formation of intramolecular short contacts within the histidine residue. Some of the intra-and intermolecular short contacts are of the hydrogen bond nature, which is also visible in the topological study within the AIM framework. We have also shown that the two imidazole nitrogen atoms, which are topologically different, are however quite similar from the point of view of the electronic (charge flow) properties. This symmetrization of the imidazole electronic properties, incurred by the PCM model, is followed by the symmetrization of the solvation shell: the imidazole nitrogen-water oxygen intermolecular distances become very similar for the ND…O and NE…O contacts. Symmetrization is also revealed in the convergence of the AIM charges of the two ring nitrogen atoms. The PCM solvation model was not able to yield the same large degree of polarization of histidine as was obtained with the explicit microsolvation and the presence of hydrogen bonding. On the other hand, even with six water molecules, the histidine residue is still additionally polarized

Conclusions
A detailed study of a histidine molecule microsolvated explicitly by up to six water molecules is here presented, and the effect of implicit PCM solvation with water as a solvent is also included. All three tautomers (neutral HID, HIE and protonated HIP) are considered in our study. The structural investigations show that the explicit solvation with a low number of water molecules leads to the formation of intramolecular short contacts within the histidine residue. Some of the intra-and intermolecular short contacts are of the hydrogen bond nature, which is also visible in the topological study within the AIM framework. We have also shown that the two imidazole nitrogen atoms, which are topologically different, are however quite similar from the point of view of the electronic (charge flow) properties. This symmetrization of the imidazole electronic properties, incurred by the PCM model, is followed by the symmetrization of the solvation shell: the imidazole nitrogen-water oxygen intermolecular distances become very similar for the ND . . . O and NE . . . O contacts. Symmetrization is also revealed in the convergence of the AIM charges of the two ring nitrogen atoms. The PCM solvation model was not able to yield the same large degree of polarization of histidine as was obtained with the explicit microsolvation and the presence of hydrogen bonding. On the other hand, even with six water molecules, the histidine residue is still additionally polarized when the PCM solvation is included. This suggests the use of a combination of explicit microsolvation and the PCM implicit solvent in further studies, and we consider it the most important conclusion of the current study.
The nature of the interactions within the microsolvated clusters was investigated using SAPT partitioning of the interaction energy. A high correlation was found between the results obtained with the aug-cc-pVDZ and aug-cc-pVTZ basis sets, which is important due to the quickly growing memory requirements of the SAPT scheme. We show that the aug-cc-pVDZ basis set can be used to model the larger systems adequately. The origin of the microsolvation cluster's stability was also investigated. While the largest absolute contributions to the interaction energy are formed by the electrostatics and Pauli repulsion, these two effects mostly cancel each other out. This makes polarization the most important factor for the actual interaction energy, especially with regards to the protonated HIP residue. The increased stability of the charged complexes of HIP-1 to HIP-6 was put down to the ability of the HIP cation to polarize the explicit solvent molecules included in the cluster.