The free energy of pore formation in lipid bilayers has been previously calculated using a variety of reaction coordinates. Here, we use free energy perturbation of a cylindrical lipid exclusion restraint to compute the free energy profile as a function of pore radius in dimyristoylphosphatidylcholine (DMPC) and dioleoylphosphatidylcholine (DOPC) bilayers. Additionally restraining the headgroups to lie on the membrane surface allows us to also calculate the free energy profile of hydrophobic pores, i.e., cylindrical pores lined by acyl chains. For certain pore radii, the free energy of wetting of hydrophobic pores is calculated using the density bias method. It is found that wetting of hydrophobic pores becomes thermodynamically favorable at 5.0 Å for DMPC and 6.5 Å for DOPC, although significant barriers prevent spontaneous wetting of the latter on a nanosecond time scale. The free energy of transformation of hydrophilic pores to hydrophobic ones is also calculated using free energy perturbation of headgroup restraints along the bilayer normal. This quantity, along with wetting and pore growth free energies, provides complete free energy profiles as a function of radius. Pore line tension values for the hydrophilic pores obtained from the slope of the free energy profiles are 37.6 pN for DMPC and 53.7 pN for DOPC. The free energy profiles for the hydrophobic pores are analyzed in terms of elementary interfacial tensions. It is found that a positive three-phase line tension is required to explain the results. The estimated value for this three-phase line tension (51.2 pN) lies within the expected range.

Biological membranes self-assemble, thanks to the amphiphilic nature of their lipid constituents. They serve as a barrier between the interior and the exterior of a cell and as intracellular compartment separators. Unregulated formation of pores in biological membranes can be detrimental to a cell. This is a strategy employed by many pathogens1 and by the innate immune system against invaders.2 On the other hand, membrane pores need to be harmlessly induced in a variety of technological applications such as drug delivery, genetic transformation, sensing, and nanoelectronics.3,4 Thus, an understanding of the mechanism and knowledge of the quantitative thermodynamics of this process are essential for both fundamental and practical reasons.

Pore formation in lipid membranes is commonly induced experimentally by electric fields.5–7 In this process, termed electroporation, short pulses of high voltage cause a temporary breakdown of the cell membrane followed by complete resealing. Under certain conditions, however, the breakdown can be irreversible.8 The experimental results have been interpreted based on the formation of hydrophobic and hydrophilic defects,5 but the exact mechanism is still being investigated.9 Pores can also be induced by increasing membrane tension, either by osmotic pressure10 or by micropipette aspiration.11,12 In the former experiment, an osmotic pressure difference is created by diluting the salt concentration outside the lipid vesicles. In the latter experiment, a micropipette is used to exert suction on giant unilamellar vesicles until the vesicle breaks. Normally, the pores are short-lived, but applying tension to electroporated liposomes resulted in long-lived metastable pores.13 Local intermediate structures (nonconductive or hydrophobic prepores) before the formation of pores have been inferred in various experiments.12,14,15

The classical nucleation model describes the free energy of a circular pore of radius R in a two-dimensional surface as a competition between an unfavorable edge energy term and a favorable surface tension term,10,16

ΔGR> =2π >R> γπR2σ.
(1)

The first term in Eq. (1) is the free energy cost of having an edge of length 2πR, the circumference of a circular pore, and γ is the pore line tension. The second term is the gain in free energy from reducing the membrane surface area by πR2, σ being the surface tension. Equation (1) exhibits a maximum at R* = γ/σ. If the pore radius is less than this critical radius, the pore tends to close; otherwise, the pore grows in size indefinitely causing membrane rupture. To model electroporation, Eq. (1) is supplemented with a term proportional to the square of the applied voltage, representing the change in the energy of a capacitor when the low-dielectric material is replaced by the high-dielectric solvent.5 Equation (1) cannot be entirely valid down to zero radius, because a zero radius hydrophilic pore (i.e., a toroidal pore where the lipid headgroups meet at the center of the membrane) differs from an intact membrane. More detailed models predicted a metastable “nascent” hydrophilic pore, separated by a barrier from the intact membrane.17 

The above model is phenomenological, hiding all microscopic details into the line tension γ, assumed to be independent of R. Attempts to predict the line tension from the molecular properties of the lipids have been made based on continuum or mean field theories.18–24 For example, the line tension has been calculated from an extended “opposing forces” model that includes the contribution from the conformational free energy of the hydrocarbon chains.18 A lattice self-consistent field theory has provided the edge energy as a function of thickness and bending modulus of the bilayer and the spontaneous curvature of the monolayer.23 Self-consistent field theory has also been used to examine how the constitution of a lipid component affects the line tension.21,22 Ting et al.20 combined the dynamic self-consistent field theory with the string model to calculate the minimum energy path for pore formation and rupture in the lipid membrane. Finally, Akimov et al.24 have recently estimated the pore energy as a function of hydrophobic thickness of the lipid membrane at different pore radii by using continuum elasticity theory.

More molecular details have been obtained by molecular dynamics or Monte Carlo simulations of coarse-grained lipid representations in which the lipids typically consist of 3–5 particles, and solvent is taken into account either implicitly25–27 or explicitly.28–30 In these studies, pores may appear spontaneously25,29 or may be induced by stretching26,28 or other appropriate biases.27,30 These studies provided a testing ground of analytical theories and offered valuable insights into the relationship between pore formation energetics and the molecular characteristics of the lipids. The ability of more detailed coarse-grained lipid representations to provide a realistic view of membrane pores has also been examined.31 

Other studies employed an atomistic representation of lipids and water. The pore was induced by applying tension32,33 or electric field32,34–37 or by forcing a lipid headgroup into the membrane center.31,38,39 Efforts have been made to calculate the free energy of pore formation using a variety of reaction coordinates. One is a collective radial coordinate proposed by Tolpekina et al.,28 first applied to a coarse-grained model and then to an atomistic model.40 The free energy profile of the pore was found to be quadratic for small and linear for large pore radii. Two other reaction coordinates used were the distance between the single phosphate group and the bilayer center38 and the average water density inside a cylinder parallel to the bilayer normal.41 The free energy was obtained by either the mean constraint force approach28 or umbrella sampling.42 A comparative study of the above three reaction coordinates found that they exhibit hysteresis, to different extents, in the free energy profiles of pore opening and pore closing processes due to slow relaxation along the degrees of freedom that are orthogonal to the reaction coordinate.43 The collective radial coordinate perturbed the membrane more strongly than was required for pore formation and gave much higher pore-free energy than the other two coordinates. In all three reaction coordinates, the free energy increased monotonically, and no metastable pore was observed. More recently, Hub and Awasthi44 defined a new reaction coordinate to overcome the issue of hysteresis and observed a shallow metastable prepore state (a local minimum in the free energy profile) under certain conditions.

One limitation of the above studies is that they apply best to pore generation, rather than pore growth. The pore radius is not a natural coordinate in these methods but at most a parameter chosen at the outset. In the present work, we use an alternative approach to compute the free energy of a pore as a function of pore radius using a harmonic restraint that excludes all lipid atoms from a cylinder of a given radius. The free energy perturbation method45 is then used to compute the free energy change upon gradual increase in the cylinder radius in dimyristoylphosphatidylcholine (DMPC) and dioleoylphosphatidylcholine (DOPC) lipid bilayers. By employing additional restraints on the headgroups, we also estimate the free energy profile of hydrophobic pores, i.e., pores lined by lipid acyl chains. Additionally, to obtain complete free energy profiles, we calculate the free energy of wetting of hydrophobic pores using the density bias method41 and the free energy of converting a hydrophilic pore to a hydrophobic one using free energy perturbation of a planar exclusion restraint on the lipid headgroups. This approach provides data on the growth of hydrophilic pores, which can be used to obtain the pore line tension, and the free energy of idealized hydrophobic pores, which can be used to obtain fundamental interfacial tensions. In addition to providing fundamental insights, the data could be useful for building simplified models for peptide-decorated membrane pores.

In this work, we apply to all lipid atoms an asymmetric harmonic restraint that excludes them from a cylinder of radius R parallel to the membrane normal (z-axis),

UR=i=1NURi,URi=kriR2 >when> ri<R,> URi=0 >when> ri>R,
(2)

