1. Introduction

Studies on vesicular membranes composed of bipolar tetraether lipids (BTLs) and diester (or diether) phospholipids are of great biophysical interest. First, the plasma membranes of crenarchaeota, which include thermoacidophilic archaea, contain primarily BTLs and a small amount of diether lipids ^[1]. Archaeal BTLs are able to form vesicular membranes where the BTL molecules adopt an up-right monomolecular configuration spanning the entire membrane ^[1]. Diether lipids, on the other hand, form bilayer membrane structures. At present, it is not clear how the monolayer structures of BTLs and the bilayer structures of diethers are organized in the archaea plasma membrane. Further, the ratio of diethers to tetraethers in crenarchaeota increases with increasing pressure and decreasing temperature ^[2]. This implies that variations in growth conditions can change the amount of bilayer structures relative to monolayer structures in archaeal cell membranes. How such a change in membrane organization affects the functions of the archaea plasma membranes remains elusive. To tackle these biologically relevant yet complex issues, it would be helpful to conduct a biophysical study on well-defined tetraether/diether hybrid liposomes.

Second, BTLs can be used as a membrane stabilizing agent to enhance the usefulness of diester liposomes. For example, tetraether/diester hybrid liposomes have been used as non-viral transfection vectors for gene delivery ^[3,4], nano-particles for oral, skin, and intravenous delivery of peptides or other agents ^[5,6], and a vehicle for photodynamic suppression of local infection ^[7]. It is conceivable that the applications of tetraether/diester (or diether) hybrid liposomes will grow tremendously as researchers gain greater understanding of the role of the monomolecular structures of BTLs in the hybrid liposomes.

Many archaeal BTLs are asymmetric macrocyclic molecules carrying one negative charge at one polar end while the other polar end has no net charge at neutral pH ^[8,9,10,11]. One of such a BTL molecule (GDNT-0) is illustrated in Figure 1a. In the present study, we consider a tetraether/diester (e.g., GDNT-0/dimyristoylphosphatidylcholine (DMPC)) unilamellar vesicle, in which the tetraether lipids (e.g., GDNT-0) adopt an up-right monomolecular structure spanning the membrane ^[12], whereas the zwitterionic diester phospholipids (e.g., DMPC) form a bilayer structure (Figure 1b). Hence, at the neutral pH, the negatively charged phosphate moiety on GDNT-0 could be located close to the outer membrane surface or near the inner surface.

Charged lipid vesicles are usually modelled by using Poisson-Boltzmann equation ^[13]. The approximate solution of the Poisson Boltzmann equation for the electric potential in the outer-vesicular space is ^[14]: $\varphi (R) = {\varphi _0}\frac{{{R_2}}}{R}\exp \left[ {-\frac{{R-{R_2}}}{{{\lambda _D}}}} \right]$ where R₂ is the outer radius of the vesicle, ${\varphi _0}$ is the electric potential at the outer surface of the vesicle and ${\lambda _D}$ is the Debye length in the extra-vesicular ion solution. This approximation is valid at low potentials when $\sinh \left ({\frac{{e\varphi }}{{kT}}} \right) \approx \frac{{e\varphi }}{{kT}}$ where e, k and T are the elementary charge, the Boltzmann constant and the absolute temperature, respectively. Recently an exact analytical solution of the Poisson-Boltzmann equation in cases of spherical, axial, and planar geometry has been obtained in the form of the logarithm of a power series ^[15]. According to the shell theorem ^[16] in the intra-vesicular space the electric potential is constant. In these continuous models the membrane charges both on the extra- and intra-vesicular surface are spread evenly and described by two constants, ${\sigma _{om}}$ and ${\sigma _{im}}$ (called surface charge densities) that are independent of the polar angles, while the location of ions in the intra- and extra-vesicular space are characterized by the continuous ion distribution ^[17,18]. By using these models one can calculate the ion distribution in the extra- and intra-vesicular space and get the following linear relationship between ${\sigma _{om}}$ and the electric potential (${\varphi _0}$) at the outer surface of the vesicle: ${\sigma _{om}} = {\varphi _0}{\varepsilon _0}\varepsilon \left ({1 + {R_2}/{\lambda _D}} \right)/{R_2}$, where ${\varepsilon _0}$ is the vacuum permittivity and $\varepsilon$ the relative permittivity of the extra-vesicular solution. Note, obtaining the above relationships the electric coupling between the outer- and inner-layer of membrane charges was neglected ^[18]. Regarding the proportion of the surface charge densities Israelachvili derived the following relationship ^[19]: $\frac{{{\sigma _{om}}}}{{{\sigma _{im}}}} = \frac{{{R_1}({R_2} + {\lambda _D})}}{{{R_2}({R_1}-{\lambda _D})}}$ where ${R_1}$ is the inner radius of the vesicle membrane.

