Crystal structure refinement using analytical derivatives of the energy

Multibilayer structure of dimyristoylphosphatidylcholine dihydrate as determined by energy minimization. Garret Vanderkooi. Biochemistry 1991 30 (44),...
0 downloads 0 Views 920KB Size
Download PDF
4366

J . Phys. Chem. 1990, 94. 4366-4372

Crystal Structure Refinement Using Analytical Derivatives of the Energy Function. Application to 1,2-Dilauroyl-~~-phosphatidylethanolamine Garret Vanderkooi Department of Chemistry, Northern Illinois University, DeKalb, Illinois 601 15 (Received: August 25, 1989) Methods are presented for refining crystal structures of organic molecules by minimization of the energy. The crystal energy is expressed as a sum of intramolecular and intermolecular terms and includes nonbonded, hydrogen-bonded, electrostatic, and torsional contributions. Equations are given for computing the first derivatives of the crystal energy with respect to intramolecular bond rotations, rigid-body motions, and the lattice constants. These methods were used to carry out an energy refinement of the crystal structure of 1,2-dilauroyl-~~-phosphatidylethanolamine:acetic acid (DLPE:HAc). Minimization was simultaneously performed with respect to all intramolecular bond rotations (34 in DLPE and 2 in the cocrystallized molecule of acetic acid) and all rigid-body parameters (6 each for DLPE and acetic acid). The crystal energy changed from -1 29.7 to -142.5 kcal/mol and the derivatives decreased to values close to zero as a result of minimization. Prior to energy refinement of the structure, minimization with respect to the lattice constants caused an increase of 57 A3 in the unit cell volume, but after refinement the energy-minimized lattice constants were within the error limits of the experimental values. An analysis of the total crystal energy is given in terms of the monolayer, intrabilayer, and interbilayer contributions to the total energy.

There is widespread interest in the properties of lipid bilayers and membranes, but detailed structural information has been difficult to obtain for these systems since lipids do not crystallize very readily or very well. The first successful three-dimensional X-ray structural analysis of a phospholipid was published in 1974 by Hitchcock et ai.,' for 1,2-dilauroyl-~~-phosphatidylethanolamine:acetic acid (DLPE:HAc), and refined coordinates were published in 1977.* Since that time, several other lipids or lipid-like molecules have been crystallized and their structures determined by X-ray diffraction, including dimyristoyldimyristoylpho~phatidate,~dimyristoylpho~phatidylcholine,~ ph~sphatidylglycerol,~ cerebroside,6 and several lyso compounds and other glycerol esters. In all of these crystals there is a problem with disorder, especially in the aliphatic-chain regions. In spite of this, these structures do provide a substantial experimental basis for understanding the conformations and intermolecular interactions of phospholipids in the bilayer arrangement. Conformational energy minimization methods have been highly developed by several laboratories for studies on macromolecular structures of proteins and nucleic acids, but much less effort has been devoted to energetic calculations of lipid bilayer conformations. The work that has been done can be briefly summarized. The earliest lipid calculations were published in 1973 by Vanderkooi7 and by McAlister et aL8 and dealt with the problem of the sterically permitted and energetically preferred conformations of isolated molecules. The Pullman group,9 as well as Gupta et al.,IOcarried out extensive quantum mechanical calculations on lipids in the mid- 1970s, and Frischleder and co-workers,l'Sl2 Kreissler et aI.,l3 and Brasseur et aI.l4each carried out calculations at various levels of detail on small clusters of lipid molecules. The ( I ) Hitchcock, P. B.; Mason, R.; Thomas, K. M.; Shipley, G. G. Proc. Natl. Acad. Sci. U . S . A . 1974, 71, 3036. ( 2 ) Elder, M.; Hitchcock, P.; Mason, R.; Shipley, G. G. Proc. R . Soc. London 1977. A354, 157. (3) Pearson, R. H.; Pascher, I. Nature 1979, 281, 499. (4) Harlos. K.; Eibl, H.: Pascher, I.; Sundell, S. Chem. Phys. Lipids 1984, 34, 115. ( 5 ) Pascher, I.;Sundell, S.; Harlos, K.; Eibl, H. Biochim. Biophys. Acta 1987, 896, 77. (6) Pascher, 1.; Sundell, S. Chem. Phys. Lipids 1977, 20, 175. (7) Vanderkooi, G. Chem. Phys. Lipids 1973, 11, 148. (8) McAlister, J.; Yathindra, N.; Sundaralingam, M. Biochemistry 1973, 12. 1189,

(9) Pullman, B.; Saran, A. lnr. J . Quantum Chem.: Quantum Biol. Symp. 1975. No. 2. 7 1 ,

(IO) Gupta, S. P.; Govil, G.; Mishra, R. K. J . Theor. Biol. 1975, 51, 13. ( I 1 ) Frischleder, H. Chem. Phys. Lipids 1978, 22, 343. ( 1 2) Frischleder, H.; Lochmann, R. Int. J . Quantum Chem. 1979, 16, 203. ( I 3) Kreissler. M.: Lemaire. B.; Bothorel, P. Biochim. Biophys. Acta 1983, 735, 12. (14) Brasseur, R.; Goormaghtigh, E.; Ruysschaert, J. M. Biochem. Biophys. Res. Commun. 1981, 103. 301.

