*Article* **Dissecting Bonding Interactions in Cysteine Dimers**

**Santiago Gómez <sup>1</sup> , Sara Gómez <sup>2</sup> , Jorge David <sup>3</sup> , Doris Guerra <sup>1</sup> , Chiara Cappelli 2,\* and Albeiro Restrepo 1,\***


**Abstract:** Neutral (*n*) and zwitterionic (*z*) forms of cysteine monomers are combined in this work to extensively explore the potential energy surfaces for the formation of cysteine dimers in aqueous environments represented by a continuum. A simulated annealing search followed by optimization and characterization of the candidate structures afforded a total of 746 structurally different dimers held together via 80 different types of intermolecular contacts in 2894 individual non-covalent interactions as concluded from Natural Bond Orbitals (NBO), Quantum Theory of Atoms in Molecules (QTAIM) and Non-Covalent Interactions (NCI) analyses. This large pool of interaction possibilities includes the traditional primary hydrogen bonds and salt bridges which actually dictate the structures of the dimers, as well as the less common secondary hydrogen bonds, exotic X· · · Y (X = C, N, O, S) contacts, and H· · · H dihydrogen bonds. These interactions are not homogeneous but have rather complex distributions of strengths, interfragment distances and overall stabilities. Judging by their Gibbs bonding energies, most of the structures located here are suitable for experimental detection at room conditions.

**Keywords:** cysteine dimers; NBO; NCI; QTAIM; stochastic optimization; hydrogen bonding; salt bridges; non-covalent interactions

#### **1. Introduction**

Cysteine, HOOCCH(NH<sup>2</sup> )CH<sup>2</sup> SH is the only amino acid among the unique list of 20 found in proteins that possesses a thiol functional group [1]. This thiol group, which is comparatively a weaker Brønsted–Lowry acid than O–H in carboxylic acids, is extremely important in biochemistry. Among a large number of known functionalities, the S–H group is responsible for nucleophilic additions to *α*, *β*−unsaturated carbonyl compounds via Michael reactions [2–5], serves as a deprotonation agent [6], and its strong nucleophilicity renders cysteine a key component of the active sites of several protease enzymes [7,8]. In addition to S–H, depending on the conditions, sulfur atoms in cysteine engage in S–S disulfide bonds, which are a central element determining secondary and tertiary structure in proteins [9–11] and are relevant in physiological redox activity. According to the SwissProt databank [12], six percent of all proteins contain at least one disulfide bridge, and the median number of disulfide bonds is two.

Protein· · · protein interaction is one of the central problems in molecular biology. Unfortunately, with present days computational methods, a thorough understanding from a molecular perspective is unattainable because the number of explicit contacts grows exponentially with the size of the protein. For example, insulin is one of the smallest biologically active proteins containing a primary sequence of just 51 amino acids (six cysteines among them), for the dimer of this protein, not counting salt bridges and other intermolecular interactions, there are at least 6 <sup>×</sup> <sup>10</sup><sup>4</sup> possibilities for hydrogen bonding in the classical X *<sup>δ</sup>*<sup>−</sup> · · · <sup>+</sup>*δ*H–Y*δ*<sup>−</sup> description [13]. Typical proteins and other biomolecules contain in excess

**Citation:** Gómez, S.; Gómez, S.; David, J.; Guerra, D.; Capelli, C.; Restrepo, A. Dissecting Bonding Interactions in the Cysteine Dimers. *Molecules* **2022**, *27*, 8665. https:// doi.org/10.3390/molecules27248665

Academic Editors: Qingzhong Li, Steve Scheiner and Zhiwu Yu

Received: 11 November 2022 Accepted: 2 December 2022 Published: 7 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of 1000 amino acids and it is not uncommon to find very large proteins such as titin, which depending on the splice isoform, contains between 27,000 to 35,000 amino acids [14]; thus, the number of specific amino acid· · · amino acid contacts quickly becomes intractable.

