Asymmetrical calcium ions induced stress and remodeling in lipid bilayer membranes

Abstract

Ca²⁺ ions play crucial roles in regulating many chemical and biological processes, but their impact on lipid bilayer membranes remains elusive, especially when the impacts on the two leaflets are asymmetrical. Using a recently developed multisite Ca²⁺ model, we performed molecular dynamics simulations to study the impact of Ca²⁺ on the properties of membranes composed of POPC and POPS and observed that both the structure and fluidity of the membranes were significantly affected. In particular, we examined the influence of asymmetrically distributed Ca²⁺ on asymmetric lipid bilayers and found that imbalanced stress in the two leaflets was generated, with the negatively charged leaflet on the Ca²⁺-rich side becoming more condensed, which in turn induced membrane curvature that bent the membrane away from the Ca²⁺-rich side. We employed continuum mechanics to study the large-scale deformations of a spherical vesicle and found that the vesicle can go through vesiculation to form a multi-spherical shape in which a number of spheres are connected with infinitesimal necks, depending on the specific Ca²⁺ distributions. These results provide new insights into the underlying mechanisms of many biological phenomena involving Ca²⁺-membrane interactions and may lead to new methods for manipulating the membrane curvature of vesicles in chemical, biological, and nanosystems.

---

Introduction

Biological membranes are essential for the integrity of cells, providing a barrier between the internal and external environments. They also accommodate a great diversity of proteins, which participate in a wide range of critical biological processes, such as energy or signal transduction, solute transport, cell–cell recognition, and trafficking.¹ The physical properties of membranes are important for their barrier function, which can be regulated by the surrounding environmental factors, such as the presence of various molecules and ions. For example, previous studies have shown that the adsorption of Ca²⁺, one of the most important secondary messengers in cells, can influence many physical properties of membranes, such as the thicknesses,² head group orientations of lipids,^2–5 ordering of acyl chains,^3,5–10 lipid dehydration,^2,4–6,8 and lateral diffusion.^10–16

Phosphatidylcholine (PC) is the most abundant zwitterionic phospholipid in eukaryotic cell membranes and subcellular organelles.¹⁷ Phosphatidylserine (PS) is an anionic phospholipid present in most inner leaflets of cellular organelles, including the plasma membrane, where it accounts for as much as 28% of all lipids.^18–21 Therefore, biological membranes are intrinsically asymmetric with different lipid compositions in the inner and outer leaflets, and a mixture of PC and PS is often used as a general model for membrane studies. Moreover, the concentration of calcium ions across membranes is often distributed asymmetrically. Intracellularly, the free calcium ion concentration typically ranges between 10 and 100 nM. However, within organelles such as the endoplasmic reticulum (ER) and sarcoplasmic reticulum (SR), the concentration can reach significantly higher levels, ranging from 60 to 500 μM, and in some cases, even up to 1 mM.^22,23 In contrast, the extracellular free calcium ion concentration is approximately 1.8 mM.²⁴ In addition, PS lipids were found to play additional roles by interacting with signaling proteins,²⁵ regulating surface charges and protein localization,²⁶ and inducing protein aggregation.^27,28 Some proteins attract PS lipids by nonspecific electrostatic interactions, and the binding can be regulated by Ca²⁺.²⁵ Therefore, the interaction between Ca²⁺ and the negatively charged PS is of particular interest.

It is well documented that the dehydrated head groups of PS and PIP₂, induced by Ca²⁺ adsorption, could lead to separation of lipid components, aggregation of certain types of lipids, phase transitions of membranes, and even collapse of the membrane bilayer structure in experiments.^{13–16,29,30} Notably, there were two contradictory experimental results with respect to the Ca²⁺ induced membrane curvature, one showing that the presence of Ca²⁺ on one side of a membrane can induce a negative curvature^31,32 such that the membrane bends away from Ca²⁺, and the other showing the opposite.³³ To fully understand the experimental observations and the underlying mechanisms, the atomistic interactions between Ca²⁺ and membranes, particularly the impact of asymmetrically distributed Ca²⁺ on asymmetric membranes, need to be further investigated.

Molecular dynamics (MD) simulation is a powerful tool to reveal molecular interactions and has also been extensively utilized to study the interactions between Ca²⁺ and membranes.^{2–5,9–12,34,35} However, previous studies showed that the Ca²⁺ models used in classical nonpolarizable force fields were often inaccurate in describing the Ca²⁺-biomolecule interactions.³⁶ Recently, we developed a multisite Ca²⁺ model, which significantly improved the calculation accuracy for the Ca²⁺–protein interaction.³⁷ Considering that the dielectric constant of membranes is close to proteins, we expect that the new model would generate more reliable simulation results in the study of Ca²⁺-membrane interactions as well.

Another limitation of MD simulation in the membrane study is the use of periodic boundary conditions, which makes it difficult to observe large-scale deformations and remodeling of membranes. To solve this problem, a strategy of combining multiscale methods is helpful. Therefore, in this study, we combined MD simulations and continuum mechanics calculations to investigate the physical properties and lateral pressure profiles of POPC/POPS membranes in the presence of Ca²⁺, as well as the resulting curvature and steady shapes of membranes. Our results showed that the mechanical properties of lipid bilayers, especially PS-rich bilayers, can be significantly influenced by the presence of Ca²⁺, and asymmetrically distributed Ca²⁺ ions lead to asymmetrical lateral stress in the two leaflets of lipid bilayers. With this information adopted, continuum mechanics calculations showed that asymmetric stress can induce a negative curvature on membranes with respect to excessive Ca²⁺ ions, and the resulting steady shape of vesicle membranes, either locally or globally spherical or multi-spherical, depends on the specific Ca²⁺ distribution around the membrane. These results can improve our understanding of Ca²⁺-induced membrane remodeling, which may regulate many biological processes, such as endocytosis, exocytosis, membrane fusion and fission, and perhaps those activities mediated by mechanosensitive membrane proteins as well.^18,38

Results

Adsorption and distribution of ions around lipid bilayer membranes

We performed MD simulations for membrane systems in a solution of CaCl₂. The membranes were composed of either pure POPC or a mixture of POPC : POPS (3 : 1) (Fig. 1(A)), and both the recently developed multisite Ca²⁺ model³⁷ and the default Ca²⁺ model of the CHARMM force field were used in the simulations for comparison, which were denoted Cam and Cal, respectively. The Ca²⁺ concentration was raised stepwise from 0 to 50, 100, 200, 500 and 1,000 mM. The majority of the membrane-bound Ca²⁺ ions bind to the phosphate groups of POPC and POPS, or the carboxylate groups of POPS lipids in simulations (Fig. 1(B)). Then, we analyzed the distribution of various groups along the membrane's normal direction in the equilibrated simulation systems using the SuAVE software.^39,40 The density profiles were similar with various CaCl₂ concentrations, so we only showed the results obtained for two representative cases: 1 M CaCl₂ – a high-concentration regime (Fig. 1(C) and (D)) and 200 mM CaCl₂ – a low-concentration regime (Fig. 1(E) and (F)). Simulations with the two Ca²⁺ 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. 1(C)–(F) pink lines). The density peak of Ca²⁺ roughly overlapped with that of phosphate groups, indicating the adsorption and accumulation of Ca²⁺ under the membrane surface. Notably, more Ca²⁺ was observed around this region in the POPC/POPS system than in the pure POPC system (Fig. 1(C)–(F), yellow lines), which is reasonable since POPS is negatively charged and has a stronger binding affinity to Ca²⁺.

