Understanding the role of intramolecular ion-pair interactions in conformational stability using an ab initio thermodynamic cycle

Intramolecular ion-pair interactions yield shape and functionality to many molecules. With proper orientation, these interactions overcome steric factors and are responsible for the compact structures of several peptides. In this study, we present a thermodynamic cycle based on isoelectronic and alchemical mutation to estimate intramolecular ion-pair interaction energy. We determine these energies for 26 benchmark molecules with common ion-pair combinations and compare them with results obtained using intramolecular symmetry-adapted perturbation theory. For systems with long linkers, the ion-pair energies evaluated using both approaches deviate by less than 2.5% in vacuum phase. The thermodynamic cycle based on density functional theory facilitates calculations of salt-bridge interactions in model tripeptides with continuum/microsolvation modeling, and four large peptides: 1EJG (crambin), 1BDK (bradykinin), 1L2Y (a mini-protein with a tryptophan cage), and 1SCO (a toxin from the scorpion venom).

Intramolecular ion-pair interactions yield shape and functionality to many molecules. With proper orientation, these interactions overcome steric factors and are responsible for the compact structures of several peptides. In this study, we present a thermodynamic cycle based on isoelectronic and alchemical mutation to estimate intramolecular ion-pair interaction energy. We determine these energies for 26 benchmark molecules with common ion-pair combinations and compare them with results obtained using intramolecular symmetry-adapted perturbation theory. For systems with long linkers, the ion-pair energies evaluated using both approaches deviate by less than 2.5% in vacuum phase. The thermodynamic cycle based on density functional theory facilitates calculations of saltbridge interactions in model tripeptides with continuum/microsolvation modeling, and four large peptides: 1EJG (crambin), 1BDK (bradykinin), 1L2Y (a mini-protein with a tryptophan cage), and 1SCO (a toxin from the scorpion venom).