In an effort to understand the intricacies of protein· · · protein interactions, the astonishingly large number of specific contacts calls for the use of reduced molecular models, often as gas phase isolated dimers of individual amino acids [13]. In this context, we attempt a detailed study of the cysteine dimers. This is a very complicated issue in its own right: First, there are the two enantiomeric forms of the amino acid. Second, there is the possibility of neutral and zwitterionic forms. Third, there is a large number of conformers of the monomer in a small energy window amenable to form dimers, which, for just the neutral form, has been estimated via ab initio computations to be 42, 51 or 71 depending on the sophistication of the model chemistry [1,9,15]. Notice that six well defined conformers have been experimentally identified via IR and MW spectroscopies [16,17]. Fourth, cysteine contains seven hydrogen atoms and four electronegative atoms; thus, ignoring salt bridges, dispersive dihydrogen interactions and other exotic contacts [13], from the classical X *<sup>δ</sup>*<sup>−</sup> · · · <sup>+</sup>*δ*H–Y*δ*<sup>−</sup> perspective, a total of 28 individual primary plus secondary hydrogen bonds are possible for each dimer. The number of possibilities is reduced to 20, distributed as 12 primary and 8 secondary HBs if the two H–C*<sup>β</sup>* bonds are grouped into just one type and if the two N–H bonds are considered as another type. Fifth, as seen for example in alanine [13], several dimers are attached by more than one contact. Sixth, S–H leads to considerably weaker interactions than O–H, then, the potential energy surface (PES) for the dimers is expected to be considerably richer in weakly bound pairs and thus high levels of electron correlation are needed to correctly describe the intermolecular interactions. Seventh, the environment sensibly impacts the ability of biomolecules to interact in biological settings; thus, using gas phase dimers as a reduced model for protein· · · protein interactions does not seem enough and at least solvent effects must be included.

Cysteine has been thoroughly studied through experiments and computations. Besides the above mentioned publications dealing with the conformations of the monomer [1,9,15–17], Kaminski et al. [18] and Sadlej et al. [19] undertook somewhat exhaustive explorations of the conformational PES for the monomers to rationalize Raman Optical Activity (ROA) and Vibrational Circular Dichroism (VCD) spectra. For the dimers, early studies focused on gas phase and implicitly solvated models with limited explorations of the PES using a few hand constructed configurations [20], later studies considered both explicit water molecules and the neutral and zwitterionic forms [21,22]. There are reports dealing with the formation of the dimers, their stability and bonding (via density differences) when adsorbed in gold surfaces [23,24]. Group IA cations bonded to cysteine dimers have also been studied [25].

In view of the expected complexity arising from the multiple classical donor and acceptor hydrogen bonding sites of cysteine, which as a reduced molecular model has profound implications in the protein interaction problem, the brief summary of the scientific literature dealing with cysteine dimers just exposed reveals an unsatisfactory level of understanding not only of the potential energy surface but of the nature of the intermolecular bonding interactions for cysteine· · · cysteine. The present work is an attempt to remedy this situation. To that end, we undertake systematic explorations of the neutral (*n*) and zwitterionic (*z*) pairs in *n* · · · *n*, *n* · · · *z* and *z* · · · *z* combinations of low lying energy monomers via stochastic samplings of the corresponding PES, and dissected the nature of the interactions using formal quantum descriptors of bonding as provided by the Quantum Theory of Atoms in Molecules [26–29] (QTAIM), the Natural Bond Orbitals [30–33] (NBO) and the Non Covalent Interactions [34,35] (NCI), as discussed in the Section 2.

#### **2. Computer Methods**

Sampling the potential energy surfaces for all possible neutral (72 computed, 6 experimentally found) and zwitterionic forms (12 computed) is not only impossible but unnecessary under the premise that a few representative pairs would capture the vast majority of the specific contacts and thus would provide a sound picture applicable to all