Fig. 1 (A) Single-membrane simulation system. The lipid bilayer is shown in gray. Ca²⁺ and Cl⁻ ions are shown in blue and purple, respectively. Water is shown with a transparent cyan surface. (B) Representative structures of the Ca²⁺–lipid complexes. A Ca²⁺ ion was bound to the PO₄⁻ of POPC, PO₄⁻ of POPS, and COO⁻ of POPS. (C) and (D) Density profiles of selected groups calculated along the bilayer normal direction for pure POPC or mixed POPC/POPS bilayers with 1 M CaCl₂. (E) and (F) Density profiles of selected groups calculated along the bilayer normal direction for pure POPC or mixed POPC/POPS bilayers with 200 mM CaCl₂. The solid lines were obtained using the multisite Ca²⁺ model (Cam) and the dashed lines were obtained using the default Ca²⁺ model (Cal).

There were some differences in the distributions of Ca²⁺ around the membranes when using the two Ca²⁺ models. As shown in Fig. 1(C)–(F), more Cam Ca²⁺ was adsorbed into the POPC bilayer than Cal, although the difference in the mixed POPC/POPS was less significant. We analyzed the radial distribution function (RDF) of Ca²⁺ around the negatively charged groups of the lipid molecules (phosphate and carboxylate oxygen), which showed that the first density peak is within 3.2 Å (Fig. S1, ESI†). Therefore, we used 3.2 Å as the cutoff to determine the lipid-bound Ca²⁺ ions and the coordinated groups around them. We calculated the coordination numbers (Tabel S1, ESI†) by integrating the radial distribution function of its first peak. The RDFs of oxygen atoms of lipid molecules around the bound Ca²⁺ (Fig. S2, ESI†) 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²⁺ ions bound to a mixed POPC/POPS bilayer in a 200-mM CaCl₂ 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 S1, ESI†). This might indicate an advantage of the multisite Cam model over the default Cal model in describing the Ca²⁺-membrane interactions, since previous experimental measurements showed that Ca²⁺ may have a stronger binding affinity to phosphate than to carboxylate groups in a POPS bilayer.^2,41,42

Calcium ions alter the structure and fluidity of lipid bilayers

Area per lipid (APL) is a fundamental parameter describing the molecular packing of lipid bilayers and can reach stable values on a 100 nanosecond timescale in MD simulations. As shown in Fig. 2, APL gradually decreased in the process of adding more CaCl₂ to the system (Fig. 2(A) and (B)), a phenomenon that is more significant in the mixed POPC/POPS membrane. When the Ca²⁺ concentration exceeded 200 mM, the change in APL became less pronounced. It appears that the effect of Ca²⁺ on the packing of membranes containing negatively charged PS is more significant. With saturated CaCl₂ (200 mM), the APL was reduced by approximately 4.2% for the pure POPC bilayer (from ∼0.72 nm² to ∼0.69 nm²) and 23.8% (from ∼0.88 nm² to ∼0.67 nm²) for the mixed POPC/POPS bilayer. The reduction in APL was accompanied by an increase in bilayer thickness (Fig. 2(C) and (D)), which was also more significant in the mixed POPC/POPS bilayer. There was no significant difference in the APL and membrane thickness when using two different Ca²⁺ models, indicating that the observed impact of Ca²⁺ on the overall membrane size was rather robust and independent of the specific Ca²⁺ model used for simulations.

Fig. 2 Structural properties of the simulated lipid bilayers. Area per lipid (APL) for pure POPC (A) and mixed POPC/POPS (B) lipid bilayers with increasing CaCl₂ concentrations. The APL and bilayer thickness were calculated by the SuAVE software^39,40(C) and (D).

We also examined the structural properties of the lipid tails by analyzing the acyl chain C–H bond order parameters (eqn (1) and Fig. 3). Evidently, the lipid tails in the mixed POPC/POPS bilayer systems became more ordered with increasing concentrations of Ca²⁺ ions, which is consistent with previous findings,^5,9,10 and the influence became less significant when the concentration of Ca²⁺ exceeded 200 mM. In contrast, the changes in the order parameters in the pure POPC bilayer systems were minor in all the simulations. The two Ca²⁺ models showed consistent effects on the acyl chain order parameters as well.

Fig. 3 The acyl chain order parameter of POPC in the pure POPC membrane systems (A) and mixed POPC/POPS membrane systems (B), and the acyl chain order parameter of POPS in the mixed membrane systems (C). The solid lines were obtained using the multisite Ca²⁺ model (Cam) and the dashed lines were obtained using the default Ca²⁺ model (Cal).

To further evaluate the performance of the Ca²⁺ models used in this study, we compared our simulation results with experimental data. The C–H bond order parameters obtained from NMR experiments serve as valuable metrics for validating the consistency of lipid structures in MD simulations with experimental data. These measures are particularly relevant as they are associated with key bilayer properties such as area per lipid, thickness, and the conformational variability of individual lipids.⁴³ Previous MD studies used these metrics to evaluate the Ca²⁺ model with electronic continuum correction (ECC).^34,35 Here, we also analyzed the changes of the headgroup C–H order parameters in our simulations with increasing CaCl₂ concentrations and compared with experimental values.^44,45 As can be seen, the change tendencies of the β, α, g₃ order parameters were consistent with the experimental results when both the Cam and Cal models were used (Fig. S3 and S4, ESI†). However, it is also noticeable that most of the simulation results with the Cam model showed better consistency with the experimental data than the Cal model, except for the β order parameter in the pure POPC bilayer (Fig. S3B, ESI†). This indicates that our Cam model overall performed better than the default Cal model in simulating the lipid bilayer systems.

The lateral diffusion coefficient of lipids is closely related to the structure and phase of membranes, which is one of the most important dynamic parameters of biomembranes. From our atomistic MD simulations, we calculated the lateral diffusion coefficients of the lipid molecules in our model membranes from the mean square displacement of lipid phosphorus atoms (Fig. S5, ESI†). The diffusion constant was calculated by fitting a straight line from 100 ns to 150 ns according to eqn (2). The values were approximately 0.4 × 10⁻⁷ cm² s⁻¹ to 1.4 × 10⁻⁷ cm² s⁻¹ depending on the Ca²⁺ concentration (Fig. 4). As a comparison, previous experimental measurements from fluorescence recovery after photobleaching,⁴⁶ pulsed NMR,⁴⁷ single particle tracking,⁴⁸ and fluorescence correlation spectroscopy⁴⁹ yielded the similar magnitude of the diffusion coefficients of approximately 0.5 × 10⁻⁷ cm² s⁻¹ to 1.0 × 10⁻⁷ cm² s⁻¹. Interestingly, the presence of Ca²⁺ significantly reduced the mobility of lipid molecules in the mixed POPC/POPS bilayers, but not to the same extent in the pure POPC bilayers (Fig. 4). Therefore, our results showed that membranes containing negatively charged lipids can become significantly less diffusive in the presence of Ca²⁺. Again, the two Ca²⁺ models did not show significant differences in affecting the lipid diffusion coefficient.

Fig. 4 The lateral diffusion coefficients of lipid molecules in the pure POPC bilayers (A) and mixed POPC/POPS bilayers (B), with increasing CaCl₂ concentrations.

Asymmetrical calcium ions induced asymmetric pressure profiles across lipid bilayers