where ri is the distance of the lipid atom i from the z-axis (center of the pore), N is the total number of lipid atoms, and k is the force constant.

The free energy difference for increasing the pore radius from R to R + ΔR is computed by the free energy perturbation formula,45 

ΔGRR+ΔR=kBT> >lnexpUR+ΔRUR/kBTR,
(3)

where kB is Boltzmann’s constant and T is the absolute temperature. The average is computed with the system simulated at a cylinder radius R. We have performed the free energy perturbation simulations for four values of the force constant (k): 5.0 kcal mol−1 Å−2, 10.0 kcal mol−1 Å−2, 20.0 kcal mol−1 Å−2, and 30.0 kcal mol−1 Å−2. For each value, we computed the free energy for increasing the radius from 5.0 Å to 5.5 Å. The calculated values are similar (Table S1 in the supplementary material), so we chose k = 10 kcal mol−1 Å−2 for the remaining simulations. The radial increment ΔR should be the largest possible that gives accurate results. We found that a ΔR of 1 Å gives a result too high relative to two steps of 0.5 Å, whereas a ΔR of 0.5 Å gives a result similar to two steps of 0.25 Å (Table S2 in the supplementary material). Thus, we chose ΔR = 0.5 Å. The miscellaneous mean field potential (MMFP) module of the CHARMM program was used to implement this geometrical restraint. The cylinder radius R varied from 0 Å to 10.0 Å in 0.5-Å increments. The free energy changes were calculated using the PERT module of CHARMM,46 which allows the MMFP restraint energies to be included in the perturbation calculation. For each of the 20 windows, we performed 40 ns of simulations for DMPC and 48 ns for DOPC. 10 ns and 14 ns, respectively, were discarded as equilibration, and the remainder was used to calculate averages. Statistical uncertainties were estimated as standard deviations of block averages.

It is important to note that Eq. (3) can only be applied in the direction of increasing the pore radius. This is best-appreciated in the case of a hard sphere potential: performing a simulation with UR+ΔR pushes lipids out of a cylinder of radius R + ΔR, which of course also excludes them from a cylinder of radius R. Thus, the calculated free energy from R + ΔR to R evaluated at R + ΔR would be zero, whereas that from R to R + ΔR would be finite and positive. The situation is analogous to computing the chemical potential of a hard sphere solute. Applying the perturbation formula in the direction of annihilating a hard sphere would give a zero free energy change. This issue has been discussed in the literature, and possible remedies have been proposed.47–49 For soft potentials, this asymmetry is less sharp, but the bias still exists, leading to the realization that creation of a solute is always better than annihilation.50 For this reason, doing the usual “double wide sampling” (conducting a simulation at R and perturbing to R − ΔR and R + ΔR) would give erroneous results. Therefore, we do all perturbations in the direction of increasing R.

To generate hydrophobic pores, i.e., pores lined by acyl chains, we restrained the lipid headgroups to lie close to the bilayer surface. The average z distance between lipid P atoms in the upper and lower leaflets is 37.0 Å ± 4 Å (DMPC) and 39.5 Å ± 4 Å (DOPC), respectively (Fig. S1 in the supplementary material). Thus, the phosphates were restrained to lie at z = ±18.5 Å for DMPC and z = ±19.75 for DOPC, with an allowance of 4 Å. Beyond that, a harmonic force was applied with a force constant (kpo4) of 1.00 kcal mol−1 Å−2. Again, the MMFP module of CHARMM was used for these restraints. For both hydrophobic and hydrophilic pores, the system at each value of R was equilibrated for 15 ns.

To calculate the free energy of converting a hydrophilic, toroidal pore to a wet hydrophobic pore in DMPC, we also computed the free energy profiles of moving the headgroups from the center of the pore to the surface of the lipid bilayer by using the free energy perturbation of a planar restraint that excludes these groups from a slab of certain thickness. This calculation was done in increments of 1.0 Å from z = 0 Å to 11 Å for DMPC and z = 0 Å to 14 Å for DOPC. After that, the free energy of changing from the slab exclusion restraint to the planar position restraint that was used in the simulations of the hydrophobic pores was calculated. This calculation was done for R = 3 Å in DMPC and R = 5 Å in DOPC and included the phosphate atoms within 8 Å and 10 Å from the pore axis, respectively. A 5-ns simulation was carried out at each z value. Statistical uncertainties were estimated by splitting the sample into two blocks. The MMFP and PERT modules of CHARMM were also used for this calculation.

The molecular structures of DMPC and DOPC, a schematic picture of a lipid bilayer showing the total thickness, hydrophobic thickness, and headgroup thickness, and snapshots of a hydrophilic pore and a hydrophobic pore are shown in Fig. 1. Hydrophilic pores are toroidal, lined by lipid headgroups, whereas hydrophobic pores are cylindrical, lined by lipid acyl chains. They can be “wet,” as in Fig. 1, or they can be “dry” if water does not enter the pore.

FIG. 1.

(a) The molecular structure of DMPC and DOPC lipids. Blue, red, yellow, and gray colors are for nitrogen, oxygen, phosphorus, and carbon atoms, respectively. (b) Schematic picture of a lipid bilayer showing the total thickness, headgroup thickness, and hydrophobic thickness. [(c) and (d)] Snapshots of a hydrophilic pore (c) and a wet hydrophobic pore (d). Water is shown as red lines and lipid headgroup phospate atoms as blue spheres. The rest of the lipids is omitted for clarity.

FIG. 1.

(a) The molecular structure of DMPC and DOPC lipids. Blue, red, yellow, and gray colors are for nitrogen, oxygen, phosphorus, and carbon atoms, respectively. (b) Schematic picture of a lipid bilayer showing the total thickness, headgroup thickness, and hydrophobic thickness. [(c) and (d)] Snapshots of a hydrophilic pore (c) and a wet hydrophobic pore (d). Water is shown as red lines and lipid headgroup phospate atoms as blue spheres. The rest of the lipids is omitted for clarity.

Close modal

The density bias method41 with umbrella sampling42 was used to obtain free energy profiles between the dry and the wet state of hydrophobic pores for selected pore radii. These calculations employed the same cylindrical restraints as the free energy simulations. The water atom density within a cylinder covering the pore was controlled and restrained at specific values using a force constant of 1 kcal mol−1 nm−6. The cylinder extended from z = −13.5 Å to +13.5 Å for DMPC and from −14.5 to +14.5 for DOPC and had a radius equal to the pore radius. The water density was varied in windows spaced 3 nm−3 apart, and each window was simulated for 1 ns. In these calculations, the radial switching distance (RW) had to be large (10 Å or 5 Å); otherwise, the pore radius would increase as the algorithm pushed water radially out of the designated cylinder. Potentials of mean force were computed by the weighted histogram analysis method51 using primarily Grossfield’s code,52 double-checking the results with other programs.53,54 With Grossfield’s code, it was necessary to use a large number of bins (100 or more) to obtain reliable results. Statistical uncertainties were estimated by a Monte Carlo bootstrap method,52 but the noise in the profile is sometimes larger than these estimates.

Classical molecular dynamics simulations and free energy calculations were performed using the CHARMM program.46 Initial structures for the dimyristoylphosphatidylcholine (DMPC) and dioleoylphosphatidylcholine (DOPC) lipid bilayers were generated using the membrane builder55 module of the CHARMM–GUI web server.56 The CHARMM36 force field57 was used for the lipids. 140 DMPC molecules were hydrated with 5920 TIP3P58 water molecules in a 66 Å × 66 Å × 83 Å cell. In the case of DOPC, 140 lipid molecules were hydrated with 7084 TIP3P water molecules in a 70 Å × 70 Å × 90 Å cell. The leapfrog integrator59 with a 2-fs time step was used to integrate the equations of motion. The particle mesh Ewald method60 was used to compute the long range interactions with a cutoff of 10 Å, a k value of 0.34, an order of B-spline interpolation of 6, and a grid spacing of 1 Å. For the Lennard-Jones interactions, a switch function was used in the range 10 Å–16 Å. The temperature of the system was maintained at 303 K by using the Hoover thermostat.61 The pressure was maintained at 1 bar with the Langevin piston method,62 with piston mass of ∼4000 amu and collision frequency of 20 ps−1. Bonds involving hydrogen atoms were constrained using SHAKE.63 