cases. Accordingly, we took two of the experimentally detected neutral monomers [17] (*n*1, *n*<sup>2</sup> in Table 1) and two of the computed lowest energy zwitterions [9] (*z*1, *z*<sup>2</sup> in Table 1) and exhausted all *x*(*x* − 1)/2 + *x* = 10 possible dimeric combinations for *x* = 4. Each pair was superimposed at the center of a cubic box of 512 Å<sup>3</sup> (8 Å side) and was allowed to evolve under simulated annealing conditions [36–38] as implemented in the ASCEC program [39]. Superimposing the interacting system at the center of the box gives the algorithm the worst possible starting point (we call this the big bang initial conditions) and guarantees that the located stationary points within the corresponding PES are free of any structural bias. ASCEC [40,41], after its Spanish acronym Annealing Simulado Con Energía Cuántica, randomly explores the quantum energy landscape for the dimeric interaction, subjects the generated structures to a modified Metropolis acceptance test, and delivers a set of candidate structures that undergo further optimization via gradient following techniques and characterization as true minima via harmonic vibrational analysis. Each one of the 10 possible dimeric combinations was treated to duplicate ASCEC runs. All ASCEC runs and geometry optimizations were carried out in an aqueous environment represented by a continuum according to the PCM (ASCEC) and IEFPCM (optimization) models. [42–44].

**Table 1.** Structures of the B3LYP–D3/6–311++G(*d*, *p*) monomers of neutral (*n*<sup>1</sup> , *n*2) and zwitterionic (*z*<sup>1</sup> , *z*2) L–cysteine. In each case, ∆∆*G* are the corresponding differences in Gibbs energies at room conditions with respect to *n*<sup>1</sup> , *z*<sup>1</sup> , the lowest energy monomers. Descriptors of intramolecular bonding derived from QTAIM and NBO are included along with the specific NBO orbitals responsible for the interactions. ∆∆*G*, *E* (2) *d*→*a* in kcal mol−<sup>1</sup> , all other descriptors in a.u. *n*<sup>1</sup> and *n*<sup>2</sup> have been experimentally detected [17].


Final equilibrium geometries and Gibbs energies for every located dimer computed with the Gaussian09 suite of programs [45] are reported here using the dispersion corrected B3LYP–D3/6–311++G(*d*, *p*) model chemistry. Binding using the Gibbs energies at room conditions (1 atm, 298.16 K) are calculated as the negative difference between the energy of the cluster and the energy of the fragments, *BE* = −(*Ecluster* − ∑ *Ef ragments*), in this way, positive binding energies indicate strongly bonded clusters.

Once the molecular wavefunctions and electron densities for the optimized geometries are recovered by the procedure just stated, we use them to gain insight into the nature of intermolecular bonding interactions using the tools provided by QTAIM, NBO, and NCI following strategies described elsewhere [46–53]. At this point, we state that we use those methods as well established analysis tools, the interested reader is directed to the specialized literature for detailed discussions of their merits and shortcomings and for a description of how the calculated descriptors are related to bonding [6,33,49]. In short, use the Multiwfn suite [54] to find the bond critical points (BCPs, **r***c*) corresponding to

intermolecular interactions, analyze their properties, i.e., the electron density *ρ*(**r***c*), its Laplacian <sup>∇</sup>2*ρ*(**r***c*), the total, kinetic, and potential energy densities <sup>H</sup>(**r***c*) <sup>=</sup> <sup>G</sup>(**r***c*) <sup>+</sup> <sup>V</sup>(**r***c*), and the virial ratio |V(**r***c*)|/G(**r***c*). With the same program, we calculate the Laplacian of the electron density, a scalar field that gives direct information about the most probable regions to find the electrons. Then, we use NBO6 [55] to pinpoint the specific orbitals involved in the intermolecular interactions associated with each BCP and estimate the strength of the interaction via second order perturbation energy for the interaction between the donor and acceptor orbitals, *E* (2) *d*→*a* . The NCIPLOT program [56] was used to derive the non-covalent interaction surfaces. Jmol and VMD [57,58] were used to visualize the molecular structures, and their related surfaces and orbital interactions.

#### **3. Results**