Brasseur group has also published a series of papers on the preferred arrangements of lipid molecules surrounding amphipathic drug molecules or membrane peptides (see, e.g., ref 15). For fluid systems such as micelles or lipid bilayers at temperatures above the thermal phase transition, the energy minimization approach is not appropriate on account of the many equivalent energy states and the consequent importance of the entropy contribution to the free energy. For such systems, molecular dynamics methods may be used, and a report has recently appeared in which a molecular dynamics simulation was carried out on a hydrated micelle of lysophosphatidylethanolamine.I6 The approach taken in this work is to use the known crystal structures of lipids as starting points for energy calculations on lipids in the gel state. In this paper, general methods are presented for refinement of crystal structures by energy minimization, and these methods are used to refine the X-ray structure of DLPE:HAc2 The earlier literature on the use of energy calculations in crystal structure analysis has been extensively reviewed by Timofeeva et aI.,l7 and the mathematical methods developed by Williams's-20 have been It is clear from these reviews that while a large number of papers have appeared in the general area of energy calculations on crystal structures, the majority deal with structures of relatively small molecules with few internal degrees of freedom. A notable exception is the recent paper by Glasser and ScheragaZZin which the crystal packing of Leu-enkephalin was studied, involving a total of 332 atoms in the asymmetric unit and 176 independent variables. At least three computer programs have previously been described for calculating and minimizing the energy of crystals. Glasser and Scheraga** adapted the WMIN program of BusingZ3 for use in their Leu-enkephalin studies. This program uses the accelerated convergence procedure19 and employs minimization methods that do not require the calculation of analytical gradients. Williams has released a program called PCK8324 in which the first and second derivatives of his energy function are computed with ( I 5) Brasseur, R.; Cabiaux, V.; Killian, J. A,; de Kruijff, B.; Ruysschaert. J . M. Biochim. Biophys. Acta 1986, 855, 317. (16) Wendoloski, J. J.; Kimatian, S. J.; Schutt, C. E.; Salemme, F. R. Science 1989, 243, 636. (17) Timofeeva, T. V.; Chernikova, N. Yu.; Zorkii, P. M. Russ. Chem. Reo. 1980, 49, 509. (18) Williams, D. E. Acta Crystallogr. 1969, ,425, 464. (19) Williams, D. E. Acta Crystallogr. 1971, A27, 452. (20) Williams, D. E. Acta Crystallogr. 1972, A28, 629. (21) Williams, D. E. Top. Curr. Phys. 1981, 26, 3. (22) Glasser, L.; Scheraga, H. A. J . Mol. Biol. 1988, 199, 513. (23) Busing, W. R. ORNL-5747; Oak Ridge National Laboratory: Oak Ridge, TN, 1981. (24) Williams, D. E. QCPE Bull. 1984, 4 (3), 82 (QCPE Prog. No. 481).

0 1990 American Chemical Society

Phosphatidylethanolamine Crystal Refinement

The Journal of Physical Chemistry, Vol. 94, No. 10, 1990 4361

N ,a6

Figure 1. Diagram of DLPE and acetic acid showing the atom and bond numbering system empl~yed.~

respect to the lattice constants, molecular rotations and translations, and internal rotations, and uses either the Newton-Raphson method or the steepest descent method to minimize the energy. A program has also been described by Pietila and R a ~ m u s s e n ~ ~ that has many of the characteristics of the Williams program. The program used in the present work was developed with the explicit application in mind of being used for calculations on crystal structures of molecules with many internal degrees of freedom and with the possibility of several nonsymmetry-related molecules within the asymmetric unit. Analytical first derivatives of the energy function are used in the minimization procedure, but the accelerated convergence method has not been implemented. The efficient derivative factorization method described by Levitt26is used, with the result that the time required for a complete gradient plus energy calculation is only about twice that of an energy calculation alone.

Computational Methods Geometry. Molecular Structure. DLPE cocrystallizes with acetic acid in the P 2 , / c space group, with 4 DLPE and 4 acetic acid molecules per unit cell and 97 atoms in the asymmetric unit. An inversion center relates pairs of the D and L molecules of the racemic mixture of DLPE. The lipid molecules are in a stacked multibilayer arrangement, and the acetic acid molecules are associated with the polar head groups of DLPE in the spaces between the bilayers. The DLPE structure is shown in Figure 1. The atom and bond labeling scheme follows that used in the original crystal and correponds to the nomenclature proposed by S~ndaralingam.~’ The fractional cell coordinates given by Elder et aL2 for the non-hydrogen atoms provided the starting point for the calculations. All hydrogen atoms are explicitly included, with their coordinates being generated from standard bond lengths and angles.28 Bond lengths and angles were not permitted to vary during the energy-refinement procedures, but torsional rotations about all bonds were permitted. Symmetry Operations. Fractional cell coordinates are transformed to Cartesian coordinates by using the D matrix:’* ri = Dx, (1) ri denotes the position vector of the ith atom in Cartesian space and xi the same atom in unit cell space. Symmetry-related molecules are generated in unit cell space by the equation

si is the jth point group symmetry operation, and t k is the kth cell translation operation. An analogous equation applies in Cartesian space: rijk = DsiD-lrj + Dt, (3) The molecular geometry must be maintained constant while the lattice constants are varied during energy minimization. This may be done in two ways: i n the first, which was used in the (25) Pietila, L. 0.;Rasmussen, K . J . Compuf. Chem. 1984, 5 , 252. (26) Levitt, M. J. Mol. Biol. 1983, 170, 723. (27) Sundaralingam, M. Ann. N . Y . Acad. Sci. 1972, 195, 324. (28) Vanderkooi, G . J. Phys. Chem. 1983, 87, 5121.

