J. Phys. Chem. 1995, 99, 5724-5731
5724
Orientation and Diffusion of a Drug Analogue in Biomembranes: Molecular Dynamics Simulations Howard E. Alper and Terry R. Stouch” Bristol-Myers Squibb Pharmaceutical Research Institute, P. 0. Box 4000, Princeton, New Jersey 08543-4000 Received: October 25, 1994@
Four nanoseconds of molecular dynamics simulation of a nifedipine analogue in a phospholipid bilayer was performed in order to gain insight into the interactions of a drug molecule with the membrane and the mechanism of drug diffusion through membranes. An all-atom representation of a fully hydrated DMPC bilayer in the L a phase was employed. Such simulations can provide detail inaccessible experimentally. It was found that the rate of diffusion did not vary with location in the bilayer, in contrast to that of the much smaller benzene, which was 3 times faster in the bilayer center than the interfacial region. Benzene diffusion is accelerated by its ability to jump between voids in the bilayer that are unavailable to nifedipine due to its larger size. Toward the water interface the nifedipine analogue experienced a diverse environment including interactions with lipid carbonyl and “headgroup” atoms as well as with water. On average one lipid carbonyl group and 4-6 water molecules interacted strongly with this molecule. Due to its motion, this environment varied during the simulation. Near the water interface, the preferred orientation of the drug analogue was significantly tilted relative to the bilayer normal. This differs from the orientation often assumed for the location of drug molecules in membranes during interpretation of X-ray scattering studies. This tilt evolved through concerted changes in side-chain torsions of the solute and the orientation of the entire solute which resulted in a significant increase in the hydrogen bonding of the solute with the lipids and (especially) water. This implies that an important contribution to the behavior of drugs in membranes is optimization of hydrogen bonding between the drug and its environment. Correlation between the movement of some of lipids and the solute suggests strong interactions or transient binding although not all neighboring lipids experienced this correlation. The nifedipine analogue did not introduce significant perturbations into the gross properties of the bilayer.
Introduction Many important biochemical processes take place in membranes. Of particular importance to drug bioavailability is the permeation of membranes by drugs. Once a molecule has partitioned in to a membrane, its rate of permeation is larely dependent on its rate of diffusion within the membrane. Although this is an important process, little is known concerning the atomic level details of the mechanism of permeation. Early theories related solute permeability to the partition coefficient of the molecule into a hydrophobic solvent.’ However, low molecular weight molecules were found to violate this rule and permeated faster than predicted by Overton’s r ~ l e . ~Walter .~ and Gutknecht studied a series of molecules of different size, weight, and complexity and verified Overton’s rule for solutes in the series with molecular weights between 50 and 300, while significant deviations from Overtons rule were found for solutes with molecular weight ‘50. For the latter it was found that permeability is inversely proportional to molecular size. However, this study employed only a limited number of different solutes, and in addition no attempt was made to determine the effects separately of molecular weight, size, and complexity on permeability. A proposed explanation for the deviations from Overton’s rule was the “soft-polymer’’ hypothesis.4~~ This states that small solutes diffuse by jumping between spontaneously arising, mobile voids, the size of the solute or larger, which are separated by low-energy barriers. Other aspects of solute diffusion through membranes have been probed by modem experimental methods. ESR and I3C spin relaxation experiments suggested that in dipalmitoyllethi-
* To whom correspondence should be addressed. Abstract published in Advance ACS Abstracts, March 1, 1995.
0022-365419512099-5724$09.00/0
cin membranes the translational motion of nitroxide solutes of the size of benzene is faster in the bilayer center than near the hydrocarbodwater interface? On this basis, permeation through the headgroup region was proposed as the rate-limiting step in transport of small molecules. NMR studies’ indicated that partitioning of benzene into DMPC and DPPC bilayer was dependent on the surface density of the lipid hydrocarbon chains. These two studies demonstrated that bilayers are interfacial in nature and cannot be treated as a homogeneous bulk phase. Although a number of simulations have been performed to understand the structure and dynamics of phospholipid bilayers,*-18 only a few have addressed the diffusion of drugs or other penetrants in lipid bilayers at the atomic l e ~ e l . ’ ~ , ~ ~ Recently we presented the results of simulations of benzene in a bilayer of dimyristoylphosphatidylcholine(DMPC) in the L a phase.l9 Simulations employing different numbers of benzene molecules were completed to determine the effect of solute concentration and location on the behavior of the benzene. Calculated diffusion coefficients and rotational correlation times were found to be in excellent agreement with experiment. Although the benzene moved freely throughout the bilayer, it was found that the rate of diffusion was determined in part by the position of the benzene molecule within the bilayer, in agreement with the experiments mentioned above.6 This shows that the bilayer and even its hydrocarbon interior has an internal structure and so cannot be treated as a uniform bulk phase in accord with DeYoung and Dill’s hypothesis. It was also found that the mechanism of benzene diffusion was similar to the “hopping” mechanism for diffusion of gases through polymers21 and lends support to the “soft-polymer’’ hypothesis of Lieb and Stein. 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99,No. 15, 1995 5725
Orientation and Diffusion of a Drug Analogue Benzene is a fairly small, uncomplicated molecule. Its molecular weight of 78 places it at the bottom to center of the series studied by Walter and Gutknecht to determine the validity for Overton’s rule.3 The present study expands our previous work by investigating the effect of increased solute size, polarity, and hydrogen-bonding capability on the mechaqism and rate of solute diffusion. The solute we study here is an analogue of nifedipine, a 1,4-dihydropyridine calcium channel antagonist (DHPs). Of particular interest is a comparison of the mechanism of diffusion of this molecule to that of a smaller nonpolar solute such as benzene (which has one-fourth the mass and one-third the volume of the nifedipine analogue) in order to clarify the relationship between the solute composition (size, weight, and complexity) on the rate and mechanism of diffusion and on agreement with Overton’s rule. In particular, the present study focuses on the consequences of the location of this solute within biomembranes. It has been proposed that drugs diffuse to receptors by first penetrating the membrane to a certain location and then laterally diffusing to the receptor22 and so the preferred location of drugs in membranes might be important. However, considerable uncertainty exists in the literature concerning the location of DHPs in membranes. For example, X-ray studies of charged and uncharged DHPs suggest that they are located at the hydrocarbod water interface in the membrax~e.*~-~~ However, since the resolution of such studies is low, DHPs with equal validity could be placed either in the hydrocarbon region or at the hydrocarbod water interface which is only 10-15 A distant. Furthermore, N M R experiments suggest that uncharged DHPs are distributed uniformly throughout the bilayer.27 Even should these approaches determine the approximate location of solute, they haven’t sufficient precision to determine solute orientation. Yet the orientation of the drug, in part, determines its interactions with membrane components, and a knowledge of these interactions is necessary for an atomic-level understanding of drug penetration and diffusion, as well as action at membrane-bound receptors. In this study, we present two 2-ns simulations of the diffusion of the nifedipine analogue within an all-atom, fully hydrated phospholipid membrane (over 7100 atoms) in order to determine the effect of the solute environment on the diffusion of the drug, as well as its orientation and interaction with the bilayer. In one, the solute was in the bilayer center in accord with the experiments that indicate that uncharged drug molecules can occupy the center of the bilayer.27 In the other, the nifedipine analogue was close to the hydrocarbodwater interface in accord with suppositions based on the X-ray scattering data. Methods The details of the methodology are described in accounts of our studies of monolayers, l6,I7 bilayers, and benzene diffusion in bi1a~ers.l~The system consisted of 36 dimyristoylphosphatidylcholine (DMPC) lipid molecules arranged into two monolayers (18 in each layer) with lateral dimensions of 34.5 A by 34.5 A in the x-z plane. Each monolayer was solvated with 481 water molecules (which represents full hydration). Two dimensional periodic boundary conditions were employed in the x-z plane, and repulsive walls were used to restrict the water molecules in the third dimension and maintain the appropriate density. The system contained a single solute, an analogue of nifedipine (Figure 1). An all-atom representation of all molecules, including hydrocarbon hydrogen atoms, was employed, because it has been s h o ~ n ~ that * - ~using ~ a united atom model for studies of the diffusion of small solutes in polymers resulted in artificially 15918
H4
I
Figure 1. Nifedipine analogue. C13, C14, C17, C18, C19 are methyl
groups. Phenyl hydrogen atoms are not shown. Hydrogens were numbered according to the number of the non-hydrogen atom to which they were attached (e.g., H191 was attached to carbon 19). A methyl group replaced the nitro group of nifedipine. For these simulations, a pyridine rather than the dihydropyridine ring of nifedipine was used. TABLE 1: Solute Partial Charge# c1 c2 c3 N4 H4 c5 C6 c7 C8 c9 H9 c10 H10 c11 H11 c12 H12 C13 H131 H132 H133 C14
0.0000 0.1792 0.1100 -0.5000 0.2800
0.1100 0.1792 0.000 0.0000 -0.1000 0.1000
-0.1000 0.1000 -0.1000 0.1000 -0.1000 0.1000 -0.3000 0.1000 0.1000 0.1000 -0.3000
H141 H142 H143 C15 01 C16 02 03 04 C17 H171 H172 H173 C18 H181 H182 H183 C19 H191 H192 H193
0.1000 0.1000 0.1000 0.2308 -0.3800 0.2308 -0.3800 -0.1800 -0.1800 -0.1500 0.1000 0.1000 0.1000 -0.1500 0.1000 0.1000 0.1000 -0.3000
0.1000 0.1000 0.1000
a See Figure 1 for atom names. Hydrogens were numbered according to the number of the non-hydrogen atom to which they were attached (e.g., H191 was attached to carbon 19).
high diffusion rates. This was also the conclusion reached from studies of n - a l k a n e ~and ~ ~ solvents such as toluene and tetrah ~ d r o f u r a n .The ~ ~ force field for the lipids is the same as that employed in the work of Stouch et al.33 except that the hydrocarbon chain 1,4 interactions were scaled by 0.6. The water potential employed was a flexible version of the SPClike potential34 which has been shown to be in good agreement with both experiment and the computed properties of other water models.35 For the nifedipine analogue the CVFF force field36137was employed and the partial atomic charges (Table 1) were obtained by a toplogy-based “bond-increment” method calibrated on electrostatic potential derived point
charge^.^^,^^ The Verlet algorithm was used to integrate the equations of motim40 using a time step of 1 fs. The temperature of the entire system was maintained at 320 f 3 K (approximately 22 K above the gel L a phase transition temperature for DMPC bilayers) for all simulations by coupling the system to an external temperature bath?l The temperature of the lipids and solute were scaled separately from the water. A further check on the thermal equilibration was performed by computing the temperature of the system in 1 A slabs parallel perpendicular to the
Alper and Stouch
5726 J. Phys. Chem., Vol. 99,No. 15, 1995
bilayer plane. The temperature of all slabs was within k5 K of the target temperature 320 K, showing that the simulation did not produce any “hot spots” in the system which could lead to artifacts in computed properties of the system. One of the most important components of the simulation methodology is the choice of the nonbonded cutoffs. Simulations of a variety of biomolecular systems16J7*42-43 have shown that truncation of the potential and forces can lead to significant perturbations in the structure and dynamics of the system and that these perturbations can be eliminated only through the use of long cutoffs or effective means such as the Ewald sum or reaction field. Therefore the following protocol was employed for the cutoffs. Interactions between water, the lipid hydrocarbon chains, and the glycerol moieties were treated using “group-based” cutoffs as described previously16(neutral groups of a small number of atoms are used, where one atom in each group is chosen to evaluate group-group distances.) For all interactions involving the zwitterionic phosphocholine headgroup of the lipids (except for interactions with the nifedipine analogue) a “re~idue-based”~*’~ cutoff was used (neutral groups composed of larger numbers of atoms are used, where residueresidue and residue-group interactions are based on the closest :tom-atom pair). If at least one pair of atoms was within 10 A, all interactions between the residues were calculated. Because of the larger size of the “residues”, this implies an “effective” cutoff of -16 8, as described in previous work.16J7 For all interactions involving the nifedipine analogue a “residuebased” cutoff of 20 8,,corresponding to an “effective” cutoff of about 30 A, was used. All interactions were truncated at the cutoff; no switching of the potentials and forces was employed since studies of a variety of systems have shown that switching functions introduce artifacts into the properties computed from simulation^^^^^^^^^. The nonbonded neighbor list was evaluated to a distance 2 8, longer than the potential and force cutoffs. This “buffer” ensured that atoms did not frequently move in and out of the cutoff. This allowed a list update frequency of 50 steps to be employed. The starting bilayer configuration was taken from a wellequilibrated simulation of a fully hydrated DMPC bi1a~er.l~ The solute was incorporated into the bilayer as follows. Making use of the substantial amount of free space within the bilayer,18 first the rings were manually inserted into free spaces in the bilayer then the substituents were “grown” onto the rings rapidly, within 10 ps. The systems were minimized to remove any residual unfavorable interactions, heated to 320 K in increments of 50 K ps, and equilibrated for another 100 ps. During the 2 ns production runs configurations were saved every 0.1 ps for analysis. Two simulations were completed. In one the solute was wholely placed into the center of the bilayer. In the other, the solute was placed with the dihydropyridine ring closer to the region of the hydrocarbon/water interface such that its amine and carbonyl groups could interact with water and lipid headgroups and with the phenyl ring pointing toward the bilayer center. The behavior of the drug analogue-membrane system was analysed by calculating various properties of the solute and the lipids and water, as before.19 The diffusion coefficient of the nifedipine analogue was obtained from the mean-squared displacement of the solute center of mass. To further elucidate the details of the molecules motion, the distance moved by the solute center of mass in Consecutive 5 ps intervals was studied. To determine the extent to which jumps between voids in the bialyer contributed to the diffusion of the drug analogue, the voids present in a neat DMPC bilayer were analyzed. The location and environment of the solute was studied relative to
TABLE 2: X, Y, and 2 Ranges of the Solute simulation Ax AY center interfacial
12.59 10.92
5.38 5.63
A2 17.83 25.34
the distributions of the lipid atoms along the normal to the bilayer plane. The immediate environment of the nifedipine was established by tabulating which lipid and water atoms approached within 5 8, of any atom of the solute. Other properties are described in the body of the Results. The response of the bilayer to the solute was studied as before19 through determination of bilayer thickness and hydrocarbon chain order parameters. ‘ All simulations were performed using a heavily modified version of the DISCOVER simulation package.44 The simulations took 1.8 s/step on a Cray YMP/2E, for a total of 2000 cpu hours for the 4 ns of simulation time. A prototype version of the molecular analysis program Decipher was used for the bulk of the analysis.
Results Bilayer Center. During the 2 ns trajectory, the solute translated substantially throughout the hydrocarbon chain region and rotated slowly. The environment of the solute consisted solely of the hydrocarbon chains, resulting in an essentially homogeneous, hydrophobic environment. The probability of neighboring the solute was highest for the terminal CH2 or CH3 groups (C12-C14), smaller for the methylene groups in the chain center (C8-Cll), and essentially zero for the upper portions of the chains (Cl-C3). The translational motion of the solute is shown in Table 2. During the course of the simulation, the solute oscillated perpendicular to the bilayer plane by only f 3 8, (center of mass), and so the solute always remained near the bilayer center. Note that this is in contrast to our simulations of the smaller, more isotropic benzene molecule, which moved 2-4 times this amount in half the time. Parallel to the bilayer plane the solute explored space more extensively, ranging 13 and 18 A in the two directions, ranges equivalent to those traversed by benzene in half the time. The diffusion coefficient D resulting from this motion was 1.3 x cm2/s. (Comparison of the values for D at 1, 1.5, and 2 ns showed that this quantity was well converged). This is consistent with experimental studies which showed that the diffusion constant in phospholipid membranes of solutes of a size similar to that employed here is on the order of cm2/ s6 and that the diffusion of larger, more complex solutes in lo-* cm2/s.24,45.46 This phospholipid bilayers is of order coefficient is about one-third that of benzene at the same location (3-4 x cm2/s).19 Note, however, that based on molecular weight alone DnifedipineanaloguJn~ene should be closer to onehalf since the nifedipine analogue is about 4 times as massive as benzene. That the observed ratio is smaller is no doubt partly due to the shape of the present solute. In addition, as will be discussed below, it also appears that its increased size over that of benzene alters its mechanism of diffusion. The orientation and rotation of the solute were investigated by the behavior of the vector formed by the amine nitrogen (N) and the opposing carbon atom (Cl) of the dihydropyridine ring. The angle formed by this vector with the bilayer normal as a function of time is shown in Figure 2a, where it can be seen that this angle explores most of the possible range of values, from almost parallel to almost antiparallel to the bilayer normal. The autocorrelation function of this vector (Figure 2B) demonstrates the slow decorrelation over the simulation and suggests that the time scale of this motion is on the order of nanoseconds.
J. Phys. Chem., Vol. 99, No. 15, 1995 5727
Orientation and Diffusion of a Drug Analogue
35
g
50
55
60
Figure 3. Atom probability distributions along the bilayer normal. Light lines, lipid: solid (P), phosphorus; dotted (N), nitrogen; long carbonyl oxygens; solid (M) methyl groups. Heavy lines, dashes (0)
-
solute: solid, amine nitrogen; dotted, 4-phenyl carbon; dashed, 4-pyridyl carbon.
00-
4.J
45
Posilion(angstroms)
IO -
0.1
40
-
0
So0
IWO
1x4
2”
Time@)
Figure 2. (a, top) Angle formed by the C1-N vector with the bilayer normal as a function of time, for the solute in the bilayer center. (b, bottom) Autocorrelation function of the C1-N vector, obtained as the expectation value of the normalized first Legendre polynomial, for the solute in the bilayer center.
The torsion angles of the COOCH3 side chains were also investigated. For the first dihedral angle out from the dihydropyridine ring (formed by the atoms Cl-C2-C15-04 and Cl-C6-C1-603) few interconversions (less than 3) between different conformers were observed, and “residence times” for different conformational states of as small as 200 ps and as large as 1300 ps were observed. The second torsions from the dihydropyridine ring (formed by the atoms C2-Cl5-04-Cl8 and C6-Cl6-03-Cl7) remained in the trans state for the entire 2 ns simulation. The outermost torsions from the pyridine ring, involving rotation of the methyl groups, experienced frequent interconversion between rotameric states; as many as 40-45 interconversions during the course of the 2 ns simulation resulting in “residence times” as small as 25 ps and as large as 200 ps. This rapid conformational interconversion is consistent with the low energetic barriers to rotation about terminal methyl groups. “Interfacial” Region. A significant body of research indicates that some drug molecules situate in the hydrocarbodwater interface, suggesting that this might position them to bind to membrane-bound receptor^.^^-*^ However, as stated above, these experiments lack the resolution to determine the exact location of the drug, its orientation, and the details of its intermolecular interactions. During this simulation, the solute moved substantially parallel to the bialyer plane and “bobbed” perpendicular to the plane. On average, the solute as a whole was tilted (defined by the angle between the C 1-N vector and the bilayer normal), and it rotated to explore a range of orientations with respect to the
bilayer normal. However, this reorientation of the solute was substantially slower than when the solute was in the bilayer center. Within this region, the solute experienced a rich and fluctuating environment, which included hydrogen bonding with the headgroups and carbonyl groups of the lipids as well as with water. During the course of the simulation the center of mass of the solute moved perpendicular to the bilayer plane over a 6 8, range (Table 2). The range of locations of the solute in the membrane is shown in Figure 3 relative to the distributions for the locations of various lipid atoms. The width of these distributions results from the significant freedom of motion experienced within the L a phase. The amine group of the solute overlaps significantly with the region covered by the lipid carbonyl groups and to a lesser degree (covering about a 4 8, range) with the region spanned by headgroup atoms such as phosphorous and nitrogen. This is in agreement with the results of X-ray studies of the location of various nifedipine analogues in phospholipid membra ne^,^^-^^ which show that, on averge, the dihydropyridine rings of many 1A-dihydropyridine calcium channel antagonists (1,4-DHPs) situate themselves at the “hydrocarbon core/water interface”,just below the phosphocholine headgroups. The solute’s phenyl ring penetrated significantly into the region of the terminal methyl groups and, with small probability, advanced into the region of the lipid carbonyl groups. These results are also consistent with NMR experiments27of uncharged 1,4-DHP drug molecules in membranes of selectively deuterated lipids that indicate that such solutes explore different regions of the bilayer. Increases in the quadropole splittings of all deuterated segments with increasing drug concentration implied a homogeneous distribution of the drug throughout the bilayer. These distributions are a reflection of the motion of the solute as well as the lipids. The solute also translated extensively parallel to the bilayer plane to an extent even greater than that in the bilayer center. The range of motion was again about the same as benzene experienced in one-half the time except in the z direction, where it was somewhat greater than this. The diffusion coefficient in this region was 1.3 x cm2/s (analysis of the diffusion constant as a function of simulation time showed that this quantity was well converged by the full simulation time of 2 ns), essentially the same as in the bilayer center. This is in contrast to the results we obtained for the smaller benzene
Alper and Stouch
J. Phys. Chem., Vol. 99, No. 15, 1995
m1
sn2
Figure 5. Frequency with which different lipid atoms approached within 5 8, of any atom of the solute when near the interface. (CGl, glycerol carbon atoms; C1 and 0 1 , carbonyl carbon and oxygen; C14, terminal methyl carbons; 0 2 , ether oxygens; 03-06, phosphate oxygens).
-1
8
00
0
JW
Im,
ISm
ma
nm(ps)
Figure 4. (a, top) Angle formed by the C1-N vector with the bilayer normal as a function of time, for the solute near the interface. (b, bottom) The autocorrelation function of the C1-N vector, obtained as
in Figure 2b, for the solute near the interface. molecule which diffuses at different rates in the two regions;lg a mechanism to explain this different behavior will be proposed below. The orientation and rotation of the solute were determined, as before, using the C1-N vector. Its angle with the bilayer normal vs time (Figure 4a) explored a more limited range in the interfacial region than in the bilayer center and ranged from roughly perpendicular to the bilayer normal (-90") to approximately parallel to the normal (-180"). On average the nifedipine analogue was tilted with respect to the bilayer normal. During the 2 ns simulation, the autoconelation function of the C1-N vector (Figure 4b) never decreases below 0.5, showing that the rate of rotation and reorientation is much slower in this region than in the bilayer center and that the time scale of this motion is at least on the order of nanoseconds. The nifedipine experienced a rich and diverse environment in the hydrocarbodwater interfacial region in contrast to that in the bilayer center. The total number of lipid molecules that had any atom within 5 8, of the solute varied from 4 to 10. The frequency with which different lipid atoms were within 5 8, of any atom of the solute is plotted in Figure 5 (accumulated for 40 configurations taken at 50 ps intervals). The environment of the solute primarily consisted of the carbonyl groups, the glycerol backbone, and nearby methylene units (C2-C5) of the hydrocarbon chains, in agreement with Figure 3. The solute also neighbors the lipid headgroups, though to a smaller degree. The number of water molecules within 5 8, of the solute varied from 6 to 22 during the simulation. On average the N-H group neighbored one lipid carbonyl group, though the number of such interactions ranged from 0 to 2 during the simulation. In particular, one lipid contributed
its carbonyl group for 1500 ps of the simulation, while another lipid contributed its carbonyl group for 500 ps. This shows that the interaction between the N-H group of the drug analogue and the lipid carbonyl group is strong and of long duration. The N-H group also neighbored the phosphoether oxygen atom of one lipid but for a much shorter time period, about 200 ps. The short duration of the interaction may be related to the relatively small degree of penetration of the dihydropyridine ring into the region occupied by the lipid headgroups (see Figure 3) and also to the fact that the phosphoether oxygen is sterically crowded, being bonded to both the phosphorus and the glycerol. That the amine does not hydrogen bond to the less crowded and more strongly charged phosphoryl oxygens is probably mostly a reflection of the relative positions of those atoms. The phosphoryl oxygens are almost always oriented away from the bulk of the bilayer and toward the water layer, with which it strongly hydrogen bonds. On average 2-3 water molecules were within 5 8, of the amine moiety, although the number ranged from 0 to 5 during the simulation. Many of these neighboring water molecules exchanged within 5 ps, though some water molecules resided about the amine moiety for as long as 30 ps. This is similar to the residence times computed for water about the polar phosphocholine headgroups of the lipids in a DMPC monolayer,17 where residence times of 10 and 20 ps were observed for water about the phosphate and choline moieties, respectively. The two carbonyl groups of the solute also had significant interactions with water, though the number of neighboring water molecules was substantially different for the two because of the overall tilt of the solute. One carbonyl group had on average 2-3 neighboring water molecules (though the number varied between 0 and 6), and these water molecules exhibited frequent exchange, as about the N-H grouf. The other carbonyl group had no water molecules within 5 A for most of the simulation, but for a 200 ps interval there was one neighboring water molecule. The identity of this water also varied due to frequent exchange. Interestingly the torsions of the ester groups behaved differently in this location than in the bilayer center. The dihedrals of the two ester groups immediately off of the pyridine ring (rotation about the single bond to the carbonyl carbon) were highly correlated during the simulation such that when one changed from trans to cis (the carbonyl group stayed coplanar with the pyridine ring) the other changed oppositely at the same time. No such correlation was observed when it was in the center. One possible explanation for this behavior is that it maximizes the separation between the carbonyl groups and so
Orientation and Diffusion of a Drug Analogue
J. Phys. Chem., Vol. 99, No. 15, 1995 5729 0.6
0.5
0.4
m
E
5
P
0.3
B
e
0
0.2
0.1
10
15
25
20
30
X
Figure 6. Position in the bilayer plane of the solute (left side, represented by the C7 atom coordinates) and one lipid (right side, coordinates of the phosphorus atom) for the simulation of the solute near the interfacial region. The data points are 150 ps apart and are numbered sequentially in time for clarity.
allows for the optimum contact of these groups with water. These torsions underwent fewer interconversions between trans and cis conformations in the interfacial region than when it was in the bilayer center and residence times were correspondingly larger, 700 and 1500 ps, respectively. The next torsions out from the ring (about the bond linking the carbonyl carbon to ether oxygen) stayed in the trans state throughout the 2 ns run, as in the simulation of the analogue in the bilayer center. The rotational behavior about the ester methyl groups was similar for the solute in both locations, and the rotamer residence times ranged from 25 to 200 ps. However, in the interfacial region the number of interconversions for the two terminal methyl groups were unequal (-35 and 65 interconversions, respectively) unlike in the center (-40 and 45, respectively). In particular, the greater number of torsional interconversions occurred for the terminal methyl group associated with that solute carbonyl group which, because of the overall molecular tilt, had more neighboring water molecules. Dynamical Events at the Interface. The diverse environment about the solute in the interfacial region resulted in several interesting dynamic events. At about one-third through the simulation, the interactions of the nifedipine analogue with its surrounding altered significantly due to two simultaneous changes. The orientation of the solute as a whole changed suddenly and moved one of its ester groups closer to the water. At the same time, the C Tst torsions off the dihydropyridine ring in the ester groups underwent concerted conformational changes which pointed the carbonyl oxygen of the above-mentioned ester directly into the water and increased the number of neighboring water molecules by 33-100%. This transition did not change the number of adjacent lipid carbonyl and phosphate oxygen atoms. The overall binding energy of the solute became 50% more favorable. The resulting orientation of the molecule was at 100” from the bilayer normal with one ester group pointing into the water and another pointing toward the bilayer center. Also interesting was the lipid motion relative to the solute motion. The motion of several neighboring lipids appeared to be correlated with the soIute motion to varying degrees. Figure 6 shows a high correlation between the positions within the bilayer plane of one such lipid (phosphorus position) and the solute (atom C7 which is near the center of the solute.) for the entire trajectory. The xosition of another lipid close to the solute was similarly correla’ed with that of the solute, but for only portions of the trajectory. These differing degrees of correlation between the solute and lipid motion no doubt reflect different lipid-solute binding arrangements. This also implies that the
0.0
.
.
.
.
.
.
.
.
1
2
3
4
3
6
7
8
. 9
.
.
10
11
I2
Carbon Number
Figure 7. Order parameters for the methylene units of the lipid hydrocarbon chains: simulation of the neat bilayer (solid line). simulation of the solute in the interfacial region of the bilayer (dots), and simulation of the solute in the bilayer center (dashed).
nifedipine analogue may affect the rate of diffusion of lipids that are strongly bound to it. Effect of the Solute on the Membrane. Similar to the results of our previous studies of benzene in DMPC bilayers,lg the nifedipine analogue did not affect gross changes in the properties of the bilayer, compared to the neat bilayer, though there may have been small perturbations in some properties. For example, the bilayer thickness was found to be 25.2 & 0.2 8, (sd) and 25.4 ?C 0.2 8, for nifedipine in the center and interfacial region respectively, which agrees with the 25.3 f 0.3 (sd) A thickness of the neat bi1a~er.l~ The number of gauche torsions per hydrocarbon chain was 3.0 k 0.2 (sd) (center) and 3.2 f 0.2 (interface), within error of the range 2.7-3.3 for neat DMPC bi1a~ers.l~The distributions along the bilayer normal of the lipid atoms (Figure 3) were essentially identical to each other and those of the neat bi1a~er.l~ As with benzene,19 there appeared to be small changes in the order parameters between the two simulations of this study and the neat bilayer (Figure 7). The solute in either region decreased the order parameters near the interface relative to those of the neat bilayer. In the middle of the chains (C4-C8) the results depended on the location of the solute; when it was near the interface, the solute increased the order and when near the center, it decreased the order. When the solute is at the interface, the C4-C8 region is coincident with at least the solute’s phenyl ring. “Packing” with the solute could prevent the neighboring lipid hydrocarbon chains from assuming conformations which would otherwise pass through the volume occupied by the solute. This could cause them to straighten and assume fewer gauche angles in their dihedral angles which would result in increased order parameters of their methylene groups. When at the bilayer center, the solute would prevent neighboring lipids’ hydrocarbon chains from assuming straightened conformations which would otherwise allow them to penetrate the volume occupied by the solute. To compensate, the chains would have to tilt or assume greater numbers of gauche torsions, both of which would lower their order parameters. The order at the chain ends (C9-C12) were essentially unchanged by the presence of the solute, even when the solute was in the center. This suggests that the amount of disorder in the center, due to rapid torsional isomerization and substantial free volume, is great enough to be unperturbed even by a solute the size of nifedipine. The average time between rotameric interconversions of the hydrocarbon torsion angles for the solute in the center, the interfacial region, and
5730 J. Phys. Chem., Vol. 99, No. 15, 1995
the neat bilayer were of 56.1, 55.9, and 55 ps, which might suggest a decrease in interconversion rate when the solute was in the center. Discussion. Perhaps the most important result of the present work coneems the relationship between solute size, the variation of the solute diffusion constant with position in the bilayer, and the mechanism of diffusion. In particular, although the diffusion of benzene in the bilayer center was found to be 3 times faster than diffusion in the interfacial region,lg the nifedipine analogue diffused at the equivalent rates in either region. While early theories of solute permeation related it to the partition coefficient of the molecule into a hydrophobic solvent,* low molecular weight molecules were found to diffuse faster than predicted by Overton’s The “soft polymer” theory states that in a matrix enhanced rates of diffusion for small solutes are due to jumps between voids in the matrix of sufficient size to contain the solute^.^*^ This theory seems to apply to the rate of diffusion of benzene in the bilayer, where the enhanced rates of diffusion could be correlated to increased numbers of voids and increased numbers of jumps. Studies of the location and size of voids in a neat bilayer18 found a high probability of voids on the order of the volume of a benzene molecule (> 80 A3) in the bilayer center, whereas the interfacial region contained little free volume and no voids approaching the volume of benzene. However, voids seldom approached the size of the nifedipine analogue even in the center. Therefore, unlike benzene, enhanced rates of diffusion via the hole-jumping mechanism were unavailable to the nifedipine analogue in any region, and we conclude that benzene and the nifedipine analogue diffuse within the membrane by different mechanisms. The present results are in agreement with another recent computer simulation study of the diffusion of water through a DPPC bilayer?O which showed Dcen&Dhtcfiace 10. Increasing solute size, from water to benzene to the nifedipine analogue, decreases the diffusion gradient, so that, finally, for solutes of sufficient size, diffusion will occur at a similar rate independent of the solute location. The atomic-level resolution of these studies can help in the interpretation of experimental data. For example, although other configurations were sampled, the most energetically favorable and longest lasting configuration sampled during this simulation, that with the greatest hydrogen bonding between the solute, the lipids, and the water, was tilted substantially along its long axis relative to the bilayer normal. This contrasts with the assumptions often employed in X-ray studies of the location of drug molecules in membra ne^.^^-^^ Although such studies can approximate location, orientation cannot be determined and can only be assumed. The assumption most commonly employed is that the C 1-N vector (or the equivalent) in a drug molecule aligns parallel to the bilayer normal, presumably because such an orientation is most energetically favored. The present results show that the solute orientation can differ significantly from this assumption. The amine group of the solute is high enough in the interfacial region to neighbor lipid carbonyl groups and water for a variety of solute orientations. However, the ability of the solute carbonyl groups to neighbor water is more sensitive to solute orientation since they are 2-3 8, lower than the amine group and are located on side chains of the solute. When the C1-N vector of the solute is parallel to the bilayer normal the solute carbonyl groups are too low to neighbor water and so little interaction with water is seen. When the solute has a significant tilt one of the carbonyl groups is raised sufficiently high that it can interact with water. This means that a greater number of energetically favorable interactions with water are
-
Alper and Stouch possible when the solute is tilted than when the C1-N vector is parallel to the bilayer normal.
Conclusions Drudmembrane interactions were probed through over 4 ns of molecular dynamics simulations of a nifedipine analogue in an all-atom, fully hydrated, DMPC bilayer in the L a phase. Long-range electrostatic interactions were included to a distance of 30 8, to minimize artifacts of the simulation. The calculated rate of diffusion was quite close to that expected. Unlike the smaller benzene molecule, this rate did not change with the drug analogue’s, position within the bilayer. This difference could be rationalized based on the volumes of the solutes and those of spontaneously arising voids within the bilayer which are often of the size of benzene but essentially never the size of nifedipine. The rate of benzene diffusion increases when it can make large jumps between these voids. These voids are predominantly in the bilayer center which gives rise to the different rates of diffusion. The analogue explored regions of the bilayer by “bobbing” in a perpendicular direction and translating even more extensively later to the bilayer plane. When situated near the bilayer/ water interface, the amine moiety of the analogue hydrogen bonded extensively with the carbonyl oxygens of the fatty acid ester linkages and water and transiently with phosphoester oxygens. When the long axis of the solute was perpendicular to the bilayer plane, its carbonyl esters made hydrogen bonds with water only with difficulty. However, it tilted to make a substantial number of such interactions with one of the ester groups. This suggests that a preferred orientation of nifedipinelike compounds orientation within biological membranes will be tilted relative to the bilayer normal. Although when near the interface it essentially always formed hydrogen bonds with water, the waters exchanged freely and had residence times of only several tens of picoseconds. Several different relationships were seen with the lipid molecules. Although some of them exchanged freely, the motion of one was correlated with the solute through most of the simulation. Although the solute experienced considerable translational and conformational motion during the two 2-ns simulations, the time scale of whole-body rotation of the solute was on the order of nanoseconds, regardless of location. Future studies will attempt to probe the nature of these longer time motions and their effects on drug/membrane interactions.
Acknowledgment. The authors would like to thank Dr. Donna Bassolino for many helpful discussions. The authors also thank Richard Shaginaw, John Stringer (Cray Research Inc.), Richard Gopstein and Gregory Burnham of the BristolMyers Squibb High Performance Computing Center, and Stewart Samuels (Bristol-Myers Squibb Macromolecular Modeling Department), who play key roles in our work by orchestrating, updating, and maintaining our Cray-Y/MP and Silicon Graphics computing network. References and Notes (1) Overton, E. Vierteljahrsschr. Natulforsch. Ges. Zuerich 1899,44,
88. (2) Cohen, B. E. J . Membr. B i d . 1975,20, 235. (3) Walter, J.; Gutknecht, J. J . Membr. B i d . 1986,90, 207. (4) Cohen, M.H.; Tumbull, D. J . Chem. Phys. 1959,31, 1164. (5) Lieb, W. R.; Stein, W. D. Nature 1969,224, 240. (6) Dix, J. A.; Kivelson, D.; Diamond, J. M. J . Membr. B i d . 1978, 40, 315. (7) DeYoung, L. R.; Dill, K. A. Biochemistry 1988,27, 5281. (8) Jonsson, B.; Edholm, 0.;Teleman, 0. J. Chem. Phys. 1986,85, 2263.
Orientation and Diffusion of a Drug Analogue (9) Wendoloski. I. J.: Kimatian. S. J.: Schutt. C. E.: Salemme. F. R. Science 1989, 243, 636. (10) Eeberts. E. Ph.D. Thesis. Riiksuniversiteit. Gronineen. 1988. (11) fatanabe, K.; Ferrario, M.; klein, M. L. J . Chem. Fhjrys. 1988,92, 818. (12) Egberts, E.; Berendsen, H. J. C. J . Chem. Phys. 1988, 89, 3718. (13) Berkowitz, M. L.; Raghavan, K. Lungmuir 1991, 7, 1042. (14) Raghavan, K.; Rami Reddy, M.; Berkowitz, M. L. Longmuir 1992, 8, 233. (15) Stouch, T. R. Mol. Sim. 1993, IO, 335. (16) Alper, H. E.; Bassolino, D.; Stouch, T. R. J . Chem. Phys. 1993, 98, 9798. (17) Alper, H. E.; Bassolino, D.; Stouch, T. R. J . Chem. Phys. 1993, 99, 5547. (18) Stouch, T. R.; Alper, H. E.; Bassolino, D. Int. J. Supercomput.App. 1994, 8, 6-23. (19) Bassolino, D.; Alper, H. E.; Stouch, T. R. Biochemistry 1993, 32, 12624. (20) Marrink, S.-J.; Berendsen, H. J. C. J . Phys. Chem. 1994,98,41554168. (21) Sok, R. M.; Berendsen, H. J. C.; van Gunsteren, W. F. J. Chem. Phys. 1992, 96, 4699. (22) Rhodes, D. G.; Sarmiento, J. G.; Herbette, L. G. Mol. Pham. 1985, 27, 612. (23) Herbette, L. G.; Trumbore, M.; Chester, D. W.; Katz, A. M. J . Mol. Cell. Cardiol. 1988, 20, 373. (24) Mason, R. P.; Chester. D. W. BioDhvs. J. 1989, 56. 1193-1201. (25) Herbette, L. G.; Vant Ewe, Y. M.'H.; Rhodes, D. G. J . Mol. Cell. Cardiol. 1989, 21, 187. (26) Mason, P. R.; Moring, J.; Herbette, L. G. Nucl. Med. Biol. 1990, 17, 13. (27) Bauerle, H.-D.; Seelig, J. Biochemistry 1991, 30, 7203. (28) Muller-Plathe, F.; Rogers, S. C.; van Gunsteren, W. F. Chem. Phys. Lett. 1992, 199, 237-243. (29) Yoon, D. Y.; Smith, G. D.; Matsuda, T. J . Chem. Phys. 1993,98, 10037- 10043.
J. Phys. Chem., Vol. 99,No. 15, 1995 5731 (30) Pant, P. V. K.; Boyd, R. H. Macromolecules 1993, 26, 679-686. (31) Ryckaert, J. P.; Klein, M. L. J . Chem. Phys. 1986, 85, 1613. (32) Bareman, J. P.; Reid, R. I.; Hrymak, A. N.; Kavassalis, T. A. Mol. Sim. 1993, 11, 243-250. (33) Stouch, T. R.; Ward, K. B.; Altieri, A,; Hagler, A. T. J . Comput. Chem. 1991, 12, 1033-1046. Williams, D. E.; Stouch, T. R. J. Comput. Chem. 1993, 14 (9), 1066-1076. (34) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Inremolecular Forces; Pullman, B., Ed.; Reidel; Dordrecht, Holland, 1981; p 331. (35) Lau, K.; Alper, H. E.; Thacher, T.; Stouch, T. R. J. Phys. Chem. 1994, 98, 8785-8792. (36) Hagler, A. T.; Lifson, S.; Dauber, P. J . Am. Chem. SOC.1979, 101, 5122. (37) Dauber-Osguthorpe, P.; Roberts, V. A,; Osguthorpe, D. A,; Wolfe, J.; Genest, M.; Hagler, A. T. Proreins: Structure, Function, Genetics 1988, 4, 31. (38) Maple, J. R.; Hwang, M.-J.; Stockfkch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J . Comput. Chem. 1994, 15, 162. (39) Oie, T.; Maggiora, G. M.; Christoffersen, R. E.; Duchamp, D. J. Int. J . Quantum Chem. Quantum Biol. Symp. 1981, 8, 1. (40) Verlet, L. Phys. Rev. 1967, 159, 98-103. (41) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Di Nola, A.; Haak, J. R. J . Chem. Phys. 1984, 81, 3684. (42) Alper, H. E.; Levy, R. M. J . Chem. Phys. 1989, 91, 1242. (43) Smith, P. E.; Pettitt, B. M. J . Chem. Phys. 1991, 95, 8430. (44) Discover biomolecular simulation package v 2.60, BIOSYM Technologies, Inc., 10065 Barnes Canyon Road, Suite A, San Diego, CA 92121. (45) Wu, E.-S.; Jacobson, K.; Papahadjopoulos, D. Biochemistry 1977, 16, 3936. (46) V u , W. L. C.; Clegg, R. M.; Hallman, D. Biochemistry 1985,24, 781. JP942977U
