To further characterize the effect of Ca 2+ ions on membrane properties, we computed the lateral pressure proiles (LPPs) across the lipid bilayers. Considering that PS lipids are mostly located in the inner lealet of biological membranes, we set up a double-membrane and asymmetric simulation system as shown in Fig. , in which an asymmetric bilayer system was duplicated along the z direction. Please note that such a double-membrane system resulted in two isolated water compartments, as the top and bottom water regions were actually connected into one compartment due to the periodic boundary condition used in the simulations. K + ions were used to neutralize the negative charges generated from PS and evenly distributed in the two water compartments, while different amounts of Ca 2+ ions were added to the two water compartments to generate an asym- metrical Ca 2+ distribution. Cl -ions were added to ensure that no charge imbalance was generated between the two water compartments. Thus, both the ion distribution and composition of lipid bilayers were asymmetrical, which can help to characterize the impact of asymmetrical Ca 2+ distribution on biomembrane properties. For instance, the sarcoplasmic reticulum (SR) membrane faces a higher inner Ca 2+ concentration than the outer, which resembles the upper lipid bilayer in Fig. , while the outer lealet of plasma membranes may face a higher Ca 2+ concentration than the inner lealet that enriches POPS, which resembles the lower lipid bilayer in Fig. . Note that in our MD simulations, we used an NPT-ensemble such that the normal pressure p zz is kept almost constant along the z-axis (Fig. , black curve). This was achieved by adjusting the density difference of water across the membrane such that the osmotic pressure caused by Ca 2+ density gradient is just balanced off by the opposite water density gradient. The same strategy was also used in the previous study. The LPPs of the asymmetrical Ca 2+ treated bilayer lipids are shown in Fig. . As can be seen, moving from the bulk solvent toward the bilayer center, LPPs displayed a moderate positive pressure peak (I) attributed to the repulsion of the lipid headgroups. A large negative peak (II) was observed in each lealet, indicative of the lipid-water interfacial tension, followed by another positive peak (III) close to the center of the bilayer, which was demonstrated to align well with the location of the acyl chain double bonds and caused by the relative incompressibility of the double bonds in the unsaturated lipids. In addition, LPP had a positive peak at the bilayer center (IV), which was related to the geometries of the speciic lipid types. We calculated the differential stress for each lealet by integrating the lateral pressure π(z) along z from the bulk water to the bilayer center (Fig. ). The differential stress of the two lealets was then obtained and denoted as ∆ σ 1 and ∆ σ 2 (Fig. ). We analyzed the lateral pressure proiles π(z) in simulation systems without Ca 2+ to show the effect of asymmetric lipid compositions (Fig S6 ). The LPP was also aligned to the bilayer centers to give an easier visual comparison across bilayers (Fig. and). The ∆ σ 1 and ∆ σ 2 were 0.06 ± 0.68 mN/m and -1.05 ± 0.72 mN/m in this system, indicating that the lipid composition asymmetry in our simulations system did not lead to a signiicant differential stress between the two lealets. Therefore, the differential stress shown in Fig. was mostly caused by the addition Ca 2+ . As shown in Fig. , the addition of Ca 2+ had a minor effect on ∆ σ 1 (black dotted line), but the impact was signiicant for ∆ σ 2 (black solid line), which was increased up to -6.37 mN/m while increasing Ca 2+ concentration on the POPSrich side. This indicates that the addition of Ca 2+ to the negatively charged lealet of membranes will lead to a signiicant change in the differential stress in this lealet, and consequently, a stress imbalance will be generated in the two lealets of the membrane. Such a stress imbalance would cause a spontaneous curvature of the membrane. We calculated the curvature order parameter using the SuAVE software to explore the effects of Ca 2+ on membrane shapes. The results showed no signiicant difference in the average order parameter across the simulation systems with various Ca 2+ concentrations, and the values were consistent with those calculated for planar bilayers (Fig. ). This indicates that no signiicant curvature was generated by the presence of Ca 2+ in the MD simulations, probably owing to the inluence of the periodic boundary condition and the limited size of the simulation systems. Therefore, we had to use other methods to analyze the shape change of a continuous membrane.