Effects of lipid composition on membrane permeation†

Passive permeation through lipid membranes is an essential process in biology. In vivo membranes typically consist of mixtures of lamellar and nonlamellar lipids. Lamellar lipids are characterized by their tendency to form lamellar sheet-like structures, which are predominant in nature. Nonlamellar lipids, when isolated, instead form more geometrically complex nonlamellar phases. While mixed lamellar/nonlamellar lipid membranes tend to adopt the ubiquitous lamellar bilayer structure, the presence of nonlamellar lipids is known to have profound effects on key membrane properties, such as internal distributions of stress and elastic properties, which in turn may alter related biological processes. This work focuses on one such process, i.e., permeation, by utilising atomistic molecular dynamics simulations in order to obtain transfer free energy profiles, diffusion profiles and permeation coefficients for a series of thirteen small molecules and drugs. Each permeant is tested on two bilayer membranes of different lipid composition, i.e., purely lamellar and mixed lamellar/nonlamellar. Our results indicate that the presence of nonlamellar lipids reduces permeation for smaller molecules (molecular weight < 100) but facilitates it for the largest ones (molecular weight > 100). This work represents an advancement towards the development of more realistic in silico permeability assays, which may have a substantial future impact in the area of rational drug design.

more complex geometries, depending on the specific lipid types and thermodynamic conditions, such as inverse hexagonal structures characterized by lipid-lined water channels. Nonlamellar phases are uncommon, and typically only appear in transient processes such as membrane fusion. Notably, in vivo cell membranes are mixtures comprising both lamellar and nonlamellar lipids, even though the overall structures formed are almonst invariably lamellar [2][3][4][5] . In this study, we consider the representative lamellar lipid dioleoylglycerophosphocholine (DOPC) and the representative nonlamellar lipid dioleyolglycerophosphoethanolamine (DOPE), whose chemical structures are reported in figure 1. Under biological conditions, when dispersed in an aqueous environment, pure DOPC forms lamellar phases, while pure DOPE forms nonlamellar (inverse hexagonal) phases 6 . Structurally, DOPC and DOPE are identical apart from the head terminal group, which is choline for DOPC and ethanolamine for DOPE (see figure 1). This difference in the aminoalcohol affects the size and shape of the polar head, making DOPC lipids bulkier.
The effect of different lipids on the mechanical properties and dynamic behaviour of bilayer membranes has been studied extensively in the past using experiments 7,8 , analytical theory 9,10 and molecular dynamics (MD) simulations [11][12][13][14][15][16][17][18][19][20][21][22] . In particular, Ding et al. 21 have employed atomistic simulations to quantify the effects of changes in the ratio of lamellar vs. nonlamellar lipids on the physical properties of bilayer membranes. A key finding was that the addition of DOPE lipids to DOPC bilayers had a negligible effect on the structure of the bilayers but induced substantial changes to the lateral pressure profile, in agreement with the available qualitative experimental evidence 7 . The lateral pressure profile, arguably the most fundamental physical property of lipid bilayers 2 , characterizes the net stresses across membranes, thus quantifying the pressure that membrane inclusions (such as proteins) or permeants "feel" inside the bilayer 23 . This pressure varies over a vast range (hundreds of atmospheres) with the bilayer depth, and three distinct regions can be identified depending on the nature of the forces present; a region of high repulsion corresponding to the upper lipid headgroups, a region of high attraction in the polar/apolar interface, and a region of varying repulsion in the bilayer hydrophobic core. The lateral pressure profile underlies many fundamental membrane phenomena, such as phase transition 24,25 , water penetration 26 , drug transport 27,28 , anaesthesia [29][30][31][32] and in general membrane protein function 23,[33][34][35][36][37][38] . For instance, changes in the lateral pressure profile are believed to control the opening and closing of ion channels 33 .
With regards to lamellar/nonlamellar mixtures, we are aware of only two relevant previous studies, both experimental. Huster et al. 56 observed that the addition of nonlamellar DOPE lipids reduced the permeation of water by 40 % in comparison to a pure DOPC bilayer and by 18 % when mixed with a pure polyunsaturated 18:1-22:6 PC bilayer. The results were explained based on lipid order increase, which is more prominent in saturated and monounsaturated lipids, leading to tighter packing and smaller area per lipid. Purushothaman et al. 70 examined the permeation of the antibiotic norfloxacin through pure and mixed DOPC, DOPE and DOPG bilayers. The addition of DOPE lipids generally reduced the permeation coefficient, although the reduction was higher for smaller concentrations of DOPE.
In this work, we use molecular dynamics simulations to characterize the effects of a change in the lamellar/nonlamellar lipid composition on passive permeation by simulating bilayers comprising either pure DOPC (lamellar) or a 1 to 3 mixture of DOPC to DOPE (nonlamellar) lipid molecules, respectively. Bilayers comprising different relative concentrations of these two lipid were shown previously to differ substantially in their lateral pressure profile 7,21,71 . In particular, an increase in repulsive forces was observed in the hydrophobic part of the bilayer when DOPE (nonlamellar) lipids were present. We hypothesise that the increase in repulsive forces should lead to an increase in the resistance to permeation through the membrane. In the remainder of the paper, we report a series of results produced with atomistic molecular dynamics simulations for 13 small molecules and drugs, namely ammonia, water, fluoromethane, carbon dioxide, propane, ethanol, urea, isopropanol, glycine, phenol, benzoic acid, coumarin and paracetamol. The chemical structure of all permeants is shown in figure 2 while their physical properties are reported in the supplementary material (table S1 ).

