Impact of Lipid Composition on Membrane Partitioning and Permeability of Gas Molecules

1. Introduction

Lipid membranes provide physical barriers that control the exchange and flow of substances in and out of living cells [1,2]. Unlike polar and charged molecules such as ions, which are energetically unfavorable to dissolve in the hydrophobic core of a lipid membrane [3,4], small, nonpolar gas species, such as O₂, NO, and CO₂, may diffuse readily through pure lipid bilayers. Although this mode has long been accepted as the main mechanism for gas transport across the cellular membranes  [5,6,7], some biological membranes exhibit surprisingly low gas permeabilities [8,9]. Examples include the gastric parietal chief cells [10], erythrocytes [11,12], and fiber cells of the eye lenses [13]. These reports evoke the question as to what degree the lipid composition can determine the gas permeability of a membrane. Distinct gas permeation rates have been reported for membranes composed of different lipid constituents [13,14,15,16,17,18,19,20]. Membrane channels, e.g., aquaporins and Rhesus (Rh) glycoproteins, are also documented to facilitate gas transport across the membrane [21,22,23,24,25,26,27,28], but their physiological significance in this regard remains disputed [29,30], mainly due to experimental limitations [31,32]. Demonstrating that the lipid composition can modulate gas permeability is the first step in establishing the physiological significance of membrane proteins in the overall gas exchange in living cells.

Given its microscopic results, molecular dynamics (MD) simulation can offer insights into the delivery and transport of molecules across the membranes at the atomic level. Physical properties, such as molecular partitioning and permeability coefficients can also be predicted using MD [3,33,34,35]. For years, lipid membranes have been simulated with only one glycerophospholipid type, such as POPC, DOPC or POPE, which are highly fluid under ambient conditions. In recent years, all-atom simulations of membranes with heterogeneous lipid compositions containing multiple glycerophospholipid species, as well as cholesterol (CHL) and sphingomyelins (SM) have been performed [36,37,38].

CHL and SM (Figure 1) are key lipids for structural integrity and morphology of mammalian membranes, and perform other regulatory functions [2,39]. They are abundant in typical mammalian plasma membranes, including erythrocyte and lens membranes, which are reported to have low gas permeability [40]. In those membranes, CHL:phospholipid ratios can range from 1:1 to 4:1, where SM constitutes $>15$% of the lipid content [41,42]. CHL is enriched in the liquid-ordered phase of the membrane [43,44] and also known to induce lipid domain formations [45,46,47]. The clustering between SM with CHL has been implicated to the formation of liquid, highly-ordered domains, known as rafts, which have been proposed to be involved in cell signaling [48,49,50,51].

The effects of CHL on membrane permeability have been probed by MD simulations. Previous studies with extensive simulations and analysis of membranes with high CHL content (>30 mole%) predicted a reduced water permeability [52,53]. Other studies also found apparent reduction in the partitioning of small molecules and gases in membranes with CHL relative to CHL-free membranes [54,55,56,57].

The present study examines the effects of lipid composition heterogeneity, particularly the presence of high CHL and SM, on gas permeability by explicitly simulating O₂ or CO₂ molecules in membranes. The sets of independent simulations comprise CHL-free (pure) bilayers of POPC, SM or DPPC, binary compositions of POPC:CHL, SM:CHL, DPPC:CHL, POPC:SM and POPC:DPPC, and ternary compositions of POPC:PSM:CHL and POPC:DPPC:CHL. SM is represented by the palmitoylsphingomyelin (PSM) (Figure 1). The fully saturated DPPC (dipalmitoylphosphatidylcholine) lipids account for the effects of lipid saturation on gas permeability. We explicitly calculated membrane partitioning free energy profiles and permeability coefficients of the gas molecules from the equilibrium MD trajectories. The results show clearly the dependence of the calculated gas permeability on the lipid composition; they demonstrate that high CHL content reduces solubility of gas molecules in the membrane. Reduced gas permeability is pronounced in the simulated POPC-free PSM membranes, due to slow lipid dynamics and gel-liquid phase behavior under physiological body temperature (310 K).

