Curvature effects in charge-regulated lipid bilayers

The effects of electrostatic interactions Markovich et al. (2021); Bossa et al. (2021) on properties of lipid bilayer membranes are significant, ubiquitous and biologically Hong and Brown (2008) and biomedically Ewert et al. (2021) relevant, since naturally occurring as well as technologically relevant lipids usually carry numerous dissociable or polar groups Cevc (1990), and lipid membrane charge density has been identified as a universal parameter that quantifies the transfection efficiency of lamellar cationic lipid–DNA complexes Ewert et al. (2021). The particular relevance of electrostatic interactions in lipidic systems very clearly transpired from the structural principles of the two vaccines approved against COVID-19 in December 2020 and based on the cationic lipid vesicle vector and the anionic mRNA payload developed by Pfizer/BioNTech Mulligan et al. (2020) and Moderna Baden et al. (2021). In general, charged lipids capable of modulating their charge depending on the bathing environmental conditions, have been recognized as a key component of ionizable lipid nanoparticles which are currently acknowledged as one of the most advanced non-viral vectors for the efficient delivery of nucleic acids Schlich et al. (2021).

While electrostatic interactions in bio-soft matter are universal Holm et al. (2001), there is a fundamental difference between the standard colloid electrostatics and membrane electrostatics Borkovec et al. (2001), in the sense that the membrane charge, just as the protein charge Lund and Jönsson (2005); Barroso da Silva and Jönsson (2009); da Silva and Dias (2017), depends on the solution environment of the lipid Avni et al. (2018) and its changes can engender also changes in the shape of the lipids and concomitant structure of lipid assemblies. One example of this solution environment effect would be changes in protonation/deprotonation equilibria of the dissociable phospholipid moieties depending on the solution pH, that in general affect the charge of lipids’ headgroup, and another would be the modification of the strength of electrostatic interactions wrought by the ionic strength of the bathing aqueous solution that affects the ionic screening. In particular, for therapeutic gene delivery Podgornik et al. (2008) ionizable lipids are designed to have positive charges at acidic pH in the production stage to ensure a strong interaction with nucleic acids, after which they should become almost neutral under physiological conditions ($\mathrm{pH}~{}7.4$) to prevent sequestration in the liver, and should finally turn positive again upon being introduced into a typical acidic environment of the endosomes, promoting the delivery of the genetic cargo, followed finally by release of the genetic cargo at neutral $\mathrm{pH}$ enabling the cargo to interact with the cell machinery Schlich et al. (2021).

To connect the dissociation equilibria of dissociable surface groups and electrostatic fields, Ninham and Parsegian (NP) Ninham and Parsegian (1971) introduced the charge regulation (CR) mechanism, that couples the Langmuir isotherm model of the local charge association/dissociation process with the local electrostatic potential Markovich et al. (2021), resulting in an electrostatic self-consistent boundary condition. The NP model can be generalized to include other adsorption isotherms, depending on the detailed nature of the dissociation process Podgornik (2018).

The protonation/deprotonation equilibria at the dissociable surface groups in phospholipid membranes involve local $\mathrm{pH}$ and local bathing solution ion concentrations, which in general differ from the bulk conditions Nap et al. (2014). This implies that the changes in the bathing solution properties will have a pronounced effect on the effective charge of the membrane, thus modifying $\mathrm{pH}$ sensing and $\mathrm{pH}$ response of lipid membranes Angelova et al. (2018). In fact, $\mathrm{pH}$-dependence of the lipid charging state can directly enable an ionizable lipid with a dissociation constant $\mathrm{p}K_{a}\sim 6.5$ to be neutral in the blood circulation, thus preserving a bilayer structure, and then revert to its charged protonated form at an endo-lysosomal $\mathrm{pH}$, consequently reverting to an hexagonal phase upon contact with anionic membrane lipids Jayaraman et al. (2012).