Permeability calculations
A simple framework to predict membrane permeability, first introduced over a hundred years ago by Meyer and Overton 72,73 , is based on the octanol-water partition coefficient logP oct/water of a permeant, which indicates the permeant's solubility preference between an octanol and a water phase. On this basis, the bulk solubility-diffusion model of permeability was proposed 74 in which lipid membranes were considered as homogeneous bulk bodies. Diamond et al. 75 later accounted for the heterogeneity of membranes by developing the inhomogeneous solubilitydiffusion model, whereby the permeation coefficient P of a solute through a membrane can be predicted as 76,77 : where T is the simulation temperature, R is the universal gas constant (which is equal to the product of Boltzmann's constant k B with Avogadro's number N A , R = k B · N A ) and z is the position normal to the bilayer surface, with z 1 and z 2 representing the bulk water regions on the two sides of the membrane. Also, R(z) is the local resistance to permeation, D(z) is the local diffusion coefficient of the solute and ∆G(z) is the Gibbs free energy difference between the thermodynamic states of the permeant in bulk water and at position z. To obtain ∆G(z) and D(z) from molecular dynamics simulations, several enhanced sampling approaches have been developed, such as the z-constraint or the z-restraint methods, which are discussed in detail elsewhere [78][79][80] . In this work we utilise the z-restraint method. The free energy difference can be obtained by using the umbrella sampling scheme 81 , whereby a harmonic potential restrains the movement of the permeant in a small "window" around each position along the reaction coordinate path, which is the z-direction normal to the bilayer plane in our work. The free energy difference is then calculated as: where V b (z) is the biasing potential and P b (z) is the permeant's spatial distribution along z positions. Finally, in order to obtain the unbiased free energy difference ∆G, the weighted histogram analysis method (WHAM) is used 82,83 .
Regarding the diffusion coefficient (D(z) in equation 1), it should be noted that, in general, computing reliable local diffusion coefficients from restrained simulations remains a very active research field [84][85][86][87] . One of the most popular methods is the one introduced by Hummer 88 (in turn based on the previous works of Woolf and Roux 89 and Berne et al. 90 ). In this method, when umbrella sampling simulations are performed with a harmonic bias along a reaction coordinate, the diffusion coefficient can be computed as: where var(z) = z 2 − z 2 is the variance of the z-distance between the centres-of-mass of the permeant and membrane and τ is the characteristic time of the z-distance autocorrelation function: where according to the definition of the autocorrelation function δ z(t) = z(t) − z . The integral can be computed with the method proposed by O'Neill et al. 91 , where the integration domain is taken from the beginning until the first time that the autocorrelation function becomes zero. Finally, permeability coefficients are computed by direct substitution into equation 1.
It must be noted that permeation through water pores has also been proposed as an alternative or complimentary pathway to the solubility-diffusion mechanism. 92 In particular, the importance of pore mediated permeation has been recently established for large, charged molecules such as cell penetrating peptides. 93,94 However, for smaller and neutral molecules, such as those considered in this work, the solubility-diffusion mechanism is regarded as predominant. 95-97