present work, the Cartesian coordinates of all of the atoms in the reference asymmetric unit are held fixed while the lattice changes, while in the second method, the unit cell coordinates of a particular atom are held fixed, while the unit cell coordinates of the remaining atoms change so as to maintain the same relative geometry in Cartesian space. The latter method was used by Williamsz1and should be employed if an atom lies at a special position, but it gives problems in energy minimization if the rigid-body translations and the lattice constants are treated simultaneously as variables. Changes in the lattice constants result in molecular translations in Cartesian space, and hence the lattice constants and rigid-body translations are not linearly independent. This problem is avoided, without any loss of generality, by using the first method in which the Cartesian coordinates of the reference asymmetric unit are not altered by lattice constant changes. By this method, the lattice variables affect the translation vectors between symmetry-related groups, while the rigid-body variables determine the positions and orientations of the one or more molecules within the asymmetric unit. Rigid-Body Geometry. Specification of the spatial orientation and position of a rigid body requires six parameters: three translational and three rotational. A molecule-based coordinate system was defined by using three atoms whose coordinates do not change as a result of any bond rotations. One of these atoms is the “center atom”, which defines the center about which rigid-body rotation will occur. For DLPE, C2 was placed at the origin; the C2-C3 bond defined the z axis and 021, together with the C2 and C3 atoms, defined the x z plane. A right-handed Cartesian coordinate system was constructed on the basis of these atoms. For acetic acid, the origin was on C41 (carbonyl carbon), the z axis was along the C41-C42 bond, and 0 4 2 completed the coordinate description. Rigid-body rotations were defined in terms of the moleculebased coordinate systems using the Eulerian angles a , 6, and y, as defined by Margenau and Murphy:29 a corresponds to rotation about the z axis, 6about the x axis, and y about z’, which is the orientation of the z axis after the a and 6 rotations have been applied. These rotations must be transformed so as to be expressed in the fixed coordinate system which coincides with the orthogonalized crystal lattice. If there is more than one molecule in the crystallographic asymmetric unit, the molecules are numbered with the principal molecule being first. Rigid-body motions of the second and higher numbered molecules occur in local coordinate systems based on those molecules, but these are in turn related to the local coordinate system based on the first molecule. These other molecules are therefore carried along with the first molecule when the first molecule undergoes rotations or translations. In the present case, DLPE is identified as the first molecule and acetic acid is associated with it as the second molecule. Let El be the Eulerian rotation matrix for molecule 1, and let M I be the matrix that brings the coordinate system based on molecule 1 into alignment with the fixed coordinate system. A similarity transform (MIEIMI-I)can then be constructed that will (29) Margenau, H.; Murphy, G. M . The Mathematics of Physics and Chemistry, 2nd ed.; Van Nostrand: Princeton, NJ, 1956; pp 286-289.

4368 The Journal of Physical Chemistry, Vol. 94, No. 10, 1990 have the effect of carrying out the El rotation in the fixed coordinate system. Translations of molecule 1 (TI)are taken parallel to the fixed axes and may be carried out in the same operation as the rotation: rill = MlE,Ml-l(rj1- rOl)+ rol + T,

(4)

roI are the coordinates of the center for rigid-body rotation, and ril and rill are the coordinates of atom i before and after the rotation and translation, respectively, expressed in the fixed coordinate system. When more than one molecule is present in the asymmetric unit, the molecule-based rotations and translations of the second and higher numbered molecules are transformed back first into the molecule-based coordinate system of the first molecule and second to the fixed system (In practice, these operations are combined into a single step.):

+ roz + T, rj2” = MlE,Ml-1(ri2’- rol) + rol + Tl ri2’ = M2EZM2-l(ri2 - r,,)

(5)

M2 transforms the coordinate system based on molecule 2 into the system based on molecule 1, r i i are the rotated and translated coordinates of molecule 2 expressed in the molecule 1 system, and r,2” are the transformed coordinates of molecule 2 in the fixed system. Bond Rotations. The effects of bond rotations are computed by setting up a local coordinate system on the bonded atoms which define the dihedral angle and using a similarity transform to relate the effects of the bond rotation to the fixed system: rl’ = MRM-l(r, - r,)

+ ro

(7)

M is a matrix that transforms the bond-based coordinate system into the fixed system, R is a rotation matrix for rotating about the local axis that parallels the bond, ri and r,’ are the atomic coordinates of atoms affected by the bond rotation, and ro are the coordinates of the first atom of the bond about which rotation occurs. The standard dihedral angle convention is employed; Le., when sighting down the bond in question, Oo corresponds to the cis orientation of the first and fourth defining atoms and clockwise rotation of the far atom relative to the near atom is the positive sense. Energy Calculations. Energy Functions and Parameters. The empirical energy was computed as a sum of nonbonded, electrostatic, hydrogen-bonded, and intrinsic torsional energy terms with the parameters and equations of Momany et aL30 The nonbonded energy is calculated as the sum of all intramolecular and intermolecular atom-atom pairwise interactions, using a Lennard-Jones “6-1 2” potential function. Intramolecular interactions between atoms separated by more than three bonds are treated in the same manner as intermolecular interactions, but the strength of the repulsive part of the nonbonded function is decreased by half for “1-4” interactions between atoms separated by three bonds. Hydrogen bonds were represented with the “10-1 2” hydrogen bond function developed by McGuire et al.31 In this approach, the 6-1 2 nonbonded interaction between the hydrogen-bonding H.-0 or He-N atoms is replaced with a special 10-12 hydrogen bond function, but all other nonbonded and electrostatic energy terms between the hydrogen-bonded groups are evaluated in the usual manner. No explicit angular dependence is employed. The crystal-refined hydrogen bond parameters of Vanderkooi28were used for interactions involving the phosphate group, while those of Momany et aL30 with revisions as noted by Sippl et al.,32were used for other types of hydrogen bonds. The electrostatic energy was calculated in the monopole approximation using the Coulomb equation, by assigning partial (30) Momany, F. A,; Carruthers, L. M.; McGuire, R. F.; Scheraga, H. A. 1595. (31) McGuire, R . F.; Momany, F. A.; Scheraga, H. A . J. Phys. Chem. 1912, 76, 315. (32) Sippl, M. J.; NBmethy, G.; Scheraga, H. A. J. Phys. Chem. 1984,88, 623 1 ,

