After applying 1500 steps of energy minimization, the simulation system was equilibrated for about ~2 ns with increasing time steps (from 1 fs to 20 fs, see Table for more details). The production run for each replica was 20 µs. During the equilibration and production runs, LINCS 185 constraints were used for ring systems and stiff bond constraints; the LINCS order and iteration were set to 12 and 2, respectively. The temperatures of lipids and solvent (NA, CL, and W beads) were separately kept at 310 K with a time constant for coupling (τT) of 1.0 ps, using the Berendsen 187 and v-rescale 198 thermostats during the equilibration and production, respectively. The pressure was kept at 1 bar initially using Berendsen 187 and Parrinello-Rahman 188 barostats during the equilibration and production steps, respectively. The compressibility was set to 3 × 10 -4 bar -1 . The cut-off distance for the short-range neighbor list was set to 1.35 nm for proper neighbor list updates, as suggested by z-coordinate positional restraints were applied to the phosphate beads of POPC lipids only in the upper leaflet, with a force constant of 2 kJ/mol/nm 2 , to maintain a flat bilayer as explained in earlier work .