In what follows we will show that even in the case of a chemically symmetric curved membrane, i.e. with the same lipid composition on the inner and outer leaflets, the CR process introduces a charge asymmetry and an out-of-plane membrane polarization vector in a highly non-linear fashion, characterized by a sudden symmetry breaking transition involving in all other respects chemically identical outer and inner leaflet surfaces of the membrane. This phenomenology closely parallels the recently discovered symmetry breaking charge transitions in interactions between pairs Majee et al. (2018) or stacks Majee et al. (2019) of charged membranes as well as the complex coacervation driven by the charge symmetry broken states in solutions of macroions with dissociable surface moieties Majee et al. (2020).

Charged lipids. Different naturally occurring charged lipid species, an important fraction of which is anionic Galassi and Wilke (2021), represent a notable part of biological membranes, and consequently their role in the structural, dynamical and mechanical properties of model lipid bilayers as well as their effect on the membrane proteins is of great biological interest Pöyry and Vattulainen (2016). The less commonly occuring cationic lipids, on the other hand, have played an important role as the main charged component of the cationic liposome vectors for nucleic acid therapeutics Podgornik et al. (2008).

The charge of phospholipid polar heads originates in deprotonated phosphate groups, protonated amine group, and deprotonated carboxylate group, with the corresponding dissociation constants expected to lie in the region of $0\leq\mathrm{p}K_{a}\leq 2$ for the primary phosphate group, in the region $6\leq\mathrm{p}K_{a}\leq 7$ for the secondary phosphate group, in the region $3\leq\mathrm{p}K_{a}\leq 5$ for the carboxylate groups, and $9\leq\mathrm{p}K_{a}\leq 11$ for the primary amine group Marsh (2013). Since in the subsequent analysis the effect of the electrostatic interactions - including the local polarity, the ionic strength and the distribution of neighboring dissociable groups - will be considered explicitly, we consider only the intrinsic $\mathrm{p}K_{a}$, noting that even in the absence of CR certain structural and continuum electrostatic effects are sometimes included in the calculation of the dissociation constants when determined by e.g. PROPKA Olsson et al. (2011).

Of the different phospholipids, phosphatidylcholines have very low $\mathrm{p}K_{a}$ and can be considered as undissociated, but in particular phosphatidylserine, phosphatidylinositol, as well as phosphatidylglycerol, phosphatidylethanolamine, cardiolipin and phosphatidic acid all contain at least one dissociable moiety at relevant pH conditions and should be considered as potentially charged Cevc (2018). Among the cationic lipids DOTAP (1,2-dioleoyl-3-trimethylammonium-propane) remains charged irrespective of pH, while the amine groups of other cationic lipids acquire their charge by protonation only below a certain pH Leventis and Silvius (1990). As an extreme case one could mention the custom-synthesized cationic lipid MVLBG2 with a headgroup that bears no less than 16 positive charges at full protonation Ewert et al. (2021). In ionizable lipids the quaternary ammonium head of cationic lipids can be in addition substituted with a titratable molecular moiety with an engineered dissociation constant that would ensure the charge would be mouldable by the environmental $\mathrm{pH}$ Schlich et al. (2021). Among the ionizable lipids that respond strongly to pH in the relevant range one can list DODAP (1,2-dioleoyl-3-dimethylammonium propane), DLinDMA, DLin-KC2-DMA and finally DLin-MC3-DMA, all of them used for proper lipid vector formation in cytosolic delivery of therapeutic cargo as reviewed recently in Ref. Schlich et al. (2021).

Protonation and pH effects. Detection and response of the cells to either global pH environment or local pH gradients is crucial for optimal cellular metabolism, growth and proliferation Casey et al. (2010). The sensing of pH heterogeneities is important for initiating migration, polarization and deformations of lipid bilayer assemblies Angelova et al. (2018) and can be driven by protonation changes of phospholipids that can furthermore directly affect the membrane curvature Ben-Dov and Korenstein (2012). In fact pH regulation of phospholipid charges have been implicated in pH-induced migration, pH-induced bilayer local or global deformation, as well as pH-induced vesicle chemical polarization Angelova et al. (2018). The fact that the charge state of many phospholipids depends on the solution pH implicates by and large a charge regulation mechanism as indeed first proposed for surfaces containing dissociable groups by Ninham and Parsegian Ninham and Parsegian (1971).

