CXC and CC Chemokines Form Mixed Heterodimers

CXC and CC chemokines are involved in numerous biological processes, and their function in situ may be significantly influenced by heterodimer formation, as was recently reported, for example, for CXC chemokines CXCL4/PF4 and CXCL8/IL8 that interact to form heterodimers that modulate chemotactic and cell proliferation activities. Here we used molecular dynamics simulations to determine relative association free energies (overall average and per residue) for homo- and heterodimer pairs of CXC (CXCL4/PF4, CXCL8/IL8, CXCL1/Gro-α, and CXCL7/NAP-2) and CC (CCL5/RANTES, CCL2/MCP-1, and CCL8/MCP-2) chemokines. Even though structural homology among monomer folds of all CXC and CC chemokines permits heterodimer assembly, our calculated association free energies depend upon the particular pair of chemokines in terms of the net electrostatic and nonelectrostatic forces involved, as well as (for CC/CXC mixed chemokines) the selection of dimer type (CC or CXC). These relative free energies indicate that association of some pairs of chemokines is more favorable than others. Our approach is validated by correlation of calculated and experimentally determined free energies. Results are discussed in terms of CXC and CC chemokine function and have significant biological implications.

kine CCL5/RANTES (19), as well as CC secondary lymphoid tissue chemokine (CCL21/SLC) and CXC B cell attracting chemokine-1 (CXCL13/BCA-1) (20). The functional result is that heterodimerization dramatically modulates the biological activities of these chemokines. For example, the presence of angiogenic CXCL8 in solution with anti-angiogenic CXCL4 increases the anti-proliferative activity of CXCL4 against endothelial cells (16). The co-presence of CXCL4 and CXCL8, in turn, attenuates the CXCL8-mediated rise in intracellular calcium in a myeloid progenitor cell line and enhances CXCL8-induced migration of bone marrow-derived pro-B-cells (Baf/3) (16,17).
This investigation, which was motivated by these studies (16,17), explores the energetic basis for heterodimerization of CXC and CC chemokines. Here we used the molecular mechanics and Poisson-Boltzmann surface area (MM-PBSA) approach to calculate binding free energies (21)(22)(23) and to predict which pairs of CXC and CC chemokines would likely form in solution. Free energy calculations from MD simulations have been used in various biological applications (24 -28), thus establishing their usefulness to better understand protein-protein associations (29). The MM-PBSA approach has emerged recently as a rapid computational approach that is broadly applicable to molecular systems that differ substantially in structure and/or are of comparable size, such as protein-protein complexes (21,23,28,29). Here we investigated formation of chemokine com-plexes by evaluating absolute binding free energies and comparing them to experimentally determined values.