Simulation Protocol
Each simulation conducted in this study comprised 4300 water molecules, 128 lipid molecules (64 per leaflet) and 1 permeating molecule (see figure 2 for the full set of permeants considered). Each permeant was simulated in two membrane systems, pure (comprising DOPC only) and mixed (DOPC:DOPE (1:3) mixture). A typical snapshot for the mixed system is reported in figure 3. Both starting configurations for the membrane systems were taken pre-equilibrated from our previous work 21 , and the permeant was manually placed in the required position along the bilayer. Initially, an energy minimisation was performed to remove any high energy overlaps and then a short constant temperature and pressure (NPT) 100 ps equilibration was run. Both during the minimisation and the equilibration, the distance between the centres of mass of the membrane and the permeant was constrained, in order to ensure the correct distance between them in the beginning of the production simulation. Overall, 28 z-positions were examined for each permeant-membrane combination, and for each such case a 100 ns NPT production simulation was performed. Molecular dynamics simulations were conducted with the GROMACS software 98-103 version 5.1.1. Van der Waals forces were approximated with a Lennard-Jones potential with a switching cutoff from 1 nm to 1.2 nm, short range electrostatics were approximated with a Coulomb potential with a cutoff at 1.2 nm and long range electrostatics were treated with the Smooth Particle-Mesh Ewald (SPME) 104 method. Lipid molecules were modelled with the CHARMM36 (August 2015 version) force field 105,106 , permeant molecules were modelled with CHARMM36 or compatible CGenFF 107,108 parameters and water molecules were modelled with the CHARMM implementation of TIP3P 109 .
For temperature coupling during the equilibration and production simulations, the V-rescale algorithm 110 was used. All the systems were kept at 300 K and the coupling time constant was set to 1 ps. Pressure was kept at 1 bar with the Berendsen 111 and Parrinello-Rahman 112-114 barostats during equilibration and production simulations, respectively; the coupling time constant was 5 ps and the coupling type was semi-isotropic, i.e., isotropic for the x and y directions but independently from the z direction (as is common practice in simulations of interfacial systems 34,94,115 ). A 2 fs timestep was used; covalent bonds with hydrogen atoms were constrained with the SETTLE algorithm 116 for water molecules and with the LINCS algorithm 117,118 for all other molecules.

Analysis
The restraining forces and z-axis fluctuations of the permeant were sampled every 1 ps resulting in two timeseries of 10 5 points each. The first 30 ns (30 000 points) were discarded as extra equilibration time and all final results were produced using the last 70 ns. Standard errors for the calculated properties were computed with the block averaging method 119,120 , where the 70 ns timeseries were separated in 7 blocks of 10 ns. In order to increase the computational efficiency and since the monolayers' composition was the same for each examined membrane, the permeants were positioned only across one leaflet of the bilayer. Afterwards, the position-dependent results were treated as symmetrical to cover the entire z-dimension of the bilayer.
The free energy profiles ∆G(z) were computed with the GRO-MACS implementation 121 of the weighted histogram analysis method 83,122 . Each permeant was restrained by a virtual spring with a harmonic force constant of 1000 kJ mol −1 nm −2 along the z-axis, for 28 positions, every 0.1 nm from the water slab to the bilayer core. Local diffusion coefficients D(z) were computed according to equation 3; extra analysis on the numerical integration of the autocorrelation function and the handling of oscillatory profiles with the application of filters is reported in the supplementary material. The resistance profiles and the permeability coefficients were computed by direct substitution of ∆G(z) and D(z) into equation 1. Statistical significance testing of the differences between permeation coefficients for the pure and mixed membranes was carried out with paired t-tests. Hydrogen bonds were computed between the permeants and the lipid-water molecules. The permeant lateral mobility was investigated qualitatively through the corresponding lateral trajectories; representative traces are reported in the supplementary information.