In this paper we model BTL/diester large unilamellar vesicle, at low BTL mole fractions. In our model of charged unilamellar vesicle we take into consideration all the discrete charge-charge interactions between intra-vesicular ions, extra-vesicular ions, and membrane associated charges. By using Monte Carlo simulation techniques, we determine the distribution of intra- and extra-vesicular ions and membrane associated charges and compare the results of our discrete charge model with that of the above mentioned continuous models.

---

2. Materials and Methods

Here, we develop a simple model of charged lipid membranes containing a mixture of bipolar tetraether lipid (e.g., GDNT-0) and diester phospholipid (e.g., DMPC) molecules (see Figure 1c). Our model considers only one large unilamellar vesicle that is located in the center of a sphere of radius R₃ = 600 nm. The inner and outer radius of the lipid vesicle is R₁ = 100 nm and R₂ = 105 nm, respectively. Note that for better visualization in Figure 1c we use a higher proportion of R₂/R₁. Five groups of charges are distinguished in the model: (group 1) intra-vesicular charges of positive ions, (group 2) intra-vesicular charges of negative ions, (group 3) negative charges of membrane lipids, (group 4) extra-vesicular charges of positive ions, (group 5) extra-vesicular charges of negative ions.

In order to reduce the time of the Monte Carlo simulations we consider the case of low ion concentration in the intra- and extra-vesicular space, i.e., a total of 5000 positive and 5000 negative monovalent ions are located in these spaces. It is assumed that there is no ion transport through the membrane and thus the number of intra vesicular ions (either positive or negative) is constant and it is N₁ = 23 for the positive and N₂ = 23 for the negative ions. While the number of the remaining ions, located in the extra vesicular space, is N₄ = 4977 for the positive and N₅ = 4977 for the negative ions. With these choices of ion numbers the average intra- and extra-vesicular ion concentrations are equal. So far we ran simulations with three different numbers of charged membrane lipid molecules: N₃ = 1000, 5000 and 10000. When N₃ = 10000, the tetraether content in the liposome is estimated to be 2.3 mol %.

---

2.1. Initial configuration

At the beginning of the Monte Carlo simulation, the ions in the intra- and extra-vesicular space are homogeneously distributed. Also the charges of the tetraether lipids are homogeneously distributed on the inner and outer surface of the vesicle membrane, and their number on these surfaces is proportional with the surface area of the inner and outer surface of the vesicle, respectively. Within each group the charges are numbered.

---

2.2. Generating trial configuration

We pick the s-th charge from the j-th group of charges and randomly select a trial position for the selected charge within the region of the respective group. The energy difference between the trial and original configuration, ${E_{trial}}-{E_{orig}}$, is:

${{E}_{trial}}-{{E}_{orig}}=\frac{{{k}_{e}}}{{{\varepsilon }_{w}}}\sum\limits_{k=1}^{5}{\sum\limits_{\begin{align} & i=1 \\ & (i\ne s \\ & at \\ & k=j) \\ \end{align}}^{{{N}_{k}}}{\left[ \frac{{{q}^{j}}{{q}^{k}}}{\bar{r}_{si}^{jk}}-\frac{{{q}^{j}}{{q}^{k}}}{r_{si}^{jk}} \right]}}$

(1)

where the Coulomb's constant ${k_e} = \frac{1}{{4\pi {\varepsilon _o}}} = 8.98755 \times {10^9}Jm{C^{-2}}$, the absolute charge of a charge from the k-th group $|{q^k}| = 1.6 \times {10^{-19}}C$, and the dielectric constant of water ${\varepsilon _w} = 80$, while $\bar r_{si}^{jk}$ and $r_{si}^{jk}$are the trial and original distance, respectively, between the s-th charge picked from the j-th group of charges and the i-th charge of the k-th group. When an intra-vesicular charge is selected the trial spherical angles are unrestricted while the trial radius should be between zero and R₁. When a membrane charge is selected, the trial spherical angles are unrestricted while the trial radius can be either R₁ or R₂. When an extra-vesicular charge is selected the trial spherical angles are unrestricted while the trial radius should be between R₂ and R₃. Generating the homogeneous initial configuration and the trial configuration in the intra- and extra-vesicular space we calculated the spherical coordinates, R, $\varphi$ and θ, from the random number, Ran, as follows:

From

$Ran=\frac{\left( 4\pi {{R}^{3}}/3 \right)-\left( 4\pi R_{i}^{3}/3 \right)}{\left( 4\pi R_{i+1}^{3}/3 \right)-\left( 4\pi R_{i}^{3}/3 \right)}$

(2)

we obtain

$R={{\left[ R_{i}^{3}+Ran\cdot \left( R_{i+1}^{3}-R_{i}^{3} \right) \right]}^{1/3}}.$

(3)

From

$1-2\cdot Ran=\frac{z}{{{R}_{i+1}}}=\cos \varphi$

(4)

where z is one of the x, y, z coordinates of the charge, we obtain

$\phi =\arccos \left( 1-2\cdot Ran \right)$

(5)

Finally, θ can be calculated from

$\theta =2\pi \cdot Ran$

(6)

In Eqs. 2–4, in the case of extra-vesicular charges index i = 2, and in the case of intra-vesicular charges index i = 0 (where R₀ = 0).

Generating the homogeneous initial configuration and the trial configuration for the membrane bound charges (located at either R = R₁ or = R₂) we calculated the spherical coordinate, R, from the random number, Ran, as follows:

$R={{R}_{1}}\ \ \text{if}\ \ Ran<\frac{4\pi {{\left( {{R}_{1}} \right)}^{2}}}{4\pi {{\left( {{R}_{1}} \right)}^{2}}+4\pi {{\left( {{R}_{2}} \right)}^{2}}}$

(7)

otherwise R = R₂, while the $\varphi$ and θ spherical coordinates of the membrane bound charges are generated according to Eq.5 and Eq.6 respectively.

---

2.3. Decision making

The trial configuration is accepted as new original configuration if the following inequality holds:

$Ran<\exp -\left( {{E}_{trial}}-{{E}_{orig}} \right)/kT$

(8)

Otherwise we do not change the configuration. Here $k=1.38\times {{10}^{-23}}J/K$ and $T=300K$.

---

2.4. Monte Carlo cycle

Within a Monte Carlo cycle, every charge in each charge group is picked once, one by one, according to a certain order. At N₃ = 5000 we ran N_MC = 8000 Monte Carlo cycles, while for cases with N₃ = 1000 and N₃ = 10000 we ran N_MC = 10000 Monte Carlo cycles.

---

3. Results

At the end of each Monte Carlo cycle we counted separately the number of positive and negative ions at every shell of the intra- and extra-vesicular space (where the thickness of each shell is 5 nm) and calculated the ion concentration in each shell as follows:

$C_{i}^{+}(K)=\frac{N_{i}^{+}(K)}{\frac{4\pi }{3}\left[ {{\left( {{R}_{i+1}} \right)}^{3}}-{{\left( {{R}_{i}} \right)}^{3}} \right]}$

(9)

where $N_{i}^{+}\left (K \right)$ is the number of positive ions in the i-th shell after the K-th Monte Carlo cycle while R_i and R_i+1 are the lower and upper radius of the i-th shell.

In the case of our model parameters the equilibrium distribution of the charges has been formed by the N_eq = 4000-th Monte Carlo cycle. After the equilibrium is reached we started to calculate the average concentration of positive ions in the i-th shell as follows:

$\bar{C}_{i}^{+}=\frac{1}{{{N}_{MC}}-{{N}_{eq}}}\sum\limits_{K={{N}_{eq}}}^{{{N}_{MC}}}{C_{i}^{+}(K)}$

(10)

Similarly, we also calculated the average number of negative membrane bound charges at radius R₁ as follows:

$\bar{N}_{{{R}_{1}}}^{-}=\frac{1}{{{N}_{MC}}-{{N}_{eq}}}\sum\limits_{K={{N}_{eq}}}^{{{N}_{MC}}}{N_{{{R}_{1}}}^{-}(K)}$

(11)