I. INTRODUCTION
Ion-pair interactions are ubiquitous and are the fifth most prevalent in biological systems as salt-bridges [1]. There are many definitions of salt-bridge based on the constitution and spatial orientation of the interacting moieties [2][3][4][5][6]. In peptides, such interactions are predominantly found in guanidinium-carboxylate or ammoniumcarboxylate ion pairs [7]. Presence of only a few saltbridges can have a significant impact on both the structure and the function of bio-macromolecules [8][9][10][11][12][13][14][15][16][17][18][19]. The stability of a salt-bridge depends on its local environment and exposure to polar solvents. Depending on these factors, they can be stabilising [20][21][22][23][24][25], insignificant [26][27][28][29], or even destabilising [30]. Experimentally, the presence of a salt-bridge can be inferred through nuclear magnetic resonance spectroscopy [31] or infrared vibrational spectroscopy [32,33]. While vibrational spectroscopy also detects different salt-bridging patterns (such as endon/side-on or monodentate/polydentate [32]), they cannot provide any insight on their relative thermodynamic stabilities. The absolute value of the salt-bridge interaction energy of a folded protein is not an experimental observable. Between a folded and an unfolded state of a macromolecule, the Gibbs free energy of folding, ∆G folding , can be experimentally determined in terms of stability constants. This quantity is essentially the sum of contributions due to structural relaxation or strain, ∆G strain , salt-bridge (or ion-pair) interaction, ∆G ip , and reorganization due to interactions with solvent molecules. The contribution due to salt-bridge, ∆G ip , can compensate for the other two terms resulting in a net negative ∆G folding . Experimentally, ∆G ip is determined using denaturation titrations of the wild-type protein and its single and double mutants [31,34]. These experiments have been extremely useful to shed light on the strength * ramakrishnan@tifrh.res.in of different salt-bridges across many different proteins. The two salt-bridges (E6/R92 and E62/K46) in thermophilic ribosomal protein L30e stabilize the system by 0.5-1.2 kcal/mol [35]. A double mutant cycle analysis revealed the K11/E34 surface salt-bridge in ubiquitin to have a strength of 0.24 kcal/mol [36], while the conserved salt-bridges (R69/D93 and R83/D75) in barnase revealed an interaction energy of 3.0-3.5 kcal/mol [37]. White et al. used molecular dynamics simulations to investigate charged amino-acids in water, when interaction strength (∆G int ) was found to vary from 0.5-2.0 kcal/mol [38]. Typically, the stabilization due to salt-bridges in most proteins are within the 4-5 kcal/mol range [7]. However, studies have shown that this value can be influenced by interactions with osmolytes [39]. In general, engineered surface salt-bridges are destabilizing but conserved saltbridges present inside the protein are stabilizing [40].
We note that in denaturation experiments the interacting moieties in ion-pairs are in two very different paradigms: In the folded state, the ion-pair is not exposed to any solvent molecules, while in the unfolded state, it is exposed to both solvent molecules and ions. Furthermore, when ion-pairs are exposed to the solvent media, the interaction strength is significantly quenched due to the dielectric constant of the solvent. Perhaps, the ideal scenario of understanding unfolding must be one where contributions due to solvent can be completely ignored, but at this stage it is experimentally unfeasible. Again, in double mutant experiments, the single mutants involve replacing one of the charged residues with a neutral one. The assumption is that the effects due to mutation will exactly cancels out. However, the wild-type and mutated peptides may exhibit different conformational preferences resulting in non-cancellation of energy terms due to structural relaxation. Another computational protocol would be to access a conformer of the biomolecule where the interacting groups are well-separated [41]. Not only does this exercise become increasingly challenging when it comes to estimating salt-bridges in large peptides, but more fundamentally, the strength of the salt-bridge is now dependent on the type of elongated conformer accessed.
Focusing only on the enthalpic contributions, a theoretical estimation of intramolecular ion-pair interaction energy, E intra ip , is more challenging than that of the intermolecular counterpart, E inter ip . This is because the latter is a well-defined quantity related to the dissociation energy, while such analogy is not suitable for intramolecular bonding. E inter ip can be estimated using the supermolecule approach as the difference between the total energies of the complex and its isolated monomers. E inter ip determined for model salt-bridges provides insights into E intra ip of proteins.
Besides the supermolecule approach, intermolecular interaction energies can also be estimated using energy decomposition analysis (EDA) [42][43][44][45][46] or symmetryadapted perturbation theory (SAPT) [47][48][49][50][51][52][53][54][55][56][57]. Both EDA and SAPT provide a breakdown of interaction energies quantifying various physical effects. These methods have been benchmarked using results from the supermolecule approach based on highly accurate wavefunction methods as references. The accuracy of SAPT is limited by the reference wavefunction used to model the interacting fragments and the order of perturbation corrections. In the intramolecular analog of SAPT (I-SAPT), Hartree-Fock (HF) molecular orbitals (MOs) are localized on the interacting fragments and the linker. Collectively, this procedure is denoted as F/ISAPT0, where 0 indicates the use of HF states for the fragments prior to the application of perturbation theory [58]. The initial applications of F/ISAPT0 to organic molecules such as pentanediol have yielded physically meaningful partitioning of the interaction energy, of E intra ip [53]. In general, the SAPT formalism delivers excellent predictions for weak non-covalent interactions in charge-neutral systems, while they underperform for ionic systems [59][60][61]. Hence, for zwitterions stabilized by intramolecular ion-pair interactions, it is desirable to compare F/ISAPT0's predictions to results from other approaches.
In this study, we present a thermodynamic cycle to estimate E intra ip by introducing isoelectronic atomic mutations. For model dimers with intermolecular ionic bonding, we benchmark SAPT0 with two basis sets against CCSD(T) references. We use the best performing basis set in combination with F/ISAPT0 method, and determine intramolecular ion-pair interactions for a set of 26 molecules. These results are compared with E intra ip predictions using the thermodynamic cycle based on total energies from the domain-based local pair natural orbital coupled cluster with singles doubles and perturbative triples method, DLPNO-CCSD(T) [62,63] and various density funcitonal theory (DFT) approximations. We probe the convergence of E intra ip with different fragmentation procedures in vacuum and implicit solvation regimes. We explore the effect of microsolvation on structural preferences, salt-bridge interaction energies, and conformer stabilities of four model tripeptides. Finally, we estimate salt-bridge interaction energies in four large peptides and rationalize their experimentally noted geometries.

