Physical Chemistry Chemical Physics

Thermodynamics of stacking interactions in proteins

Simone Marsili^ab, Riccardo Chelli^ab, Vincenzo Schettino^ab and Piero Procacci*^ab
aDipartimento di Chimica, Università di Firenze, Via della Lastruccia 3, I-50019, Sesto Fiorentino, Italy
^bEuropean Laboratory for Nonlinear Spectroscopy (LENS), Via Nello Carrara 1, I-50019, Sesto Fiorentino, Italy. E-mail: procacci@chim.unifi.it

Using a database of 6166 experimental structures taken from the Protein Data Bank, we have studied pair interactions between planar residues (Phe, Tyr, His, Arg, Glu and Asp) in proteins, known as – interactions. On the basis of appropriate coordinates defining the mutual arrangement of two residues, we have calculated 2-D potentials of mean force aimed at determining the stability of the most probable structures for aromatic–aromatic, aromatic–cation and aromatic–anion bound pairs. Our analysis reveals the thermodynamic relevance and the ubiquity of stacked complexes in proteins.

Stacking interactions between systems play a major role in molecular recognition,^1–3 in shaping the structure and the functionality of nucleic acids⁴ and in determining the tertiary and quaternary structure of proteins.^5,6 The physical nature of these interactions is now quite well understood. Vertical stacking between planar molecular systems is favored in polar environment, mainly due to entropic reasons, while electrostatic forces are in general determinant for planar or T-shaped structures in non polar media. In this respect, the traditional terminology – stacking is somewhat misleading since it is suggestive of some kind of special interaction between molecular orbitals. This speculation, still very popular today, is supported by the fact that vertical stacking is indeed quite frequent in proteins,^5–7 suggesting in effect the intervening of a sort of – bonding interaction.

Recent accurate ab initio studies^8–10 have actually shown that electron density reorganization upon formation of vertical stacked complexes of neutral monomers is negligible and even less important than that found in T-shaped structures. One fortunate consequence of this fact is that stacked interactions of aromatic systems can effectively be modeled^10–12 using molecular mechanics approaches with non polarizable force fields. In this respect, molecular dynamics simulations of aromatic systems in polar solvents have consistently shown that the hydrophobic effect together with the compactness of stacked planar geometry plays an important role in stabilizing these structures.^10,12,13 In other words, because of the planar geometry, the vertical stacked configuration is the most effective way of reducing the surface area of the solute cavity, allowing for a surprisingly high entropy gain of the surrounding solvent molecules. Indeed, according to a recent molecular dynamics study,¹³ already two water molecules are sufficient to change the coplanar arrangement of an H-bonded adenine–thymine base pair into the stacked one. Even more impressively, computer simulations^14,15 have demonstrated that vertical stacked complexes made up of two monomers with like charges are stabilized by entropy gain in polar media. While both experimental and theoretical studies have by now well characterized the chemical physics of the stacked interactions of planar and aromatic systems, their precise role in determining the protein tertiary and quaternary structure is still a matter of discussion and intensive research.

systems such as Phe or Tyr have been experimentally found to promote amyloid fibril aggregation^16–18 in small peptides^19,20 as well as in proteins.¹⁶ The role of – stacking in this promotion is, however, still controversial.¹⁶ In model systems^17,21 it has been also recently observed that – and –cation (either Lys or Arg) interactions may act as cooperative driving forces in the formation of salt-bridges on the protein surface. Based on this and other experimental evidence on sequence-designed polypeptides, it has been further speculated that such concurrent interactions are likely to be important in amyloidogenesis.^17,22

A possible additional role of -cation and -anion³ stacking in the stability of tertiary and quaternary structures of proteins can be envisaged in the light of some recent experimental data concerning the precise nature of the hydration structure of chaotropic planar denaturant ions such as the guanidinium ion²³ or the SCN⁻ ion and on the anionic penetration of lipid monolayers.^23,24 According to these studies, the Hofmeister series²⁵ can be rationalized in terms of direct ion–macromolecule interactions rather than in terms of disruption of the H-bond network of the surrounding water.²⁶ In the specific case of the guanidinium and SCN⁻ ions,²³ an entropic force combined with the possibility of formation of compact vertical stacked structures aimed at exposing partially buried planar arrangements could be important in the unfolding process of a macromolecule.