2. Materials and Methods

2.1. Simulation Systems

The simulated membrane systems are summarized in Table 1. They comprise mixtures of POPC:PSM:CHL or POPC:DPPC:CHL lipids. The membrane patches were prepared using CHARMM-GUI [58]. The membranes were solvated with TIP3P water and ionized with 0.2 M NaCl. Structures and interactions of lipids (POPC, PSM and CHL), water, ions and gases were described by the CHARMM36 force field parameter set [59,60,61,62]. The initial dimension of each membrane plane was 100 Å × 100 Å.

Depending on the membrane system, a simulation was performed to equilibrate the system in the absence of gas molecules for tens to a few hundred nanoseconds. The area of the $xy$ plane of each system was monitored and the simulation was continued until it became steady. To explicitly describe the partitioning and permeation of gases, O₂ or CO₂, 125 copies of each species were added to the equilibrated systems. The molecules were initially placed in the aqueous solution. We refer to this type of simulations as “flooding” simulations.

2.2. Simulation Protocols

All simulations were performed using NAMD2 [63,64] with the CHARMM36 force field [59,65] and a time step of 2 fs. All of the simulations were performed under NPT ensemble with the flexible cell under the constant $x/y$ ratio, with the z axis representing the membrane normal. The periodic boundary condition (PBC) was used throughout the simulations. To evaluate long-range electrostatic interactions in PBC without truncation, the particle mesh Ewald (PME) method [66] with a grid density of 1/ų was used. The temperature was maintained at 310 K by Langevin dynamics [67] with a damping coefficient $\gamma$ of 1/ps. The modified Nosé-Hoover method [67,68], in which Langevin dynamics was used to maintain the pressure at 1 atm with a piston period of 200 fs. All bonds involving hydrogen atoms were kept rigid using the SHAKE algorithm [69]. The cutoff for van der Waals interactions was set at 12 Å.

2.3. Membrane Partitioning and Permeability of Gases

The partitioning free energy profiles ($\Delta G$) and membrane permeability coefficients ($p_{\mathrm{d},\mathrm{mem}}$) of O₂ or CO₂ were calculated using the last 175 ns part of trajectories for the 100% POPC membrane systems, and the last 400 ns parts of trajectories for the pure PSM, 0:50:50 POPC:PSM:CHL, and 0:50:50 POPC:DPPC:CHL systems. For the other membrane systems (50:0:50, 35:0:65, 50:50:0, 33:33:33 and 25:25:50 POPC:PSM:CHL, and 50:50:0 and 33:33:33 POPC:DPPC:CHL), which were simulated for 325 ns, we used the last 250 ns for the analysis. $\Delta G$ profiles representing the occupancy of the gas molecule was expressed as:

$$\Delta G=-\mathrm{R}Tln{\displaystyle \frac{P_i}{P_{\mathrm{bulk}}}},$$

with $P_i$ and $P_{\mathrm{bulk}}$ representing occupancy of the gas molecule at position i and in the aqueous solution, respectively. The VolMap plugin of VMD was used to calculate the occupancy coefficients, which were then projected onto a one-dimensional (z) axis (membrane normal) using the Volutil plugin.

$p_{\mathrm{d},\mathrm{mem}}$ was determined from the cumulative profiles of the number of permeation events (N) from either the extracellular side or the intracellular side over a time period ($\delta t$) using linear regression. It is in the unit of cm/s and expressed as:

$$p_{\mathrm{d},\mathrm{mem}}={\displaystyle \frac{m}{N_A{C}_{b}A}},$$

where m is the slope of N over $\delta t$, $C_b$ is the concentration of the gas molecule in the aqueous solution, and A is the area of the $xy$ plane of the membrane. $C_b$ was determined from the average number of the gas molecules outside the membrane normalized by the number of water molecules in the same region.

