Quantum Chemical Modeling of the Cardiolipin Headgroup - The

J. Phys. Chem. A 2010, 114, 4375–4387

4375

Quantum Chemical Modeling of the Cardiolipin Headgroup Martin Dahlberg,*,† Alberto Marini,‡ Benedetta Mennucci,‡ and Arnold Maliniak† DiVision of Physical Chemistry, Arrhenius Laboratory, Stockholm UniVersity, S-106 91 Stockholm, Sweden, and Department of Chemistry, UniVersity of Pisa, Via Risorgimento 35, 56126 Pisa, Italy ReceiVed: NoVember 19, 2009; ReVised Manuscript ReceiVed: February 5, 2010

Cardiolipin is a key lipid component in many biological membranes. Proton conduction and proton-lipid interactions on the membrane surface are thought to be central to mitochondrial energy production. However, details on the cardiolipin headgroup structure are lacking and the protonation state of this lipid at physiological pH is not fully established. Here we present ab initio DFT calculations of the cardiolipin (CL) headgroup and its 2′-deoxy derivative (dCL), with the aim of establishing a connection between structure and acid-base equilibrium in CL. Furthermore, we investigate the effects of solvation on the molecular conformations. In our model, both CL and dCL showed a significant gap between the two pKa values, with pKa2 above the physiological range, and intramolecular hydrogen bonds were found to play a central role in the conformations of both molecules. This behavior was also observed experimentally in CL. Structures derived from the DFT calculations were compared with those obtained experimentally, collected for CL in the Protein Data Bank, and conformations from previous as well as new molecular dynamics simulations of cardiolipin bilayers. Transition states for proton transfer in CL were investigated, and we estimate that protons can exchange between phosphate groups with an approximate 4-5 kcal/mol barrier. Computed NMR and IR spectral properties were found to be in reasonable agreement with experimental results available in the literature. Introduction Transmembrane proton gradients are central to the production of energy in the eukaryotic cell.1 A proton gradient can be created by active proton pumping from the inside of the matrix in the mitochondrion into the space between the inner and outer mitochondrial membranes. The spatial distribution of the protons has been the subject of disputeswith bulk or surface solvation being the main contending views.2,3 It has been proposed that anionic lipids at the membrane surface could work as proton sinks,4 and rapid lateral proton conduction at the lipid-water interface has been shown to be significantly faster than equilibration with bulk water.5–9 In parallel with these observations, one of the main lipid components of the inner mitochondrial membrane (IMM), cardiolipin (CL), has been found to have a peculiar pH-dissociation curve, indicating the ability to bind protons at a significantly higher pH than previously reported. In initial investigations, a double pKa value in the pH range 1-4 was reported for cardiolipin.10–12 A second pKa, between 7.5 and 9.5, was later proposed based on titration measurements.13 Recent Fourier transform infrared spectroscopy (FT-IR) data also support the presence of a second pKa at 7.9.14 Because cardiolipin was found at especially high concentrations in the IMM, and pKa2 in the physiological range, it was suggested that the lipid headgroup could serve as the IMM proton sink.15 In particular, intramolecular hydrogen bonding (in an intramolecular hydrogen bonded network, IHBN) was proposed to enable cardiolipin to bind protons at high pH.16 This notion was supported by simple models with the 2′hydroxyl group of the central glycerol group engaged in hydrogen bonding with the phosphate groups in the 1′ and 3′ * To whom correspondence should be addressed, [email protected]. † Division of Physical Chemistry, Arrhenius Laboratory, Stockholm University. ‡ Department of Chemistry, University of Pisa.

Figure 1. (left) Cardiolipin with chain connectivity. In black the subset of atoms used for the headgroup model. (right) Nomenclature for the headgroup model.

positions (Figure 1). Similar results were also found in lipids with analogous headgroups, such as the diphytanoyl phospholipids of extreme halophile bacteria.17–20 In the titration and proton conduction experiments mentioned above, a variant of the cardiolipin headgroup was also investigated: by changing the 2′-hydroxyl group for a hydrogen (deoxy-CL (dCL), 1′,3′propanediol) the possibility of intramolecular bonding was reduced, with significant effects on titration curves and conduction velocities. From a functional point of view, it is well established that cardiolipin-depleted mitochondria have decreased capacity for energy production.21,22 Although a changed cardiolipin density could contribute directly to a decrease in proton conduction, it is likely that other changes also contribute. Cardiolipin is associated with all major components of the oxidative phosphorylation chain and is involved in protein complex formation (for a review see, e.g., Schlame et al.23). There is also evidence that cardiolipin has effects on the mechanical properties of membranes.24

10.1021/jp9110019  2010 American Chemical Society Published on Web 02/26/2010

4376

J. Phys. Chem. A, Vol. 114, No. 12, 2010

The available physicochemical data were critically reviewed by Lewis and McElhaney.25 Although none of the available data show the existence of the intramolecular bond directly, a competing hypothesis must account for the differences in the pKa values of cardiolipin and its deoxy variant and further account for differences in proton conduction. The status of cardiolipin charge in the bilayer is still unclear. Experiments on mixed bilayers suggest that charge neutralization is achieved when monocationic surfactant to cardiolipin ratio is 2:1,26 where a ratio of 1:1 would be expected for a monoprotonated cardiolipin. To access the molecular details of acid-base reactions at surfaces, the experimental results require sophisticated analysis. The pH in the aqueous solution is not the same as at the surface, because changes in hydration and the existence of hydrogen bonds at the surface need to be considered.27–29 An additional complication is that the experimental conditions and techniques that have been used to determine the pKa of cardiolipin differ considerably. Cardiolipin, which has four acyl chains and a relatively small headgroup, has a propensity to go from the lamellar into the inverse hexagonal phase at low pH, high monovalent salt, or in the presence of divalent cations.12 Such phase transitions can be assumed to influence several of the parameters associated with the protonation state of the lipids and, therefore, increase the difficulty in finding the true cardiolipin protonation state. In this study we investigate the molecular structure of the cardiolipin headgroup and its 2′-deoxy derivative. We aim to give a description of the molecular requirements for an IHBN, which, according to the current hypothesis, is different in CL and dCL. Previous modeling work carried out in our group has involved coarse-grained properties of cardiolipin systems,30 and molecular dynamics at the atomistic scale,31 also recently done by others.32,33 Such methods are not suitable for investigating hydrogen bond formation or proton transfer in detail and have indeed not modeled the proposed IHBN. Instead, we turn to ab initio density functional theory (DFT) calculations to better describe the details of the cardiolipin headgroup. Recent developments in solvation continuum models34 have allowed accurate ab initio descriptions of acid-base equilibria.35–38 Others have carried out work on protons interacting with lipids with hybrid quantum mechanical/molecular dynamics models,39,40 but we are not aware of any similar work at the ab initio level on the cardiolipin headgroup. Early modeling investigations of cardiolipin were based on molecular mechanics force fields and did not focus on the headgroup structure, and both phosphate groups were assumed deprotonated.41,42 Recent theoretical work on phosphate acid-base equilibria have shown good agreement with experimentally determined pKa values,43–45 which is a requirement for an accurate description of acid-base equilibria in cardiolipin. We construct a model of the cardiolipin headgroup and its deoxy analogue where the two diacylglycerol groups have been replaced with methyl groups (see Figure 1) to make the computational problem more tractable. We use this truncated model as a first approximation, including only the most central structural features of the lipid headgroup thought to be involved in the IHBN. However, our model retains the characteristic phosphate diester groups found in CL and other lipids. Reliable environmental properties at the lipid-water interface are not available at the present, but we consider constraints of the phosphorus-phosphorus distance as an approximation of the effect of the lipid chains. We also describe the sequential protonation of the headgroup phosphates and determine the pKa

Dahlberg et al. SCHEME 1: Thermodynamic Cycles for the Free Energy Calculation of the Sequential Deprotonation of CL and dCL

of the process from the thermodynamic cycle shown in Scheme 1. We investigate transfer of protons between different sites in the cardiolipin headgroup, by means of identifying transition states, where the central glycerol hydroxyl group mediates proton exchange. The intramolecular barriers to proton transfer in the cardiolipin headgroup could limit proton exchange on the surface, and based on the transition state energies we show that such shuttling of protons is possible on the nanosecond time scale. Furthermore, the geometries of the CL and dCL models are compared with available structure data deposited in the Protein Data Bank, where lipids have been cocrystallized with proteins. Recent developments in ab initio methods and computer hardware have made accurate evaluation of NMR chemical shieldings and infrared spectra of many molecular systems possible.46–48 We compare available experimental NMR and FT-IR spectroscopic data for the cardiolipin headgroup with their computed equivalents. Finally, we revisit previous MD simulations of fully charged cardiolipin, and investigate the intramolecular and intermolecular effects of adding a proton to one of the phosphate groups in the headgroup. Methods DFT calculations using the B3LYP functional49,50 and the 6-31G(d) basis set have been performed to obtain the most populated conformational states of the polar headgroups of the CL and dCL systems. Gaussian0351 (build numbers E.01 and B.05) was used for all ab initio calculations, and initial conformations were built with GaussView 3.0. We use the notation CL0H/dCL0H, CL1H/dCL1H, and CL2H/dCL2H for the unprotonated, monoprotonated, and diprotonated species, respectively, and CL/dCL to address the molecules independent of protonation state, so that CL denotes molecules with a glycerol backbone, and dCL molecules with a 1,3-propanediol backbone. We further used letters (a-z) to label conformations according to the energy order in vacuum, with subscripts 0, 1, or 2 depending on the number of protons, and with a prime to denote the deoxy form. As an example, b1′ is the second (b) lowest energy state of the mono(1)-protonated deoxy(′) cardiolipin. PES and Conformation Generation. Conformations of CL and dCL models were generated from the variations of the dihedral angles φA, φB, ψA, ψB, χΑ, and χΒ (see Figure 1). For all the investigated systems, the main conformational search was performed in the B region employing a step by step relaxation scan procedure where the scanned dihedral angles were frozen and all the other structural parameters remained free. This procedure gave the potential energy surfaces PES(ψB,φB) and PES(χB,φB) in the range -180 to 180°, with a resolution of 20° (18 × 18 matrix) for the CL species. With the information gained from the CL systems, we used a lower resolution (60°) for the dCL system, which was sufficient to identify the significant minima. The final sets of minimum energy conformers were then obtained by using the structural information acquired in the region B, and the optimized dihedral angles for region A (φA, ψA, χA) to build conformations which were not tested in the dihedral scan. Finally, the geometry of each conformer was optimized without imposing any geo-