Whatever the role of – and –ion stacked interactions in proteins is, either in connecting distinct secondary structure elements, or promoting amyloid aggregation, salt bridge formation or in functioning as denaturing hot spots, there should be detectable thermodynamic signatures of these interactions in the available protein structures. With regard to this, several systematic studies on the incidence and thermodynamics of stacked interactions have been conducted in the past on limited sets of non-homologous proteins.^27–29 Due to the poor sampling, conclusions on the relevance of – and –cation interactions and assessment of stabilization free energies were, in general, quite contradictory. The results from statistical analysis of structural data appeared to have finally settled some time ago when various studies on the Protein Data Bank³⁰ (PDB) involving 505⁵ and 1671¹⁰ non-homologous structures came to the same conclusions, i.e., stacking is stable in proteins and the stabilizing energy was consistently found in the order of a few kcal mol⁻¹ for most analyzed pairs. A further very recent investigation³¹ on about 2000 non-homologous proteins taken from the PDB was focused mainly on the geometrical rather than the thermodynamical aspect of non-bonded interactions involving planar groups. The 1998 study by McGaughey et al.⁵ dealt with a PDB including a total of 400 folds, 650 superfamilies, and 950 families.³² By the beginning of 2002, the year of the survey of ref. 10 on – interactions, the numbers in the SCOP³² classification statistics had almost doubled with 650 folds, 1100 superfamilies, and 1700 families. On June 2005,³¹ the PDB again almost doubled, numbering about 1000 folds, 1600 superfamilies, and 3000 families.

In the present article we construct a new expanded database of structures based on the PDB as of June 2007 with emphasis not only on intraprotein contacts (non-homologous proteins), but also on protein–protein interactions for examining the thermodynamics of – stacking in determining the protein tertiary and quaternary structure. In a Boltzmann-like approach, the expanded protein database used in the present study allows for a high grade of sequence homology in such a way that the most populated protein families have a correspondingly larger weight. Concerning the statistics of stacking interactions, our protein database and two of the most popular databases of non-homologous proteins^33,34 are in excellent agreement (see the ESI). In order to have a complete and clear thermodynamic picture, we also systematically integrate the energetics of – interactions with data relative to the mean solvent exposure of stable complexes.

The outline of the article follows. In section II, we discuss the methodology adopted in the present study, by characterizing the analyzed protein database and by describing the theoretical aspects related to the calculation of 2-D potentials of mean force (PMFs) of complexes and of the solvent exposure of stable pair configurations. In section III, we present the main results, namely the 2-D PMFs of aromatic–aromatic, aromatic–cation, and aromatic–anion residue pairs expressing quantitatively the thermodynamics of stacking in proteins. Conclusive remarks are given in section IV.

The PDB, as of 18 June 2007, contained 44095 structures. Since we are interested in the overall thermodynamics of – interactions, i.e. in their relevance for determining both the tertiary and quaternary structures in proteins, contacts between residues belonging to the same and to different subunits (namely amino acid chains with distinct chain-id) have been taken into account in the statistical analysis. Therefore, in order to construct an extended and representative set of protein structures, we sorted all available structures in the PDB eliminating structural degeneracy. In particular, two structures (PDB files) have been considered identical or degenerate if all following criteria are fulfilled: (a) the structures have the same name (the PDB keyword HEADER); (b) the structures have the same number of subunits; (c) in the two structures, the numbers of amino acids of two corresponding subunits (selected following the subunit length in descending order) differ by less than four. If structural degeneracy occurs, then only one structure, chosen randomly, is retained in the database. Application of the above criteria reduces the database to 9280 structures. These structures have then been processed in order to discard PDB files with only C-traces and PDB files containing defective amino acids or non-standard residues, finally leading to a database of 6166 structures (the full list of the proteins in the database is available upon request). The average number of subunits per protein is 2.9, the average number of amino acids per protein is 874, and the average number of amino acids per subunit is 297. The composition of the database in terms of secondary structure yields 28.9% -helix and 16.5% -sheet.

In this article we examine the thermodynamics of pair interactions of planar residues in proteins, with emphasis on stacking. The – interactions that we have analyzed involve the residues R = Phe, Tyr, His, Arg, Glu, Asp. There are many possible choices for studying the incidence of – arrangements in proteins.^5,31 Here we follow a well established procedure¹⁰ based on the definition of a PMF for a restricted set of coordinates, assuming that the collection of pairs of planar residues in the database reflects an equilibrium distribution at ambient temperature and pressure.