The free energy profiles of growing a pore in a DMPC bilayer with and without headgroup restraints are shown in Fig. 2. Below R = 3.0 Å, the pore is hydrophobic in both cases, and the free energy curves are essentially identical and quadratic up to about 1 Å. Without headgroup restraints, a rearrangement occurs at R = 3.0 Å. At that point, water enters the pore [Fig. 3(b)], dragging along the lipid headgroups to create a toroidal, hydrophilic pore [Figs. 3(c) and 3(d)]. The pore remains fully hydrated at all larger radii, and the slope of the curve decreases.

FIG. 2.

The free energy profile as a function of pore radius for a DMPC lipid bilayer. In the red curve, the lipid headgroups are restrained to lie near one of the two membrane surfaces, whereas in the blue and black curves, no such restraints are applied. All pores along the red and blue curves are hydrophobic, whereas those along the black curve are hydrophilic.

FIG. 2.

The free energy profile as a function of pore radius for a DMPC lipid bilayer. In the red curve, the lipid headgroups are restrained to lie near one of the two membrane surfaces, whereas in the blue and black curves, no such restraints are applied. All pores along the red and blue curves are hydrophobic, whereas those along the black curve are hydrophilic.

Close modal
FIG. 3.

Snapshots of the simulation of DMPC at a pore radius of 3.0 Å without restraints on the headgroups at different times in the simulation. (a) 8 ns, (b) 8.345 ns, (c) 10 ns, and (d) 30 ns.

FIG. 3.

Snapshots of the simulation of DMPC at a pore radius of 3.0 Å without restraints on the headgroups at different times in the simulation. (a) 8 ns, (b) 8.345 ns, (c) 10 ns, and (d) 30 ns.

Close modal

This transformation of the pore into a wet, toroidal pore incurs a drop in free energy that involves degrees of freedom unrelated to the exclusion cylinder radius. This free energy must be calculated using a different approach. We split the process into two steps: wetting of a hydrophobic pore and relaxation of the headgroups from a wet hydrophobic pore to a wet toroidal pore. The free energy of the first step was calculated by the density bias method41 as +12.3 kcal/mol ± 0.3 kcal/mol (see Sec. III C). The free energy of the second step was calculated by free energy perturbation of restraints on the headgroups starting from a toroidal pore and gradually excluding lipid headgroups from a slab of increasing thickness. The free energy of this step was calculated as 21.7 kcal/mol ± 2 kcal/mol, or in the opposite direction as −21.7 kcal/mol ± 2 kcal/mol. The sum of the two steps gives ΔG = −9.4 kcal/mol ± 2 kcal/mol. Thus, a drop in this magnitude is incorporated at R = 3 Å between the red and the black curves of Fig. 2.

Taking the 3-Å hydrophilic pore and removing the radial restraints for 4 ns, we form hydrophilic pores of smaller radii. By applying then free energy perturbation from R = 0 Å to 3 Å in 0.5-Å increments, as described above, we can extend the wet hydrophilic pore curve down to 0 Å. The free energy of the zero radius (nascent) hydrophilic pore was calculated as 6.8 kcal/mol.

With restraints to keep headgroups on the membrane surface, the pore remains dry until R = 7 Å, and the free energy rises steeply. Pore wetting starts after 12.5 ns in the R = 7.0 Å simulation (Fig. 4) and causes a decrease in the slope of the free energy profile. In Sec. III C, we compute the free energy of wetting the DMPC hydrophobic pore at 7 Å as −12.4 kcal/mol ± 0.3 kcal/mol. This drop in free energy is included in the red curve in Fig. 2.

FIG. 4.

Snapshots of the simulation of DMPC at a pore radius of 7.0 Å restraining headgroups at the membrane surface.

FIG. 4.

Snapshots of the simulation of DMPC at a pore radius of 7.0 Å restraining headgroups at the membrane surface.

Close modal

The free energy profiles of pore growth in DOPC with and without headgroup restraints are given in Fig. 5. Without headgroup restraints, pore wetting starts after 12 ns in the R = 5.0-Å simulation (Fig. 6), and a fully toroidal pore is formed at R = 5.5 Å (Fig. S2). As in DMPC, the slope of the curve after formation of the toroidal pore decreases. The free energy of wetting a 5-Å DOPC pore was calculated to be +16.3 kcal/mol ± 0.2 kcal/mol, and the free energy of relaxing the headgroups was calculated to be −17 kcal/mol ± 1 kcal/mol, smaller than that for DMPC. Combining the two numbers, we obtained a drop of only about −0.7 kcal/mol ± 1 kcal/mol from the red hydrophobic pore curve to the black hydrophilic pore curve in Fig. 5. We have not extended the hydrophilic pore curve down to zero radius, but linear extrapolation gives ∼25 kcal/mol for a nascent hydrophilic pore at zero radius.

FIG. 5.

The free energy profile as a function of pore radius for a DOPC lipid bilayer. In the red curve, the lipid headgroups are restrained to lie near one of the two membrane surfaces, whereas in the black and blue curves, no such restraints are applied. The hydrophilic pore curve is extrapolated back to zero radius.

FIG. 5.

The free energy profile as a function of pore radius for a DOPC lipid bilayer. In the red curve, the lipid headgroups are restrained to lie near one of the two membrane surfaces, whereas in the black and blue curves, no such restraints are applied. The hydrophilic pore curve is extrapolated back to zero radius.

Close modal
FIG. 6.

Snapshots of the hydrophilic pore of DOPC at a pore radius of 5.0 Å. (a) 12.119 ns, (b) 12.125 ns, (c) 14 ns, and (d) 30 ns.

FIG. 6.

Snapshots of the hydrophilic pore of DOPC at a pore radius of 5.0 Å. (a) 12.119 ns, (b) 12.125 ns, (c) 14 ns, and (d) 30 ns.

Close modal

With headgroup restraints to maintain a hydrophobic pore, there is a slight increase in the slope as the radius is increased. The hydrophobic pore is never fully hydrated in these simulations; only a partial entry of water molecules near the entrances is observed. In Sec. III C below, we determine that wet hydrophobic pores are more stable than dry ones for R > 6.5 Å. Thus, the pores here remain dry for kinetic reasons.

During the growth of hydrophobic pores, we observed spontaneous wetting in DMPC at R = 7 Å, whereas in DOPC, no wetting was observed up to 10 Å. Because the wetting process is independent from our control variable R, the free energy of wetting needs to be estimated by a different method. We have used the density bias method,41 which restrains the water density in a user-specified cylinder using radial and axial switching functions to provide a continuous control variable. The results of this calculation are very sensitive to the type of restraints used for the lipids (see the  Appendix). To obtain results compatible and consistent with the pore growth free energy profiles, we used the same type of restraints used in those calculations (one-sided harmonic exclusion from a cylinder of a given radius).

Figures 7 and 8 show free energy profiles of wetting of hydrophobic pores in DMPC and DOPC, respectively. The x-axis in these plots is the atom number density (VDEN), which for bulk water is 100 nm−3 (here, this variable can be larger than 100 because of the large RW switching distance; see Sec. II) For DMPC at 7 Å, where spontaneous wetting is observed within 15 ns, the ΔG is −12.4 kcal/mol ± 0.3 kcal/mol, and at 5 Å, ΔG is +1.2 kcal/mol ± 0.3 kcal/mol. At R = 3 Å, there is no local minimum in the wet state. The pore fills at about VDEN = 110, and then, additional water expands the pore. As the free energy of wetting the 3-Å pore, we select the value at 110, which is +12.3 kcal/mol ± 0.3 kcal/mol. For DOPC, the values are −5.2 kcal/mol ± 0.1 kcal/mol at 7 Å and +16.3 kcal/mol ± 0.2 kcal/mol at 5 Å. By linear interpolation, we find that ΔG = 0 at about R = 6.55 Å. Even though wetting is thermodynamically favorable at 7 Å, the barrier is too large to be traversed in a few nanoseconds. Significant kinetic barriers to capillary condensation and evaporation have been observed in previous studies of hydrophobic wetting/dewetting.64–66 

