How Ethanolic Disinfectants Disintegrate Coronavirus Model Membranes: A Dissipative Particle Dynamics Simulation Study.

We have developed dissipative particle dynamics models for pure dipalmitoylphosphatidylcholine (DPPC), dioleoylphosphatidylcholine (DOPC), and dimyristoylphosphatidylcholine (DMPC) as well as their binary and ternary mixed membranes, as coronavirus model membranes. The stabilities of pure and mixed membranes, surrounded by aqueous solutions containing up to 70 mol % ethanol (alcoholic disinfectants), have been investigated at room temperature. We found that aqueous solutions containing 5-10 mol % ethanol already have a significant weakening effect on the pure and mixed membranes. The magnitude of the effect depends on the membrane composition and the ethanol concentration. Ethanol permeabilizes the membrane, causing its lateral swelling and thickness shrinking and reducing the orientational order of the hydrocarbon tail of the bilayer. The free energy barrier for the permeation of ethanol in the bilayers is considerably reduced by the ethanol uptake. The rupture-critical ethanol concentrations causing the membrane failure are 20.7, 27.5, and 31.7 mol % in the aqueous phase surrounding pure DMPC, DOPC, and DPPC membranes, respectively. Characterizing the failure of lipid membranes by a machine-learning neural network framework, we found that all mixed binary and/or ternary membranes disrupt when immersed in an aqueous solution containing a rupture-critical ethanol concentration, ranging from 20.7 to 31.7 mol %, depending on the composition of the membrane; the DPPC-rich membranes are more intact, while the DMPC-rich membranes are least intact. Due to the tight packing of long, saturated hydrocarbon tails in DPPC, increasing the DPPC content of the mixed membrane increases its stability against the disinfectant. At high DPPC concentrations, where the DOPC and DMPC molecules are confined between the DPPC lipids, the ordered hydrocarbon tails of DPPC also induce order in the DOPC and DMPC molecules and, hence, stabilize the membrane more. Our simulations on pure and mixed membranes of a diversity of compositions reveal that a maximum ethanol concentration of 32 mol % (55 wt %) in the alcohol-based disinfectants is enough to disintegrate any membrane composed of these three lipids.

Introduction

The infectious respiratory coronavirus disease 2019 (COVID-19) caused by the 2019 novel coronavirus (2019-nCoV), also known as SARS-CoV-2 and HcoV-19, has spread throughout the whole world.[1] The outer layer of the coronavirus envelope is composed of the membrane (M), spike (S), and envelope (E) proteins and the host-derived lipid bilayer, which gives the virus its distinctive shape and structure and protects its RNA from the surrounding environment.[2−5]

Recently, there has been rapid progress in the development of safe and effective vaccines against the coronavirus and development of potential therapies for SARS-CoV-2.[6−9] However, still a considerable fraction of the world’s population is unvaccinated, the new variants of the virus may spread, and there is no established treatment for SARS-CoV-2 infection. Washing one’s hands with soap and water or hand sanitizer that contains at least 60–70 wt % alcohol (usually ethanol, n-propanol, isopropyl alcohol, or a mixture of them) is still one important way to prevent the spread of COVID-19 and other corona viruses.[10,11] Since the viral membrane acts as a barrier to the penetration of small molecules through it, the deactivation of the virus is primarily controlled by its membrane permeability to alcohol. By using different techniques, investigators have studied the alcohol-induced changes on various lipid bilayer systems.[11−14] Concentrated alcohol solutions increase the area per lipid molecule, accompanied by a decrease in the bilayer thickness along with disordering and enhanced interdigitation of lipid acyl chains.[15,16] These changes cause a loss of membrane integrity and make it permeable to the passage of alcohol molecules, water, and other species, a process that eventually leads to membrane rupture.[17] In a recent study, we have reported that the stability of the pure dipalmitoylphosphatidylcholine (DPPC) membrane against alcoholic disinfectants strongly depends on the phase of the membrane.[15] We reliably observed the disintegration of the DPPC membrane in its liquid crystalline phase (323 K), at ethanol concentrations ≈15 mol % in the aqueous phase surrounding the membrane, while its gel phase (298 K) remained intact even at higher ethanol concentrations (up to about 20 mol % examined).

Although the exact composition of the viral membrane is unknown and presumably changes between individual virus particles, we know at least that it contains a mixture of lipids that mechanically anchor the S- and E-proteins. Among all lipid types, phosphatidylcholines (PCs) are the most important lipid components of living organisms. Specifically, they are the main components of the endoplasmic reticulum Golgi intermediate compartment (ERGIC), where coronaviruses are replicated and assembled.[18,19] Additionally, the lung, the primary organ affected by the coronavirus, mostly uses dipalmitoylphosphatidylcholine (DPPC), i.e., one of the PCs, as the abundant constituent of its surfactants.[20] As well as PCs (∼50%), the ERGIC of a mammalian cell contains smaller amounts of (∼15–25%) phosphatidylethanolamines (PEs) and (∼10–15%) phosphatidylinositols (PIs).[21] The PCs-, PEs-, and PIs-molecules have different head groups but the same hydrophobic tails,[21,22] which consist of saturated and/or unsaturated acyl chains of various lengths. Previous experimental[23−26] and simulation[15] reports show that the phase of the membrane, which in turn depends on the hydrocarbon tail length and its degree of saturation, is the main factor determining its stability. Our results in this work also confirm such a trend for the stabilities of membranes made up of PCs as well of PC mixtures (see the section Membrane Failure). We have compared experimental data[27] on the phases (stabilities) of PC-, PE-, and PI-membranes, of different hydrocarbon tail lengths with a different number of unsaturated bonds, in Supporting Information Figure S1. Based on these data, we would argue that the lipid headgroup does not have a dominant role in the stability of the membrane. Besides, our previous atomistic simulation results show that the largest free-energy barrier for the passage of small molecules is observed very close to the center of the membrane at both low and high ethanol concentrations.[15] Its height depends on the membrane thickness, i.e., the gap that a penetrant needs to cross.[15] It has been found experimentally[28] that the thickness of the membrane depends nearly linearly on the length of hydrocarbon tail. In contrast, the type of the headgroup does not have a noticeable influence on the membrane thickness, which is related to the stability. The headgroup is found to act only as a secondary barrier to membrane penetration, at low ethanol concentrations. Thus, the type of headgroup determines how the first few ethanol molecules partition into the membrane and their accumulation around the headgroup region. At high ethanol concentrations, penetration to the tail region causes a substantial lateral swelling of the membrane, preceding its rupture. Therefore, it is the hydrocarbon tails and their composition, which ultimately controls the membrane stability. This argument is also in line with previous simulations[29,30] of the stability of membranes of pure palmitoyloleoylphosphatidylcholine (POPC, PCs) and phosphatidylethanolamine (POPE, PEs) at high ethanol concentrations. Due to the large number of lipid constituents of ERGIC membranes (whose composition resembles that of the viral coating membrane, since the viral genome does not provide for lipid manufacture), we are forced to only simulate a subset. However, the arguments above allow us to reasonably justify simulating mixed-PC membranes as coronavirus model membranes, to investigate their stability against disinfectants. Hence, we selected three PCs with different classes of hydrocarbon chains, namely dipalmitoylphosphatidylcholine (DPPC, consisting of long saturated hydrocarbon chains), dioleoylphosphatidylcholine (DOPC, consisting of long unsaturated hydrocarbon chains), and dimyristoylphosphatidylcholine (DMPC, consisting of short saturated hydrocarbon chains) as the components of mixed lipid bilayers. Besides, the coronavirus membrane anchors S-proteins, which fuse with the host cell membrane and facilitate virus entrance to the cell, as well as E- and M-proteins. The presence of these proteins in the structure of the model membrane could possibly influence its stability against damage by ethanol. However, our recent atomistic study on the effect of the E-protein on the stability of a mixed lipid bilayer, palmitoylsphingomyelin (PSM) and POPC, immersed in ethanol–water mixtures, has shown that the E-protein has a negligible effect on the partitioning of water and ethanol from the aqueous phase to the lipid phase of the membrane.[31] In other words, the E-peptide has no appreciable effect on ethanol-induced viral membrane failure in the range of alcohol concentrations studied. Therefore, in this study, we focus our attention on the stability of protein-free DPPC-DMPC-DOPC mixed membranes, as models of the coronavirus membrane, immersed in water–ethanol solutions as disinfectants. It is worth mentioning that the stabilities of pure and mixed membranes, immersed in pure water, have been studied by several investigators.[32,33] However, reports on the stability of mixed membranes, immersed in aqueous solutions containing alcohol, are scarce. This is particularly true at high alcohol concentrations, i.e., common concentrations in alcohol-based disinfectants, in which membrane disintegration becomes an issue of importance.