To this end, we represent each residue by two vectors in space, r^R_i and n^R_i, where R refers to the amino acid type and 1 i N_R, with N_R being the number of residues of type R in the protein database. For R = Phe, Tyr, His, the vector r^R_i is the position of the geometric ring center. For the other residues, r^R_i is the vector position of the carbon atom of the guanidinium (Arg) and of the carboxy group (Glu and Asp). n^R_i corresponds to the normal to the molecular plane of the ith residue of type R. For a pair of residues of type R and R, labeled i and j, respectively, we define the following coordinates:

(2.1)

(2.2)

(2.3)

	Fig. 1 Left: visualization of the coordinates defining the mutual position and orientation of the planar moieties (see eqns (2.1), (2.2) and (2.3)). Right: schematic representation of the most common stacked and T-shaped pairing arrangements of the planar moieties in the protein database.

In using the definitions of eqns (2.1), (2.2), and (2.3), not only do we assume a D_h planar symmetry of the monomers, but for 0 we also choose not to distinguish between the non equivalent structures obtained when the normal to the plane of one monomer precedes around the normal to the plane of the other monomer. Nonetheless, for short distances r, these non equivalent structures all bearing the same coordinates are dominated by few attractors at specified , angles whose structure may often be deduced judiciously, taking into account the chemical topology of the dimers. Moreover, as approaches zero, i.e., for stacked structures, all non equivalent structures become degenerate given the assumed D_h symmetry of the monomers. Schematic pictures of these stacked arrangements with related classification are shown in Fig. 1.

We define an unnormalized pair correlation function in the space of the restricted set of variables r, , and

(2.4)

_RR

_RR

_RR

_RR

RR(r,,) = gRR(r,,) + gRR(r,,).

(2.5)

_RR

_RR

_RR

The PMF in the r, , space can be calculated (up to an addictive constant) from the pair correlation function of eqn (2.5) as follows

wRR(r,,) = −RT ln RR(r,,),

(2.6)

(2.7)

_RR

_RR

_RR

_RR

For a more detailed analysis of the structural features of pair stacked arrangements, we define (up to addictive constants) two 2-D PMFs as

(2.8)

(2.9)

_RR

_RR

We stress that the thermodynamic potentials w_RR(r,,), (r,) and (r,) reflect any sort of deviation from the perfectly random probability. Therefore, their topology will retain the characteristics of the residue pair RR of enthalpic and entropic origin (which are the main object of the present study), but also those related to general topological features due, for example, to the primary and secondary structure of proteins.

In order to highlight the connection between the function g_RR(r,,) and the so-called centroid–centroid contact potential,^6,35,36 we point out that the conditional probability that two residues are in contact given that they are of type R and R is

(2.10)

_RR

_RR

The solvent accessible surface of a generic residue R is calculated on the planar moieties of the side chains using the analytical method proposed in ref. 37, using van der Waals radii from ref. 38 and a solvent particle radius³⁷ of 1.4 . A residue or a bound dimer is considered to be exposed if a surface of at least 12 ² is exposed, independently of the size of the residue or dimer. This area corresponds approximately to the size of the largest section of a sphere with a volume corresponding to the mean volume of a water molecule in pure liquid at normal conditions. Given a class or ensemble of molecular systems denoted with C (pairs of residues, single residues or atoms), the probability P^(C)_E of any member of the class C being exposed is given by P^(C)_E = _E/, where _E and are the number of exposed members and the total number of members (exposed or not) in the class C, respectively. The P^(C)_E values for the residues considered in the present study are the following: 0.36 (Phe), 0.57 (Tyr), 0.68 (His), 0.80 (Asp), 0.81 (Arg), and 0.86 (Glu).

In the next section we will report and discuss the solvent exposure of bound residues, focusing on the differences arising from their pair arrangement, i.e., stacked or T-shaped. For each pair of residues RR, solvent exposure for stacked arrangement is evaluated over the class made up of pairs falling into the S domain, while for the T-shaped arrangement, pairs belonging to the T domain are considered (see above for the definitions of the S and T domains).