FIG. 7.

Free energy of wetting DMPC hydrophobic pores as a function of atom number density. R = 7 Å (black curve), R = 5 Å (red curve), and R = 3 Å (blue curve).

FIG. 7.

Free energy of wetting DMPC hydrophobic pores as a function of atom number density. R = 7 Å (black curve), R = 5 Å (red curve), and R = 3 Å (blue curve).

Close modal
FIG. 8.

Free energy of wetting DOPC hydrophobic pores as a function of atom number density. R = 7 Å (black curve) and R = 5 Å (red curve).

FIG. 8.

Free energy of wetting DOPC hydrophobic pores as a function of atom number density. R = 7 Å (black curve) and R = 5 Å (red curve).

Close modal

From Eq. (1) (without the surface tension term), the slope of the free energy profile vs R is equal to 2πγ. Experiments probe hydrophilic (toroidal) pores, as they are clearly more stable. Our free energy profiles for hydrophilic pores are very close to linear. Only in DMPC, there is a slight upward curvature, implying a slightly lower line tension for small pores. Ignoring that, we have fitted the free energy profile of the hydrophilic pore in the range R = 3.0 Å–10.5 Å using a linear regression. To assess the convergence of this result, we performed this calculation over 2-ns segments of the 40-ns simulations and also cumulatively after the first 10 ns (Fig. S3 in the supplementary material). The initial values are high, but the calculation appears to have converged after 15 ns of simulation. The average value of the line tension coefficient for 10 ns–40 ns is 37.6 pN. There are no experimental data for DMPC, but a value of 31 pN was estimated from elasticity theory based on the bending modulus and the thickness.67,68 For the (unrealistic) hydrophobic pores, the calculated value of the line tension is 54.5 pN in the range of R of 7.0 Å–10.5 Å (where the pores are fully hydrated) and 105.3 pN in the range of R of 0.0 Å–7.0 Å (where the pores are dry).

The line tension of the DOPC bilayer was obtained by fitting the free energy profile of the hydrophilic pore in the range of pore radii of R = 5.0 Å–10.5 Å. The plot of line tension vs simulation time is given in Fig. S4 in the supplementary material. Again, the calculation appears converged after about 15 ns. The average value of the line tension is 53.7 pN for 10 ns–48 ns. Experimental values are in the range of 7 pN–28 pN,69 whereas previous theoretical calculations gave a value of 45 pN.68 The computed line tension of the hydrophobic pore of DOPC is 127 pN in the range of pore radii of 0.0 Å–8.0 Å.

The slopes observed in the free energy profiles of the hydrophobic pores can be rationalized in terms of interfacial tensions. A schematic diagram of a hydrophobic pore before and after wetting is shown in Fig. 9. When a dry hydrophobic pore grows, four interfaces increase in size: the acyl chain–vacuum interface (thickness t, interfacial tension σa), the headgroup–water interface (thickness h, interfacial tension σh), and the two water–vacuum interfaces (surface tension σw) at the entrances of the pore. In addition, there is a contribution from three-phase line tension τ at the locus of points where vapor (vacuum), water, and lipids meet.70–72 Thus, the free energy of a dry hydrophobic pore of radius R is

Gdry=2π >R t> σa+2π >R hσh+2πR2σw+22π >R> τ.
(4)

The presence of the quadratic term in Eq. (4) could explain why the free energy profile for the DOPC dry pore increases gradually in slope. When a wet hydrophobic pore grows, only the alkane–water and headgroup–water interfaces increase,

Gwet=2π >R t> σaw+2π >R hσh=2π >R> γ,
(5)

where σaw is the alkane–water interfacial tension. By analyzing the average z positions of choline N and carbonyl C atoms, we estimate t to be 22 Å for DMPC and 27 Å for DOPC and h to be about 12 Å for both.

FIG. 9.

Schematic diagram of a hydrophobic pore before wetting (left) and after wetting (right). The red ellipses are the polar headgroups, and the yellow lines are the acyl chains.

FIG. 9.

Schematic diagram of a hydrophobic pore before wetting (left) and after wetting (right). The red ellipses are the polar headgroups, and the yellow lines are the acyl chains.

Close modal

Using experimental values for the interfacial and surface tensions to rationalize our data presents some difficulties: (a) the force field predictions for these properties often deviate significantly from experimental values;73 (b) the microscopic interfaces in our simulations deviate from macroscopic interfaces; (c) to our knowledge, no experimental values are available for the lipid headgroup–water interfacial tension and for the specific three-phase line tension in the present systems.

We will try to rationalize the observed slopes and extract a value for τ using the calculated values of surface and interfacial tensions for the force field we are using.73 From Eq. (5), for DMPC, using σaw = 44.9 mN/m,73 we get σh = −51.1 mN/m. A negative interfacial tension is reasonable for headgroups, because phosphocholine and glycerophosphocholine by themselves are water-soluble. Because the quadratic term in Eq. (4) is the water surface term, σw could be determined from the curvature of the dry pore profiles,

d2G/dR2=4πσw.
(6)

Fitting the DOPC dry pore curve (which occurs over a larger range of R) to a second degree polynomial, we obtain σw = 55.5 mN/m. For comparison, the value calculated for a free water surface for the TIP3P model using a similar nonpolar truncation scheme as here is 55.3 mN/m.74 

Assuming the same surface tension for the acyl chains of DMPC and DOPC, the slopes of the dry pore profiles at R = 6.5 Å, divided by 2π, are

DMPC> 22  σa+12  σh+2*6.5  σw+2τ1440pN/10,
(7)
DOPC> 27  σa+12  σh+2*6.5  σw+2τ1510pN/10.
(8)

Subtracting these two, we find σa = 14 mN/m, which is close to the literature value.73 Using this value for σa and values for σh and σw derived above in Eq. (7), we get τ = 51.2 pN. This value is within the commonly reported range of the three-phase line tension in various systems.72,75 While some studies find a positive three-phase line tension,65 others find a negative value.66,76,77 The origin of the discrepancy is not clear.

The difference between the dry and wet pores is

GdryGwet=t> σaσaw >R> +2τ >R> +σwR2.
(9)

The first term on the right is negative, and the others are positive. The values derived above predict that the energy of the dry pore is always higher than that of the wet pore. However, the wetting free energy calculations showed that pores smaller than 5 Å in DMPC and 6.5 Å in DOPC have lower free energy when dry. One possible explanation for this is the breakdown of the macroscopic approach or, equivalently, a significant increase in the acyl chain–water interfacial tension at very small pores. This makes physical sense, as water in very narrow spaces loses a lot of hydrogen bonding interactions.78 Further work is needed to prove this point.

In this work, we have used free energy perturbation of geometric restraints to obtain the free energy of pore growth and pore transformation from hydrophobic to hydrophilic in DMPC and DOPC bilayers. Combining these with wetting free energies, we obtained complete free energy profiles of hydrophobic and hydrophilic pores as functions of radius. We essentially use three progress variables: pore exclusion radius, headgroup exclusion zone, and pore water density. Pore radius is sufficient at small radii because the pores are spontaneously dry and the headgroups remain on the surface. At larger radii, the other two variables come into play and are treated by separate calculations.

The pore formation process followed in the present work is similar to that imagined in previous theoretical models.5,17 That is, a hydrophobic pore forms first, and then, when it expands sufficiently, the headgroups rearrange to make a hydrophilic pore. Glaser et al.17 predicted this to happen at R ∼ 3 Å–5 Å. Indeed, in our simulations, it happens at 3 Å for DMPC and 5 Å for DOPC. One departure from their assumptions is that we find that narrow hydrophobic pores are unlikely to be wet. Transformation of a hydrophobic pore to a hydrophilic one is coupled with the wetting of the pore. Fitting their model to electroporation data, Glaser et al. estimated a barrier of 45 kBT for the formation of a hydrophilic pore and assumed it involved the rearrangement of the headgroups. We find that this rearrangement is essentially downhill, and the barrier involves the actual opening of the hydrophobic pore. In the specific path we are following, this barrier is ∼21 kcal/mol for DMPC and ∼48 kcal/mol for DOPC.

