Molecular dynamics simulations of ethanol permeation through single and double-lipid bilayers

Transport of small molecules across phospholipid membranes is an important biological process in drug design, toxicology, and production of biofuels such as ethanol, a widely used chemical in industry. Ethanol, due to its amphipathic character, is an important representative of molecules, which can modulate the properties of the membrane.¹ Recently, bacteria and yeast have been exploited to produce various chemicals such as ethanol. Yeast cells (Saccharomyces cerevisiae) can produce ethanol for use as a biofuel.² The yeast plasma membrane contains different lipids including inositol sphingolipid, phosphatidylserine (PS), phosphatidylinositol (PI), phosphatidylethanolamine (PE), and phosphatidylcholine (PC). Sterol ratio to phospholipid in yeast is 0.5 for the plasma membrane.³ It is known that most biofuels such as ethanol can perturb the membrane of the cells causing toxicity, which is a serious limiting factor in biofuel production.^4,5 Membrane permeation is of great importance in biofuel production where a higher yield requires tolerance to high concentrations of the biofuel produced. Information on membrane permeation of these biofuels enables engineering microbial cell membranes with higher resistance and thereby an increased yield of biofuels.^6,6 Therefore, having tools that help us to quantitatively predict the permeability of different lipid membranes to ethanol will prove useful in engineering lipid membranes with improved tolerance to the production of biofuels.

Passive permeation of small molecules through the membrane is also an essential process for the activity of produced small molecule drugs. High throughput screening of compounds based on their membrane permeation potential is an active field of pharmaceutical research, and computational approaches to calculate permeabilities are gaining increased attention due to the progress in supercomputers and computational speed.^7,8 The ultimate goal is to provide a tool for a reliable and efficient prediction of membrane permeability of the desired drug compounds. One example of a computational method to predict permeability is the quantitative structure–activity relationships (QSARs), which are based on mathematical relationships between the chemical properties of compounds and their biological activities.^9,10 However, these models lack precision due to limitations on data size and quality.¹¹ Molecular dynamics (MD) simulations and the inhomogeneous solubility-diffusion model (ISD) provide a direct measurement and physical detail about permeation of small molecules through membranes.¹² 

Membrane permeability calculation started in 1901 when Overton and Meyer observed independently that permeability of different solutes through membranes is correlated with their oil/water partition coefficient. This rule leads to a relation between membrane permeability and partition coefficient called solubility-diffusion model (or Meyer–Overton rule),

where K is the oil/water partition coefficient, D is the solute’s diffusion coefficient, and h is the membrane thickness.^13,14

A more realistic description of membrane permeation is provided by the ISD model where the bilayer is considered to be an inhomogeneous medium,¹⁵ 

Equation (2) is derived from the Smoluchowski diffusion equation at the steady state.¹² In this equation, the membrane is centered around z = 0, F_ref is the reference free energy in water, D_⊥(z) is the position-dependent diffusivity along the membrane normal, and β is 1/k_BT (k_B is Boltzmann’s constant and T is the temperature).

Accurate estimates of permeability coefficients from simulation in the ISD model depend on the reliability of free energy and diffusion profiles. From the exponential dependence of permeation on free energy, it is evident that an error as small as 2.3k_BT in F(z) can translate into changes in permeability by an order of magnitude. For permeants that show poor sampling across the bilayer, conventional MD does not give an accurate estimate of potential of mean force (PMF) and diffusion profiles. Therefore, various enhanced sampling methods have been used to compute free energy profiles. Marrink and Berensden¹² combined insertion of fictitious particles with the computation of the average force felt by molecules constrained to specific positions along the bilayer normal. Bennion et al.¹⁶ used umbrella sampling (US) to compute the PMF of a series of small-molecule drugs across the membrane. In this approach, permeants are restrained to different positions along the bilayer normal. The effect of bias is then removed from probability distributions using the Weighted Histogram Analysis Method (WHAM). Comer et al.^17,18 used an Adaptive Biasing Force (ABF) algorithm to calculate the PMF underlying the translocation of alcohols through the membrane. Sun et al.¹⁹ used transition-tempered metadynamics to obtain a reliable estimate of the PMF with two different reaction coordinates. They calculated permeability as a line integral along the minimum free energy path of the PMF. Ghaemi et al.²⁰ used bias exchange metadynamics to compute the free energy profile and built a multidimensional rate model to calculate the mean permeation time through a membrane. Considering these various approaches for PMFs, accurate estimation is difficult due to the slow convergence of degrees of freedom and abrupt transitions as a function of solute insertion depth, which is underestimated in most free energy calculations.²¹ Moreover, these methods use low concentrations of permeants, which is less relevant for designing tolerant membranes. Lee et al.²² investigated different methods for calculating PMFs including umbrella sampling (US), Adaptive Biasing Force (ABF), Replica Exchange Umbrella Sampling (REUS), and Multiple Walker ABF (MW-ABF) using three different compounds: codeine, benzoic acid, and urea. They found that all these methods converged to the same PMF.

A reliable estimation of normal diffusivity D(z) is also complicated due to the rapid change in the PMF as a function of position, which requires averaging over short timescales. On the other hand, short time dynamics include contributions from inertial and solvent memory effects, which requires averaging over long timescales to get optimal diffusivity. There are multiple methods for the estimation of normal diffusivity. The generalized Langevin equation uses restrained MD simulations to estimate the position-dependent diffusivities from autocorrelation functions. In this method, the solute is described as a harmonic oscillator, and the remainder of the system acts as the frictional bath for the solute.²³ Gaalswyk et al.²⁴ used a generalized Langevin equation for the calculation of position-dependent diffusivities where the solute is restrained to a series of positions by a harmonic potential, and diffusivity is then calculated from the position autocorrelation function (PACF) or the velocity autocorrelation function (VACF). They assessed the effect of thermostat on diffusivities calculated from the generalized Langevin equation for self-diffusion of the TIP3P-water model and a friction coefficient ranging from γ = (1–10) ps⁻¹ and compared the results to the Lowe–Anderson thermostat. The PACF and VACF methods showed a decline in diffusion coefficients as the friction constant increased, which is due to the damping of dynamics of the solute by the frictional force of the thermostat. This effect was smallest at the lowest friction coefficient of γ = 1 ps⁻¹.