II. THEORY
Here, we focus on the ion-pair interaction energy, E ip , and its contribution to the stability of conformers. An idealized potential energy profile of an ion-pair system connecting a compact conformer and a longer one is presented in FIG.1a. The net ion-pair bonding in the compact conformer, ∆E ip , overcomes the steric strain associated with folding.
Hence, on the potential energy surface (PES) of the non-polar mutant, the longer conformer is energetically favored. In principle, it is possible to apply the schematics shown in FIG.1a to estimate the difference in E ip between a compact and a long conformer (∆E ip ). In practice, the interest lies in estimating E ip for an experimentally observed folded geometry such as those collected in the protein data bank (PDB) [64], without depending on a reference unfolded geometry. While determining an unfolded conformer is possible by exploring the torsional PES (as routinely explored in molecular dynamics studies), our thermodynamic scheme discussed below, avoids this step.
In FIG.1b we present a scheme to estimate E intra ip of an ion-pair system using a thermodynamic cycle. The specific case of ammonium-carboxylate ion-pair is used as an example. Here, wt denotes the system with an intramolecular ion-pair bonding, and mwt i is obtained by combinatorially mutating the polar groups with isoelectronic substitutions. Thus, N atoms in the positively charged polar moiety in wt are substituted with C atoms while in the negatively charged moiety, O atoms are replaced by F. The difference between the total energies of wt and mwt quantifies the ion-pair interaction that is present in wt while absent in mwt, and contributions due to the net change in the external potential and nuclear repulsion going from wt to mwt: The alchemical contribution, E al. , can be estimated separately by using the total energies of the fragmented monomers of wt and mwt, where the Cartesian coordinates of atoms of the polar groups are kept frozen. Hence, Fragmenting wt to infinitely separated units comprising the charged moieties and a limited degree of their local chemical environment yields two fragmented systems, denoted as fr pos , and fr neg (commonly denoted as fr ). These polar moieties are saturated with terminal groups of varying length whose performances are later investigated. Further, for all possible mutation paths in the wild type system, wt−→mwt i , there are corresponding mutations in the fragments, fr pos/neg −→mfr pos,j/neg,k . Our scheme can be applied to determine E ip with microsolvation to capture solvent effects. In this regime, fr mfr Fragmentation Figure 1. Development of the thermodynamic scheme. a) Schematic illustration of typical potential energy profiles of a molecule in the zwitterionic form (bottom) and in its alchemically mutated non-polar form (top). The horizontal axis is denoted by θ, a general reaction coordinate. For simplicity, effect of hydration is not considered. Net stabilization of the compact structure of an ion-pair can be estimated by assuming a cancellation of alchemical changes due to introduction of charge neutral atoms or groups, E alch. ≈ E alch. . b) Thermodynamic cycle for estimating intramolecular ion-pair interaction energy, Eip. As an example, we consider the salt-bridging interaction between the guanidinium and carboxylate moieties (highlighted with a red ellipse).
Here, wt denotes the wild-type variant of the molecule. Atomistic mutations to remove charges on moieties leads to mutated wild-type variants denoted by mwti, where i runs from 1 to 6, spanning all possible mutation paths. We now fragment wt to two separate units, one containing the positively charged moiety and the other containing the negatively charged moiety, denoted as frpos and frneg respectively. Dangling bonds in both moieties are saturated. Thus, mutating frpos, and frneg leads to mfrpos,j, and mfr neg,k where j covers two mutation possibilities in the carboxylate, while k covers the three possibilities in guanidinium. Cartesian coordinates are preserved during the two mutation processes, namely, wt−→mwti and fr pos/neg −→mfr pos,j/neg,k . c) Distribution of Eip,i for symmetric and asymmetric salt-bridges, E intra ip is taken as the mean.
we hydrate our target molecules in a polarizable continuum medium and the geometries are optimized to capture reorganization energies. When a molecule with intramolecular ion-pair interactions is relaxed with a fixed number of water molecules, the resulting fragments will also contain the same number of water molecules with same orientation as in wt/mwt.
Combining (1) and (2) we arrive at for the i-th mutation path. The approximation holds as long as E al. ≈ E al. which can be expected when the interacting moieties are appropriately selected to adequately capture polarization effects on the linker atoms. The final value of E intra ip for a given geometry of wt is estimated by averaging over all possible alchemical muta- The standard deviation with respect to the mean indicates asymmetry of the ion-pair bonding pattern as shown schematically in FIG.1c. Hence, an asymmetric bonding pattern will result in distinct E ip values whose standard deviation corresponds to the method's uncertainty. For symmetric bonding patterns, the standard deviation vanishes as all mutation paths are equivalent.

