Atomistic Simulations of Pore Formation and Closure in Lipid Bilayers

3 downloads 0 Views 2MB Size Report
0006-3495/14/01/0210/10 $2.00 http://dx.doi.org/10.1016/j.bpj.2013.11.4486 ...... Soft Matter. 5:3295–3302. 15. Kornberg, R. D., and H. M. McConnell. 1971.
210

Biophysical Journal

Volume 106

January 2014

210–219

Atomistic Simulations of Pore Formation and Closure in Lipid Bilayers W. F. Drew Bennett, Nicolas Sapay, and D. Peter Tieleman* University of Calgary, Department of Biological Sciences and Centre for Molecular Simulation, Calgary, Alberta, Canada

ABSTRACT Cellular membranes separate distinct aqueous compartments, but can be breached by transient hydrophilic pores. A large energetic cost prevents pore formation, which is largely dependent on the composition and structure of the lipid bilayer. The softness of bilayers and the disordered structure of pores make their characterization difficult. We use moleculardynamics simulations with atomistic detail to study the thermodynamics, kinetics, and mechanism of pore formation and closure in DLPC, DMPC, and DPPC bilayers, with pore formation free energies of 17, 45, and 78 kJ/mol, respectively. By using atomistic computer simulations, we are able to determine not only the free energy for pore formation, but also the enthalpy and entropy, which yields what is believed to be significant new insights in the molecular driving forces behind membrane defects. The free energy cost for pore formation is due to a large unfavorable entropic contribution and a favorable change in enthalpy. Changes in hydrogen bonding patterns occur, with increased lipid-water interactions, and fewer water-water hydrogen bonds, but the total number of overall hydrogen bonds is constant. Equilibrium pore formation is directly observed in the thin DLPC lipid bilayer. Multiple long timescale simulations of pore closure are used to predict pore lifetimes. Our results are important for biological applications, including the activity of antimicrobial peptides and a better understanding of membrane protein folding, and improve our understanding of the fundamental physicochemical nature of membranes.

INTRODUCTION Pores are important but transient structures in membrane biology and in biotechnology. Membranes function to separate cells and cellular compartments, which is necessary for energy transduction, metabolism, cellular signaling, and life in general. Disruption of the membrane can cause cell death, and is a possible mechanism for antimicrobial peptides (1–3). Most drugs (many charged and/or polar) must be designed to cross lipid membranes to reach intracellular targets, and might cross more quickly through membrane pores. Applying a large potential difference across membranes can cause pore formation, with potential therapeutic applications (4). Cell-penetrating peptides are short, polycationic peptides that can translocate cargo molecules through membranes and into cells, possibly through a pore-mediated mechanism (5). Amyloid peptides have been shown to induce pore formation in model membranes (6). The release of cytochrome c from the mitochondria is a key step in apoptosis. Cytochrome c formed pores in model membranes free of other protein, but that contained cardiolipin (7). A lipocentric view of all of these processes suggests a common