Table 1 shows the structures of the *n*1, *n*<sup>2</sup> neutral and *z*1, *z*<sup>2</sup> zwitterionic forms chosen in this work to study the dimers of cysteine. The four non-covalent intramolecular interactions, as derived from QTAIM are displayed as dotted lines along with the involved NBOs. QTAIM and NBO descriptors are included as well. Only *n*<sup>2</sup> has a structure free from intramolecular hydrogen bonds while *z*<sup>1</sup> exhibits two intramolecular contacts. Except for *n*1, all intramolecular interactions are characterized as weak, long range contacts because of the positive Laplacians, relatively small accumulation of electron densities at BCPs, virial ratios smaller than 1, and positive bond degree parameters. However, the *n*<sup>N</sup> → *σ* ∗ O−H interaction in *n*<sup>1</sup> is uncharacteristically strong, with values for the bonding descriptors that in every case surpass those of the archetypal hydrogen bond in the water dimer [30,59]. These intramolecular contacts are quite important because the formation of the dimers will usually involve investing energy to eliminate those interactions in favor of dimeric contacts. The electrostatic potentials in Figure S1 of the Supplementary Information show the blue and red regions which are more susceptible to the formation of intermolecular contacts according to the classic electrostatic view of hydrogen bonding.

#### *3.1. Structures and Energies*

Complex and rich potential energy surfaces are uncovered by our stochastic searches in every case. A total of 746 distinct well defined dimers were located in the 10 monomer + monomer possible combinations. Table 2 lists the number of structural isomers for each PES and also shows that the vicinities of the putative global minima are populated with other close energy dimers; thus, all structures accounting for populations larger than 1% are within 2.1 kcal/mol of the *n* · · · *n* lowest energy structure, and so on. This point is emphasized by the results shown in Tables S1–S3 in the Supplementary information and in Figure 1, which clearly show that there are no dominant isomers.

**Table 2.** Summary of structural and energetical properties of the cysteine dimers. ∆*G* range: Gibbs energy difference between the highest and lowest energy structure in kcal/mol. %*x<sup>i</sup>* : isomer populations listed in Tables S1–S3 in the Supplementary Material.


Although 746 is a very large number of structures and is considerably higher than the numbers reported in any of the previous studies, we recognize that given the complexity of our problem, no stochastic or analytic search algorithm is able to locate all possible geometries. A representative set including only those dimers with populations exceeding 5% within each PES is shown in Figure 1, along with the NBOs responsible for the inter-

molecular interactions. Cartesian coordinates for all 746 structures located in this work are provided in the Supplementary Material.

**Figure 1.** Lowest energy structures and the NBOs responsible for the strongest intermolecular interactions in the neutral *n* · · · *n* (top), neutral + zwitterionic *n* · · · *z* (middle) and zwitterionic *z* · · · *z* (bottom) B3LYP–D3/6–311++g(*d*, *p*) potential energy surfaces of the cysteine dimers under the continuum IEFPCM solvent model for water. Solid/meshed surfaces correspond to charge donor/acceptor orbitals, respectively. *BE*: binding energies in kcal/mol calculated using the Gibbs free energies at room conditions. See Table 1 for the structures of *n*<sup>1</sup> , *n*2, *z*<sup>1</sup> , *z*2, the isolated monomers. Only those structures with populations (%*x<sup>i</sup>* ) higher than 5% within each PES are included. Energetics for the entire set of 746 dimers is provided in Tables S1–S3 of the supplementary material.

On the basis of purely ZPE–corrected electronic energies (Tables S1–S3), all cysteine dimers are stable towards fragmentation into the corresponding monomers, however, Figure 2 shows that consideration of temperature and entropy leads to 235 clusters (80 *n* · · · *n*, 154 *n* · · · *z* and 1 *z* · · · *z*) having negative binding energies as calculated from the Gibbs energies; thus, those particular structures correspond to unstable dimers and are not amenable to experimental detection at room conditions in aqueous environments, a fact that is emphasized by their %*x<sup>i</sup>* ≈ 0 populations. Notice the contrast with the 416 *n* · · · *n* gas phase equilibrium structures reported for the Alanine dimers [13], which are all strongly bonded. Binding energies show a clear *BEn*···*<sup>n</sup>* < *BEn*···*<sup>z</sup>* < *BEz*···*<sup>z</sup>* ordering; thus, there is a marked preference for charged cysteine dimers in aqueous environments.