The standard deviation of the average concentration of the positive ions in the i-th shell, ${{s}_{\overline{C}_{i}^{+}}}$, is calculated from the standard deviation of the average of the number of positive ions in the i-th shell, ${{s}_{\overline{N}_{i}^{+}}}$, by

${{s}_{\overline{C}_{i}^{+}}}=\frac{{{s}_{\overline{N}_{i}^{+}}}}{(4\pi /3)\left\{ {{\left( {{R}_{i+1}} \right)}^{3}}-{{\left( {{R}_{i}} \right)}^{3}} \right\}}$

(12)

Where ${{s}_{\overline{N}_{i}^{+}}}$ is calculated from the standard deviation of the numbers of positive ions in the i-th shell, ${{s}_{N_{i}^{+}}}$ obtained from ${{N}_{MC}}-{{N}_{eq}}$ Monte Carlo cycles as follows

${{s}_{\overline{N}_{i}^{+}}}=\frac{{{s}_{N_{i}^{+}}}}{\sqrt{{{N}_{MC}}-{{N}_{eq}}}}$

(13)

By assuming Poisson distribution the standard deviation (or the width of the distribution) of the numbers of positive ions in the i-th shell is ${{s}_{N_{i}^{+}}}=\sqrt{\overline{N}_{i}^{+}}$ where $\overline{N}_{i}^{+}$ is the average of the number of positive ions in the i-th shell obtained from ${{N}_{MC}}-{{N}_{eq}}$ Monte Carlo cycles.

In Figures 2, 3, 4 and 6, the average concentrations of the positive and negative ions, $\bar{C}_{i}^{+}$ and $\bar{C}_{i}^{-}$ are plotted as a function of the medium radius of the shell in the case of 1000, 5000, 10000, and 2624 membrane bound charges, respectively. The error bar at each point was calculated by Eq. 12.

In Figures 2, 3, 4 and 6 we made weighted fit based on the Levenberg-Marquard algorithm ^[20] of the following double exponential curve to the average extra-vesicular ion concentrations obtained from our Monte Carlo simulations:

${{\overline{C}}^{+/-}}=c\cdot \exp \{\pm \exp [b\cdot ({{R}_{i}}-{{R}_{2}})+a]\}$

(14)

where the values of the a, b and c parameters obtained from the weighted curve fits are listed in Table 1. In the curve fit the weight on each point is the square of the reciprocal of the standard deviation of the respective average concentration (see Eq.12). We choose to fit this double exponential curve to the extra-vesicular average concentrations because the plot of the double logarithm of ${{\overline{C}}^{+}}({{R}_{i}})/{{\overline{C}}^{+}}({{R}_{3}})$ and ${{\overline{C}}^{-}}({{R}_{i}})/{{\overline{C}}^{-}}({{R}_{3}})$ (in Figure 4D and Figure 4B, respectively) resulted in a linear arrangement of the points (see Figures S1C, S2C in the Supplemental Material). On the other hand to the average intra-vesicular concentrations we made weighted fit of a horizontal line.

In Figure 5 the average number of negative charges on the inner- and outer surface of the vesicle membrane are plotted against the total number of membrane charges. The standard deviations of these averages, listed in Table 2, can be calculated similarly to the calculation of the standard deviation of the average number of ions located in a given shell (see Eq. 13).

In Figure 6, the average concentrations of the positive and negative ions, $\bar{C}_{i}^{+}$ and $\bar{C}_{i}^{-}$ are plotted as a function of the medium radius of the shell in the case of 2624 membrane bound charges. In this case all the membrane bound charges are located on the outer surface of the vesicle membrane.

---

4. Discussion

---

4.1. On the extra- and intra-vesicular charge distribution

The simulated average ion concentrations calculated from N_MC-N_eq Monte Carlo cycles severely scatter especially in the case of Figure 2. The inhomogeneous extra-vesicular charge distribution is the consequence of membrane associated charges. Particularly in the case of low amount of membrane associated charges the thermal noise is able to disturb the effect of membrane associated charges on the distribution of intra- and extra-vesicular charges.