Atomistic simulations for the compositionally complex lipid mixtures are computationally expensive. Here, we employed a relatively affordable computer simulation method, dissipative particle dynamics (DPD), which, if carefully parametrized, can provide meaningful results. The DPD interaction parameters were chosen from the well developed four-to-one coarse-grained (CG) mapping scheme of the MARTINI-like models,[34,35] in which four heavy atoms and their attached hydrogen atoms were mapped to one DPD bead. To validate our DPD models, we compared the structural and thermodynamic properties of lipid membranes of different compositions, surrounded by pure water, with those available from simulation and experiment. For the ethanol-containing systems, we further examined these properties to capture the effect of ethanol on the membrane stability. In addition to the elucidation of the mechanism of a disinfectant’s influence on virus deactivation, we believe that the present study provides insight into numerous other applications such as drug delivery, anesthesia, and cryopreservation, where high concentrations of alcohols are used to modulate functions of biological membranes. Comparing the DPD results with our previous atomistic studies of coronavirus model membranes, we also assess the predictive ability of DPD as a CG model, allowing the achievement of longer time and length scales, which are inaccessible by atomistic simulations.

Simulation Details

Model

The DPD method has been well described in the literature.[36] For its details, we refer the reader to excellent reviews in the literature.[34,37] Here, we restrict ourselves to a brief explanation of the method. The nonbonded conservative force between two DPD beads i and j, F = α(1 – r/r)e(r < r), separated by a distance r (e is the corresponding unit vector), is purely repulsive. The repulsion parameter α controls the magnitude of repulsion. The cutoff distance (r), bead mass (m), and thermal energy per one bead (kBT) are chosen as reduced units of length, mass, and energy in DPD simulation. The reduced time unit is defined accordingly as .

In this work, we have used the DPD model based on the four-to-one mapping scheme of the MARTINI model.[34,35] For water at a reduced density ρ = 3, we map four water molecules into a single DPD bead. The volume of four water molecules is 0.12 nm3. As a cube of volume r3 contains 3 water beads, the cutoff distance corresponds to 0.71 nm. At 298 K, each ethanol molecule has a volume of 0.097 nm3; therefore, the volume of 1.24 ethanol molecules (modeled as a single DPD bead) corresponds to the volume of four water molecules. For the lipid molecules, the volume of the DPD beads (see Figure ) corresponds to the volume of a single water bead. Setting the DPD repulsion parameter for water, a = 100, accurately reproduces the compressibility of water at room temperature, k–1 = 16.[36,37] The repulsion strength between the similarly charged beads of head groups was increased to compensate for the electrostatic repulsion (α = α = 110).[34] All other repulsion parameters for beads of the same type, α, were set to 100. The repulsion parameter for water–ethanol interaction, 102, was taken from the literature.[34,37] All other repulsion parameters between unlike beads, α, were determined from the Flory–Huggins χ-parameter, according to the following expression by Li et al.[34]where Δα = α – α. The Flory–Huggins χ-parameters between water and hydrocarbons and between the lipid head groups and hydrocarbons have been well discussed in the literature.[34,37] Li et al.[34] reported that these parameters can well describe the compressibility and bending rigidity of the real membranes. In this work, some modifications have been made to the parameter set reported in the literature.[34,37] The DPD repulsion parameter α for the ethanol-tail (hydrocarbon) interaction was tuned, by scanning it over the range 104 to 110 (based on the reports in the literature),[34,37] against the free energy barrier for the passage of a single ethanol molecule across the membrane, according to our previous atomistic simulation results (see the section Validation of the Model for Ethanol-Containing Systems).[15] We have summarized the final DPD repulsion parameters for interactions between all bead types in Figure (f).

Figure 1

Models of (a) DMPC (dimyristoylphosphatidylcholine), (b) DPPC (dipalmitoylphosphatidylcholine), (c) DOPC (dioleoylphosphatidylcholine), (d) water, and (e) ethanol in this work. Each bead has a mass and volume comparable to four realistic water molecules. (f) The DPD repulsion parameters for all bead–bead interactions.

The DPD beads in phospholipids are connected by harmonic bondswhere U is the bond potential, and k and r0 are the spring constant and equilibrium bond length, respectively. The bending angle potential, Uθ, is defined aswhere kθ is the force constant, and θ0 is the equilibrium angle. We have reported the values of k, r0, kθ, and θ0 in Table .[35,38] It should be noted that hydrocarbon chains of DPPC and DOPC are modeled by four connected tail beads, while those of DMPC are modeled as three tail beads. In the case of DOPC, the equilibrium 4-5-6 and 8-9-10 angles are set to 120° to mimic the unsaturated hydrocarbon chains.[35,38]

Table 1

Equilibrium Bond Lengths and Angles and Their Corresponding Force Constants for DMPC (Dimyristoylphosphatidylcholine), DPPC (Dipalmitoylphosphatidylcholine), and DOPC (Dioleoylphosphatidylcholine)

bond potential			angle potential
type	r0 (rc)	kb (kBT/rc2)	type	θ0 (deg)	kθ (kBT)
h1 – h2	0.47	512	h2 – h3 – t	180	6
h2 – h3	0.47	512	h2 – h3 – h3	120	6
h3 – h3	0.31	512	h3 – t – t	180	6
h3 – t	0.59	512	t – t – ta	120, 180	6
t – t	0.59	512

Except for the 4-5-6 and 8-9-10 angles in DOPC, for which the equilibrium t – t – t angle is 120°, the rest of the t – t – t equilibrium angles are 180°.

In the ethanol-free systems, a total number of 1152 lipid molecules were placed in the center of the xy plane of an xyz periodic box. The dimensions of the simulation box along x and y directions, L and L, were set according to the known area per lipid molecule for the phospholipids.[35,38] For the mixed membranes, L and L were set based on the additivity of the area per lipid for different components. We have tested different initial values of L = L, around the calculated area per lipid, and found that simulations of a constant number, N, of lipid molecules at constant temperature, T, and constant pressure, P, (the NPT ensemble) did not change the lateral dimensions of the simulation box noticeably. The initial box size along the z direction was fixed at 24 r for all systems. The box size in our simulations was beyond the limit to which finite size effects have been reported to influence the properties studied.[39] Three types of systems, namely pure systems (DPPC, DOPC, and DOPC), two-component systems (DPPC-DOPC, DPPC-DMPC, and DMPC-DOPC), and three-component systems (DPPC-DOPC-DMPC), were simulated. In all mixed-membrane systems, the lipids were randomly placed in each leaflet. Therefore, the compositions of both leaflets were the same, but the leaflets were not symmetric (see Figure (a)). The compositions of all ethanol-free lipid bilayers simulated in this work are summarized in Figure (b). We have simulated the afore-cited lipid bilayers, immersed initially in pure water, in the NPT ensemble to obtain relaxed planar bilayers. To simulate a tensionless membrane, the sizes of the simulation box in the lateral (xy) and normal (z) directions were allowed to change independently, by coupling them to a Berendsen barostat[40,41] (the time constants for pressure couplings were 10 tdpd). The lateral and normal components of the pressure were fixed at 89 kBT/r3, which is the same as that for bulk water at a reduced density of 3. Simulations were done for 5 × 105 steps to achieve equilibrium and for another 5 × 105 steps for data collection.