Modeling of the Cardiolipin Headgroup metrical constraint. Different proton positions were systematically generated for the protonated species. The solvation free energy for all conformers was calculated with the integral equation formalism-polarized continuum model method52 with the dielectric constant, ε, set to 4, 15, and 80, and the solvation cavity was created using default atom radii for the PCM method as implemented in Gaussian03. The intermediate dielectric constant, ε ) 15, was chosen to model the lipid-water interface. Additionally, geometry optimization with solvation (ε ) 80) was performed starting from gas phase optimized structures. pKa Calculations. The pKa of an acid/base reaction is directly proportional to the free energy difference between the reactants and products. This free energy difference was estimated from the thermodynamic cycles in Scheme 1. We calculated the pKa values according to

pKa ) 1/(2.303RT)[Gg°(A-) + Gsol(A-) - Gg°(HA) Gsol(HA) - 270.29] where the proton free energies in water and vacuum have been assumed to be -265.953 and -6.28 kcal/mol,54 respectively. Due to the change in standard states between 1 atm and 1 M, a constant factor of 1.89 kcal/mol has been added. Acids and bases for the respective reactions are denoted HA and A-, respectively. Free energies in vacuum, Gg°, were calculated with the zeropoint energies (ZPE) and thermal corrections at 298.15 K from vibrational analysis, and Gsol denotes the solvation free energy. The ensemble pKa can be calculated from the site-specific pKaij according to

pKa ) pKaij + log fi - log fj′ with fi and fj′ the populations of the unprotonated state i and protonated state j, respectively.55 The populations were calculated using the free energy of the corresponding state. Cardiolipin conformations were also extracted from Protein Data Bank entries. Residues with the cardiolipin identifier “CDL” or “CDN” (which both represent cardiolipin) were considered, and the coordinates of the heavy atoms in the headgroup (Figure 1) were used for interatomic distance calculations and geometry optimizations. Only “CDL” residues with all the heavy atoms of our model (see Figure 1) were used for optimizations. We found 73 such conformations. The considerably fewer “CDN” residues were not used for the optimizations because they had significant overlap in P-P distance with the “CDL” residues and represent relatively few protein structures. Hydrogens were added using the Reduce software.56 For CL1H and dCL1H all four protonation positions on the phosphate groups were considered. Transition State Calculations. For the CL1H system we located in vacuo transition states (TS) between reactants (R) and products (P) with the QST2 method.57,58 Structures P were constructed from R by moving hydrogen bonded protons along their respective bond vectors to new positions, changing the role of donor and acceptor with corresponding bond lengths. Before locating the TS, all structures in P were geometry optimized in vacuo. Vibrational frequencies at the TS were calculated, and one large negative mode was found for each TS structure. Additional stability for a subset of conformations was assessed at the 6-311G+(2d,p) level in the CL1H system. NMR Properties. The 13C and 31P chemical shielding (CS) tensors were calculated for all the minimum-energy conformers

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4377 found for CL and dCL model systems. The computations were carried out following the recommendations for chemical shielding calculations of organic molecules given by Cheeseman et al.59 For chemical shielding prediction we used the MPW1PW91 functional,60,61 which is particularly accurate in predicting the chemical shifts. Furthermore, because NMR chemical shielding calculations require accurate wave functions, the CS calculations were performed using the 6-311+G(d,p) basis set, which increases the accuracy and seems to be the best compromise between accuracy and CPU time. The gauge including atomic orbital (GIAO) method62 was used for CS calculations. All calculations were performed in vacuo and in solution as described above. We calculated 13C and 31P chemical shift parameters δ11, δ22, δ33, δiso, and ∆δ, where δ11, δ22, δ33, are the principal components of the chemical shift tensor (with δ11 g δ22 g δ33) and δiso, ∆δ are the isotropic chemical shift and the chemical shift anisotropy, defined as δiso ) 1/3(δ11 + δ22 + δ33) and ∆δ ) δ11 - 1/2(δ22 + δ33). Finally, the chemical shift tensors were evaluated relative to reference compounds, calculated using the same method and level of theory. Absolute chemical shieldings obtained by DFT are 181.7 ppm for 13C in tetramethylsilane (TMS) and 328.4 ppm for 31P in H3PO4.63,64 IR Properties. Infrared spectral frequencies and intensities were obtained from vibrational mode analysis of the optimized geometries in Gaussian03.51 The spectra were visualized by plotting each mode above a threshold intensity as a delta function. The threshold was chosen globally so that the hydroxyl and phosphate symmetric and asymmetric modes were included. The frequencies were not corrected to account for the harmonic approximation. MD Simulations. Molecular dynamics (MD) simulations of monoprotonated cardiolipin, i.e., corresponding to CL1H, were performed based on our previous work.31 A more detailed description of the technical details of these simulations is given in the Supporting Information. Results Description of Conformers. We use two levels of detail to describe the main findings in this study. First, we present a detailed geometric description of the conformations found to be important in vacuum and in solution. The energies of the optimized conformations in vacuum and with IEF-PCM (ε ) 80, opt) are presented in Figure 2, together with single-point calculations on vacuum optimized conformations with IEF-PCM (ε ) 4, 15, 80). In Figure 3 we show some of the important conformations of CL and dCL obtained in vacuum and in solution. Second, we look at the overall geometry of the molecules and the number of hydrogen bonds. Phosphate oxygen-hydroxyl oxygen (O-O) and phosphorus-phosphorus (P-P) distances are used as the two characteristic parameters for the molecular description at this level. All the studied molecules share a common backbone, which serves as a basis for comparison between different protonation states and between oxy and deoxy forms. We use the IUPAC conventions to describe the backbone conformations and focus on the φA/φB pair, which describes the central O-C-C-C-O fragment (see Figure 1). We also use the [φ,ψ] pair (see Figure 1) with bracket notation (for either side, A or B) to give a basis for hydrogen bonding between the central hydroxyl group and the phosphate groups. CL0H. Preliminary tests showed that hydrogen bonds between the central hydroxyl group and oxygens of one of the phosphate groups in CL0H produced a decrease in energy of

4378

J. Phys. Chem. A, Vol. 114, No. 12, 2010

Dahlberg et al.

Figure 2. Energies (kcal/mol) of the CL and dCL conformers, relative to lowest energy conformer in vacuum. From top to bottom: vacuum, ε ) 4, 15, 80 with vacuum geometry, and ε ) 80 with PCM optimization. Sequential naming (a-z) based on the vacuum energies.

Figure 3. Lowest energy conformers for protonation states 0H (charge ) -2), 1H (-1), and 2H (0) in CL and dCL: vacuum and optimized with PCM (ε ) 80). Several energy minimized PDB structures are included. Large images represent the lowest energy in the corresponding group. Energies are relative to the lowest state within the group, and selected P-P distances and backbone conformations are shown.

about 4-5 kcal/mol relative to nonbonded conformations. We then identified three minima with such hydrogen bonds (a0, b0, c0), which share the backbone torsion pair -sc/-sc, but each with a bond to a different oxygen atom of one of the phosphate groups. Conformations a0 and b0 exhibited, as expected, lower energy than c0 (bond to O-Me), with energy differences between c0 and the lowest energy conformer being 3.6 and 2.5 kcal/mol in vacuum and in solution, respectively. Conformers a0 and b0 had very similar energies, with a0 being the minimum by 0.44 kcal/mol in vacuum, and b0 the minimum by 0.56 kcal/ mol with PCM. Thus, the energies were insensitive to methyl group position, except when the hydroxyl bonded directly to the oxygen attached to the methyl group. All conformations showed a slight decrease of the distance between the charged phosphate groups when optimized with solvation, which is expected for a flexible molecule with negative charges when a dielectric medium is introduced. CL1H. A large variety of conformations were investigated for this species, which compared to CL0H has been protonated

and has a single negative charge. This is the protonation state that was hypothesized as the most stable at physiological pH. The first observation in vacuum was that in analogy with the CL0H energy map, formation of the first hydrogen bond resulted in a decrease of about 4-5 kcal/mol in the torsion scan. The second hydrogen bond lowered the energy by additional ∼6-7 kcal/mol relative to energy minima with one hydrogen bond. We found that conformers with the hypothesized intramolecular hydrogen bond network with bonds from phosphate to hydroxyl to phosphate indeed had low energies but that these were not the global minima. Instead, conformers where the charged phosphate group acted as an acceptor of hydrogen bonds from both the hydroxyl and the protonated phosphate group were the most stable. Interestingly, such conformations were of two kinds: with +sc/-sc (a1) and -sc/-sc (b1) torsion angles. The minimum energy conformation (a1) in vacuum was closely followed (0.35 kcal/mol) by the b1 conformer. In contrast, the main contributions from the hypothesized IHBN motif were 2.4 kcal/mol (c1) or higher above the minimum. The most significant