However, we should note that the paths in configuration space traced by our free energy profiles are not necessarily realistic. Pore formation likely occurs via more complex paths, such as those depicted in previous computational works that employ collective reaction coordinates.40,41,44 Still, our profiles are useful because they refer to well-defined, idealized pore structures. For example, large-radius hydrophobic pores in pure bilayers are clearly unrealistic, but their free energy could be useful in building energy models for peptide-lined pores. In addition, wetting of hydrophobic channels is of fundamental interest.65,76,79 The process of wetting we observe here is similar to that observed in artificial hydrophobic channels.66 

Glaser et al.17 argued that there should be a local free energy minimum in the free energy profile right after the first formation of a hydrophilic pore. In their model, this metastable “prepore” state is a result of steric repulsion and hydration effects as a hydrophilic pore shrinks. A metastable state was found in a recent elasticity-based model at the point where the radius of the hydrophilic pore is equal to the hydrophobic thickness of the lipid bilayer.24 This metastable state is usually not observed in simulations,40,41 except in a recent work that used an elaborate reaction coordinate44,80 where a shallow local minimum of at most 2 kBT was found under certain conditions. In our work, the depth of this local minimum is likely exaggerated in DMPC due to the path taken, i.e., a path through a high-energy dry hydrophobic pore. Despite that, no local minimum is found in DOPC, which seems consistent with the conclusion of Ting et al. that unsaturated chains destabilize the metastable prepore minimum.80 However, our finding applies to R = 5 Å, while the nascent DOPC pore is extrapolated to be much lower in free energy. So, it is still possible that a barrier to closing exists. The “silent,” nonconductive “prepores” that some electroporation experiments suggested14 are likely nascent, zero radius hydrophilic pores. Their persistence on the millisecond time scale could perhaps be explained by further stabilization of these local minima by voltage. More work is needed to prove this. The very long (multisecond) resealing times seen in cell and planar bilayer experiments even in the absence of voltage need to be explained by a different mechanism. One explanation has been offered for pipette experiments,81 but it does not apply to other experimental setups. Our hydrophilic pore curves are reminiscent of Karal and Yamazaki’s modification of classical nucleation theory.12 However, their estimate of the energy of the nascent hydrophilic pore (4.7 kBT) is much smaller than ours and other theoretical estimates.

A number of approaches have been used to compute the free energy of a pore.39–41,43,44 Comparison of the free energy values obtained by each of them is not simple because the pore states they refer to are in principle different. However, most of them compute the free energy of a nascent pore, i.e., a pore with the smallest possible opening. For such a pore in DMPC, Hub and Awasthi44 calculated a free energy of 8.4 kcal/mol at 300 K. 11 kcal/mol–19 kcal/mol, depending on the force field, was obtained from phosphate placement at the bilayer center.82 More recently, a slightly lower value of 6 kcal/mol was obtained by Ting et al.80 We obtain a free energy of 6.8 kcal/mol for the DMPC nascent pore and 14 kcal/mol at R = 3 Å, values comparable with those obtained by the previous studies. For a nascent DOPC pore, Ting et al. calculated a free energy of 22.8 kcal/mol, and we get (by extrapolation) 25 kcal/mol. Similar values are obtained for the barrier to phospholipid flip flop using the Berger united atom force field.31 Coarse-grained simulations of DOPC estimated a much higher pore nucleation free energy of 42.3 kcal/mol.83 

At very small radii, the profile appears quadratic, as anticipated by continuum theory17 and observed in other simulation studies.40 However, the quadratic region here extends only up to 1 Å for both DMPC and DOPC, whereas in previous studies, it extended up to almost 3 Å40 or even longer.17 The quadratic dependence on R has been attributed to a linear increase in the alkane surface tension with radius17 or to density fluctuations.40 A simpler explanation can be offered based on the concepts of scaled particle theory.84 The free energy of inserting a small cavity in a homogeneous solvent is ΔGc = −kBT ln(1 − 4/3πR3ρ), where R is the distance of the closest approach and ρ is the solvent density. For inserting a cylinder in a membrane of thickness h, the formula becomes ΔGc = −kBT ln(1 − hπR2ρ) ∼ R2. This formula is only valid for R less than the radius of a spherical solvent molecule. Here, the “solvent” is polyatomic lipid molecules and that may complicate the picture, but if we consider the solvent consisting of 2.5-Å diameter groups, the range of validity of the quadratic behavior observed is reasonable. The difference in the shape of the free energy profiles between the present work and that of Wohlert et al. is likely due to the different pore configurations sampled in the two works (no hydrophobic pores were observed using the collective coordinate of that work).

The hydrophilic pore free energy profiles are close to linear with respect to pore radius, suggesting a constant line tension. The calculated values of the line tension of hydrophilic pores in DMPC and DOPC gradually decrease as simulation time increases but appear to converge by about 15 ns–20 ns per window. Our best estimates for the line tension are 37.6 pN for DMPC (Fig. S3) and 53.7 pN for DOPC (Fig. S4). No experimental value has been determined for DMPC, but a value of 31 pN has been estimated based on the bending modulus and the thickness.68 The experimentally reported values of the line tension of DOPC are in the range of 7 pN–28 pN.69 Kindt and co-workers68 calculated via the pressure tensor a line tension of 45 pN for DOPC using the same force field (Charmm 36) and 12 pN–35 pN for DMPC using the Gromos force field.85 It would be interesting to calculate the full free energy profile of DOPC hydrophilic pores down to zero radius to see if the slope, and, thus, the line tension, will be lower at smaller radii, as predicted by continuum elasticity theory.24 It is noted that some studies found a strong dependence of pore formation free energies on the force field82 but others did not.43,68 This apparent discrepancy seems to be due to the hysteresis and sampling problems present when the phosphate position reaction coordinate is used.44,86

Compared to previous work, our approach has advantages and disadvantages. One main advantage is that the pore radius is the most natural variable to characterize a pore and appears in classical nucleation theory so that direct comparisons are possible. This choice guarantees the presence of a true pore, whereas other collective coordinates based on water density or lipid lateral position cannot distinguish between water density that is contiguous in space or dispersed in small packets. Some can also lead to unphysical negative values for the radius.40 Furthermore, to obtain the free energy as a function of pore radius using the other collective variables, one has to repeat the calculation from scratch for every radius, whereas here, one can simply start from an existing pore and increase its radius. Because the radius is the pre-eminent characteristic of an aqueous pore, lack of control over this parameter is a serious deficiency. One disadvantage of the present method is that the calculation can only be done in the direction of increasing pore radius. In addition, it is likely that the radius, as defined here, is not the true reaction coordinate for pore nucleation in the highly multidimensional space of lipid and water degrees of freedom. For that, other approaches with emphasis on pore nucleation are more useful; the strength of the present approach is in pore growth.

A model based on fundamental interfacial tensions can reproduce the hydrophobic pore results only with a positive three-phase line tension. Other simulation studies of nanopore wetting are also consistent with the positive value of this parameter,65 while others predict a negative value.66 Whether this discrepancy is caused by differences in force field parameters or system characteristics needs to be further investigated. Most experimental studies suggest a positive value for the three-phase line tension,87 while most theoretical studies indicate a magnitude from 1 pN to 100 pN.72 

The present approach for computing the free energy as a function of pore radius is conceptually simple and easily implemented. Calculated pore line tensions could provide additional targets for force field parameterization, although the large variability in experimental values is somewhat disheartening. Free energy perturbation of restraints could potentially be useful in obtaining the dependence of pore free energies on other geometrical characteristics such as curvature or lipid composition of the pore edge.

See the supplementary material for Tables S1 and S2, which show free energy changes for different values of the force constant and different choices of the pore radius increment, respectively. The average distances between the phosphate groups of the upper leaflet and the lower leaflet for DMPC and DOPC are given in Fig. S1. A snapshot of the hydrophilic pore of DOPC is given in Fig. S2. The convergence of line tension calculations is shown in Figs. S3 and S4.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the National Institutes of Health (Grant No. GM117146). Computational resources were provided by XSEDE, which is supported by the National Science Foundation (Grant No. TG-MCB170009), and the high performance computing facilities of The City University of New York, supported by the NSF. Infrastructure support was provided in part by Research Centers in Minority Institutions (Grant No. 8G12MD007603) from the National Institutes of Health. The authors thank Professor Alan Grossfield for clarifications on his WHAM code.