In addition, the phospholipid protonation state and the membrane charge distribution can modify the interaction energy with a protein on approach to the membrane Lund and Jönsson (2013). This effect too is due to perturbed protonation equilibrium at the membrane titratable sites, with phospholipid charge regulation determining the strength and even the sign of the protein-membrane interaction in variable chemical environments Javidpour et al. (2021). Charge regulation can thus be seen as a biological $\mathrm{pH}$-sensitive switch for protein binding to phospholipid membranes and thus conferring signaling functions to different types of lipids Tanguy et al. (2018).

Curvature and lipid dissociation effects. Global and local changes in pH, inducing changes in the charged states of lipids, can affect lipid vesicle shape deformations as well as phase-separated domain formation Angelova et al. (2018). Usually these effects are framed within the electrostatic contribution to the mechano-elastic properties of membranes, such as surface tension and bending rigidity, which have been reviewed in detail, see e.g. Refs. Andelman (1995), Fogden and Ninham (1999) and Galassi and Wilke (2021). Electrostatic coupling of membrane charges leads to a salt concentration dependent increase in the bending rigidity of the charged membrane compared to a charge-free membrane Faizi et al. (2019). Usually the details of lipid dissociation mechanism are not considered in these works and analyses going beyond the fixed charge assumption, incorporating direct lipid charge dissociation models, are less common but have been nevertheless invoked when modelling the flexoelectric properties of membranes Petrov (1999). In fact Langmuir adsorption isotherm was used to estimate the change in the surface density of counterions on the inner and outer leaflet surface of the phospholipid membrane induced by the coupling between curvature and electrostatics Hristova et al. (1991).

In what follows we will present a detailed analysis of the charge regulation effects in a system, see Fig. 1, comprised of a single spherical unilamellar lipid vesicle of a fixed phospholipid membrane thickness $w\simeq 4\,\mathrm{nm}$, dielectric constant $\varepsilon_{p}\simeq 5$ and a variable inner radius $R$, immersed in a simple univalent aqueous electrolyte solution, of dielectric constant $\varepsilon_{w}\simeq 80$. We will use the full Poisson-Boltzmann (PB) theory to evaluate the electrostatic part of the free energy, as well as the often invoked and easier to implement linearized Debye-Hückel (DH) theory in the second order curvature expansion approximation. This will allow us to ascertain which effects are non-linear in nature as far as electrostatics is concerned. We will show that charge regulation of an otherwise symmetric bilayer membrane induces three important phenomena: charge symmetry breaking, non-linear flexoelectricity and anomalous curvature dependence of free energy and that in general the pH effects go beyond the paradigm of electrostatic renormalization of the mechano-elastic properties of membranes.

The outline of the paper is as follows: In Sec. II, we describe the details of the charge regulation model and formalism for this study. Section III is devoted for the results of the study. Finally, the discussion of the results of the curvature effects in charge-regulated lipid bilayers is given in Sec. IV.

For all the numerical evaluations presented herein, we have used typical system parameters such as the Bjerrum length $\ell_{B}=0.74\,\mathrm{nm}$, the membrane thickness $w=4~{}{\rm nm}$, the dielectric constant of water $\varepsilon_{w}=80$ and that of the lipid as $\varepsilon_{p}=5$. For the dimensionless lipid density at the neutral plane, $\tilde{n}_{0}$, defined in Eq. (17), we took a reasonable value $\tilde{n}_{0}=7.65$ which for an aqueous uni-univalent electrolyte solution of $140\,\mathrm{mM}$ salt concentration (inverse Debye length $\kappa_{D}=1.215\,\mathrm{nm}^{-1}$, or equivalently, screening length $\lambda_{D}=0.823\,\mathrm{nm}$), corresponds to $n_{0}=1\,\mathrm{nm}^{-2}$. With fixed $\tilde{n}_{0}$, because of the scaling relations, Eq. (17), changing the screening length implies also changing the lipid density at the neutral plane, $n_{0}$.

In making sense of the numerical results we need to take due cognizance of the fact that our continuous electrostatic formulation remains valid only for radii of curvature much larger then the thickness of the membrane, $R\gg w$. Because all of the distances are scaled by the Debye length this furthermore implies that once the screening length is chosen, the numerical results are consistent only for $h\times(\kappa_{D}w)\ll 1$.