To further characterize the effect of Ca²⁺ ions on membrane properties, we computed the lateral pressure profiles (LPPs) across the lipid bilayers. Considering that PS lipids are mostly located in the inner leaflet of biological membranes, we set up a double-membrane and asymmetric simulation system as shown in Fig. 5(A), 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²⁺ ions were added to the two water compartments to generate an asymmetrical Ca²⁺ 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²⁺ distribution on biomembrane properties. For instance, the sarcoplasmic reticulum (SR) membrane faces a higher inner Ca²⁺ concentration than the outer, which resembles the upper lipid bilayer in Fig. 5(A), while the outer leaflet of plasma membranes may face a higher Ca²⁺ concentration than the inner leaflet that enriches POPS, which resembles the lower lipid bilayer in Fig. 5(A). 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. 5(B), black curve). This was achieved by adjusting the density difference of water across the membrane such that the osmotic pressure caused by Ca²⁺ density gradient is just balanced off by the opposite water density gradient. The same strategy was also used in the previous study.⁵⁰

Fig. 5 The double-membrane simulation system allowed the calculation of the asymmetrical Ca²⁺-induced lateral pressure profiles and spontaneous curvature. (A) Two asymmetric lipid bilayer membranes created two solvent compartments to impose asymmetrical ionic concentrations. PC lipids were shown in gray, while PS lipids were shown in pink. (B) Lateral pressure profiles calculated from the simulation system with 1 M Ca²⁺. (C) The average differential stress between the leaflets (left y-axis), and the spontaneous curvature of the bilayers c¹₀ and c²₀ (right y-axis). The error bars represented the standard deviation of Δσ or c₀ obtained from five independent MD trajectories. (D) The sign convention of the membrane curvature.

The LPPs of the asymmetrical Ca²⁺ treated bilayer lipids are shown in Fig. 5(B). 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 leaflet, 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.^51–53 In addition, LPP had a positive peak at the bilayer center (IV), which was related to the geometries of the specific lipid types.^54,55

We calculated the differential stress for each leaflet by integrating the lateral pressure π(z) along z from the bulk water to the bilayer center (Fig. 5(B)). The differential stress of the two leaflets was then obtained and denoted as Δσ₁ and Δσ₂ (Fig. 5(A)). We analyzed the lateral pressure profiles π(z) in simulation systems without Ca²⁺ to show the effect of asymmetric lipid compositions (Fig. S6, ESI†). The LPP was also aligned to the bilayer centers to give an easier visual comparison across bilayers (Fig. S6C and S7, ESI†). The Δσ₁ and Δσ₂ 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 significant differential stress between the two leaflets. Therefore, the differential stress shown in Fig. 5(C) was mostly caused by the addition Ca²⁺. As shown in Fig. 5(C), the addition of Ca²⁺ had a minor effect on Δσ₁ (black dotted line), but the impact was significant for Δσ₂ (black solid line), which was increased up to −6.37 mN m⁻¹ while increasing Ca²⁺ concentration on the POPS-rich side. This indicates that the addition of Ca²⁺ to the negatively charged leaflet of membranes will lead to a significant change in the differential stress in this leaflet, and consequently, a stress imbalance will be generated in the two leaflets 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²⁺ on membrane shapes. The results showed no significant difference in the average order parameter across the simulation systems with various Ca²⁺ concentrations, and the values were consistent with those calculated for planar bilayers (Fig. S8, ESI†). This indicates that no significant curvature was generated by the presence of Ca²⁺ in the MD simulations, probably owing to the influence 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.

Asymmetrical stress induces membrane curvature and remodeling

As a result of the asymmetrical stress within the lipid bilayer induced by the Ca²⁺ ions, the membrane tends to bend to release the stress. This tendency for the membrane to bend can be characterized by the spontaneous curvature c₀ in the classical Helfrich model, which has been widely used to model the deformation of a membrane on a large scale.^56–59 Herein, we studied the deformation of a spherical vesicle with a diameter of 100 nm with the Helfrich model, by employing the MD simulation results as the input physical parameters. Basically, we estimated the spontaneous curvature induced by the asymmetrically distributed Ca²⁺ ions viaeqn (6) based on the integral of the first moment of LPPs obtained from MD simulation analyses, which were then input to the Helfrich model to calculate the steady shape of the vesicle (see continuum mechanics calculation in the section of Materials and methods for a detailed description).

Note that any asymmetries between the two leaflets could lead to spontaneous curvature of the membrane. Three types of asymmetries have been studied: (1) compositional asymmetry—one of the leaflets has a different lipid composition or different ratios of lipid compositions than the other; (2) number asymmetry—the two leaflets have different numbers of lipids; (3) environmental asymmetry —for example, the headgroups of the lipids in the two leaflets face different Ca²⁺ concentrations. The spontaneous curvature caused by differential stress as a result of number asymmetry (2) has been well characterized in ref. 50 and 60. In this paper, we studied both the compositional asymmetry (1) and environmental asymmetry (3), but neglected the number asymmetry by intentionally equalizing the number of lipid molecules in the two leaflets. As shown in Fig. 5(C), in the absence of Ca²⁺ concentration, compositional asymmetry alone could induce a spontaneous curvature c¹₀ ≃ 0.04 nm⁻¹ for the lower bilayer and c²₀ ≃ −0.06 nm⁻¹ for the upper bilayer. Here we chose the sign convention that a negative spontaneous curvature tended to bend away from the middle water layer while a positive spontaneous curvatured tended to bend towards the middle water layer. The addition of Ca²⁺ ions significantly enhanced the upper bilayer's tendency to bend away from the middle water layer, and could reach up to c²₀ = − 0.3 nm⁻¹ at high centration of Ca²⁺ ions. While for the lower bilayer, the presence of Ca²⁺ ions was able to reverse its bending direction, from towards the middle water layer to away from it. Combining the two results, we inferred that Ca²⁺ ions near a membrane tended to bend the membrane away from the Ca²⁺ ions. To verify that our observation was not caused by the specific double-membrane simulation setup, we calculated the curvature of another double membrane system, where the lower membrane is upside down, treated with 200 mM Ca²⁺ (Fig. S9, ESI†). We found that Ca²⁺ could still induce a negature curvature for the lower membrane (c¹₀ ≃ −0.16 nm⁻¹), which was qualitatively aligned with those of the original simulation system (Fig. 5(C)). Therefore, although the curvature values may vary depending on the specific simulation setup, the bending directions of the membranes are always consistent.

In order to calculate how a spherical vesicle of a diameter of 100 nm would deform in response to the change in spontaneous curvature c₀, we needed to specify how the surface area A of the vesicle and the volume V enclosed by the vesicle varied when the spontaneous curvature c₀ was changed. Here for simplicity, we considered two types of deformations: (i) in long timescale, the surface area A was assumed to be conserved, while the volume V was allowed to change to keep a fixed osmotic pressure p. (ii) In short timescale, the volume V was assumed to be conserved, while the surface area A was allowed to change to keep a fixed membrane tension σ. Our calculation results showed that the stable shapes under the two conditions are similar, therefore, we presented the results for type (i) deformation in the main text, while leaving the results for type (ii) deformation in the ESI† (Fig. S10 and S11).

We first studied the case in which the area impacted by the Ca²⁺ ions was localized to a small patch (1% of the total surface area) of the vesicle near the north pole. For a vesicle of 100 nm in diameter, as the spontaneous curvature c₀ was gradually increased from 0 to 0.3 nm⁻¹, which matched the maximum spontaneous curvature estimated from eqn (6) in Fig. 5(C), the impacted membrane patch protruded further outward and developed into a hemispherical cap (Fig. 6). The positive spontaneous curvature could be induced by a local accumulation of Ca²⁺ in the inner part of the membrane. Similar trends were also observed for negative c₀, in which the impacted membrane patch bent inward and developed into a small spherical vesicle with a diameter of ∼20 nm (Fig. 6(B)), as indicated by the relative height, which is the difference of the pole-to-pole distance between the deformed vesicle and the undeformed spherical vesicle (Fig. 6(C)). If the diameter of the undeformed spherical vesicle was increased to 200 nm, the impacted membrane patch (1% of the total surface area) would form a pearl-shaped structure in which multi-spherical vesicles are interconnected via narrow necks (Fig. S12, ESI†).