III. COMPUTATIONAL METHODS
For four large peptides, 1EJG, 1BDK, 1L2Y, and 1SCO, we collected geometries from PDB. For all other systems, we performed geometry optimizations using the ωB97X-D3-DFT [65,66] method with the def2-TZVP basis set. This method was used because of its ability to converge to zwitterionic molecules while preserving the intended bond connectivities [67]. Solvent effects (water medium) were modeled using the conductor-like polarizable continuum medium (CPCM) approach [68,69]. Explicit solvent interactions were modeled using the microsolvation approach where 2, 4, or 8 explicit water molecules were included in the DFT/CPCM geometry optimizations. In microsolvation investigations, different conformers were generated starting with randomly sampled arrangements of the water molecules. For a given solute, the center of the centroid-centroid distance of the polar moieties involved in non-covalent ion-pair interaction was chosen as a center of a sphere on whose surface the water molecules were randomly dispersed. This implies that in the initial geometries generated, the water molecules were close to the ion-pair containing end of the molecules. Both solutes and solvents were then allowed to undergo geometry relaxation. In vacuum phase geometry optimizations, N-H bond lengths were constrained at the cationic terminals in order to converge to an ionpair local minimum and prevent converging to a thermodynamically more stable proton-transferred structure. Such constraints were not required with CPCM and microsolvation for most of the systems studied here. All single-point calculations were performed on the ωB97X-D3/def2TZVP geometries. For the aforestated large peptides, such calculations were done using the PDB geometries. DLPNO-CCSD(T) [62,63] calculations were performed with the aug-cc-pVTZ basis set and TightPNO settings. DFT and DLPNO-CCSD(T) calculations were performed using ORCA 5.0 [70,71] with the resolution-ofidentity (RI) approximation [72,73] for Coulomb (J) and 'chain-of-spheres' (COS) algorithm for exchange integrals (RIJCOSX). For RI calculations we used the Weigend auxiliary basis sets [74]. DFT calculations were done with the default grid, defgrid2. Intramolecular ion-pair interaction energies were modeled with F/ISAPT0 [53] using aug-cc-pVTZ basis set. This basis set was selected after comparing its accuracy with that of jun-cc-pVDZ for modeling SAPT0-level intermolecular molecular ionpair interaction energies. All SAPT calculations were performed with the code Psi4 [75].

IV. RESULTS AND DISCUSSION
In FIG. 1, we discussed the influence of an ion-pair interaction on the net stability of the compact conformer over the long conformer. While the compact conformer is most stable when its polar groups are suitably oriented for maximal interaction, this comes at the cost of internal strain. In large molecules, this balance between strain and ion-pair bonding along with myriad other non-covalent interactions determine the global minimum's structure [1]. Although studies have quantified the effect of interactions between polar groups from a super-molecule approach [38,90], or as relative free energy changes using molecular dynamics simulations [91,92] it is of interest to directly estimate such intramolecular interaction energies for a given structure.
In this study, we pursue this goal with F/ISAPT0 and a thermodynamic cycle. To benchmark these methods, we design two sets of model systems. The first set comprises nine dimers stabilized by intermolecular ion-pair bonding. The interaction energies of these molecules were estimated in a supermolecule fashion with the CCSD(T)/aug-cc-pVTZ method. Using these values as reference, we benchmarked SAPT0/jun-cc-pVDZ and SAPT0/aug-cc-pVTZ methods, and found the augcc-pVTZ to be more suitable for SAPT modeling of ionpair systems. Subsequently, we modeled intramolecular ion-pair interactions in a second set with 26 molecules using the aug-cc-pVTZ basis set in combination with the F/ISAPT0 method. For this set, we also benchmark E ip (Eq. 3) determined with the thermodynamic cycle based on different DFT methods. With the best set up, we report salt-bridge interaction energies of model tripeptides and large proteins in solvent medium.

A. Model ion-pair systems
In FIG. 2 we describe the composition of model ion-pair systems considered in this study. The polar groups comprises three cationic (C n ), three anionic (A n ) moieties that are biochemically relevant [90]. The diverse resonance structures of these polar moieties facilitate various bonding patterns as seen, for example, in guanidinium-carboxylate salt-bridges [32]. The cation set includes ammonium (C 1 ), formamidinium (C 2 ), and guanidinium (C 3 ) moieties with 1, 2 and 3 resonance structures, respectively, while the anion set includes oxide (A 1 ), carboxylate (A 2 ), and sulphate (A 3 ) moieties with a similarly increasing number of resonance structures.
The fist set of model systems used for quantifying E inter ip were formed by attaching methyl groups to the polar moieties, directed away from the ion-pair bonding. This set contains 9 dimers formed by combining 3 cations with 3 anions. For each system, we estimated E inter ip for 3 different centroid-centroid distances, resulting in 27 energies. The second set of model systems to quantify E intra ip comprise these same polar moieties by covalently connected via three different linker (L n ) fragments, n = 1, 2, 3. For linkers, we chose linear alkane chains with 3 (L 1 ), 6 (L 2 ), and 9 (L 3 ) methylene units. These non-polar linkers are not expected to interfere with the salt-bridging interactions. Of the 27 possibilities, C 2 -L 1 -A 1 was ignored because during geometry optimization it cyclized to an amide.