Modeling of the Cardiolipin Headgroup contributions from these structures were from conformers with a -sc/ap torsion pair in the central fragment. The potential energy surface (PES) was quite flat above the two (a1 and b1) minimal conformers. The global minimum in vacuum was also the global minimum when conformers were optimized with PCM. There were some slight rearrangements in the order of conformations with respect to energy, most notably in that the energy for the -sc/-sc conformer b1 was significantly higher in energy (3.60 kcal/mol). Similarly, the lowest energy IHBN conformer with the -sc/ -sc motif (j1) was 3.65 kcal/mol above the minimum. Solvation flattened the energy landscape close to the energy minimum, with the second and third conformers (c1 and k1) only 1.30 and 1.35 kcal/mol from the minimum. The two latter differ in that the methyl group was attached to different oxygens on the phosphate that donate a hydrogen to the central hydroxyl group, showing, as in CL0H, that the energy was insensitive to the methyl group position. The change in phosphate positions upon solvent optimization seen for CL0H was not observed for the CL1H conformers, which is reasonable given the reduced overall charge and two stabilizing hydrogen bonds for each conformer. CL2H. The neutral CL2H species presented a challenge because of the many possible hydrogen bonds with two protons and a hydroxyl group. The conformation found to be the energy minimum in vacuum (a2) had three hydrogen bonds. The +sc/ -sc torsion pair, with phosphates in close proximity introduced the possibility of forming a cyclic hydrogen bond network with each phosphate having a hydrogen bond to hydroxyl and to the other phosphate. The five lowest energy conformers of CL2H had cyclic hydrogen bond networks, with the energy rising steeply from 2.97 kcal/mol between a2 and b2 and 8.69 kcal/ mol between a2 and c2. The first noncyclic conformer was f2, with a -sc/ap torsion pair, 11.1 kcal/mol above a2. Clearly each hydrogen bond gave a strong reduction in the energy in vacuum, and the molecule backbone was flexible enough to allow cyclic networks without significant strain. As seen in Figure 2 there were major rearrangements in the PES as solvation was introduced, and all but one of the cyclic structures were unstable upon geometry optimization in solution. Interestingly, the cyclic a2 conformer was stable and still represented the minimum energy, but now with closer competition from three main structures e2, f2, and g2, all within 1 kcal/ mol. The changes here require further comment. In the case of e2, based on a -sc/ap backbone and the hydrogen donor on phosphate PA1 also being the hydrogen acceptor of phosphate PB2, the cyclic bond network collapsed. The minimum conformation was very close in energy to a2: 0.33 kcal/mol separating the two. Only one hydrogen bond remained: donor phosphate to hydroxyl. The contrast between a2 and e2 is rather striking, with three hydrogen bonds in the former and only one in the latter (see Figure 3). The f2 conformer constitutes an intermediate case with two hydrogen bonds. Clearly, the implicit solvation successfully competes with intramolecular hydrogen bonds for this neutral molecule, making the benefit of explicit intramolecular bonds smaller. The effect of solvent was also clear for conformer f2, with a -sc/ap backbone and the noncyclic IHBN involving the hydroxyl group, which changed very little from the geometry found in vacuum. The magnitude of the solvation energy was quite different for a2 and f2, as seen in Figure 2. Whereas a2 with the vacuum conformation exhibited sequential changes in energy of -7.94, -3.29, and -1.17 kcal/ mol at ε ) 4, 15, and 80, respectively, the corresponding energies for f2 were -14.28, -5.70, and -1.99 kcal/mol. The

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4379 TABLE 1: Torsion Angles [φ,ψ] Observed for the Backbone of CL with Hydroxyl Hydrogen Bonds in Vacuum Calculationsa CL0H

[-sc,-ac]

a

CL1H

CL2H

[ap,-sc] [ap,-ac] [-sc,-sc] [-sc,-ac] [-sc,+sc] [+sc,+sc] [+sc,+ac] [+ac,+ac]

[ap,-sc] [ap,-ac] [-sc,-sc] [-sc,-ac] [+sc,+sc] [+sc,+ac]

Items in bold also found with the PCM (solution) model.

PCM optimized geometry can be seen in Figure 2, and the final energy of f2 was only 0.44 kcal/mol higher than a2. Finally, conformer g2 had solvation energy similar to that of f2, here ending at 0.84 kcal/mol above a2. The difference between g2 and f2 was the connectivity of the methyl group at the phosphate donating the hydrogen in the bond to the central hydroxyl group, showing as in CL0H and CL1H that the energy was only weakly dependent on methyl connectivity. For the three lowest energy conformers the number of hydrogen bonds varied between 1 and 3, showing the strong effect of solvation. Importantly, seven of the eight cyclic conformations were unstable in solution, indicating that hydrogen bonds in such species are in strong competition with surrounding water. Noncyclic conformations would thus likely dominate in a medium with a dielectric constant close to that of water. With a partition of the torsion space into the six IUPAC angles, and noting that the sp (synperiplanar) conformation was not found in any of the minima for the backbone conformations, there are 25 possible pairs of [φ,ψ] angles. In Table 1 we have summarized the eight torsion pairs that were observed to be consistent with a hydrogen bond between the hydroxyl and a phosphate group. We note that hydrogen bonds were possible from either major backbone φ torsion angle [ap,x], [+sc,x], and [-sc,x], due to the flexibility afforded by the ψ torsion angle (here denoted x). For one conformer (w1), we observed a hydrogen bond with a +ac φ angle. The conformation was however only 2.6° from +sc, and does not constitute a different group of structures but should rather be seen as an extreme of the [+sc,+ac] pair. To gain more insight into the conformation requirements for hydrogen bonds, we also calculated the φ′ torsion angle (O-C-C-OH). With this definition of φ′, we found that hydrogen bonds were only possible when φ′ was +sc or -sc, but not in the ap position (i.e., φ elements in Table 1 transform to φ′ elements +sc or -sc). Because all ψ torsions were in the synclinal or anticlincal conformations, we note that only pairs of synclinal φ′ and synclinal or anticlinal ψ allowed hydrogen bonds with the central hydroxyl group. In the Supporting Information the torsion angles of previous MD simulations31 of cardiolipin bilayers (corresponding to CL0H) are presented. The only observed ψ in that study were in the ap basin, and the hydroxyl to phosphate group hydrogen bond was not observed, which is consistent with the results observed here. dCL0H. In analogy with CL0H, this system has two negative charges. With no possibility of forming hydrogen bonds, the conformation search was focused on the carbohydrate backbone. In vacuum, the fully extended, all trans (ap/ap) a′, 0 conformation had the lowest energy, followed by the +sc/ap conformer b0′ (+1.68 kcal/mol) and the -sc/ap conformer c0′ (+2.42 kcal/

4380

J. Phys. Chem. A, Vol. 114, No. 12, 2010

mol). Comparing the backbone conformation of CL0H, -sc/ -sc for all the studied minima, we see that the hydroxyl group in CL0H favors the ap to -sc conformation change on the side of the molecule not involved in hydrogen bonding. Upon geometry optimization in solution, the phosphate group distance decreased slightly, but less so than in CL0H. The order of the conformations was unchanged with respect to energy, but the energy differences were reduced significantly, with b0′ and c0′ only 0.02 and 0.28 kcal/mol above a0′, respectively. We interpret these changes in the energy landscape as an effective charge screening by the solvent model. Because the energies with PCM were close to indistinguishable, it is likely that the torsion energy of the backbone plays little role also in vacuum, the energy being dominated by phosphate-phosphate repulsion. dCL1H. The most stable dCL0H conformations were protonated as a first guess of the dCL1H conformations. We found that the resulting conformations were approximately 20 kcal/ mol higher in energy than the energy minimum, which was instead found to include phosphate-to-phosphate hydrogen bonding. The 10-membered cycles thus formed were flexible enough to give minimal strain energy. Four different cyclic conformers (a1′-d1′) were found within 4.5 kcal/mol. The solution energy changed the order within the set of cyclic conformers, with b1′ becoming the most stable, and reduced the spread in energies to 1.3 kcal/mol. The open conformers were 6-7 kcal/mol higher in energy in solution. dCL2H. The important energy minima for the neutral dCL2H species were cyclic in vacuum. The global minimum was found to have a +sc/-sc backbone torsion pair a′, 2 similar to the CL2H minimum a2. In vacuum, the closest conformer was 3.3 kcal/ mol above a′2 with open conformers approximately 10 kcal/mol above the minimum. The solution energy differences were reduced significantly, but unlike CL2H, cyclic arrangements were stable during geometry optimization. The open conformations f2′ and g2′ were 1.6 and 1.7 kcal/mol above the global minimum. Taken together these results suggest that although the phosphate-phosphate hydrogen bonds were destabilized relative to the open conformations, the strain the bonds imposed was less than that in CL2H, where the majority of cyclic conformers opened as solvent was introduced. Hydrogen Bonds and the Constraints of the Headgroup. Intramolecular hydrogen bonds constrain the possible donor and acceptor atom distances. We use the hydroxyl oxygen (O1) to phosphate oxygen (OA3, OA4, OA5, OB3, OB4, OB5) distance as a criterion for possible hydrogen bonds in our model of the CL headgroup. In Figure 4, the two shortest such oxygen-oxygen (O-O) distances are shown for different protonation states. Also shown are the corresponding distances for cardiolipin structures deposited in the Protein Data Bank. The PDB names for structures with potential hydrogen bonds are shown. The CL1H conformations showed the overall shortest first and second O-O distances, consistent with strong hydrogen bonding, and the CL2H and CL0H systems exhibited slightly longer O-O distances. Notably, structures 2H(HK/IT/JS: Reaction Centre from R. sphaeroides), 1SQP (cytochrome b-c1 complex), and 2C3E (Bovine Mitochondrial ADP-ATP Carrier) were localized close to a cluster of ab initio O-O distances composed of all CL0H structures and the a1 and b1 structures of CL1H. The 1SQP cardiolipin with short O-O distance is shown in Figure 4, with the minimum energy conformation CL0H a0 for comparison. The majority of the PDB structures were not compatible with intramolecular hydrogen bonds involving the central hydroxyl group, with distances for all potential acceptors and donors well