Figure 2, showing distribution plots of the Gibbs binding energies leads to a few relevant observations: Dashed vertical lines indicate the expected values of the binding energies using the Boltzmann populations of the Gibbs energies within each PES as weighing factors. 14.3, 16.6 and 20.9 kcal/mol are obtained for *n* · · · *n*, *n* · · · *z*, *z* · · · *z*, again showing a preference for charged dimers in aqueous environments. To put these binding energies in context, they are larger than the gas phase Gibbs binding energies of acetamide and acetic acid, which are 2.1 and 3.8 kcal/mol respectively, according to Copeland et al. [60] Notice that the same authors reported substantially higher binding energies when using only the ZPE-corrected electronic energies: 12.3 kcal/mol for acetamide and 14.7 kcal/mol for acetic acid. High ZPE-corrected binding energies have also been reported for the lowest

energy structures in similar systems: 16.6 kcal/mol for the dimers of formic acid according to Kalescky et al. [61] and 12.7 according to Farfán et al. [62], 19.0 kcal/mol for the dimers of carbonic acid [63], and 20.9 kcal/mol for the alanine dimers [13]. Tables S1–S3 in the Supplementary material show exactly the same trend for all cysteine dimers calculated here, that is, comparatively much higher binding energies are obtained when only ZPE-corrected are considered with expected values of 25.9, 28.9, 33.7 kcal/mol for the *n* · · · *n*, *n* · · · *z*, and *z* · · · *z* cases, respectively. Notice that these numbers are up to over 6 times larger than the 5.0 kcal/mol binding energy arising because of a single hydrogen bond in the archetypal water dimer [59]. Finally, notice that those structures being unstable towards fragmentation (*BE* < 0) have minimal populations and thus do not contribute to the expected value of the binding energy. The role of dispersive interactions is clearly seen in the fact that when the D3 correction is removed from B3LYP, all strongly bound isomers become unstable towards fragmentation (values within parentheses in Tables S1–S3). For the cysteine dimers with positive binding energies, Tables S1–S3 show that the structures with the largest populations are strongly bonded.

**Figure 2.** Distribution of binding energies of the cysteine dimers using the Gibbs energies. Dashed vertical lines mark the expected value for each potential energy surface using the Boltzmann populations at room conditions as weighing factors.

As general structural features of the cysteine dimers, we point out that in all cases where neutral monomers are involved, *n*<sup>2</sup> (no intramolecular HB, Table 1) leads to lower energy dimers. Additionally, except for *Dnz* 1 , in all structures that contain *n*1, the intramolecular HB in *n*<sup>1</sup> remains in the dimer. A surprising result is that contrary to the well known structures of the dimers of carboxylic acids, out of the 244 well characterized *n* · · · *n* local minima, only two (*Dnn* 5 , %*x<sup>i</sup>* = 1.7 and *Dnn* 7 , %*x<sup>i</sup>* = 1.3) exhibit the traditional eight atom, cyclic double C=O· · · H–O stabilizing network. We attribute this to two factors: one, the influence of the solvent which favors other configurations, and two, the intramolecular hydrogen bond occupying the O–H bond in *n*<sup>1</sup> remains in all but one *n* · · · *n* dimer; thus, this bond is not available for intermolecular bonding (see the dissection of intermolecular bonding interactions below).

#### *3.2. Bonding*

The configurational space for the cysteine dimers is complex and rich. We located and characterized a total of 746 structures and there might as well be many more. This geometrical variety arises because of the large number of possible interactions discussed above. Our stochastic search and subsequent dissection of bonding interactions (see below) uncovers an astonishing total of 80 well characterized physically different types of direct intermolecular contacts listed in Tables 3 and 4. Gratifyingly, the found structures account for every single one of the 20 possible hydrogen bonds among the monomers as exposed in