Vanderkooi charges to all of the atoms in the polar parts of the molecules. An overlap-normalized C N D 0 / 2 calculation was carried out to determine the partial charges. Atoms beyond C24 and C34 in the acyl chains were given zero charges. The dielectric constant was held at a fixed value of 2.0, in agreement with the value employed in the energy parameter refinement calculation^.^^^^^ An intrinsic torsional potential was computed for each bond about which rotation could occur, using the single-term cosine function and parameters as given by Momany et al.33for all except the 0-P bonds. In the latter case, a two-term function which included 2-fold and 3-fold cosine terms was used, with barrier heights as given by Weiner et al.,34in order to reproduce the known tendency of phosphate ester bonds to assume a gauche configuration. Computational Procedures. A system involving electrically neutral “sphere groups” was adopted to permit the use of cutoff distances for the intermolecular electrostatic and nonbonded energies without creating spurious unbalanced electrical charges. Each molecule is divided into one or more clusters of atoms for which the sum of the partial charges can be made to equal zero by minor numerical adjustments of the computed charges. A dummy atom is placed at the approximate center of each cluster, or “sphere group”, and the radius of the group is taken as the distance from the dummy atom to the farthest real atom of the group. The dummy atom thereafter transforms along with the real atoms. Initial distance checks between these clusters of atoms are done between their respective dummy atoms. If the dummy atom-dummy atom distance is less than the sum of the respective radii of the two sphere groups involved, plus the specified cutoff distance, then electrostatic calculations are carried out for all pairs of atoms between the two groups, but nonbonded calculations are only done between pairs whose actual interatomic distance is less than the specified cutoff distance. Provision was also made that during any given iteration or cycle of energy minimization, the same set of sphere-group interactions was maintained throughout the cycle. In this way the abrupt small steps in the energy were avoided which would otherwise occur as groups moved into or out of the region of computation, and which cause havoc with derivative-based minimization methods. This procedure is an alternative solution to the cutoff problem from that employed by Levitt,26who used a modified potential function that asymptotically approaches zero more rapidly than the unmodified potential function. For interactions between the reference asymmetric unit and its symmetry-related environment, an explicit list of space group symmetry operations is developed that will generate all molecules out to the cutoff distance. Symmetry operations are applied in a trial manner to the sphere group dummy atoms, and the distances between the symmetry-related dummy atoms thus generated and those in the reference asymmetric unit are calculated. If any of the distances between the dummy atoms in the reference unit and the dummy atoms generated by a trial symmetry operation are less than the sum of the respective sphere group radii plus the cutoff distance, then the trial symmetry operation is provisionally added to the list of operations to be subsequently used in energy calculations. Before actually adding a new operation to the symmetry list, however, it is tested to see if it is the exact inverse of an operation already in the list. If the answer to this test is positive, then the new operation is not added to the list, but its inverse already in the list is marked so that energies calculated by it will be doubled. The total energy, expressed as energy per mole of asymmetric units, equals the sum of the intramolecular energy, plus the intermolecular pairwise energies between molecules within the asymmetric unit, plus half of the sum of all interactions between the reference asymmetric unit and its surrounding symmetryrelated environment.

J . Phys. Chem. 1974, 78,

(33) Momany, F. A,; McGuire, R. F.; Burgess, A. W.; Scheraga, H . A. J. Phys. Chem. 1975, 79, 2361. (34) Weiner, S. J.; Kollman, P. A.; Case, D. A,; Singh, U.C.; Ghio, C.; Alagona, G.: Profeta. S.; Weiner. P. J. Am. Chem. SOC.1984, 106, 765.

The Journal of Physical Chemistry, Vol. 94, No. 10, 1990 4369

Phosphatidylethanolamine Crystal Refinement

Energy Minimization. Analytical derivatives were computed simultaneously with the evaluation of the energy. Energy minimization was carried out by using the DUMING and DUMIDH programs in the IMSL Math Library, version 1.0 (1987). The DUMlNG program uses a quasi-Newton method,3s and DUMIDH35*36 uses a modified Newton method; both employ analytical first derivatives and numerically generate the Hessian of second derivatives. Differentiation. The total energy function ( E ) may conveniently be written as a sum of intrinsic torsional potentials (I) and the combined electrostatic, nonbonded, and hydrogen-bonded terms which depend upon the pairwise interatomic distances (V): E=I+V

(8)

The partial derivative of E with respect to an arbitrary geometrical variable 0 is then given as aE

ae

-

a i + -av ae ae

Rigid-Body Translations. For rigid-body translations (0,) which are defined in the fixed coordinate system, ari/ae, simply equals the unit vectors along the x, y , or z axis (ex, e,,, ez). Translations of molecule 1 are defined in the fixed system, but for higher numbered molecules, the local molecule-based coordinate system rotates along with molecule 1, and hence the derivative vectors must also be transformed back to the fixed system using the transformation matrix for molecule 1:

(9)

The intrinsic torsional potentials are dependent only on the dihedral angles of the bond rotations, and hence all80 may be evaluated simply and directly. Evaluation of aV/aO, on the other hand, requires use of the chain rule for differentiation. Let wijk be the pairwise potential energy function for the interaction between atom i on the reference molecule and atom j on a molecule generated by the kth symmetry operation. (Let the identity symmetry operation correspond to k = 1; this term will then include intramolecular and inter-rigid-body interactions within the reference asymmetric unit.) Differentiation of the pairwise potential function with respect to 0 can be written as

Rigid-Body Rotations. The Eulerian angle rotation matrix must be differentiated term by term with respect to each of the Eulerian angles in order to get the corresponding set of derivative matrices (aE/aO,, where Be are the Eulerian angles and E is the Eulerian rotation matrix). These derivative matrices are defined in the local, molecule-based coordinate systems and must be transformed into the fixed coordinate system. Partial derivatives with respect to the Eulerian angles of molecule 1 are obtained as follows:

a rill a El Mi-l(ril - rol) - - MI aeei

see,

The symbols have the same meanings as in eq 4, and rol and ril correspond to the unrotated coordinates in the fixed system. Partial derivatives with respect to the Eulerian angles of molecule 2 (ee2) are obtained in the fixed system by using the combined effects of the transformation given in eqs 5 and 6: ari2/’

--

