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.

## I. INTRODUCTION

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 pathogens^{1} 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 pressure^{10} 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}

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 πR^{2}, σ 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 implicitly^{25–27} or explicitly.^{28–30} In these studies, pores may appear spontaneously^{25,29} or may be induced by stretching^{26,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 tension^{32,33} or electric field^{32,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 center^{38} 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 approach^{28} 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 Awasthi^{44} 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 method^{45} 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 method^{41} 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.

## II. METHODS

### A. Free energy perturbation of restraints

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),

where r_{i} 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}

where k_{B} 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 U_{R+Δ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 (k_{po4}) 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.

### B. Free energy of wetting by the density bias method

The density bias method^{41} with umbrella sampling^{42} 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 method^{51} 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.

### C. Computational details

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 builder^{55} module of the CHARMM–GUI web server.^{56} The CHARMM36 force field^{57} was used for the lipids. 140 DMPC molecules were hydrated with 5920 TIP3P^{58} 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 integrator^{59} with a 2-fs time step was used to integrate the equations of motion. The particle mesh Ewald method^{60} 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}

## III. RESULTS

### A. Free energy of pore growth in the DMPC bilayer

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.

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 method^{41} 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.

### B. Free energy of pore growth in the DOPC bilayer

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.

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.

### C. Free energy of wetting of hydrophobic pores

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}

### D. Line tension of DMPC and DOPC

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 Å.

### E. Analysis of hydrophobic pore results based on interfacial tensions

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

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,

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.

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,

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

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

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.

## IV. DISCUSSION

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 k_{B}T 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 coordinate^{44,80} where a shallow local minimum of at most 2 k_{B}T 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 suggested^{14} 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 k_{B}T) 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 Awasthi^{44} 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 theory^{17} 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 radius^{17} 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 ΔG_{c} = −k_{B}T ln(1 − 4/3πR^{3}ρ), 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 ΔG_{c} = −k_{B}T ln(1 − hπR^{2}ρ) ∼ R^{2}. 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-workers^{68} 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 field^{82} 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.

## SUPPLEMENTARY MATERIAL

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.

## DATA AVAILABILITY

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

## ACKNOWLEDGMENTS

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.

### APPENDIX: INFLUENCE OF RESTRAINTS ON THE FREE ENERGY OF WETTING A HYDROPHOBIC PORE

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!.

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.