Here f ′ denoted the derivative of function f (u) with respect to u. With these geometrical expressions, the total free energy Eq. ( ) became a functional of F tot [r, z, ψ, h, α, β], where α and β were functional Lagrangian multipliers to impose the geometric relations Eq. (11). For type (i) deformation, the pressure difference p was given while the membrane tension σ was a Lagrangian multiplier to impose a ixed surface area A. For type (ii) deformation, the membrane tension σ was given while the pressure difference p was Lagrangian multiplier to impose a ixed volume V . By performing variations against these functions, we obtained a set of ordinary differential equations that can be solved with the bvp4c solver in MATLAB. The detailed derivations of these equations were not carried out in this paper, but could be found in our previous work. 6 Author Contributions Changes of the order parameters for the β, α, and g 3 carbon segments in pure POPC bilayer systems were shown in (B), (C), and (D), respectively. The changes of the headgroup order parameters were obtained by subtracting the order parameters of the system without Ca 2+ from those with Ca 2+ (∆S C H ). We calculated the time-averaged order parameter for each lipid, and then calculated the mean value as well as the standard error of the mean for all the lipid molecules in the simulation system. The simulation results using the Cal and Cam models at 313 K, as well as the experimental data, were shown for comparison. The red lines or dots were calculated for the membrane with POPC/POPS facing Ca 2+ (the upper membrane in Fig. ). The grey lines or dots were calculated for the membrane with POPC facing Ca 2+ (the lower membrane in Fig. ).