Figure 2

(a) Schematic of the arrangement of lipid chains in each leaflet in a ternary membrane. Blue, green, and red colors represent DPPC (dipalmitoylphosphatidylcholine), DMPC (dimyristoylphosphatidylcholine), and DOPC (dioleoylphosphatidylcholine) lipid molecules, respectively. The compositions of both leaflets are the same, but the membrane is not symmetric. (b) Schematic illustration of systems simulated in this work. The points represent compositions of the systems simulated in this work.

In the ethanol-containing systems, the number of water and lipid molecules in the system was the same as that of the ethanol-free systems, but ethanol molecules were added to the aqueous phase to reach the desired concentration of ethanol in water. According to the adopted mapping scheme, the mole fraction of ethanol in the aqueous phase surrounding the membrane is expressed aswhere nethanol and nwater are the number of DPD beads of ethanol and water, respectively. The factors 1.24 and 4 in eq account for the fact that in our mapping one DPD bead represents 1.24 real ethanol molecules but 4 water molecules. Thus, the ethanol mole fraction in this paper corresponds to the experimental mole fraction, not to the mole fraction of DPD beads used. The latter is much closer to the experimental volume fraction, since all DPD beads have the same size. In contrast to the solvent phase, no adjustment is necessary for the mole fractions describing the lipid compositions of the membranes. We have summarized the compositions of systems simulated in this work in Table .

Table 2

Compositions of the Water–Ethanol Phase in the Systems Simulated in This Work

mol % of ethanol	wt % of ethanol
0.00	0.00
5.19	12.3
10.3	22.7
13.3	28.1
15.4	31.8
17.1	34.6
20.2	39.4
23.7	44.3
27.5	49.3
31.7	54.4
42.0	64.9
55.4	76.1
73.6	87.7

Based on the mapping scheme adopted in this work,[34] the cutoff distance and the time step were 0.71 nm and 1.43 ps (0.01 tdpd), respectively. In order to compare our results with the reports in the literature, we use physical, rather than reduced DPD, units from here on.

Analysis Method

We use four parameters, including the area per lipid molecule, bilayer thickness, orientational order of the hydrocarbon chain, and the bending modulus, to characterize the degree of molecular perturbation caused by the ethanol disinfectant. Specifically, the area per lipid molecule, bilayer thickness, and bending modulus are useful only for the membranes that are still intact, i.e., without obvious holes. The membrane thickness l is defined as the average distance between the choline groups (h1 in Figure ) in two leaflets and the area per lipid, a, is defined aswhere N is the number of lipid molecules in each leaflet. The orientational order is defined in terms of the following second-Legendre polynomialwhere θ is the angle between a unit vector along the hydrocarbon chain and the bilayer normal unit vector (z axis), and the brackets denote ensemble average. Two hydrocarbon chains of each lipid molecule (Figure (a)–(c)) are calculated separately to the connection points (h3 groups). Goetz et al.[42] found that the value of bending modulus, κ, deduced from the analysis of the shape fluctuations of the bilayer membranes (in the tensionless state), can be expressed as the following equationwhere K is the area compressibility. In order to calculate K, we measure fluctuations of the interfacial tension, γ. The area per lipid molecule, a, is varied by modifying the lateral size in the NVT ensemble (V being the volume). Then, K can be obtained as the zero-tension limit of the slope of the interfacial tension versus (a – a0)/a0 as[42]where a0 is the area per lipid molecule in a tensionless membrane (membrane immersed in the pure water). It should be noted that this procedure for calculation of κ is only suitable for the membranes with negligible spontaneous curvature.[34] The interfacial tension γ can be obtained from the pressure anisotropy aswhere P, P, and P are the diagonal components of the pressure tensor.

Results and Discussion

Validation of the Bilayer Models for Ethanol-Free Systems

To validate the lipid bilayer model, immersed in water, we have calculated the area per lipid molecule, a0, the membrane thickness, l0, the orientation order of hydrocarbon groups, S0, and the bending modulus, κ0. As summarized in Table , our calculated a0 and l0 of DPPC, DMPC, and DOPC are in good agreement with experiments[43−46] at temperatures close to 300 K. Poghosyan and Gharabekyan[33] have investigated the DPPC-DMPC membranes of various compositions at 325 K by performing atomistic simulations. Their observed trend of decreasing a0 versus xDPPC is compatible with the trend found at 298 K here (see Table ). For DPPC, closeness of value of S0 to 1 (S0 = 0.894) indicates that the bilayer is in the gel phase, in which the hydrocarbon chains are in a highly ordered state (Figure (a)), while DMPC (S0 = 0.467) and DOPC (S0 = 0.426) bilayers are in the fluid phase, as illustrated in Figure (b),(c)). The deviations of bending moduli obtained in this work from experimental values of DPPC (∼29% to 76%),[47] DMPC (∼−16%),[48] and DOPC (∼−30% to −75%)[47] can partly be due to the softness of DPD beads in our simulations and partly due to experimental uncertainties in measuring the bending modulus. Chaurasia et al.[49] have compared several methods of evaluating the bilayer bending modulus in simulations and reported that the bending modulus of PC lipid membranes in the fluid state (325 K) was within the range of κ0 ∼ (0.4–2.1) × 10–19 J, which matches well with our calculated results (∼(0.4–0.9) × 10–19 J) for PC membranes in the liquid phase at 298 K. Because of the lack of experimental data, we cannot compare our calculated bending moduli for mixed membranes with the experiment. Based on the agreement of our calculated values of a0, l0, S0, and κ0, for membranes immersed in pure water, with the corresponding experimental values, we conclude that our DPD models are well parametrized against experimental measurements.

Table 3

Calculated Equilibrium Thermodynamic and Mechanical Properties of Lipid Bilayers, Compared with Previous Simulations and Experimentsa

ratio			this work 298 K				previous simulations			previous experiments
xDPPC	xDMPC	xDOPC	l0 (nm)	a0 (nm2)	S0	κ0 (10–19J)	l0 (nm)	a0 (nm2)	κ0 (10–19J)	l0 (nm)	a0 (nm2)	κ0 (10–19J)
1	0	0	4.76	0.576	0.894	12.496	4.71[43]a	0.487[43]a	7.1–8.9[47]a	4.35[15]a	0.50[50]a
0	1	0	3.52	0.651	0.467	0.471	3.53[51]a	0.597[52],a 0.606[51]a	0.56[48]a	3.47[53]a	0.616[53]a
0	0	1	3.89	0.706	0.426	0.557	3.67[44],a 3.71[46]a	0.724[44],a 0.674[45]a	0.8–2.2[47]a	3.82[53]a	0.695[53]a
0.25	0.75	0	3.57	0.647	0.476	0.454					0.74[33]b
0.5	0.5	0	4.38	0.627	0.556	0.851					0.72[33]b
0.75	0.25	0	4.53	0.601	0.636	5.24					0.70[33]b
0.333	0.333	0.333	4.12	0.671	0.516	0.62

Represents the reference value obtained at 298–303 K.