The average extra-vesicular positive and negative charge concentration close to the vesicle surface (i.e. close to R₂) is given by $c\cdot \text{exp}\left (exp\left[ a \right] \right)$ and $c\cdot \text{exp}\left (-exp\left[ a \right] \right)$ (shown in Table 1), respectively. The concentration of these positive and negative charges show almost perfect positive and negative linear correlation with the total number of membrane charges, N₃ (see Figure 7). The slopes of the fitted red and blue straight lines are $2.4053\cdot {{10}^{-11}}\pm 2.3958\cdot {{10}^{-13}}\left[ n{{m}^{-3}} \right]$ and $-2.3435\cdot {{10}^{-11}}\pm 2.672\cdot {{10}^{-13}}\left[ n{{m}^{-3}} \right]$, respectively. The absolute values of the slopes are close to each other but the one belonging to the red line is significantly larger. This slight asymmetry may be attributed to the fact that close to the membrane the concentration of the negative ions cannot be less than zero while there is no upper bound for the concentration of the positive ions.

In Table 1 the fitted parameter c estimates the concentration of positive and negative ions at $R=\infty [nm]$. Although these concentrations for positive ions and negative ions are very close to ${{N}_{extra}}/{{V}_{extra}}=5.53\cdot {{10}^{-6}}[n{{m}^{-3}}]$ (the concentration of the homogeneously distributed ions, either positive or negative, in the extra-vesicular space), the deviation increases with increasing the number of membrane associated charges, N₃ and it becomes larger than the error range of parameter c. With increasing N₃, the estimated concentration of the positive and negative ions at $R=\infty [nm]$ is increasingly lower and higher, respectively, than N_extra/V_extra. This is the case because the negatively charged membrane attracts positive and repels negative ions causing higher and lower ion concentrations than N_extra/V_extra, respectively, close to the membrane (and far from the membrane lower and higher ion concentrations than N_extra/V_extra, respectively).

Within the error range the average intra-vesicular ion concentrations obtained from the simulations (and given in the legends of Figures 2, 3, 4 and 6) always agree with the concentration ${{N}_{\operatorname{int}ra}}/{{V}_{\operatorname{int}ra}}=5.491\cdot {{10}^{-6}}[n{{m}^{-3}}]$. Note that close to $R=0[nm]$ the scattering from the average is particularly high. This is the case because during the simulation the number of charges found in the very small sphere, of 5 nm radius, around the center of the vesicle scatters a lot. With increasing distance from the center of the vesicle the sampling volume, i.e. the shell volume, increases quadratically and the average ion concentrations in these shells closely approach ${{N}_{\operatorname{int}ra}}/{{V}_{\operatorname{int}ra}}$.

In summary, according to our simulations, the extra-vesicular concentration of the positive and negative ions decreases and increases double-exponentially with increasing distance from the vesicle surface. On the other hand, the charge distribution is homogeneous in the intra-vesicular space.

---

4.2. On the energy of an extra-vesicular charge

In thermodynamic equilibrium, the charges in the extra-vesicular space follow Boltzmann distribution ^[13]:

${{\overline{C}}^{{}}}({{R}_{i}})\tilde{\ }\exp (-\frac{ze\varphi ({{R}_{i}})}{kT})$

(15)

where e is the elementary charge, z is the charge number, and $\varphi ({{R}_{i}})$ is the electric potential in the i-th shell. Consequently the found double-exponential change of the average concentration (see Eq. 14) means that the average electric energy of an extra-vesicular ion of charge ze changes according to a single-exponential function:

$-\frac{ze}{kT}\varphi ({{R}_{i}})=z\exp \left[ b\cdot ({{R}_{i}}-{{R}_{2}})+a \right]$

(16)

Thus in Figure 8 we calculated the distance dependence of the average electric energy of a positive and negative ion (in thermal energy and mV units) at the extra-vesicular space at different numbers of negative membrane charges.

The rate of the exponential decay of the ion energy is characterized by the Debye length, ${{\lambda }_{D}}={{(-b)}^{-1}}\left[ nm \right]$ (see Table 3).

---

4.3. Proportion of charges on the outer to the inner surface of the vesicle membrane