Dahlberg et al.

Figure 4. Two shortest opposite side (A-B in Figure 1) phosphate oxygen-hydroxyl oxygen distances in ab initio calculations: CL0H (dark green), CL1H (green), and CL2H (blue) represented with filled circles (vacuum) and crosses (PCM, ε ) 80). Arrowheads: lowest vacuum energy conformers. Images of a0 and 1SQP are shown for reference. Data from PDB (black filled circles). MD data31 on selected O-O distances (see text) shown as distribution (gray).

above 3 Å. This observation was similar to previous results derived from molecular dynamics (MD) simulations of pure cardiolipin bilayers,31 the hydroxyl to phosphate O-O distribution of which is also shown in Figure 4 as a distribution projected onto the horizontal axis. The peak at 6 Å from the simulations is an effect of including four phosphate oxygens (OA3, OA4, OB3, OB4) in the analysis, rather than just the two closest oxygens (the OA5 and OB5 oxygens are not shown but represent no deviation from the general behavior). The MD model did not show intramolecular hydrogen bonding, but its O-O distances agreed reasonably well with the majority of the structures from the PDB. There was also no evidence for phosphate-phosphate hydrogen bonds in the PDB set: the minimum O-O distance for such phosphates was 4.7 Å, well above the limit for hydrogen bonding. So far the discussion has been limited to the unconstrained molecule with an isotropic solvation model, and we have not considered the effect of chains and environment anisotropy. The solvation energies indicate that phosphate-phosphate bonds and ring structures such as a2 are in strong competition with surrounding water. It is clear that conformations with P-P distances at around 4-5 Å are cycles constrained by the hydrogen bonds. As a test of the importance of these bonds, and to estimate the effect of separating the phosphate groups from each other, we optimized CL0H and CL1H conformers from the PDB with the following constraint: the P-P distance was fixed, and the other degrees of freedom were optimized from the coordinates found in the PDB. We chose to constrain the phosphorus atoms rather than glycerol or acyl chain atoms because the former define the position of the charged or strongly polar phosphate groups that have well-defined contact areas with the proteins. Acyl chains in the known cardiolipin-protein binding cases, in contrast, have less-well-defined interactions with the protein surface.65 We stress that the optimization with P-P constraints does not, in general, preserve the geometry of

Modeling of the Cardiolipin Headgroup

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4381 TABLE 2: Calculated pKa for the CL and dCL Systemsa pKa1 CL, min/min CL, weighted n)1 n)2 n)3 dCL, min/min dCL, weighted n)1 n)2 n)3

1.80 2.41 1.80 1.83 2.42 3.42 5.78 4.10 3.88 3.91

(1.80) (1.81) (1.71) (3.66) (4.02) (4.96)

pKa2 11.38 11.57 11.05 11.12 11.37 10.97 11.28 11.72 12.23 11.31

(12.12) (11.02) (11.11) (11.23) (11.18) (11.31)

a Key: min/min, from minimum energies in solvation and vacuum; weighted, from all members of the ensemble; n, number of included members, sorted on PCM (vacuum) energies.

Figure 5. Energy as a function of P-P distance in the CL and dCL sets. Vacuum (black) and PCM (red, ε ) 80, optimized) for the a-z conformers, and vacuum (gray) and PCM (light red, ε ) 80) for the PDB conformers. Dotted lines show Coulomb energy of two unit charges as a function of distance relative to the lowest energy conformer.

the PDB conformations but rather that our method presents a set of realistic starting structures where exhaustive conformation sampling would be expensive. On the basis of the observed stability of intramolecular hydrogen bonds, and the lack of conformations where such bonds are possible in the PDB conformers, we expect the unoptimized PDB coordinates to be high-energy conformers. The energy dependence on the P-P distance for CL and dCL is shown in Figure 5. CL0H conformations from the PDB set reproduced our main conformer set within 0.5 kcal/mol, with very similar P-P distances. There was a clear inverse correlation between the P-P distance and the energy. This can be attributed to the repulsion between charged phosphate groups. The energy difference of two negative charges at a distance of 5.5 and 7 Å in vacuum with the Coulomb law gives 13 kcal/mol, which is in reasonable agreement with the ∼10 kcal/mol change seen for our conformer set. The balance between the two competing effects, phosphate repulsion and hydrogen bonding, was predictably shifted in solution, giving shorter P-P distances for the hydrogen bonded low-energy conformations in CL0H, and large overall energy reductions for dCL0H. The main difference between CL0H and dCL0H was in the positions of the PDB conformer set relative to the other optimized geometries. Here we saw an approximate 5 kcal/mol energy decrease due to hydrogen bonding in CL0H. One of the dCL0H PDB conformers with a +sc/+sc backbone (see Figure 3) had a significantly lower energy in solution than in our previously optimized set of solvated conformers and was included in the pKa calculation below. This conformation was not stable in vacuum and was therefore not found in the initial conformer search. For CL1H, the energy dependence on the P-P distance was much weaker than that in CL0H. Conformations with the hydroxyl bridge were found for the low-energy conformers with P-P distances of near 6 and 7 Å. As a test of the effect of the constraint on hydrogen bonding, the 6.9 Å conformer was optimized with PCM and found to be 4.3 kcal/mol above the global minimum. This we interpret as an effect of strain in the hydrogen bonds. The longest P-P distance in the PDB set was 7.35 Å, which was similar to the CL0H conformations with a single hydrogen bond between the hydroxyl group and the

charged phosphate group. This conformation was 10.2 kcal/ mol above the global CL1H minimum, (10.0 kcal/mol after PCM optimization). From these calculations we tentatively place the maximal P-P distance compatible with the IHBN between 7 and 7.35 Å. For dCL1H there was a positive correlation between distance and energy, due to charge neutralization and the possibility for hydrogen bonding with very short P-P distances. The solvated energies without such hydrogen bonds were 6 kcal/mol or higher above the minimum energy and quite uniform in the 5.5 - 7.8 Å P-P range covered by our test set. The 6 Å region, which had several low-energy CL1H conformers, is thus qualitatively different in dCL1H and CL1H. We did not investigate the CL2H/dCL2H case for the PDB conformers because we were mainly interested in the effect of the first protonation reaction. However, a few words on the energy-distance relations for CL2H and dCL2H are in place. The a2 conformation (see Figure 3) and similar cyclic structures with short P-P distances dominated in vacuum. Solvation effectively destroyed this advantage, and the energy profile was very similar to that of CL1H, with significant contributions at 4 and 6 Å. The situation in dCL2H was similar, but the energy of the extended conformations (∼6.5-7.8 Å) was 10 and 1.6-2 kcal/mol above minimum, as opposed to 20 and 6 kcal/mol in dCL1H in vacuum and PCM, respectively. Partial charge neutralization for both phosphate groups can (only) be attained in the cyclic dCL1H conformations, which may explain this effect. Upon adding an additional proton, charge neutralization becomes independent of the P-P distance, and the remaining differences can be seen as conformation and hydrogen bond energies. pKa Calculations. The calculated pKa values for the two systems CL and dCL are shown in Table 2. Due to the flexibility of the molecule, and different sites for the proton, we weighted the site-specific pKa values with their Boltzmann population fractions, calculated from the free energies of the corresponding conformers. In the pKa calculation, we included the dCL0H conformer from the PDB set that was the PCM energy minimum. In order to estimate the magnitude of the contribution from multiple states, we also calculated the pKa values by comparing the lowest energy conformation in vacuum and in solution, respectively (denoted min/min). This approach gave pKa1 ) 1.8 and pKa2 ) 11.4 for CL, pKa1 ) 3.4 and pKa2 ) 11.0 for dCL, which corresponded to a pKa gap of 9.6 for CL and 7.6 for dCL. The values for pKa2 were relatively insensitive to the choice of method, but pKa1 for dCL showed a shift of +2.4 pH units when the weighted method was used instead of the min/min method. Because the number of conformations for each protonation state was different, the comparison of free

4382

J. Phys. Chem. A, Vol. 114, No. 12, 2010

Figure 6. Reaction coordinates for proton transfer in the CL1H system. From top to bottom: vacuum, ε ) 4, 15, and 80. Left to right in order of increasing vacuum energy. Energies relative to the transition state are given in kcal/mol. Inset images show the reaction coordinate for a representative reaction coordinate.