In this section we address the stabilization free energies of stacked and T-shaped configurations in the protein database, characterizing at the same time these pair interactions in terms of solvent exposure. Stable dimeric structures are identified through the analysis of the 2-D (r,) and (r,) PMFs (eqns (2.8) and (2.9), respectively). A further analysis is conducted by calculating the excess free energy of stacked () and T-shaped () structures in the space of the r, , coordinates (see eqn (2.7) and related discussion). We divide – interactions in three classes which will be reviewed separately. The first class, aromatic–aromatic, includes the Phe–Phe, Phe–Tyr, Phe–His, Tyr–Tyr, Tyr–His and His–His pairs. The second class, aromatic–cation (or aromatic–Arg), includes Phe–Arg, Tyr–Arg and His–Arg pairs. For reasons that will be apparent later on, in this class we also include the Arg–Arg pair. The last class, aromatic–anion, includes Phe–Glu, Tyr–Glu, His–Glu, Phe–Asp, Tyr–Asp and His–Asp pairs.

In Fig. 2 and 3 we show the (r,) and (r,) PMFs for pairs of the aromatic–aromatic class. Overall, (r,) shows two competing minima corresponding to parallel and orthogonal arrangements. Parallel configurations are localized at low values (see for example the series of minima going from 3.5 to 5 in the w_FF(r,) PMF), while the minimum concerning the orthogonal structures, if any, is present at 90° and r 5 . In the (r,) PMF, i.e., for parallel structures, two features are common to all aromatic–aromatic pairs: (i) the presence of an extended banana shaped minimum region corresponding, for < 10°, to vertical stacked configurations, while for greater values, to a manifold of parallel displaced structures; (ii) the existence of an energetically forbidden region at small values ( < 30° and 5.5 < r < 7 ). This latter feature is a consequence of the tight packing in proteins: if two planar residues were at r and in the high free energy region, then no space would be left to allocate a third residue and a cavity ought to be formed. These features are suggestive of a minimum free energy path for stacking formation involving an approach at large as in an intercalation process.

	Fig. 2 Contour plots of the (r,) (eqn (2.8)) and (r,) (eqn (2.9)) PMFs for Phe–Phe, Phe–Tyr and Phe–His pairs (left and right sides, respectively). On the right of each picture, the chromatic scale refers to the value of the PMF expressed in kJ mol−1. The bin size of the and angles is 10°, while that of the distance r is 0.2 .

	Fig. 3 Contour plots of the (r,) (eqn (2.8)) and (r,) (eqn (2.9)) PMFs for Tyr–Tyr, Tyr–His and His–His pairs (left and right sides, respectively). On the right of each picture, the chromatic scale refers to the value of the PMF expressed in kJ mol−1. The bin size of the and angles is 10°, while that of the distance r is 0.2 .

The salient data of all possible aromatic–aromatic T-shaped and stacked dimers (see definitions of T and S domains in section IIB) are reported in Table 1. In general, the relative importance of T-shaped arrangements with respect to stacked ones decreases with an increase in the product of the solvent exposure of the two monomers (P^(R)_E from Table 1). This result is consistent with (and at the same time validates) the picture that emerges from theoretical studies of models of planar amino acid interactions,¹⁰ according to which T-shaped configurations are stabilized in non polar or weakly polar environments (e.g., the core in globular proteins), while stacked arrangements are prevalent in a highly polar medium (the protein surface). Also, the presence of a stable T-shaped structure, which is characterized by the largest degeneracy, makes the associated contact potential (see eqn (2.10)) more stabilizing⁶ with respect to other dimers that show only the stacked minimum, as the aromatic–cation and aromatic–anion pairs (see sections IIIB and IIIC).

RR					P(R)E
FF	−6.8	−7.7	0.25	0.31	0.36	0.36
FY	−5.7	−6.9	0.53	0.53	0.36	0.57
FH	−7.9	−6.0	0.76	0.66	0.36	0.68
YY	−7.4	−6.0	0.79	0.70	0.57	0.57
YH	−9.1	−5.2	0.91	0.73	0.57	0.68
HH	−9.7	−3.2	0.81	0.89	0.68	0.68

