13638
J. Phys. Chem. B 2007, 111, 13638-13644
Molecular Conformations in a Phospholipid Bilayer Extracted from Dipolar Couplings: A Computer Simulation Study Johan Thaning, Carl-Johan Ho1 gberg, Baltzar Stevensson, Alexander P. Lyubartsev, and Arnold Maliniak* DiVision of Physical Chemistry, Arrhenius Laboratory, Stockholm UniVersity, SE-106 91 Stockholm, Sweden ReceiVed: July 6, 2007; In Final Form: September 14, 2007
This paper describes an analysis of NMR dipolar couplings in a bilayer formed by dimyristoylphosphatidylcholine (DMPC). The couplings are calculated from a trajectory generated in a molecular dynamics (MD) simulation based on a realistic atom-atom interaction potential. The analysis is carried out employing a recently developed approach that focuses on the construction of the conformational distribution function. This approach is a combination of two models, the additive potential (AP) model and the maximum entropy (ME) method, and is therefore called APME. In contrast to the AP model, the APME procedure does not require an intuition-based choice of the functional form of the torsional potential and is, unlike the ME method, applicable to weakly ordered systems. The conformational distribution function for the glycerol moiety of the DMPC molecule derived from the APME analysis of the dipolar couplings is in reasonable agreement with the “true” distributions calculated from the trajectory. Analyses of dipolar couplings derived from MD trajectories can, in general, serve as guidelines for experimental investigations of bilayers and other complex biological systems.
Introduction Lipids are the main components of biomembranes, which in turn regulate the bioactivities and transport of various molecular species in living cells. These amphiphilic molecules in water also display a variety of complex mesophases such as bicelles, ribbons, vesicles, and different lamellar structures to name a few.1 This rich polymorphism originates from the hydrophobic effect, which is closely connected to the amphiphilic nature of the lipid molecule. The structures and the temperature region of the phases formed by lipids are governed by a delicate balance between the chemical composition and the molecular properties such as the number of lipid tails and the polar nature of the headgroup. The most studied lipid systems, sometimes called the workhorses in bilayer research, are dipalmitoylphosphatidylcholine (DPPC) and dimyristoylphosphatidylcholine (DMPC) bilayers. In the present work, we focus on the bilayers formed by DMPC (see Figure 1), which consists of a nonionic but polar phosphatidylcholine headgroup, whereas the hydrophobic part is made of two tails, each being 14 carbons long. During the past few years, a dramatic increase in computer power has led to a significant increase of interest in computer simulations of complex chemical systems. Many investigations of bilayers have been performed with various degrees of approximation in the interaction model.2-12 Lipid bilayers have been studied extensively using experimental methods, and therefore a large data set is available for comparisons with the results from computer simulations. Nuclear magnetic resonance (NMR) spectroscopy is a powerful experimental method for investigations of liquid crystals in general and lyotropic mesophases in particular. Deuterium (2H) NMR has been employed extensively for studies of the * Author for correspondence. E-mail:
[email protected].
Figure 1. Schematic of DMPC.
molecular order, structure, and dynamics in lipid bilayers.13-19 Recently, several investigations of molecular self-diffusion have been reported.20-24 Deuterium spectra, though easy to interpret, reflect only the motion of the C-D vector and are therefore relatively insensitive to molecular conformations. Carbon-13 NMR spectroscopy has several advantages for studying lipid bilayers, because anisotropic spin interactions, such as chemical shifts or dipolar couplings, provide important information on molecular structure and ordering.25-34 In particular, the throughspace magnetic dipole-dipole couplings have proven to be of great importance. These interactions depend on the spin-spin distances and on the orientations of the internuclear vectors with respect to the external magnetic field. This means that the dipolar coupling is a valuable probe of long-range order and molecular structure. The measurement of 1H-13C and 1H-1H dipolar couplings in macromolecules dissolved in aqueous dilute liquid crystals is an emerging area in structural biology and biochemistry. The dipole-dipole interactions are used together with the traditional observables (nuclear Overhauser effects and Jcouplings) in the structure determination of proteins,35-37 nucleic acids,38,39 and carbohydrates.40-43 To extract useful information about conformation and order from the experimental dipolar couplings in a flexible molecule, we need a model that takes into account both molecular tumbling and internal degrees of freedom. Two approaches that have been frequently used for interpretations of dipolar couplings in bulk
10.1021/jp075278t CCC: $37.00 © 2007 American Chemical Society Published on Web 11/13/2007
Dipolar Couplings in Lipid Bilayers liquid crystals are the additive potential (AP) model44 and the maximum entropy (ME) approach.45 Both models have been successfully applied in experimental and computational studies of the orientational order and molecular conformations in mesogenic molecules.46-50 These models suffer, however, from serious limitations. The AP approach requires a priori knowledge of the functional form of the torsion potential relevant for the molecular fragment. The ME procedure gives the flattest possible conformational distribution consistent with the experimental data set, which results in an incorrect description of systems with low orientational order.45,48 These two limitations are obviously relevant for investigations of complex macromolecules in liquid crystals: the torsion potential function is frequently unknown, and the orientational order is usually low. In fact, the low order problem in the ME approach has been partially solved by introducing conformational constrains,51-54 which results in a distribution with the functional form corresponding to the AP model. Recently, we proposed a new approach for the analysis of NMR parameters in general and dipolar couplings in particular.55 This procedure, derived for the low-order limit, is a combination of the AP and ME methods and is denoted APME. In analogy with the AP model, the APME approach also requires the knowledge of the internal potential function. It is, however, always consistent with the ME principle and therefore no intuition-based assumptions are necessary. The APME approach was demonstrated and used for interpretation of the dipolar interactions determined experimentally,43,55,56 and calculated from a trajectory generated in a molecular dynamics (MD) computer simulation of a nematic liquid crystal.57 The APME approach was subsequently generalized,58 leading to extended limits of validity for the model. In the present study, the APME analyses were performed on a set of dipolar couplings calculated from an MD trajectory11 of a phospholipid bilayer formed by DMPC. Using computer simulations to generate dipole-dipole interactions has obvious strengths: (i) rigorous tests of theoretical approaches for analyzing couplings are rewarding to perform since the distribution functions derived from the analyses can be immediately compared with the simulated counterparts, (ii) all possible dipolar interactions may be calculated from the trajectory, and (iii) the minimal required set and the most sensitive couplings can be selected. Three sets of dipolar interactions were used for the analysis: 13C-13C, 1H-13C, and 1H-1H couplings. In addition, 17O-17O, 13C-17O, and 1H-17O were employed to check the consistency of the procedure. In the conformational analysis, we have considered the two torsion angles (φ1 and φ2 in Figure 1) to be related to the glycerol fragment of the lipid molecule. The structure of this molecular fragment has been previously investigated in the lamellar phase by means of NMR and computer modeling.19,27 Computational Methods Simulation of the Bilayer. Details of the computer simulation, the parameters of the force fields, and the properties of the lamellar liquid crystal are reported in a recent publication;11 here we only summarize the essential technical details. The MD simulation of the bilayer was performed on a system consisting of 128 (64 × 2) DMPC lipids and 3655 water molecules.11 The initial configuration was taken from a well-equilibrated (200 ns) simulation in the liquid crystalline phase carried out at T ) 303 K.11 In the present work, the temperature was set to T ) 310 K, and the simulation continued for 50 ns, from which the last 45 ns were used for the analysis.
J. Phys. Chem. B, Vol. 111, No. 48, 2007 13639 The lipid force field parameters for bonded and nonbonded interactions and atomic partial charges are based on the GROMOS force field.59,60 The united atom model was used for the CH, CH2, and CH3 groups in DMPC lipids. Previous studies of phospholipid bilayers have shown that the united-atom and all-atom force fields yield essentially identical results, whereas the united-atom model results in a 3-fold saving of the computer time. The bond lengths were kept constant at their equilibrium values, and a time step of 5 fs was used in the simulation. For water, the simple point charge (SPC) model was employed. The simulation was carried out using Gromacs simulation package, version 3.2.61 The temperature of the liquid crystalline phase of the DMPC bilayer was regulated separately for lipids and water using the Nose-Hoover thermostat.62 The pressure, set to 1 atm, was regulated anisotropically, using the Parrinello-Rahman barostat,63 allowing the x-, y-, and z-dimensions of the simulation box to change independently. The long-range electrostatic forces were treated by the particle mesh Ewald summation (PME) algorithm.64 Analysis of the Dipolar Couplings. The through-space magnetic dipolar coupling between spins k and l, with magnetogyric ratios γk and γl, is (in Hz) given by
dkl ) -
µ0 16π
2
γkγlp〈(3 cos2 θkl - 1)rkl-3〉
(1)
where rkl is the spin-spin distance, and θkl is the angle between the spin-spin vector and the external field. In the analysis, we assume that the director (defined as the normal to the bilayer plane) and the magnetic field coincide. The angular bracket in eq 1 indicates that the dipolar couplings are averaged over both molecular tumbling and internal bond rotations. In order to derive the molecular conformation, the dipolar interaction for a fixed conformational state Φ is expressed as65
dkl(Φ) ) bkl(Φ)〈D0,02(ΩLP(Φ))〉 ) bkl(Φ)D0,02(ΩLD) × 2
〈D0,m2(ΩDM)〉Dm,02(ΩMP(Φ)) ∑ m)-2
(2)
where the dipole-dipole coupling constant is defined as bkl(Φ) ) -µ0γkγlp/8π2rkl3(Φ) and Dn,m2(Ω) is the second-rank Wigner function.66 The transformation between the principal axis system (P) of the dipolar interaction and the laboratory frame (L) is performed explicitly using three successive rotations. Transformation between frame P and the molecular frame M is characterized by ΩMP(Φ), which, together with bkl(Φ), reflects the direct conformational dependence of the coupling. The Euler angles β and γ that specify relative orientations of molecular and director (D) frames are denoted ΩDM. Note that the averaging in eq 2 refers only to molecular orientations, and the order parameters, 〈D0,m2(ΩDM)〉, are conformation dependent. The relative orientation of the director and the magnetic field is given by ΩLD. The assumption of the director and magnetic field (laboratory frame L) being parallel results in D0,02(ΩLD) ) 1. The motionally averaged dipolar couplings in a flexible molecule, such as DMPC, are given by
dkl )
∫dkl(Φ)P(Φ)dΦ
(3)
where Φ ) {φ1,φ2} denotes, in our case, a conformational state defined by the two torsion angles related to the glycerol part of the DMPC molecule. The torsion angle probability, P(Φ), is
13640 J. Phys. Chem. B, Vol. 111, No. 48, 2007
Thaning et al.
obtained from the ensemble average over the singlet orientational distribution function (ODF) P(β,γ,Φ):
P(Φ) ) Z-1
∫∫P(β,γ,Φ) sin βdβdγ
(4)
where β and γ are Euler angles specifying the relative orientations of molecular and director frames, and Z is a normalization factor. In eq 2, dkl(Φ) is the dipolar coupling for a conformation with probability P(Φ). The central task in the interpretation of the dipolar interactions is therefore to determine the torsion angle distribution function. The expression for the ODF and thus P(Φ) is derived by combining the AP model44 and the ME approach,45 resulting in PAPME(Φ). The details of the derivation of the APME approach are collected in the Appendix. The final expression for the distribution function takes the following form:57
PAPME(φ1,φ2) ) Z′-1 exp{-
λkldkl(φ1,φ2)} ∑ kl
(5)
where the conformational parameter Φ was replaced by the explicit torsion angles φ1 and φ2, Z′ is a normalization constant, and λkl are adjustable parameters that can be determined by numerical fitting of calculated couplings with those obtained from the MD trajectory. The conformation-dependent but orientationally averaged dipolar coupling dkl(φ1,φ2) is given by eq A8 (Appendix). In contrast to the ME procedure, the APME approach exhibits a very weak dependence of the molecular conformations on the orientational order (see Figure 2 in ref 57). It is so because the decreased order, described by the parameters in eq A8, is compensated for by a corresponding increase of the λ parameters. In general, a torsion angle is defined by four atoms (A, B, C, and D) and can be seen as an angle between two planes, ABC and BCD, respectively. For the APME analysis of the glycerol moiety in the DMPC molecule, local coordinate systems in the three rigid subunits (Cg1, Cg2, and Cg3) related to the two torsion angles φ1 and φ2 are required. The local coordinate frames are located with the z axes parallel to the Cg1-Cg2, Cg2Cg3, and Cg3-Cg2 bonds, whereas the y axes are normal to the planes Og1-Cg1-Cg2, Cg1-Cg2-Cg3, and Og3-Cg3-Cg2. In principle, each rigid fragment requires five linearly independent dipolar couplings (order parameters) to determine the local order tensor. If, however, the z axes of the two rigid subunits centered on Cg2 and Cg3 atoms coincide, the number of order parameters required for the analysis can be reduced by one. In the APME analysis, the local orientational order is described by the iab parameters (eq A8). The potential model employed in the MD simulation is based on united atoms, which means that no hydrogen atoms were included. Thus, the hydrogen positions, required for 1H-1H, 1H-13C, and 1H-17O dipolar couplings, were calculated from the coordinates of the heavy atoms using valence angles from the trajectory and assuming the C-H distances to be 1.095 Å. The absence of protons in the MD simulation implies that vibrational contributions to the couplings are partially neglected. The strategy used in the analysis is to fit the dipole-dipole interactions collected in Table 1 to eq 3. The APME distribution function PAPME(φ1,φ2) is constructed using i and λ parameters (eqs 5 and A8), which describe the orientational order and molecular conformations, respectively. We will return to the physical meaning of these parameters and the choice of relevant couplings in the following section. The fitting was performed using a computer code written in-house based on the MATLAB
TABLE 1: Dipolar Couplings (in Hz) Calculated from the MD Trajectory (dtraj kl ) and Derived from the APME Analysis (dAPME )a kl coupling
dtraj kl
dAPME kl
dtraj kl
coupling
dAPME kl
Cg1-Cg2 571 (13) 571 Hg1a-Hg2 1117 (58) 1117 Cg2-Cg3 -1144 (54) -1144 Hg1a-Hg3a 396 (38) 396 Cg1-Hg1a -5021 (354) -5021 Hg1a-Hg3b 216 (44) 216 Cg1-Hg1b -139 (731) -139 Hg1b-Hg2 1526 (18) 1526 Cg2-Hg2 4320 (241) 4320 Hg1b-Hg3a -408 (53) -408 Cg3-Hg3a 6090 (186) 6090 Hg1b-Hg3b -276 (16) -276 Cg3-Hg3b 4344 (275) 4344 Hg2-Hg3b -4000 (212) -4000 Cg1-Hg2 843 (24) 843 Hg2-Hg3a -1021 (193) -1021 Cg2-Hg1a 852 (26) 852 Cg2-Hg3b -1112 (35) -1112 Cg3-Hg2 -740 (106) -740 Hg1a-Hg1b -5958 (480) -5958 Hg3a-Hg3b 6450 (170) 6450 Cg1-Cg3 -136 (7) -129 Og3- Cg2 112 (4) 110 Cg1-Hg3a -100 (25) -123 Og3- Cg3 -33 (21) -25 Cg1-Hg3b -218(10) -223 Og1- Hg1a 547 (68) 535 Cg2-Hg1b 456 (34) 457 Og1- Hg1b -282 (18) -277 Cg2-Hg3a -752 (92) -722 Og1- Hg2 -161 (9) -177 Cg3-Hg1a 230 (18) 233 Og1- Hg3a 52 (9) 67 Cg3-Hg1b -190 (10) -179 Og1- Hg3b 82 (6) 57 Og1-Og2 21 (1) 16 Og2- Hg1a -161 (6) -176 Og1-Og3 -23 (1) -17 Og2- Hg1b -284 (6) -285 Og2-Og3 -18 (1) -19 Og2- Hg2 -388 (69) -378 Og1- Cg1 64 (1) 64 Og2- Hg3a 472 (18) 469 Og1- Cg2 -21 (4) -20 Og2- Hg3b 207 (12) 198 Og1- Cg3 27 (2) 20 Og3- Hg1a 600 (22) 705 Og2- Cg1 -99 (1) -95 Og3- Hg1b 282 (11) 253 Og2- Cg2 -120 (33) -104 Og3- Hg2 111 (16) 120 Og2- Cg3 107 (5) 108 Og3- Hg3a -393 (40) -366 Og3- Cg1 126 (3) 121 Og3- Hg3b -463 (290) -450 a Numbers in parentheses are estimated errors given as σtraj (eq 7). kl Note that the dipolar couplings in the upper part of the table were used in the AMPE analysis, whereas the couplings in the lower part provided a test of the distribution function (see text for details).
subroutine fminu.67 The fitting program minimizes the quality factor
Q)
∑ kl
(
)
APME dtraj kl - dkl
σtraj kl
2
(6)
APME where dtraj are dipolar couplings calculated from the kl and dkl trajectory and derived from the APME analysis, respectively. The standard deviation of each coupling is defined as
2 1/2 traj traj σtraj kl ) 〈(〈dkl 〉N - 〈dkl 〉N,t) 〉t
(7)
where N is the number of molecules, and t is time (length of the trajectory). In practice, the trajectory was divided into nine 5 ns blocks, which were averaged separately. Note that only the couplings in the upper part of the table are used in the construction of the distribution function (yielding Q ) 0); the other couplings are employed for the check of the procedure. Results and Discussion The conformational analysis of the dipolar couplings was restricted to the glycerol moiety of the DMPC molecule. Two torsion angles are relevant for this segment (see Figure 1): φ1 defined by atoms Cg1-Cg2-Cg3-Og3, and φ2 given by Og1Cg1-Cg2-Cg3. The normalized conformational distribution function derived from the MD trajectory Ptraj(φ1,φ2) is shown in Figure 2a, while the contour plot of the distribution is displayed in Figure 3. It is obvious that two regions of
Dipolar Couplings in Lipid Bilayers
J. Phys. Chem. B, Vol. 111, No. 48, 2007 13641 TABLE 2: Local Order Parameters Derived from Experiments and Calculated from the Trajectory, 2 Straj CH ) 〈3 cos φ - 1〉/2, Where φ Is the Angle between the CH Vector and the Normal to the Bilayer (Director) vector
a Straj CH
|SCH|b
|SCD|c
Cg1Hg1a Cg1Hg1b Cg2Hg2 Cg3Hg3a Cg3Hg3b
0.003 0.22 -0.19 -0.26 -0.19
0 0.16 0.20 0.23
0 0.16 0.19 0.21
a The errors are on the order of 1/xN, where N ) 128. b Taken from ref 29. c Taken from ref 69.
TABLE 3: Parameters Obtained from the APME Analysis (iab, where i ) g1, g2, g3 and λkl) of the Dipolar Couplings in DMPC (Table 1)a
Figure 2. The normalized conformational distribution functions (a) calculated from the trajectory generated in an MD computer simulation (Ptraj(φ1,φ2)), and (b) derived from the APME analysis of the dipolar couplings (PAPME(φ1,φ2)).
parameter
APME analysis
g1 zz g1 xx-yy g1 xy g1 xz g1 yz g2 zz g2 xx-yy g2 xy g2 xz g2 yz g3 zz g3 xx-yy g3 xy g3 xz g3 yz λg1a-g2
-0.3491 3.3338 1.5477 -0.6476 -1.6434 0.9446 -1.0298 0.3159 -0.9282 1.7506 0.9446 1.0741 -0.7853 -1.8798 1.1460 1.930 0.890 -8.914 -4.352 0.166 -1.342 1.834 0.449
λg1a-g3a λg1a-g3b λg1b-g2 λg1b-g3a λg1b-g3b λg2-g3a λg2-g3b Figure 3. The contour plots of the conformational distribution functions (Figure 2): Ptraj(φ1,φ2) - gray, and PAPME(φ1,φ2) - black.
Ptraj(φ1,φ2) are highly populated. These conformational regions are in agreement with results obtained in a study where deuterium NMR was combined with molecular mechanics calculations19 and derived from X-ray data.68 We start the analysis of the dipolar interactions by considering the distance (conformation)-independent dipolar couplings (left columns of the upper part of Table 1), which only reflect the orientational order of the three rigid fragments. Strictly speaking, no dipolar couplings are completely conformation independent because the ordering matrix is unique for every conformational state. Therefore, the one-bond 1H-13C couplings also contain indirect structural information. The order parameters relevant for one-bond 1H-13C dipolar couplings, Straj CH, calculated from the trajectory and derived from 2D MAS (R-PDLF)29 and deuterium NMR experiments,69 are collected in Table 2. The agreement is reasonable, which, considering all the approximations made in the simulation procedure, indicates that our force field adequately describes the structure and order of the lipid bilayer. Note that the two order parameters related to the Cg3 group are indistinguishable in experiments, but slightly different in the simulation (with an average corresponding to the experimental value).
a Note that all the λ parameters were derived from the 1H-1H kl couplings and were scaled with a factor of 103.
The 13 dipolar couplings in the first three columns of the upper section of Table 1 are used for determination of the orientational order, whereas the eight couplings in the following columns have strong conformation-dependent spin-spin distances and are therefore used for extraction of the λkl parameters. The number of conformation-dependent dipolar couplings (and thus λkl parameters) was chosen on the basis of several initial tests that guided our final choice of the minimal set required for the analysis. The normalized conformational distribution function derived from the APME analysis, PAPME(φ1,φ2), is shown in Figure 2b, with the contour plot displayed in Figure 3. The adjustable parameters obtained from the APME procedure and used for the calculation of PAPME(φ1,φ2) are collected in g3 Table 3. Note, that g2 zz and zz are identical because the two z axes coincide; furthermore, g1 zz can be expressed as a linear APME ) combination of the five elements: g2 ab. The calculated (dkl dipolar couplings were obtained using eq 3, substituting P(Φ) APME in with PAPME(φ1,φ2). The agreement between dtraj kl and dkl the upper part of Table 1 is perfect, which is expected for a situation where the number of couplings is equal to the number of fitting parameters. In order to check the consistency of the APME procedure, 34 additional dipole-dipole interactions were calculated from the trajectory (dtraj kl ) and compared with those
13642 J. Phys. Chem. B, Vol. 111, No. 48, 2007
Thaning et al. of valence angles, that results from an analysis carried out using a single value of the θ deviates from the true conformational distribution, Ptraj(φ1,φ2). Before closing this section, we comment on the validity limit of the APME procedure. Previously, we estimated57 that the APME approach is valid up to a limit of the orientational order parameter of S ≈ 0.6. It is so because in the derivation of APME, a low-order approximation was used (see eqs A6 and A7 in the Appendix). The molecular order of the lipids in our bilayer system is S ) 0.51 and therefore close to the validity limit. However, in view of the extension of these limits suggested in ref 58, we believe that the APME procedure can be used in the analysis of the dipolar couplings in DMPC. Conclusions
Figure 4. The probability distribution of the valence angle θ (Cg1Cg2-Cg3) derived from the MD trajectory.
Figure 5. The quality factor Q (eq 6) dependence of the valence angle θ (Cg1-Cg2-Cg3) obtained using the APME analysis of the dipolar couplings. The line is drawn to guide the eye. Note that all (55) couplings collected in Table 1 are used for calculation of the quality factor.
calculated (dAPME ) using the distribution function PAPME(φ1,φ2) kl in eq 3. These couplings are collected in the lower part of Table 1, and the agreement between the dipolar couplings calculated from the trajectory and those derived from the APME analysis is very good. We note that the agreement between the two distributions Ptraj(φ1,φ2) and PAPME(φ1,φ2) in Figures 2 and 3 is reasonable but not perfect. In particular, the conformational probabilities exhibit two highly populated regions with similar relative volumes and positions. An obvious observation is that the Ptraj(φ1,φ2) distribution is significantly broader. One possible explanation for the discrepancy between the two distributions may be the fact that we did not attempt to optimize the geometrical parameters of the DMPC molecule. Yet, another source of error can be the neglect of the vibrational contributions to the couplings. It is known that, in addition to torsion angles, dipolar interactions are sensitive to bending angles.70,71 Preliminary tests indicated that the couplings are particularly sensitive to one bending angle, namely, Cg1-Cg2-Cg3, denoted θ. The probability distribution for that angle was calculated from the trajectory and is shown in Figure 4. The distribution is relatively broad, with full width at half-maximum of fwhm ) 8°, and a mean value of θ ) 106°. The APME analysis was therefore carried out assuming θ ) 106°. The dependence of quality factor Q (eq 6) on θ displayed in Figure 5 exhibits a minimum at 106-107° and rises steeply for smaller and larger angles. Thus, it is not surprising, given the broad distribution
We have demonstrated a recently developed method for determination of conformational distribution functions P(φ1,φ2) from NMR dipolar couplings. The method is a hybrid of the AP model and the ME technique and is therefore referred to as the APME approach. The NMR dipole-dipole interactions were derived from a trajectory generated in an MD computer simulation of a bilayer formed by DMPC. The simulation was carried out employing a detailed atom-atom interaction model. The dipolar couplings were subsequently used for investigations of the two torsion angles, φ1 and φ2, in the DMPC molecule related to the glycerol moiety. The aim of this study was to estimate the molecular structure by constructing torsion angle distribution functions and to compare these with the “true” distributions calculated from the trajectory. In order to give guidelines for experimental investigations, a minimal set of the couplings was employed in the analysis. The conformational distribution derived from the APME analysis, PAPME(φ1,φ2), is in reasonable agreement with the distribution calculated directly from the trajectory, Ptraj(φ1,φ2). A strong dependence of the conformational distributions on the central valence angle of the glycerol fragment was observed. In a forthcoming study, we intend to modify the APME approach to include the vibrational degrees of freedom. Our results indicate that the APME procedure is a possible and attractive tool for the interpretation of experimental results in bilayers, liquid crystals, and other complex chemical systems. In particular, several sets of experimental dipolar couplings determined in lipid bilayers were recently reported in the literature.29,30,34,72 Acknowledgment. This work was supported by the Swedish Research Council (VR), and the Carl Trygger Foundation. Appendix APME Procedure: A Hybrid Model Based on ME and Molecular Field Theory. The crucial step in the analysis of dipolar couplings in flexible molecules involves construction of the torsion angle distribution function for the molecular fragment under consideration. The dipolar interaction for a fixed conformational state Φ is given by
dkl(Φ) ) bkl(Φ)〈D0,02(ΩLP(Φ))〉 ) bkl(Φ)D0,02(ΩLD) × 2
〈D0,m2(ΩDM)〉Dm,02(ΩMP) ∑ m)-2
(A1)
where the dipole-dipole coupling constant is defined as bkl(Φ) ) -µ0γkγlp/8π2rkl3(Φ), and Dn,m2(Ω) is the second-rank Wigner function.66 Note that the averaging in eq A1 refers only to molecular orientations, which implies that the order parameters
Dipolar Couplings in Lipid Bilayers
J. Phys. Chem. B, Vol. 111, No. 48, 2007 13643
〈D0,m2(ΩDM)〉 are conformation dependent. The transformation between the principal axis system (P) of the dipolar interaction and the laboratory frame (L) is performed explicitly using three successive rotations. Transformation between frame P and the molecular frame is characterized by ΩMP, whereas the Euler angles β and γ that specify relative orientations of molecular and director frames are denoted ΩDM. Finally, the orientation of the director in the magnetic field is given by ΩLD. The AP Model. The AP model rests on equilibrium statistical mechanics, and the mean potential U(β,γ,Φ) is related to the singlet ODF by44
PAP(β,γ,Φ) ) Z-1 exp{-U(β,γ,Φ)}
(A3)
The conformational potential Uint(Φ) is independent of molecular orientation, and it is usually expressed as a cosine series. The potential of mean torque, Uext(β,γ,Φ), is written as 2
Uext(β,γ,Φ) ) -
∑ n)-2
2,n(Φ)D0,n2*(0,β,γ)
(A4)
where the interaction parameters 2,n(Φ) depend on the segmental, anisotropic interactions. The conformation dependence of these parameters is achieved by expressing them in terms of Φ-independent coefficients, i2,p, which represent the interactions of the ith rigid subunit with the liquid crystalline field. Thus, 2
2,n(Φ) )
∑i p)-2 ∑
i 2,p Dp,n2(ΩiM)
(A5)
where ΩiM defines the orientation of the ith subunit in the molecular frame. We will now consider a situation where the orientational order is low, which in turn implies that the potential of mean torque, Uext(β,γ,Φ), is small. In this case, the ODF can be Taylor expanded and truncated after the second term. The analytical averaging over the molecular orientations gives the distribution function corresponding to the low-order limit of the AP model, P′AP(Φ)55-57
P′AP(Φ) ) Z′-1
∫∫ exp{-Uint(Φ)}[1 Uext(β,γ,Φ)] sin βdβdγ
) Z′-1
∫∫ exp{-Uint(Φ)}[1 +
∫∫ exp{-Uint(Φ)}[1 Uext(β,γ,Φ)]D0,m2(0,β,γ) sin βdβdγ
)
1
[D0,m2(0,β,γ) + ∫∫ 4π
2
2,n(Φ)D0,n2*(0,β,γ)D0,m2(0,β,γ)] sin βdβdγ ∑ n)-2 )
1 4π
(A2)
where U(β,γ,Φ), in units of RT, is written as the sum
U(β,γ,Φ) ) Uint(Φ) + Uext(β,γ,Φ)
〈D0,m2(ΩDM)〉 ≈ Q(Φ)-1
[
2
0+
∑ 2,n(Φ) n)-2
dkl(Φ) ) )
bkl(Φ) 5 bkl(Φ) 5
)
bkl(Φ) 5
2
2,n(Φ)D0,n2*(0,β,γ)] sin βdβdγ
) Z′-14π exp{-Uint(Φ)} ) Piso(Φ)
2,m(Φ) 5
(A7)
2
i Dn,m2(ΩiM)Dm,02(ΩMP) ∑i n)-2 ∑ m)-2 ∑ 2,n 2
i Dn,02(ΩiP) ∑i n)-2 ∑ 2,n
cos θakl cos θbkliab ∑i ∑ a,b
(A8)
where the second equality was obtained using the closure properties of the Wigner functions,66 and D0,02(ΩLD) ) 1 was assumed. In the last equality, Cartesian coordinates were used, a,b ) (xi,yi,zi) refers to the coordinate frame fixed in the ith rigid unit, and θakl is the conformation-dependent angle between the internuclear vector and the a axis. Note that the orientational order in eq A8 is described by the i coefficients. Combination of the AP and ME Approaches: The APME Model. Assuming low orientational order, the conformational distribution function P′AP(Φ) was related, in eq A6, to the internal potential Uint(Φ). The conformational distribution function for the glycerol fragment in the DMPC molecule depends on the torsion angles φ1 and φ2. The choice of the functional form for the potential functions for these angles is not obvious. The least biased approach that requires minimum a priori information is the ME method.45 The ODF, PME(β,γ,φ1,φ2), can be expressed as
2
∑
5
)
The Kronecker delta δn,m is a consequence of the orthogonality of the Wigner functions,66 and Q(Φ) ) 4π exp{-Uint(Φ)} is the conformation-dependent partition function. The glycerol fragment in the DMPC molecule consists of three rigid subunits corresponding to Cg1, Cg2, and Cg3. Thus, using eqs A1, A5, and A7, we can define a conformationdependent but orientationally averaged dipolar coupling:
PME(β,γ,φ1,φ2) ) Z-1 exp{-
n)-2
]
4πδn,m
λkldkl(β,γ,φ1,φ2)} ∑ kl
(A9)
where dkl(β,γ,φ1,φ2) is the orientation- and conformationdependent dipolar coupling
(A6)
where Piso(Φ) is related to the internal potential Uint(Φ), and Z′ and Z′/4π are the normalization constants for P′AP(Φ) and Piso(Φ), respectively. Thus, for a weakly ordered system, P′AP(Φ) ) Piso(Φ). The integral over the second term in eq A6 vanishes for all n due to properties of the Wigner functions.66 Using the low-order approximation and eq A4, we can express the conformation-dependent order parameters:
dkl(β,γ,φ1,φ2) ) bkl(φ1,φ2)D0,02(ΩLD) × 2
D0,m2(0,β,γ)Dm,02(ΩMP) ∑ m)-2
(A10)
We use the functional form of the ME approach to describe the conformational dependence of the dipolar couplings. Combination of the low orientational order approximation of the AP
13644 J. Phys. Chem. B, Vol. 111, No. 48, 2007
Thaning et al.
model with the ME distribution function yields the APME distribution function57 given by
PAPME(φ1,φ2) ) Z-1 exp{-
λkldkl(φ1,φ2)} ∑ kl
(A11)
where Z is a normalization constant, and the conformationdependent but orientationally averaged dipolar coupling dkl(φ1,φ2) is defined in eq A8. References and Notes (1) Evans, D. F., Wennerstro¨m, H. The Colloidal Domain where Physics, Chemistry, and Biology Meet; Wiley-VCH: New York, 1999. (2) Marrink, S. J.; Lindahl, E.; Edholm, O.; Mark, A. E. J. Am. Chem. Soc. 2001, 123, 8638-8639. (3) Bandyopadhyay, S.; Shelley, J. C.; Klein, M. L. J. Phys. Chem. B 2001, 105, 5979-5986. (4) den Otter, W. K. J. Chem. Phys. 2005, 123, 214906-1-214906-9. (5) Patra, M.; Salonen, E.; Terama, E.; Vattulainen, I.; Faller, R.; Lee, B. W.; Holopainen, J.; Karttunen, M. Biophys. J. 2006, 90, 1121-1135. (6) Houndonougbo, Y.; Kuczera, K.; Jas, G. S. Biochemistry 2005, 44, 1780-1792. (7) de Vries, A. H.; Chandrasekhar, I.; van Gunsteren, W. F.; Hu¨nenberger, P. H. J. Phys. Chem. B 2005, 109, 11643-11652. (8) Kranenburg, M.; Smit, B. J. Phys. Chem. B 2005, 109, 65536563. (9) MacCallum, J. L.; Tieleman, D. P. J. Am. Chem. Soc. 2006, 128, 125-130. (10) Knecht, V.; Mark, A. E.; Marrink, S. J. J. Am. Chem. Soc. 2006, 128, 2030-2034. (11) Ho¨gberg, C.-J.; Lyubartsev, A. P. J. Phys. Chem. B 2006, 110, 14326-14336. (12) Kasson, P. M.; Kelley, N. W.; Singhal, N.; Vrljic, M.; Brunger, A. T.; Pande, V. S. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 11916-11921. (13) Davis, J. H. AdV. Magn. Reson. 1989, 13, 195-223. (14) Seelig, J. Q. ReV. Biophys. 1977, 10, 353-418. (15) Bonev, B. B.; Chan, W. C.; Bycroft, B. W.; Roberts, G. C. K.; Watts, A. Biochemistry 2000, 39, 11425-11433. (16) Lehnert, R.; Eibl, H. J.; Mu¨ller, K. J. Phys. Chem. B 2004, 108, 12141-12150. (17) Shaikh, S. R.; Cherezov, V.; Caffrey, M.; Soni, S. P.; LoCascio, D.; Stillwell, W.; Wassall, S. R. J. Am. Chem. Soc. 2006, 128, 53755383. (18) Aussenac, F.; Lavigne, B.; Dufourc, E. J. Langmuir 2005, 21, 7129-7135. (19) Aussenac, F.; Laguerre, M.; Schmitter, J. M.; Dufourc, E. J. Langmuir 2003, 19, 10468-10479. (20) Filippov, A.; Ora¨dd, G.; Lindblom, G. Biophys. J. 2004, 86, 891896. (21) Ora¨dd, G.; Westerman, P. W.; Lindblom, G. Biophys. J. 2005, 89, 315-320. (22) Filippov, A.; Ora¨dd, G.; Lindblom, G. Biophys. J. 2003, 84, 30793086. (23) Ora¨dd, G.; Lindblom, G. Magn. Reson. Chem. 2004, 42, 123131. (24) Rajamoorthi, K.; Petrache, H. I.; McIntosh, T. J.; Brown, M. F. J. Am. Chem. Soc. 2005, 127, 1576-1588. (25) Koenig, B. W.; Gawrisch, K. J. Phys. Chem. B 2005, 109, 75407547. (26) Soubias, O.; Reat, V.; Saurel, O.; Milon, A. J. Magn. Reson. 2002, 158, 143-148. (27) Semchyschyn, D. J.; Macdonald, P. M. Magn. Reson. Chem. 2004, 42, 89-104. (28) Dvinskikh, S. V.; Castro, V.; Sandstro¨m, D. Magn. Reson. Chem. 2004, 42, 875-881. (29) Dvinskikh, S. V.; Castro, V.; Sandstro¨m, D. Phys. Chem. Chem. Phys. 2005, 7, 607-613. (30) Dvinskikh, S. V.; Castro, V.; Sandstro¨m, D. Phys. Chem. Chem. Phys. 2005, 7, 3255-3257. (31) Dvinskikh, S. V.; Yamamoto, K.; Ramamoorthy, A. Chem. Phys. Lett. 2006, 419, 168-173. (32) Gross, J. D.; Warschawski, D. E.; Griffin, R. G. J. Am. Chem. Soc. 1997, 119, 796-802. (33) Lu, J. X.; Damodaran, K.; Lorigan, G. A. J. Magn. Reson. 2006, 178, 283-287. (34) Dvinskikh, S.; Du¨rr, U.; Yamamoto, K.; Ramamoorthy, A. J. Am. Chem. Soc. 2006, 128, 6326-6327.
(35) Tjandra, N.; Bax, A. Science 1997, 278, 1111-1114. (36) Meiler, J.; Peti, W.; Griesinger, C. J. Am. Chem. Soc. 2003, 125, 8072-8073. (37) Barbieri, R.; Bertini, I.; Cavallaro, G.; Lee, Y. M.; Luchinat, C.; Rosato, A. J. Am. Chem. Soc. 2002, 124, 5581-5587. (38) Tjandra, N.; Tate, S.; Ono, A.; Kainosho, M.; Bax, A. J. Am. Chem. Soc. 2000, 122, 6190-6200. (39) Hansen, M. R.; Rance, M.; Pardi, A. J. Am. Chem. Soc. 1998, 120, 11210-11211. (40) Tian, F.; Al Hashimi, H. M.; Craighead, J. L.; Prestegard, J. H. J. Am. Chem. Soc. 2001, 123, 485-492. (41) Staaf, M.; Ho¨o¨g, C.; Stevensson, B.; Maliniak, A.; Widmalm, G. Biochemistry 2001, 40, 3623-3628. (42) Landersjo¨, C.; Jansson, J. L. M.; Maliniak, A.; Widmalm, G. J. Phys. Chem. B 2005, 109, 17320-17326. (43) Landersjo¨, C.; Stevensson, B.; Eklund, R.; O ¨ stervall, J.; So¨derman, P.; Widmalm, G.; Maliniak, A. J. Biomol. NMR 2006, 35, 89-101. (44) Emsley, J. W.; Luckhurst, G. R.; Stockley, C. P. Proc. R. Soc. London, Ser. A 1982, 381, 117-138. (45) Catalano, D.; Di Bari, L.; Veracini, C. A.; Shilstone, G. N.; Zannoni, C. J. Chem. Phys. 1991, 94, 3928-3935. (46) Sandstro¨m, D.; Levitt, M. H. J. Am. Chem. Soc. 1996, 118, 69666974. (47) Emsley, J. W. In Encyclopedia of NMR; Grant, D. M., Harris, D. M., Eds.; Wiley: New York, 1996; Chapter 2781. (48) Stevensson, B.; Komolkin, A. V.; Sandstro¨m, D.; Maliniak, A. J. Chem. Phys. 2001, 114, 2332-2339. (49) Alejandre, J.; Emsley, J. W.; Tildesley, D. J.; Carlson, P. J. Chem. Phys. 1994, 101, 7027-7036. (50) Palke, W. E.; Catalano, D.; Celebre, G.; Emsley, J. W. J. Chem. Phys. 1996, 105, 7026-7033. (51) Berardi, R.; Spinozzi, F.; Zannoni, C. Chem. Phys. Lett. 1996, 260, 633-638. (52) Berardi, R.; Spinozzi, F.; Zannoni, C. J. Chem. Phys. 1998, 109, 3742-3759. (53) Catalano, D.; Emsley, J. W.; La Penna, G.; Veracini, C. A. J. Chem. Phys. 1996, 105, 10595-10605. (54) Cinacchi, G.; Longeri, M.; Veracini, C. A. Phys. Chem. Chem. Phys. 2002, 4, 5582-5589. (55) Stevensson, B.; Landersjo¨, C.; Widmalm, G.; Maliniak, A. J. Am. Chem. Soc. 2002, 124, 5946-5947. (56) Thaning, J.; Stevensson, B.; Maliniak, A. J. Chem. Phys. 2005, 123, 044507-1-044507-6. (57) Stevensson, B.; Sandstro¨m, D.; Maliniak, A. J. Chem. Phys. 2003, 119, 2738-2746. (58) Celebre, G.; Cinacchi, G. J. Chem. Phys. 2006, 124, 17601-117601-2. (59) Tieleman, D. P.; Berendsen, H. J. C. J. Chem. Phys. 1996, 105, 4871-4880. (60) Berger, O.; Edholm, O.; Jahnig, F. Biophys. J. 1997, 72, 20022013. (61) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306-317. (62) Hoover, W. G. Phys. ReV. A 1986, 34, 2499-2500. (63) Parinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196-1199. (64) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 1008910092. (65) Dong, R. Y. Nuclear Magnetic Resonance of Liquid Crystals; Springer: New York, 1994. (66) Brink, D. M., Satchler, G. R. Angular Momentum; Clarendon Press: Oxford, 1993. (67) MATLAB; The MathWorks: Natick, MA, 1999. (68) Pearson, R. H.; Pascher, I. Nature 1979, 281, 499-501. (69) Gally, H. U.; Pluschke, G.; Overath, P.; Seelig, J. Biochemistry 1981, 20, 1826-1831. (70) Celebre, G.; De Luca, G.; Emsley, J. W.; Foord, E. K.; Longeri, M.; Lucchesini, F.; Pileio, G. J. Chem. Phys. 2003, 118, 6417-6426. (71) Emsley, J. W.; Longeri, M.; Merlet, D.; Pileio, G.; Suryaprakash, N. J. Magn. Reson. 2006, 180, 245-255. (72) Dvinskikh, S. V.; Du¨rr, U. H. N.; Yamamoto, K.; Ramamoorthy, A. J. Am. Chem. Soc. 2007, 129, 794-802.
