We performed MD simulations for membrane systems in a solution of CaCl 2 . The membranes were composed of either pure POPC or a mixture of POPC:POPS (3:1) (Fig. ), and both the recently developed multisite Ca 2+ model and the default Ca 2+ model of the CHARMM force ield were used in the simulations for comparison, which were denoted Cam and Cal, respectively. The Ca 2+ concentration was raised stepwise from 0 to 50, 100, 200, 500 and 1,000 mM. The majority of the membranebound Ca 2+ ions bind to the phosphate groups of POPC and POPS, or the carboxylate groups of POPS lipids in simulations (Fig. ). Then, we analyzed the distribution of various groups along the membrane's normal direction in the equilibrated simulation systems using the SuAVE software. The density proiles were similar with various CaCl 2 concentrations, so we only showed the results obtained for two representative cases: 1 M CaCl 2 ± a high-concentration regime (Fig. ) and 200 mM CaCl 2 ± a low-concentration regime (Fig. ). Simulations with the two Ca 2+ models both showed that the density peak of phosphate groups, often viewed as the membrane boundary, was located approximately 2 nm away from the membrane center (Fig. pink lines). The density peak of Ca 2+ roughly overlapped with that of phosphate groups, indicating the adsorption and accumulation of Ca 2+ under the membrane surface. Notably, more Ca 2+ was observed around this region in the POPC/POPS system than in the pure POPC system (Fig. , yellow lines), which is reasonable since POPS is negatively charged and has a stronger binding afinity to Ca 2+ . There were some differences in the distributions of Ca 2+ around the membranes when using the two Ca 2+ models. As shown in Fig. , more Cam Ca 2+ was adsorbed into the POPC bilayer than Cal, although the difference in the mixed POPC/POPS was less signiicant. We analyzed the radial distribution function (RDF) of Ca 2+ around the negatively charged groups of the lipid molecules (phosphate and carboxylate oxygen), which showed that the irst density peak is within 3.2 Å (Fig. ). Therefore, we used 3.2 Å as the cutoff to determine the lipid-bound Ca 2+ ions and the coordinated groups around them. We calculated the coordination numbers (Tabel S1) by integrating the radial distribution function of its irst peak. The RDFs of oxygen atoms of lipid molecules around the bound Ca 2+ (Fig. ) showed that, compared to Cal, the Cam model greatly reduced the interaction with PS carboxylate groups, and increased the interaction with the phosphate groups in the mixed POPC/POPS bilayer systems. For example, for the Ca 2+ ions bound to a mixed POPC/POPS bilayer in a 200-mM CaCl 2 solution, there were approximately 0.99 phosphate (0.64 from POPC and 0.35 from POPS) and 0.28 carboxylate (from POPS) groups coordinated with Cam, while there were approximately 0.29 phosphate (0.19 from POPC and 0.10 from POPS) and 1.13 carboxylate (from POPS) groups coordinated with Cal (Table ). This might indicate an advantage of the multisite Cam model over the default Cal model in describing the Ca 2+ -membrane interactions, since previous experimental measurements showed that Ca 2+ may have a stronger binding afinity to phosphate than to carboxylate groups in a POPS bilayer.