B. Basis set for SAPT modeling of ion-pair interactions
In Table. I we report SAPT0-level E inter ip using jun-cc-pVDZ and aug-cc-pVTZ basis sets for 9 model ion pairs with three different centroid-centroid distances. In these complexes, the polar moieties faced each other while the methyl groups were pointed outwards. Further, for C 1 -A 1/2/3 and C 2/3 -A 1 type complexes, we accessed the 3.0 Å, 5.0 Å, and 7.0 Å centroid-centroid distances, while the remaining complexes were placed at 4.0 Å, 5.0 Å, and 7.0 Å to avoid small interatomic distances. We compared the SAPT0 results with those from a supermolecule scheme based on CCSD(T)/aug-cc-pVTZ energies. For short centroid-centroid distances, the interaction energies are high. On an average, we note that the performance of SAPT0/jun-cc-pVDZ to be inferior to SAPT0/aug-cc-pVTZ, with the largest deviations noted at short contacts. Hence, in the remainder of the study, we use the aug-cc-pVTZ basis set for modeling E intra ip at the F/ISAPT0 level.

C. Intramolecular ion-pair interactions in model systems
For 26 model systems, we determine E intra ip using the thermodynamic cycle (Eq.3) and compare the results with F/ISAPT0/aug-cc-pVTZ values ( E int ISAPT ). The increasing size of the system with linker lengths rendered CCSD(T)/aug-cc-pVTZ level calculations prohibitively expensive. Hence, we resorted to the DLPNO variant of CCSD(T) for the thermodynamic cycle, and pare. Since the inherent limitations of both methods stem from different sources, at the limit of their mutual convergence Table I. Deviation of intermolecular ion-pair interaction energies of 9 dimer model systems estimated using SAPT0/augcc-pVTZ and SAPT0/jun-cc-pVDZ from CCSD(T)/aug-cc-pVTZ reference. The mean absolute deviation of both methods are reported in the last row. Values in parenthesis in the first column the centroid-centroid distances between the monomers in Å. All other values are reported in kcal/mol.