EXPERIMENTAL PROCEDURES
Molecular Dynamics Simulations-Molecular dynamics (MD) simulations were performed for homodimers (CXCL8, CXCL4, CXCL1, CXCL7, CCL5, CCL2, and CCL8), heterodimers (CXCL4/CXCL8, PF4M2/CXCL8, CXCL4/ CXCL1, CXCL4/CXCL7, CXCL8/CXCL1, CXCL8/CXCL7, and CXCL1/CXCL7; two types of CCL2/CCL5, CCL2/CCL8, CCL2/CXCL8, CCL5/CXCL8, CCL2/CXCL4, and CCL5/ CXCL4; and two complexes with CCL5 mutants, CCL5(E26A)/CXCL4 and CCL5(44AANA47)/CXCL4) and homotetramers of CXCL4 and PF4M2. Starting protein structures for molecular dynamics simulations were built based on x-ray or NMR structures taken from the Protein Data Bank (30) without change. The corresponding Protein Data Bank entries were 1RHP, 1PFM, 1IL8, 1MSG, 1NAP, 1RTO, 1DOM, and 1ESR. The tetramer of CXCL4 is formed as a sandwich of dimers (AB and CD). Therefore, the homodimer of CXCL4 for the simulation was obtained by deleting one of the dimers (CD) from x-ray-derived tetramer. Heterodimers of CXCL4/CXCL8 (or other CXC chemokines) were formed by replacing one of the monomer subunits from an AB-type dimer of native CXCL4 with a monomer subunit from CXCL8 after FIGURE 1. Structures and amino acid sequences of CXC and CC chemokines. The folded structures of homodimers of CXC chemokine CXCL8 (11) and CC chemokine CCL2 (12) are illustrated to exemplify three-dimensional structures of all CXC and CC chemokine homo-oligomers. The structure of the tetramer of CXCL4 (10) is also shown, along with the CXCL4 AC-type dimers as explained in the text. The amino acid sequences of CXC (PF4M2, CXCL8, CXCL4, CXCL1, and CXCL7) and CC (CCL5 and CCL2) chemokines are shown at the bottom of the figure. superimposing the CXCL4 and CXCL8 homodimer. For CXCL4/CXCL8, a heterodimer of AC-type (adjacent monomers from dimers forming a sandwich) was also built by replacing C subunit of CXCL4 with a monomer subunit from CXCL8. For CCL2/CXCL8 and CCL5/CXCL8, two types of heterodimers were built, CXC-type (using CXCL8 homodimer as a template) and CC-type (using CC chemokine partner homodimers as a template). Fig. 1 shows the three-dimensional structures of CXCL4 tetramer, CXC (or AB)-type and AC-type dimers of CXCL4 derived from tetramer, and CC-type dimer exemplified by CCL2. All initial structures were built using the Insight II program (Biosym Technologies Inc., San Diego). Hydrogen atoms were added to the crystal structure using HBUILD module of the CHARMM program (31). The ionization state of the system was set at pH 5.0, the pH value used in our experimental studies (16). At this pH, the net charge on monomers was ϩ9 for CXCL4, ϩ6 for CXCL8, ϩ8 for CXCL1, ϩ7 for CXCL7, ϩ6 for CCL5, ϩ7 for CCL2, and ϩ5 for CCL8. Amino acid sequences of these chemokines are shown in Fig. 1. Although native CXCL4 and CXCL7 have 70 amino acid residues, coordinates for the first 6 residues in CXCL4 and first 4 residues in CXCL7 are not provided in the PDB files. Therefore, truncated CXCL7 has been used in this study, and absent N-terminal residues in CXCL4, which are important in forming CC-type heterodimers, have been added based on homology modeling using PF4M2.
MD simulations were performed using the c29b2 version of CHARMM (31). After an initial equilibration of 0.005 ns, we performed a 1-ns trajectory simulation for each hetero-and homodimer. The time step in the simulations was 1 fs. Coordinates were saved at 1-ps time intervals, resulting in a total of 1000 configurations for analysis. All complexes were simulated in boxes of 74 ϫ 63 ϫ 63 Å 3 or 70 ϫ 65 ϫ 65 Å 3 with explicit solvent molecules described by the TIP3P model, along with periodic boundary conditions. To make the total charge of the box 0, chloride ions were added to neutralize the system. Simulations were carried out using the CHARMM 22 all-hydrogen force field (32) with a dielectric constant of 1. 2000 steps of steepest descent minimization, followed by gradual heating to 300 K, and 5000 steps of system equilibration preceded each simulation run. The temperature during runs was maintained at 300 K. van der Waals (vdW) interactions were truncated at 13 Å using a shifted smoothing function, whereas electrostatic interactions were calculated using the Particle Mesh Ewald method (33). The SHAKE method was employed to constrain bond lengths involving hydrogen atoms (34).
Free Energy Calculations-The MM-PBSA method combines an explicit molecular mechanical (MM) model for the solute with a continuum Poisson-Boltzmann method for the solvation free energy. In this method, a molecular dynamics simulation is initially performed on a protein complex solvated explicitly using periodic boundary conditions. Then the solvent is removed, and the binding free energy, ⌬G, is calculated for each time increment according to Equation 1, where ⌬E MM is the change in molecular mechanical energy upon binding; ⌬G PBSA is the change in solvation energy, and T⌬S is the entropic contribution to binding. An estimate of the vibrational entropic contribution to the binding energy can be done by using normal mode (NM) analysis (35).
Strictly speaking, when oligomerization is considered, separate trajectories for each monomer and oligomer should be acquired and analyzed. However, we assumed that the conformation of the monomer subunit in both free and bound states is the same and estimated ⌬G from snapshots over the trajectory of the complex only. Therefore, energetic contributions from bonds, angles, and torsion angles are essentially not considered.
Chemokine dimerization or tetramerization was described by the following set of equilibria shown in Equation 2, M, D, and T indicate monomer, dimer (homo-or hetero-), and tetramer (CXCL4 only) states, respectively. The following thermodynamic cycle (Equation 3) was used to evaluate the free energy of dimer or tetramer formation in water, ⌬G w b , Here ⌬G g b is the free energy of dimer or tetramer formation in the gas phase. ⌬G sM1 , ⌬G sM2 , and ⌬G sC are the free energies of solvation for molecule 1 (M or D), molecule 2 (M or D), and for the complex (D or T), respectively. ⌬G w b was determined as the sum shown in Equation 4, All energy calculations were performed using the MM-PBSA approach (21)(22)(23) and the CHARMM program (31). Comparisons were then made between homo-oligomers and heterooligomers.
In the MM-PBSA method, ⌬G g b is calculated from molecular mechanics interaction energies shown in Equation 5, Here ⌬G g el is the electrostatic free energy, and ⌬G g vdW is the vdW interaction free energy between two molecules in the gas phase. Only complexes, and not monomers, were simulated; therefore, only nonbonded contributions to the free energy have been considered.
The solvation free energy was determined as the sum of electrostatic (⌬G g el ) and non-electrostatic (⌬G g non-el ), terms as indicated in Equation 6, The electrostatic contribution to the solvation free energy, ⌬G el , was calculated using the finite difference Poisson-Boltzman method with the PBEQ module in CHARMM (31). Both ⌬G s el and ⌬G g el were calculated simultaneously in two steps. For the initial Poisson-Boltzmann calculation, the grid size was set to 1.0 Å. For the second step, the grid size was decreased to 0.45 Å, and the box was filled 85%. The dielectric constant for the interior of protein is usually considered to be in the range from 2 to 4. For this study it was set to 2, and the dielectric constant of water was set to 80.
The non-electrostatic contribution to the solvation free energy includes vdW interactions between solute and solvent, ⌬G vdW , and is the free energy required to create the solute cavity in the solvent, ⌬G cav . The vdW and cavity terms were represented by a linear function of the total solvent-accessible surface area shown in Equation 7, where ␥ ϭ 0.00542 kcal/(mol⅐Å 2 ), and b ϭ 0.92 kcal/mol (36). For the solvent-accessible surface area calculation, the radius of the probe sphere was set to 1.4 Å.
Contributions to the binding free energy per residue were calculated as the difference in energy of a given residue in the monomer from that in the dimer (dimerization) or in the dimer and in the tetramer (tetramerization). Electrostatic and nonelectrostatic contributions were considered separately, using the approach described above.
Vibrational entropy was estimated from normal mode analysis using VIBRAN module of CHARMM program (31). The estimation was based on 40 snapshots taken from the last 800 ps of the trajectory (500 ps for tetramers) with an even interval.