For the Phe–Phe (FF) pair, the T-shaped structure is found to be more stable than the stacked one by 1 kJ mol⁻¹ (see W^S_FF and W^T_FF in Table 1). This is expected on the basis of the solvent exposure of the corresponding monomeric units, Phe being the most hydrophobic (i.e., the less exposed) among the residues considered here. The stability of the T-shaped geometry makes the contact potential (see eqn (2.10)) for Phe–Phe interaction very stabilizing with respect to other aromatic–aromatic contact potentials, in agreement with the contact potential reported in ref. 6. Among the parallel structures, the parallel displaced one (20° < < 30°) is the absolute minimum (see w_FF(r,) in Fig. 2), with a significant presence of the vertical stacked minimum. Interestingly, we note a non-negligible contribution of a parallel minimum at longer distances (r 4.8 ) and larger values (40° < < 50°). These parallel structures, many of them having the C carbon of one Phe residue interacting with the ring of the partner (see structure A in Fig. 4), are enthalpically stabilized in a non-polar environment (note, for example, the structures FF1 and FF2 in Fig. 2 of ref. 10).

	Fig. 4 Representation of some relevant parallel configurations discussed in the text: Phe–Phe (A), Phe–His (B) and Phe–Glu (C). The Phe–Phe structures have 40° < < 50°, the Phe–His structures have 30° < < 40° and the Phe–Glu structures have 20° < < 30°. Phe–Phe is taken from the PDB file 1a76 (residues F8 and F235), Phe–His is taken from the PDB file 1dgr (residues H102C and F176C) and Phe–Glu is taken from the PDB file 1ey3 (residues F213B and E221B).

A very similar pattern, both in the PMF and in the solvent exposure, is observed for the Phe–Tyr (FY) pair, which also shows competition between a stacked minimum (W^S_FY = −5.7 kJ mol⁻¹) and a more stable T-shaped complex (W^T_FY = −6.9 kJ mol⁻¹). In the w_FY(r,) PMF, the parallel displaced minimum (r 3.9 , 20° < < 30°) prevails, while there is no clear sign of a vertical stacked minimum. Like in the Phe–Phe case, also for Phe–Tyr the driving factor in the overall exposure of the stacked bound complexes is the small exposure of the Phe residue (i.e., entropic forces do not appear to be dominant in the stabilization of the bound stacked Phe–Tyr complex).

For the Phe–His (FH) complex, the parallel displaced minimum is again the most significant. The stacked structure is the most stable in the aromatic–aromatic class, with a stabilization free energy W^S_FH = −7.9 kJ mol⁻¹. The w_FH(r,) PMF shows an extended region of stability with a typical banana shape although somewhat narrower than Phe–Phe and Phe–Tyr. Other secondary minima can be noted at r 4.3 and 30° < < 40°, where the C atom of the His residue lies over the Phe ring (see structure B in Fig. 4). The very shallow T-shaped minimum (which still contains more pairs than the deep stacked minimum due to the solid angle degeneracy; see Table 1) makes the contact potential (eqn (2.10)) for Phe–His interaction much less stabilizing with respect to Phe–Phe, in agreement with the contact potential reported in ref. 6. The solvent exposure of the bound complex appears to be dictated by the exposure of His rather than that of Phe. The His residue is in fact significantly more exposed than both Phe and Tyr and the stacked complexes are systematically more exposed than the T-shaped structures (see Table 1), indicating for the former and the latter a stabilization of entropic and enthalpic nature, respectively.

For the Tyr–Tyr (YY) pair, the excess free energy of stacking (W^S_YY = −7.4 kJ mol⁻¹) is higher than that of Phe–Phe and Phe–Tyr pairs, while we note a significantly smaller stabilization of T-shaped structures (W^T_YY = −6 kJ mol⁻¹) with respect to the same pairs. This observation is consistent with the weak contact potential for this dimer.⁶ In the w_YY(r,) PMF of Fig. 3 we notice that vertical stacked and parallel displaced structures have similar stabilities. Moreover, parallel minima with large are much less important for this pair than for pairs involving Phe. This fact, along with the stability of the vertical stacked minimum, suggests that the minimization of the exposed surface of the stacked dimer is a key factor in the stabilization of the complex. Accordingly, the stacked Tyr–Tyr structures are found on the average much more solvent exposed than the unbound pairs (see Table 1), indicating that stabilization of stacking is mainly of entropic origin.