Fig. 6 Membrane deformation induced by asymmetric Ca²⁺ concentration between the exterior and interior of a spherical vesicle. The impacted area was localized to a neighborhood of the north pole of the vesicle. The total surface area of the vesicle was conserved and the osmotic pressure was fixed when varying the spontaneous curvature. (A) Membrane shapes in the minimal energy state for negative spontaneous curvature c₀. (B) Membrane shapes in the minimal energy state for positive spontaneous curvature c₀. (C) The relative height of the deformed vesicle as a function of c₀. The dots in (C) correspond to shapes presented in (A) and (B) of the same color.

When the Ca²⁺ ions were uniformly distributed in the exterior or interior of the spherical vesicle, i.e., the spontaneous curvature of the entire vesicle was impacted by the asymmetric Ca²⁺ ions, the shape transformation showed dramatically different trends between negative and positive values of c₀. For negative spontaneous curvature c₀ < 0 (Ca²⁺ outside the vesicle), the spherical shape remained the most stable one (data not shown). However, for positive spontaneous curvature c₀ > 0 (Ca²⁺ inside the vesicle), we found branches of solutions such that the vesicle exhibited a multi-spherical shape in which several smaller vesicles are connected with narrow necks (Fig. 7(A)). There existed critical values of , where n is an integer, at which the neck is infinitesimally narrow (Fig. 7(B) and (C), indicated by the vertical lines) and the large vesicle would break into a number of n smaller vesicles, so-called vesiculation. The physics behind the vesiculation phenomenon was investigated in reference.61 For a spherical vesicle of R₀, when the spontaneous curvature was increased to c^crit₀, the vesicle would assume a shape of n connected spheres each with a radius of such that the total surface area of the n small vesicles equaled to the surface area of the original vesicle. The neck would assume a saddle shape whose mean curvature just matched c^crit₀. In this configuration, the bending energy of the membrane was almost zero and, therefore was energetically more favorable.

Fig. 7 Membrane deformation induced by asymmetric Ca²⁺ concentration between the exterior and interior of a spherical vesicle. The impacted area covered the entire membrane. The total surface area of the vesicle was conserved and the osmotic pressure was fixed when varying the spontaneous curvature. (A) Membrane shapes in the minimum energy state with increasing positive spontaneous curvature c₀. (B) The height of a deformed vesicle as a function of the spontaneous curvature c₀. Branches of singlet (black), duplet (blue), triplet (red), and tetraplet (green) are shown respectively. (C) The free energy of the vesicles in different branches. For each branch, the lower energy part was shown in solid lines and the higher energy part was shown in dash-dotted lines. The vertical lines in (B) and (C) indicate the critical spontaneous curvature , where n = 2, 3, 4, respectively. The dots in (B) and (C) correspond to shapes presented in (A) of the same color.

Note that our calculation of the membrane shape here assumes the cell is able to secret or expel other osmolytes to keep the osmotic pressure at a constant level via the osmoregulation process in the presence of Ca²⁺ gradient. Therefore, the major impact of the Ca²⁺ gradient on the membrane is to alter the spontaneous curvature of the membrane, while the osmotic pressure is kept constant. In fact, many organisms are able to keep the osmotic pressure at a relatively stable level.⁶² For instance, kidney cells play an important role in osmoregulation in humans by actively controlling the amount of hormones.

Discussion

The electrostatic properties of the membrane interface are of paramount importance in numerous situations. Ca²⁺ is particularly interesting in this respect due to its divalent nature. Although sharing common properties with other divalent cations, Ca²⁺ is one of the most critical secondary messengers and functions by inducing fine-tuned electrostatic effects. For example, it was shown that blood clotting proteins associated with membrane surfaces become activated during the coagulation cascade, and Ca²⁺ enhances protein penetration into the lipid monolayer and induces flocculation of surface lipid aggregates.⁶³ The high Ca²⁺ levels inside presynaptic vesicles can prevent even the most lipophilic transmitters, such as dopamine from adhering to the inner leaflet surface, thus rendering unhindered neurotransmitter release feasible.⁶⁴ Ca²⁺ also regulates the interactions between annexins and membranes, which may influence cell membrane-related domain organization, endocytic, exocytic transport, and trafficking.^65,66 When Ca²⁺ ions were pumped into cells by ion channels or transporters, they can generate large clusters by linking PIP₂.²⁹ In contrast, monovalent cations such as Na⁺ and K⁺ form smaller and less stable clusters. Moreover, a synergy exists between Ca²⁺ and K⁺, facilitating the growth of PIP₂ clusters and the formation of more complex structures within the cell.³⁰ Therefore, there have been extensive studies seeking to elucidate the effects of Ca²⁺ on membrane properties.

As early as 2006, Pedersen et al. performed MD simulations to study the interactions between cations and negatively charged lipid bilayers,⁹ and found that Ca²⁺ ions exert a pronounced effect on bilayers by binding to anionic groups of lipid molecules. Later, further studies showed that the binding of Ca²⁺ made the membrane bilayers more condensed, with decreasing APL, increasing thickness, and dehydration of the membrane interface. Notably, the specific lipid composition and chemistry can also influence these physical–chemical properties of membranes.^2,4,10,67,68 Overall, the binding affinity between Ca²⁺ ions and lipids was significantly stronger than that of Na⁺,⁴ K⁺ and Mg²⁺ ions, regardless of the composition of the lipid bilayers.⁵ Meanwhile, we speculate that divalent cations such as Zn²⁺ would produce similar results to Ca²⁺, while tri-valent ions may induce even stronger effects. Overall, our simulations obtained results consistent with the previous findings. The binding of Ca²⁺ to anionic groups of lipid molecules leads to a reduction in membrane surface area and an increase in membrane thickness. Notably, the effect is more significant for bilayers containing negatively charged lipid molecules. As a consequence, the acyl chains become more ordered, which may facilitate phase transitions.²⁹ Our results also showed that these findings are rather robust, irrespective of the specific ion model and parameters adopted in the MD simulations.

There are still some details that are dependent on the ion model and parameters adopted in the MD simulations. With the new multisite ion model, Ca²⁺ ions were found to be stably bound to the anionic groups of lipid molecules, with similar affinity to the phosphate groups and carboxylate groups (Table S1, ESI†). This is different from the previous simulation results that revealed a strong binding preference of Ca²⁺ to carboxylate groups,^2,4,9,11,12 as also observed in our simulations with the default Ca²⁺ model. Notably, previous infrared absorption spectroscopy experiments suggested that Ca²⁺ ions probably do not stably bind to the carboxylate groups in PS, but can bind to phosphate groups to induce dehydration.^41,42 There were also experiments showing that Ca²⁺ can bind to both carboxylate and phosphate groups of PS.^2,3 In this aspect, our new multisite Ca²⁺ model may have the advantage of reducing the binding bias to carboxylate groups. Nonetheless, although these local differences resulting from specific ion models are appreciable, they do not generate qualitative differences in the global properties of the bilayers discussed in this work.