RESULTS
CXC Heterodimers of CXCL4 and CXCL8-To validate our MD-based approach, we first investigated CXC chemokines CXCL4 and CXCL8, because heterodimerization between them has been demonstrated experimentally, and structural evidence as to how they associate as heterodimers is available (16). Three observations from these calculations on CXCL4 and CXCL8 support our approach as follows: 1) the time dependence of root mean square deviations (r.m.s.d.) for backbone heavy atoms; 2) time-averaged fluctuations per residue; and 3) expected free energy variations of residues at the homo-and heterodimer interfaces. Fig. 2A illustrates the time dependence of r.m.s.d. for backbone heavy atoms (averaged) of CXCL4 subunit in homodimer (black) and heterodimer with CXCL8 (green). Because the first six N-terminal residues (and last four C-terminal residues) in each dimer deviated substantially from their initial positions compared with other residues during MD simulations, these more flexible residues were not included in the r.m.s.d. values reported. Initial r.m.s.d. values for all structures ( Fig. 2A) leveled out at about 200 -300 ps (500 ps for CXCL4 tetramer, data not shown), and then fluctuated within less than 0.5 Å over the remainder of the simulation. The relative constancy of r.m.s.d.
values over the last 700 -800 ps of the trajectory suggests formation of stable homo-and heterodimers.
Structural stability is further supported by viewing time-averaged fluctuations per residue, as plotted in Fig. 2B for CXCL4 and CXCL8. Here the amplitude of fluctuations correlates well with elements in the folded structures of these CXC chemokines. For example, residues experiencing the largest fluctuations are located within the more flexible N and C termini and the loops, whereas residues within the less flexible, well structured regions (three-stranded ␤-sheet and C-terminal ␣-helix) show considerably smaller r.m.s.d. values (less than 1 Å) (37)(38)(39)(40)(41). This correlation supports the idea that our calculations reflect actual structural effects. Moreover, time-averaged fluctuations reflect essentially the same trends as observed experimentally (11,(37)(38)(39)(40)(41). Additionally, fluctuational amplitudes are essentially the same for homodimers and heterodimers, further supporting thermodynamic stability of the CXCL4/ CXCL8 heterodimer.
Calculated free energies are given in Table 1, and electrostatic and non-electrostatic contributions to these free energies are plotted per residue in Fig. 3, e.g. for CXCL4 and CXCL8 homodimers and heterodimers. For comparison, experimentally determined free energies are provided in Table 2. Although these calculated free energies are comparable with those determined using MM-PBSA approach with other proteins (28,29), it is readily apparent that the calculated energies far exceed those determined experimentally. Clearly this discrepancy lies in the limits of the calculated energies, as we will discuss later under "Discussion." Nevertheless, trends in these energies can be compared. In this regard, although hydrophobic interactions contribute the most and about equally to the energetics of CXCL4 and CXCL8 homodimer formation, electrostatic interactions contribute more to the energetics of homodimer formation for CXCL8 than for CXCL4. In fact, homodimerization of CXCL4 is electrostatically unfavorable, with a total electrostatic contribution to the binding free energy of ϩ5.54 kcal/mol for the homodimer and ϩ12.3 kcal/mol for the homotetramer (Table 1).
For CXCL4 homodimer formation, it was originally proposed that proximity of the two interfacial Glu-28 glutamates would be electrostatically unfavorable (10). Proximity of the two glutamates is illustrated in the CXCL4 CXC-type dimer structure inserts to Fig. 3. However, this is not the case, as positively charged residues are distributed around both glutamates. In fact, our results show that CXCL4 homodimerization is electrostatically unfavorable because of proximity of positively charged residues His-35, Lys-46, Lys-61, Lys-65, and Lys-31 which oppose each other at the homodimer interface. The side chain conformation of Lys-31, for example, changes significantly during the simulation because of electrostatic repulsion from Lys-46 on the opposing subunit. On the other hand, favorable electrostatic interactions occur between Glu-69 on one subunit and helix residues Lys-61 and Lys-65 on the other. Nevertheless, hydrophobic interactions among interfacial ␤-strand residues Ile-24, Leu-27, and Val-29 and C-terminal helix residues Tyr-60, Ile-64, Leu-67, and Leu-68, which become less accessible to water molecules when in the dimer state, drive CXCL4 homodimer formation.
In contrast to the situation with CXCL4, electrostatic interactions are favorable for CXCL8 homodimer formation (Fig. 3, top right panel), with a total electrostatic contribution of Ϫ16.5 kcal/mol (Table 1) compared with ϩ5.5 kcal/mol for the CXCL4 homodimer. The net charge on the CXCL8 monomer is ϩ6 at pH 5. However, the total number of charged amino acid residues is 26 (10 acidic and 16 basic), unlike CXCL4 which has 17 (4 acidic and 13 basic residues). The distribution of these residues in CXCL8 is such that all charged groups at the interface between subunits (Lys-23, Glu-24, Arg-26, and Glu-29) benefit energetically from dimerization. This is visualized in the CXCL8 structure inserts to Fig. 3. Even though two arginines (Arg-26) are brought close together in the dimer, the symmetric distribution of positively and negatively charged residues creates the favorable surroundings. In addition, C-terminal helix residue Arg-68 in the homodimer is positioned in a more electrostatically favorable environment, proximal to Glu-29 and Glu-E37 from two ␤-strands of the other subunit.
For CXCL4/CXCL8 heterodimer formation, there are significant changes in energetic contributions. In this case, CXCL4 gains electrostatically and loses non-electrostatically (in the tetramer), whereas CXCL8 loses electrostatically and gains nonelectrostatically (Table 1). These changes drive heterodimerization. The per residue free energies for CXCL4/CXCL8 heterodimer formation are plotted for individual CXCL4 and CXCL8 subunits in the lower part of Fig. 3, and three different views of the structure of this heterodimer are illustrated below these plots. As with homodimer formation, the most significant changes occur at the interface between subunits, because overall conformations of each chemokine monomer remain essentially the same. From the CXCL4 subunit side, Lys-46 and Lys-61 show the largest gain in electrostatic energy, because in the CXCL4 homodimer Lys-46 is situated near Lys-31, whereas in the heterodimer it is proximal to Glu-29 and Glu-37 of the CXCL8 subunit (Fig. 3, bottom, left-most structure) in the CXCL4/CXCL8 heterodimer. This now favorable electrostatic interaction promotes CXCL4/CXCL8 heterodimerization. A similar case presents itself for Lys-61, which in the CXCL4 homodimer is grouped with three other lysines (Lys-62, Lys-65, and Lys-66) and a glutamate (Glu-69) on the adjacent monomer. In the CXCL4/CXCL8 heterodimer, Lys-61 is now proximal to additional negatively charged group (Glu-63) from the  (16). From the CXCL8 side, C-terminal residues Glu-63 and Glu-70 of CXCL8 also favor heterodimer formation, because they are brought into proximity with four lysines and only one glutamate of the CXCL4 C-terminal helix (Lys-61, Lys-62, Lys-65, Lys-66, and Glu-69; Fig. 3, bottom, middle structure). On the other hand, CXCL8 residues Lys-23, Arg-26, and Glu-9 (Fig. 3, right-most structure) do not contribute favorably to heterodimer formation because of a loss in the balance of electrostatic charge in the CXCL8 homodimer. Thus, in the heterodimer, Lys-23 (CXCL8) is near Lys-31 (CXCL4), rather than Glu-29 in the CXCL8 homodimer, and even though Glu-29 still has a positively charged neighbor in the heterodimer (Lys-46 of CXCL4), Lys-46 is at a greater distance apart. In addition, Arg-68 of CXCL8 in the heterodimer loses its neighbor Glu-29 to the less favorable interaction with Lys-31 of CXCL4.
These environmental changes among residues of CXCL4 and CXCL8 subunits in the heterodimer result in quaternary structural shifts. In the heterodimer, CXCL4 and CXCL8 subunits are closer together and are oriented slightly differently from those in each homodimer. In particular, the C-terminal helix of the CXCL8 subunit is closer to the C-terminal helix from the CXCL4 subunit, because of more favorable non-electrostatic contributions from residues at the contact interface (in particular, Pro-53, Gln-59, Val-62, and Leu-66). These quaternary structural shifts, however, have little, if any, effect on the overall FIGURE 3. Binding free energies per residue for CXCL4 and CXCL8 homo-and heterodimers. Calculated free energies are plotted per residue for CXCL4 homodimers, CXCL8 homodimers, CXCL4 in CXCL4/CXCL8 heterodimers, and CXCL8 in CXCL4/CXCL8 heterodimers. Electrostatic (solid bars) and non-electrostatic (open bars) contributions to the free energy are shown. These contributions are calculated as the difference between contributions to the free energy of heterodimer and homodimer formation. Consequently, residues that have positive contributions are energetically unfavorable, whereas those residues that demonstrate a negative contribution are energetically favorable. Some residues most affected by subunit-subunit interactions in each respective dimer are labeled. Structures of CXCL4 and CXCL8 homo-and heterodimers are shown as insets below each respective panel, with two views each for CXCL4 and CXCL8 homodimers, and three views for the CXCL4/CXCL8 heterodimer. CXCL4 subunit is shown in blue and CXCL8 subunit is shown in red in each complex.
conformation of each monomer subunit. As mentioned above, NMR experimental investigations have provided evidence only for formation of the CXC-type CXCL4/CXCL8 heterodimer (16,17). In fact, if we calculate energies for another CXCL4 homodimer type that could form in solution, namely a ␤-sheet sandwich AC-type dimer (10) (see Fig. 1), time-averaged fluctuations for residues located at the inter-subunit interface are substantially larger than those for AB-type (CXC-type) homodimers, whereas fluctuations for residues located within the loops/turns and termini are essentially the same. Larger fluctuations at the interface indicate that AC-type dimers would be less stable than AB-type dimers, which is consistent with these NMR data, as well as with previous conclusions based on the x-ray structure of native CXCL4 (10).
CXC Heterodimers with CXCL1/Gro-␣ and CXCL7/NAP-2-Structural homology among all CXC chemokines suggests that other combinations will also form heterodimers. Having validated our approach with CXCL4 and CXCL8, we extended the study to two other CXC chemokines, CXCL1 and CXCL7. Binding free energies for CXCL1 and CXCL7 homodimer formation are given in Table 1, along with those for heterodimer formation with each other, as well as with CXCL4 and CXCL8. As for CXCL4 and CXCL8 discussed above, binding energies are high compared with those determined experimentally ( Table 2). Reasons for this lie in limitations to our approach and will be detailed under "Discussion." For now, one should consider that although calculated values are not absolute, comparisons or trends between/among calculated free energies do provide a realistic picture. Ideally, formation of heterodimers is deemed likely when the binding energy of heterodimerization is less than binding energies for each or at least one of corresponding homodimers, i.e. ⌬G w hetero b Յ ⌬G w homo b . In this regard, CXCL1 should readily form heterodimers with CXCL4 and CXCL8, whereas it has a lower probability to form heterodimers with CXCL7. The least likely heterodimers to form among this group of CXC chemokines are CXCL8/CXCL7 and CXCL4/CXCL7. Because CXCL7 and CXCL1 have 70 -80% sequence homology with CXCL4, similar per residue conclusions regarding homodimer formation of CXCL7 and CXCL1 can be made, as was discussed above for CXCL4. The primary difference is that electrostatic contributions to the binding free energy are more favorable for dimerization of CXCL7 and CXCL1 than for CXCL4. In CXCL7 like CXCL4, repulsive effects from negatively charged glutamates positioned at the inter-subunit interface are quelled by surrounding positively charged residues, whereas in CXCL1 that interface is composed entirely of nonpolar residues that allow closer packing of side chains.
Structural insight into why CXCL1 and CXCL7 form heterodimers with CXCL4 and CXCL8 comes from analysis of per residue contributions to the free energy, as illustrated in supplemental Fig. 1 for heterodimers CXCL4/CXCL1 (A), CXCL4/ CXCL7 (B), CXCL8/CXCL1 (C), and CXCL8/CXCL7 (D). When CXCL8 is paired with CXCL1 or CXCL7, charged residues Lys-23, Arg-26, Glu-24, and Glu-29 of CXCL8 (positioned at the inter-subunit interface) oppose heterodimer formation. This follows from the electrostatic nature of CXCL8 homodimerization, as discussed above. By placing these charged residues of CXCL8 in more unfavorable electrostatic surroundings, any CXC chemokine partner perturbs the sym- and Glu-69 from CXCL8 in the heterodimer, compensate partially for this. Alternatively, charged residues from both CXCL4 and CXCL1 (or CXCL7), although being less influential electrostatically, do contribute mostly positively to CXCL4/CXCL1 (or CXCL4/CXCL7) heterodimer formation. Although less important to CXCL4/CXCL1 heterodimer formation, nonelectrostatic contributions to the free energy from CXCL8 are significant, favoring CXCL8/CXCL1 heterodimerization. CC Heterodimers-As with CXC chemokines, structural homology among all CC-chemokine monomers makes CC heterodimer formation likely. To explore this, we performed the same free energy calculations on three CC chemokines, CCL5, CCL2, and CCL8. CCL2 and CCL8 form homodimers solely via contacts among N-terminal hydrophobic residues that form a short anti-parallel ␤-sheet (see Fig. 1). In CCL5, in addition to hydrophobic residues between Thr-7 and Tyr-14, charged residues Asp-6 and Arg-47 also contribute to homodimerization. Our calculations indicate that formation of a CCL5/CCL2 CCtype heterodimer is energetically favorable, with a ⌬G w ϭ Ϫ36.2 kcal/mol for the heterodimer, compared with Ϫ31.5 kcal/mol for CCL2 homodimer and Ϫ34.9 kcal/mol for CCL5 homodimer (Table 1). On the other hand, formation of a CXCtype CCL5/CCL2 heterodimer is, as expected, energetically much less favorable (Ϫ24.7 kcal/mol; see Table 1), due primarily to the presence of two proximal arginines (Arg-29 and Arg-30) in the middle of first ␤-strand that would be involved in intersubunit contacts in the CXC-type dimer, as well as in the vicinity of Lys-44, Lys-49, and Arg-24. A high content of positively charged residues on their C-terminal helices also makes formation of the CXC-type heterodimer of these two CC chemokines electrostatically unfavorable. Based on calculated free energies, CCL2 would also form a heterodimer with CCL8. ⌬G w for the heterodimer is Ϫ40.3 kcal/mol, compared with Ϫ31.5 kcal/mol for the CCL2 homodimer and Ϫ43.7 kcal/mol for the CCL8 homodimer. Formation of the CCL2/CCL8 heterodimer has been confirmed experimentally (42). However, although our calculated values indicate that either type of heterodimer, CC or CXC, could form, Crown et al. (43) concluded that formation of the CC-type heterodimer is more probable.
CXC and CC Mixed Chemokine Heterodimers-Because the protein backbone fold of any CC chemokine monomer is overall the same as that of any CXC chemokine monomer, we hypothesized that CC and CXC chemokines could form CXC/CC mixed heterodimers. If so, a related question would be which type of dimer, CXC-type or CC-type, would be more energetically favored. To address these questions, MD simulations and free energy calculations were performed on the following CXC/CC mixed heterodimers for both CXC-type and CC-type dimer motifs: CXCL4/CCL5, CXCL4/CCL2, CXCL8/ CCL5, and CXCL8/CCL2. Binding free energies (Table 1) are favorable for formation of CXC-CC mixed heterodimers of CXCL4/CCL5, CXCL4/ CCL2, and CXCL8/CCL2 but not of CXCL8/CCL5. Furthermore, for these pairs of CXC-CC mixed heterodimers, the CXC-type dimer is favored for CXCL4/CCL2 and CXCL8/ CCL2, whereas the CC-type dimer is favored for CXCL4/CCL5. Selection of a CXC-type or CC-type dimer depends on the particular residues at each respective inter-subunit interface.
Per residue energies are exemplified for CXC-type and CCtype heterodimers of CXCL8/CCL2 and CXCL8/CCL5 in supplemental Fig. 2. In the energetically favored CXC-type CXCL8/CCL2 heterodimer ( Fig. 2A), N-terminal residues of CCL2 (that normally interact at the CC-type dimer interface) lost energy, whereas CCL2 ␤-strand residues Ala-26 to Ser-33 at the CXC-type dimer interface gained energy. For the CXCL8 subunit, the N terminus was neutral, whereas interfacial ␤-strand residues lost energy (primarily electrostatic). In the CC-type CXCL8/CCL2 heterodimer (Fig. 2B), N-terminal residues of CXCL8 actually gained considerable energy, whereas ␤-strand residues Lys-23 to Glu-29 of CXCL8 lost even more energy (electrostatic and non-electrostatic) than when in the CXC-type heterodimer. Favorable per residue energies for the CCL2 subunit in the CXC-type heterodimer, were lost in the CC-type heterodimer. Overall, energetically the formation of CC-type heterodimers from both the side of CXCL8 and CCL2 was less favorable than the formation of CXC-type heterodimer. This was even obvious from the trajectories themselves ( Fig. 2A). In the CC-type CXCL8/CCL2 heterodimer, residues within the well structured ␤-sheet and C-terminal ␣-helix of CXCL8 fluctuated significantly, whereas in the CXCtype heterodimer, they were essentially the same as those in the CXCL8 homodimer. For the CCL2 subunit, per residue fluctuations did not differ significantly for either CC-type or CXCtype heterodimers, although they were generally somewhat lower in the CXC-type. Therefore, selection of heterodimer type was influenced more by CXCL8.
As mentioned above, the CXCL8/CCL5 heterodimers would likely not form as either CXC-or CC-type dimers. This is because the CCL5 homodimer is stabilized primarily by interactions among noncharged N-terminal residues that form a short intersubunit ␤-sheet and by electrostatic interactions between Asp-6 of one monomer and Arg-47 (third ␤-strand) of the other. These favorable interactions are disrupted when CXCL8 replaces one of the CCL5 subunits in either CXC-or CC-type heterodimer. Although CXCL8 Asp-4 might replace CCL5 Asp-6, the CXCL8 N terminus is shorter, leaving CXCL8 Asp-4 more distant from CCL5 Arg-47 in a CC-type heterodimer. In addition, proximity of CXCL8 Lys-3 and Arg-6 to CCL5 Arg-47 would electrostatically oppose heterodimerization. Furthermore, charged residues Lys-23, Glu-24, Arg-26, and Glu-29 at the CXCL8 homodimer interface lose their favorable electrostatic environment when in either CXC-or CC-type CXCL8/CCL5 heterodimer. In particular, close proximity of CCL5 Glu-26 and CXCL8 Glu-29 in a CXC-type heterodimer would be significantly electrostatically repulsive.
Although it is improbable that CXCL8 and CCL5 would form CXC-or CC-type heterodimers, we cannot rule out formation of an CXCL8/CCL5 heterodimer that has some novel dimeric structure. During the course of MD simulations for both CCand CXC-type CXCL8/CCL5 heterodimers, we noted that the mutual orientation of these two monomer subunits tended to deviate significantly, unlike another pair of chemokines we investigated. For the CC-type heterodimer, the subunits reoriented themselves in a way that Asp-4 of CXCL8 moved closer to Arg-47 of CCL5, and Asp-6 of CCL5 moved closer to Lys-11 of CXCL8. For the CXC-type heterodimer, Asp-4 of CXCL8 also moved toward Arg-47 of CCL5, whereas the two subunits moved further apart because of electrostatic repulsion between Glu-26 of CCL5 and Glu-29 of CXCL8. A much longer simulation and/or experimental evidence would be necessary for further insight into what heterodimer structure, if any, would form in this case.
Comparison with Experimentally Determined Free Energies-For a number of CXC and CC chemokines, experimentally determined equilibrium dissociation constants (and therefore free energies) are available in the literature. These K d values and related free energies are listed in Table 2, along with the solution conditions and experimental techniques used to determine them. For some of these chemokines, multiple K d values have been reported (all those known are listed in Table 2), and these can vary considerably depending upon such experimental variables as temperature, pH, and salt concentration. We chose to perform our calculations at what would approximate low salt and pH 5 with acidic residues being negatively charged and basic residues being positively charged (histidine residues, assumed to all have a pK a of 6, would be mostly positively charged). Many experimental K d values for chemokines ( Table  2) were determined either at pH 5 and low salt, or at about pH 7 and 200 -250 mM ionic strength.
In general, experimental K d values increase significantly with increasing pH and ionic strength. However, because free energies vary as the natural logarithm of K d , a 10-fold difference in K d is about a 2.3-fold difference in free energy. In Table 2, there are only a few instances where fair comparisons (i.e. using the same experimental technique, under otherwise the same solution conditions) can be made, and in these cases, either changing the pH from 5 to 7 or the ionic strength from about 20 to 200 mM results in about a 0.5-1.5 kcal/mol decrease in free energy. This is not a major factor when considering trends in our calculated free energies.
Experimental K d values listed in Table 2 are perhaps even more varied because of the technique used to determine them. The most reliable and consistent K d values seem to be those derived from NMR, calorimetry, and ultracentrifugation, whereas those derived from Biacore and chemical cross-linking appear to be overestimated in comparison.
Nevertheless, if we plot all experimentally determined free energies shown in Table 2 versus the respective calculated free energies (Table  1), we have the correlation plot shown in Fig. 4A that has a linear regression coefficient R ϭ 0.55. Although this R value seems low, it reflects the relatively large variation in experimentally determined free energies for reasons discussed above, i.e. pH, ionic strength, and experimental technique. For example, free energies reported for CXCL8 (IL8) dimer formation range from Ϫ6.0 to Ϫ9.4 kcal/mol (Table 2).
For a fairer and more reasonable comparison between calculated and experimental values, one should use experimental values determined at the same, or at least similar, solution conditions. In this regard, the largest subset of the experimentally determined free energies (Table 2) is that determined at low salt (Յ20 mM equivalents of NaCl) and pH 5, i.e. CXCL4 (PF4) and PF4-M2 (both dimer and tetramer), CXCL8 (IL8), CXCL1 (Gro-␣), CXCL7 (NAP-2), and mixed chemokines CXCL4/ CXCL8 and PF4M2/CXCL8. For IL8, however, the lowest salt concentration for which experimental data are available is 50 mM phosphate buffer; therefore, we adjusted the reported free energy from Ϫ9.4 to Ϫ10.4 kcal/mol because it has been reported that adding this much salt to solution will lower the free energy of CXC chemokine dimerization by about 1 kcal/ mol (13). Fig. 4B gives the resulting correlation plot, which has R ϭ 0.87, indicating a relatively good correlation between experiment and theory. However, even if we do not adjust this value, the R value remains relatively high at 0.83.
To expand upon this point, we calculated free energies as a function of ionic strength for CXCL1 (Gro-␣, Fig. 5) for which experimental free energies have been reported at several salt concentrations ( Table 2). As expected due to electrostatic screening effects, free energies decrease rapidly at lower salt concentrations and varied less at higher salt concentrations. In addition, only the Poisson-Boltzmann term (not unexpectedly) was affected by changes in ionic strength, whereas the van der Waals term was relatively independent of ionic strength. Experimental free energies for CXCL1 are plotted at the top of Fig. 5, and the inset shows the correlation between calculated and experimental values. Even though at lower ionic strength free energies are more variable, the correlation is relatively good with R ϭ 0.73.
To address the issue of pH variations, we also calculated free energies for CXCL1 at what would approximate the charge state of the protein at pH 7 by changing the ionization state of histidine residues (CXCL1 dimer has four histidines). In doing so, the calculated free energy decreased by 1.7 kcal/mol. Regardless of the accuracy of the calculation, the trend again parallels that observed experimentally (Table 2), where the free energy goes from Ϫ5.9 kcal/mol (pH 5) to Ϫ7.2 kcal/mol (pH 7). The same is true for CXCL7 (NAP-2), where the difference in calculated free energy is a net negative increase of about 1 kcal/ mol to be compared with an experimental free energy change from Ϫ5.4 kcal/mol at pH 5 to Ϫ5.7 kcal/mol at pH 7 (Table 2).
This having been said, CCL5 is somewhat exceptional in that it forms larger oligomers (tetramers and above) at higher pH and low salt (42). This is the reason why experimental free energies (as well as the NMR solution structure) were determined at around pH 4. However, in the absence of any high resolution structure of these higher oligomer states, we could not calculate free energies for any higher order hetero-oligomer state. The same can be said for CCL5 mutant E26A, which can form tetramers of unknown high resolution structure (42). The most that our calculations say about CCL5 is that this CC chemokine is favored or is not favored to form heterodimers with another CC or CXC chemokine.

DISCUSSION
Our previous work in this area demonstrated experimentally that CXCL4 and CXCL8 associate as heterodimers as the cause for functional modulation of their individual activities (16,17). The present results extend our knowledge on chemokine heterodimerization and indicate that heterodimer formation is thermodynamically favored among many other pairs of CXC chemokines (CXCL4, CXCL8, CXCL1, and CXCL7) as well as CC chemokines (CCL5, CCL2, and CCL8), and mixed CXC/CC chemokines. These findings are consistent with various reports in the literature as follows: 1) Biacore data from von Hundelshausen et al. (19) who showed that CXCL4 forms a heteromeric complex with CCL5; 2) immunoprecipitation results from Paoletti et al. (20) who proposed that CXCL13 and CCL21 might form a heteromeric complex; and 3) recent studies of Crown et al. (43) who studied heterodimerization between CCR2 receptor ligands and found that especially strong heterodimerization is observed between CCL2 and CCL8.
Conservation among the folded monomer structures of all of these CXC and CC chemokines generally allows heterodimerization, and alignment of specific pairs of amino acid residues at the dimer interface of individual monomers appears to primarily dictate the level of thermodynamic stability and selection of dimer type (CXC or CC). This is especially true for CXC-CC mixed heterodimers, where placement of specific amino acid residues (positively/negatively charged, polar, or hydrophobic) within ␤-strand 1 and/or within the N terminus determines selection of CXC-type or CC-type heterodimers. This study also provides insight into which specific residue interactions (i.e. electrostatic and non-electrostatic contributions to the free energy on a per residue basis) promote homodimer and heterodimer formation. This information cannot easily be determined experimentally. In general, heterodimerization is mediated primarily by non-electrostatic interactions, with the exception of CXCL8 and CCL5 where electrostatic forces contribute significantly to the binding free energy. Per residue free energies range up to about Ϫ3 kcal/mol, and although these values are relatively large on a per residue basis, they do correlate structurally, i.e. the largest values are associated with those residues directly involved at the inter-dimer interface. This, in turn, increases our level of confidence in the predictive values from these calculations.
Nevertheless, even though two of the highly energetically favorable heterodimers (CXCL4/CXCL8 and CXCL4/CCL5) have been observed experimentally (16,19) consistent with trends in our calculated binding free energies (Table 1), our calculated values are unrealistically large in magnitude compared with those determined experimentally ( Table 2). The main reason for these differences lies in the limitations of the computational approach. Use of a direct computational approach to correlate the experimental and calculated absolute binding free energies of macromolecular complexes remains quite challenging, even today. Although rigorous approaches, such as free energy perturbation and thermodynamic integration, have been shown to be applicable in the study of ligand binding (44,45), they are not very well suited for the study of large protein complex formation. Application of the end point free energy approach, such as the linear interaction model (46) within the implicit solvation model framework (21), has improved both efficiency and accuracy of these calculations. However, universal application of these computational methods remains limited by our understanding of the balance among various energy contributions to the total free energy of the system studied (21,47). These works have demonstrated that the balance among various energy contributions depends on the nature of the system, and the successful application of these methods requires regressionary scaling of these energy contributions to obtain accurate correlation with experiments. For protein association studies where the corresponding electrostatic and non-electrostatic energy terms can be significantly different (as is the case at hand), a regression type approach becomes difficult. In our present study of chemokine heterodimer formation, we evaluated absolute binding free energies based on the MM-PBSA method (21,48). Even though recent models have attempted to estimate various entropic contributions (49 -51), extension of these conformational sampling based methods toward macromolecular assembly becomes problematic. Our calculation of the absolute binding free energy was based on simulation of the protein complex alone. As such, the conformational entropic contribution of the monomeric protein was not taken fully into account, leading to higher-than-expected values.
Regardless of the limitation of the theoretical framework, calculated binding free energies are useful to predict formation of chemokine heterodimers when they are compared within this series of chemokine complexes. From the perspective of CXCL4, for example, heterodimers with CXCL8 are more energetically favored than CXCL4 homo-oligomers. On the other hand, the CXCL4/CXCL8 heterodimer is less energetically favorable than the CXCL8 homodimer. Nevertheless, CXCL4/ CXCL8 heterodimers will still form in equilibrium with homodimers, and given the fact that serum concentrations of CXCL4 are generally much higher than those of CXCL8, especially when platelets are activated as in the case of platelet aggregation and wound healing, mass action will shift the association equilibrium much in favor of the CXCL4/CXCL8 heterodimer. Therefore, biology should respond to the relative concentrations of respective pairs of CXC and CC (and possibly other) chemokine monomers.
This expanded view of CXC and CC chemokine heterodimerization presents a new way to think about chemokines in situ and about how they may function (52,53). A number of reports have found that various combinations of CXC and CC chemokines, although not all of them, can act synergistically (19,20,(52)(53)(54)(55). Overall, observed functional effects from mixed chemokines depend not only on the particular combination of chemokines but also on the biological assay and cell type used in that assay. In certain instances, we can correlate modulation of biological activities with the potential to form heterodimers, i.e. our relative free energies (Table 1). In this regard, we recently showed that heterodimerization of CXCL8 and CXCL4 increases the anti-proliferative effect from CXCL4 on endothelial cells in culture, as well as modulates the chemotactic effect from CXCL8 on Baf3 cells (16). The highly favorable free energy for formation of CXCL4/CXCL8 heterodimers (Table 1) is consistent with our previous experimental observation that CXCL4/CXCL8 heterodimers form in solution (16). Furthermore, combination of CXCL4 and CCL5 causes significant enhancement of monocyte arrest on the endothelium (19), and the CXCL4/CCL5 heterodimer is also highly thermodynamically favored. Furthermore, the chemotactic activity of CXCL8 on neutrophils is enhanced upon addition of chemokines CXCL4 and CCL2 (53), and the calculated association free energies for combination of this pair is also relatively highly favorable. On the other hand, addition of CXCL7 to CXCL8 exhibits no effect on CXCL8-mediated cell migration (29,34), and formation of CXCL7/CXCL8 heterodimers is calculated to be much less energetically favorable (Table 1).
Correlations with experiment (e.g. Fig. 4) increase the level of confidence in the predictive value of our calculated free energies, and even though further experimental evidence, both structural and functional, is required to validate all of our predictive findings and relate these to actual in vivo situations, our results do suggest significant biological implications. CXC and CC chemokines are involved in numerous biological processes where specific chemokines may be up-regulated or down-regulated depending on various biological stimuli. At the very least, the change in concentration of any one chemokine (due e.g. to inflammation) certainly would result in re-shuffling via mass action of most, if not all, chemokine homodimer and heterodimer populations that are co-localized at a particular site in situ, and this could have dramatic consequences biologically. Because of this, the potential for heterodimerization should be taken into account when presenting/discussing any experimental results involving CXC and CC chemokines.