doe2 awijk/alruk! is evaluated directly by differentiating the potential function with respect to the distance. alrijkl/a(rjk-ri) yields the unit vector upon differentiation, eijk,the elements of which are the direction cosines of the vector from atom i to atom j k . a(rjk-ri)/aO may be written out as (arjk/ae - ari/ae; the two terms correspond to the effects of 6’ on the atoms of the symmetry-related molecules and on the atoms in the reference asymmetric unit, respectively. arjk/&9is evaluated by applying the rotational part of the kth symmetry operation to the partial derivative for the corresponding jth atom of the reference molecule:

Lattice Constants. The D transformation matrix must be differentiated term by term in order to obtain the partial derivatives with respect to the lattice constants. In this case, the value of the combined expression shown in eq 10, Le., a(rjk - ri)/aOl, is obtained directly when 0, is any one of the six lattice constants. If, as is done here, the procedure is used in which the Cartesian coordinates of the atoms in the reference asymmetric unit are held constant while the lattice is varied, then the derivatives for all pairs of atoms between the reference asymmetric unit and symmetryrelated asymmetric units are simply equal to the derivative of the translation part of the symmetry operation:

Equation 10 can be factored, as shown by Levitt:26 The terms in each of the six aD/aOi matrices have been given by Williams.20

The inner sum in the first term (over i ) is the sum over all atoms i in the reference asymmetric unit interacting with a t o m j in the kth symmetry-related asymmetric unit, and likewise the inner sum over j in the second term is over all atoms j in the kth symmetry-related asymmetric unit interacting with atom i in the reference unit. These two inner sums can be accumulated at the same time as the energy summations are being performed. Evaluation of ari/a6: Bond Rotations. Partial derivatives for the atoms affected by bond rotations (6,) are obtained by taking the cross product of the unit vector lying along the bond (Le., eAB, where the bond extends from atom A to B) and the vector from atom A to atom i:36 (35) Dennis,J. E., Jr.; Schnabel, R. B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations; Prentice-Hall: Englewood Cliffs, NJ, 1983; Appendix A. ( 3 6 ) Gay, D. M. ACM Trans. Math. Soft. 1983, 9, 503-524. (37) Katz, H.; Walter, R.; Somorjai, R. L. Comput. Chem. 1979, 3, 25.

Results Structure Refinement. The conformation and packing of the DLPE and acetic acid molecules in the DLPE:HAc crystal were refined by using the derivative-based minimization methods described in the preceding sections. The lattice constants were held fixed at their crystallographic values during the structure refinement, but all other variables (bond rotations and rigid-body motions) were permitted to change. A preliminary calculation step was carried out to provide favorable orientations for the terminal CH3, NH3, and COH groups whose hydrogen coordinates were generated using standard geometry; this yielded what is later referred to as the “initial structure”. This structure was refined in several steps. In the earlier steps, a cutoff distance of 7 8, was employed, and minimization was carried out with respect to various subsets of the total list of variables. In the final stages of refinement, a 12-A cutoff distance was used, and minimization was carried out with respect to all variables except the lattice constants (34 bond rotations in DLPE, 2 bond rotations in acetic acid, and 6 rigid body variables for each of the molecules, giving a total of 48 variables.)

The Journal of Physical Chemistry, Vol. 94, No. 10, 1990

4370

TABLE I: Rigid-Body Parameters and Torsional Angles for Bond Rotations of the Initial and Energy-Refined Crystal Structures initial energy-refined structure structure change Rigid-Body Parameters DLPE translation, A X 0.000 -0.046 -0.046 0.036 Y 0.000 0.036 -0.019 Z 0.000 -0.019 rotation, deg Ly 0.000 -68.058 -68.058 0.210 3 0.000 0.2 10 66.864 -r 0.000 66.864 translation, A X

L' Z

rotation, deg a

P Y

Acetic Acid

0.000 0.000 0.000

0.244 -0.262 -0.043

0.244 -0.262 -0.043

0.000 0.000 0.000

14.284 9.019 -8.185

14.284 9.019 -8.185

Bond Rotations: Torsional Angles," deg

DLPE 307.905 187.657 205.649 58.489 66.180 106.347 67.008 96.61 I 179.315 241.475 64.806 18 1.902 185.182 175.905 186.595 178.205 180.61 1 182.510 180.477 181.744 172.679 178.907 188.570 187.359 174.043 181.289 183.043 182.056 182.198 172.597 160.887

302.812 182.151 203.651 61.51 1 58.049 99.232 69.538 170.836 97.751 176.743 239.750 66.002 174.852 18 1.443 174.701 18 1.330 177.513 18 1.433 179.349 180.781 178.003 172.293 177.475 179.972 178.829 179.640 179.728 179.177 181.953 177.022 187.426 178.357 181.782 177.91 1

-5.093 -5.506 -1.999 3.022 -8.131 -7.1 I5 2.529 1.140 -2.572 -1.725 1.197 -7.050 -3.740 -1.204 -5.264 -0.692 0.822 -3.161 0.304 -9.451 4.796 1.064 -9.741 -7.719 5.686 -2.112 -1.090 -5.034 5.227 5.760 20.895

Acetic Acid 81 62

175.544 177.276

"Initial values are not given for angles whose definitions involve generated hydrogen atom coordinates.

Vanderkooi The dihedral angles for the bond rotations and the rigid-body motion parameters are given in Table I for the initial structure and for the final energy-refined structure. The translations and rotations for the initial structure are all given as zero since the initial structure was used to define the zero positions for rigid body motions. The x,y,z translation values given in Table I corres ond to rigid-body translations of 0.06 8, for DLPE and of 0.36 for acetic acid relative to DLPE. The numerical values of the a and y Eulerian rotations for DLPE are deceptively large, but since they are nearly equal in magnitude but opposite in sign (-68.1' and 66.9"), they will largely cancel each other, since the p rotation is small. A further way in which the degree of change in the structures can be evaluated is in terms of the mean atomic displacement, which is defined as the mean of the distances between the initial and energy-refined positions of the non-hydrogen atoms, and includes the effects of both bond rotations and rigid-body movements. For DLPE, the mean displacement is 0.21 A, and the mean displacement for the acetic acid is 0.37 A. The four carbon atoms at the end of the y-acyl chain of DLPE have a very large displacement, with the largest value being 1.17 A for C3 12, and the mean for these four atoms is 0.71 A. If these atoms are omitted, the remaining DLPE atoms show a mean displacement of only 0.16 A. It is noteworthy that these four atoms which have the largest displacement also have the largest estimated standard deviation in the experimental crystal structure,* which indicates that the crystals employed were disordered to some extent in the central region of the bilayer. Most of the torsional angles in the linear portions of the p- and y-acyl chains of the refined structure are appreciably closer to 180' than they are in the initial structure. The largest change in angle occurred in 712, which changed from 160.9' to 181.8'. This angle controls the position of C312, which, as noted above, also had the largest atomic displacement. Significant changes in the dihedral angles occurred in the polar head group region, but these were to some extent mutually compensating effects, since the corresponding atomic displacements are not large. The a 3 and a4 angles chan ed by -8.1" and -7.1 ", respectively, but N1 moved by only 0.13 from the initial to the final structure, and P moved by 0.17 A. The maximum atomic displacement in the a chain was at C11, which moved by 0.28

1

1

A.

Lattice Constants. Minimization of the energy with respect to the lattice constants was not used as a method to refine the structure, since the experimentally determined lattice values are accepted as reliable results of the X-ray structure determination. Rather, lattice constant minimization was used as a test to see how well the molecules could fit within the confines of the known lattice. For these calculations, the molecular structure (as defined by the bond rotations) was held constant while the energy was minimized with respect to the 12 rigid body parameters and the 4 lattice constants not specified by the space group. One might expect the lattice to expand upon energy minimization if there are atoms that are displaced from their correct positions. This is indeed what happens when the lattice constant minimization is carried out using the initial DLPE:HAc structure; the unit cell volume increased by 57 A3 (see Table 11). By contrast, when the energy-refined structure was used, minimization of the lattice constants caused an increase of only 16 A), which is within the experimental error. The intermolecular energy also changed by a much greater amount upon lattice constant minimization before

TABLE 11: Energy-Minimized Lattice Constants and Unit Cell Volumes for the Initial and Energy-Refined Structures lattice constants a, A b, A c, A 6 , deg cell vol, A) A(vol), A3 init struct, exptl lattice const 47.7 7.77 9.95 92.0 3685.5 0.0 init struct, energy-minimized lattice const' 48.025 7.678 10.161 92.71 3742.4 56.9 47.949 energy-refined struct, energy-minimized lattice const 7.769 9.951 3701.8 16.3 92.98

intermol energy, kcal/mol -84.5 -89.6 -90.8'

"The energy was minimized with respect to the a, b, c, and p lattice constants, and also with respect to the 12 rigid-body parameters, while holding the bond rotations fixed at the initial or energy-refined values. bThe energy of the energy-refined structure, calculated at the experimental lattice constants, was -90.6 kcal/mol (see Table V ) .

Phosphatidylethanolamine Crystal Refinement

The Journal of Physical Chemistry, Vol. 94, No. 10, 1990 4371

TABLE 111: Partial Derivatives of the Energy with Respect to the Rigid-Body Parameters and the Bond Rotationsa

intramoIb

initial structure intermolC totald Rigid-Body Derivatives