Rrepresents the reference value obtained at 325 K.

Figure 3

Snapshots of the simulation box, indicating pure(a) DPPC (dipalmitoylphosphatidylcholine), (b) DOPC (dioleoylphosphatidylcholine), and (c) DMPC (dimyristoylphosphatidylcholine) lipid bilayers. Blue, green, and red colors represent the hydrocarbon tails of DPPC, DMPC, and DOPC lipids, respectively.

Effect of Composition on the Properties of Membranes in the Absence of Ethanol

Figure depicts the calculated values of a0 and l0 as a function of composition for different membranes.

Figure 4

Two-dimensional plots of variation of (a) the area per lipid molecule and (b) the thickness of the bilayer as a function of composition for the ethanol-free systems at 298 K. The color bars represent the area per lipid molecule and the membrane thickness in panels a and b, respectively. (c) Δa0,mix = a0 – a0,ideal for mixed membranes. (d–f) Composition-dependence of the area per lipid molecule for DPPC (dipalmitoylphosphatidylcholine)-DMPC (dimyristoylphosphatidylcholine), DPPC-DOPC (dipalmitoylphosphatidylcholine), and DMPC-DOPC binary membranes. The dashed lines indicate the ideal mixing behavior (see eq ).

We found that the area per lipid molecule, a0, generally decreases with increasing the mole fraction of DPPC (Figure (a)). The trend of the increase in a0 is a0(DOPC-dominated bilayer) > a0(DMPC-dominated bilayer) > a0(DPPC-dominated bilayer). The reverse is true for the bilayer thickness, i.e, l0(DPPC-dominated bilayer) > l0(DOPC-dominated bilayer) > l0(DMPC-dominated bilayer). We may conclude that in mixed PC membranes, two factors determine the degree of packing of hydrocarbon tails and, hence, the surface area per lipid molecules: the hydrocarbon tail length and its degree of saturation. Membranes consisting of longer hydrocarbon chains (like DPPC) occupy less surface area per lipid molecule than those consisting of shorter hydrocarbon tails (like DMPC). The presence of unsaturated bonds in the hydrocarbon tail acts as a defect for chain packing. Therefore, bilayers containing bends due to unsaturated bonds in their hydrocarbon tails (like DOPC) are less ordered and occupy a larger area per lipid molecule than those with saturated hydrocarbon tails of the same length (like DPPC). There is the same trend for the surface area per lipid molecule in the mixed binary and ternary bilayers. Mixed bilayers composed of DPPC (with its longer saturated hydrocarbon tail) as the major component occupy less surface area per lipid molecule than DMPC-rich (with shorter saturated hydrocarbon tail) bilayers. The maximum surface area per lipid molecule, however, belongs to mixed membranes with DOPC (with a longer unsaturated hydrocarbon tail) as the major component. As the hydrocarbon chain lengths in DOPC and DMPC do not differ considerably, the existence of an unsaturated bond in DOPC plays the dominant role in surface area per lipid molecule and bilayer thickness. We have checked if the parameters a0 and l0 for mixed membranes can be expressed as the sum of contributions due to their constituents (ideal mixing), i.e.,andin Figure . In eqs and 11x is the mole fraction of the component i, and subscript “ideal” stands for the ideal mixing. We have quantified the excess area Δa0,mix = a0 – a0,ideal in Figure (c); the same is done for the excess thickness Δl0,mix = l0 – l0,ideal, for which the results are reported in Figure S2 of the Supporting Information. The close agreement between calculated and predicted values indicates that in mixed membranes the ideal mixing rule can be regarded as a fairly good approximation for estimation of the area per lipid molecule and the membrane thickness. Shown in Figure (d),(e) are the area per lipid molecule of the binary membranes as a function of composition and their deviations from ideality. This near-ideal behavior of mixed membranes is in complete agreement with experimental measurement on DMPC-DSP (distearoylphosphatidylcholine)[54] and atomistic simulation results on DPPC-DLPC (dilauroylphosphatidylcholine)[55] membranes.

Another structural parameter to characterize the order of the lipid chain in membranes is the orientational order parameter, S0. S0 = 1 would indicate perfect alignment of the hydrocarbon tails with the chain normal, whereas S0 = 0 corresponds to random orientations. The correlation between S0 and xDPPC for ternary lipid bilayers (in which xDMPC = xDOPC) is shown as an example in Figure (a). We observe that S0 increases with increasing the mole fraction of DPPC. Snapshots of three ternary lipid bilayers with xDPPC = 0.1, 0.5, and 0.9 are shown in Figure (b)–(d). At low DPPC mole fractions (Figure (b)), the DPPC molecules are disordered by the DMPC and DOPC molecules, resulting in a low value of S0. On the other hand, the better ordered DPPC molecules induce ordering in the DOPC and DMPC molecules at higher DPPC mole fractions (Figure (d)). This induced ordering effect is more evident for DMPC than for DOPC molecules. Figure (e) shows S0 as a function of composition for mixed membranes of various compositions. We found that S0 generally increases with the DPPC content. Additionally, the DMPC-rich membranes are better ordered than DOPC-rich membranes. This can be interpreted in terms of a higher order of hydrocarbon chains in DMPC than those in DOPC (see Table ).

Figure 5

(a) Dependence of the orientational order parameter, S0, on the composition for DOPC (dipalmitoylphosphatidylcholine)-DMPC (dimyristoylphosphatidylcholine)-DPPC (dioleoylphosphatidylcholine) ternary bilayers with xDOPC = xDMPC. (b–d) From top to bottom, snapshots of the simulation box indicating mixed bilayers with (xDPPC, xDOPC, xDMPC) = (0.1, 0.45, 0.45), (0.5, 0.25, 0.25), and (0.9, 0.05, 0.05), respectively. Blue, green, and red colors represent the hydrocarbon tails of DPPC, DMPC, and DOPC lipids, respectively. (e) Dependence of the orientational order of mixed bilayers on the composition for bilayers immersed in pure water at 298 K.

We have also examined the dependence of the bending modulus κ0 on the composition of the lipid bilayer in Figure . Our findings indicate that the bending modulus is more sensitive to the composition than the order parameter, the area per lipid molecule, and the membrane thickness. In this case, the DPPC-rich (xDPPC > 0.7) membranes have much higher bending moduli than the others. At the same time, increasing the DPPC mole fraction in binary and/or ternary membranes to a regime with xDPPC ∼ 0.7 has only a marginal effect on the bending modulus, and a sharp increase in the bending modulus is seen at xDPPC > 0.7. The bending modulus of lipid membranes with a higher fraction of DPPC (xDPPC > 0.8) is much higher than those with lower DPPC content. Our findings are indicative of the diverse change in the structural properties of lipid bilayers with varying compositions. Examination of the composition-dependence of the bending modulus and its sharp transition at high DPPC mole fractions may raise a question; would it be possible to design a mixed membrane with tunable stability. This is important for the investigation of the effect of ethanol on the stability of membrane, to be discussed.

Figure 6

Dependence of the bending modulus of mixed membranes on the composition for membranes immersed in pure water at 298 K.

Validation of the Model for Ethanol-Containing Systems