Free-energy profiles
The permeation free energy profiles ∆G(z) of the selected thirteen permeant molecules, simulated in both the pure DOPC and DOPC:DOPE(1:3) bilayers, are shown in the first column of figure  4.
The convergence study of the profiles is presented in detail in the supplementary material; while general features are analyzed as in previous studies, 123,124 we report a novel approach to measure relative convergence.
The free energy difference represents the thermodynamic forces (entropic and enthalpic) that drive the process of perme-Distance from bilayer centre [nm]  ation, quantifying the spontaneity of the process and the preferred partitioning position of the permeant 125 . Considering our results, permeants can be separated into three distinct categories depending on where they exhibit the global minimum of their free energy difference, which predicts the location where the molecules preferentially partition.
In the first category, the free energy difference is always positive along the bilayer and the minimum position corresponds to the reference point for ∆G(z) = 0 in the water phase. This behaviour is observed for ammonia, glycine, urea and water; these molecules therefore are not predicted to partition inside the bilayer. Such a behaviour is consistent with the well-known polar hydrophilic nature of these compounds.
In the second category, the profiles have a small positive peak in the lipid head region, then exhibit a negative global minimum in the polar/apolar interface and finally a larger positive peak in the hydrophobic core. Coumarin for the DOPC membrane, ethanol and isopropanol, clearly belong to this category, with five more molecules exhibiting minor deviations from this behaviour. Benzoic acid has a peak in the hydrophobic core but is negative in value. Also, coumarin for the DOPC:DOPE(1:3) membrane, as well as paracetamol and phenol for both membranes, do not have a positive barrier in the head region. Fluoromethane has a positive peak in the head region, a global minimum in the interface, and a second positive peak close to the hydrophobic core; inside the tail region, the profile is marginally positive for the DOPC membrane, while marginally negative for the DOPC:DOPE(1:3) membrane. All these molecules have an amphiphatic nature and under physiological conditions they are known to partition in the polar/apolar interface 126,127 , therefore the free energy profiles are in agreement with the expected behaviour.
The third category includes carbon dioxide and propane, whose profiles feature a small positive peak in the lipid head region while being negative elsewhere, with global minima located in the bilayer centre, corresponding to their preferred partitioning position. Such behaviour is consistent with both permeants being known hydrophobic molecules.
In general, all free energy profiles presented here are in good qualitative agreement with previous studies of the same permeants and the same or different PC and PE lipids [45][46][47][48]62,[128][129][130] . Furthermore, they all fit in the categories that were introduced by Neale et al. 131 who classified over 200 free energy profiles of more than 100 small molecules from previous studies.
Permeants can also be classified by examining the effect of DOPE lipids on the transfer process. In the first group, comprising ammonia, glycine, urea and water, the free energy profile of the DOPC:DOPE mixture is the same or higher than the DOPC profile, across the whole bilayer depth. This is an indication that the presence of DOPE lipids increases the barrier to permeation, in parts or across the entire bilayer. For the second group, including carbon dioxide, ethanol, paracetamol and phenol, a lower DOPC:DOPE mixture profile is observed in the head group area but a higher DOPC:DOPE profile is observed in the interface and hydrophobic core area. Isopropanol and fluoromethane show a slight deviation from this behaviour by having a marginally lower peak in the DOPC:DOPE bilayer core than in the pure DOPC mem-brane. Benzoic acid, coumarin and propane form the third group, whereby the DOPC:DOPE mixture profile is the same or lower than the DOPC profile all across the bilayer. A summary of the permeants classification in relation to their free energy difference profiles is reported in Table 1.

Local diffusion profiles
The local diffusion coefficients for all the permeants and membranes investigated are reported in the second column of figure  4. All cases exhibit the same qualitative behaviour; diffusion is largest in the water region, then it drops quickly as permeants approach the hydrophobic lipid tails, and finally increases again in the bilayer centre. Quantitatively, for both membranes, diffusion of heavier molecules is considerably slower, especially in the water region; for example, paracetamol has ten times slower diffusion in water than ammonia. With regards to the effect of DOPE lipids, there are no significant differences between the two examined membranes. For benzoic acid, fluoromethane, paracetamol and phenol, the peak in the bilayer core is marginally higher for the DOPC:DOPE mixture. Also, for all molecules, the diffusion profiles for the DOPC systems tend to be slightly higher than DOPC:DOPE profiles, especially in the water and headgroup regions. It can also be seen that some profiles exhibit pronounced fluctuations, especially in the water and headgroup regions; extra analysis of these effects is reported in the supplementary material. Overall, while diffusion coefficients are inherently noisy, it should be stressed that their contribution to the overall permeation model is marginal compared to the free energy term, which appears as argument of an exponential function (equation 1).