Submitted September 20, 2013, and accepted for publication November 25, 2013. *Correspondence: [email protected] This is an Open Access article distributed under the terms of the Creative Commons-Attribution Noncommercial License (http://creativecommons. org/licenses/by-nc/2.0/), which permits unrestricted noncommercial use, distribution, and reproduction in any medium, provided the original work is properly cited. W. F. Drew Bennett’s current address is University of Waterloo, Department of Chemistry, 200 University Avenue West, Waterloo ON, N2L 2G1. Nicolas Sapay’s current address is BIOASTER Technology Research Institute, baˆtiment Domilyon, 321 avenue Jean Jaure`s, 140169007 Lyon, France.

mechanism and underlying energetics. Understanding pore and defect formation in membranes is also important for interpreting the results from experiments, such as drug permeation, hydrophobic mismatch of lipids and proteins, and determining amino-acid hydrophobicity scales (8–10). Biological membranes display a large degree of heterogeneity in their lipid composition and resulting structure (11). In model systems, lipids and surfactants have diverse phase behavior and aggregate morphology. Micelles, bilayers, hexagonal phases, cubic phases, and many others are observed, and they depend on many factors, such as the monomers’ chemical structure, temperature, concentration, and aqueous environment. It is clear that the structure of the membrane will affect its deformability, but the relationships among the membrane composition, membrane structure, and energy cost for pore formation are not understood. Defining pore formation at an atomistic level of detail can provide crucial insight into designing strategies for selectively porating specific cells, such as bacteria or cancerous cells (12,13). Pores will form in lipid bilayers spontaneously, but with low probability due to a large free energy cost. Once a pore forms, lipids, ions, and polar molecules are able to cross the membrane. The free energy cost for pore formation is lower than expected from the partitioning of a charge from water to a bulk hydrocarbon due to the formation of water defects and pores in the membrane to keep the charged species hydrated (14–16). This is in contrast to water and other small, polar molecules that are able to cross lipid membranes on a reasonable timescale, and are thought to translocate in a solubility-diffusion mechanism (17). Computer simulations have provided molecular level details regarding hydrophilic pore formation in model

Editor: Bert de Groot. Ó 2014 The Authors 0006-3495/14/01/0210/10 $2.00

http://dx.doi.org/10.1016/j.bpj.2013.11.4486

Pores in Lipid Bilayers

membranes (18,19). Pores have been observed by applying an electrochemical gradient across the membrane (20–26). Antimicrobial peptides and cell-penetrating peptides have been shown to induce pores in atomistic simulations (27–30). Applying a mechanical stress in the plane of the membrane caused pores to form (21,31). Despite differences in the details of pore formation, the underlying structure of the pores is surprisingly conserved, with a disordered toroidal shape. We showed that energy for moving a second arginine to the center of a membrane is much lower than the energetic cost to place the first one at the center (8). The first arginine formed a water defect in the lipid bilayer, allowing the second arginine to share the defect, reducing substantially the free energy cost for the second. This nonadditivity for the free energy of arginine partitioning was also shown experimentally (9). We also showed that arginine would remain charged at the bilayer, in contrast to lysine and the other titratable residues (32). Experiments subsequently confirmed that lysine would likely deprotonate at the center of a DOPC bilayer in contrast to arginine (33). We previously determined the free energy barriers for lipid flip-flop in pure phosphatidylcholine (PC) lipid bilayers for DLPC (dilauroylphosphatidylcholine), DMPC (dimyristoylphosphatidylcholine), and DPPC (dipalmitoylphosphatidylcholine), where we observed pore formation in the thinner bilayers (14). Here we have significantly extended the umbrella sampling simulations to decompose the free energy for DMPC pore formation into the thermodynamic driving forces. In new simulations, we determined pore closure times in the three PC membranes by starting simulations from a preformed pore state and monitoring the pore evolution. We also observed multiple pore formation events in equilibrium simulations of the DLPC bilayer. Here we analyze the energetics and atomistic details of pore formation and discuss implications on biological processes such as antimicrobial peptides, lipid translocation, and transmembrane protein structure and stability.

211 TABLE 1

Simulations performed

Typea DLPC

Close Open DMPC Umb Close DPPC

Close

Replicates/number of umbrellas 100 100 4/20 280 30

Simulation length (ns) 500 50–500 ns, 50–800 nsb 50 50–600 nsb, 150–100 ns, 80–180 nsc 20

Number of close/open 6 4 — 105 30

Total of 155 ms of simulations. a The type of simulation: Umb, umbrella sampling of pore formation; Close, pore closure; and Open, equilibrium simulations to observe pore opening. b Run with a 4-fs time step (see Methods). c Run using GROMACS 3 (see Methods).

Berendsen (40) for temperature coupling, with coupling constants of 0.1 ps. Semiisotropic weak pressure coupling was used (40), with a coupling constant of 2.5 ps, and reference pressures of 1 bar in the plane of the bilayer and normal to the bilayer. Water bonds and angles were constrained with the SETTLE algorithm (41), and all other bonds were constrained with LINCS (42).

Pore formation by umbrella sampling Using a similar protocol as in Sapay et al. (14), we determined the free energy for moving the phosphate of a single DMPC lipid into the center of a DMPC bilayer. A series of 21 independent simulations were run with the phosphate of a single lipid restrained at different depths in the bilayer, with a 0.1-nm spacing between adjacent simulations, from the lipid’s equilibrium position (~2 nm) to the bilayer center (0 nm). A harmonic restraining potential with a force constant of 3000 kJ mol–1 nm–2 was used in the direction normal to the plane of the bilayer. The weighted histogram analysis method was used to determine the free energy profile (43). In this study, we ran simulations with one umbrella potential per simulation, as opposed to two (one in either leaflet), so we could determine the change in enthalpy of the system. Previous tests on DPPC showed no significant difference in the potential of mean force (PMF) between one and two umbrella potentials per simulation (44). The PMFs for DLPC and DPPC were taken from our previous work (14). Table 1 lists the simulation time for each umbrella window and the number of independent PMFs for each system.

Pore closure and rates METHODS Simulation protocol Simulations were run with the software GROMACS, versions 3 (34) and 4 (35). GROMACS 4 was needed to reach hundreds of microseconds of simulations, whereas GROMACS 3 had been used for simulations before its release. For DLPC, DMPC, and DPPC we used the parameters from Berger et al. (36). The bilayers each contained 64 lipids with the initial structures from Sapay et al. (14). Water was modeled with the SPC model (37). We used either a 2-fs or 4-fs time step, updating the neighbor search every 10 steps (2-fs time step) or five steps (4-fs time step). Simulations run using GROMACS 3 or with a 4-fs time step are indicated in Table 1. A larger time step is possible because the PC lipids contain no hydrogens, while the properties of SPC are still reasonable at a time step of 4 fs. Electrostatic interactions were calculated with the particle-mesh Ewald method (38,39), a real-space cutoff of 0.9 nm, a fourth-order spline, and a 0.12-nm grid spacing. Lennard-Jones interactions were cut off at 0.9 nm. We used the

We estimate pore closure time by releasing the restraint on the phosphate that was held at the center of the bilayer, a preformed pore state, and then by observing the time for closure. The pore was from the umbrella sampling simulations with the restraining potential at the center of the bilayer. From the rate of pore closure (kclose) and the free energy difference between the two states, we can estimate the rate for pore opening (kopen) using

  kopen DGpore ¼ exp  ; kclose RT

where kclose is from the simulations of pore closure. The value DGpore is the free energy difference for a single lipid to move from equilibrium to the center of the membrane, which we assume is the free energy for pore formation, given that a pore is observed at this state. The free energy for pore formation in a DPPC bilayer using the pore radius as a collective reaction coordinate was found to be 75–100 kJ/mol (45). The good Biophysical Journal 106(1) 210–219

212

Bennett et al.

agreement with our free energy suggests a common process and supports our assumption that we determine the free energy for pore formation. This gives the rate per lipid, which could be converted to area by dividing by the area per lipid. We can determine the equilibrium pore density (rpore) of the three bilayers, from the free energy for pore formation (DGpore) and the area per lipid (Alip):



rpore

DGpore ¼ exp  RT



Alip :

The unidirectional flux of ions, water, lipids, or any other solute across the bilayer through a transient hydrophilic pore can then be estimated as

J ¼ ji  rpore ; where ji is the flux through the pore, and rpore is the pore density. The permeability coefficient P can then be estimated using

Pi ¼

Ji ; DC

where DC is the concentration difference across the membrane. We assume unidirectional flux, where there is an equal concentration on either leaflet. For unidirectional flux, DC is 55 M for water, and 1/(0.5dAlip) for the lipids, where d is the bilayer thickness. The concentrations are first converted to molecules per cm3. The permeation calculations are based on work in Tieleman and Marrink (46).

RESULTS Pores formation Pore structures are shown in Fig. 1, along with the structure of the unperturbed bilayers (no-pore state). We used umbrella sampling to move the phosphate of a single lipid with respect to the center of mass of the bilayer. Because the phosphate of a single lipid is moved into the bilayer interior, there is a steep slope in the free energy profile (Fig. 1 E). This corresponds to the formation of a water defect: water and lipid headgroups move toward the bilayer interior. At a certain distance from the bilayer center, pore formation becomes favorable on the timescale of the umbrella sampling simulations. Water and phosphates enter from both leaflets, and the bilayer bends into a toroidal shape. This position is the free energy plateau—once a pore forms, the free energy flattens. There is an increase in the free energy cost for pore formation for DLPC (17 5 2 kJ/mol), DMPC (45 5 4 kJ/mol), and DPPC (78 5 2 kJ/mol), as the bilayer thickness increases (Table 2). For DPPC, pore formation occurs only when the phosphate is right at the bilayer center. The PMFs for pore formation are shown in Fig. 1 E, in agreement with our previous results (14). The pores are hydrophilic with water and lipid headgroups pulled into the bilayer interior, forming a disordered toroidal shape. The lamellar state is characterized by nearzero density of water and PC headgroups within the bilayer interior, with water penetrating slightly deeper than the Biophysical Journal 106(1) 210–219

FIGURE 1 Pore structure and free energy profiles for pore formation. (A) Snapshots of the bilayers at equilibrium. (B) Bilayers with a hydrophilic pore. (Gray lines) Water; (silver lines) lipid tails; (blue lines) phosphates; and (red lines) cholines. (C) Partial density profiles for the pore (thick lines) and nonpore (thin lines) states. (Red) Water density; (black) phosphorous density. (D) Same as panel C, but the y axis is scaled. (E) Free energy profiles for moving a single lipid headgroup along the normal to the plane of the bilayer. The error bars are the difference in PMF from two independent runs for DLPC and DPPC, and the mean 5 SE from four independent PMFs for DMPC. To see this figure in color, go online.

phosphates (Fig. 1 C). The partial densities in the pore state show the phosphate density and water density increase in the bilayer interior, indicating the pore is hydrophilic. The amount of water and phosphates in the pore decreases as the bilayer becomes thicker; the pore is smaller in the thicker DPPC bilayer than DMPC and DLPC (Fig. 1 B). We determined how system properties change along our reaction coordinate, as a single lipid moves from equilibrium to the center of the DMPC bilayer (Fig. 2). The changes in the number of hydrogen bonds are plotted in Fig. 2 B. There is an average of 6186 5 0.3 hydrogen bonds in the whole system, and large fluctuations (from 6093 to 6284), but the average total number of hydrogen bonds remains constant along the reaction coordinate. As the single lipid headgroup is moved toward the bilayer center, the average number of water-water hydrogen bonds decreases by ~27 and the number of lipid-water hydrogen bonds increases by ~27. The change in hydrogen bonds coincides with an increase in solvent-accessible surface area as the pore forms (Fig. 2 C). The water-accessible surface area for the toroidal pore is greater than for the lamellar state, which allows more water-lipid hydrogen bonds to form, at the expense of water-water hydrogen bonds.

Pores in Lipid Bilayers TABLE 2

DLPC DMPC DPPC

213

Thermodynamic and kinetic data for membrane pores Thickness (nm)

DGpore (kJ/mol)

tclose (ns)

Water flux (ns1)a

Lipid flux (ns1)b

kopen (ns1)

r (pores/cm2)

P lipid (cm/s)c

P water (cm/s)c

1.45 1.58 1.86

17 45 78

þ500 123d 6.6

27.3 13.5 6.2

0.07 0.03 —e

3.6  106 4.3  1010 3.7  1014

2  1011 7.4  106 35.2

1.7  102 2.0  107 2.7  1012

0.2 3.0  106 6.6  1012

a

The average number of water molecules that move across the pore every 1 ns. The average number of lipid PC headgroups that crossed the pore every 1 ns. c The permeation coefficient for water and lipids in the bilayers through pores. d Estimated from the exponential fit from Fig. 5. e Two flips were observed in 40 simulations, which is not enough to calculate an average. b

To address the underlying driving forces for pore formation we decomposed the free energy profile into enthalpic and entropic contributions (Fig. 2 A). The enthalpic contribution was determined from the total system energy

(kinetic and potential) for each umbrella window, ignoring the PV component, because this is a negligible contribution. As the lipid is moved into the hydrophobic interior, a water defect forms to solvate the lipid; at the same time, there is a modest rise in the enthalpy of the system, and a relatively constant entropy (albeit with large error bars). The transition from initial water defect to pore formation causes a sudden increase in unfavorable entropy, and a favorable change in enthalpy of ~50 kJ/mol. This shows that pore formation in a DMPC bilayer is unfavorable due to entropy and is actually enthalpically favorable, compared to the lamellar state. Equilibrium pore formation

FIGURE 2 Structural and energetic properties during pore formation in a DMPC bilayer. (A) Free energy decomposition for DMPC pore formation: the change in enthalpy and entropy along with the free energy. We set all the curves equal to 0 at 1.7 nm from the bilayer center, which is the free energy trough. (B) The change in hydrogen bonds between different components of the system. We set all the curves equal to 0 at 1.7 nm from the bilayer center. (C) The average solvent-accessible surface area for the DMPC bilayer. For panels B and C, we plot the average of four independent systems and the resulting standard error of the mean. To see this figure in color, go online.

Starting from a lamellar DLPC bilayer, we ran 50 simulations for 500 ns and 50 simulations for 800 ns and observed four pore opening events. Fig. 3 shows the progression of pore formation for these events. In all cases, there is a prepore state that corresponds to a water wire across the membrane. After the water wire forms, lipid headgroups move into the bilayer interior, usually within 400 ps, and cause the pore to expand and stabilize. Once a stable pore forms it lasts for the rest of the simulation, as expected from the slow rate of pore closure (addressed below). Fig. 4 shows the progression of the system’s potential energy during pore formation and pore closure. The formation and closure results are two simulations out of six independent simulations for DLPC pore closing and four of pore opening, chosen to illustrate the processes. Pore opening is clearly indicated by an abrupt drop in potential energy and water-water hydrogen bonds and a gain in lipid-water hydrogen bonds. This process happens quickly with no long-lived transition states that are clearly observed. Pore closure basically mirrors pore opening and does not show long-lived transitions either, although we note that the x axis for formation and closure are on different timescales. Pore closure Starting from preformed pores we ran unrestrained simulations to determine the time required for pore dissipation in Biophysical Journal 106(1) 210–219

214

Bennett et al.

FIGURE 3 Mechanism of pore formation in four independent simulations of a DLPC bilayer. The color scheme is the same as in Fig 1. (Leftmost panels) The simulation time, and the change in time for the consecutive panels. To see this figure in color, go online.

the different bilayers. We defined pore closure as when the number of phosphates within the bilayer interior became zero for 1 ns. Occasionally, the number of phosphates within the bilayer interior can drop to zero, but the pore does not close. Table 2 lists the average time required for pore closure in the three model bilayers. For DLPC, we did not observe enough pore closure events to estimate the rate of closure, with only six closures out of 100 simulations of ~500 ns. All the 30 simulations of the DPPC pores resulted in pore closure with an average time of 6 ns. Out of 280 simulations we found only 105 pore closures for DMPC, which gave a closure time of 123 ns, from an exponential fit to the pore closure times (Fig. 5). We determined the number of water molecules and lipid molecules that cross the pore before closure (see Methods). For DLPC, on average 27.3 water molecules cross the pore in one direction per nanosecond, and 0.07 lipids. DMPC has roughly half the flux, with 13.5 water and 0.03 lipids per nanosecond. The pores close too quickly in the DPPC bilayer to obtain accurate statistics for lipid flux, but the rate for water flux is approximately half of the DMPC bilayer, with 6.2 waters per nanosecond. Biophysical Journal 106(1) 210–219

Table 2 lists the estimated rate for pore formation and the permeation of water and lipids through DLPC, DMPC, and DPPC pores. We predict orders-of-magnitude morefrequent pore formation in the thinner DLPC bilayer than the thicker DPPC bilayer.

DISCUSSION Pore formation and closure Pore formation occurs spontaneously when the phosphate of a single lipid is moved to the center of the bilayer. In a previous study on PC flip-flop, formation of pores was a corollary of the flip-flop process (14). In this work, we have used similar methods and extended simulations to investigate pore formation as well as pore closure. We find differences in both the structure and energetics of pore formation in DLPC, DMPC, and DPPC bilayers. Pore formation is much more unfavorable in the thicker bilayers, as expected from the longer hydrophobic tails. This agrees with recent simulations of arginine partitioning across saturated lipid bilayers with different tail lengths (47). The reduced free energy

Pores in Lipid Bilayers

FIGURE 4 Independent simulations of pore formation and closure in a DLPC bilayer. (Top panels) Total system energy; (bottom panels) change in the number of hydrogen bonds during the simulations. Both the energies and change in hydrogen bonds were set to 0 at the average value for the nonpore state. Note the timescale of the x axis is different for formation and dissipation. To see this figure in color, go online.

difference results in orders-of-magnitude more pores forming in the thinner DLPC bilayer than the thicker DPPC bilayer. Based on the low free energy cost for pore formation in the DLPC bilayer, we ran multiple long timescale simulations to directly observe pore opening in equilibrium simulations. We directly observed multiple spontaneous pore openings in the DLPC bilayer at 323 K that allowed us to examine the mechanism of pore formation in detail. A water

FIGURE 5 DMPC pore closure times and probabilities. (Red bars) Histograms for the number of DMPC pores that closed within that time. (Green curve) Probability the pore will still be open at a given time, or 1 – P(close), which was calculated by integrating the closure times’ normalized histogram. (Black curve) Exponential fit to the green curve. To see this figure in color, go online.

215

wire spanning the hydrophobic DLPC bilayer was observed immediately preceding pore formation. Pores formed quickly—within 2 ns—once the water wire was observed, although defining when a pore is formed is somewhat arbitrary. Similar to pore opening, where a water-wire precedes pore formation, closure occurs by the headgroups first leaving the bilayer interior, leaving a small number of waters within the hydrophobic interior, which then quickly diffuse out into bulk solution. The water-wire mechanism for pore formation and closure that we find is similar to the mechanism shown from previous atomistic simulations of electroporation in lipid bilayers (20,23,26). Our calculations provide valuable insight into the underlying energetics and mechanism of pore formation and closure. As expected, pore formation results in an increase in the solvent-accessible surface area of the bilayer. This creates more interactions and hydrogen bonds between the lipids and water, at the expense of water-water hydrogen bonds. For DLPC and DMPC, this results in a favorable change in enthalpy for the pore state compared to the nopore state. A possible explanation for this data is that lipid-water hydrogen bonds are more favorable than water-water hydrogen bonds. There is a strong tendency for the system to maintain the total number of hydrogen bonds. This is somewhat surprising given the radical structural rearrangement between the pore and no-pore states (Fig. 1), and the large changes in the average number of lipid-water and water-water hydrogen bonds (~27 for DMPC). PMF decompositions The source of the free energy cost for pore formation is not obvious, given the complex interactions at interfaces. We are unaware of a previous free energy decomposition for the process of pore formation in lipid bilayers using atomistic simulations. We found that for DMPC the free energy cost for pore formation is due to an unfavorable change in entropy. For DMPC, pore formation is favorable enthalpically compared to the lamellar bilayer state, which could be due to the increased surface area of the water-lipid headgroup interface and the increased number of water-lipid hydrogen bonds. A possible explanation for the change in entropy is that water loses translation and orientational freedom within the pore, due to the increased interactions with the membrane interphase in the pore state. It was shown that there is a favorable gain in entropy of 3.5 kJ/mol at 300 K to move a single water molecule from a DPPC interface to bulk water (48). This is also supported by the large literature on water confinement in nanotubes (for example, Pascal et al. (49)). One of the driving forces for lipid aggregation is minimizing the surface area formed with water. By increasing the surface area, moving from a flat bilayer surface, to a curved toroidal pore shape, we find an unfavorable entropic change. The strong dependence of lipid-water and water-water interactions on the Biophysical Journal 106(1) 210–219

216

energetics of pore formation, as well as the role of entropy is clearly demonstrated from our simulations. Rates Using the free energy cost for pore formation and the rate for pore closure, we are able to estimate the rate of pore opening in the different bilayers. We find a rate for pore formation in a small bilayer patch of DLPC of 4  106 ns1 (assuming a pore closure time of 500 ns, which is the length of our simulations and therefore the lower bound for the pore closure time). From the equilibrium simulations of pore formation in DLPC bilayers, we observed four openings in ~65 ms of simulation, which translates into a rate of 6  105 ns1. This order-of-magnitude difference between the predicted rate and observed rate is likely due to our assumed closure time of 500 ns, as well as the limited number of only four observed pore openings. For DMPC, we predict a rate of pore formation of 4  1010 ns1, and 4  1014 ns1 for DPPC in small bilayer patches. Pores will form orders-of-magnitude faster in thinner bilayers, in agreement with previous experimental work (50,51). From the free energy cost for pore formation, we calculate the pore density (see Methods). The permeation rate for small molecules that cross through a pore-mediated mechanism can then be calculated from pore density and flux of molecules across the pore (46). If this rate is much slower than the permeation rate via another mechanism, such as the solubility-diffusion mechanism, then the overall permeation rate will not depend strongly on the rate of pore formation. For example, water has been shown to permeate model membranes orders-of-magnitude faster than the rate for pore formation, whereas for some ions, such as chloride, the rate for bilayer permeation closely matches the rate for pore formation (16). The rate of water permeation we find for DMPC (3  106 cm/s) and DPPC (6  1012 cm/s) is orders-ofmagnitude slower than the experimental permeation rate, suggesting water crosses in a pore-independent mechanism. For DLPC, we find a permeation rate coefficient of 0.2 cm/s, which is faster than the rate found experimentally (10 103 cm/s) (52). It is possible that our model for DLPC underestimates the free energy cost for pore formation (discussed below). However, in Mathai et al. (52), it is noted that for DLPC the method for determining the water permeation rate using a carboxyfluorescein self-quenching assay could not be used, because the DLPC liposomes were too leaky. Future directions Obtaining adequate sampling is a concern for free energy calculations involving lipid bilayers (53). Because we do not observe pore closure in the umbrella simulations (the closing time is longer than the simulation time), which we would expect near the transition between water defect and pore, we might be underestimating the error in some of Biophysical Journal 106(1) 210–219

Bennett et al.

the PMFs. Our simulations of equilibrium pore formation in the DLPC bilayer qualitatively agree with our prediction of rapid pore formation from the low energy cost. Our results can aid the development of coarse-grained force fields and continuum membrane models that include pore formation. The atomistic PMF decomposition and the importance of entropy and hydrogen bonding are difficult to incorporate in CG simulations and likely difficult to model accurately other lattice-based and continuum computer simulations as well. We recently showed that the MARTINI model was deficient at modeling water defects and pore formation in lipid bilayers during lipid flip-flop (44). Whether these details can be incorporated into a future CG model remains to be seen, but might be necessary to accurately investigate many interesting membrane-deforming processes. There are a number of areas of research to pursue given the results, and limitations, of this work. Increased sampling with increased computer power will make many of these possible in the near future. We do not observe pore widening in any of our simulations, which would be important for processes that involve membrane rupture. The use of periodic boundary conditions prevents pore widening, and could also stabilize pores from closure. Using a method similar to Wohlert et al. (45) where the reaction coordinate is the pore radius, would allow one to calculate this energy. We suggest that using DLPC for a model bilayer for studies on a variety of pore-related problems, such as antimicrobial peptides, could be useful, given that we observe pore formation during equilibrium simulations. This means that complicated reaction coordinates, or highly nonequilibrium conditions, would not be needed to observe pore formation. Biological implications Pores form spontaneously in lipid bilayer systems, but how frequently this occurs is not well known, although critical for many processes. As well, the size, atomistic properties, and lifetime of the pores are important, but relatively uncharacterized. The overall structural rearrangements we observe in the porated membranes are similar to many simulations on electrochemical gradients, antimicrobial peptides, cellpenetrating peptides, and lateral stress within the plane of the membrane. A common structure suggests similar energetics—the propensity of a bilayer to form pores depends on the lipid properties, as well as the interaction with water. Fundamentally, these agents must lower the free energy cost for pore formation, either by further reducing the enthalpy or decreasing the unfavorable entropic cost. If the major cost is water getting stuck in the pore, then this presents interesting potential design ideas for future membrane porating schemes and suggests a novel mechanism for pore-forming peptides. For example, antimicrobial peptides and cell-penetrating peptides may decrease the number of waters in the pore, thus reducing the unfavorable entropic component and free energy cost. Alternatively, these peptides may bind

Pores in Lipid Bilayers

water before reaching the membrane, before pore formation, thereby reducing the change in bound water between pore and no-pore states, and again reducing the entropic penalty for pore formation. Understanding this balance may help designing specific peptides for enhanced pore formation. Lipids can move from one leaflet to the other, which is important for cellular membranes to create and maintain asymmetric lipid distributions and for the equal growth of bilayer leaflets. How cells maintain the proper distribution of lipids is not completely clear, although proteins are known to play a critical role in many aspects (11). From our calculations, we can estimate the rates for passive PC lipid translocation in model bilayers. Lipids can flip-flop across the thinner DLPC bilayer orders-of-magnitude faster than the DPPC bilayer. This is mostly due to the lower free energy cost for pore formation in the thinner bilayer, but also because the thinner bilayer pores are more stable and larger than the thicker bilayers. Experiments predict a wide range of estimates for the rate of lipid flip-flop depending on the experimental conditions and the assay (54–57). We showed previously that the rate of DPPC (58) and cholesterol (59) flip-flop in DPPC bilayers was much slower with higher concentrations of cholesterol. The presence of model hydrophobic peptides did not have a significant affect on the free energy barrier for DOPC flip-flop (60). It will be interesting to investigate the role of other membrane factors on pore formation in atomistic simulations, such as anionic headgroups, polyunsaturated tails, asymmetric leaflets, and near domain interfaces. Defect and pore formation in lipid bilayers and the collective behavior of lipids are important for interpreting the results from experiments and for understanding transmembrane protein structure and stability. A good example is reconciling the different amino-acid hydrophobicity scales derived from model peptide partitioning (61), the reversible insertion using the Sec translocon machinery (62), the folding and insertion of the OmpLA b-barrel protein (9), and computational free energies of partitioning (32). Although the scales provide qualitatively similar and correlated results, understanding the atomistic details of each experiment is needed to explain differences in the magnitudes of the free energies in each scale (63). For example, the OmpLA experiments were conducted using a DLPC vesicle, as this was the only lipid in which reversible folding of this outer membrane protein was observed (9). The higher likelihood of pore or defect formation in DLPC may have been essential for the reversible insertion of the protein due to the deformability of the thin bilayer. The physicochemical environment in the thin DLPC bilayer may be different than a thicker membrane, and the possibility of pore formation may affect the resulting free energies of transfer. Simulations of the OmpLA system helped reconcile the low free energy for the partitioning of arginine in the thin DLPC bilayer (64) compared to previous estimates for arginine on a polyleucine helix in thicker DPPC bilayers (65), using

217

the same force field. Short chain lipids were shown to increase the folding and aggregation of OpaI and OpaA b-barrel proteins (66). The insertion and translocation of peptides would also be affected by the bilayer thickness and its propensity to deform. The mechanism for antimicrobial pore formation was recently studied using DLPC and DLPG mixed bilayers (67). The relatively slow insertion of low pH-low insertion peptides into a POPC bilayer at acidic pH values and fast exit at basic conditions (68) may be related to the high cost for pore and defect formation of POPC (14). We predict that using DLPC or DMPC as a model membrane would drastically increase the rate of insertion. Because the timescales involved in defect formation in model PC membranes vary over orders of magnitude depending on membrane thickness, defect formation may be a key determinant in many macroscopically slow processes involving membrane proteins and membrane-active peptides.

CONCLUSIONS Atomistic resolution molecular dynamics simulations allowed us to observe and characterize pore formation in model membranes at an unprecedented level of detail. In agreement with previous results, we find a much reduced free energy cost for pore formation in thinner PC bilayers. Based on the free energy cost for pore formation and the rate of pore closure, we predict orders-of-magnitude more-frequent pore formation in the thinner DLPC bilayer than the DPPC bilayer. Pore formation is directly observed in thin DLPC bilayers agreeing qualitatively with previous experimental data and our rate estimates from the free energy cost. We found that pore formation is enthalpically favorable in a DMPC bilayer, and the free energy cost was due to a large unfavorable entropy. The total number of hydrogen bonds in the system was maintained during pore formation, although the number of water-water hydrogen bonds decreased at the expense of increased lipid-water hydrogen bonds. We speculate that the entropic component may be mechanistically important for a large number of applications, including antimicrobial peptides and electroporation. Membrane processes involving macroscopically long timescales that depend strongly on lipid composition likely involve defects and could be explained by the probability to form defects and pores. These results present what we believe to be a novel view of pore formation in lipid bilayers, and open many avenues for further research. Calculations were carried out on WestGrid/Compute Canada facilities. This work was supported by the Natural Sciences and Engineering Research Council. W.F.D.B. is supported by a Banting Postdoctoral Fellowship through Natural Sciences and Engineering Research Council. D.P.T. is an Alberta Innovates Health Solutions Scientist and Alberta Innovates Technology Futures Strategic Chair in (Bio)Molecular Simulation. Biophysical Journal 106(1) 210–219

218

REFERENCES 1. Epand, R. M., and H. J. Vogel. 1999. Diversity of antimicrobial peptides and their mechanisms of action. Biochim. Biophys. Acta. 1462:11–28.

Bennett et al. 23. Levine, Z. A., and P. T. Vernier. 2010. Life cycle of an electropore: field-dependent and field-independent steps in pore creation and annihilation. J. Membr. Biol. 236:27–36. 24. Delemotte, L., and M. Tarek. 2012. Molecular dynamics simulations of lipid membrane electroporation. J. Membr. Biol. 245:531–543.

2. Nguyen, L. T., E. F. Haney, and H. J. Vogel. 2011. The expanding scope of antimicrobial peptide structures and their modes of action. Trends Biotechnol. 29:464–472.

25. Gurtovenko, A. A., and J. Anwar. 2007. Ion transport through chemically induced pores in protein-free phospholipid membranes. J. Phys. Chem. B. 111:13379–13382.

3. Brogden, K. A. 2005. Antimicrobial peptides: pore formers or metabolic inhibitors in bacteria? Nat. Rev. Microbiol. 3:238–250.

26. Tieleman, D. P. 2004. The molecular basis of electroporation. BMC Biochem. 5:10.

4. Vernier, P. T., Y. Sun, ., M. A. Gundersen. 2004. Nanoelectropulseinduced phosphatidylserine translocation. Biophys. J. 86:4040–4048. 5. Zorko, M., and U. Langel. 2005. Cell-penetrating peptides: mechanism and kinetics of cargo delivery. Adv. Drug Deliv. Rev. 57:529–545. 6. Last, N. B., and A. D. Miranker. 2013. Common mechanism unites membrane poration by amyloid and antimicrobial peptides. Proc. Natl. Acad. Sci. USA. 110:6382–6387. 7. Bergstrom, C. L., P. A. Beales, ., J. T. Groves. 2013. Cytochrome c causes pore formation in cardiolipin-containing membranes. Proc. Natl. Acad. Sci. USA. 110:6269–6274. 8. MacCallum, J. L., W. F. D. Bennett, and D. P. Tieleman. 2011. Transfer of arginine into lipid bilayers is nonadditive. Biophys. J. 101:110–117. 9. Moon, C. P., and K. G. Fleming. 2011. Side-chain hydrophobicity scale derived from transmembrane protein folding into lipid bilayers. Proc. Natl. Acad. Sci. USA. 108:10174–10177.

27. Herce, H. D., and A. E. Garcia. 2007. Molecular dynamics simulations suggest a mechanism for translocation of the HIV-1 TAT peptide across lipid membranes. Proc. Natl. Acad. Sci. USA. 104:20805–20810. 28. Yesylevskyy, S., S. J. Marrink, and A. E. Mark. 2009. Alternative mechanisms for the interaction of the cell-penetrating peptides penetratin and the TAT peptide with lipid bilayers. Biophys. J. 97:40–49. 29. Huang, K., and A. E. Garcı´a. 2013. Free energy of translocating an arginine-rich cell-penetrating peptide across a lipid bilayer suggests pore formation. Biophys. J. 104:412–420. 30. Leontiadou, H., A. E. Mark, and S. J. Marrink. 2006. Antimicrobial peptides in action. J. Am. Chem. Soc. 128:12156–12161. 31. Leontiadou, H., A. E. Mark, and S. J. Marrink. 2004. Molecular dynamics simulations of hydrophilic pores in lipid bilayers. Biophys. J. 86:2156–2164.

10. White, S. H., and G. von Heijne. 2008. How translocons select transmembrane helices. Annu. Rev. Biophys. 37:23–42.

32. MacCallum, J. L., W. F. D. Bennett, and D. P. Tieleman. 2008. Distribution of amino acids in a lipid bilayer from computer simulations. Biophys. J. 94:3393–3404.

11. van Meer, G., D. R. Voelker, and G. W. Feigenson. 2008. Membrane lipids: where they are and how they behave. Nat. Rev. Mol. Cell Biol. 9:112–124.

33. Gleason, N. J., V. V. Vostrikov, ., R. E. Koeppe, 2nd. 2013. Buried lysine, but not arginine, titrates and alters transmembrane helix tilt. Proc. Natl. Acad. Sci. USA. 110:1692–1695.

12. Vernier, P. T., M. J. Ziegler, ., D. P. Tieleman. 2006. Nanopore formation and phosphatidylserine externalization in a phospholipid bilayer at high transmembrane potential. J. Am. Chem. Soc. 128:6288–6289.

34. van der Spoel, D., E. Lindahl, ., H. J. Berendsen. 2005. GROMACS: fast, flexible, and free. J. Comput. Chem. 26:1701–1718.

13. Hoskin, D. W., and A. Ramamoorthy. 2008. Studies on anticancer activities of antimicrobial peptides. Biochim. Biophys. Acta. 1778:357–375. 14. Sapay, N., W. F. D. Bennett, and D. P. Tieleman. 2009. Thermodynamics of flip-flop and desorption for a systematic series of phosphatidylcholine lipids. Soft Matter. 5:3295–3302. 15. Kornberg, R. D., and H. M. McConnell. 1971. Inside-outside transitions of phospholipids in vesicle membranes. Biochemistry. 10:1111–1120. 16. Toyoshima, Y., and T. E. Thompson. 1975. Chloride flux in bilayer membranes: chloride permeability in aqueous dispersions of singlewalled, bilayer vesicles. Biochemistry. 14:1525–1531.

35. Hess, B., C. Kutzner, ., E. Lindahl. 2008. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 4:435–447. 36. Berger, O., O. Edholm, and F. Ja¨hnig. 1997. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 72:2002–2013. 37. Berendsen, H. J. C., J. P. M. Postma, ., J. Hermans. 1981 Interaction models for water in relation to protein hydration. In Book Interaction Models for Water in Relation to Protein Hydration. B. Pullman, editor. D. Reidel, Dordrecht, The Netherlands. 38. Essmann, U., L. Perera, ., L. G. Pedersen. 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103:8577–8593.

17. Missner, A., and P. Pohl. 2009. 110 years of the Meyer-Overton rule: predicting membrane permeability of gases and other small compounds. ChemPhysChem. 10:1405–1414.

39. Darden, T., D. York, and L. Pedersen. 1993. Particle mesh Ewald—an N.Log(N) method for Ewald sums in large systems. J. Chem. Phys. 98:10089–10092.

18. Gurtovenko, A. A., J. Anwar, and I. Vattulainen. 2010. Defect-mediated trafficking across cell membranes: insights from in silico modeling. Chem. Rev. 110:6077–6103.

40. Berendsen, H. J. C., J. P. M. Postma, ., J. R. Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684–3690.

19. Marrink, S. J., A. H. de Vries, and D. P. Tieleman. 2009. Lipids on the move: simulations of membrane pores, domains, stalks and curves. Biochim. Biophys. Acta. 1788:149–168.

41. Miyamoto, S., and P. A. Kollman. 1992. SETTLE—an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 13:952–962.

20. Bo¨ckmann, R. A., B. L. de Groot, ., H. Grubmu¨ller. 2008. Kinetics, statistics, and energetics of lipid membrane electroporation studied by molecular dynamics simulations. Biophys. J. 95:1837–1850.

42. Hess, B., H. Bekker, ., J. G. E. M. Fraaije. 1997. LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 18:1463–1472.

21. Tieleman, D. P., H. Leontiadou, ., S. J. Marrink. 2003. Simulation of pore formation in lipid bilayers by mechanical stress and electric fields. J. Am. Chem. Soc. 125:6382–6383.

43. Kumar, S., D. Bouzida, ., J. M. Rosenberg. 1992. The weighted histogram analysis method for free-energy calculations on biomolecules. 1. The method. J. Comput. Chem. 13:1011–1021.

22. Vernier, P. T., M. J. Ziegler, ., D. P. Tieleman. 2006. Nanopore-facilitated, voltage-driven phosphatidylserine translocation in lipid bilayers—in cells and in silico. Phys. Biol. 3:233–247.

44. Bennett, W. F. D., and D. P. Tieleman. 2011. Water defect and pore formation in atomistic and coarse-grained lipid membranes: pushing the limits of coarse graining. J. Chem. Theory Comput. 7:2981–2988.

Biophysical Journal 106(1) 210–219

Pores in Lipid Bilayers

219

45. Wohlert, J., W. K. den Otter, ., W. J. Briels. 2006. Free energy of a trans-membrane pore calculated from atomistic molecular dynamics simulations. J. Chem. Phys. 124:154905.

57. Moreno, M. J., L. M. B. B. Estronca, and W. L. C. Vaz. 2006. Translocation of phospholipids and dithionite permeability in liquid-ordered and liquid-disordered membranes. Biophys. J. 91:873–881.

46. Tieleman, D. P., and S. J. Marrink. 2006. Lipids out of equilibrium: energetics of desorption and pore mediated flip-flop. J. Am. Chem. Soc. 128:12462–12467.

58. Bennett, W. F. D., J. L. MacCallum, and D. P. Tieleman. 2009. Thermodynamic analysis of the effect of cholesterol on dipalmitoylphosphatidylcholine lipid membranes. J. Am. Chem. Soc. 131:1972–1978.

47. Li, L. B., I. Vorobyov, and T. W. Allen. 2012. The role of membrane thickness in charged protein-lipid interactions. Biochim. Biophys. Acta. 1818:135–145.

59. Bennett, W. F. D., J. L. MacCallum, ., D. P. Tieleman. 2009. Molecular view of cholesterol flip-flop and chemical potential in different membrane environments. J. Am. Chem. Soc. 131:12714–12720.

48. Debnath, A., B. Mukherjee, ., S.-T. Lin. 2010. Entropy and dynamics of water in hydration layers of a bilayer. J. Chem. Phys. 133:174704.

60. Sapay, N., W. F. D. Bennett, and D. P. Tieleman. 2010. Molecular simulations of lipid flip-flop in the presence of model transmembrane helices. Biochemistry. 49:7665–7673.

49. Pascal, T. A., W. A. Goddard, and Y. Jung. 2011. Entropy and the driving force for the filling of carbon nanotubes with water. Proc. Natl. Acad. Sci. USA. 108:11794–11798. 50. Blok, M. C., E. C. M. van der Neut-Kok, ., J. de Gier. 1975. The effect of chain length and lipid phase transitions on the selective permeability properties of liposomes. Biochim. Biophys. Acta. 406:187–196. 51. Paula, S., A. G. Volkov, ., D. W. Deamer. 1996. Permeation of protons, potassium ions, and small polar molecules through phospholipid bilayers as a function of membrane thickness. Biophys. J. 70:339–348.

61. Wimley, W. C., and S. H. White. 1996. Experimentally determined hydrophobicity scale for proteins at membrane interfaces. Nat. Struct. Biol. 3:842–848. 62. Hessa, T., H. Kim, ., G. von Heijne. 2005. Recognition of transmembrane helices by the endoplasmic reticulum translocon. Nature. 433:377–381. 63. MacCallum, J. L., and D. P. Tieleman. 2011. Hydrophobicity scales: a thermodynamic looking glass into lipid-protein interactions. Trends Biochem. Sci. 36:653–662.

52. Mathai, J. C., S. Tristram-Nagle, ., M. L. Zeidel. 2008. Structural determinants of water permeability through the lipid membrane. J. Gen. Physiol. 131:69–76.

64. Gumbart, J., and B. Roux. 2012. Determination of membrane-insertion free energies by molecular dynamics simulations. Biophys. J. 102:795–801.

53. Neale, C., W. F. D. Bennett, ., R. Pomes. 2011. Statistical convergence of equilibrium properties in simulations of molecular solutes embedded in lipid bilayers. J. Chem. Theory Comput. 7:4175–4188.

65. Dorairaj, S., and T. W. Allen. 2007. On the thermodynamic stability of a charged arginine side chain in a transmembrane helix. Proc. Natl. Acad. Sci. USA. 104:4943–4948.

54. Liu, J., and J. C. Conboy. 2005. 1,2-diacyl-phosphatidylcholine flipflop measured directly by sum-frequency vibrational spectroscopy. Biophys. J. 89:2522–2532.

66. Dewald, A. H., J. C. Hodges, and L. Columbus. 2011. Physical determinants of b-barrel membrane protein folding in lipid vesicles. Biophys. J. 100:2131–2140.

55. Wimley, W. C., and T. E. Thompson. 1991. Transbilayer and interbilayer phospholipid exchange in dimyristoylphosphatidylcholine/dimyristoylphosphatidylethanolamine large unilamellar vesicles. Biochemistry. 30:1702–1709.

67. Rakowska, P. D., H. Jiang, ., M. G. Ryadnov. 2013. Nanoscale imaging reveals laterally expanding antimicrobial pores in lipid bilayers. Proc. Natl. Acad. Sci. USA. 110:8918–8923.

56. Nakano, M., M. Fukuda, ., T. Handa. 2009. Flip-flop of phospholipids in vesicles: kinetic analysis with time-resolved small-angle neutron scattering. J. Phys. Chem. B. 113:6745–6748.

68. Andreev, O. A., A. G. Karabadzhak, ., Y. K. Reshetnyak. 2010. pH (low) insertion peptide (pHLIP) inserts across a lipid bilayer as a helix and exits by a different path. Proc. Natl. Acad. Sci. USA. 107:4081–4086.

Biophysical Journal 106(1) 210–219