energies should be corrected. We did this by including only the 1-3 lowest energy members of each protonation state in the weighted method, and we sorted the conformers by their energies with PCM and in vacuum, respectively. The results of this calculation are also shown in Table 2. The pKa2 values were similar for CL and dCL, but a slightly higher pKa1 was found for dCL. We will return to considerations of pKa in the discussion section. Transition States and Proton Conduction. The barriers to intramolecular proton transfer in CL1H are shown in Figure 6, with a typical example of the transition state structure. We considered proton transfer in the CL1H conformations assuming that this process could occur with only minor geometry changes along the reaction coordinate. Whereas the energy in the reaction coordinate in many cases showed significant asymmetry in the gas phase, the solvated structures exhibited very similar barriers to reaction with respect to the reactant and product conformations. The activation energies in vacuum for R to TS were 1.86-3.52-5.85 kcal/mol (min-mean-max), and for P to TS 0.93-3.23-5.59 kcal/mol. The corresponding energies obtained with ε ) 80 were 3.43-4.14-4.86 and 3.04-4.27-4.69 kcal/ mol, showing a mean activation energy increase. Intermediate solvent strengths, ε ) 4 and 15, showed corresponding increases, 3.90 and 4.13 kcal/mol, respectively. As a test of stability, the activation energies for g1-TS-e1, chosen as a representative of low energy conformations with a symmetric reaction barrier, were also calculated with the larger basis set 6-311G+(2d,p). The resulting mean barrier height for the vacuum and ε ) 4, 15, and 80 calculations decreased by 0.91, 0.62, 0.54, and 0.52 kcal/mol, respectively. Reactant (g1), product (e1), and TS were also optimized with solvation and showed only minor changes in geometry and slight decreases in activation energy from the reactant side (3.85 vs 3.43 kcal/mol) and the product side (4.38 vs 3.85 kcal/mol). The activation energies found here were comparable to proton transfer barriers in clusters of phosphoric acid.45 Assuming a reaction rate k ) kBT/h e-Ea/kBT, with kBT the thermal energy at 298.15 K, h is Planck’s constant, Ea the activation energy, i.e., neglecting changes in entropy between the transition state and the reactant, we calculated characteristic reaction times (1/k) from the activation energies in vacuum and

Dahlberg et al. water to be 48 ps (Ea ) 3.37 kcal/mol) and 197 ps (Ea ) 4.21 kcal/mol), respectively. Averaging reaction times over conformations instead, we found that both vacuum and solvent reaction times were in the 0.2-0.4 ns range, with significantly larger variance in vacuum. Computed NMR and IR Spectra. The 13C isotropic shifts, shown for the minimum energy conformations in Figure 7, for both CL and dCL were comparable with those commonly measured in organic molecules having phosphate groups,66 and in good agreement with those measured in CL species. Specifically, we found that the isotropic chemical shift for the C1 central glycerol (CHO) was in the range of 69-72 ppm, while the other central glycerol carbons CA2 and CB2 were shifted to lower frequencies, 57-69 ppm, in good agreement with those measured by Kraffe et al.67 (70.5 and 66.6 ppm for C1 and CA2/ CB2, respectively). In a glycosylated CL species investigated by Rˇezanka et al., the resonances found for C1 and for CA2/ CB2 were 81.6 and 64.2/63.1, again in reasonable agreement with our results. Additionally, we noticed a significant shift (from 67-72 to 28-33 ppm) in the C1 resonance going from CL to dCL, attributable to the absence of the central OH group dCL.68 Variations in the isotropic shift on the order of 5 ppm were found for smaller basis sets and the B3LYP functional. The calculated isotropic shifts for 31P were in the 15-40 ppm range for both CL and dCL systems, which is significantly larger than shifts recorded for the full cardiolipin in chloroformmethanol (two resonances at -0.8 and -1.0 ppm).69 Phospholipid 31P chemical shifts are typically within a few ppm around zero,70 and thus the calculated properties are not very accurate. We found a large basis set and method dependency for the 31P isotropic shifts, unlike the quite precise 13C shifts, and thus the predictive value of the calculations for 31P is very limited. The 31 P properties for the MPW1PW91/6-311+G(d,p) combination are however shown in Figure 7 for completeness. In contrast to 13C shifts, 31P nuclei were sensitive to the solvation state, with anisotropic shifts in particular showing variations of more than 30 ppm, but generally in the range commonly found for phospholipids (180-190 ppm),70 and closer to that range with increasing protonation. Given the poor accuracy of the 31P isotropic shifts, the agreement can however be considered fortuitous, and we are hesitant to draw any conclusions as to the protonation state in the experiments from these data. Methyl group carbons, CA3/CB3, showed the least sensitivity to solvent effects, with variations in the 13C δiso within 4 ppm over the range of dielectric constants. In general, the C1, CA2, and CB2 carbon chemical shifts were more sensitive to conformation changes in the backbone of CL and dCL, than to the solvation state. The large changes in δiso and ∆δ for a given conformation occurred when passing from vacuum to PCM (ε ) 4), and relatively minor changes occurred between ε ) 4, 15, and 80. The most significant systematic change in the isotropic 13C shift with protonation state was found for the CL and dCL terminal methyl groups, which increased their chemical shift with phosphate protonation. Other carbons and the phosphate atoms showed quite diverse chemical shift patterns, and their intraspecies variance was comparable to the changes with protonation state. This behavior is reasonable considering the local structure being highly dependent on the existence of hydrogen bonds. Calculated infrared (IR) spectra for all solvated conformations are shown in Figure 8. The two main effects in the calculated IR frequencies originate from solvation and the protonation state.

Modeling of the Cardiolipin Headgroup

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4383

Figure 7. 13C and 31P NMR chemical shifts calculated for CL and dCL. Each protonation state contains three low-energy structures, i.e., PA1/PB2, CA2/CB2, and CA3/C3 correspond to six structures (two atoms, three conformations) each.

Figure 8. Calculated IR spectra from the vibrational analysis. Each mode represented by a line, conformations a-z from bottom to top for each protonation state. A global cutoff in the IR intensity was used to filter out weak signals. The dotted line and arrow serve as a guide for the shift in the asymmetric (∼1250 cm-1) mode between the CL1H and CL2H states.

From vacuum to PCM the main effect in the hydroxyl region was to lower the frequencies upon solvation. However, for a majority of the CL1H conformers the opposite effect was found, the increase in frequency likely an effect of reducing the strength of the intramolecular hydrogen bond network. The phosphate symmetric and asymmetric vibrations were shifted by ∼-25 cm-1 and ∼-50 cm-1, respectively, by introducing solvation (data not shown). The major difference between protonation states was observed for the hydroxyl frequencies. Monoprotonated species (CL1H and dCL1H) had the lowest frequencies, with vibrations in the

2200-3000 cm-1 region, mainly for the IHBN conformations, and some vibrations in the 3000-3500 cm-1 range, mainly conformations with phosphate-phosphate hydrogen bonds. Doubly protonated CL and dCL, and CL0H, had hydroxyl frequencies in the 3000-3700 cm-1 range. There were slight shifts in the symmetric (at approximately 1050 cm-1) and asymmetric (1250 cm-1) phosphate vibrations going from CL1H to CL2H, with the asymmetric modes shifting up more systematically on the order of ∼10 cm-1 (shown in Figure 8, with large intraspecies variances). Frequencies of the symmetric phosphate vibrations in dCL1H and dCL2H were generally lower than in CL1H and CL2H. It is worth mentioning that vibrational modes in the 1050 cm-1 region also included vibrations such as the COPOC modes, and pure symmetric phosphate modes for all conformers were not found. Frequencies agreed reasonably well with experimental IR spectra for cardiolipin and 2′-deoxycardiolipin, where the symmetric phosphate stretch was located at ∼1090 cm-1, and the asymmetric phosphate stretch at 1220-1230 cm-1.16 The differences between cardiolipin and deoxycardiolipin obtained experimentally were smaller than the variance in peak positions for our set of conformations.14 In the aqueous phase it is impossible to discriminate between frequency shifts originating from intramolecular bonds and those related to hydrogen bonds with water. For the dry films there was a slight decrease in the symmetric and asymmetric phosphate frequencies, interpreted to be an effect of intramolecular bonding in the CL headgroup. Recent FT-IR data showed a shift in the asymmetric mode from 1227 cm-1 at pH 2 to 1238 cm-1 at pH 5, and 1230 cm-1 for pH 8, which was interpreted as corresponding to the protonation states CL2H, CL1H, and CL0H, respectively.14 Our data show the opposite trend for the CL1H-CL2H transition, but with similar magnitude, and no significant change in the CL0H-CL1H transition. We note that the frequencies were not rescaled to account for anharmonicity, which for our choice of method would tend to decrease our calculated frequencies by a couple of percent.48 MD Simulations of CL1H Cardiolipin. To better understand the impact of protonating cardiolipins in a membrane environment, we performed MD simulations of bilayers consisting of CL1H protonated at either free oxygen atoms at one of the phosphate groups. The results from these simulationssdistributions of intramolecular and intermolecular distances and the overall patterns of hydrogen bondingsare collected in the Supporting Information and are discussed below.

4384

J. Phys. Chem. A, Vol. 114, No. 12, 2010