The gas molecules were considered to be in the extracellular solution when they were at least 3 Å above the center of mass of the choline moiety of the lipid head groups, and within the membrane when they were at least 3 Å below the center of mass of the choline groups.

3. Results

3.1. Variation of Membrane Structure by Lipid Composition

3.1.1. Occupancy Profiles

We first monitored the changes in membrane structure by calculating the atomic profiles of its constituents including water and lipids (Figure 2).

In all the simulated membranes, the midpoint of the lipid bilayer ($z=0$) corresponds to the lowest occupancy region. In the pure POPC bilayer, as well as in the 50:50 POPC:PSM and POPC:DPPC bilayers (Figure 2A, G and J), the highest occupancy regions correspond to the head group regions located at $z=\pm 22$ Å, inferring a membrane thickness of ∼44 Å. Coincidentally, these points are at the intersection of the water and lipid occupancy profiles, also marking the lipid-water interfaces.

Having simulated the POPC membranes with a high content of CHL (50 and 65 mole%) and without it, we found that CHL significantly alters the membrane profiles (Figure 2A vs. B and C). In the presence of high CHL, the maximum density along the membrane normal spans from the head group regions to the lipid tails centered at $z=\pm 10$ Å. CHL molecules are localized beneath the phospholipid head groups. Their positioning decreases the depth of water penetration into the membrane, as indicated by the observed dehydration between $z=-15$ and $z=-10$ Å, and between $z=10$ and $z=15$ Å.

For the pure PSM bilayer, the lipid-water interfaces center at $z=\pm 24$ Å, thickening the membrane by ∼4 Å (Figure 2D). The membrane appears to be more condensed than the CHL-free POPC, POPC:PSM and POPC:DPPC bilayers. Its average area-per-lipid is ∼51 Ų, whereas those of the CHL-free POPC, POPC:PSM and POPC:DPPC bilayers are ∼65, 59, and 62 Ų, respectively. The maximum lipid density of the pure PSM bilayer plateaus to values between $z=-20$ and $-5$ Å, and between $z=5$ and 20 Å, respectively (Figure 2D). For the 50:50 PSM:CHL and DPPC:CHL bilayers, the distance between the two lipid-water interfaces increases to about $z=$50 Å. However, the PSM:CHL one is more condensed than the DPPC:CHL one by ∼15% (Figure 2E-F). Similar to the 50:50 and 35:65 POPC:CHL bilayers, a high CHL content decreases the hydration in both PSM:CHL and DPPC:CHL bilayers. This same behavior is also apparent in the 33:33:33 POPC:PSM:CHL and POPC:DPPC:CHL bilayers (Figure 2H, I and K).

3.1.2. Lipid Ordering

The lipid molecules in the pure POPC bilayer as well as those in the other CHL-free bilayers appear disordered and highly fluid (Figure 3). Those in the simulated CHL-containing bilayers, on the other hand, are conformationally ordered and structurally extended. Uniquely, the pure PSM bilayer resembled more like a ripple phase (between gel and liquid phases). Its shape visibly appears undulated and its lipid molecules appear clustered into small domains (Figure 3).

To quantify the conformational dynamics of the lipid molecules, we calculated the deuterium order parameters ($-S_{\mathrm{CD}}$) for the aliphatic chains of POPC and DPPC, the sphingosine and acyl chains of PSM (Figure 4).

An $-S_{\mathrm{CD}}$ value of 0.5 indicates that the lipid molecules are ordered and that their tails are conformationally extended (not tilted). A value of 0 means that the lipids are highly disordered. Among the simulated membranes, the pure POPC bilayer is the most disordered membrane with the maximum $-S_{\mathrm{CD}}$ values of ∼0.2 for both palmitoyl and oleoyl chains (Figure 4A). Lipids in the simulated CHL-containing bilayers are highly ordered with maximum $-S_{\mathrm{CD}}$ values approaching 0.4 for POPC lipids and above 0.4 for PSM and DPPC lipids. In the pure PSM bilayer, $-S_{\mathrm{CD}}$ values are slightly above ∼0.3, suggesting appreciable tilting of some lipid molecules, forming an undulated membrane. The lipid order is reduced in the 50:50 POPC:PSM bilayer with the drop of the maximum $-S_{\mathrm{CD}}$ values of the PSM chains to 0.275 (Figure 4B). In the 50:50 POPC:DPPC bilayer, the $-S_{\mathrm{CD}}$ values are ∼0.2 (Figure 4C).

