Efficient Simulations of Solvent Asymmetry Across Lipid Membranes Using Flat-Bottom Restraints

The routinely employed periodic boundary conditions complicate molecular simulations of physiologically relevant asymmetric lipid membranes together with their distinct solvent environments. Therefore, separating the extracellular fluid from its cytosolic counterpart has often been performed using a costly double-bilayer setup. Here, we demonstrate that the lipid membrane and solvent asymmetry can be efficiently modeled with a single lipid bilayer by applying an inverted flat-bottom potential to ions and other solute molecules, thereby restraining them to only interact with the relevant leaflet. We carefully optimized the parameters of the suggested method so that the results obtained using the flat-bottom and double-bilayer approaches become mutually indistinguishable. Then, we apply the flat-bottom approach to lipid bilayers with various compositions and solvent environments, covering ions and cationic peptides to validate the approach in a realistic use case. We also discuss the possible limitations of the method as well as its computational efficiency and provide a step-by-step guide on how to set up such simulations in a straightforward manner.


Estimation of Relative Speed
The speed gain obtained using the flat-bottom approach compared to the double-bilayer one was estimated using both a workstation computer and a supercomputer.The workstation was equipped with an AMD Ryzen 9 5950X 16-core processor, an nVidia GeForce RTX3080 Ti GPU, and 128 GB of RAM.On the supercomputer, four nodes were used, and each was equipped with two AMD Rome 7H12 64-core processors and 256 GB memory.No GPUs were used in the calculation on the supercomputer.For the benchmark, we ran all the systems listed in Table 1 in three replicates and used their mean as a representative speed estimate.
The first 10,000 steps were discarded due to the PME optimization, after which the speed was recorded for 500,000 steps (1 ns).Finally, we calculated the ratios of the simulation speeds obtained with the flat-bottom and double-bilayer approaches, which are reported in Table 1.For the larger test system containing 2,304 DMPC lipids, 576,000 waters (for the flat-bottom approach, an additional water layer of 0.6 nm was included), and 150 mM NaCl, we used 10 nodes on the supercomputer for our tests.

Supplementary Results
We provide additional data that compare the results from flat-bottom and double-bilayer simulations and confirm the viability of the former approach.We also additionally discuss how the flat-bottom simulations behave if chloride anions are free of restraints.

Chloride Ions Free of Restraints
Since chloride ions do not really interact with lipid membranes, one tempting approach is to include them on both sides of the plasma membrane unrestrained, therefore allowing them to move freely between two solvent environments.We tested this idea on some of our systems and identified several issues with this approach.First, the ionic concentrations of cations and anions become imbalanced since their accessible volume is different, being larger for chloride ions, see density profiles in Fig. S7.This leads to a lower effective concentration of chloride compared to cations.Second, unrestrained chlorides cause artificial water orientation in the vicinity of the flat-bottom excluded volume, Fig. S8.This artifact is the consequence of a little accumulation of cations in this region, which is not balanced by the same accumulation of chloride ions that do not feel the flat-bottom potential.Therefore, we do not recommend using this hybrid approach in flat-bottom restraints, although all membrane and interfacial properties are unaffected by the free motion of chloride ions, see Figs.S7 and S9, and Table S4.
Table S1: Lipid compositions of the inner and outer leaflets of the realistic membrane S1 used in simulations of "Plasma-NaK", "'Plasma-Asym", and "Plasma-Pept" setups.The composition is slightly different from the original one due to the availability of corresponding lipids in CHARMM-GUI    Table S3: Comparison of the membrane properties (area per lipid (APL), membrane thickness D P−P , and the tilt angle of the P-N vector) calculated from R9 peptides-containing simulations with flat-bottom (FB) and double-bilayer (DB).The two values of APL are given for the inner and outer leaflets, respectively.The P-N tilt angle is calculated for phosphatidylcholine lipids only and also reported for the two leaflets separately.The error estimate shows the standard error obtained using block averaging.Table S4: Comparison of the membrane properties (area per lipid (APL), membrane thickness D P−P , and the tilt angle of the P-N vector) calculated from simulations with flat-bottom with free chloride ions (FB-freeCl) and double-bilayer (DB) setups.For asymmetric membranes, the two values of APL are given for the inner and outer leaflets, respectively.The P-N tilt angle is calculated for phosphatidylcholine lipids only and also reported for the two leaflets separately.The error estimate shows the standard error obtained using block averaging.

Figure S1 :
Figure S1: Density profiles for Na + with the optimal parameters for (r=0.3 nm, k =10,000 kJ•mol −1 •nm −2 ) and with various amounts of extra water added.The solid black line refers to a simulation without flat-bottom restraints.

Figure S3 :
FigureS3: The z coordinate of the center of mass of the membrane compared to the z dimension of the box size as a function of simulation time from flat-bottom simulations of (A) "POPC-NaK-", (B) "POPC-Xtra-", and (C) "Plasma-NaK-" systems.

Figure S4 :
FigureS4: Normalized water orientation from flat-bottom (markers) and double-bilayer (solid lines) simulations of (A) "Plasma-Asym-" and (B) "Plasma-Pept-" systems calculated as the cosine of the angle between water dipole and z axis multiplied by the number density of water oxygens.

Figure S6 :
FigureS6: The z coordinate of the center of mass of the membrane compared to the z dimension of the half box size as a function of simulation time from flat-bottom and doublebilayer simulations of(A) "Plasma-Asym-" and (B) "Plasma-Pept-" systems.

Figure S8 :
FigureS8: Normalized water orientation from flat-bottom with free chloride ions ("f.-b.& free Cl − ") and double-bilayer simulations of (A) "POPC-NaK-", (B) "POPC-Xtra-", and (C) "Plasma-NaK-" systems calculated as the cosine of the angle between water dipole and z -axis multiplied by the number density of water oxygens.

Table S2 :
The number of ion crossings across the flat-bottom potential in the salt solution simulations with different values of k and r.