SAPT0 System
Reference aug-cc-pVTZ jun-cc-pVDZ their predictions can be expected to be closer to the reality. In Table. II, we report on the deviation in interaction energies (E int ISAPT ) determined by F/ISAPT0 and E ip using Eq. 3 with total energies from DLPNO-CCSD(T) and ωB97X-D3-DFT levels for various choices of terminal groups: -CH 3 , -CH 2 CH 3 , and -CH 2 CH 2 CH 3 .
In general, for all choices of terminal groups, and across all methods, the deviation between E int ISAPT and E ip [DLPNO-CCSD(T)] largely follows the order L 1 > L 2 > L 3 . This trend is reflected better in the standard deviation (RMSD) and mean percentage absolute deviation (MPAD) than the mean absolute deviation (MAD). In L 1 systems, the net ion-pair interaction is not only mediated through space but also through polarization of the short linker. This through-bond contribution is not accounted for in E int ISAPT . However, in the case of longer linkers (L 2 and L 3 ), the through-bond contribution to E int ISAPT becomes negligible justifying the MO localiza-tion scheme adopted. Hence, the F/ISAPT0 modeling of ion-pair interaction is most appropriate when the polar groups are separated by sufficiently long linkers.
In the thermodynamic cycle, we do not employ localized orbitals in order to describe the electron density as realistically as possible. Hence, when the spatial extent of the polar group is confined only until the nearest neighbor group (terminal group: -CH 3 ) and further through-bond polarization is not captured. Therefore, increasing the length of the terminal group from -CH 3 to -CH 2 CH 2 CH 3 results in improved agreement between E int ISAPT and DLPNO-CCSD(T)-based E ip with a mean deviation of 2.73 ± 2.49 kcal/mol (deviation of ≈ 2%) for the 9 C x L 3 A y molecules. The remaining discrepancies are presumably within the uncertainty of the model arising from averaging over mutation pathways (see FIG.1b).  Scatterplot comparison of intramolecular ionpair interaction energies, estimated with ISAPT, E int , and Eip, estimated using the thermodynamic cycle (Eq. 3) based on ωB97X-D3/def2-SVPD total energies. All values are in kcal/mol. Red triangles, green squares, and blue circles correspond to CxL1Ay, CxL2Ay, and CxL3Ay systems, respectively. For each case, Pearson correlation coefficients are also provided.
When replacing the DLPNO-coupled-cluster energies in Eq. 3 with DFT values, we see good agreement for basis sets with a diffuse function: def2-SVPD and def2-TZVPD. For the L 3 systems with the longest terminal group explored, the mean deviations of the ωB97X-D3-DFT based estimations of E ip from the ISAPT values are 3.25 and 2.97 for def2-SVPD and def2-TZVPD basis sets, respectively. Hence in further applications, we use the def2-SVPD basis set due to its favorable cost-accuracy trade-off.
In FIG. 3, we compare the absolute values of E int ISAPT and ωB97X-D3/def2-SVPD-based E ip for 26 bench-   mark ion-pairs. E ip values were estimated using the −CH 2 CH 2 CH 3 terminal group. Interaction energies computed using both methods lie in the range of -160 to -80 kcal/mol. Such large values are expected in the vacuum phase as noted previously in the EDA calculations of model ion-pair systems [93]. As noted in Table. II, with increasing linker length, the agreement improves as reflected by the Pearson correlation coefficients. The correlation is 0.98 when using a sufficiently long linker. In general, the magnitude of interaction energy is large when the polar groups are connected by a long and flexible linker resulting in a more compact structure. The best agreement between F/ISAPT0 and the thermodynamic cycle is seen for the 9 C x L 3 A y systems with the X/Y=CH 2 CH 2 CH 3 terminal groups, in Table. III. Hence, we benchmark the performance of 10 DFT approximations only for these systems. Several of the range-separated and long-range corrected hybrid-DFT methods show good agreement with the reference F/ISAPT0 values. All methods (barring LC-PBE) differ amongst each other within 1 kcal/mol suggesting the uncertainty due to the choice of the DFT method to be smaller than the inherent uncertainty in the thermodynamic cycle. Since ωB97X-D3 was used for geometry optimization, we continue with the same method in combination with the def2-SVPD basis set for the rest of the study.
Motivated by the good agreement with the DFT-based thermodynamic cycle and F/ISAPT0 gas-phase interaction energies, we now move on to see the effect of solvent on E ip . This is relevant because salt-bridging are mostly found in proteins, which in turn are often experimentally isolated in the aqueous phase. Alas, the ISAPTformalism has not been extended to solvent phase. However, there is no limitation in applying the thermodynamic cycle in the solvent phase. The importance of accurately modeling the aqueous phase to study noncovalent interactions [94,95], reactions [96], and ligandmetal interactions [97] has been well established. The magnitudes of E ip is expected to decrease in the aqueous phase where the ion-pair interaction is screened. To this end, using the vacuum phase geometries, we performed CPCM total energy calculations for all 26 benchmark systems. FIG. 4 presents E ip calculated at ωB97X-D3/def2-SVPD and ωB97X-D3(CPCM; water)/def2-SVPD levels. While the vacuum phase values are < −80 kcal/mol, in the CPCM phase even the strongest interaction has an E ip of about -40 kcal/mol. In general, the magnitude of vacuum phase E ip decreases by a factor of 3.6-9.3 in the CPCM phase. Further, the effect of CPCM on E ip preserves the qualitative trends in ion-pair strengths as noted in the vacuum phase. In both media, we note the interactions due to monodentate polar moieties (ammonium or oxide) resulting in the strongest ion-pair interaction. We explored the ion-pair interactions with longer terminal groups (−CH 2 CH 3 and −CH 2 CH 2 CH 3 ), and noted their influence on E ip to be minimal in the CPCM phase. FIG. 4 ascertains that our estimation of E ip is efficient in the continuum solvation paradigm capturing commonly expected physical trends. This motivated us to investigate E ip for some biologically relevant model systems with CPCM and with microsolvation modeling to account for explicit solute-solvent interactions.