The two interaction parameters $\alpha$ and $\chi$ can be varied within relevant ranges, $-30\leq\alpha\leq 30$ and $0\leq\chi\leq 40$ as argued before. In order to convert the dimensionless paramaters to physical ones, we note from Eq. (7) that $(\mathrm{p}K_{a}-\mathrm{pH})=\alpha/\ln 10$ and $\chi$ is given in thermal units. In addition, we also present the differences between the full non-linear PB theory and the linearized DH theory expanded to the second order in the curvature of the membrane. This comparison is important since the linearized DH theory expanded to second order in curvature is often used to quantify the electrostatic and curvature effects in membranes Safinya and Radler (2021). Clearly both calculations exhibit similar qualitative features, but can be quantitatively quite different. We scan the parameter space in order to identify the important phenomena connected with charge regulation, which is our primary focus here, but do not specifically apply our model to any particular lipid membrane system.

From our previous works Majee et al. (2018, 2019) on two CR surfaces interacting across an electrolyte solution we know that depending upon the values of the parameters $\alpha$ and $\chi$ the Frumkin-Fowler-Guggenheim adsorption isotherm can lead to a charge symmetry breaking transition, corresponding to unequal equilibrium charge of two surfaces. This happens as a result of a competition among the three major interactions present in the system: the adsorption energy of ions or protons onto the charge-regulated surfaces, interaction of neighboring protonated sites on the surface and the electrostatic interaction between two charge-regulated surfaces. While the last one is at play only for a pair of surfaces, the former two are relevant even for a single CR surface and it can be shown that they lead to equally deep minima of the grand potential for charge densitites differing in magnitude as well as in sign for each CR surface in the absence of the other. A pair of interacting surfaces then automatically adopts an asymmetric charge distribution owing to an electrostatic-attraction-mediated reduction of the system energy. The charge asymmetry was found to be highest around the line $\chi=-2\alpha$ and the charge symmetry broken region persisted also in the neighborhood of this line in the earlier works Majee et al. (2018, 2019). The present case, describing a dielectric layer sandwiched between electrolyte layers, differs from the model of Ref. Majee et al. (2018), pertinent to an electrolyte layer sandwiched between dielectric layers in the sense that the region between the two charged phospholipid leaflets is a dielectric, impermeable to electrolyte ions, while the electrolyte solution fills the rest of the space. The salient features of the charging behavior thus display a rather different behavior.

In fact while in our system we still find a symmetry breaking transition in the charging fractions, $\eta_{1,2}$, the symmetry broken charge region and the symmetric charge regions occupy switched places in the $(\alpha,\chi)$ "phase diagram" if compared to the case of Ref. Majee et al. (2018); the line $\chi=-2\alpha$ and its neighborhood thus corresponds to the symmetric charge region, while the rest of the "phase diagram" is symmetry broken. In addition, the extent of the symmetric region also depends on the curvature of the membrane, $h$, in such a way that the larger the curvature the larger is the extent of the charge symmetry broken region. This trend starts already at very small values of curvature. We now elucidate these statements in all the relevant detail.

We note that the parameters relevant for the plots are the dimensionless dissociation energy $\alpha$, the dimensionless lateral pair interaction energy of occupied neighboring sites $\chi$, and the dimensionless curvature $h$, see Table 1 for definitions. As already stated, when converting from dimensionless to physical quantitites, once the screening length is chosen, one can only consider numerical results for $h\times(\kappa_{D}w)\ll 1$.