the Introduction and also reveal additional salt bridges, dihydrogen bonds, and a number of exotic X· · · Y (X, Y = O, S, N, C) and C· · · H–C, C· · · H–S contacts. Notice that the *n* · · · *n* dimers exhibit a well balanced field of all non-charged interactions while the *z* · · · *z* dimers favor the salt bridges by a long shot (159 appearances) and, to a lesser extent, other interactions where only one of the fragments is charged. What should be clear is that the largest contributors to the stabilization of the dimers are N· · · H–O interactions in the *<sup>n</sup>* · · · *<sup>n</sup>* dimers, C=O<sup>−</sup> · · · H–O in *<sup>n</sup>* · · · *<sup>z</sup>*, and C=O<sup>−</sup> · · · H–N<sup>+</sup> salt bridges in *<sup>z</sup>* · · · *<sup>z</sup>*. We think it is important to point out that, as a general rule, due to the comparatively larger interaction strength, it is the primary neutral, charged, or salt bridges forms of HBs that determine the molecular geometry of the dimers while secondary HBs and exotic contacts are a consequence of the structure (*vide infra*), however, the collective action of the multiple weak interactions on the stabilization energy of each cluster should not be ignored.

Our topological analysis of the electron densities of the 746 equilibrium structures located in this work affords a total of 2894 intermolecular contacts, which are collected into 80 different types in Tables 3 and 4. Without a single exception, positive Laplacians at bond critical points (see Figure S2 in the supplementary material) indicate that bonding in the *n* · · · *n*, *n* · · · *z* and *z* · · · *z* cysteine dimers occurs via closed shell interactions, in the form of either ionic bonding or long range weak interactions. We dissect the nature of intermolecular interactions next.

**Table 3.** Properties of the 20 types of primary hydrogen bonds, 10 types of secondary hydrogen bonds, and 12 types of dihydrogen bonds stabilizing the cysteine dimers. *NA*···*<sup>B</sup> i* is the number of times that the interaction appears in the corresponding type of dimers. *φ<sup>d</sup>* , *φ<sup>a</sup>* are the charge donor and acceptor orbitals as identified from NBO. An example of a dimer containing each particular interaction is given in the rightmost column.



**Table 3.** *Cont.*

**Table 4.** Properties of the "exotic" intermolecular contacts found in the cysteine dimers. *NA*···*<sup>B</sup> i* is the number of times that the interaction appears in the corresponding type of dimers. *φ<sup>d</sup>* , *φ<sup>a</sup>* are the charge donor and acceptor orbitals as identified from NBO. An example of a dimer containing each particular interaction is given in the rightmost column.



#### 3.2.1. Interaction Distances

Figure 3 shows the distribution of the distances associated with individual intermolecular contacts separated by interaction type, that is, primary and secondary hydrogen bonds, dihydrogen bonds, and exotic contacts for all dimers. Remarkably, the spectrum of A· · · B distances for direct intermolecular contacts covers a wide range, from the very short (1.50 Å for a C=O<sup>−</sup> · · · H–O in *<sup>D</sup>nz* <sup>38</sup> ) to the very large (4.19 Å for the exotic S· · · S in *<sup>D</sup>nn* <sup>77</sup> ), which sensibly departs in both directions from the reference 1.98 Å in the isolated gas phase water dimer. Notice that regardless of the constituting monomers, only primary hydrogen bonds fall below 1.98 Å. In a classical sense, a zwitterion may be conveniently seen as two remote charge islands within the same molecule, in this view, the effect of the charges in the structural complexity of the cysteine dimers is clear: on one hand, intermolecular distances are reduced for the dimers with more charge islands, i.e., *r nn AB* > *r nz AB* > *r zz AB*, on the other, the structural complexity is also sensibly reduced for the more charge-separated species because the *n* · · · *n* distributions have more peaks than *n* · · · *z* which in turn have more peaks than *z* · · · *z*. In addition, it may be argued that among all the types of interactions stabilizing the cysteine dimers, salt bridges should be the strongest and thus the most important structural determining factor whenever they occur. Indeed, the lowest energy *z* · · · *z* dimers with populations larger than 5% shown in Figure 1, are actually stabilized by two salt bridges. Notice that the center of the peak for the distribution of C=O<sup>−</sup> · · · H-N<sup>+</sup> distances (1.66 Å, Figure 3C) is actually larger than 1.57 Å, the center of the peak for the distribution of C=O<sup>−</sup> · · · H-O interactions (Figure 3B), which are *a priori* not as strong as the salt bridges but which dictate the structures of the *n* · · · *z* dimers, the reason for this apparent contradiction is that formation of the two salt bridges confers structural rigidity to the clusters.