The effect of Ca²⁺ ions on the dynamic properties of lipid bilayers has not been systematically studied before. Our analysis of the MD trajectories showed that the lipid molecules become less diffusive upon Ca²⁺ binding, for both the zwitterionic and negatively charged lipid bilayers. However, the effect is significant for bilayers containing negatively charged lipids, but only minor for zwitterionic lipid bilayers. We attribute this to the electrostatic “glue” effect of Ca²⁺ on the negatively charged lipids, which makes those lipids more sticky and less diffusive. It appears that when the concentration of Ca²⁺ reached a certain threshold (200 mM in our simulations), the effect became saturated, probably because the binding between Ca²⁺ and lipids was saturated above the threshold and therefore the bilayer did not adsorb more Ca²⁺ to become stickier.

It should be noted that the Ca²⁺ distribution can often be asymmetrical on the two sides of membranes. For example, the plasma membrane often faces a much higher Ca²⁺ concentration on the outer leaflet than on the inner leaflet, while the SR membrane usually experiences a higher Ca²⁺ concentration on the inner leaflet than on the outer leaflet. In fact, membranes themselves are often asymmetric. For example, the inner leaflet of plasma membranes contains much more negatively charged lipid molecules than the outer leaflet. Therefore, studying the effect of asymmetrical Ca²⁺ on asymmetric membranes is of particular interest. In line with previous findings, our atomistic simulations showed that the Ca²⁺-exposed leaflet of membranes becomes more packed. One would expect such a scenario to generate imbalanced stress in the two leaflets of membranes. Indeed, our pressure profiles analysis showed that the addition of Ca²⁺ to the negatively charged leaflet leads to a significant imbalance of the lateral pressure profiles across the membrane Fig. 5(C), while the effect is much less pronounced in the case of the Ca²⁺ ions added to the leaflet with only zwitterionic lipids. As a consequence, this imbalanced lateral pressure distribution can induce spontaneous curvature of the membrane, as indicated by our continuum mechanics calculation and other recent studies.^31,32,69 There were two experimental results regarding the effect of Ca²⁺ on the membrane's bending: one suggests a negative curvature,³²i.e., bending away from the Ca²⁺, while the other suggests the opposite.³³ Our simulation results are in line with Graber et al.'s experimental observation that Ca²⁺ can induce a negative curvature of membranes.³² We suspect that the observation of the positive curvature was due to the different experimental setup,³³ which may introduce other factors that were not accounted for in simulations.

Interestingly, our continuum mechanics calculations showed that, when the impact of Ca²⁺ is focused on a local surface area of a vesicle, a pearl chain extending toward the other side of the membrane can be generated (Fig. 6 and Fig. S10, S12, ESI†). This is similar to recent observations that vesicle chain-growth from the GUV membrane can be maintained as long as there is a supply of membrane material and Ca²⁺, analogous to native secretory vesicles.^31,32 This kind of remodeling can only be generated by an asymmetric effect in the two leaflets of membranes, which may be related to membrane budding and fusion,^70,71 morphology of ER exits,⁷² as well as the generation of interconnected vesicles.⁷³ For example, dynamic membrane remodeling may be achieved by the interplay between lipids and Ca²⁺ to initiate the formation of tubules.^74–76 In contrast, when the impact of Ca²⁺ is uniform on the membrane surface and Ca²⁺ ions are added to the outside, the vesicle tends to remain spherical. When low-concentration Ca²⁺ ions are uniformly added to the inside, vesicles tend to remain spherical as well. Intriguingly, the vesicle will become multi-spherical while more Ca²⁺ ions are added to the inside (larger spontaneous curvature). This suggests that a significant reduction in stress within the inner leaflet of membranes may facilitate vesicle fission, and could potentially aid cell division as well. A significant increase in stress within the outer leaflet of membranes should have a similar effect, and so does the combination of the two. It should be noted though, we did not explicitly consider the asymmetry of the lipid composition of membranes in our continuum mechanics calculations. Moreover, the actual membrane compositions and environments are much more complex than the ideal cases studied here, so the eventual shape may be the result of many factors, including the distribution of proteins.⁷⁷ Nonetheless, the distribution of Ca²⁺ and negatively charged lipid molecules and their interplay may play important roles here.

Our MD simulations and continuum mechanics calculations provided additional atomistic details of Ca²⁺-membrane interactions, and further clarified the question about the controversial membrane curvature caused by Ca²⁺. Although the Ca²⁺ concentration used here was much greater than the average physiological Ca²⁺ concentration measured experimentally, the electrolytic environment around the biomembranes was by no means a static average distribution of ions. Significant localized concentration spikes of Ca²⁺ ions occur along Ca²⁺ signaling pathways, which are always challenging to measure experimentally and can span in excess of three orders of magnitude, and the Ca²⁺ concentration near membranes are known to be higher than that in bulk.^78,79 Therefore, we believe that our results are not only quantitatively informative for chemical systems but also qualitatively meaningful for biological systems, which would be helpful for the understanding of many membrane stress/curvature-related processes.

A potential issue in our MD simulation was that we used the same number of lipids for the asymmetric bilayers to remove the effects of lipid numbers, which may lead to a somewhat frustrated state of the bilayer in the rectangular simulation box due to the size difference of the asymmetric lipids. To alleviate this concern, we implemented a correction by subtracting the lateral pressure of the system in the absence of calcium ions from the lateral pressure of each simulation scenario. This adjustment can effectively neutralize the effects of the potential frustration induced by the asymmetric lipid composition.

Considering the intrinsically asymmetric nature of plasma membranes, with mostly zwitterionic lipids on the outer leaflet and abundant negatively charged lipids on the inner leaflet, it may be reasonable to think that evolution has made the outer surface of cells more inert to external cation distribution in order to maintain a stable shape, while the inner surface of cells is more sensitive to Ca²⁺ for responding to internal signaling. In addition, by influencing the mechanical properties of membranes, Ca²⁺ may also regulate membrane-bound proteins, especially mechanosensitive proteins. For example, mechanosensitive ion channels are involved in many physiological processes, such as touch, hearing, balance, and the sensation of pain,^80,81 which can be activated by forces exerted on membranes and the resulting surface tension, imbalanced stress, and membrane curvature. Therefore, the imbalanced stress across membranes and the resulting membrane curvature induced by Ca²⁺ might provide an additional regulatory mechanism for the activation of mechanosensitive channels. The addition of Ca²⁺ ions to vesicles containing negatively charged lipids may provide a straightforward way to manipulate membrane curvature, which can be useful for regulating many membrane-involved biological, chemical, and nanoscale processes, such as membrane fusion, drug delivery, and even capturing the activated structures of mechanosensitive ion channels. In addition, Ca²⁺ may also have indirect effects on membrane structures, for example through calcium-activated scramblases. When those proteins are activated by Ca²⁺, PS lipids would flip from the inner leaflet to the outer leaflet, which might also change the membrane structures.⁵²

Conclusion

In this work, we combined MD simulations and continuum mechanics to study how Ca²⁺ ions influence the microscopic properties and macroscopic deformation of lipid membranes composed of POPC and POPS. We found that both the structure and fluidity of the lipid membranes were significantly affected by the presence of Ca²⁺ ions. A gradient of Ca²⁺ across the lipid bilayer could create imbalanced stress in the two leaflets, with the negatively charged leaflet on the Ca²⁺-rich side becoming more condensed, thus inducing membrane curvature that bent the membrane away from the Ca²⁺-rich side. For a spherical vesicle, the asymmetric distribution of Ca²⁺ on the two sides of the membrane could deform the vesicle into a locally or a globally multi-spherical shape in which a number of spheres are connected with infinitesimal necks. These results are summarized in Fig. 8, which can help clarify the question about the controversial membrane curvature caused by Ca²⁺ and may lead to new methods for manipulating the membrane curvature of vesicles in chemical, biological, and nanosystems.