Discussion Although the experimental titration results suggest that intramolecular hydrogen bonds are important for the acid-base reactions of the cardiolipin headgroup, there is little direct evidence of the actual bonds. A detailed molecular understanding is restricted by experimental difficulties, which for instance limit the usefulness of IR spectroscopy due to headgroup exposure to water, and further insights into the properties of the cardiolipin headgroup can be aided by computational models. We introduced simplifications in the cardiolipin models by replacing the diacylglycerol fragments with methyl groups. Although this approximation ignores the effects of, e.g., carbonyl dipoles, it still represents all the atoms necessary for the proposed IHBN. Because the whole cardiolipin molecule was not modeled, some comparisons with similar compounds are in place. The structurally related compounds glycerol-1,3bisphosphate (CL analogue) and 1,3-bisphospho-1,3-propanediol (dCL analogue), were previously shown to have pKa1/pKa2 at 2.25/2.35, and pKa3/pKa4 at 6.3/7.1, respectively,13 i.e., with no sign of anomalous titration behavior, and behaving much like phosphoric acid in aqueous solution (pKa1 ) 2.2, pKa2 ) 7.2).71 A slight decrease in the stability of the protonated species is expected for methyl-substituted phosphates, an example being dimethyl phosphate with an experimental pKa of 1.29.72 Lipid analogues with phosphatidylglycerolphosphate (PGP) headgroups have also been investigated experimentally. TomoaiaCotisel et al. titrated the methylated PGP-Me (pK1 ) 3.00, pK2 ) 3.61), PGP (pK1, pK2 as PGP-Me, pK3 ) 11.12), and dPGP (pK1 ) 2.72, pK2 ) 4.09, pK3 ) 8.43) lipids, which showed that the gap between the second and third deprotonation step was considerably widened due to the presence of the hydroxyl group. This was taken as support for the notion of the intramolecular hydrogen bond in the hydroxylated species.20 In contrast to our calculations, the first and second deprotonation steps for all these model compounds are at low pH. The only experimental observation consistent with our pKa results for CL is the titrations of the full lipids of Kates et al.,13 who did not, however, find a similar pKa gap in the full dCL lipid. Further experimental work, e.g., titrations, on the analogues we use could show if the IHBN can indeed be present in solution or if the constraints are necessary. In the present study we found that the CL and dCL species both had a large gap between pKa1 and pKa2, which can be attributed to intramolecular hydrogen bonds. Notably, intramolecular bonds were important not only in CL, which was suggested previously, but also in dCL. For CL we found small energy differences between intramolecular bonds directly between phosphate groups and bonds involving the IHBN with the central hydroxyl group. Because dCL cannot form the hydroxyl bonds, the high pKa2 can be attributed only to the phosphate-phosphate hydrogen bond. With reduction of the charge on one or both phosphate groups, the energy surface flattened considerably and allowed more compact geometries with several hydrogen bonds. It is notable, however, that the gains in hydrogen bond energy decreased upon protonation, which can be seen in the case of CL2H. There, the number of potential hydrogen bonds was three (counting only single donors and acceptors), but most conformations with three hydrogen bonds were unstable when solvated. Equivalently, the added benefit of a hydrogen bond was reduced in the neutral state. We did not try explicit representations of water, which might influence the strength of the hydrogen bonds due to specific interactions. In a previous work on myo-inositol-2-monophosphate, intramolecular hydrogen bonds between phosphate and

Dahlberg et al. hydroxyl groups were not disrupted by adding explicit water molecules, which indicates that an implicit solvation model is reasonable for this type of system.43 The most stable ab initio conformations of CL1H and CL2H have no counterpart in the lipids found in the PDB, as seen, e.g., by the P-P distance. A few PDB lipids were close to the geometry of CL0H, with hydrogen bonds between the central hydroxyl group and one phosphate group. This implies that proton trapping by the IHBN mechanism suggested for membrane cardiolipins is not valid for the available protein bound cardiolipins. It is important to note that the lipid-protein interaction potentially changes the pKa of the lipid, and we should not extrapolate to the free lipid situation from these considerations of geometry and energy. On the basis of the representatives of cardiolipin-protein interactions deposited in the PDB set, we are led to conclude that interactions are primarily between protein and the doubly anionic form of cardiolipin at physiological pH. We found activation energies for the proton transfer reaction to be on the order of 4 kcal/mol. It is known that proton barriers are underestimated with DFT methods as compared to postHartree-Fock methods.73,74 For malonaldehyde the difference between CCSD(T) and B3LYP, the latter of which was used in the present study, was approximately 1 kcal/mol. Compensating for the effect of using DFT, an increase in the activation energy of 1 kcal/mol corresponds to a decrease in the reaction rate with a factor ∼5 (5.41). Thus we expect reactions on the time scale of a few nanoseconds. This can be compared to the typical time scale of proton release from surface to bulk, which is on the hundreds of microseconds time scale8 allowing many intramolecular proton transfer events along the coordinate that we described. We did not investigate proton transfer in the noncyclic CL2H structures, because it would lead to charge separation. The cyclic conformation a2 potentially has a similar proton transfer circuit, but this was not pursued further as it would not lead to a net transfer of charge. The proton transfer barriers found here were similar in magnitude to the strength of the hydrogen bonds for CL0H and CL1H, as judged from the initial PES scans and the constrained optimizations. We therefore predict headgroup dynamics (breaking and forming of hydrogen bonds) to be on roughly the same time scale as proton transfer. The MD simulations of CL1H cardiolipin indicate that although hydrogen bond formation is sensitive to the LJ parameters of hydrogen, presence or absence of hydrogen bonding tends not to affect the intramolecular conformations and the area per lipid strongly, at least on the 5-30 ns time scale. Although phosphate-phosphate hydrogen bonds were observed in the MD simulations, these were between adjacent lipids, and the intramolecular P-P distance distribution shows that intramolecular phosphate-protonated phosphate bonding did not occur, indicating that such conformations are not compatible with cardiolipin in bilayers. We did not observe the IHBN in the MD simulations, but the significant number of intermolecular hydrogen bonds in the no-LJ simulations indicates that the lack of intramolecular bonds might be due to the LJ parameters on the central hydroxyl group hydrogen in the original simulations. However, setting the LJ parameters to zero, which is the standard treatment for hydrogens, creates excessively stable, unphysical, bonds between the hydroxyl group and OA2/OB2 phosphate ester oxygens. This shortcoming of the MD force field highlights the sensitive balance of forces in the headgroup, and reinforces the need for a detailed quantum chemical description of the cardiolipin headgroup.