Local resistance profiles
The resistance profiles for both membrane systems and all permeants studied are reported in the third column of figure 4. Carbon dioxide and propane, both hydrophobic molecules, experience higher permeation resistance in the hydrophilic region of the bilayer, while hydrophilic molecules (ammonia, glycine, urea and water) experience the highest resistance in the bilayer core, as expected. Amphiphilic molecules are characterised by two peaks in their resistance profiles. The first features in the hydrophilic area of the bilayer, indicating a permeation barrier of the hydrophobic part of the permeant as it dissolves inside the hydrophilic region of the bilayer. The second can be seen in the lower chains region and the bilayer core, indicating a resistance to permeation of the hydrophilic part of the permeant in the hydrophobic region of the bilayer. The effect of the lipid composition depends on the hydrophilicity of the permeants. The DOPC:DOPE (1:3) bilayer induces higher resistance for hydrophilic permeants than the pure DOPC bilayer, especially in the lipid chains region. For hydrophobic molecules the resistance is marginally higher in the polar region of the DOPC:DOPE (1:3) bilayer but it is unaffected elsewhere. Finally, for the amphiphilic permeants, the DOPC:DOPE (1:3) membrane has lower resistance in the bilayer centre, higher resistance in the polar part and the 0.5 nm to 1.5 nm regions, and the same resistance in the other regions. Coumarin is an exception as the resistance is lower along the whole apolar part of the bilayer.
Overall it can be seen that, except for minor deviations, resistance profiles are qualitatively similar to the free energy profiles (compare first and third column of figure 4). Any deviation can be ascribed to differences in the diffusion behaviour, as reported also in previous permeation studies 46,47,132,133 . From our results, it can be seen that differences in the diffusion profiles can alter the relative resistance between the DOPC and DOPC:DOPE membrane for some permeants. In particular, resistance to permeation of both membranes to ammonia, urea, glycine and phenol becomes equal or higher for DOPC in the bilayer centre, although the free energy profile is clearly higher in the same region for the DOPC:DOPE mixed system. In other cases, e.g., for water, fluoromethane and benzoic acid, the differences between the membranes are amplified.