To confirm the tilting of the lipid molecules, we also calculated the tilt angles ($\theta$) of the lipid tails (using carbon atoms) with respect to the membrane normal. In the absence of CHL, lipid tails adopt a broad range of orientations with $\theta$ spanning from 0 to 90°; the distribution peaks are ∼25° for POPC and DPPC lipids and ∼20° for PSM lipids (Figure 5).

This reflects the conformational heterogeneity of the lipid molecules, showing some are highly tilted or structurally bent. In the simulated bilayers with CHL, lipid tails become significantly less tilted; the $\theta$ values range between 0 and 90° with the peaks at 10°, consistent with increasing lipid order. For the pure PSM bilayers, although the distribution peaks of $\theta$ are ∼10°, $\theta$ spreads from 0 to 70° (Figure 5B), confirming high degrees of tilting in some lipid molecules.

3.2. Membrane Partitioning and Permeability of Gas Molecules

3.2.1. Pure POPC and Binary POPC:CHL Bilayers

The changes in membrane structure and dynamics by lipid compositions can have profound effects on the membrane partitioning and permeability of gas molecules. Using the equilibrium flooding MD trajectories, we directly calculated the permeability coefficients ($p_{d,mem}$) and partitioning $\Delta G$ profiles of O₂ and CO₂ (Table 2, Figure 6 and Figure 7). For the pure POPC bilayers, the last 175 ns of the 225-ns simulations were used for the calculations. The calculated $p_{d,mem}$ values are ∼12 and 16.6 cm/s for O₂ and CO₂, respectively (Figure 6A-B). The calculated $\Delta G$ profiles indicate the highest accumulation of O₂ and CO₂ to be in the midplane of the membrane, corresponding to free energy minima ($\Delta G_{min}$) of $-2$ kcal/mol for O₂ and $-1$ kcal/mol for CO₂ with respect to the bulk solution. The free energy maxima ($\Delta G_{max}$) of ∼0.5 kcal/mol, indicating the lowest occupancy or solubility regions of the gas molecules, are located in the headgroup regions. Thus, according to the $\Delta G$ profiles, the free-energy barriers of diffusing from the center of the membrane to the aqueous solution ($\Delta G^{\mathrm{cen}-\mathrm{bulk}}$) are 2.5 kcal/mol for O₂ and 1.25 kcal/mol for CO₂ (Figure 7A and C).

For the simulated bilayers with 50 or 65 CHL mole%, the last 250 ns of the 325-ns simulations were used in $p_{d,mem}$ and $\Delta G$ calculations. O₂ and CO₂ remain accumulated in the membrane relative to the aqueous solution, and $\Delta G_{min}$ also remains located in the midplane of the membrane. The incorporation of CHL molecules in the membrane decreases the occupancy of the gas molecules in the tail regions, as indicated by higher $\Delta G_{\mathrm{tail}}$ values, peaked at $z=\pm 10$ Å (Figure 7A). In the pure POPC bilayers, the $\Delta G$ value drastically decreases as the gas molecule diffuses from the head group region towards the center of the membrane. This decrease in membrane solubility, however, does not decrease O₂ permeability. At 50 CHL mole%, $\Delta G_{\mathrm{tail}}$ plateaus at ∼0 kcal/mol, so $\Delta G^{\mathrm{cen}-\mathrm{bulk}}$ for O₂ is ∼0.5 kcal/mol lower than in the pure POPC bilayer, which results in a higher $p_{\mathrm{d},\mathrm{mem}}$ value of 16 cm/s. At 65 mole%, $\Delta G_{\mathrm{tail}}$ increases to 0.5 kcal/mol, forming new local $\Delta G_{min}$ at ±18 Å. Still, $\Delta G^{\mathrm{cen}-\mathrm{bulk}}$ is about the same as that in the pure POPC bilayer (Figure 6A), as reflected by a $p_{\mathrm{d},\mathrm{mem}}$ value of 10 cm/s, which corresponds to a 16 % decrease.