Modeling of the Cardiolipin Headgroup In principle, the free energy differences for the different protonation states should be corrected for symmetry and invariance of proton position. In the pKa calculations we have only corrected for the number of low-energy conformations, which will in theory give a bias to CL0H and dCL0H, relative to the protonated species. Because many different conformations are accessible, and we are not able to exhaustively search the entire conformational space, free energy estimates given here are associated with some uncertainty. The different methods tested for pKa calculations show this, with variations of approximately (1 pH unit (corresponding to 1.4 kcal/mol). The observed differences between pKa1 and pKa2 for both the CL and dCL systems are much larger (7-9 pH units) than the variance of the different methods, which is an indication of the strength of the hydrogen bonds formed. Importantly, the magnitude of the pK gap is influenced by the size of the solvent cavity (see, e.g., the work of Verdolino et al.38), which was not varied in the present study. Thus, the significance of the 7-9 pH unit gap that we found is in that it was similar in CL and dCL. The dielectric constants used in the solvation model were chosen to span the different environments close to the membrane interface, with ε approximately 4 in the hydrocarbon core, 80 in the aqueous phase, and intermediate values, here represented by ε ) 15, at the interface. Importantly, the interfacial electrostatic interactions are highly anisotropic,75 which is not included in the present model. We also conclude, however, that optimization with the solvation model is crucial for accurate representations of the energy landscape. This was most clearly seen in the CL1H system, where the majority of the cyclic structures were stable in vacuum, but highly unstable in solution. Solvated conformers also gave spectroscopic parameters in reasonable agreement with experimental data. The main physical insights presented here are that multiple glycerol backbone geometries are compatible with the IHBN and give similar transition state energies, 4-5 kcal/mol, for the proton transfer. These energies correspond to reaction rates on the nanosecond time scale and were not sensitive to decreasing the dielectric constant to values relevant for a lipid-water interface (ε ) 15). We also present limits on the geometry which still allow the IHBN to persist: P-P distances up to approximately 7 Å with energy penalties of approximately 4 kcal/ mol from the minimum and 3 kcal/mol from the minimum IHBN conformation. This is relevant for considering the headgroup constrained by, e.g., acyl chains in the bilayer or cardiolipin binding to a protein surface, and will assist in future MD models of cardiolipin, which to date have not modeled the IHBN. In the case of the cardiolipin model membrane, intermolecular hydrogen bonds are relatively scarce because the 2′-hydroxyl group is the only hydrogen bond donor. Competing hydrogen bonds from water, which is only included in the continuum approximation, and electrostatic interactions from ions and other lipids can be considered in the perspective of the approximated hydrogen bond strengths in the IHBN found here. We have shown that intramolecular hydrogen bonds in CL and dCL are consistent with a significant gap in pKa values, and that two types of cyclic conformations are represented. The cyclic structures that involved hydrogen bonds between the phosphate groups have not been reported before. These structures were sensitive to solvation effects, but some conformations persisted also with a dielectric constant of 80. In CL these conformations were within 1.3 kcal/mol from the IHBN conformations and probably are not present in bilayer lipids in the monoprotonated cardiolipin. In dCL, the phosphate-phosphate

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4385 hydrogen bonds in the monoprotonated form were more stable than extended conformations, which was unlike the unprotonated form where P-P distances of 6-8 Å were without significant energy differences. Imposing constraints on P-P should therefore lead to dissociation of the proton at neutral pH, in line with the experimental titration curves. It is not clear how significant the phosphate-phosphate bond is in aqueous solution, where the DFT energies suggest both CL and dCL might have such forms. In dCL we note that it is only stabilized by a single hydrogen bond and that these bonds are dynamic. Future work on the full lipid or with realistic treatment of the membrane environment is necessary to determine the role of the IHBN in the presence of other charged lipids and, importantly, in the presence of salt solutions.76 Acknowledgment. We thank Per Siegbahn and Margareta Blomberg for fruitful discussions. This work was supported by grants from the Swedish Research Council (VR), the Carl Trygger Foundation, and the Swedish National Infrastructure for Computing (SNIC 001-09-49) via PDC (www.pdc.kth.se), within the research project “Quantum Chemical Modeling of the Cardiolipin Lipid Headgroup”. Supporting Information Available: Details of methods and results of previous and new MD simulations. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Mitchell, P. Coupling of Phosphorylation to Electron and Hydrogen Transfer by a Chemi-Osmotic Type of Mechanism. Nature 1961, 191, 144– 148. (2) Williams, R. J. The Multifarious Couplings of Energy Transduction. Biochim. Biophys. Acta 1978, 505, 1–44. (3) Kocherginsky, N. Acidic Lipids, H(+)-ATPases, and Mechanism of Oxidative Phosphorylation. Physico-Chemical Ideas 30 Years After P. Mitchell’s Nobel Prize Award. Prog. Biophys. Mol. Biol. 2009, 99, 20–41. (4) Haines, T. H. Anionic Lipid Headgroups as a Proton-Conducting Pathway Along the Surface of Membranes: A Hypothesis. Proc. Natl. Acad. Sci. U.S.A. 1983, 80, 160–164. (5) Teissie, J.; Prats, M.; Soucaille, P.; Tocanne, J. F. Evidence for Conduction of Protons Along the Interface between Water and a Polar Lipid Monolayer. Proc. Natl. Acad. Sci. U.S.A. 1985, 82, 3217–3221. (6) Prats, M.; Teissie, J.; Tocanne, J. Lateral Proton Conduction at Lipid-Water Interfaces and its Implications for the Chemiosmotic-Coupling Hypothesis. Nature 1986, 322, 756–758. (7) Tocanne, J.; Teissie´, J. Ionization of Phospholipids and Phospholipid-Supported Interfacial Lateral Diffusion of Protons in Membrane Model Systems. Biochim. Biophys. Acta 1990, 1031, 111–142. (8) Heberle, J.; Riesle, J.; Thiedemann, G.; Oesterhelt, D.; Dencher, N. A. Proton Migration Along the Membrane Surface and Retarded Surface to Bulk Transfer. Nature 1994, 370, 379–382. (9) Branden, M.; Sanden, T.; Brzezinski, P.; Widengren, J. Localized Proton Microcircuits at the Biological Membrane-Water Interface. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 19766–19770. (10) Few, A. V.; Gilby, A. R.; Seaman, G. V. F. An Electrophoretic Study of Structural Components of Micrococcus Lysodeikticus. Biochim. Biophys. Acta 1960, 38, 130–136. (11) Coulon-Morelec, M. J.; Faure, M.; Marechal, J. Controlled Degradation of Diphosphatidylglycerol (Cardiolipid) in an Acid Medium. Study of the Phosphatide Derivatives obtained. Bull. Soc. Chim. Biol. (Paris) 1962, 44, 171–183. (12) Seddon, J. M.; Kaye, R. D.; Marsh, D. Induction of the LamellarInverted Hexagonal Phase Transition in Cardiolipin by Protons and Monovalent Cations. Biochim. Biophys. Acta 1983, 734, 347–352. (13) Kates, M.; Syz, J. Y.; Gosser, D.; Haines, T. H. PH-Dissociation Characteristics of Cardiolipin and its 2′-Deoxy Analogue. Lipids 1993, 28, 877–882. (14) Hielscher, R.; Wenz, T.; Hunte, C.; Hellwig, P. Monitoring the Redox and Protonation Dependent Contributions of Cardiolipin in Electrochemically Induced FTIR Difference Spectra of the Cytochrome Bc(1) Complex from Yeast. Biochim. Biophys. Acta 2009, 1787, 617–625. (15) Haines, T. H.; Dencher, N. A. Cardiolipin: A Proton Trap for Oxidative Phosphorylation. FEBS Lett. 2002, 528, 35–39. (16) Hubner, W.; Mantsch, H. H.; Kates, M. Intramolecular Hydrogen Bonding in Cardiolipin. Biochim. Biophys. Acta 1991, 1066, 166–174.

4386

J. Phys. Chem. A, Vol. 114, No. 12, 2010

(17) Stewart, L. C.; Kates, M.; Smith, I. C. Synthesis and Characterization of Deoxy Analogues of Diphytanylglycerol Phospholipids. Chem. Phys. Lipids 1988, 48, 177–188. (18) Quinn, P. J.; Kates, M.; Tocanne, J. F.; Tomoaia-Cotisel, M. Surface Characteristics of Phosphatidylglycerol Phosphate from the Extreme Halophile Halobacterium Cutirubrum Compared with those of its Deoxy Analogue, at the air/water Interface. Biochem. J. 1989, 261, 377–381. (19) Teissie, J.; Prats, M.; LeMassu, A.; Stewart, L. C.; Kates, M. Lateral Proton Conduction in Monolayers of Phospholipids from Extreme Halophiles. Biochemistry 1990, 29, 59–65. (20) Tomoaia-Cotisel, M.; Stewart, L. C.; Kates, M.; Zsako, J.; Chifu, E.; Mocanu, A.; Frangopol, P. T.; Noe, L. J.; Quinn, P. J. Acid Dissociation Constants of Diphytanylglycerolphosphorylglycerol-Methylphosphate, and Diphytanylglycerolphosphorylglycerophosphate and its Deoxy Analog. Chem. Phys. Lipids 1999, 100, 41–54. (21) Jiang, F.; Ryan, M. T.; Schlame, M.; Zhao, M.; Gu, Z.; Klingenberg, M.; Pfanner, N.; Greenberg, M. L. Absence of Cardiolipin in the crd1 Null Mutant Results in Decreased Mitochondrial Membrane Potential and Reduced Mitochondrial Function. J. Biol. Chem. 2000, 275, 22387–22394. (22) Hoch, F. L. Cardiolipins and Mitochondrial Proton-Selective Leakage. J. Bioenerg. Biomembr. 1998, 30, 511–532. (23) Schlame, M.; Rua, D.; Greenberg, M. L. The Biosynthesis and Functional Role of Cardiolipin. Prog. Lipid Res. 2000, 39, 257–288. (24) Nichols-Smith, S.; Teh, S. Y.; Kuhl, T. L. Thermodynamic and Mechanical Properties of Model Mitochondrial Membranes. Biochim. Biophys. Acta 2004, 1663, 82–88. (25) Lewis, R. N. A. H.; McElhaney, R. N. The Physicochemical Properties of Cardiolipin Bilayers and Cardiolipin-Containing Lipid Membranes. Biochim. Biophys. Acta 2009, 1788, 2069–2079. (26) Lewis, R. N.; McElhaney, R. N. Surface Charge Markedly Attenuates the Nonlamellar Phase-Forming Propensities of Lipid Bilayer Membranes: Calorimetric and (31)P-Nuclear Magnetic Resonance Studies of Mixtures of Cationic, Anionic, and Zwitterionic Lipids. Biophys. J. 2000, 79, 1455–1464. (27) Egorova, E. M. Determination of the Apparent pK Values of Ionizable Groups in Highly Charged Lipid Membranes. Colloids Surf., A 1996, 110, 47–53. (28) Egorova, E. M. Dissociation Constants of Lipid Ionizable Groups I. Corrected Values for Two Anionic Lipids. Colloids Surf., A 1998, 131, 7–18. (29) Egorova, E. M. Dissociation Constants of Lipid Ionizable Groups II. Changes in Surface pK at Low Ionic Strengths. Colloids Surf., A 1998, 131, 19–31. (30) Dahlberg, M. Polymorphic Phase Behavior of Cardiolipin Derivatives Studied by Coarse-Grained Molecular Dynamics. J. Phys. Chem. B 2007, 111, 7194–7200. (31) Dahlberg, M.; Maliniak, A. Molecular Dynamics Simulations of Cardiolipin Bilayers. J. Phys. Chem. B 2008, 112, 11655–11663. (32) Rog, T.; Martinez-Seara, H.; Munck, N.; Oresic, M.; Karttunen, M.; Vattulainen, I. Role of Cardiolipins in the Inner Mitochondrial Membrane: Insight Gained through Atom-Scale Simulations. J. Phys. Chem. B 2009, 113, 3413–3422. (33) Po¨yry, S.; Ro´g, T.; Karttunen, M.; Vattulainen, I. Mitochondrial Membranes with Mono- and Divalent Salt: Changes Induced by Salt Ions on Structure and Dynamics. J. Phys. Chem. B 2009, 113, 15513–15521. (34) Mennucci, B.; Cammi, R. In Continuum solVation models in chemical physics: from theory to applications; John Wiley & Sons: Chichester, England, and Hoboken, NJ, 2007; pp 619. (35) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. ReV. 1999, 99, 2161–2200. (36) Jang, Y. H.; Goddard, W. A.; Noyes, K. T.; Sowers, L. C.; Hwang, S.; Chung, D. S. PK(a) Values of Guanine in Water: Density Functional Theory Calculations Combined with Poisson-Boltzmann ContinuumSolvation Model. J. Phys. Chem. B 2003, 107, 344–357. (37) Sadlej-Sosnowska, N. Calculation of Acidic Dissociation Constants in Water: Solvation Free Energy Terms. their Accuracy and Impact. Theor. Chem. Acc. 2007, 118, 281–293. (38) Verdolino, V.; Cammi, R.; Munk, B. H.; Schlegel, H. B. Calculation of pK(a) Values of Nucleobases and the Guanine Oxidation Products Guanidinohydantoin and Spiroiminodihydantoin using Density Functional Theory and a Polarizable Continuum Model. J. Phys. Chem. B 2008, 112, 16860–16873. (39) Zahn, D.; Brickmann, J. Quantum-Classical Simulation of Proton Transport Via a Phospholipid Bilayer. Phys. Chem. Chem. Phys. 2001, 3, 848–852. (40) Smondyrev, A. M.; Voth, G. A. Molecular Dynamics Simulation of Proton Transport Near the Surface of a Phospholipid Membrane. Biophys. J. 2002, 82, 1460–1468. (41) Schlame, M.; Brody, S.; Hostetler, K. Y. Mitochondrial Cardiolipin in Diverse Eukaryotes. Comparison of Biosynthetic Reactions and Molecular Acyl Species. Eur. J. Biochem. 1993, 212, 727–735. (42) Goormaghtigh, E.; Huart, P.; Praet, M.; Brasseur, R.; Ruysschaert, J.-. Structure of the Adriamycin-Cardiolipin Complex: Role in Mitochondrial Toxicity. Biophys. Chem. 1990, 35, 247–257.