A fundamentally different approach employs Bayesian inference. In this approach, the bilayer is divided into multiple bins along the transition coordinate z, and transitions between these bins are monitored at finite time intervals. A likelihood function, which is defined as the probability of observing a set of transitions between bins in MD trajectories, is used to connect the diffusive model to the observed simulation data. Diffusivity and free energy profiles are varied randomly with Monte Carlo moves, and the acceptance of the profile is determined by the likelihood function. Ghysels et al. extended this approach to determine the position-dependent diffusion coefficients parallel and normal to the bilayer.²⁵ They sampled free energy and diffusion profiles randomly with a Monte Carlo procedure to find profiles that optimally explain their observed transition counts. Krämer et al.²⁶ employed the maximum likelihood estimation (MLE) instead of Bayesian Monte Carlo sampling using a derivative free optimization method, the Covariance Matrix Adaptation-Evolution Strategy (CMA-ES).²⁷ This MLE approach is used to calculate diffusion profiles in this paper.

Another approach to calculate membrane permeability is through counting the number of crossings through the bilayer. There are two different counting methods: flux-based (FB) counting where one needs to simulate a double bilayer and calculate the concentration of the permeant in the solvent compartment vs time and a transition-based (TB) counting method where the total number of bilayer crossings are evaluated.²⁸ A double bilayer has two compartments with different concentrations of the permeant and can be used directly with Fick’s first law to calculate permeability. In contrast, single bilayer systems do not allow different concentrations across the bilayer due to the periodic boundaries, which neutralizes the concentration difference across the single bilayer. Double bilayer simulations have been used to calculate permeabilities of ions and small molecules. Venable et al.²⁸ used a double bilayer simulation to calculate the water permeability of DLPC bilayers. In addition, researchers have used double bilayer membrane systems to study ion gradients and pore formation.^29–31 For simulations with only a single water compartment, permeability can instead be calculated by counting crossings in both directions using the TB counting method for related permeation events. These types of approaches have been used in various studies of water, O₂, and CO₂ permeation through membranes and channel proteins.^28,32–37 However, the TB counting method is based on the assumption that the crossing rate is only determined by the equilibrium concentration of the permeant in water rather than the actual concentration gradient.²⁸ This assumption might break down for membrane-modulating permeants such as ethanol, especially at high concentrations where the outer leaflet is maximally perturbed by the alcohol, while the inner leaflet gets gradually saturated by the permeating molecules. Therefore, the present work studies permeation of ethanol using both double and single bilayer simulations.

Venable et al.²⁸ compared the ISD model and counting method to calculate the water permeability of different lipid membranes. However, due to the small number of crossings and insufficient sampling in the bilayer core, permeability estimates from both methods carried large uncertainties. However, for POPC, the TB counting method produced smaller error bars compared to the ISD model.

In this paper, we first use the MLE approach to calculate diffusion and free energy profiles from an unbiased MD simulation and the ISD model. Next, we present two different counting methods to calculate permeability for single and double bilayers. In Sec. IV, we first present the findings of the ISD based MLE approach for POPC, POPS, and POPE membranes and their permeabilities followed by the results of TB counting for single bilayers. Next, we present the results for double bilayer permeability calculations. At the end, we compare permeabilities calculated from different approaches and also predict the permeabilities of membranes with asymmetric concentration of ethanol from single bilayer simulations and compare to the results from double bilayers.

In this paper, three different lipids were chosen for this study: PS, PE, and PC with the palmitoyl oleoyl (PO, 18:1/16:0) chains used for all lipids. All bilayers were built using Membrane Builder in CHARMM-GUI.^38–43 Single bilayers contained 50 lipids per leaflet and 50 water molecules per lipid. For POPS membranes, K⁺ counterions were added to the system to neutralize the charges. Double bilayer systems were built by replicating single bilayers in the z direction. Scripts from CHARMM-GUI were used to minimize the system and equilibrate the system. The lipid force field in all the simulations was CHARMM36,⁴⁴ and the water model was the CHARMM-compatible TIP3P.⁴⁵ A Langevin thermostat with a friction coefficient of γ = 1 ps⁻¹ was used as the integrator following the previous work of Gaalswyk et al.⁴⁶ 30 ns post-equilibration was done for all single and double bilayer systems, and ethanol was added to the double bilayer systems in the bulk phase outside the double bilayer at a molar concentration of 20 mol. % with respect to the bulk phase using Packmol.⁴⁷ The CHARMM general force field (CGenFF) was used to represent the parameters for ethanol.⁴⁸ In single bilayer systems, ethanol was added after 30 ns post-equilibration to the top and bottom of the lipid membrane at six different concentrations. All equilibration and production simulations were performed under the NPT ensemble at 298.15 K. A summary of the systems is shown in Table I. Since the POPS and POPE single bilayer systems took longer to equilibrate, we extended their simulation time to 600 ns. For single bilayers, the last 200 ns of simulation was used for analysis and permeability calculations.