For the Tyr–His (YH) pair, the excess free energy of stacking W^S_YH is −9.1 kJ mol⁻¹, mainly due to parallel displaced arrangements (see w_YH(r,) PMF of Fig. 3). The w_YH(r,) PMF shows a shallow T-shaped minimum (W^T_YH = −5.2 kJ mol⁻¹), which implies a less stabilizing contact potential. Overall, the stacked structures are found much more exposed than the T-shaped ones, suggesting that entropic thermodynamic forces are the main stabilizing factor.

The His–His (HH) pair shows the deepest stacked minimum in the aromatic–aromatic class with excess free energy W^S_HH of −9.7 kJ mol⁻¹. At the same time, the His–His pair exhibits a very shallow and much less stabilizing T-shaped region (W^T_HH = −3.2 kJ mol⁻¹) with respect to the pairs considered so far. Table 1 shows that stacked configurations are in general more exposed than the T-shaped ones, indicating for the former a stabilization mainly from entropic forces (like for Phe–His, Tyr–Tyr and Tyr–His pairs). The entropic stabilization of the stacked configurations is reflected in the w_HH(r,) PMF by the contraction of the banana-shaped region typical of aromatic–aromatic interactions.

In Fig. 5 and 6 we show the (r,) and (r,) PMFs for the aromatic–Arg class and Arg–Arg pair, respectively. The excess free energies and the solvent-exposure probabilities for stacked and T-shaped dimers are summarized in Table 2. As a general consideration on the (r,) PMF, we notice that this type of pair interactions is consistently characterized by a very stable stacked minimum. Accordingly, for the stacked configurations, the excess free energy is always lower than −8.7 kJ mol⁻¹, while for the T-shaped configurations it falls above −3.5 kJ mol⁻¹. From Table 2, we see that the solvent exposure of the complex is in general driven by the Arg exposure.

	Fig. 5 Contour plots of the (r,) (eqn (2.8)) and (r,) (eqn (2.9)) PMFs for Phe–Arg, Tyr–Arg and His–Arg pairs (left and right sides, respectively). On the right of each picture, the chromatic scale refers to the value of the PMF expressed in kJ mol−1. The bin size of the and angles is 10°, while that of the distance r is 0.2 .

	Fig. 6 Contour plots of the (r,) (eqn (2.8)) and (r,) (eqn (2.9)) PMFs for Arg–Arg pair (left and right sides, respectively). On the right of each picture, the chromatic scale refers to the value of the PMF expressed in kJ mol−1. The bin size of the and angles is 10°, while that of the distance r is 0.2 .

RR					P(R)E
FR	−8.7	−3.5	0.87	0.79	0.36	0.81
YR	−9.7	−2.9	0.89	0.83	0.57	0.81
HR	−9.1	−1.5	0.86	0.60	0.68	0.81
RR	−9.6	0.0	0.91	0.97	0.81	0.81

The Phe–Arg (FR) complex shows a w_FR(r,) PMF dominated by a deep minimum corresponding to stacked arrangements whose excess free energy is −8.7 kJ mol⁻¹. The w_FR(r,) PMF is characterized by a minimum including both vertical stacked (with the guanidinium carbon atom insisting on the centroid of the Phe ring) and parallel displaced configurations. The stacked minimum includes structures that are significantly more exposed than those of T-shaped bound dimers, again indicating a stabilization of entropic origin. A quick survey on a sample of vertical stacked structures reveals that the Arg residue systematically protects Phe from being totally exposed to the solvent, resulting in significantly larger solvent exposure of the complex with respect to that of the unbound monomers.

For the Tyr–Arg (YR) and His–Arg (HR) pairs, the pattern of the PMFs is very similar to that of the Phe–Arg complex exhibiting a neat minimum for vertical stacked and parallel displaced configurations. It is nonetheless worthwhile mentioning that the excess free energy of T-shaped structures (see Table 2) decreases in the order His–Arg > Tyr–Arg > Phe–Arg, i.e., as the partner residue of Arg becomes more hydrophobic. For all the three pairs, the solvent exposure of a stacked structure is always greater than that of a T-shaped structure. This suggests that the T-shaped configuration is stabilized basically by enthalpic forces, while stabilization is regulated by entropic factors in the case of stacking.