In contrast to O₂, the presence of high CHL contents decreases CO₂ permeability (Figure 6B). The calculated $p_{\mathrm{d},\mathrm{mem}}$ values are ∼10 cm/s at 50 mole%, and ∼6 cm/s at 65 mole%. $\Delta G_{\mathrm{tail}}$ are ∼0.25 kcal/mol at 50 mole%, and ∼1 kcal/mol at 65 mole%, becoming the new global $\Delta G_{max}$ (Figure 7C).

3.2.2. Special Cases of Pure PSM and Binary PSM:CHL Bilayers

For the PSM bilayers with 0 and 50 CHL mole%, the last 400 ns of the 525-ns simulation trajectories were used to calculate the $p_{\mathrm{d},\mathrm{mem}}$ and $\Delta G$ profiles. The solubility of the gas molecules remain highest at the midplane of the membrane. In the pure PSM bilayers, no local $\Delta G_{max}$ is formed in the tail regions despite substantial alterations of membrane structure and lipid dynamics compared to the pure POPC bilayers. $\Delta G^{\mathrm{cen}-\mathrm{bulk}}$ remains ∼2.5 kcal/mol for O₂, and 1.5 kcal/mol for CO₂, as in the pure POPC bilayers (Figure 7B and D). However, gas permeability in the pure PSM bilayers becomes appreciably smaller (by a factor of 10) than in the pure POPC and binary POPC:CHL bilayers. The calculated $p_{\mathrm{d},\mathrm{mem}}$ values are ∼1 cm/s for both O₂ and CO₂ (Figure 6A-B).

At 50 CHL mole%, the calculated $p_{\mathrm{d},\mathrm{mem}}$ for O₂ is ∼2 cm/s which is two times higher than CHL-free membranes. It is still at least five times lower than any of the simulated POPC bilayers. With the $\Delta G_{max}$ at the tail regions,

$\Delta G^{\mathrm{cen}-\mathrm{bulk}}$ is ∼3 kcal/mol (Figure 7B). For CO₂, the calculated $p_{\mathrm{d},\mathrm{mem}}$ value is ∼1 cm/s (Figure 6B), while $\Delta G^{\mathrm{cen}-\mathrm{bulk}}$ increases to ∼2.75 kcal/mol (Figure 7C).

3.2.3. POPC:PSM:CHL and POPC:DPPC:CHL Bilayers

POPC contains an unsaturated fatty chain and a saturated one, whereas PSM contains a saturated fatty chain and a sphingosine chain in the trans conformer. The saturation of lipid tail could modulate the permeation of gas molecules, as the pure PSM bilayers appear to resemble a gel phase under the simulated temperature of 310 K. We explored this potential effect by simulating CO₂ dynamics in a new set of membrane systems constituted of either POPC and PSM or POPC and DPPC (fully saturated). This new set included 50:50:0, 33:33:33 and 25:25:50 POPC:PSM:CHL, and 50:50:0 and 33:33:33 POPC:DPPC:CHL bilayers.