Fig. 8 Schematic diagram of the effect of Ca²⁺ (blue spheres) on membrane curvature and remodeling.

Materials and nethods

Simulation systems

The molecular dynamics simulation systems were composed of pure POPC or mixed POPC/POPS in an aqueous CaCl₂ solution. Each single-membrane system was composed of 128 lipids and approximately 4368 water molecules, with varying concentrations of CaCl₂ (Table 1). The bilayers, composed of either POPC or POPC/POPS, were symmetric in these single-membrane systems (Fig. 1(A)).

Table 1 List of the simulation systems

No. of bilayers	Bilayer composition	Lipid amounts	[Ca^2+] (mM)	Duration (ns)	Repeats
Single	POPC	128	[0, 50, 100, 200, 500, 1000]	400	1
	POPC/POPS	96/32		400	1
Double	POPC : POPC/POPS	128 : 96/32		500	5

---

To study the effect of asymmetrical Ca²⁺ on the lateral stress of asymmetric membranes, double-membrane systems were generated by duplicating a single-membrane system along the z direction (Fig. 5(A)). The outer leaflet and inner leaflet of each asymmetric lipid bilayer consisted of POPC (64) and POPC : POPS (48 : 16), respectively. In these double-membrane systems, K⁺ ions were added to neutralize the negative charges of PS lipids. Ca²⁺ ions were added to only one water compartment. Cl⁻ ions were added to neutralize the systems and eliminate the transmembrane potential caused by ion distributions. The details of all the simulation systems are listed in Table 1.

Molecular dynamics simulations

All MD simulations were performed using GROMACS software 2018.6.⁸² The CHARMM36 force field and the CHARMM implementation of the TIP3 water model were used for the simulations.^83–85 In addition to the default ion models (Ca²⁺, K⁺, Cl⁻) of the CHARMM force field, we utilized the recently developed multisite Ca²⁺ model that was optimized to achieve more accurate binding affinity with proteins.³⁷ The temperature was kept constant at 310 K using the Nose–Hoover thermostat.^86,87 The pressure was maintained at 1 bar using the semi-isotropic Parinello–Rahman barostat,^88,89 with a time constant of 1 ps and compressibility of 4.5 × 10⁻⁵ bar⁻¹. A cutoff of 1.2 nm was used for the Lennard-Jones interaction, with a smooth shift function from 1.0 to 1.2 nm. The particle mesh Ewald (PME) method was used for the long-range electrostatics.⁹⁰ After energy minimization, we ran 200-ps NVT and 500-ps NPT simulations to preliminarily equilibrate the systems, and then we carried out 400-ns production simulations for the single-membrane systems and 500-ns production simulations for the double-membrane systems (Table 1).

Analysis of MD trajectories

If not otherwise stated, the last 300 ns trajectories were used for the analysis of the single-membrane systems; while for the double-membrane systems, the last 200 ns trajectories were used for analysis.

The order parameters of lipid tails were analyzed using the following definition:

where θ_z is the angle between the C–H bond and the membrane normal (z direction).

The lateral diffusion coefficient of lipids was calculated with the following equation:

〈r²(t)〉 = 〈(r(t) −r(0))²〉 = 4Dt (2)

where 〈r²(t)〉 is the two-dimensional mean squared displacement (MSD) of phosphorus atoms, D is the diffusion coefficient, and t is the simulation time. A time interval between 100 ns and 150 ns was used for the linear fit to obtain D. The error estimate was done by calculating the difference of the diffusion coefficients obtained from the two halves of the fit interval, as implemented in the “GROMACS msd” tool.

The lateral pressure profiles (LPPs) were calculated using the GROMACS-LS tool.⁵¹ The calculation of the LPPs was based on the last 200 ns of the simulation data (Table 1), with the coordinates and velocities saved every 10 ps. An electrostatic cutoff of 2.0 nm was used for the LPP calculation.⁵¹ The profiles were calculated with about 140 slabs, corresponding to an approximate slab width of 1 Å. The outputs of GROMACS-LS were the components of the stress tensor, σ, as a function of the z coordinate (bilayer normal). LPP was calculated using the formula

π(z) = P_L(z) − P_N (3)

where P_L(z) is the lateral component of the pressure tensor (P_L(z) = 0.5[P_XX(z) + P_YY(z)]), and P_N is the normal component (P_N = P_ZZ).

Continuum mechanics calculations

We first explained how to estimate the spontaneous curvature caused by asymmetric Ca²⁺ concentration from the MD simulations. When two lipid monolayers were combined to form a bilayer, we used the mid-plane between the two monolayers to represent the shape of the bilayer, which was parameterized as X. This position was denoted as z = 0 along the normal direction of the bilayer. Assume the thickness of the two monolayers were d⁺ and d⁻, the differential stress Δσ is defined as

We calculated the spontaneous curvature of each monolayer via the integral^50,60

where κ^± and c₀^± denoted the bending rigidity and the spontaneous curvature of the two monalayers respectively, and π(z) was the LPP defined in (3). In practice, we perform the integration from the midplane of the bilayer to bulk water. The spontaneous curvature of the bilayer then read

For both the bilayer and the monolayer, the first moment of the LLP π(z) gave the product of the bending rigidity and the spontaneous curvature. In order to obtain the spontaneous curvature, the values for the bending rigidity κ and κ^± were needed which in principle could be obtained by a buckling simulation.⁹¹ Here for simplicity we assumed κ^± = 12.5 k_BT and κ = κ⁺ + κ⁻ = 25 k_BT, which were the typical values for the plasma membrane made of POPC.^92–94 Here we neglected the effect that the bilayer bending rigidity κ could be different for the two monolayers due to compositional asymmetry. In fact, the bending rigidity can be influenced by the lipid composition, as well as ions in the external environment. If the membrane is made of a mixture of POPC and POPS (70 : 30), the surface charge might stiffen the membrane and increase the bending rigidity from 25 k_B to 30 k_BT.⁹³ However, the presence of ions can counteract the effect of surface charge and soften the membrane.⁹⁵ As a result of the complicated nature of the bending rigidity, the approximation κ^± = 12.5 k_BT may affect the quantitative accuracy of the calculated spontaneous curvature, but not the sign of it (bending direction of membranes).

We next explained how to calculate the deformation of the membrane as a result of the spontaneous curvature c₀ induced by a local asymmetric Ca⁺ concentration between the inside and outside of the membrane. Assuming that only a small patch of the membrane of area a₀ near the north pole changed its spontaneous curvature to a value of C₀, while the other part of the membrane has a zero spontaneous curvature, c₀ is modeled as a function of the surface area a calculated from the north pole

which exhibited a sharp drop from C₀ to 0 when the area a exceeds a₀. Here β was a parameter to control the sharpness of the drop in the spontaneous curvature. In the case of the entire membrane being influenced by the Ca²⁺ ions, a₀ is set to be a value greater than the total surface area of the membrane. The total free energy of the membrane isHere we neglected the Gaussian bending energy term because we only considered a closed vesicle of spherical topology and the Gaussian bending energy is invariant under deformations of the membrane that retains the topology. The pressure difference between the interior and exterior of the vesicle was assumed to be a constant and denoted by p, and the volume enclosed by the vesicle was denoted by V. For an axisymmetric membrane, the shape is expressed as