We have tuned the DPD repulsion parameter for the ethanol-tail interaction against our recent atomistic simulation results[56] for the partitioning of ethanol between the aqueous and membrane phases for DPPC immersed in aqueous solutions containing ethanol. The experimental gel to liquid crystalline phase transition temperature of DPPC is 315 K.[57] To characterize the thermodynamic state of our DPPC model, we have shown the temperature-dependence of the order parameter in Figure S2. The sudden change in the order parameter at T = 320 K is indicative of the gel to liquid crystalline phase transition. We performed DPD simulations for DPPC immersed in aqueous solutions containing 5 and 10 mol % ethanol at 323 K. We have calculated the Gibbs free energy (molar) profiles for translocation of ethanol molecules across the membrane from the equilibrium density profiles, i.e.,where μ is the excess chemical potential, i.e., the difference between the chemical potential and that of an ideal gas, ρ is the number density, kB is the Boltzmann constant, T is the temperature, and μ(aqueous) and ρ(aqueous) are the chemical potential and the density in the aqueous phase (surrounding the membrane), respectively. In eq , the ΔG(z) is the molar Gibbs free energy change for transferring a solute molecule i (ethanol) from the bulk aqueous phase to a position z. The free energy profiles for the permeation of ethanol (calculated by scanning the DPD interaction parameter for ethanol-tail, α, between 104 and 110) through a DPPC bilayer immersed in an aqueous solution containing 5 and 10 mol % ethanol are shown in Figure (a). The repulsion parameter α = 105 best predicts the calculated barrier heights from atomistic simulations. While our previous atomistic simulations show that ethanol introduces a big hole in an DPPC membrane immersed in an aqueous solution containing (17.5 mol %) ethanol, the present DPD simulations predict the DPPC membrane to remain intact at similar ethanol concentrations of 17.7 mol % (up to 1 μs). To search for the minimum ethanol concentration in the aqueous phase necessary to cause membrane rupture, we have simulated a number of systems in which DPPC is immersed in aqueous solutions containing 17 to 22 mol % ethanol. Our DPD simulations confirm that the membrane undergoes rupture at a slightly higher ethanol concentration (20.2 instead of 17.5 mol %) (see Figure ). This is a second check for the agreement of our DPD simulation results with atomistic simulation results. Furthermore, in agreement with our atomistic simulation results, we found that the DPPC membrane remains intact at an ethanol concentration of 15.4 mol % at 298 K. This is due to the fact that the DPPC exists in the more intact gel phase at 298 K. The tight packing of lipid molecules in DPPC is weakened in the presence of 15.4 mol % ethanol (Figure (b)), but the membrane remains intact. We show the number density profiles for ethanol and water, partitioned between the aqueous and membrane phases, in Figure (c),(d). We observed that increasing the ethanol concentration in the aqueous phase causes accumulation of both ethanol and water molecules near the membrane head groups, and they penetrate further into the region of the hydrocarbon-headgroup interface. This is in complete agreement with the results of our previous atomistic simulations. Because of the fact that the DPD beads are much softer than the atomistic sites, we observe that water and ethanol molecules are exchanged several times between the aqueous phase and the membrane during the time scale of our simulations.

Figure 7

(a) Free energy barrier for the permeation of ethanol in the DPPC (dipalmitoylphosphatidylcholine) lipid bilayer, immersed in an aqueous solution containing ∼5 and 10 mol % ethanol, at 323 K. The position z = 0 corresponds to the center of the bilayer. The DPD repulsion parameter α for ethanol-tail (hydrocarbon) beads is varied from 104 to 110 to find the best match between barrier heights calculated from atomistic[15] and DPD simulations. (b) and (c) Snapshots (at 1 μs) of a DPPC lipid bilayer, surrounded in aqueous solutions containing (b) 17.7 mol % and (c) 20.2 mol % ethanol, at 323 K. The blue and red spheres show head groups and tails, respectively. The water and ethanol molecules are not displayed for clarity.

Figure 8

Snapshots (at 1 μs) of DPPC (dipalmitoylphosphatidylcholine) lipid bilayers, surrounded by (a) pure water and by (b) an aqueous solution containing 15.4 mol % ethanol at 298 K. The blue and red spheres show head groups and tails, respectively. The water and ethanol molecules are not displayed for clarity. Comparison of the number density profiles, calculated from DPD (full curves) and atomistic (dashed curves) simulations, for (c) ethanol and (d) water across the DPPC membrane surrounded by aqueous solutions containing 0, 5, 10, and 15 mol % ethanol at 298 K.

Effect of Ethanol on the Structure of the Lipid Bilayer

We have summarized the variation of the area per lipid molecule, a, and the membrane thickness, l, for mixed membranes immersed in aqueous solutions containing 5 mol % ethanol and compared them with the corresponding values for ethanol-free water in Figure . We found ethanol uptake expands all pure and mixed lipid bilayers laterally (Figure (a)) and shrinks them vertically (Figure (b)). Among all membranes, the area per lipid and the membrane thickness of DPPC-rich (xDPPC > 0.8) membranes were not significantly affected by the presence of ethanol. The largest changes of the area per lipid molecule and the membrane thickness between membranes immersed in water–ethanol solutions and in ethanol-free solutions are seen for lipid bilayers with 0 < xDPPC < 0.8. This indicates that the effect of ethanol on the stability of a membrane depends on its composition. It is worth mentioning that all pure and mixed membranes examined in this work, surrounded by an aqueous solution containing 5 mol % ethanol, are intact (based on a visual inspection of the simulation and density profiles).

Figure 9

(a) Dependence of the ratio of surface area per lipid molecule for membranes surrounded by an aqueous solution containing 5 mol % ethanol to the corresponding value in the absence of ethanol (a/a0) on the composition of membrane. (b) The same as (a) for the membrane thickness. a0 and l0 are the area per lipid and membrane thickness, respectively, in the ethanol-free systems.

Membrane Failure

Characterization in Terms of Area Per Lipid, Thickness, and Orientational Order

We have searched for the rupture-critical ethanol concentration in the aqueous phase, needed to cause the rupture of pure DPPC, DMPC, and DOPC membranes, immersed in such a solution. In addition to the ethanol concentrations tabulated in Table , we have simulated pure membranes, each immersed in 12 aqueous solutions containing ethanol of varying concentrations from 20 to 32 mol %. Here, the membrane failure was characterized in terms of visual observation of long-lived holes in the membrane. The rupture-critical ethanol concentrations for membrane failure at 298 K are 20.7, 27.5, and 31.7 mol % in the aqueous phase surrounding pure DMPC, DOPC, and DPPC, respectively. In order to establish a link between membrane stability and the parameters (area per lipid molecule a, membrane thickness l, and orientational order of hydrocarbon chains s) discussed above, we have calculated the relative change of each parameter upon transferring the pure membrane (DPPC, DMPC, and DOPC) surrounded by water to an aqueous solution containing ethanol of a given concentration. In Figure , we show the ratios of the area per lipid molecule, membrane thickness, and orientational order of the hydrocarbon tail for pure DPPC, DMPC, and DOPC membranes (immersed in aqueous solutions containing ethanol) to the corresponding value in the absence of ethanol. Because of the ambiguity in calculating the area per lipid molecule and membrane thickness for disrupted membranes, we have reported the ratios a/a0 and l/l0 only for intact membranes (up to 1 μs simulation time) in Figure (a),(b), but the ratio S/S0 (Figure (c)) was reported for both intact and ruptured membranes. For all three membranes, the a/a0 increases (by a factor up to ≈1.6 for DPPC), and l/l0 decreases (by a factor ≈0.6 for DPPC) with the increasing ethanol concentration in the aqueous phase surrounding the membrane. We have also shown the ratio of (a·l)/(a0·l0) as a function of the ethanol concentration (see Supporting Information Figure S4). We found that the membrane volume is less sensitive to the presence of ethanol in the aqueous phase as compared to the surface area and the membrane thickness. Larger fluctuations in the area per lipid molecule and membrane thickness at high ethanol concentrations (which can be regarded as a signature of the stability-failure transition) are the results of frequent opening and closing of holes in the membrane, induced by ethanol. This indicates that ethanol weakens the membrane by fluidizing it, which increases the area per lipid molecule and decreases the membrane thickness, and eventually introduces holes in the membrane, leading to its rupture. The fluidizing effect of ethanol on the membrane obviously decreases the orientational order of the hydrocarbon tail (see Figure (d)). Compared to the ratios a/a0 and l/l0, larger deviations are seen for the ratio S/S0 upon transferring the membrane from pure water to water–ethanol solutions. In this case, the ruptured states can be identified by small values of S/S0 (close to zero), where S = 0 corresponds to the complete random orientation of tails. However, no sharp transition from the intact to the ruptured state, and vice versa, is observed in terms of the orientational order parameter. In other words, we cannot quantitatively predict the location of the phase transition point based on the area per lipid molecule, the membrane thickness, and the orientational order parameter. In order to discriminate between intact and ruptured membranes, we show snapshots (taken at 1 μs) of the structure of DOPC, immersed in aqueous solutions whose ethanol concentration is close to (but below) the rupture-critical ethanol concentration needed to disrupt the membrane (Figure (d)). At the highest ethanol concentrations, the membrane undulates, and even a few lipid molecules are extracted from it. At 26.7 mol % ethanol concentration (just below the phase transition point), the membrane thickness is nonuniform, but it still remains intact (up to 1 μs). Upon further increase of the ethanol concentration, some parts of the membrane become thin enough to allow free (barrierless) passage of ethanol across the membrane, i.e, a big hole is formed in the membrane.