energy-refined structure intramolb intermol'

totald

DLPE

translation X

Y z

0.000 0.000 0.000

7.782 17.911 24.844

7.782 17.911 24.844

0.000 0.000 0.000

0.006 0.019 -0.006

0.006 0.019 -0.006

0.000 0.000 0.000

114.398 -25.849 114.398

114.398 -25.849 114.398

0.000 0.000 0.000

0.062 0.027 0.056

0.062 0.027 0.056

rotation cy

P Y

Acetic Acid translation X

Y z

0.000 0.000 0.000

-9.115 9.939 8.539

-9.115 9.939 8.539

0.000 0.000 0.000

0.010 0.02s 0.013

0.010 0.025 0.013

0.000 0.000 0.000

6.224 -8.419 6.224

6.224 -8.419 6.224

0.000 0.000 0.000

0.002 0.007 0.008

0.002 0.007 0.008

2.107 8.581 7.548 -4.94s -6.035 5.31 1 6.635 -3.651 -1.430 14.777 -3.845 6.531 -3.008 4.559 -3.270 3.613 -1.483 1.884 -0.651 0.942 -0.493 -2.450 4.519 -1.678 2.915 -0.901 1.986 -0.531 1.691 -1.546 3.192 -1.194 1.178 -1.018

-2.240 -8.562 -7.917 4.823 6.054 -5.234 -6.657 3.653 1.423 -14.766 3.878 -6.511 2.982 -4.553 3.259 -3.609 1.478 -1.889 0.646 -0.949 0.487 2.455 -4.523 1.686 -2.927 0.919 -2.007 0.542 -1.690 1.545 -3.202 1.192 -1.179 1.024

-0.133 0.019 -0.069 -0.122 0.019 0.076 -0.023 0.003 -0.007 0.01 1 0.033 0.019 -0.026 0.006 -0.01 1 0.004

-0.950 -0.571

0.966 0.571

0.015 -0.001

rotation a

P Y

61

83 01

n2 a3

a4 a5

4.220 -0.451 6.890 -4.630 -2.388 5.151 6.340

Bond Rotation Derivatives DLPE -41.796 -37.576 0.532 0.081 18.195 25.085 -0.747 -5.377 19.666 17.278 23.247 28.398 -2.572 3.768

a6 37.845 -7.250 36.999 2.136 3.780 3.311 1.181 5.100 1.641 1.886 1.641 1.280

-61.632 -104.338 -62.838 -48.460 46.487 -41.235 37.484 -43.648 26.098 -31.104 16.169 -16.768

-23.787 -1 11.589 -25.839 -46.324 50.267 -37.924 38.665 -38.548 27.738 -29.218 17.810 -15.488

710

8.625 -1.133 4.033 3.163 7.039 -2.992 2.806 1.008 1.284 1.537

Yll Y12

-2.231 -7.019

12.966 4.273 4.614 5.817 6.202 3.977 1.575 4.063 0.142 2.556 -0.553 5.353

21.590 3.141 8.646 8.979 13.241 0.985 4.380 5.071 1.426 4.094 -2.784 -1.667

PI

P2 P3 P4 PS P6 07 P8 P9 PI0 PI 1

PI2 PI3 Yl Y2 73 Y4 Y5 Y6 Y7 Y8 Y9

713

-0.005 -0.005 -0.005

-0.007 -0.006 0.005

-0.004 0.0u8 -0.012 0.018 -0.021 0.010

0.001 -0.001 -0.0I O -0.002 -0.001 0.005

Acetic Acid 81

62

a The derivatives are expressed in units of kcal/(mol.A) for translations and kcal/(mol.rad) for angular measure. Intramolecular components of the partial derivatives. Intermolecular components of the partial derivatives. dSum of the intramolecular and intermolecular components.

refinement than after refinement (from -84.5 to -89.6 kcal/mol before refinement versus -90.6 to -90.8 kcal/mol after refinement). Derioatiues. Table 111 gives the values of the partial derivatives of the energy with respect to all of the rigid-body and bond rotation variables, before and after the energy refinement. The intramolecular and intermolecular contributions were summed separately and are given along with the total values. The small total values for the derivatives after refinement is evidence that the net forces acting in each of these dimensions were reduced essentially to zero. The breakdown into the intramolecular and intermolecular

components, on the other hand, gives an indication of the extent to which the energy-minimized bond rotations are the result of a compensation effect between the intramolecular and intermolecular forces. The magnitudes of the derivative values do not indicate how much the dihedral angles would change if the intermolecular interactions were absent, however. Computer Times. The times required for energy calculations and for energy plus derivative calculations are given in Table IV, for calculations using cutoff distances of 7, 12, and 15 A. The total energy is also given for each cutoff distance. These calculations were carried out on the Amdahl 5870 Mainframe computer

4372

Vanderkooi

The Journal of Physical Chemistry, Vol. 94, No. 10, 1990

TABLE IV: Comouter Times as a Function of the Cutoff Distance cutoff energy energy and gradient total energy,' distance, k, calcn, s calcn (48 variables), s kcal/mol 7 .O 0.56 1.20 -135.2 12.0 1.48 2.84 -142.5 15.0 2.34 4.39 -143.4 ~

~

~

~~~~~~

"Energy as a function of cutoff distance, for the energy-refined structure.

TABLE V: Contributions to the Total Energy" initial energy-refined structure structure intramolecular energy, kcal/mol electrostatic -38.94 -38.88 -14.22 -17.45 nonbonded hydrogen bonded -0.16 -0.17 8.16 4.56 torsional