To find the diffusion and free energy profiles of single bilayers, we used the diffusioncma package.²⁶ The code is written in python and C++ and deposited in gitlab (https://gitlab.com/Olllom/diffusioncma). Diffusion profiles were generated at 10 different lag times t_i (10 ps–100 ps with 10 ps intervals). Diffusivities D that capture diffusive dynamics are then captured by a linear fit of MLE estimates D_i vs 1/t_i by extrapolating the profiles to infinite lag times.²⁵ 

In the counting methods, the dividing surface between water and membrane is defined as the average z-coordinate of phosphate groups with respect to the center of mass of the membrane. The permeability is calculated from the number of crossing events. The code to calculate permeability from TB counting is free available on gitlab (https://gitlab.com/Olllom/rickflow). A sample CSV file of permeation events is given in Table S1 of the supplementary material. For FB counting, we used a CHARMM script to count the number of ethanol molecules in the outer and inner water compartments, as defined by the dividing surface. A snapshot of the single and double bilayer system is shown in Fig. 1.

Membrane permeabilities are calculated with three methods: MLE for the ISD model, FB counting for double bilayers, and TB counting for single bilayers.

The free energy and diffusion profiles are calculated using the MLE approach in Ref. 26, which closely follows a Bayesian analysis that was originally developed by Hummer.⁴⁹ In this approach, the simulation domain is discretized into multiple equidistant bins, and transitions between these bins are counted vs time. The propagator p(z, t|z₀, t₀) is defined as the probability of the permeant being in position z (normal to the bilayer) at time t, given that it was at position z₀ at time t₀. The time evolution of this propagator is described by the Smoluchowski equation. In the discretized version of the diffusion equation, a transition rate matrix R is defined, which describes the transition rates between neighboring bins,

The propagator in this model is given by

where p(i, t|j, 0) is the conditional probability that a permeant starting from bin i is in bin j after time t. This theoretical propagator can be used to extract free energy and diffusion profiles from MD trajectories. The free energy and diffusion profiles are optimized with an efficient optimization method called the Covariance Matrix Adaptation-Evolution Strategy (CMA-ES).²⁷ This method is combined with a Bayesian approach to calculate free energy and diffusion profiles and gives more robust profiles than the Monte Carlo method. The plausibility of each pair of profiles is then assessed with the log-likelihood function that is defined as

The sum is over all possible transitions between bins, and N_ij is the observed number of such transitions at lag time in the simulated trajectory. The free energy and diffusion profiles are represented by linear combinations of basis functions (cosine series) to lower the amount of sampling required. Ten basis functions for free energy and six basis functions for diffusion profiles are used. Another difficulty in the Bayesian approach is choosing the right lag time, which should be long enough to be sensitive to the diffusion process for the permeant. However, a very long lag time makes the matrix exponential in Eq. (4) indifferent to local diffusion, which can blur the local shape of the free energy profile.²⁵ Having the optimum free energy and diffusion profiles, the permeability is then calculated from the ISD model [Eq. (2)].

Another method for directly calculating membrane permeability is through counting the number of permeants crossing through the bilayer in a double bilayer system, as depicted in Fig. 2. In this system, the only way for the permeant molecules to move into another compartment is to cross a bilayer, and this approach allows a direct measurement of membrane permeability. Fick’s first law relates the crossing rate to the membrane permeability and describes diffusion between compartments across a uniform slab of thickness h,

n₁ and c₁ refer to the number of ethanol molecules and its concentration in compartment 1 (Fig. 2). In these equations, J is the molar flux and the cross-sectional area is multiplied by 2 because of the system being a double bilayer. Solving for the number of molecules results in

where N is the total number of permeants. These equations are more elaborated by Krämer et al.²⁶ 

FB counting can be used to calculate the permeability of double bilayer simulations. Venable et al. simulated water permeation in a DLPC double bilayer setup.²⁸ Based on the experimental data, they estimated the equilibrium time (for two compartments to be reasonably mixed) to be 20 µs, which is too long to generate through MD simulations on conventional clusters. On the other hand, a short time expansion of exponential in Eq. (8) results in the following equation:

Permeability from double bilayer membranes is then calculated from the slope of the number of molecules in compartment 2 vs time. Although the short time expansion is reasonable for water, solutes with higher permeability such as ethanol require a modification. To take into account the dependence of slope dn₁/dt, the short time expansion can be replaced by the following equation:

In the above equation, J is the total flux of the permeant into the inner compartment, b is the length of the outer water compartment, and l is the length of the inner water compartment, as defined by the dividing surfaces located at the two phosphate planes. The numerator in Eq. (10) is the slope of the number of permeants in compartment 1 vs time. The number in the denominator is the average over a part of the simulation. Figure 3 shows the average number of ethanol molecules in the inner compartment vs time for POPE, POPS, and POPC. Figure 3 shows that for POPC after 1000 ns, about 250 ethanol molecules crossed the phosphate planes and entered the inner compartment in between the two bilayers. This number is about 100 for POPS and 125 for POPE membranes. Noting the total number of ethanol molecules in the system, which is 1250, the concentration of ethanol for POPC is 16% in the outer compartment and 4% for the inner compartment after 1 µs.

FB counting is not always practical to extract permeability information of bilayers from simulations. After the equilibrium is reached, the flux will be zero, and the linear fit is no longer applicable to fit the model. Moreover, Venable et al. demonstrated that for permeants that are highly soluble in the membrane and have a low concentration in the solvent, one would need a very large system to simulate permeation of the molecules.²⁸ In their paper, Venable et al. showed that these problems are solved by determining the crossing rate that is directly related to the permeability. They showed that the permeability can be calculated from

where r is the transition rate through the membrane in both directions and c_w is the average concentration of the permeant in the inner and outer compartments. This method assumes that permeation events for each particle are equally probable, crossings are independent, and the rate of crossing is constant.

The free energy profiles of different membranes at different ethanol concentrations were derived from MLE. Figure 4 shows the free energy and diffusion profiles at 8% ethanol concentration for the POPC simulation. The free energy and diffusion profiles for the POPC membrane at different ethanol concentrations are presented in Fig. 5. Figures S1 and S2 of the supplementary material show the profiles for POPE and POPS membranes.

The PMF for all the membranes show a barrier near the phosphate region. There is a local maximum and minimum in the free energy profiles below the phosphate headgroups around 20 Å from the membrane midplane. This means that ethanol has a favorable location (minimum) below the negatively charged phosphate plane of the membrane and an unfavorable location (maximum) in the membrane midplane. The barrier height near the headgroup region as well as the membrane midplane for three different headgroups follows the general order POPS > POPE > POPC, although for lower ethanol concentrations, the free energies for the POPS and POPE membrane are similar. POPS lipids possess a negatively charged headgroup, and the interfacial peak is about 30% greater than the POPC bilayer. Furthermore, the two local minima in the headgroup region are further apart in the case of POPS than POPE and POPC. This reflects that the POPC membrane is thinner with a higher surface area than POPE and POPS membranes, likely contributing to greater permeability.^50–53 

The free energy barrier for the transfer of the permeant from water into the membrane (ΔG_{water→lipid}) shows strong similarity to the previous estimated values from simulations. This can be estimated by the highest peak in the PMF and can be compared with experimental transfer free energies of short-tail alcohols from water to n-hexane. The experimental value for transferring ethanol from water to n-hexane is 2.73 kcal/mol (4.6k_BT), which shows similarity for the calculated free energies at 6% ethanol concentration for the POPC membrane.^54,55 Whitting et al. measured the thermodynamic of transfer of solutes from water to hexadecane. Their value for ethanol transfer from water to hexadecane is 2.99 kcal/mol, which is close to the measured value for n-hexane.⁵⁶ Chipot and Comer⁵⁷ adopted a fractional Smoluchowski equation to account for sub-diffusion and the ABF method to calculate the PMF of ethanol across the POPC bilayer and the barrier height was 2.9 ± 0.2 kcal/mol, which is very similar to the value calculated here 3 kcal/mol. The shape of the free energy profile here is also similar to the profiles calculated using the fractional Smoluchowski equation and ABF method for computing the free energy profile.⁵⁹ 

Although the free energy profiles from MLE approach at low ethanol concentration are similar to the experimental values and other computational works, at high ethanol concentrations the error in free energy and diffusion profiles is high, and the ISD model is not accurate for the calculation of permeability. Figure S3 compares the free energy and diffusion profiles for the POPC membrane at 2% and 18% ethanol concentrations. At the higher concentration, there is about 1 k_BT error in estimation of free energy. This huge error makes the ISD model inappropriate at higher concentrations. The reason for this failure at a higher permeant concentration is the breakdown of the Markovian assumption for ethanol. This effect is elaborated in Ref. 26. In summary, ethanol molecules fill up the space below the lipid headgroups where they form a hydrogen-bond network with themselves as well as the lipid headgroups. Permeation involves breaking out of this hydrogen-bond network, which is energetically unfavorable. This effect is stronger at higher ethanol concentrations, where there is a stronger hydrogen-bond network below the lipid headgroups, and the Markovian assumption is no longer valid. Krämer et al.²⁶ showed this effect in their Fig. 4. The free energy profiles are also calculated from the histogram count in the TB counting method. Figure 6 shows the free energy profiles for POPC and POPE membranes at different concentrations from the distributions in the counting method. The free energy profiles from the counting method are less accurate in lower concentrations due to the small number of membrane crossings.

It is observed that upon increasing the ethanol concentration for the POPC membrane, the barrier height for transfer to the center of the bilayer decreases from 6.3 ± 0.1k_BT for 1% ethanol to 3.3 ± 1.0k_BT for 18% ethanol concentration, as shown in Fig. 5(a). Moreover, at high ethanol concentrations, the dip near the headgroup region of POPC and POPE disappears. Electron density profiles (EDPs) for single bilayers of POPS, POPE, and POPC are shown in Figs. S4–S12 of the supplementary material. Figures S4, S7, and S10 of the supplementary material show the EDP for ethanol in POPS, POPE, and POPC single bilayers, respectively. The hydrophobic core of the bilayer has negligible ethanol density. There is a shoulder in the EDP for ethanol around 15 Å, which corresponds to the dip in the free energy profile near the headgroup region, which confirms that this is a favorable location for ethanol molecules. At a high ethanol concentration (18%), it is observed that the headgroup region is saturated and the dip in the EDP near the headgroup region disappears for the POPS and POPE membranes but not for the POPC membrane (see Fig. S10 of the supplementary material). The headgroup of the POPC membrane has a higher capacity for ethanol, which relates to its higher surface area per lipid. Figure S5 of the supplementary material shows the effect of the ethanol concentration on the EDP of phosphates in POPS. It is evident that a higher ethanol concentration causes the peaks for phosphates to become wider and closer to each other. This means that ethanol is thinning the membrane, resulting in closer phosphate headgroups. Similar effects are seen for POPE and POPC membranes in Figs. S8 and S11 of the supplementary material, respectively. The total EDPs for POPS, POPE, and POPC single bilayers are shown in Figs. S6, S9, and S12 of the supplementary material, respectively. At the highest ethanol concentration (18%), the shape of the total EDP changes and the shoulder near the carbonyl group disappears. Figure S13 of the supplementary material shows the component EDP for POPC, POPE, and POPS double bilayer systems. The phosphate peak of the inner leaflet is narrower compared to the outer leaflet. This is because the outer leaflet is exposed to a higher alcohol concentration than the inner leaflet. There is also a local minimum in the EDP of ethanol near the phosphate headgroup for POPS and POPC lipids, which disappears for the POPE membrane. This demonstrates the lower affinity of the headgroup of POPS and POPC than POPE for ethanol and more exclusion near the phosphate group. Figure S14 of the supplementary material shows the surface area per lipid for double bilayer systems. POPC has the largest surface area, and POPE and POPS lipid membranes have a similar surface area per lipid, and this correlates with the EDPs and membrane thickness that is more influenced for the POPC membrane.

The reason for this difference in these three membrane systems is the difference in the surface area per lipid in these systems. POPC has the highest surface area, making the free energy barrier for permeation the smallest. POPS is also a charged lipid that causes a high barrier near the headgroup region for ethanol molecules.⁵ The choline group in POPC is important for its permeability. It sterically hinders the ethanol from leaving the membrane, allowing ethanol to cluster within the membrane. This clustering perturbs the membrane, exposing the hydrophobic core to solvent and increasing the permeability. Similar effects were observed by Konas et al.⁵ 

MLE is used to calculate the diffusion profiles of membranes. Figures 4(b) and 5(b) plot the normal diffusion profiles for the POPC membrane. The profiles were generated for 10 different lag times, and diffusivity was obtained by extrapolating the profiles to infinite lag time.⁴⁶ While some modeling errors are to be expected from the collective behavior of ethanol, the diffusion profiles still illuminate the kinetics in different regions along the z-coordinate. Figure 4(b) shows that diffusivity drops near the headgroup region and has the minimum value in the membrane midplane. Figure 5(b) shows the diffusivity profile for POPC at different ethanol concentrations. Ethanol diffuses faster in solution compared to the membrane. The self-diffusion coefficient of ethanol is lower than the diffusion coefficient of ethanol in water. Moreover, the ethanol diffusion coefficient in the water phase was faster than ethanol self-diffusion and decreased as a function of ethanol concentration. This finding is consistent with a previous MD study of water/ethanol mixtures in the absence of a bilayer.⁵⁸ Therefore, at higher ethanol concentrations, self-diffusion of ethanol molecules causes the peak in the headgroup region to decrease. Seydel et al.⁵⁹ used quasi-elastic neutron spectroscopy (QENS) and measured the ensemble-averaged self-diffusion coefficients of water and ethanol as a function of water–ethanol mixing ratio on the picosecond timescale. The diffusion coefficients in these experiments were decreasing from 1.6 × 10⁻⁵ cm²/s in pure water to about 0.6 × 10⁻⁵ cm²/s in 20 mol. % ethanol at 285 K. The diffusion coefficient calculated from simulation data for the lowest ethanol concentration in the TIP3P water is 2.1 ± 0.1 × 10⁻⁵ cm²/s at 298.15 K, which is higher compared to the experimental diffusivity in the water phase (∼1.5 × 10⁻⁵ cm²/s). At the highest ethanol concentration in our simulations, the computational diffusivity of ethanol in the water phase is more than two times greater than that of the experimental value (∼0.7 × 10⁻⁵ cm²/s). The difference between the experimental and simulated values is because of errors in viscosity (force field) and temperature. In Ref. 26, the effect of the friction constant γ = 1 ps⁻¹ in the Langevin thermostat was studied and compared to the Nosé–Hoover thermostat for computing diffusivity of ethanol in the membrane and solvent. The difference in the membrane was small suggesting that simulations with 1 ps⁻¹ friction coefficient produce acceptable results. However, it is more advisable to use a thermostat that does not affect the diffusive behavior of the solute. Moreover, the water baseline for diffusivity is too high because of the parameterization of the TIP3P-water model, resulted in an overestimate in the diffusivity.⁶⁰ 

Figure 7 shows permeabilities of three different membranes (POPC, POPE, and POPS) studied at six different ethanol concentrations ranging from 1% to 18%. An assumption for the ISD model is that mass transport is described by the Smoluchowski equation, which is assumed that the permeant follows a pure diffusive behavior across the membrane. The ISD model assumes that the permeant obeys classical diffusion such that the probability to find a permeant at position z₂ at time t₂ given that it was at position z₁ at time t₁ follows a Gaussian distribution in an isotropic environment.¹² Anomalous diffusion or non-diffusivity can come from memory effects, inertial effects, and a complex medium. To describe these effects accurately, one needs to extend the Smoluchowski equation to the generalized Fokker–Plank equation⁶¹ or the fractional Smoluchowski equation.

Figure 7 shows that the calculated ethanol permeabilities follow the order POPC > POPE > POPS. Ethanol permeability of membranes increases with ethanol concentration, and the trend is linear for POPC and POPE membranes. However, for the POPS membrane, although a positive trend is generally observed, the concentration dependence is less pronounced due to the membrane’s resistance to permeation. This is due to the low surface area of the POPS lipid and the charge on the headgroup, which makes it less permeable to ethanol. The error bars in the estimation of permeability are larger at a higher ethanol concentration where the MLE approach fails to capture the right free energy and diffusion profiles due to incompleteness of the reaction coordinate, which arises from the non-Markovian behavior of ethanol described before.

Figure 8 shows permeabilities for different concentrations of ethanol in the POPE, POPC, and POPS membranes based on the TB counting method. The error bars for permeability are large at low concentrations and small at higher concentrations. This is because the number of transitions increase at higher concentrations, which gives better statistics for the permeability calculation. The POPC membrane has the largest area fluctuation among the three lipid bilayers studied here, which is the reason for larger error bars in the permeability calculation for POPC. Table S2 of the supplementary material shows the result of TB counting for the POPE membrane. In each concentration, the permeability is calculated from the number of crossings based on the last 200 ns of simulation. Error bars are from 95% confidence interval, which are calculated assuming a Poisson distribution. Permeabilities calculated from the ISD model are about three times larger than the ones calculated from the counting method at higher concentrations. However, at low concentrations, they are within the same error range as the counting method. The reason for this is in the assumptions in the ISD model. The reaction coordinate for the calculation of free energy and diffusion profiles in the ISD model is the z-distance normal to the bilayer plane. This is the natural transition coordinate at low concentrations, where the ethanol kinetics are largely decoupled, which means they transition the bilayer independent of one another. The kinetics become less independent at higher concentrations due to the formation of hydrogen bonds. In effect, the transition coordinate becomes highly nontrivial and can no longer be described in terms of single molecules.

Figure 9 shows the probability distribution of angles and demonstrates the rotation of ethanol when crossing the membrane. In the bulk phase, the angle distribution is wide from 60° to 120°. When approaching the headgroup region, ethanol rotates to have an angle of 30°–50° near the headgroup of the membrane. Due to less sampling in the hydrophobic core of the membrane, the distribution in this region is noisy. However, it is evident that ethanol molecules rotate 180° in the hydrophobic core region. The angle distribution does not change with the ethanol concentration. Moreover, the hydration number of ethanol molecules when permeating the membrane is shown in Fig. S15 of the supplementary material. The hydration number for ethanol in the bulk phase decreases at higher ethanol concentrations. Inside the membrane, different ethanol concentrations follow the same behavior with the highest hydration number at the phosphate region and the lowest at the bilayer center.

In another approach, we simulated double bilayer membranes, which allows an estimation of membrane permeability from counting the number of ethanol molecules crossing the membrane and moving into the inner water compartment. Fick’s law of diffusion is used to calculate the membrane permeability. The results of calculating permeabilities for three different bilayers at a starting ethanol concentration of 20% in the bulk phase outside the bilayer are shown in Table II, where permeabilities are calculated using the short time expansion and from the definition of permeability in Eq. (10) selecting all the data and the first and last 500 ns blocks. The choice of the dividing surface for the calculation of permeability is of great importance and can impact the permeability calculations. Based on the definition of permeability, we decided to define the dividing surface as the phosphate plane of the inner leaflet. It is observed that the permeability for the second block (500 ns–1000 ns) is greater than the first block (0 ns–500 ns) for POPE and POPC lipids, which means that the permeability is increasing with the simulation time. In the POPS membrane system, however, permeabilities between the time ranges are statistically equivalent, which is due to the low permeability of the POPS membrane.

For single bilayers, the transition rate, r, is the number of crossings in both directions divided by the length of simulation and cross-sectional area of the bilayer. TB counting can also be used with double bilayers to calculate the permeability with an anisotropic concentration of the permeant on either side of the membrane. For double bilayers, r is the number of crossings over any bilayer in any direction divided by the length of simulation and the total membrane area multiplied by two. At short timescales, TB counting is equivalent to FB counting. Table III shows the results for TB counting in the double bilayer systems. Permeabilities calculated from TB counting are slightly greater than FB counting, which is due to a higher number of crossings. In TB counting, the simulation is divided into regions of 10 ns, and the permeability is calculated by counting the number of crossings in a 10 ns window and accumulating the data for subsequent windows. Equation (11) shows that the permeability is affected by both the surface area and concentration of the permeant in the water compartment, and therefore, permeability can change through time. Figure 10(a) shows that c_w (average concentration of the permeant in the inner and outer compartments) increases linearly with time. However, permeability also increases with time and the reason is the increase in transition rate [see Fig. 10(c)] with the surface area of the membrane. Figures S16 and S17 of the supplementary material show the same behavior for POPE and POPS membranes. As the permeants cross the outer leaflet, the concentration of permeants below the headgroups of the inner leaflet increases, but since the membrane does not reach equilibrium, the concentration of the permeant in the outer leaflet is higher than the inner leaflet. Therefore, the defect and deformations are more pronounced in the outer leaflet. This effect cannot be seen in a single bilayer setup since both leaflets feel the same concentration of the permeant.

In real biological systems, the concentration of the permeant inside and outside the membrane is different and cells can experience an anisotropic effect of ethanol across their membrane. The results from single bilayer simulations are used to calculate the permeability of membranes where the concentration on the two sides of membrane is different. This cannot be done in single bilayer simulations because of the periodic boundary conditions. In double bilayer systems, one can have different concentrations of the permeant on different sides of the bilayer and calculate the permeability in this anisotropic system. This can be challenging due to the computational cost of simulating a double bilayer membrane, and the timescale of simulation should be long enough for permeants to cross the bilayer. Another approach is to simulate single bilayer membranes at different concentrations and use the solubility-diffusion model to calculate free energy and diffusion profiles. These profiles can further be combined in a compartmental model with two compartments to calculate the permeability of a membrane with different concentrations on either side. The permeability in this method can be calculated by summing up the resistivities of the two leaflets or, equivalently,

Table IV shows the calculated permeabilities of a lipid membrane of POPE using permeabilities from TB counting in single bilayer membranes at different concentrations. Comparisons with double bilayer results will be made in Sec. VI.

Experimental membrane permeability measurements are done using the micropipette aspiration technique. In this experiment, area expansion measurement of vesicles exposed to flowing alcohol solutions provide a direct measurement of membrane permeability.⁵⁰ However, there is a huge difference between theoretical and experimental membrane permeability to ethanol. The experimental permeability of the SOPC membrane to ethanol at 1.70M concentration is reported to be 3.8 ± 3 × 10⁻⁵ cm/s, whereas the smallest theoretically calculated permeability for the POPS membrane to ethanol is on the order of 10⁻² cm/s. Experimental membrane permeabilities to ethanol are scarce, and theoretical membrane permeabilities are on the same order of magnitude as in this work. The permeabilities calculated in this paper for the POPC membrane at 1%–18% ethanol from the counting method ranges from 0.170 cm/s to 0.358 cm/s. Ghaemi et al. used bias exchange metadynamics simulation to estimate the ethanol permeability of the POPC bilayer and inferred permeability from a flux-based rate model with the GROMACS force field to be 6 × 10⁻² cm/s. Chipot and Comer used the fractional Smoluchowski equation and adaptive biasing force (ABF) to calculate the permeability of the POPC membrane to methanol with the CHARMM36 force field and obtained a value of 0.158 cm/s.

To explain the discrepancy between experimental and theoretical estimates of membrane permeability to ethanol, one might consider several reasons. The difference in the lipid chain used in experiment (SOPC) and simulation, either in its chemical nature or its length, is unlikely to explain three orders of magnitude difference between theoretical and experimental permeability. In the experimental setting, an ethanol-free vesicle is placed in a constant flow of ethanol, and the area expansion of the vesicle is monitored through time. Permeability is then calculated from the area expansion of the vesicle vs time. Ghaemi et al. calculated the partition coefficient of ethanol to be 6.5 × 10⁻³, and by setting the experimental values of K, P, and h in the equation P = (K · D)/h, they calculated a value on the order of 10⁻⁸ cm² s⁻¹ for the diffusion coefficient of ethanol in the membrane. This means that ethanol must diffuse about 1000 times faster in water than in the membrane, which is shown to be unlikely for molecules at the same size. Ghaemi et al. showed that to reach the stationary condition when the concentration inside and outside the liposome is similar, ethanol molecules must not only traverse the membrane but also diffuse inside the liposome. They calculated the relaxation time of ethanol concentration at the center of a liposome modeled as a spherical container having the same radius as in experiment. This relaxation time is the timescale of diffusion of ethanol molecules in water until equilibrium. The diffusion timescale is calculated to be a few seconds consistent with experiments, which suggest that the area expansion in experiment does not only reflect the translocation time across the membrane and they should also consider the time required for the permeant to equilibrate within the vesicle.²⁰ On the other hand, they argued that in the experimental permeability calculation, one needs to consider not only crossing of the permeant inside the membrane but also diffusion inside the vesicle. Since the diffusion and relaxation timescales for ethanol are very close, interpretation of experimental results is difficult due to uncertainty in the rate limiting step for permeation. In addition, it was shown that the unstirred layer effect (USL) plays a major role in permeation of small molecules through the membrane.⁶² USL consists of a few layers of ordered water molecules. Transport through the USL imposes an additional resistance to diffusion of small molecules through the membrane. The size of the USL is comparable to the size of giant unilamellar vesicles (GUVs) that are used for experimental permeability calculations for ethanol, and therefore, the time for diffusion into these layers cannot be ignored in experimental permeability calculations. Hannesschlaeger et al.⁶² calculated the permeation across USL for CO₂ considering a planar bilayer with a diameter of 150 µm, a USL thickness of 100 µm, and a solute diffusion coefficient of 2 × 10⁻⁵ cm/s, which resulted in CO₂ permeability of ∼10⁻³ cm/s for the USL. The membrane permeability to CO₂ is shown to be orders of magnitude higher than the USL permeability with a value of about 20 cm/s.⁶² This value is roughly equal to the CO₂ permeability of only a 10 nm USL. Therefore, diffusion across the USL is the rate limiting step for CO₂ permeation across the membrane and cannot be neglected. However, the situation is less dramatic for ethanol since the hydrophobic core of the membrane has the highest energy barrier for permeation of ethanol. On the other hand, CO₂ does not experience a high barrier in the membrane midplane, and the phosphate headgroup region presents the highest free energy barrier for permeation. This infers that even a 100 μm USL does not completely explain the permeability on the order of 10⁻⁵ cm/s for ethanol. However, the USL still plays a significant role, and its effect should be carefully considered in experimental permeability calculations.

The results of our simulation show that the free energy barrier for permeation of ethanol is very close to the experimental solvation free energy of ethanol in water and n-hexane, which rules out the inaccuracy of the force field in reproducing this thermodynamic property. 5⁵¹ The TIP3P model of water exaggerates self-diffusion, which could be the reason for overestimation of diffusivity in the bulk water phase.⁵⁰ The inaccuracy of the reaction coordinate in the ISD model at higher ethanol concentrations is another reason for the difference between permeabilities calculated from simulation and experiment using the ISD model. A small system size can adversely affect the permeability calculations because deformations due to permeation are precluded in such small system sizes.⁶³ Moreover, the effect of periodic boundary conditions is huge for lipid and protein diffusion in a small system size. However, it was shown that the effect of system size is small for diffusion of small molecules such as water and O_2.

In this paper, we compared two different frameworks for the calculation of membrane permeation: (1) an ISD model using a MLE approach to get free energy and diffusion profiles and (2) two different counting (flux-based and transition-based) methods. The free energy and diffusion profiles are extracted from MD simulation trajectories described by the Smoluchowski equation for diffusion of molecules in a free energy surface. The free energy barrier against permeation of ethanol for a series of lipid bilayers of POPE, POPS, and POPC is obtained. Consistent with their area per lipid, the ethanol permeation of membranes follows the trend of area per lipid with POPC > POPE > POPS. The ISD model is used to calculate the permeability of single bilayer membranes, and the results show that the membrane permeability increases with increasing ethanol concentration. The TB counting method is used for single bilayer membranes where the number of transitions is measured through time to estimate the membrane permeability. The ISD model overestimates the permeability for ethanol, which is due to the breakdown of the Markovian assumption in diffusive modeling of ethanol (see Ref. 26 for a more detailed discussion). This effect is more pronounced at higher ethanol concentrations where the network of hydrogen bonds between permeating ethanol is stronger. However, this effect is small in lower permeant concentrations, and the permeabilities obtained in this regime have lower error bars for the POPC membrane and are comparable to the results from the TB counting method. Krämer et al.²⁶ showed that water and oxygen molecules do not have this collective behavior, and the ISD model is accurate.

Considering the work here and in Ref. 26, the non-Markovian permeant behavior resulting in an inaccurate ISD model prediction is caused by two properties: first, permeating the membrane with a relatively low energy barrier and second, being able to form a network of hydrogen-bonds with itself and with membrane headgroup. Water and oxygen do not satisfy these two properties, and the ISD model is accurate. When considering the permeation of other produced biofuels and chemicals or drugs, the use of ISD should consider if non-Markovian behavior will exist due to the high permeability and ability to form a hydrogen-bond network. The TB counting method should not be used at low permeant concentrations due to insufficient sampling, while the ISD model should not be used for ethanol-like molecules. When the ISD model fails, the TB counting method can be used, as it does not rely on the definition of a reaction coordinate.

To represent permeation in a cell-like environment, we used microsecond simulations of double bilayer membranes to calculate the permeability using Fick’s law of diffusion in the FB-counting method where the concentration of the permeant in the inner and outer compartments is evaluated vs time. Due to limitations in FB counting, we used transition-based counting to calculate the permeability of single and double bilayer membranes. Figure 10 shows that due to the increase in the surface area of the membrane vs time, the transition rate and subsequently membrane permeability also increase with time when considering cumulative blocks of 10 ns for calculating permeability. The short time expansion of Fick’s first law is replaced by the definition of permeability in Eq. (10), which allows using compartments of data for the calculation of permeability from Fick’s law. By simulation of the double bilayer membrane in a regime before equilibrium is reached, one can divide the exponential function of the number of ethanol molecules vs time by multiple linear equations that correspond to the concentration difference at the particular time interval and calculate permeability using Eq. (10).

Subsequently, we used the results of single bilayer permeability calculation from the TB counting method in a compartmental model to calculate the permeability of membranes with the anisotropic concentration of permeants. For the POPC double bilayer, the second block of data that corresponds to 16% ethanol concentration in the outer compartment and 4% in the inner compartment results in a permeability of 0.233 ± 0.021 cm/s, and the result from the compartmental model for the same concentration is 0.236 ± 0.029 cm/s. As seen, the results closely match for permeation in the POPC membrane. For the POPE membrane, the permeability from double bilayer FB counting is 0.120 ± 0.010 cm/s, which is very close to the result from TB counting and shows remarkable similarity to the permeability from the compartmental model at the same concentrations to be 0.124 ± 0.019 cm/s. Moreover, the speed of a simulation scales roughly as N log N with the number of atoms in the simulation, and doubling the size of the system from 30 000 atoms to 60 000 can increase the simulation time by a factor of about 2.1. Therefore, to avoid computational cost of the simulating double bilayer membranes, one can replace them by a single bilayer at different permeant concentrations. Although simulating double bilayer membranes for the permeability calculation is expensive, there are still benefits to do these simulations. For a single bilayer, c_w does not change with time. In contrast, Fig. 10 shows that c_w for a double bilayer increases with time. On the other hand, in a double bilayer, simulation inner and outer leaflets are exposed to different concentrations of the permeant, which cannot be done in a single bilayer setup. Furthermore, if one is interested in calculating the permeability of membranes with asymmetric lipid composition on leaflets, one needs to simulate a double bilayer membrane where the inner and outer leaflets have different lipid compositions.

The calculation of membrane permeability to small molecules is an inherently challenging task due to the complex and inhomogeneous environment of membranes. The MLE approach and ISD model can give inaccurate permeabilities if they are used for a high concentration of amphiphilic permeants due to the breakdown of the Markovian assumption for diffusive modeling. One can instead use the TB counting method if there are considerable crossings over the bilayer. By combining permeabilities from single bilayers, we showed that the permeability in this compartmental model is similar to the permeability calculated from a double bilayer at the same concentration. Permeabilities calculated from the counting method for the POPC membrane ranges from 0.170 cm/s to 0.358 cm/s for 1%–18% ethanol. The calculated permeabilities in this work and others from simulation are orders of magnitude higher than the past experimental values but are likely due to the poor interpretation of experimental results based on the general convergence of permeation values from simulation with varying force fields. However, there are still opportunities for the development of permeability calculation using the counting method and ISD model. The present work is useful for the general community interested in permeability calculation for choosing which method is best for their system of interest.

See the supplementary material for additional tables and figures. Multiple tables are presented for permeabilities of POPC, POPE, and POPS single bilayers for different concentrations of ethanol. Free energy and diffusion profiles are given at different ethanol concentrations for the three bilayers. Structural information about the bilayers are presented as electron density profiles and surface area per lipid. A profile for the hydration number of ethanols based on the position relative to the membrane center is presented. Several plots for the permeability of double bilayers are given.

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

This research was supported by the NSF (Grant No. CBET-1604576; J.B.K. and M.G.) and Intramural Research at NIH (A.K., E.W., and M.G.). The high-performance computational resource used for this research is Deepthought2, which is maintained by the Division of Information Technology at the University of Maryland. Production runs used the Extreme Science and Engineering Discovery Environment (XSEDE) allocation under Grant No. MCB-100139 on Stampede2 and Bridges, and the computational resources at the Maryland Advanced Research Computing Center (MARCC).