Figure 10

Dependence of (a) area per lipid molecule, (b) membrane thickness, and (c) the orientational order parameter of pure DMPC (dimyristoylphosphatidylcholine), DOPC (dioleoylphosphatidylcholine), and DPPC (dipalmitoylphosphatidylcholine) membranes on the ethanol concentration in the aqueous phase surrounding the membrane. Open (panel c) and filled markers represent the ruptured and intact bilayers, respectively. (d) Snapshots (at 1 μs) of DOPC lipid bilayers at 298 K, surrounded by aqueous solutions containing 24.4–26.7 mol % ethanol, i.e., close to the rupture-critical ethanol concentration (27.5 mol %) needed for membrane rupture. The blue and red spheres show head groups and the lipid tails, respectively.

Ethanol Penetration: Partition Coefficient and Permeation Dynamics

The partitioning of ethanol between the lipid membrane and the aqueous phase can be characterized in terms of the partition coefficientwhere cethanollipid and cethanolaq are the equilibrium mole concentrations of ethanol in the lipid membrane and in the aqueous phase, respectively. It is noted that cethanolaq is defined as the equilibrium ethanol concentration, which is different from the initial ethanol concentration (before equilibrium). In order to calculate K, the number of ethanol molecules in the lipid membrane (nethanollipid) and aqueous phase (nethanolaq) and the respective volumes of these two regions (Vlipid and Vaq) are computed. These parameters, however, cannot be calculated unambiguously, especially at the high ethanol concentrations, where the interface between the regions becomes blurred. Terama et al.[58] have calculated cethanolaq by assigning the ethanol molecules in the region away from the membrane boundary to the aqueous phase (where the density profiles almost converge to a constant value) and assigning the rest of the ethanol molecules (bounded) to the membrane. These ”bounded” ethanol molecules are partly influenced by the membrane, although not completely partitioned into it. Inspired by a recent work,[59] we have adopted a method based on the assumption that any deviations in the Gibbs free energy profiles from the corresponding values in the aqueous phase (ΔGbulk ≈ 0) must be due to interactions with the membrane. We have partitioned the simulation along the z direction into a number of bins of width 0.305 nm. The boundary between the lipid membrane and aqueous phases is defined as the layer with |ΔG(z) – ΔG(z)bulk| = 0.1 kJ/mol. The results for K of three single-lipid membranes in the aqueous phases containing ethanol of varying concentrations from 5 mol % to the concentrations just before rupture are shown in Figure . The partition coefficients of the pure membranes increase with the increasing ethanol concentration in the aqueous phase surrounding the membrane. This means that the amount of ethanol in the membrane increases disproportionately with the ethanol concentration in the aqueous phase. This is indicative of a synergistic effect. Moreover, the values of the partition coefficients close to the rupture points are nearly the same for the three membranes.

Figure 11

Dependence of the partition coefficient K of pure DMPC (dimyristoylphosphatidylcholine), DOPC (dioleoylphosphatidylcholine), and DPPC (dipalmitoylphosphatidylcholine) membranes on the ethanol concentration, in the aqueous phase surrounding the membrane, at 298 K.

To study the dynamics of the permeation of ethanol through the membrane, we applied an external force to drag two selected ethanol molecules through the membrane (z direction).[60] We used pure DPPC, DMPC, and DOPC membranes immersed in pure water (where the membrane is stable) and in aqueous solutions containing 20.7, 27.5, and 31.7 mol % ethanol (corresponding to the state point at which the membrane disrupts), respectively. For membranes in pure water, two ethanol molecules were added to the aqueous phase for this purpose. In ethanol-containing systems, two ethanol molecules were selected near the box boundary. Simulations were performed from the equilibrated systems (at 1 μs). Constant external forces of equal magnitude but opposite directions were imposed on the two ethanol molecules (Figure (a)). We tried several drag forces with magnitude ranging from 0.58 to 17.4 pN for all systems. No crossing event has been observed for the DPPC membrane in pure water with an external force up to 5.8 pN within 143 ns. For DOPC and DMPC lipid membranes, however, we have seen at least one crossing event with the same force. For the DPPC membrane, increasing the external force, a full cycle (an ethanol molecule traversing the entire box length along the z direction) was observed at 8.8 pN within 143 ns. Hence, we have fixed the external force at 11.6 pN, to be able to observe enough crossing events and also not so fast crossing the membrane in the high ethanol concentration systems. As shown in Figure (b), for membranes immersed in pure water, the crossing period for DPPC is longer than those of DPMC and DOPC, and the limitation of the crossing is due to the time spent in the membrane. At high ethanol concentrations, however, the crossing periods for all three membranes are comparable, and crossing the membrane is as fast as crossing the aqueous phase. In order to quantitatively describe the crossing behavior, we have run 5 independent simulations for each system for a long time (715 ns). The average velocity of ethanol molecules crossing the lipid membrane is defined based on the average time it takes them to cross the membrane. The results have been summarized in Table . For all three membranes, the number of crossing events and the crossing rate are similar when the membranes are in the ruptured state.

Figure 12

(a) Snapshots (at 1 μs) of the DPPC (dipalmitoylphosphatidylcholine) lipid bilayer, surrounded by pure water with the addition of two ethanol molecules (purple colored spheres), experiencing external forces along the z axis in different directions at 298 K. The blue and red spheres show head groups and lipid tails, respectively, and the green spheres show water molecules. Tracking the position of ethanol molecules experiencing external forces along the z axis (11.6 pN) in pure DMPC (dimyristoylphosphatidylcholine), DOPC (dioleoylphosphatidylcholine), and DPPC membranes surrounded by (b) pure water and (c) aqueous solutions containing 20.7, 27.5, and 31.7 mol % ethanol, respectively, at 298 K.

Table 4

Number of Crossing Events within 715 ns and the Average Velocity of Ethanol Molecules Crossing the Lipid Membrane Experiencing an External Force of 11.6 pN at 298 Ka

system	ethanol (mol %)	ncross	|vcross| (m/s)
DPPC in intact state	0	6	0.06
DMPC in intact state	0	14	0.10
DOPC in intact state	0	10	0.12
DPPC in ruptured state	31.7	17	0.41
DMPC in ruptured state	20.7	19	0.44
DOPC in ruptured state	27.5	18	0.43