change 0.06 -3.23 -0.01 -3.60

total intramolecular

-45.16

-5 I .94

-6.78

intermolecular energy, kcal/mol electrostatic non bonded hydrogen bonded total intermolecular

-22.58 -49.29 -1 2.66 -84.53

-22.66 -54.40 -13.50

-0.08 -5.1 1 -0.84 -6.03

-129.69

-142.50

intramol

+ intermol, kcal/mol

TABLE VI: Monolayer, Intrabilayer, and Interbilayer Energy Components of the Energy-Refined Structure hydrogen electrostatic nonbonded bonded monolayer" -4.67 DLPE + HAcb -2 1.48 -45.03 -3.60 -19.02 -39.45 DLPE' -1.07 -2.46 -5.58 HAC~ intra bilayer' 0.0 DLPE + HAC* 0.0 -2.57 -2.57 0.0 DLPE' 0.0 0.0 HAC~ 0.0 0.0 interbilayerf DLPE + HAcb -1.18 -8.83 -6.80 0.47 DLPE' -0.41 0.0 -6.39 -1.65 -8.23 HAC~

totalg -71.18 -62.07 -9.1 1 -2.57 -2.57 0.0 -16.81 0.06 -16.87

Energy calculated from the reference unit to other molecules within the same monolayer (lateral interactions). *DLPE and acetic acid both included in calculation. cAcetic acid molecule omitted. dAcetic acid contribution to the energy; calculated as the difference between (DLPE HAC) and DLPE energies. 'Energy from the reference asymmetric unit to the molecules in the other half of the same bilayer (tail-tail interactions). /Energy from the reference asymmetric unit to the molecules in the succeeding bilayer (head-head interactions). gThe sum of the DPLE + HAC energies (monolayer, intrabilayer, and interbilayer) equals the intermolecular crystal energy (-90.56 kcal/mol). as given in Table V .

+

-90.56

-12.81

'These calculations were carried out using the experimental values of the lattice constants. A cutoff distance of 12 k, was employed.

at Northern Illinois University, operating under MVS/SP 1.3.6, using VS Fortran compiled under OPT(2).

Discussion Components of the Total Energy. The nonbonded, electrostatic, hydrogen-bonded, and torsional contributions to the intramolecular and intermolecular energies are listed in Table V. These are given for the initial and energy-refined structures, which enables one to observe which energy contributions are mainly responsible for the change in the total energy upon minimization. It can be seen that energy minimization causes decreases of a similar magnitude in both the intramolecular and intermolecular components (-6.7 and -6.0 kcal/mol, respectively). The intramolecular energy decrease is almost equally divided between the nonbonded and torsional energy terms, whereas for the intermolecular energy, the major decrease takes place in the nonbonded energy, with a smaller contribution from the hydrogen bond energy. In neither case is there a significant change in the electrostatic energy. This is consistent with the nature of the electrostatic energy function, which varies much more slowly with distance than does the nonbonded or hydrogen-bonded energy function. Bilayer Interactions. The DLPEHAc crystal consists of a stack of planar bilayers. The crystal energy is a sum of the interactions between the molecules in the same monolayer (lateral interactions), between the two monolayers of the same bilayer (tail-tail interactions), and between succeeding bilayers (head-head interactions). Each of these components of the intermolecular energy was determined for the energy-refined structure. A similar calculation was done while omitting the acetic acid molecules from the structure, so that the energy contributions which involve the acetic acid could be determined by difference. Table VI gives the results obtained. For lateral interactions between molecules within the same monolayer, the largest energy contribution is the nonbonded term, but the electrostatic and hydrogen-bonded terms also contribute a fair amount. Acetic acid has a relatively small effect on the lateral interaction energy.

Only the nonbonded energy term is nonzero for interactions between the monolayers of the same bilayer. This is reasonable, since the terminal portions of the fatty acyl chains are uncharged and the polar head groups are separated by a distance which is considerably larger than the cutoff distance used in the energy calculations. The tail-tail interaction energy (-2.57 kcal/mol) is the energy computed from the reference molecule in one leaflet of the bilayer to all atoms within the 12-A cutoff distance in the other leaflet of the same bilayer. For interactions between successive bilayers in the stacked multilayer (head-head interactions) the acetic acid molecules play a major part. The interaction between DLPE molecules in opposing bilayers calculated in the absence of acetic acid shows that the total energy is actually positive, arising from a positive electrostatic term and a smaller negative nonbonded term. When acetic acid is included, however, there are large negative nonbonded and hydrogen-bonded contributions from the acetic acid interactions, which serve to stabilize the stacked multilayer. The derivative-based minimization procedures used here have proven to be very effective in finding a multidimensional local minimum of the DLPE:HAc crystal structure which is close to the structure obtained by X-ray diffraction. The same methods should be useful in refining other crystal structures where experimental data are available but are not of sufficient quality to accurately define the structure. In the present work, the bond lengths and angles implied by the published coordinates were accepted without modification (except for the hydrogen coordinates which were generated by using standard geometry), but these could also be refined, if desired, using one of the available molecular mechanics programs in a preliminary isolated-moleculecalculation. The problem of refining crystal structures is of course much simpler than predicting crystal structures a priori, since the experimental data provide a starting point which is presumed to be near the global minimum, whereas in the a priori approach, the multiple local minimum problem must be dealt with.

Acknowledgment. This work was supported in part by Grant BRSG SO7 RR07 176 awarded by the Biomedical Support Grant Program, Division of Research Resources, National Institutes of Health.