When immersed in a continuum aqueous environment, there is partial dissociation of the O–H bonds upon the formation of the dimers. Figure 4 shows the changes in the corresponding distances and Wiberg Bond Indices (WBI) compared against the reference monomer. Evidence for partial dissociation is provided by the peak centered at ≈0.58 WBI, which actually corresponds to O–H groups of the low energy, high population dimers where neutral monomers are involved.

## 3.2.2. Electron Densities at Bond Critical Points, *ρ*(**r***c*)

The relationship between electron density at bond critical points and the nature of the interaction is clear: large accumulations of electron densities at BCPs indicate that the electrons are shared between two fragments or atoms, otherwise known as covalent bonding. Conversely, small electron densities at BCPs indicate that the electrons are displaced towards the nuclei, thus signaling either ionic bonding or long range interactions. Figure 5 shows the values for the calculated electron densities at the 2894 bond critical points associated to intermolecular interactions in the 746 cysteine dimers. Electron densities at those points cover the [9.1 <sup>×</sup> <sup>10</sup>−<sup>4</sup> , 7.6 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ] a.u. interval. These values are sensibly smaller than the 0.24 and 0.35 a.u obtained for the covalent C–C and O–H bonds in *Dnn* 1 . The smaller electron densities correspond to secondary HBs and exotic contacts while among primary HBs, those with the smallest densities involve the S–H group. Only some primary HBs and salt bridges have larger densities than the reference water dimer. The

distributions plotted in Figure 5 are wide; thus, there are many possibilities for the same type of interaction. Finally, as expected [49,64–66], there seems to be an inverse correlation between interaction distance and electron density at intermolecular BCPs.

**Figure 3.** Distributions of the A· · · B distances for intermolecular contacts for all dimers of cysteine found in this work. The distributions are fitted to the actual histograms, so the center of the peaks of the distributions are statistically relevant. The top row shows only primary hydrogen bonds including salt bridges. The bottom row shows secondary hydrogen bonds, dihydrogen bonds, and all exotic interactions. The left column is reserved for the *n* · · · *n* dimers (subfigures (**A**,**D**)), the middle column for *n* · · · *z* (subfigures (**B**,**E**)) and the right column for *z* · · · *z* (subfigures (**C**,**F**)). All distances taken from the B3LYP–D3/6–311++G(*d*, *p*) potential energy surfaces with water represented as a continuum solvent. The dashed vertical lines at 1.98 Å mark the reference H· · · O distance in the gas phase water dimer.

**Figure 4.** Variation of the O–H bond properties as a consequence of intermolecular interactions. The reference values for the *n*<sup>2</sup> cysteine monomer are included as vertical dashed lines.

**Figure 5.** Electron densities at bond critical points for intermolecular contacts for all dimers of cysteine found in this work. The top row shows only primary hydrogen bonds including salt bridges. The bottom row shows secondary hydrogen bonds, dihydrogen bonds, and all exotic interactions. The left column is reserved for the *n* · · · *n* dimers (subfigures (**A**,**D**)), the middle column for *n* · · · *z* (subfigures (**B**,**E**)) and the right column for *z* · · · *z* (subfigures (**C**,**F**)). All values taken from the B3LYP– D3/6–311++G(*d*, *p*) potential energy surfaces with water represented as a continuum solvent. The dashed vertical lines mark the reference value for the gas phase water dimer.