D. Microsolvation effects on model tripeptides
Upon realizing the quenching effect of CPCM on E ip in the 26 model systems we expand our study towards biologically relevant molecules with implicit and microsolvation modeling. The functionality of most macromolecules is strongly influenced by their interaction with water molecules in the medium. For example, explicit solute-solvent interactions influence pharmacodynamic activities of ligand-protein complexes [98]. It has been shown that quantitatively accurate modeling of saltbridges realistic an adequate description of the solvent phase [99]. Since previous studies have shown that the zwitterionic form of amino acids is best described in the presence of water molecules [99,100], we modeled the four tripeptides-lysine-glycine-aspartic acid (KGD), lysine-glycine-glutamic acid (KGE), arginine-glycineaspartic acid (RGD), and arginine-glycine-glutamic acid (RGE)-with various degrees of microsolvation. Their C−terminals were amidated with methylamine to confine salt-bridging interactions to only between side-chains and to prevent the participation of main chain carboxylic groups [2].
In FIG. 5, we show the definitions of the polar terminal groups of for the four tripeptides studied here. We fragmented the positive amino acids by taking all atoms of the amino acid up to the C α atom of the adjacent residue, i.e. glycine in this case. Hence, the NH part of the K- G and R-G amide bonds were included as a part of the cationic fragments, K and R. Similarly, the CO part of the DG and EG amide bonds were included as a part of the anionic fragments, D and E. The dangling bonds were saturated by capping with an H atom. Hence, the terminal groups are now different for the two non-interacting components of fr (in FIG. 1b) and they depend on the exact sequence of amino acids running from the N to C-terminals. This implies that the terminal groups for fragments of KGD and DGK will be different. The coordinates of the atoms in these fragments were kept frozen as in the compact structures and the mutation paths were followed as discussed before in Sec.II.
In order to study the influence of water molecules on E ip , we inspected the 10 conformers of RGE in the presence of 4 water molecules in a background CPCM environment. While the influence of conformers on intramolecular H-bonding interactions has been studied before [101], here we wanted to study the conformer in-fluence on intramolecular salt-bridging interactions. We started with RGE relaxed in CPCM without any explicit water molecule as the initial structure to which we added water molecules. Hence, when the initial structure showed a bidentate side-on [32] interaction between guanidinium and carboxylate groups, it is surprising to observe in FIG. 6 a monodentate structure as the global minima. Previous studies have noted a better stabilization of solvent-separated ion-pair in lysine-glutamate dipeptide with an increasing number of explicit water molecules [99]. The remaining 9 complexes showed sideon interactions. We calculated the relative E ip values with reference to the most stable conformer with a sideon interaction.
While for all side-on conformers, the ∆E ip is very similar, for the monodentate complex it is found to be destabilizing. Inspecting this structure revealed that this is due to an unfavorable orientation that diminishes the salt-bridging interaction. The net stability of the global minimum is due to solvation energy dominating over the salt-bridging contribution. It may be realized that for a highly stable complex, the water molecules and the polar groups must involve in maximal bonding (including H−bond network among the water molecules). The role of solvation in stabilizing a system decreases when the H−bonds are formed with weakly polar groups or when the water network is disconnected. In FIG. 6, we note this effect on the conformers on the right side, where the complexes destabilize when all solvent molecules do not form extended networks with the polar groups. Therefore, it may be concluded that, while favorable orientations with water molecules may benefit the overall stability of the complex, it comes at the cost of quenched salt-bridging interactions.
Varying the number of water molecules in the microsolvation of tripeptide will provide insights into the effect of solvation in the bulk water scenario. In FIG. 7, we plotted E ip for all 4 tripeptides in the presence of 0, 2, 4, and 8 water molecules. The zero-water case corresponds to the CPCM limit with one minimum energy structure. For the case of 2, 4, and 8 water complexes, we sampled 10 conformers starting with different arrangements of water molecules, of which we report only the global minimum's E ip regardless of the bonding denticity. Thus, for all 4 peptides, increasing the number of water molecules mostly results in quenched salt-bridging interactions. This is a physically expected trend that has been captured by the model. However, the degree to which the E ip changes is different across the peptides. This can be attributed to the relatively small conformer space sampled and due to the finite degree of hydration.
A closer look at the structures sheds more light on the magnitude of the interaction energies reported in FIG. 7. For KGD, with an increasing number of water molecules, the interactions between the polar fragments diminish. While this trend is seemingly similar in KGE, in the case of 8 water molecules KGE adopts a solvent separated ion-pair form, while that of KGD is solvent-shared with a Conformers are labelled from 0-9 (0 being the global minimum). For clarity, structures have been arranged above and below the horizontal axis. The leftmost structure is with a monodentate salt-bridge bonding where only one of the two O atoms of the carboxylate moiety binds to both NH atoms (see inset) while the structure right to it shows a bidentate interaction (see inset).
stronger ion-pair interaction. For RGD, all solvated complexes prefer the monodentate arrangement. Hence with increasing degree of hydration, E ip of RGD approaches a limit of > −10 kcal/mol. In the case of RGE, 0, 2, and 8 water complexes form a side-on bidentate arrangement, while with 4 water molecules, the most stable conformer has a monodentate arrangement with an elevated E ip . It may be noted that the linker length is longer in RGE compared to RGD favoring a more stable side-on saltbridging. Hence, by comparing the bonding pattern of the minima we can qualitatively interpret E ip values estimated using the thermodynamic cycle.