We attempted to calculate the free energy of wetting of hydrophobic pores using different types of restraints and found that the type of restraint influences the results dramatically. In the main text, we report the results obtained using the same type of restraint as that in the pore growth free energy simulations (exclusion from the cylinder using a one-sided harmonic restraint with a force constant of k = 10 kcal mol−1 Å−2). This is appropriate for obtaining results compatible and consistent with the pore growth free energy profiles (Figs. 2 and 5). When we used harmonic position restraints on all lipid C atoms, we obtained results that depend dramatically on the force constant used. For example, Fig. 10 shows free energy profiles for the DMPC pore at R = 7 Å using different values for the force constant. With k = 5 kcal mol−1 Å−2, we obtain ΔG = −47 kcal/mol, but with k = 0.4 kcal mol−1 Å−2, we obtain ΔG = +17.5 kcal/mol!.

FIG. 10.

Free energy profiles for wetting a 7-Å DMPC pore for different values of the force constant used to restrain the position of the lipids: k = 5 kcal mol−1 Å−2 (black curve) and k = 0.4 kcal mol−1 Å−2 (red curve).

FIG. 10.

Free energy profiles for wetting a 7-Å DMPC pore for different values of the force constant used to restrain the position of the lipids: k = 5 kcal mol−1 Å−2 (black curve) and k = 0.4 kcal mol−1 Å−2 (red curve).

Close modal

This dependence of the wetting behavior on restraint strength has been observed in previous studies.88–90 It was partly attributed to an effective narrowing of the space available to water.89,91 However, this is not enough to explain the difference.88 The main reason for this dependence on the force constant is that the lipid–water interaction becomes more favorable, and thus, the lipid–water interfacial tension becomes lower, as the lipids become more rigid. This is readily verified by calculating average interaction energies between pore water and lipids for simulations run using different restraint force constants.

1.
M. D.
Peraro
and
F. G.
Van Der Goot
, “
Pore-forming toxins: Ancient, but never really out of fashion
,”
Nat. Rev. Microbiol.
14
(
2
),
77
92
(
2016
).
2.
K. A.
Brogden
, “
Antimicrobial peptides: Pore formers or metabolic inhibitors in bacteria?
,”
Nat. Rev. Microbiol.
3
(
3
),
238
250
(
2005
).
3.
T.
Kotnik
,
W.
Frey
,
M.
Sack
,
S.
Haberl Meglič
,
M.
Peterka
, and
D.
Miklavčič
, “
Electroporation-based applications in biotechnology
,”
Trends Biotechnol.
33
(
8
),
480
488
(
2015
).
4.
S.
Majd
,
E. C.
Yusko
,
Y. N.
Billeh
,
M. X.
Macrae
,
J.
Yang
, and
M.
Mayer
, “
Applications of biological pores in nanomedicine, sensing, and nanoelectronics
,”
Curr. Opin. Biotechnol.
21
(
4
),
439
476
(
2010
).
5.
I. G.
Abidor
,
V. B.
Arakelyan
,
L. V.
Chernomordik
,
Y. A.
Chizmadzhev
,
V. F.
Pastushenko
, and
M. P.
Tarasevich
, “
Electric breakdown of bilayer lipid membranes. I. The main experimental facts and their qualitative discussion
,”
J. Electroanal. Chem.
104
(
C
),
37
52
(
1979
).
6.
T. Y.
Tsong
, “
Electroporation of cell membranes
,”
Biophys. J.
60
(
2
),
297
306
(
1991
).
7.
J. C.
Weaver
and
Y. A.
Chizmadzhev
, “
Theory of electroporation: A review
,”
Bioelectrochem. Bioenerg.
41
(
2
),
135
160
(
1996
).
8.
C.
Wilhelm
,
M.
Winterhalter
,
U.
Zimmermann
, and
R.
Benz
, “
Kinetics of pore size during irreversible electrical breakdown of lipid bilayer membranes
,”
Biophys. J.
64
(
1
),
121
128
(
1993
).
9.
M.
Casciola
and
M.
Tarek
, “
A molecular insight into the electro-transfer of small molecules through electropores driven by electric fields
,”
Biochim. Biophys. Acta, Biomembr.
1858
(
10
),
2278
2289
(
2016
).
10.
C.
Taupin
,
M.
Dvolaitzky
, and
C.
Sauterey
, “
Osmotic-pressure induced pores in phospholipid vesicles
,”
Biochemistry
14
(
21
),
4771
4775
(
1975
).
11.
E.
Evans
,
V.
Heinrich
,
F.
Ludwig
, and
W.
Rawicz
, “
Dynamic tension spectroscopy and strength of biomembranes
,”
Biophys. J.
85
(
4
),
2342
2350
(
2003
).
12.
M. A. S.
Karal
and
M.
Yamazaki
, “
Communication: Activation energy of tension-induced pore formation in lipid membranes
,”
J. Chem. Phys.
143
(
8
),
081103
(
2015
).
13.
D. V.
Zhelev
and
D.
Needham
, “
Tension-stabilized pores in giant vesicles - determination of pore-size and pore line tension
,”
Biochim. Biophys. Acta
1147
(
1
),
89
104
(
1993
).
14.
K. C.
Melikov
,
V. A.
Frolov
,
A.
Shcherbakov
,
A. V.
Samsonov
,
Y. A.
Chizmadzhev
, and
L. V.
Chernomordik
, “
Voltage-induced nonconductive pre-pores and metastable single pores in unmodified planar lipid bilayer
,”
Biophys. J.
80
(
4
),
1829
1836
(
2001
).
15.
E. A.
Kotova
,
A. V.
Kuzevanov
,
A. A.
Pashkovskaya
, and
Y. N.
Antonenko
, “
Selective permeabilization of lipid membranes by photodynamic action via formation of hydrophobic defects or pre-pores
,”
Biochim. Biophys. Acta, Biomembr.
1808
(
9
),
2252
2257
(
2011
).
16.
J. D.
Litster
, “
Stability of lipid bilayers and red blood cell membranes
,”
Phys. Lett. A
53
(
3
),
193
194
(
1975
).
17.
R. W.
Glaser
,
S. L.
Leikin
,
L. V.
Chernomordik
,
V. F.
Pastushenko
, and
A. I.
Sokirko
, “
Formation and evolution of pores
,”
Biochim. Biophys. Acta
940
(
1988
),
275
287
(
1988
).
18.
S.
May
, “
A molecular model for the line tension of lipid membranes
,”
Eur. Phys. J. E
3
,
37
44
(
2000
).
19.
V.
Talanquer
and
D. W.
Oxtoby
, “
Nucleation of pores in amphiphile bilayers
,”
J. Chem. Phys.
118
(
2
),
872
877
(
2003
).
20.
C. L.
Ting
,
D.
Appelö
, and
Z. G.
Wang
, “
Minimum energy path to membrane pore formation and rupture
,”
Phys. Rev. Lett.
106
(
16
),
168101
(
2011
).
21.
J.
Li
,
K. A.
Pastor
,
A. C.
Shi
,
F.
Schmid
, and
J.
Zhou
, “
Elastic properties and line tension of self-assembled bilayer membranes
,”
Phys. Rev. E
88
(
1
),
012718
(
2013
).
22.
A.
Dehghan
,
K. A.
Pastor
, and
A. C.
Shi
, “
Line tension of multicomponent bilayer membranes
,”
Phys. Rev. E
91
(
2
),
022713
(
2015
).
23.
H.
Pera
,
J. M.
Kleijn
, and
F. A. M.
Leermakers
, “
On the edge energy of lipid membranes and the thermodynamic stability of pores
,”
J. Chem. Phys.
142
(
3
),
034101
(
2015
).
24.
S. A.
Akimov
,
P. E.
Volynsky
,
T. R.
Galimzyanov
,
P. I.
Kuzmin
,
K. V.
Pavlov
, and
O. V.
Batishchev
, “
Pore formation in lipid membrane I: Continuous reversible trajectory from intact bilayer through hydrophobic defect to transversal pore
,”
Sci. Rep.
7
(
1
),
12152
(
2017
).
25.
O.
Farago
, “‘
Water-free’ computer model for fluid bilayer membranes
,”
J. Chem. Phys.
119
(
1
),
596
605
(
2003
).
26.
I. R.
Cooke
,
K.
Kremer
, and
M.
Deserno
, “
Tunable generic model for fluid bilayer membranes
,”
Phys. Rev. E
72
(
1
),
011506
(
2005
).
27.
Z. J.
Wang
and
D.
Frenkel
, “
Pore nucleation in mechanically stretched bilayer membranes
,”
J. Chem. Phys.
123
(
15
),
154701
(
2005
).
28.
T. V.
Tolpekina
,
W. K.
den Otter
, and
W. J.
Briels
, “
Nucleation free energy of pore formation in an amphiphilic bilayer studied by molecular dynamics simulations
,”
J. Chem. Phys.
121
(
23
),
12060
12066
(
2004
).
29.
C.
Loison
,
M.
Mareschal
, and
F.
Schmid
, “
Pores in bilayer membranes of amphiphilic molecules: Coarse-grained molecular dynamics simulations compared with simple mesoscopic models
,”
J. Chem. Phys.
121
(
4
),
1890
1900
(
2004
).
30.
W. K.
den Otter
, “
Free energies of stable and metastable pores in lipid membranes under tension
,”
J. Chem. Phys.
131
(
20
),
205101
(
2009
).
31.
W. F. D.
Bennett
and
D. P.
Tieleman
, “
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
(
2011
).
32.
D. P.
Tieleman
,
H.
Leontiadou
,
A. E.
Mark
, and
S.-J.
Marrink
, “
Simulation of pore formation in lipid bilayers by mechanical stress and electric fields
,”
J. Am. Chem. Soc.
125
(
21
),
6382
6383
(
2003
).
33.
H.
Leontiadou
,
A. E.
Mark
, and
S. J.
Marrink
, “
Molecular dynamics simulations of hydrophilic pores in lipid bilayers
,”
Biophys. J.
86
(
4
),
2156
2164
(
2004
).
34.
D. P.
Tieleman
, “
The molecular basis of electroporation
,”
BMC Biochem.
5
,
10
(
2004
).
35.
M.
Tarek
, “
Membrane electroporation: A molecular dynamics simulation
,”
Biophys. J.
88
(
6
),
4045
4053
(
2005
).
36.
R. A.
Böckmann
,
B. L.
De Groot
,
S.
Kakorin
,
E.
Neumann
, and
H.
Grubmüller
, “
Kinetics, statistics, and energetics of lipid membrane electroporation studied by molecular dynamics simulations
,”
Biophys. J.
95
(
4
),
1837
1850
(
2008
).
37.
T. J.
Piggot
,
D. A.
Holdbrook
, and
S.
Khalid
, “
Electroporation of the E. Coli and S. Aureus membranes: Molecular dynamics simulations of complex bacterial membranes
,”
J. Phys. Chem. B
115
(
45
),
13381
13388
(
2011
).
38.
D. P.
Tieleman
and
S.-J.
Marrink
, “
Lipids out of equilibrium: Energetics of desorption and pore mediated flip-flop
,”
J. Am. Chem. Soc.
128
(
38
),
12462
12467
(
2006
).
39.
W. F. D.
Bennett
,
N.
Sapay
, and
D. P.
Tieleman
, “
Atomistic simulations of pore formation and closure in lipid bilayers
,”
Biophys. J.
106
(
1
),
210
219
(
2014
).
40.
J.
Wohlert
,
W. K.
den Otter
,
O.
Edholm
, and
W. J.
Briels
, “
Free energy of a trans-membrane pore calculated from atomistic molecular dynamics simulations
,”
J. Chem. Phys.
124
(
15
),
154905
(
2006
).
41.
V.
Mirjalili
and
M.
Feig
, “
Density-biased sampling: A robust computational method for studying pore formation in membranes
,”
J. Chem. Theory Comput.
11
(
1
),
343
350
(
2015
).
42.
G. M.
Torrie
and
J. P.
Valleau
, “
Nonphysical sampling distribution in Monte Carlo free energy estimation: Umbrella sampling
,”
J. Comput. Phys.
23
,
187
199
(
1977
).
43.
N.
Awasthi
and
J. S.
Hub
, “
Simulations of pore formation in lipid membranes: Reaction coordinates, convergence, hysteresis, and finite-size effects
,”
J. Chem. Theory Comput.
12
(
7
),
3261
3269
(
2016
).
44.
J. S.
Hub
and
N.
Awasthi
, “
Probing a continuous polar defect: A reaction coordinate for pore formation in lipid membranes
,”
J. Chem. Theory Comput.
13
(
5
),
2352
2366
(
2017
).
45.
R. W.
Zwanzig
, “
High-temperature equation of state by a perturbation method I. Nonpolar gases
,”
J. Chem. Phys.
22
,
1420
1426
(
1954
).
46.
B. R.
Brooks
,
C. L.
Brooks
 III