Similar to the other simulated systems, the highest solubility region for the gas molecules is the center of the membrane. Like the simulated CHL-free bilayers, no local $\Delta G_{max}$ is apparent in the tail regions. $p_{\mathrm{d},\mathrm{mem}}$ are 8.4 cm/s and 14.3 cm/s in the 50:50 POPC:PSM and POPC:DPPC bilayers, respectively (Figure 6B-C), in correlation with the differences in lipid order and packing density. At 33 and 50 CHL mole%, CO₂ occupancy is reduced in the tail region but to a lesser degree than the 50:50 PSM:CHL bilayer. CO₂ permeability was also higher than in the pure PSM and 50:50 PSM:CHL bilayers. $p_{\mathrm{d},\mathrm{mem}}$ in the 33:33:33 POPC:PSM:CHL and POPC:DPPC:CHL bilayer were 6.5 cm/s and 7 cm/s, respectively (Figure 6B-C), in correlation with lipid occupancy profiles (Figure 2H and K). The value in the 25:25:50 POPC:PSM:CHL bilayer was 4.3 cm/s.

We also simulated CO₂ diffusion in a DPPC bilayer with 50 CHL mole% and analyzed its last 400 ns trajectory for $p_{\mathrm{d},\mathrm{mem}}$ and $\Delta G$ calculations. The calculated $p_{\mathrm{d},\mathrm{mem}}$ value was ∼3 cm/s (Figure 6C), which is lower than those in the simulated POPC-containing bilayers but higher than the pure PSM and 50:50 PSM:CHL bilayers. $\Delta G^{\mathrm{cen}-\mathrm{bulk}}$ was ∼2.25 kcal/mol, which is 0.5 kcal/mol lower than in the 50:50 PSM:CHL bilayer (Figure 7E).

3.3. Diffusivity of Gas Molecules

Following the inhomogeneous solubility diffusion model [3,75,76,77], membrane permeability of a gas molecule is modulated by its solubility or partitioning in the membrane, its diffusivity, and the membrane thickness. $p_{\mathrm{d},\mathrm{mem}}$ is expressed as:

$${\displaystyle \frac{1}{p_{\mathrm{d},\mathrm{mem}}}}={\int }_{-d}^{+d}{\displaystyle \frac{\mathrm{d}z}{D\left(z\right){\mathrm{e}}^{-\Delta G\left(z\right)/\mathrm{R}T}}},$$

where $\Delta G\left(z\right)$ is the local partitioning free energy and $D\left(z\right)$ is the local translational diffusion coefficient in membrane segment (dz). $+d$ and $-d$ corresponds to the upper and lower membrane part for which diffusivity is calculated, respectively; the distance resembles the membrane height, denoted as $\Delta z$ in Table 2. $exp-\Delta G\left(z\right)/\mathrm{R}T$ is the partitioning or solubility coefficient of a gas molecule with respect to the aqueous solution and can be calculated from the partitioning $\Delta G$ profiles (Figure 7). With the calculated $p_{\mathrm{d},\mathrm{mem}}$, the overall diffusion coefficients of a gas molecule in the membrane ($D_{\mathrm{app}}$) can be approximated to:

$$D_{app}=p_{\mathrm{d},\mathrm{mem}}{\int }_{-d}^{+d}{\displaystyle \frac{\mathrm{d}z}{e^{-\Delta G\left(z\right)/\mathrm{R}T}}}.$$

For all of the simulated POPC containing bilayers (i.e., pure POPC, POPC:CHL binary, and POPC:PSM:CHL and POPC:DPPC:CHL ternary bilayers), as well as for 50:50 DPPC:CHL, the decrease in the permeability of O₂ and CO₂ is resulted by the decrease of their membrane solubility. The estimated $D_{\mathrm{app}}$ values for O₂ and CO₂ in these membranes are similar, ranging from $2\times {10}^{-6}$ to $3.5\times {10}^{-6}$ cm²/s (Figure 8), which is about 10 times slower than their diffusion in water ($\sim 3\times {10}^{-5}$ cm²/s).

In the 50:50 PSM:CHL bilayers, $D_{\mathrm{app}}$ are $\sim 1.5\times {10}^{-6}$ cm²/s (Figure 8A-B). Gas diffusion is even slower in the pure PSM bilayers with $D_{\mathrm{app}}$ values of $\sim 5\times {10}^{-7}$ cm²/s (Figure 8A-B), which is at least 4 times slower than those in the simulated POPC-containing bilayers.