The two charge fractions $\eta_{1,2}$ as a function of $\alpha$ are displayed for different $\chi$ and a fixed membrane curvature $h$ in Fig. 2. Clearly, for $\chi=20$ even at the small curvature $h=0.02$ (corresponding to $R=\rm 41.15~{}nm$), both the PB and DH results show a pronounced charge asymmetry, which becomes even more prominent for higher values of $h$, see Fig. 3 for details. This charge asymmetry corresponds to the charge symmetry breaking in the bilayer. In addition, the charge fractions are not only different but can exhibit a discontinuous transition as a function of $\alpha$. For $\chi=20$ only the PB results show this transition for $\eta_{2}$, whereas the DH results do not; for $\chi=35$ the PB results show a discontinuity for both $\eta_{1}$ and $\eta_{2}$, while the DH results show a discontinuity only for $\eta_{2}$ in the charge state of the lipids. As the curvature is increased the system remains in a charge symmetry broken state between the outer and inner leaflets of the membrane. We therefore conclude that for a spherical vesicle with finite, even if small, curvature the charge regulation standardly leads to a charge symmetry broken state. As a result the values for the charge fractions of the inner and the outer membrane layer differ even if the two leaflets of the lipid membrane are chemically identical. In addition, depending on the values of the charge regulation model parameters and curvature, the dependence of $\eta_{1,2}$ on either $\alpha$ or $\chi$ can be continuous or discontinuous.

In order to be able to analyze the dependence of the charge densities $\sigma_{1},\sigma_{2}$ of the two membrane leaflets on the charge regulation parameters, we present a phase diagram in Fig. 3 that shows the variation of $|\sigma_{1}-\sigma_{2}|$ as a function of $\alpha\in[-20,20]$ and $\chi\in[0,40]$ for $h=0.1,0.5,1.0,$ and $2.0$. The (white) regions, corresponding to charge symmetric states of the two leaflets, are clearly discerned and their location and extent depends on $\alpha$ and $\chi$, as well as on the the bilayer curvature $h$, which obviously has a profound effect on the charge state of the bilayer, i.e., on the value of $|\sigma_{1}-\sigma_{2}|$. We reiterate that in the previous work Majee et al. (2018, 2019) the symmetry broken region was centered on the line $\chi=-2\alpha$, while in the present case it is the symmetric state which is centered on that line, while the rest of the phase diagram corresponds to a symmetry broken state.