,
A. D.
Mackerell
, Jr.
,
L.
Nilsson
,
R. J.
Petrella
,
B.
Roux
,
Y.
Won
,
G.
Archontis
,
C.
Bartels
,
S.
Boresch
, et al, “
CHARMM: The biomolecular simulation program
,”
J. Comput. Chem.
30
(
10
),
1545
1614
(
2009
).
47.
N. G.
Parsonage
, “
Determination of the chemical potential by the particle insertion method and by its inverse
,”
J. Chem. Soc., Faraday Trans.
91
(
17
),
2971
2973
(
1995
).
48.
D. A.
Kofke
and
P. T.
Cummings
, “
Precision and accuracy of staged free-energy perturbation methods for computing the chemical potential by molecular simulation
,”
Fluid Phase Equilib.
150-151
(
151
),
41
49
(
1998
).
49.
G. C.
Boulougouris
,
I. G.
Economou
, and
D. N.
Theodorou
, “
On the calculation of the chemical potential using the particle deletion scheme
,”
Mol. Phys.
96
(
6
),
905
913
(
1999
).
50.
I.
Cabeza De Vaca
,
R.
Zarzuela
,
J.
Tirado-Rives
, and
W. L.
Jorgensen
, “
Robust free energy perturbation protocols for creating molecules in solution
,”
J. Chem. Theory Comput.
15
(
7
),
3941
3948
(
2019
).
51.
S.
Kumar
,
J. M.
Rosenberg
,
D.
Bouzida
,
R. H.
Swendsen
, and
P. A.
Kollman
, “
The weighted histogram analysis method for free-energy calculations on biomolecules .1. The method
,”
J. Comput. Chem.
13
(
8
),
1011
1021
(
1992
).
52.
A.
Grossfield
, WHAM: An implementation of the weighted histogram analysis method.
53.
J. S.
Hub
,
B. L.
De Groot
, and
D.
Van Der Spoel
, “
G-Whams-a free weighted histogram analysis implementation including robust error and autocorrelation estimates
,”
J. Chem. Theory Comput.
6
(
12
),
3713
3720
(
2010
).
54.
L.
Sun
,
J. K.
Noel
,
J. I.
Sulkowska
,
H.
Levine
, and
J. N.
Onuchic
, “
Connecting thermal and mechanical protein (Un)folding landscapes
,”
Biophys. J.
107
(
12
),
2950
2961
(
2014
).
55.
E. L.
Wu
,
X.
Cheng
,
S.
Jo
,
H.
Rui
,
K. C.
Song
,
E. M.
Dávila-Contreras
,
Y.
Qi
,
J.
Lee
,
V.
Monje-Galvan
,
R. M.
Venable
, et al, “
CHARMM-GUI membrane builder toward realistic biological membrane simulations
,”
J. Comput. Chem.
35
(
27
),
1997
2004
(
2014
).
56.
S.
Jo
,
T.
Kim
,
V. G.
Iyer
, and
W.
Im
, “
CHARMM-GUI: A web-based graphical user interface for CHARMM
,”
J. Comput. Chem.
29
(
11
),
1859
1865
(
2008
).
57.
J. B.
Klauda
,
R. M.
Venable
,
J. A.
Freites
,
J. W.
O’Connor
,
D. J.
Tobias
,
C.
Mondragon-Ramirez
,
I.
Vorobyov
,
A. D.
MacKerell
, Jr.
, and
R. W.
Pastor
, “
Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types
,”
J. Phys. Chem. B
114
(
23
),
7830
7843
(
2010
).
58.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
, “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
,
926
935
(
1983
).
59.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon Press
,
Oxford
,
1989
).
60.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An NlogN method for Ewald sums in large systems
,”
J. Chem. Phys.
98
(
12
),
10089
10092
(
1993
).
61.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
62.
S. E.
Feller
,
Y.
Zhang
,
R. W.
Pastor
, and
B. R.
Brooks
, “
Constant-pressure molecular-dynamics simulation - the Langevin piston method
,”
J. Chem. Phys.
103
(
11
),
4613
4621
(
1995
).
63.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J. C.
Berendsen
, “
Numerical-integration of Cartesian equations of motion of a system with constraints - molecular-dynamics of N-alkanes
,”
J. Comput. Phys.
23
(
3
),
327
341
(
1977
).
64.
O.
Beckstein
and
M. S. P.
Sansom
, “
Liquid-vapor oscillations of water in hydrophobic nanopores
,”
Proc. Natl. Acad. Sci. U. S. A.
100
(
12
),
7063
7068
(
2003
).
65.
S.
Sharma
and
P. G.
Debenedetti
, “
Free energy barriers to evaporation of water in hydrophobic confinement
,”
J. Phys. Chem. B
116
(
44
),
13282
13289
(
2012
).
66.
A.
Tinti
,
A.
Giacomello
,
Y.
Grosu
, and
C. M.
Casciola
, “
Intrusion and extrusion of water in hydrophobic nanopores
,”
Proc. Natl. Acad. Sci. U. S. A.
114
(
48
),
E10266
E10273
(
2017
).
67.
N.
Kučerka
,
Y.
Liu
,
N.
Chu
,
H. I.
Petrache
,
S.
Tristram-Nagle
, and
J. F.
Nagle
, “
Structure of fully hydrated fluid phase DMPC and DLPC lipid bilayers using X-ray scattering from oriented multilamellar arrays and from unilamellar vesicles
,”
Biophys. J.
88
(
4
),
2626
2637
(
2005
).
68.
A.
West
,
K.
Ma
,
J. L.
Chung
, and
J. T.
Kindt
, “
Simulation studies of structure and edge tension of lipid bilayer edges: Effects of tail structure and force-field
,”
J. Phys. Chem. A
117
(
32
),
7114
7123
(
2013
).
69.
T.
Portet
and
R.
Dimova
, “
A new method for measuring edge tensions and stability of lipid bilayers: Effect of membrane composition
,”
Biophys. J.
99
(
10
),
3264
3273
(
2010
).
70.
B. V.
Toshev
,
D.
Platikanov
, and
A.
Scheludko
, “
Line tension in three-phase equilibrium systems
,”
Langmuir
4
(
3
),
489
499
(
1988
).
71.
B.
Widom
, “
Line tension and the shape of a sessile drop
,”
J. Phys. Chem.
99
(
9
),
2803
2806
(
1995
).
72.
J.
Drelich
, “
The significance and magnitude of the line tension in three-phase (solid-liquid-fluid) systems
,”
Colloids Surf., A
116
(
1–2
),
43
54
(
1996
).
73.
S. E.
Feller
,
Computational Modelling of Membrane Bilayers
(
Academic Press
,
Amsterdam
,
2008
).
74.
J. B.
Klauda
,
X.
Wu
,
R. W.
Pastor
, and
B. R.
Brooks
, “
Long-range Lennard-Jones and electrostatic interactions in interfaces: Application of the isotropic periodic sum method
,”
J. Phys. Chem. B
111
(
17
),
4393
4400
(
2007
).
75.
J.
Drelich
and
J.
Miller
, “
The line/pseudo-line tension in three-phase systems
,”
Part. Sci. Technol.
10
,
1
20
(
1992
).
76.
L.
Guillemot
,
T.
Biben
,
A.
Galarneau
,
G.
Vigier
, and
E.
Charlaix
, “
Activated drying in hydrophobic nanopores and the line tension of water
,”
Proc. Natl. Acad. Sci. U. S. A.
109
(
48
),
19557
19562
(
2012
).
77.
M.
Ghasemi
,
S. M.
Ramsheh
, and
S.
Sharma
, “
Quantitative assessment of thermodynamic theory in elucidating the behavior of water under hydrophobic confinement
,”
J. Phys. Chem. B
122
(
50
),
12087
12096
(
2018
).
78.
J. C.
Rasaiah
,
S.
Garde
, and
G.
Hummer
, “
Water in nonpolar confinement: From nanotubes to proteins and beyond
,”
Annu. Rev. Phys. Chem.
59
(
1
),
713
740
(
2008
).
79.
R. C.
Remsing
,
E.
Xi
,
S.
Vembanur
,
S.
Sharma
,
P. G.
Debenedetti
,
S.
Garde
,
A. J.
Patel
, and
K. A.
Dill
, “
Pathways to dewetting in hydrophobic confinement
,”
Proc. Natl. Acad. Sci. U. S. A.
112
(
27
),
8181
8186
(
2015
).
80.
C. L.
Ting
,
N.
Awasthi
,
M.
Müller
, and
J. S.
Hub
, “
Metastable prepores in tension-free lipid bilayers
,”
Phys. Rev. Lett.
120
(
12
),
128103
128109
(
2018
).
81.
J. D.
Moroz
and
P.
Nelson
, “
Dynamically stabilized pores in bilayer membranes
,”
Biophys. J.
72
(
5
),
2211
2216
(
1997
).
82.
W. F. D.
Bennett
,
C. K.
Hong
,
Y.
Wang
, and
D. P.
Tieleman
, “
Antimicrobial peptide simulations and the influence of force field on the free energy for pore formation in lipid bilayers
,”
J. Chem. Theory Comput.
12
(
9
),
4524
4533
(
2016
).
83.
Y.
Hu
,
S. K.
Sinha
, and
S.
Patel
, “
Investigating hydrophilic pores in model lipid bilayers using molecular simulations: Correlating bilayer properties with pore-formation thermodynamics
,”
Langmuir
31
(
24
),
6615
6631
(
2015
).
84.
H.
Reiss
,
H. L.
Frisch
, and
J. L.
Lebowitz
, “
Statistical mechanics of rigid spheres
,”
J. Chem. Phys.
31
(
2
),
369
380
(
1959
).
85.
F. Y.
Jiang
,
Y.
Bouret
, and
J. T.
Kindt
, “
Molecular dynamics simulations of the lipid bilayer edge
,”
Biophys. J.
87
(
1
),
182
192
(
2004
).
86.
K.
Huang
and
A. E.
García
, “
Effects of truncating van Der Waals interactions in lipid bilayer simulations
,”
J. Chem. Phys.
141
(
10
),
105101
(
2014
).
87.
A.
Amirfazli
and
A. W.
Neumann
, “
Status of the three-phase line tension: A review
,”
Adv. Colloid Interface Sci.
110
(
3
),
121
141
(
2004
).
88.
O.
Beckstein
and
M. S. P.
Sansom
, “
The influence of geometry, surface character, and flexibility on the permeation of ions and water through biological pores
,”
Phys. Biol.
1
(
1
),
42
52
(
2004
).
89.
S.
Andreev
,
D.
Reichman
, and
G.
Hummer
, “
Effect of flexibility on hydrophobic behavior of nanotube water channels
,”
J. Chem. Phys.
123
(
19
),
194502
(
2005
).
90.
Y. E.
Altabet
,
A.
Haji-Akbari
, and
P. G.
Debenedetti
, “
Effect of material flexibility on the thermodynamics and kinetics of hydrophobically induced evaporation of water
,”
Proc. Natl. Acad. Sci. U. S. A.
114
(
13
),
E2548
E2555
(
2017
).
91.
Y. E.
Altabet
and
P. G.
Debenedetti
, “
The role of material flexibility on the drying transition of water between hydrophobic objects: A thermodynamic analysis
,”
J. Chem. Phys.
141
,
18C531
(
2014
).

Supplementary Material