X(u,ϕ) = [r(u)cos(ϕ),r(u)sin(ϕ),z(u)] (9)

where the parameter ϕ ∈ [0,2π] denoted the azimuthal angle, and u ∈ [0,1] denoted the position along the meridional direction, with u = 0 corresponding to the north pole and u = 1 to the south pole. The mean curvature of the membrane then readwhere ψ and h satisfied the relation

r′ = h cos ψ, z′ = −h sin ψ. (11)

Here f′ denoted the derivative of function f(u) with respect to u. With these geometrical expressions, the total free energy eqn (8) became a functional of F_tot[r,z,ψ,h,α,β], where α and β were functional Lagrangian multipliers to impose the geometric relations eqn (11). For type (i) deformation, the pressure difference p was given while the membrane tension σ was a Lagrangian multiplier to impose a fixed surface area A. For type (ii) deformation, the membrane tension σ was given while the pressure difference p was Lagrangian multiplier to impose a fixed 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.⁹⁶

Author contributions

C. S. and R. M. designed and supervised the research. C. L. performed molecular dynamics simulations and analysis. Q. Z. conducted continuum mechanics calculation. K. K. parallelized the GROMACS-LS code for the stress analysis. All the authors contributed to the manuscript writing.

Data availability

The multisite calcium ion model used in this work can be downloaded from https://github.com/ComputBiophys/Multisite_Calcium. The simulation data are available upon reasonable request.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by grants from the National Natural Science Foundation of China (32071251 and 21873006 to C. S., 12474199 to R. M.), and Fundamental Research Funds for Central Universities (20720240144 to R. M.). The MD simulations were performed on the computing platform of the Center for Life Sciences at Peking University and the National Supercomputer Center in Tianjin.

References