In addition, for $h=0.1$, we actually see not one but three regions of charge symmetric states, corresponding to $\sigma_{1}=\sigma_{2}$ (indicated by white color in Fig. 3, almost coinciding for both non-linear PB and linear DH calculations. From the phase diagram it is difficult to see what the three symmetric branches correspond to, but we will later show that they correspond to the change in sign of the charge density difference, $\sigma_{1}-\sigma_{2}$. Among the regions with charge symmetry, the one centered on $\chi=-2\alpha$ (the middle region) passes through almost the same range of $\alpha$ for both calculations. At $h=0.5$, the line $\chi=-2\alpha$ fully passes through the range $\alpha\in[-20,0]$ in the DH theory, while in the PB theory, the line is terminated at $\alpha\approx-16$. Further increasing the dimensionless curvature to $1.0$ and beyond, the line is terminated at $\alpha=-20$ for both theories. The more pronounced asymmetric states in Fig. 3 of the PB case extend over a broader region than for the DH case at $h=0.1$ and $0.5$. With increasing $h$ up to $1.0$ and higher, the phase diagram shows clearly that only the PB case can exhibit the highest asymmetry (represented by dark blue color, while the DH solution remains broadly less asymmetric. Also the PB theory exhibits a wider range of $|\sigma_{1}-\sigma_{2}|$ values than the DH solution.

We now turn our attention to the details of the curvature dependence of the charge state of the membrane. At the beginning a caveat is in order: smaller values of the curvature are only accessible in the DH approximation, whereas the solution of the full PB theory take unreasonably long time because the system size has to be increases concurrently with the decrease in curvature. It is for this reason that the curvature dependence in the PB theory is truncated at finite values of curvature.

The curvature dependence of the difference between the inner and outer charge density, $\sigma_{1}-\sigma_{2}$, for negative and positive values of $\alpha$ is fully displayed in Figs. 4 and 5, respectively. Obviously this dependence is fundamentally non-monotonic, including the changes of sign. For negative values of $\alpha$, see Fig. 4, the difference $\sigma_{1}-\sigma_{2}$ is a non-monotonic function of the curvature exhibiting regions of unbroken charge symmetry, corresponding to $\sigma_{1}=\sigma_{2}$, as well as regions of broken charge symmetry with $\sigma_{1}\neq\sigma_{2}$, with the difference between the two surface charge densities varying from positive values to negative values on increasing the curvature. In fact this behavior can be seen clearly also in Fig. 3 where it corresponds to a curvature cut through the three lines of charge symmetry for a fixed $(\alpha,\chi)$ combination in the phase diagram. Clearly for a sufficiently negative $\alpha$, e.g. for $(\alpha,\chi)$ combination $(-20,5)$, the difference between the inner and outer charge density, $\sigma_{1}-\sigma_{2}$ ceases to change sign within both the theories, but their non-monotonic behavior is still retained. The two thus represent separate features of the charging state of the curved bilayer. The behavior for the positive values of $\alpha$, see Fig. 5, differs significantly. For a variety of $(\alpha,\chi)$ cuts through the phase diagram we detect here no changes in sign for the difference between the inner and outer charge density, but we do see remarkable non-monotonicity in its dependence on curvature with very clear-cut differences between the PB and DH results. Nevertheless, $\sigma_{1}-\sigma_{2}$ seems to increase for small curvature, reaching a local maximum, then dropping and increasing again but less steeply for larger curvatures. As for the range of curvatures that we display on Fig. 5, one also needs to consider the inherent limitations of the continuum assumptions which form the basis of the present calculations and the curvature should not be extended to arbitrary large values.

The curvature dependence of the difference in surface charge density between the inner and outer leaflet, $\sigma_{1}-\sigma_{2}$, implies the existence of a dipolar moment in the direction of the normal of the bilayer, see Eq. (19), and therefore also the existence of flexoelectricity, but with a variable magnitude and sign of the flexoelectric coefficient. Clearly, the region of a simple proportionality between the bilayer polarization and curvature is limited and depends crucially on the charge regulation parameters.

For $\alpha>0$, Fig. 5, the difference $\sigma_{1}-\sigma_{2}$ does not in general change sign, but remains nevertheless non-monotonic so that the system remains in a broken charge symmetry state for all indicated values of the charge regulation parameters $\alpha$ and $\chi$. The calculations seem to indicate two separate regions of approximately linear flexoelectricity, but with flexoelectric coefficients of different magnitude: one at small curvatures (below $h\simeq 0.5$ for PB calculations and $h\simeq 0.1$ for DH calculations) and another one at larger curvatures (above $h\simeq 0.75$ for PB calculations and $h\simeq 0.2$ for DH calculations), until finally $\sigma_{1}-\sigma_{2}$ varies only weakly with curvature.

For $\alpha\ll 0$, Fig. 4, the situation seems more complicated and also the differences between the PB and DH calculations more evident. There appears a linear flexoelectric regime at very small curvatures that changes sign for larger curvatures (evident for e.g. $\alpha,\chi=-20,5$ or $-10,10$ case for PB and DH calculations, Fig. 4). As $\chi$ is increased the linear flexoelectric regime for small curvatures is eventually cut short for large enough curvatures and the system fully restores its charge symmetry with vanishing flexoelectricity (evident for e.g. $\alpha,\chi=-15,30$ case for PB and DH calculations, Fig. 4). Since the PB calculation cannot probe the regime of very small curvature because of numerical problems we can rely only on the DH results that show a vanishing flexoelectricity also for very small curvatures at not too large $\alpha<0$ and $\chi>0$. The non-linear features of flexoelectricity just described are pertinent to the charge regulation model which takes into account certain salient features of the phospholipid protonation/deprotonation or other charge dissociation reactions at the membrane-electrolyte solution boundary. They do not appear in fixed charge models of membrane electrostatics. Furthermore it is clear that the property of flexoelectricity depends crucially on the membrane dissociation properties as well as the solution conditions, and is thus far from being a universal property of the membrane composition.

We finally examine the hypothesis that the electrostatic effects in membrane can be reduced to a renormalized spontaneous curvature and renormalized bending rigidity, standardly invoked in the context of electrostatic interactions in phospholipid membranes Andelman (1995); Faizi et al. (2019); Shojaei et al. (2016). In Fig. 6 we plot the curvature dependence of the equilibrium electrostatic surface free energy density as obtained from the linearized DH or fully non-linear PB theory. If indeed the electrostatic effects could be reduced solely to renormalized values of the bending rigidity and spontaneous curvature, then the curvature dependence of the equilibrium electrostatic surface free energy density should show a parabolic dependence centered at the spontaneous curvature. What we observe in Fig. 6 corroborates this expectations but only for positive values of $\alpha$. In general, and in particular for negative values of $\alpha$, however, the curvature dependence exhibits a behavior quite different from these expectations. In this case one observes, see Fig. 6, multiple curvature minima and in some cases a vanishing value of the electrostatic free energy as both surfaces become completely neutralized. The paradigm of renormalized spontaneous curvature and renormalized bending rigidity has thus only a limited validity and its range of applicability depends again crucially on the phospholipid dissociation properties and solution conditions. The membrane composition, affecting the two charge regulation parameters, indeed plays a role in the curvature properties but so do the bathing solution conditions that collectively determine the nature and magnitude of electrostatic effects.

There are of course a plethora of works based on molecular and coarse grained simulations of lipid bilayers and their interactions Feller (2000); Wang and Deserno (2010); Bruhn et al. (2016); Joshi and Deshmukh (2021) but none of these studies exhaustively explored the contribution of electrostatic interactions to the charge states of dissociable lipid headgroups Schlich et al. (2021). Even in the most accurate parametrization of the lipid molecular force fields that take into account the dissociation state of the lipid headgroup through the partial charges, the charge regulation, i.e., the dependence of the dissociation state on the local bathing solution conditions, is usually not considered Yu et al. (2021). To include the effects discussed in our paper one would need to implement the charge dissociation reactions explicitly into the simulation scheme along the lines of recent advances in the simulation techniques for reaction ensembles Landsgesell et al. (2020); Curk et al. (2022). Inclusion of the explicit charge regulation model into the detailed mean-field electrostatic theory applied to a curved phospholipid bilayer membrane allowed us to uncover and analyse several effects that have heretofore remained outside the paradigm of electrostatic renormalization of the membrane mechanical properties Faizi et al. (2019).

The first effect depending on the detailed charge regulation mechanism is the charge symmetry breaking, leading to unequal charges of the inner and outer phospholipid membrane surface even if - and this is important - they be chemically identical, i.e. described by the same free energy parameters. This effect was first observed in the case of membranes interacting across an electrolyte solution Majee et al. (2018, 2019), however, the charge symmetry breaking for a curved bilayer differs from this case since the two charged surfaces in a membrane interact across a simple dielectric, and charge symmetry breaking is in some sense inverse to the case of interacting membranes. In the most drastic case the charge symmetry broken state can be characterized by one of the monolayers near neutral and only one charged (see the case corresponding to $\chi=20$ and $\alpha\approx-9$ in Fig. 2, for example).

The second important effect is the existence of non-linear flexoelectricity Deng et al. (2014), with flexoelectricity itself being well known and even explained by a type of charge regulation mechanism for membranes Hristova et al. (1991). While sophisticated electrostatic models have been formulated for flexoelectricity analysis Loubet et al. (2013), more detailed charge regulation description and full mean-field electrostatic theory indicate that the flexoelectric constitutive relation is in general not linear and in addition its form depends crucially on the charge regulation parameters. Only certain intervals of dissociation parameters would correspond coarsely to a linear flexoelectric constitutive relation, with dipolar moment proportional to the curvature Loubet et al. (2013). In general, however, the proportionality is non-linear or there might actually be no flexoelectricity for sufficiently large curvatures.

The last important modification in our analysis of the charge regulation effects is the dependence of the free energy on the curvature. As stated before, the prevailing paradigm is to see the electrostatic effects as commonly renormalizing the elasto-mechanical properties of the membrane, such as surface tension, curvature and bending rigidity (for a recent description see Ref. Faizi et al. (2019)). Our results indicate that this is not the complete story and that the free energy of a charge-regulated phospholipid membrane exhibits a much richer variety of curvature dependence. Only for a limited interval of the phospholipid dissociation parameters does one indeed observe a clear quadratic dependence on curvature that would be consistent with the electrostatic renormalization of the membrane elasto-mechanics.

The bilayer membrane asymmetry need not necessarily involve the obviously asymmetric protein distribution in biological membranes. It can also exhibit differences in lipid composition between the two leaflets Ingólfsson et al. (2014) or be a consequence of different solution conditions across the separating membrane Karimi et al. (2018). Based on our calculations we can state that even without any compositional asymmetry, or any solution asymmetry, the phospholipid dissociation coupled to electrostatic interactions and curvature would itself contribute an essential charge asymmetry to the otherwise completely symmetric membrane properties.

Experimentally, the very same model that we used above to describe the charge regulated membrane electrostatics was used to elucidate the liquid-liquid ($L_{\alpha}\longrightarrow L_{\alpha^{\prime}}$) phase transition observed in osmotic pressure measurements of certain charged amphiphilic membranes Harries et al. (2006). The phenomenon of charge symmetry breaking, that can induce attractions between chemically identical membranes, was argued to induce an attractive disjoining pressure in plant thylakoid membranes and photosynthetic membranes of a family of cyanobacteria Majee et al. (2019). We are thus confident that the charge regulation effects coupled to membrane curvature are real and could be detected in solutions of varying acidity and salt activity. While there is no lack in experiments probing the effects of pH on membrane properties of giant unilamellar vesicle such as a pH change induced vesicle migration and global deformation Kodama et al. (2016), vesicle polarization coupled to phase-separated membrane domains Staneva et al. (2012) and the effect of localized pH heterogeneities on membrane deformations Khalifat et al. (2011), a systematic quantitative comparison between expected pH-induced effects and observed membrane response is lacking. One of the reasons for this is that the pH-curvature coupling is complicated and non-linear, as we argued and demonstrated above. The phenomenology of this coupling presented in this work should help to elucidate the details of the pH effects and the role played by the dissociation mechanism in phospholipid membranes.

While the measured effective rigidity of lipid membranes seem to corroborate the paradigm of electrostatic renormalization of elasto-mechanical properties of membranes Faizi et al. (2019), the experimental work of Angelova and collaborators, quantifying the effects of pH changes on lipid membrane deformations, polarization and migration of lipid bilayer assemblies seem to present a more complicated picture Staneva et al. (2012); Kodama et al. (2016); Angelova et al. (2018). They demonstrated that a global or a local pH change can induce localized deformations of the membrane as well as the whole vesicle, polarization in membranes with phase separated lipid domains as well as whole vesicle migration. The models of pH effects used in this context are usually not based on explicit free energy contributions of charge dissociation, being the starting point of our work, but rather build upon assumed effective values of the membrane charges, and while the experimental studies do show the existence of different instabilities the detailed quantitative connection with our work is at present difficult to establish. Nevertheless, our elucidation of the coupling between the charging processes in lipid membranes and electrostatic interactions in general provide a solid underpinning for these type of phenomena which are clearly beyond the paradigm of electrostatic renormalization of elasto-mechanical properties of membranes.

Finally, our methodology and the model used have understandably and unavoidably many limitations. First of all, the limitations inherent in the PB continuum electrostatics, applying also to our theory, are well known and understood Markovich et al. (2021). The assumption of incompressibily is relatively standard in membrane physics but of course entails certain well recognized limitations Faizi et al. (2019). We only consider a quenched membrane state with fixed density and type of lipids, thus disregarding the annealing of the composition either by in plane diffusion or by trans-membrane flip-flops of different lipid species both associated presumably with large(er) time scales Hossein and Deserno (2020). We do not exactly specify the composition of the membrane, neither its lipid part nor, possibly even more important, the membrane protein counterpart Javidpour et al. (2021), aiming for salient characteristics of the behavior and disregarding the certainly important consequences of molecular identity Nap et al. (2014). We believe that irrespective of all these acknowledged limitations we uncover some features of membrane electrostatics which can be crucial in assessing and controlling the behavior of charged, decorated lipid vesicles.

RP and PK acknowledge funding from the Key Project No. 12034019 of the National Natural Science Foundation of China and the 1000-Talents Program of the Chinese Foreign Experts Bureau and the School of Physics, University of Chinese Academy of Sciences.