Dahlberg et al. (43) Yang, P.; Murthy, P. P.; Brown, R. E. Synergy of Intramolecular Hydrogen-Bonding Network in Myo-Inositol 2-Monophosphate: Theoretical Investigations into the Electronic Structure, Proton Transfer, and pKa. J. Am. Chem. Soc. 2005, 127, 15848–15861. (44) Naor, M. M.; Walker, M. D.; Van Brocklyn, J. R.; Tigyi, G.; Parrill, A. L. Sphingosine 1-Phosphate pKa and Binding Constants: Intramolecular and Intermolecular Influences. J. Mol. Graphics Modell. 2007, 26, 519– 528. (45) Vilciauskas, L.; Paddison, S. J.; Kreuer, K. D. Ab Initio Modeling of Proton Transfer in Phosphoric Acid Clusters. J. Phys. Chem. A 2009, 113, 9193–9201. (46) Kaupp, M.; Malkin, V. G.; Bu¨hl, M. In Calculation of NMR and EPR parameters: theory and applications; Kaupp, M., Malkin, V. G., Bu¨hl, M., Eds.; Wiley-VCH: Weinheim, 2004. (47) de Dios, A.; Pearson, J.; Oldfield, E. Secondary and Tertiary Structural Effects on Protein NMR Chemical Shifts: An Ab Initio Approach. Science 1993, 260, 1491–1496. (48) Merrick, J. P.; Moran, D.; Radom, L. An Evaluation of Harmonic Vibrational Frequency Scale Factors. J Phys Chem A 2007, 111, 11683– 11700. (49) Becke, A. D. A New Mixing of Hartree--Fock and Local DensityFunctional Theories. J. Chem. Phys. 1993, 98, 1372–1377. (50) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. ReV. B 1988, 37, 785–789. (51) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Scalmani, G.; Mennucci, B.; Barone, V.; Petersson, G. A.; Caricato, M.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Li, X.; Hratchian, H. P.; Peralta, J. E.; Izmaylov, A. F.; Kudin, K. N.; Heyd, J. J.; Brothers, E.; Staroverov, V. N.; Zheng, G.; Kobayashi, R.; Normand, J.; Sonnenberg, J. L.; Ogliaro, F.; Bearpark, M.; Parandekar, P. V.; Ferguson, G. A.; Mayhall, N. J.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Burant, J. C.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Chen, W.; Wong, M. W.; Pople J. A. ReVision B.05 and E.05, Gaussian, Inc.: Wallingford, CT, 2004. (52) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. ReV. 2005, 105, 2999–3093. (53) Lim, C.; Bashford, D.; Karplus, M. Absolute pKa Calculations with Continuum Dielectric Methods. J. Phys. Chem. 1991, 95, 5610–5620. (54) Camaioni, D. M.; Schwerdtfeger, C. A. Comment on “Accurate Experimental Values for the Free Energies of Hydration of H+, OH-, and H3O+. J. Phys. Chem. A 2005, 109, 10795–10797. (55) Jang, Y. H.; Goddard, W. A.; Noyes, K. T.; Sowers, L. C.; Hwang, S.; Chung, D. S. First Principles Calculations of the Tautomers and pKa Values of 8-Oxoguanine: Implications for Mutagenicity and Repair. Chem. Res. Toxicol. 2002, 15, 1023–1035. (56) Word, J. M.; Lovell, S. C.; Richardson, J. S.; Richardson, D. C. Asparagine and Glutamine: Using Hydrogen Atom Contacts in the Choice of Side-Chain Amide Orientation. J. Mol. Biol. 1999, 285, 1735–1747. (57) Peng, C.; Schlegel, H. B. Combining Synchronous Transit and Quasi-Newton Methods to Find Transition States. Isr. J. Chem. 1993, 33, 449. (58) Peng, C.; Ayala, P. Y.; Schlegel, H. B.; Frisch, M. J. Using Redundant Internal Coordinates to Optimize Equilibrium Geometries and Transition States. J. Comput. Chem. 1996, 17, 49–56. (59) Cheeseman, J. R.; Trucks, G. W.; Keith, T. A.; Frisch, M. J. A Comparison of Models for Calculating Nuclear Magnetic Resonance Shielding Tensors. J. Chem. Phys. 1996, 104, 5497–5509. (60) Adamo, C.; Barone, V. Exchange Functionals with Improved LongRange Behavior and Adiabatic Connection Methods without Adjustable Parameters: The mPW and mPW1PW Models. J. Chem. Phys. 1998, 108, 664–675. (61) Perdew, J. P.; Burke, K.; Wang, Y. Generalized Gradient Approximation for the Exchange-Correlation Hole of a Many-Electron System. Phys. ReV. B 1996, 54, 16533. (62) Ditchfield, R. Self-Consistent Perturbation Theory of DiamagnetismsI. A Gauge-Invariant LCAO Method for N.M.R. Chemical Shifts. Mol. Phys. 1974, 27, 789. (63) Mason, J. In Multinuclear NMR; Plenum Press: New York, 1987; pp 639. (64) Jameson, C. J.; De Dios, A.; Keith Jameson, A. Absolute Shielding Scale for 31P from Gas-Phase NMR Studies. Chem. Phys. Lett. 1990, 167, 575–582.

Modeling of the Cardiolipin Headgroup (65) Schlame, M. Cardiolipin Synthesis for the Assembly of Bacterial and Mitochondrial Membranes. J. Lipid Res. 2008, 49, 1607–1620. (66) Silverstein, R. M.; Webster, F. X.; Kiemle, D. In Spectrometric Identification of Organic Compounds; John Wiley & Sons: Southern Gate, Chichester West Sussex, 2005. (67) Kraffe, E.; Soudant, P.; Marty, Y.; Kervarec, N.; Jehan, P. Evidence of a Tetradocosahexaenoic Cardiolipin in some Marine Bivalves. Lipids 2002, 37, 507–514. (68) Rezanka, T.; Siristova, L.; Melzoch, K.; Sigler, K. Direct ESI-MS Analysis of O-Acyl Glycosylated Cardiolipins from the Thermophilic Bacterium Alicyclobacillus Acidoterrestris. Chem. Phys. Lipids 2009, 161, 115–121. (69) Henderson, T. O.; Glonek, T.; Myers, T. C. Phosphorus-31 Nuclear Magnetic Resonance Spectroscopy of Phospholipids. Biochemistry (N. Y.) 1974, 13, 623–628. (70) Seelig, J. 31P Nuclear Magnetic Resonance and the Head Group Structure of Phospholipids in Membranes. Biochim. Biophys. Acta 1978, 515, 105–140.

J. Phys. Chem. A, Vol. 114, No. 12, 2010 4387 (71) Westheimer, F. Why Nature Chose Phosphates. Science 1987, 235, 1173–1178. (72) Lopez, X.; Schaefer, M.; Dejaegere, A.; Karplus, M. Theoretical Evaluation of pK(a) in Phosphoranes: Implications for Phosphate Ester Hydrolysis. J. Am. Chem. Soc. 2002, 124, 5010–5018. (73) Sadhukhan, S.; Mun˜oz, D.; Adamo, C.; Scuseria, G. E. Predicting Proton Transfer Barriers with Density Functional Methods. Chem. Phys. Lett. 1999, 306, 83–87. (74) Barone, V.; Orlandini, L.; Adamo, C. Proton Transfer in Model Hydrogen-Bonded Systems by a Density Functional Approach. Chem. Phys. Lett. 1994, 231, 295–300. (75) Wohlert, J.; Edholm, O. The Range and Shielding of Dipole-Dipole Interactions in Phospholipid Bilayers. Biophys. J. 2004, 87, 2433–2445. (76) Egorova, E. M.; Tupolev, V. V. Electrochemical Adsorption Method for the Determination of the pK of Ionizable Groups of Lipid Membranes. J. Electroanal. Chem. 1994, 379, 479.

JP9110019