1. M. Edidin, Nat. Rev. Mol. Cell Biol., 2003, 4, 414–418.
2. A. Melcrová, S. Pokorna, S. Pullanchery, M. Kohagen, P. Jurkiewicz, M. Hof, P. Jungwirth, P. S. Cremer and L. Cwiklik, Sci. Rep., 2016, 6, 1–12.
3. J. M. Boettcher, R. L. Davis-Harrison, M. C. Clay, A. J. Nieuwkoop, Y. Z. Ohkubo, E. Tajkhorshid, J. H. Morrissey and C. M. Rienstra, Biochemistry, 2011, 50, 2264–2273.
4. P. T. Vernier, M. J. Ziegler and R. Dimova, Langmuir, 2009, 25, 1020–1027.
5. H.-H. G. Tsai, W.-X. Lai, H.-D. Lin, J.-B. Lee, W.-F. Juang and W.-H. Tseng, Biochim. Biophys. Acta, Biomembr., 2012, 1818, 2742–2755.
6. P. Garidel, A. Blume and W. Hübner, Biochim. Biophys. Acta, Biomembr., 2000, 1466, 245–259.
7. H. Akutsu and J. Seelig, Biochemistry, 1981, 20, 7366–7373.
8. H. Binder and O. Zschörnig, Chem. Phys. Lipids, 2002, 115, 39–61.
9. U. R. Pedersen, C. Leidy, P. Westh and G. H. Peters, Biochim. Biophys. Acta, Biomembr., 2006, 1758, 573–582.
10. A. Melcrová, S. Pokorna, M. Vošahlíková, J. Sýkora, P. Svoboda, M. Hof, L. Cwiklik and P. Jurkiewicz, Langmuir, 2019, 35, 11358–11368.
11. M. L. Valentine, A. E. Cardenas, R. Elber and C. R. Baiz, Biophys. J., 2018, 115, 1541–1551.
12. A. Martín-Molina, C. Rodríguez-Beas and J. Faraudo, Biophys. J., 2012, 102, 2095–2103.
13. S. Nir, J. Wilschut, J. Bentz, C. Newton, A. Portis and D. Papahadjopoulos, et al. , J. Membr. Biol., 1981, 59, 115–125.
14. G. Feigenson, Biochemistry, 1986, 25, 5819–5825.
15. D. Leckband, C. Helm and J. Israelachvili, Biochemistry, 1993, 32, 1127–1140.
16. T. Ito and S.-I. Ohnishi, Biochim. Biophys. Acta, Biomembr., 1974, 352, 29–37.
17. T. A. Lagace and N. D. Ridgway, Biochim. Biophys. Acta, Mol. Cell Res., 2013, 1833, 2499–2510.
18. L. Li, X. Shi, X. Guo, H. Li and C. Xu, Trends Biochem. Sci., 2014, 39, 130–140.
19. S. E. Gordesky and G. Marinetti, Biochem. Biophys. Res. Commun., 1973, 50, 1027–1031.
20. G. D. Fairn, N. L. Schieber, N. Ariotti, S. Murphy, L. Kuerschner, R. I. Webb, S. Grinstein and R. G. Parton, J. Cell Biol., 2011, 194, 257–275.
21. J. Lorent, K. Levental, L. Ganesan, G. Rivera-Longsworth, E. Sezgin, M. Doktorova, E. Lyman and I. Levental, Nat. Chem. Biol., 2020, 16, 644–652.
22. M. J. Phillips and G. K. Voeltz, Nat. Rev. Mol. Cell Biol., 2016, 17, 69–82.
23. A. P. Ziman, C. W. Ward, G. G. Rodney, W. J. Lederer and R. J. Bloch, Biophys. J., 2010, 99, 2705–2714.
24. A. Bogdanova, A. Makhro, J. Wang, P. Lipp and L. Kaestner, Int. J. Mol. Sci., 2013, 14, 9848–9872.
25. P. A. Leventis and S. Grinstein, Annu. Rev. Biophys., 2010, 39, 407–427.
26. T. Yeung, G. E. Gilbert, J. Shi, J. Silvius, A. Kapus and S. Grinstein, Science, 2008, 319, 210–213.
27. H. Zhao, E. K. Tuominen and P. K. Kinnunen, Biochemistry, 2004, 43, 10302–10307.
28. G. P. Gorbenko and P. K. Kinnunen, Chem. Phys. Lipids, 2006, 141, 72–82.
29. E. Bilkova, R. Pleskot, S. Rissanen, S. Sun, A. Czogalla, L. Cwiklik, T. Rog, I. Vattulainen, P. S. Cremer and P. Jungwirth, et al. , J. Am. Chem. Soc., 2017, 139, 4019–4024.
30. K. Han, S. H. Kim, R. M. Venable and R. W. Pastor, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2202647119.
31. O. Tarun, H. I. Okur, P. Rangamani and S. Roke, Commun. Chem., 2020, 3, 1–8.
32. Z. Graber, Z. Shi and T. Baumgart, Phys. Chem. Chem. Phys., 2017, 19, 15285–15295.
33. M. Simunovic, K. Y. C. Lee and P. Bassereau, Soft Matter, 2015, 11, 5030–5036.
34. J. Melcr, H. Martinez-Seara, R. Nencini, J. Kolafa, P. Jungwirth and O. S. Ollila, J. Phys. Chem. B, 2018, 122, 4546–4557.
35. J. Melcr, T. M. Ferreira, P. Jungwirth and O. S. Ollila, J. Chem. Theory Comput., 2019, 16, 738–748.
36. H. Li, V. Ngo, M. C. Da Silva, D. R. Salahub, K. Callahan, B. Roux and S. Y. Noskov, J. Phys. Chem. B, 2015, 119, 9401–9416.
37. A. Zhang, H. Yu, C. Liu and C. Song, Nat. Commun., 2020, 11, 1–10.
38. M. J. Berridge, Biochim. Biophys. Acta, Mol. Cell Res., 2009, 1793, 933–940.
39. D. E. Santos, F. J. Pontes, R. D. Lins, K. Coutinho and T. A. Soares, J. Chem. Inf. Model., 2019, 60, 473–484.
40. D. E. Santos, K. Coutinho and T. A. Soares, J. Chem. Inf. Model., 2022, 62, 4690–4701.
41. R. Dluhy, D. G. Cameron, H. H. Mantsch and R. Mendelsohn, Biochemistry, 1983, 22, 6318–6325.
42. W. Hübner and A. Blume, Chem. Phys. Lipids, 1998, 96, 99–123.
43. O. S. Ollila and G. Pabst, Biochim. Biophys. Acta, Biomembr., 2016, 1858, 2512–2528.
44. C. Altenbach and J. Seelig, Biochemistry, 1984, 23, 3913–3920.
45. M. Roux and M. Bloom, Biochemistry, 1990, 29, 7077–7089.
46. L. Song, E. Hennink, I. Young and H. Tanke, Biophys. J., 1995, 68, 2588–2600.
47. A. L. Kuo and C. G. Wade, Biochemistry, 1979, 18, 2300.
48. T. Fujiwara, K. Ritchie, H. Murakoshi and J. A. Kusumi, J. Cell Biol., 2002, 157, 1071–1081.
49. P. Schwille, J. Korlach and W. W. Webb, Cytometry, 2015, 36, 176–182.
50. B. Różycki and R. Lipowsky, J. Chem. Phys., 2015, 142, 02B601_1.
51. J. M. Vanegas, A. Torres-Sánchez and M. Arroyo, J. Chem. Theory Comput., 2014, 10, 691–702.
52. M. Doktorova, J. L. Symons and I. Levental, Nat. Chem. Biol., 2020, 16, 1321–1330.
53. M. Doktorova, M. V. LeVine, G. Khelashvili and H. Weinstein, Biophys. J., 2019, 116, 487–502.
54. K. J. Boyd, N. N. Alder and E. R. May, Biophys. J., 2018, 114, 2116–2127.
55. C.-P. Chng, Y. Sadovsky, K. J. Hsia and C. Huang, Extreme Mech. Lett., 2021, 43, 101174.
56. W. Helfrich, Z. Naturforsch., C: J. Biosci., 1973, 28, 693–703.
57. U. Seifert, K. Berndl and R. Lipowsky, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 1182.
58. U. Seifert, Adv. Phys., 1997, 46, 13–137.
59. A. Agrawal and D. J. Steigmann, Biomech. Model. Mechanobiol., 2009, 8, 371–379.
60. A. Hossein and M. Deserno, Biophys. J., 2020, 118, 624–642.
61. L. Miao, B. Fourcade, M. Rao, M. Wortis and R. K. Zia, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 43, 6843.
62. E. Solomon, L. Berg and D. Martin, Biology, Brooks/Cole Publishing, Belmont, CA, 6th edn, 2002.
63. J. Craft, K. Epps, S. Brancato and A. Serfis, J. Colloid Interface Sci., 2003, 268, 181–187.
64. S. Mokkila, P. A. Postila, S. Rissanen, H. Juhola, I. Vattulainen and T. Rog, ACS Chem. Neurosci., 2017, 8, 1242–1250.
65. V. Gerke, C. E. Creutz and S. E. Moss, Nat. Rev. Mol. Cell Biol., 2005, 6, 449–461.
66. Y.-H. Ma, B. Li, J. Yang, X. Han, Z. Chen and X. Lu, Chem. Commun., 2020, 56, 1653–1656.
67. D. E. Santos, A. De Nicola, V. F. Dos Santos, G. Milano and T. A. Soares, J. Phys. Chem. B, 2023, 127, 6694–6702.
68. A. Nascimento, F. J. Pontes, R. D. Lins and T. A. Soares, Chem. Commun., 2014, 50, 231–233.
69. C. Allolio and D. Harries, ACS Nano, 2021, 15, 12880–12887.
70. C. G. Sinn, M. Antonietti and R. Dimova, Colloids Surf., A, 2006, 282, 410–419.
71. P. Drücker, M. Pejic, H.-J. Galla and V. Gerke, J. Biol. Chem., 2013, 288, 24764–24776.
72. A. V. Weigel, C.-L. Chang, G. Shtengel, C. S. Xu, D. P. Hoffman, M. Freeman, N. Iyer, J. Aaron, S. Khuon and J. Bogovic, et al. , Cell, 2021, 184, 2412–2429.
73. B. Ali Doosti, D. Fjällborg, K. Kustanovich, A. Jesorka, A.-S. Cans and T. Lobovkina, Sci. Rep., 2020, 10, 1–7.
74. V. Kirejev, B. Ali Doosti, M. Shaali, G. D. Jeffries and T. Lobovkina, Small, 2018, 14, 1703541.
75. S. Jones, A. Huynh, Y. Gao and Y. Yu, Mater. Chem. Front., 2018, 2, 603–608.
76. A. Doosti, W. Pezeshkian, D. S. Bruhn, J. H. Ipsen, H. Khandelia, G. D. Jeffries and T. Lobovkina, Langmuir, 2017, 33, 11010–11017.
77. C. Kluge, M. Pöhnl and R. A. Böckmann, Biophys. J., 2022, 121, 671–683.
78. N. Xu, M. Francis, D. L. Cioffi and T. Stevens, Studies on the resolution of subcellular free calcium concentrations: a technological advance. Focus on “Detection of differentially regulated subsarcolemmal calcium signals activated by vasoactive agonists in rat pulmonary artery smooth muscle cells”, 2014.
79. M. J. Berridge and R. F. Irvine, Nature, 1984, 312, 315–321.
80. S. S. Ranade, R. Syeda and A. Patapoutian, Neuron, 2015, 87, 1162–1179.
81. S. Sukharev and F. Sachs, J. Cell Sci., 2012, 125, 3075–3083.
82. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1, 19–25.
83. J. B. Klauda, R. M. Venable, J. A. Freites, J. W. OConnor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. MacKerell Jr and R. W. Pastor, J. Phys. Chem. B, 2010, 114, 7830–7843.
84. A. D. MacKerell Jr, D. Bashford, M. Bellott, R. L. Dunbrack Jr, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo and S. Ha, et al. , J. Phys. Chem. B, 1998, 102, 3586–3616.
85. S. R. Durell, B. R. Brooks and A. Ben-Naim, J. Phys. Chem., 1994, 98, 2198–2202.
86. S. Nosé, Mol. Phys., 1984, 52, 255–268.
87. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695.
88. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
89. S. Nosé and M. Klein, Mol. Phys., 1983, 50, 1055–1076.
90. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
91. M. Hu, P. Diggins and M. Deserno, J. Chem. Phys., 2013, 138, 214110.
92. M. Deserno, Chem. Phys. Lipids, 2015, 185, 11–45.
93. M. Doktorova, D. Harries and G. Khelashvili, Phys. Chem. Chem. Phys., 2017, 19, 16806–16818.
94. M. A. S. Karal, M. M. Billah, M. Ahmed and M. K. Ahamed, Soft Matter, 2023, 19, 8285–8304.
95. H. A. Faizi, S. L. Frey, J. Steinkühler, R. Dimova and P. M. Vlahovska, Soft Matter, 2019, 15, 6006–6013.
96. R. Ma and J. Berro, Biophys. J., 2021, 120, 1625–1640.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp01715c

---