The number of charges on both the inner (${{N}_{im}}$) and outer surface of the membrane (${{N}_{om}}$) increases linearly with the increase of the number of the total membrane charges (${{N}_{3}}$) (Figure 5), and their ratio is ${{N}_{om}}/{{N}_{im}}=1.104\pm 0.00024$ (this is the weighted average of the ratios given in Table 2. The ratio of the outer to the inner membrane surface is $\frac{4\pi {{({{R}_{2}})}^{2}}}{4\pi {{({{R}_{1}})}^{2}}}=1.1025$ only 1.27% less than 1.104. Since this deviation is larger than the error of ${{N}_{om}}/{{N}_{im}}$ we conclude that the charge density on the outer membrane surface is slightly higher than on the inner membrane surface. What is the reason of this difference? When the intra- and extra-vesicular ion concentration is zero, the charge density at the outer- and inner membrane surface is the same. The situation changes in the presence of intra- and extra-vesicular ions. The negative membrane charges attract and repel the positive and negative extra-vesicular ions, respectively, causing a difference in the average concentration of positive and negative ions close to the vesicle surface (see Figure 7). But the interaction between the membrane charges and extra-vesicular ions is mutual. Since close to the membrane surface the concentration of the positive ions is higher than the concentration of the negative ions the membrane charges are attracted to the outer membrane surface. Thus the charge density at the outer membrane surface becomes slightly higher than at the inner membrane surface. Note that the intra-vesicular charges are homogeneously distributed (see Figures 2-4, 6) and do not affect the ${{N}_{om}}/{{N}_{im}}$ ratio.

---

4.4. Partitioning of charges between the inner and outer surface of the vesicle membrane does not affect significantly the extra-vesicular charge distribution

In the cases of Figures 2-4, the fixed numbers of membrane charges are partitioning between the inner and outer membrane surface in a way that minimizes the system's free energy. Here, in Figure 6, we investigate the effect of a forced partitioning (i.e. ${{N}_{3}}={{N}_{om}}=2624$ and ${{N}_{im}}=0$) on the membrane charges that is rather different from the above mentioned optimal one. Like in the case of optimal partitioning, here we get the proper average ion concentrations far from the membrane surface (i.e. the values of the fitted parameter c for both positive and negative ions in Table 1 are very close to ${{N}_{extra}}/{{V}_{extra}}=5.53\cdot {{10}^{-6}}[n{{m}^{-3}}]$). Also, the average intra-vesicular ion concentrations obtained from the simulations (see legend to Figure 6) agree within the error limits with the concentration ${{N}_{\operatorname{int}ra}}/{{V}_{\operatorname{int}ra}}=5.491\cdot {{10}^{-6}}[n{{m}^{-3}}]$. In this case again the extra-vesicular charge density close to ${{R}_{2}}$ is determined by the total number of membrane charges, ${{N}_{3}}$ and it does not show significant dependence on the way of the partitioning between the outer and inner membrane surface. This is the case because in Figure 7 the red and blue closed circles, referring to the forced partitioning of the membrane charges, are on the red and blue lines, respectively. Finally the value of parameter b is not affected either by the partitioning of the negative charges between the outer and inner membrane surface (see Table 1).

---

4.5. Comparing the results of the continuous models with the discrete charge model

In the discrete charge model, at the considered low ion concentration, the electric potential changes exponentially with the radial distance from the vesicle surface (Eq.16). A similar distance dependence was obtained for the solution of the Poisson-Boltzmann equation in the case of planar charged surface ^[13]. The approximate solution of the Poisson-Boltzmann equation at spherical geometry has been given in the Introduction. In the Supplemental Material this approximate solution of the Poisson-Boltzmann equation was fitted to the result of our Monte Carlo simulations (see green curves in Figure S3). The simple exponential function was also fitted to the simulated data (see black curves in Figure S3). With increasing surface charge the error bars are getting smaller and these two fitted curves get rather close to each other.

According to the discrete charge model the intra-vesicular charges are homogeneously distributed (see Figures 2-4, 6) and do not affect the ${{N}_{om}}/{{N}_{im}}$ ratio. This finding is in accordance with the shell theorem ^[16] mentioned in the Introduction.

In Table 3 there is no significant difference between the listed Debye lengths derived from the results of the Monte Carlo simulations. The weighted average of these Debye lengths is $135.75\pm 5.059\left[ nm \right]$. This finding is consonant with the definition of the Debye length derived from the Debye-Huckel equation for monovalent ions ^[13]:

${{\lambda }_{D}}={{\left( \frac{{{\varepsilon }_{0}}{{\varepsilon }_{w}}kT}{{{({{q}^{+}})}^{2}}[N_{extra}^{+}/{{V}_{extra}}]+{{({{q}^{-}})}^{2}}[N_{extra}^{-}/{{V}_{extra}}]} \right)}^{1/2}}=101.744[nm],i.e.$

its value is close to the ones in Table 3 and it does not change if the number of extra-vesicular positive and negative ions remain the same.

In Section 4.3 we concluded that the charge density on the outer membrane surface is slightly higher than on the inner membrane surface. This is in qualitative agreement with Israelachvili's result ^[19] discussed in the Introduction. However, by substituting the above value of the Debye length into Israelachvili's equation we get: $\frac{{{\sigma }_{om}}}{{{\sigma }_{im}}}=\frac{{{R}_{1}}({{R}_{2}}+{{\lambda }_{D}})}{{{R}_{2}}({{R}_{1}}-{{\lambda }_{D}})}=-112.9$, which is a physically meaningless large negative value. This is the case because the inner radius of the vesicle membrane is smaller than the Debye length. We think that the above equation is incorrect because during the derivation Israelachvili assumed that the intra-vesicular ion distribution is similar to the extra-vesicular one, i.e. it is inhomogeneous.

Finally, in Section 4.2 we found that the electric potential at the outer surface of the vesicle is proportional with the number of charged lipids (see inset to Figure 8). In the Introduction we mention a similar result as a solution of the Poisson-Boltzmann equation, i.e. ${{\sigma }_{om}}={{\phi }_{0}}{{\varepsilon }_{0}}\varepsilon \left (1+{{R}_{2}}/{{\lambda }_{D}} \right)/{{R}_{2}}$. However, in the case of the discrete charge model, the proportionality constant is $9.625\left[ \frac{C}{V{{m}^{2}}} \right]$, while in the case of the continuous model it is only $0.0137\left[ \frac{C}{V{{m}^{2}}} \right]$. Thus the electric potential close to the vesicle surface is much lower in the case of discrete surface charges than in the case of evenly spread surface charges. This is because in the case of the discrete charge model the electric energy close to the vesicle membrane surface is inhomogeneous and the nearby ions show inhomogeneous distribution too, i.e. the ions tend to locate in the energy minima.

In the future we plan to improve the precision of the values of the model parameters (a, b and c) by changing the simulation strategy. Currently there is no limitation for the selection of the trial configuration. For example, if we pick an intra-vesicular charge, its trial position can be anywhere in the intra-vesicular space. Besides the thermal noise, we believe that this is another reason of the strong scattering of the intra- and extra-vesicular charge distribution. In the future, we intend to continue the simulation by switching from unlimited to limited selection of trial configuration, i.e. the trial configuration of the picked charge should be selected in the neighborhood of the original position of the charge.

---

5. Conclusions

Here, we model a negatively charged lipid vesicle by a spherical shell that has zero ion permeability. We take into consideration all the charge-charge interactions between intra-vesicular ions, extra-vesicular ions, and membrane associated charges. Monte Carlo simulations result in homogeneous and double-exponential ion distribution in the intra- and extra-vesicular space, respectively. The extra-vesicular ion concentration close to the membrane surface is proportional to the total amount of the membrane charges and is independent of the partitioning of the membrane charges between the outer and inner membrane surface. If the partitioning of the membrane charges (N_om/N_im) is not restricted (i.e., flip-flop of GDNT-0 is permitted), then N_om/N_im is slightly, but significantly, higher than the ratio of the outer to the inner surface area of the membrane. This model also shows that one should not disregard the effect of the charges on the inner membrane surface when we calculate the ion distributions around a charged vesicle. We compared the results of our discrete charge model and continuous models based on the solutions of the Poisson-Boltzmann equation and pointed out qualitative similarities and sometimes major quantitative differences between these two types of models. Our current simple membrane model considers only charge-charge interactions between the BTL molecules. In the future we intend to develop a more realistic model of BTL/diester membrane where short-range lateral interactions between the membrane components and high BTL concentrations are also taken into consideration. The results obtained from this work may be applied to other negatively charged lipid vesicles such as phosphatidylserine (PS)-containing liposomes (^[21] and references cited therein) and to cancer or other diseased cells or cell-derived micro-vesicles where PS exposure on the outer surface is abundant (^[22] and references cited therein).

---

Conflicts of Interest

There is no conflict of interest regarding this paper.

---

Acknowledgments

PLGC acknowledges the support from NSF (CBET1350841, CBET1437930, DMR1105277). IPS acknowledges Sri Chinmoy's help.

---