DMPC (dimyristoylphosphatidylcholine), DOPC (dioleoylphosphatidylcholine), and DPPC (dipalmitoylphosphatidylcholine) membranes surrounded by aqueous solutions containing 20.7, 27.5, and 31.7 mol % are in the ruptured state.

Characterization in Terms of Machine-Learned State Variables

We have further implemented a machine-learning framework to characterize the failure of lipid membranes, namely, MembraneNN. The key component of this framework is a deep neural network (DNN), which is a nonlinear model mapping the particle coordinates of a single lipid molecule to state variables representing the order of this lipid molecule in the membrane. As depicted in Figure , the DNN is a standard feed-forward network composed of fully connected layers, in which the data flows from the input layer through the hidden layers toward the output layer. Each layer has a certain number of nodes (below), called neurons, which store information about the importance of the input and associations between the importance of combinations of inputs. In all layers, linear operations are applied to the input data din, e.g., , where w and b are the weights from layer l to the output layer k and the bias produced at layer k, respectively (Figure (b)). Additionally, the rectified linear unit (ReLU) activation functionwhere x is the input to a neuron, is applied in intermediate layers to break the linearity: .[61] In the output layer, the neurons are activated by a sigmoidal activation function to produce two state variables in the range of 0 to 1. Two state variables λ1 and λ2 are introduced to characterize the membrane failure. We only label the data at two extreme states: λ1 = 1 and λ2 = 0 (representing an intact membrane) and λ1 = 0 and λ2 = 1 (representing a ruptured membrane); see Figure (a). Overfitting is usually a serious problem in DNN models with an increasing number of layers and neurons per layer. To prevent overfitting and to improve the generalization error,[62] the dropout technique,[63] which randomly drops out neurons (and their connections) during training, is implemented in our DNN model with a dropout rate of 0.2. This dropout rate indicates that there are 20% neurons randomly deactivated in each hidden layer. Another important element of the DNN is the loss function. Essentially, the loss function defines properties that the DNN attempts to optimize. Here, we use the mean square error (MSE)where n is the number of lipid molecules for training, and y and ỹ are the manually labeled target state variables and state variables predicted by the DNN (λ1 and λ2) for a single lipid molecule, respectively.

Figure 13

(a) Schematic figure of the deep neural network model. The Cartesian coordinates of all particles in a single lipid molecule are fed as input data to the network. The network has an output layer with two values, λ1 and λ2, that represent the intact and ruptured states of a single lipid molecule, respectively. (b) Schematic figure for a single neuron.

Figure 14

Flowchart of the MembraneNN framework to characterize the failure of lipid membranes (MembraneNN) in this work. (a) Example of traning data: snapshots of the intact and ruptured DMPC (dimyristoylphosphatidylcholine) lipid membranes (h1 head groups with blue color, h2 and h3 head groups with red color, and lipid tails with cyan color). (b) Architecture of the DNN model. (c) Snapshots of DMPC lipid membranes in pure water and in aqueous solutions containing 23.7 and 27.5 mol % ethanol (from top to bottom). The blue and red colors in the snapshots show head and tail lipid groups, respectively. (d) Variation of state variables λ1 and λ2 for a DMPC bilayer, immersed in aqueous solutions containing various concentrations of ethanol at 298 K. Open and filled markers represent the ruptured and intact bilayers, respectively, as determined by the visual inspection.

The DNN model is built and trained using PyTorch[64] (version 1.8.1). Hyperparameters, whose values are used to control the learning process such as the number of layers and neurons in each layer, are usually difficult to determine.[65] In our framework, we use a grid search to determine these parameters in order to achieve efficiency and accuracy of our DNN. In the production run, we choose the network architecture: 1 input layer, 3 hidden layers, and 1 output layer with neurons in each layer 30-90-90-90-2. The ADAM (adaptive moment estimation) optimizer[66] with a learning rate of 0.001 is employed to train the DNN model. During training, data are iteratively fed into the network until a termination criterion is reached. Here, we stop the training when a reasonably low value of the loss function is reached, and it does not decrease significantly in further runs.

The workflow of MembraneNN is summarized here (Figure ): The first step is to prepare the data for training and prediction. It is basically the trajectories of lipid molecules generated from the coarse-grained simulations of lipid membranes with various ethanol concentrations. At each ethanol concentration, particle coordinates of lipid molecules from 40 trajectory frames at equilibrium are collected, which means we have totally nmole ≈ 40 × 1000 ≈ 4 × 104 samples for developing the DNN model for each type of lipid molecule. It is noted that we move the center-of-mass position of the lipid molecule to the Cartesian origin in order to make the coordinates of a single lipid molecule translation-invariant for training. As aforementioned, we only label the lipid molecules (DPPC, DOPC, and DMPC) at the intact and ruptured states to be (λ1 = 1, λ2 = 0) and (λ1 = 0, λ2 = 1), respectively. These labeled data are used separately to develop the DNN models with 70% for training and 30% for testing, resulting in one DNN model per lipid type. We use the trained models to predict state variables (λ1 and λ2) for lipid molecules at various ethanol concentrations (an example is shown in Figure (d)). The integrity of the entire membrane is then defined as the averaged state variables of all lipid molecules in the membrane at a given alcohol concentration. We observe noticeable distributions of the predicted state variables λ1 and λ2 of lipid molecules in an intact membrane (see examples in Figure S6), although we manually labeled all lipid molecules of an intact membrane to be λ1 = 1 and λ2 = 0.

We employ the MembraneNN framework to predict the integrity of all pure and mixed membranes. The intact-to-ruptured transition point is determined via the procedure described in the section Machined Learning Model of the Supporting Information. The results are summarized in Table . The rupture-critical ethanol concentration in the aqueous phase surrounding the mixed membranes, composed of DMPC, DOPC, and DPPC, varies between ≈20 mol % ethanol (for DMPC dominated membranes) and ≈31 mol % ethanol (for DPPC dominated membranes). It is lowest for pure DMPC (the least robust membrane) and highest for pure DPPC (the most resilient membrane).

Table 5

Rupture-Critical Ethanol Concentration in the Aqueous Phase, Surrounding Lipid Membranesa

xDPPC	xDMPC	ethanol (mol %)	xDPPC	xDMPC	ethanol (mol %)	xDPPC	xDMPC	ethanol (mol %)
0.0	1.0	20.7	0.3	0.4	25.6	0.1	0.1	28.2
0.1	0.9	20.8	0.0	0.2	25.9	0.2	0.0	28.3
0.0	0.9	20.9	0.1	0.45	26.0	0.45	0.1	28.8
0.0	0.8	21.5	0.0	0.1	26.1	0.3	0.0	28.9
0.05	0.9	22.1	0.45	0.45	26.2	0.4	0.0	28.9
0.2	0.8	22.5	0.4	0.2	26.8	0.6	0.2	29.2
0.1	0.8	22.5	0.3	0.3	27.0	0.5	0.0	29.4
0.0	0.6	23.0	0.4	0.4	27.1	0.8	0.0	29.6
0.0	0.4	23.4	0.333	0.333	27.4	0.8	0.2	29.8
0.0	0.5	23.6	0.0	0.0	27.5	0.9	0.1	30.2
0.4	0.6	24.1	0.25	0.25	27.5	0.8	0.1	30.4
0.5	0.5	24.4	0.4	0.3	27.6	0.9	0.0	31.0
0.2	0.6	24.6	0.1	0.0	28.0	0.9	0.05	31.1
0.2	0.4	24.8	0.5	0.25	28.0	1.0	0.0	31.6
0.6	0.4	25.2	0.2	0.2	28.0
0.25	0.5	25.3	0.05	0.05	28.2