4. Discussion

The present MD study examines the extent to which changes in the lipid compositions, such as the CHL and SM contents and the degree of lipid saturation, influence the permeability of bioactive O₂ and CO₂ gases across the membrane. Our results shows that increasing CHL content results in an increase in lipid order, packing density and membrane thickness, in correlation with a reduced gas solubility in the membrane. The degree of lipid saturation (i.e., presence of POPC, PSM and/or DPPC) also contributes to the rate of gas permeation. Among the simulated bilayers, gas permeability is slowest is a pure PSM bilayer ($p_{\mathrm{d},\mathrm{mem}}\sim 1$ cm/s), and the incorporation of monosaturated POPC lipids increases the permeability by several folds.

We note that our simulation used a temperature of 310 K (body temperature). According to an X-ray diffraction study, the ripple-fluid phase transition temperature of the PSM bilayer is at 314 K [78]. This explains the observed undulated membrane structure at the simulated temperature (Fig. Figure 3). This also results in the dynamics of lipid molecules in the pure PSM bilayer to be slower than those in the simulated POPC or CHL bilayers, which are either in a liquid-ordered or liquid-disordered phase or coexist in both phases [44,79]. Low membrane fluidity, together with tight packing of PSM lipids, rationalizes the appreciably lower calculated $p_{\mathrm{d},\mathrm{mem}}$ and $D_{\mathrm{app}}$ values for O₂ and CO₂ in the pure PSM bilayers. The apparently larger variations of the $\Delta G$ profiles in the PSM bilayer than in the CHL-free POPC, POPC:PSM and POPC:DPPC bilayers (Figure 7) are most likely due to the coexistence of solid and liquid lipid domains in the ripple phase.

We also show that the presence of CHL in the membrane limits the permeability of CO₂ more significantly than O₂ (Figure 6A-B). These findings agree with a recent MD study [57], observing only 23% reduction of O₂ permeability at 62.5 CHL mole%. CHL partitioning in the membrane was also found to strengthen the hydrophobic lipid-lipid contacts [55], thereby tightening the membrane. This may make CO₂, which has a larger volume, more difficult to diffuse than O₂. This observation also explains the lower $p_{\mathrm{d},\mathrm{mem}}$ value for CO₂ than O₂, as CO₂ shows a greater decrease in its solubility in the lipid-tail region occupied by CHL molecules (Figure 7A and C). Since the pure POPC and other simulated CHL bilayers have similar $D_{\mathrm{app}}$ vlaues (Figure 8), gas permeation through fluid bilayers still follows the Meyer-Overton rule, which predicts their permeability from their solubility [29].

In all of our simulated bilayers, gas concentration is highest in the midplane of the membrane. $\Delta G$ barriers for the entry of gas molecules from the aqueous solution into the membrane is negligibly small. It is the diffusion of gas molecules from the midplane to the solution that constitutes the main barrier of their permeation. These features could play vital roles in a living cell by creating a local concentration gradient of gases that are substrates of many biochemical reactions taking place within the membrane. A notable example is the catalysis of O₂ to water by cytochrome c oxidases of which the entrance of the O₂ access pathway is located near the center of the membrane [80].

Many aspects potentially controlling membrane permeability remain to be explored, and this study focuses only on neutral phospholipids (POPC, DPPC and PSM) and CHL. Realistic biological membranes also contain anionic and more complex lipids, as well as integral and peripheral proteins and adhesion proteins. These other components, especially proteins, are structurally more condensed and less dynamic than lipids, so their presence could add significant barriers for gas molecules to penetrate in. In such conditions, membrane channels, such as aquaporins which contain high-affinity, always-open pathways for gases [81,82], may be important for maintaining cellular homeostasis and reducing cellular toxicity.