Permeability
Permeability coefficients for each molecule were computed from the local resistance profiles according to the inhomogeneous solubility-diffusion model as defined in equation 1. The permeation values of the molecules examined are presented in ascending permeability order in table 2 along with selected representative literature results for the same permeants (an extended collection of all the available literature data for the examined molecules can be found in table S4 in the supplementary material).
In order to make comparisons between different molecules easier, the logarithms of base 10, logP, of the permeation coefficients are also presented in table 2 and figure 5. From figure 5, it can be seen that most permeants have a negative logP value, with the most negative corresponding to the four hydrophilic molecules (urea, water, glycine and ammonia). Propane and carbon dioxide, both hydrophobic, have the highest logP together with benzoic acid. Fluoromethane, ethanol, isopropanol, phenol, coumarin and paracetamol have logP values between -1 and 0.7. It can also be seen that the majority of the permeants are on the right side of the equality line, indicating a higher logP value through the DOPC membrane; only paracetamol and coumarin are on the left of the equality line, indicating higher logP value through the mixed DOPC:DOPE bilayer. To evaluate whether the membrane composition has a significant effect on the permeability coefficient, we carried out onetailed paired t-tests on the differences in logP induced by exposing a permeant to the two different membranes. Two tests were performed on two different permeant groups, one including permeants with molecular weight smaller than 100 g mol −1 (smallsized molecules) and one including permeants with molecular weight larger than 100 g mol −1 (medium-sized molecules). For molecules with molecular weight lower than 100 g mol −1 , there is strong evidence (p − value = 0.0019) that the logP values for the Overall, our results are consistent with available experimental and computational data previously reported in the literature.
In particular, Jansen and Blume 60 examined water permeation through several large unilamellar lipid vesicles. At T=343 K, DMPE and DPPE reduced permeation by ≈ 99% in comparison to pure DMPC and DPPC, respectively. The same behaviour was observed in our study, although the relative difference between the DOPC and DOPC:DOPE membrane was smaller (−27%). Also Huster et al. 56 studied water permeation, with O 17 nuclear magnetic resonance, through a pure DOPC and a DOPC:DOPE(1:1) bilayer, at T=303 K. They reported that for the latter, the permeation coefficient reduced by ≈ 39%, consistent with the present work. Seo et al. 141 performed a study on the effect of lipid composition on the passive permeation of molecules through PAMPA assays, and observed an 825% increase in the permeation of benzoic acid through a DOPE phase, in comparison to pure DOPC. Although the PAMPA system differs from the bilayers in our MD simulations, in that the lipid phase in PAMPA is much thicker (≈ 100 µm), it is noteworthy that we predicted the same qualitative effect, i.e., a permeation increase in the DOPC:DOPE mixed bilayer compared to pure DOPC (specifically, we observed a 31% increase). Hub et al. 62 examined the permeation of ammonia through a pure POPC and a pure POPE bilayer at T=300 K. They reported an ≈ 87% reduction in permeation for the latter, which is consistent with our observation of a 71% reduction for DOPC:DOPE(1:3). Wennberg et al. 69 examined the partitioning of ammonia, ethanol and propane in a pure POPC and pure POPE bilayer and they reported that the transfer free energy through the latter was increased for all solutes. Finally, Zocher et al. 48 also observed a decrease in the permeation of ammonia through a pure DOPE membrane in comparison to a pure DOPC.
Overall, our results are mostly in agreement with the permeation and logP values previously reported in the literature. Where permeation coefficients were not available, free energy profiles from the literature can be used to assess the findings. Comparisons of our logP values to literature values for each permeant are detailed in the following paragraphs.
Urea. The computed logP value of −6.17 for the DOPC membrane is close to the −5.85 for DOPC 53 at 303 K and −6.27 for DMPC 130 at 298 K found in the literature.
Water. Studies of water permeation have reported a plethora of results spread over a large range, depending on the experimental and simulation protocols and lipid compositions. In general, logP values of water are between −1.15 to −5.64 and our logP values of −3.40 for DOPC and −3.54 for DOPE belong to the lower end of this range. In particular, values of this study are two orders of magnitude smaller than the DOPC and POPC permeation coefficients of Mathai et al. 59 , Paula et al. 53 , Huster et al. 56 and Koenig et al. 136 , one order of magnitude smaller than those reported by Olbrich et al. 135 and Comer et al. 85 and very close to the values of DOPC from Carruthers et al. 134 and DMPC/DPPC of Jansen and Blume. 60 Glycine. Chakrabarti et al. 137 reported a logP of −10.7 through large unilamellar vesicles of DMPC, which is somewhat smaller than our result for DOPC.
Ammonia. Our MD simulations underestimated the permeation coefficients of ammonia in comparison to the rest of the literature. With regards to studies of similar lipid compositions, the values presented here are 1 to 2 orders of magnitude smaller than the POPC, POPE and DOPC MD simulations of Hub et al. 62 and Zocher et al. 48 The discrepancy with both studies may be ascribed to the different force fields used; specifically, the referenced studies employed the Berger united-atom model for lipids, OPLS for ammonia and TIP4P for water, whereas we used allatom CHARMM, CGenFF and TIP3P, respectively (as detailed in the Simulation Protocol section).
Ethanol. Literature values vary over 6 orders of magnitude and they are mostly negative, apart from the simulation work of Comer et al. 84 that reported a logP value of 0.30 for POPC at T=308 K. Ghaemi et al. 138 computed a logP of −1.07 for POPC at T=323 K which is very similar to our values. The experimental work of Ly and Longo 139 , with SOPC, returned a logP of −4.42 at room temperature. Further validation results were available from the works of MacCallum et al. 50 and Carpenter et al. 142 that studied the free energy of ethanol permeating through a DOPC bilayer. In both cases, our findings are almost identical to theirs, qualitatively and quantitatively. Therefore it is expected that also the permeation coefficients would be very similar.
Benzoic Acid. Our results are in good agreement with the simulation work of Lee et al. 130 who computed a logP of 0.45 for DMPC at room temperature, similar to our value of 0.80 for DOPC. These coefficients are both many orders of magnitude higher than the experimental findings for DOPC 140,141 and DOPE 141 .
Carbon dioxide. The logP obtained in our work is very similar to the one by Hub et al. 62 , also from MD simulations, for carbon dioxide permeation through a POPC:POPE bilayer.
Propane. While no reported permeation coefficients were found for propane, MacCallum et al. 50 have presented the free energy profile of permeating propane. As with ethanol, our findings are in excellent agreement with theirs, observing the same −3.2 kcal mol −1 free energy trough in the bilayer centre and overall qualitative behaviour. Therefore, it is expected that the permeation coefficient would also be identical.
Isopropanol, fluoromethane. No computational or experimental results were found regarding the permeation of isopropanol or fluoromethane. Ours appears to be the first report of the permeation coefficients for these two molecules.
Phenol, paracetamol and coumarin. While there are no literature permeation studies with similar lipid compositions, Paloncỳovà et al. 143 have reported free energy profiles of coumarin in a DOPC bilayer at 310 K that are qualitatively close to our findings, with a global minimum around 1.4 nm and an energy barrier in the centre of the bilayer. However, quantitatively the results differ considerably, possibly because of the difference 1-13 | 9 in temperature and force field, as Paloncỳovà et al. 143 used the united atom Berger model for lipids and custom parameters for coumarin based on GROMOS 53a6. In our study, the global minimum is −1.5 kcal mol −1 , while they report a global minimum of −5.7 kcal mol −1 to −6.7 kcal mol −1 . Likewise, in the bilayer centre they report a −3.5 kcal mol −1 to −4.1 kcal mol −1 negative peak in contrast to the positive 1.2 kcal mol −1 barrier that was found in this work. Therefore, the permeation coefficient for coumarin from Paloncỳovà et al. 143 would be probably a few orders of magnitude higher than the one reported here.