In the case of the Arg–Arg (RR) pair, we see a single neat and profound stacked minimum in the w_RR(r,) PMF, with an excess free energy well below −9 kJ mol⁻¹. The solvent exposure of the stacked dimer is systematic and significantly greater than that of the monomer, thus indicating the relevant role played by entropy in the stabilization of this complex. The stacked complex has a nearly perfect vertical geometry as can be seen from the w_RR(r,) PMF. Such results, although quite surprising given the net charge of the Arg residue, are fully consistent with theoretical investigations^14,15 claiming a strong entropic effect in the stabilization of the dimer. The large stabilization energy for the Arg–Arg stacked complex is also consistent with previous studies,⁶ where the Arg–Arg contact potential has been found comparable to that of the Phe–Arg and His–Arg pairs.

In contrast to aromatic–cation dimers, aromatic–anion ones have received in the past much less attention in chemical studies of non covalent supramolecular interactions, related to biological recognition.³ This may be due to the popular argument whereby a net negative charge is likely to interact unfavorably with an electron rich system. However, in recent years, aromatic–anion interactions have started to be extensively studied by means of accurate ab initio methods.^7,39,40 In particular, Kim et al.,³⁹ in their excellent study of complexes made of aromatic systems and halides or trigonal planar ions, thoroughly characterized the nature of –anion (and –cation) interactions at the MP2 level of theory with large basis sets. By applying a rigorous decomposition scheme, they found that the vertical stacked –anion configuration is as stable as –cation with a slightly larger contribution in the former arising from dispersion and polarization interactions.

In Fig. 7 and 8, we show the (r,) and (r,) PMFs for the aromatic–anion class. Salient data regarding the excess free energies and solvent exposure are reported in Table 3. As far as the PMFs are concerned, the differences with the aromatic–Arg PMFs are relevant, except for the Phe–Asp (FD) and Phe–Glu (FE) pairs. In the specific case of Phe–Asp and Phe–Glu pairs, we must stress once more the central role of stacked structures in agreement with ab initio data.³⁹

	Fig. 7 Contour plots of the (r,) (eqn (2.8)) and (r,) (eqn (2.9)) PMFs for Phe–Asp, Tyr–Asp and His–Asp pairs (left and right sides, respectively). On the right of each picture, the chromatic scale refers to the value of the PMF expressed in kJ mol−1. The bin size of the and angles is 10, while that of the distance r is 0.2 .

	Fig. 8 Contour plots of the (r,) (eqn (2.8)) and (r,) (eqn (2.9)) PMFs for Phe–Glu, Tyr–Glu and His–Glu pairs (left and right sides, respectively). On the right of each picture, the chromatic scale refers to the value of the PMF expressed in kJ mol−1. The bin size of the and angles is 10°, while that of the distance r is 0.2 .

RR					P(R)E
FD	−4.4	−2.0	0.64	0.69	0.36	0.80
YD	−4.1	−2.3	0.83	0.80	0.57	0.80
HD	−5.6	−0.7	0.82	0.72	0.68	0.80
FE	−5.5	−3.2	0.78	0.77	0.36	0.86
YE	−5.7	−2.8	0.88	0.85	0.57	0.86
HE	−5.8	−1.4	0.91	0.90	0.68	0.86

Stacked structures have an excess free energy between −5.8 and −4.1 kJ mol⁻¹, significantly higher than that observed for aromatic–Arg pairs. T-shaped configurations of aromatic–anion and aromatic–Arg pairs are instead comparable. The minimum in the (r,) PMF, seen for all pairs at < 10° and r < 4.5 , is again in all cases due to vertical or parallel displaced stacked structures (see (r,) PMF).

In a recent paper Jackson et al.⁷ have argued that the interaction of planar anions with aromatic rings occurs preferentially through coplanar structures (see Fig. 1). This speculation was supported by a statistical analysis on 946 proteins, where a clear preference for the = 0° angle between the normals of the planar anion and aromatic system was observed. Our data show that (i) stable coplanar geometries are not present for the Phe–anion pairs; (ii) weak minima corresponding to coplanar structures ( 0°, 9°) are visible in the (r,) PMF at distances around 4.5 for the His–anion pairs and at distances around 5.5 for Tyr–anion pairs. However, these latter coplanar structures present a PMF much lower than stacked structures. Such coplanar minima are actually a subensemble of a huge class of H-bonded structures clearly visible in the (r,) PMF of Tyr–anion and His–anion pairs as a wide belt with extending from 0° to 90°.