xDPPC and xDMPC indicate the mole fractions of DPPC (dipalmitoylphosphatidylcholine) and DMPC (dimyristoylphosphatidylcholine) in the membrane, and the mole fraction of DOPC (dioleoylphosphatidylcholine) can be derived as xDOPC = 1 – xDPPC – xDMPC.

We have checked whether there exists a correlation between the characteristics of mixed membranes immersed in pure water (always intact) and their stabilities against ethanol. In Figure , we correlate the orientational order parameter for the hydrocarbon tails and the bending modulus of mixed membranes, immersed in water, with the rupture-critical ethanol concentration in the aqueous phase, at which the membrane disrupts. A general trend for stabilities in the presence of ethanol is seen; a higher ethanol concentration in the aqueous phase is needed to disrupt a more ordered (higher S0) and stiffer (higher κ0) mixed membrane. The DPPC-rich membranes (xDPPC > 0.8), however, do not follow the trend of stability of the other mixed membranes; the values of both S0 and κ0 are much larger than the corresponding values for membranes with xDPPC < 0.8. On the other hand, an unusually high ethanol concentration in the aqueous phase is needed to disrupt the DPPC-rich membranes. The reason is that at T = 298 K, where we did our simulations, the DPPC exists in the gel phase, but the DMPC and DOPC are always in the fluid phase. This is in agreement with the reported gel–liquid crystalline phase transition temperature of DPPC (315 K),[57] DMPC (297.25 K),[67] and DOPC (256.65 K).[57] Therefore, the hydrocarbon tails of the DPPC are much better ordered than those of DOPC and DMPC, and hence, the DPPC has a much higher bending modulus than the DMPC and DPOC. However, when the surrounding solution contains ethanol, it penetrates into the membrane and fluidizes it. Upon increasing the ethanol concentration in the aqueous phase to ≈30 mol %, even the DPPC-rich membranes become weak enough (as a result of ethanol penetration) to rupture.

Figure 15

Correlation between (a) the orientational order parameter and (b) the bending modulus of membranes, immersed in pure water, and the rupture-critical ethanol concentration (in the aqueous phase surrounding the membrane) at which the membrane disrupts.

Conclusions

We have developed DPD models of three lipids, DPPC (dipalmitoylphosphatidylcholine), DOPC (dioleoylphosphatidylcholine), and DMPC (dimyristoylphosphatidylcholine), and their neat, binary, and ternary mixed membranes as models of the SARS-CoV-2 membrane, immersed in water and in water–ethanol solutions (disinfectants) up to concentrations where the membranes undergo rupture. It is worth noting that the composition of the SARS-CoV-2 membrane is not known; however, it is presumably composed of different phospholipids stolen from its host, which anchor several membrane proteins. Since our previous atomistic simulations[31] revealed that the E-peptide does not offer noticeable protection against the ethanol-induced failure of the viral membrane, we have simulated mixed membranes composed of DPPC, DMPC, and DOPC, as models of the SARS-CoV-2 membrane without membrane proteins. The molecular mechanisms of ethanol-induced weakening and the impact of ethanol concentration on the deactivation (membrane rupture) of SARS-CoV-2 have been investigated in detail. We have found that two factors influence the membrane stability: the lipid composition of the membrane and the concentration of ethanol in the aqueous disinfectant solution surrounding it. We looked at different properties, including the area per lipid molecule, bilayer thickness, orientational order of the lipid tails, bending modulus, and ethanol permeability, and their dependence on membrane composition and ethanol concentration.

Our DPD models for pure membranes were validated by comparing their areas per lipid molecule and thicknesses against experimental data and reported atomistic simulation results in the literature. Based on the agreement of the calculated properties for neat bilayers in pure water with experimental and previous atomistic simulation reports, we conjecture that at 298 K (where we did our DPD simulations) the DPPC exists in the gel phase, but the DMPC and DOPC exist in the fluid phase. For mixed DPPC-PMPC-DOPC bilayers in pure water, we found that structural properties depend on the length and the degree of saturation of the hydrocarbon tail of the components. Mixed bilayers composed of DPPC (with longer saturated hydrocarbon tails) as the major constituent occupy less surface area per lipid molecule than DMPC-rich (with shorter saturated hydrocarbon tails) bilayers. However, the maximum surface area per lipid molecule belongs to mixed bilayers dominated by DOPC (with unsaturated bonds in the hydrocarbon tail). Expectedly, the bilayer thickness depends inversely on the surface area per lipid molecule, i.e., the DPPC-rich mixed membranes are thicker than the DOPC- or DMPC-rich bilayers. In the absence of ethanol, the variation of the surface area per lipid a0 and the bilayer thickness l0 is about 20% between more fluidic (DMPC and DOPC) and gel-like (DPPC) membranes. Therefore, a0 and l0 cannot be usefully employed to predict the stabilities of mixed bilayers of different compositions. The orientational order of the hydrocarbon tail S0 is, however, a more sensitive parameter in this respect; it varies by a factor of 2 from DMPC/DOPC to DPPC. Interestingly, an even sharper composition-dependent change is seen for the bending modulus κ0. It varies by a factor of ≈26 from DMPC/DOPC to DPPC. Another noticeable point is that while the values of a0 and l0 for mixed membranes nearly follow an ideal mixing rule, the values of S0 and κ0 do not. In terms of the bending modulus, a marked region of stability is seen for DPPC-rich (xDPPC > 0.8) membranes.

In complete agreement with our previous atomistic simulations,[15] we found that aqueous solutions containing 5–10 mol % ethanol have a significant weakening effect on the neat and mixed membranes. With increasing ethanol concentration in the disinfectant solution, the ethanol uptake of the membrane increases disproportionately: There is a synergistic effect, as the ethanol already in the membrane fluidizes it and facilitates the absorption of even more ethanol. The dissolution of ethanol in the membrane causes lateral membrane swelling and the shrinkage of its thickness. Obviously, the ethanol uptake reduces the orientational order of the hydrocarbon tails of the lipids. However, we cannot quantitatively predict the location of the phase transition point based on the area per lipid molecule, the membrane thickness, or the orientational order parameter. Hence, we have further developed a machine-learning framework to access the integrity of lipid membranes in place of visual inspections. We found that the rupture-critical ethanol concentrations needed to induce damage in the mixed membranes vary between that for pure DMPC (the least robust membrane) to that for pure DPPC (the most robust membrane). At 30 mol % ethanol concentration, only DPPC-rich membranes (xDPPC > 0.8) remain intact. This stabilizing effect of DPPC can be correlated with a larger bending modulus for membranes in pure water. The tight packing of DPPC molecules, which forms a gel phase when pure, makes the bilayer less permeable to ethanol. It also induces orientational order among the DOPC and DMPC lipids at high DPPC concentrations. Membrane failure can not only be identified by human or machine inspection of its structure. Also the relative transmembrane permeation rates of ethanol, calculated by the nonequilibrium molecular dynamics, show a sudden increase when the membrane develops holes, jumping by as much as a factor of ≈7 for DPPC and a smaller factor of ≈4 for DMPC and DOPC.

The least robust membrane studied in this work is pure DMPC, which is punctured already by a disinfectant solution containing ≈20 mol % (corresponding to ≈40 wt %) ethanol. On the other hand, as the DPPC exists in the gel phase at room temperature, it forms one of the most resilient membranes of the coronavirus. Even such a robust membrane, however, fails in aqueous disinfectant solutions containing more than ≈32 mol % (corresponding to ≈55 wt %) ethanol. Therefore, an ethanol concentration below ≈20 mol % in the disinfectants is hardly efficient in the deactivation of the coronavirus. On the other hand, an ethanol concentration above ≈32 mol % in the disinfectants will make the viral membrane dysfunctional, regardless of which of the lipids studied here dominates in the membrane.