Hydrogen bonds
The average number of hydrogen bonds per frame (h.b.p.f.) formed by the permeants is shown in the last column of figure  4. In particular, the total number of hydrogen bonds are shown together with their decomposition into the number of hydrogen bonds formed between the permeant and either the bilayer or the solvent (water). Hydrogen bonds with the solvent are predominant in comparison to the ones formed with the bilayer, and generally increase with the distance from the bilayer centre as also observed in previous studies 79,127 . It is noteworthy however, that hydrophilic and amphiphilic molecules form hydrogen bonds with water almost right down to the bilayer core, showing the propensity of permeants to retain their hydration shell even deep into the hydrophobic region 144 . Quantitatively, as expected, hydrophilic molecules formed more hydrogen bonds with water compared to hydrophobic molecules. Hydrogen bonds with the bilayer can be split in two profile categories. For fluoromethane, carbon dioxide, coumarin and ammonia/DOPC, they are very close to zero along the bilayer. For the rest, hydrogen bonds start to form around 0.1 nm to 0.5 nm, then reach a maximum at 1.5 nm to 2.0 nm (of 0.4 h.b.p.f. for water, urea, glycine and paracetamol, of 0.2 for ethanol, isopropanol, phenol and benzoic acid, of 0.1 for ammonia/DOPC:DOPE) and then fade back to zero in the solvent region.
Overall, no significant difference was observed in the hydrogen bond formations of the molecules with the bilayer between the DOPC and DOPC:DOPE(1:3) membranes, i.e., the lipid composition did not affect the total number of hydrogen bonds formed between the permeants and their environment. Such a finding may appear counterintuitive, due to the presence of an extra hydrogen bond donor in DOPE (part of the headgroup amine moiety) compared to DOPC; however, as was observed in the past 11 , any corresponding additional DOPE hydrogen bonds are mostly formed with either other DOPE or DOPC headgroups, or with solvent molecules, so that permeants are excluded.

Conclusions
In this work we used atomistic molecular dynamics simulations to examine how changes in the lipid composition can affect the biologically crucial process of passive permeation. The free energy, local diffusion and local resistance profiles, as well as the permeation coefficients were reported for thirteen small molecules and drugs. For eleven of the molecules selected, permeation coefficients were reported for the first time in relation to DOPC and DOPC:DOPE(1:3) membranes. Our key findings are that the presence of the nonlamellar DOPE lipids reduces permeation for the smaller molecules (molecular weight < 100 g mol −1 ) while enhances permeation for the largest (molecular weight > 100 g mol −1 ). This is the first systematic investigation of the effect of changes in the lamellar vs. nonlamellar lipid composition on membrane permeation. In general, our study represents an advancement towards the development of more realistic membrane models for the in silico prediction of molecular permeation.

Conflicts of interest
There are no conflicts to declare.