Finally, we note that aromatic–Glu stacked minima are essentially parallel displaced, while aromatic–Asp pairs have a clear preference for vertical stacking (see (r,)). In this regard, we find that Phe–Glu, Tyr–Glu parallel displaced interactions are generally stabilized by a hydrophobic interaction between the aliphatic side chain of Glu and the ring of the aromatic residues, as shown in Fig. 4 (structure C). These parallel displaced stacking configurations exhibit a very high solvent exposure probability (see Table 3).

In this article we have analyzed the thermodynamics of – interactions in proteins, extracting the statistical ensemble of proteins from the Protein Data Bank. In particular, we have evaluated the PMF in a 2-D space defined by the distance between the centroids of the residues and by the angle between the normals to the monomer planes ((r,) PMF). The stacked configurations have been further analyzed defining a PMF aimed at determining the parallel displacement of the residues ((r,) PMF). The results have been supplemented with data concerning the solvent exposure of the dimers and monomers. Overall, the results show the importance and ubiquity of the vertical stacked complexes in proteins, in many cases being the only relevant thermodynamic minimum found in the free energy surface. T-shaped structures are also favored in proteins although to a less extent with respect to stacked structures. In the ESI we show that the results obtained with our protein database are consistent with the corresponding data obtained with two of the most popular databases of non-homologous proteins.^33,34 In Fig. 9 we plot the excess free energy of stacked versus T-shaped complexes (data from Tables 1, 2 and 3). With the exception of Phe–Phe and Phe–Tyr pairs, stacked complexes are found to be consistently more stable than the T-shaped ones. For the series of aromatic–aromatic and aromatic–Arg pairs, the stability of the T-shaped complexes decreases with decreasing solvent exposure of the intervening monomers. Indeed, at an end, pairs involving the mainly buried Phe residue have comparable T-shaped and stacked stabilization energy, while, at the other end, dimers involving highly polar and exposed residues are characterized by an exceedingly stable stacked structure.

	Fig. 9 Stacked versus T-shaped excess free energy ( and , respectively) for the studied pairs (the monomers are denoted with a one-letter code). The dashed line is drawn as a guide for the eyes.

The aromatic–anion pairs, involving the highly solvent exposed Glu and Asp residues, also show very low T-shaped stabilization, but they also present less stable stacked complexes with respect to the other dimers. As a result, aromatic–anion pairs do not fit in the trend of the excess free energy correlation plot observed for aromatic–aromatic and aromatic–Arg pairs in Fig. 9. Since recent accurate quantum mechanical calculations³⁹ in vacuo have shown that stacked –anion systems are as stable as –cation ones, we must conclude that the entropic stabilization of the highly exposed aromatic–anion pairs is significantly less effective than in aromatic–aromatic, aromatic–Arg and Arg–Arg pairs. Very likely this is due to the smaller planar moiety of Glu and Asp with respect to the aromatic residues Phe, Tyr and His: the disparity in the physical dimensions of the monomers results in a non-optimal superposition of the planar systems, leaving the stacked complex more vulnerable to solvent intercalation.

With the exception of pairs involving the highly buried Phe residue for which energetically stable stacked structures are possible also in a non-polar environment,^10,41 in general we can say that T-shaped structures tend to occur more frequently in the protein interior with respect to the unbound dimers, while stacked arrangements show the opposite behavior. This is consistent with the results from computational studies of model systems^10,15 thus validating the methodology exposed therein for predicting the effect of solvent polarity onto the thermodynamic properties of – interactions. Conversely, it appears that the driving forces for the interactions of – systems in proteins mainly stem from basic thermodynamics rather than from specific structural or electronic effects.

In the light of recent evidence that the ion, and not the bulk water structure, is central to the Hofmeister series,²⁴ our results indeed suggest a possible role of the –ion stacking interactions in the denaturation of proteins. The capability of stacked interactions to locally increase the entropy at the protein surfaces may indeed explain the effectiveness of urea or of chaotropic planar ions in the Hofmeister series, such as SCN⁻ and guanidinium ions, in disrupting the tertiary structure.

Footnote

Electronic supplementary information (ESI) available: Excess free energies for residue pair arrangements from different databases of protein structures. See DOI: 10.1039/b718519g