E. Application to large peptides
We collected the NMR resolved structures of crambin (1EJG) [102], bradykinin (1BDK) [103], tryptophan zipper I (1L2Y) [104], and a scorpion toxin OSK1 (1SCO) [105]. Crambin (1EJG) is a neutral plant protein with 46 amino-acid residues and 642 atoms stabilized by the salt-bridging interaction between a carboxylate group on the main chain of Asn46 and a guanidinium group of Arg10. 1BDK is 11 residues long with 187 atoms with a net charge of +3. 1L2Y has a salt-bridging interaction between Arg16 and Asp9, containing 304 atoms and a net charge of +1. 1SCO is a scorpion toxin with an overall charge of +8 and 595 atoms, presenting a i-to-i+4 salt-bridge between Arg12 and Glu16. These 4 peptides provide an opportunity to understand the significant role orientation plays in deciding the strength of a salt-bridge. Their structures are on display in FIG. 8. In these proteins, centroid-centroid distances of ion-pairs are < 5.0Å. 1EJG shows a side-on salt-bridge (Asn46-Arg10) with a short centroid-centroid distance leading to a significantly large magnitude of E ip , -24.7 kcal/mol. 1BDK is stabilized in a monodentate fashion with a relatively weaker interaction energy than 1EJG. This is followed by 1L2Y presenting an even more diminished interaction. Finally, in 1SCO we note a small positive E ip .
A closer inspection of the individual E ip values obtained through various mutation paths provide us with further insight. As noted before in FIG. 1c, a salt-bridge with a symmetric bonding pattern should show a smaller variance across mutation paths. Hence, for 1EJG, where the salt-bridge is symmetric, all E ip obtained are of similar magnitude with a small standard deviation. However, for monodentate arrangements in 1BDK and 1L2Y, the symmetry is lost and we note a larger variation in the E ip . For 1SCO, while most mutation paths yield small positive values, there is one path that yields a small negative value. Hence, the salt-bridging in this system is es-sentially non-interacting, if not mildly unfavorable. This could be speculated to be due to the unfavorable orientation of the interacting moieties. Despite the relatively short centroid-centroid distance, the polar groups are oriented away from each other, resulting in a net repulsion, unlike all the other cases explored here. The small positive values arise when we explore mutation paths that involve these unfavorably oriented atoms. This example demonstrates how a salt-bridging with a reasonably small centroid-centroid distance can be overall irrelevant.

V. CONCLUSIONS
In this study, we present a thermodynamic cycle to quantify intramolecular ion-pair interactions between two polar moieties of a molecule, without complex bifurcating interactions. We benchmark the efficiency of this approach by comparing intramolecular interaction energies obtained for a set of 26 zwitterions stabilized by ion-pair interactions against predictions from ISAPT. We found the deviations between the two approaches to be minimal for molecules containing long linkers and when the  terminal groups are adequately capped. Going from the vacuum phase to the aqueous phase, the ion-pair interactions are quenched as expected.
We modeled microsolvation of four biologically relevant tripeptides and found increasing degree of solvation to quench their E ip . In the presence of solvent molecules, the weakest interaction was found for a solvent-separated salt-bridge where the polar groups are surrounded by a cage of water molecules. The salt-bridge interaction is stronger for conformers with a bidentate side-on bonding. On the other hand, in conformers with a monodentate bonding, the salt-bridge is more exposed to solvent molecules resulting in better hydration of the biomolecule. In the aqueous environment both hydration and salt-bridging are competing stabilizing factors, the latter favoring a folded structure. We applied the thermodynamic cycle on experimental structures of 4 large proteins and found the salt-bridging interaction to be strongest for 1EJG with a side-on bonding, while 1BDK and 1L2Y with a monodentate bonding showed weak interactions. In 1SCO, even though the distance between polar groups suggests a possible salt-bridge bonding, due to the unfavorable orientation of the corresponding moieties, its E ip turned out to be insignificant. The favorable computational overhead of the proposed method promises explorations of salt-bridging interaction in other biological systems in the solvent phase. Since the present model can only account for simple, non-bifurcated saltbridges, future endeavors may concentrate on quantifying interactions for complex salt-bridges.

A. Conflicts of Interest
The authors have no conflicts of interest to disclose.

B. Data Availability
The data that support the findings of this study are openly available in Github, Refs. 106.