Bridges of Calcium Bicarbonate Tightly Couple Dipolar Lipid Membranes

The interaction between lipid membranes and ions is associated with a range of key physiological processes. Most earlier studies have focused on the interaction of lipids with cations, while the specific effects of the anions have been largely overlooked. Owing to dissolved atmospheric carbon dioxide, bicarbonate is an important ubiquitous anion in aqueous media. In this paper, we report on the effect of bicarbonate anions on the interactions between dipolar lipid membranes in the presence of previously adsorbed calcium cations. Using a combination of solution X-ray scattering, osmotic stress, and molecular dynamics simulations, we followed the interactions between 1,2-didodecanoyl-sn-glycero-3-phosphocholine (DLPC) lipid membranes that were dialyzed against CaCl2 solutions in the presence and absence of bicarbonate anions. Calcium cations adsorbed onto DLPC membranes, charge them, and lead to their swelling. In the presence of bicarbonate anions, however, the calcium cations can tightly couple one dipolar DLPC membrane to the other and form a highly condensed and dehydrated lamellar phase with a repeat distance of 3.45 ± 0.02 nm. Similar tight condensation and dehydration has only been observed between charged membranes in the presence of multivalent counterions. Bridging between bilayers by calcium bicarbonate complexes induced this arrangement. Furthermore, in this condensed phase, lipid molecules and adsorbed ions were arranged in a two-dimensional oblique lattice.

B . Figure S1: A. Solution X-ray scattering curves from a dispersion of 20 mg/ml DLPC, whose wide-angle X-ray scattering part (q > 6) is presented in Figure 1B. The solution was dialyzed at 45 • C for two months against a solution of 8.3 mM CaCl 2 and 10 wt% PEG (black curve), and at room temperature for a week against the same solution (blue curve). The series of q > 6 peaks are explained in Figure 1B. The low q data of the blue curve show a single lamellar phase, with a repeat distance of 5.5 ± 0.1 nm. The low q data of the black curve shows the two coexisting phases, the swollen phase with D of 5.7 ± 0.1 nn and the second lamellar phase with D = 3.45 ± 0.02 nm. B. Solution X-ray scattering curves from 20 mg/mL DLPC dispersion, dialyzed at 45 • C against 5 wt% PEG in various salt solutions, as indicated in the figure. The first lamellar phase had a repeat distance, D, of 5.5 ± 0.1 nm as revealed by the peaks indicated by the black arrows. The lamellar repeat distance of the second lamellar phase was D = 3.45 ± 0.02 nm, as revealed by the peaks indicated by the purple arrows. The dialysis was performed for a week (curve A), two weeks (curve B), or three weeks (curves C to G). Experiments D-G were repeated 3 times. C was repeated 5 times. A was repeated 10 times and in single case the tightly condensed phase was not observed. B was repeated 36 times. We note, however, that in 14 cases a dialysis of two weeks was insufficient for observing the tightly condensed phase. In three of those cases an additional week was required for detecting the tightly condensed phase and in the other cases it was not detected, and most likely required longer incubations. S-2 Figure S2: 2D image of X-ray scattering data from a 20 mg/mL DLPC dispersion, dialyzed at 45 • C against a solution of 8.3 mM CaCl 2 and 5 wt% PEG. The first lamellar phase had a repeat distance, D, of 5.5 ± 0.1 nm as revealed by the full scattering rings. The lamellar repeat distance of the second phase was 3.45 ± 0.02 nm, as revealed by scattered spots and scattering arcs.
S-3 Figure S3: 2D image of X-ray scattering data from a 20 mg/mL DLPC dispersion, dialyzed at 45 • C against a solution of 8.3 mM CaCl 2 , 1 mM Na 2 CO 3 , and 5 wt% PEG. The first lamellar phase had a repeat distance, D, of 5.5 ± 0.1 nm as revealed by the full scattering rings. The lamellar repeat distance of the second phase was 3.45 ± 0.02 nm, as revealed by scattered spots and partial scattering arcs. S-4 1 . 2 I n t e n s i t y ( a . u ) q ( n m -1 ) Figure S4: Solution X-ray scattering curves from silver behenate. When the sample was smaller than the cross section of the beam (black curve), the full width at half maximum, σ, of the principle correlation peak was 0.019 nm −1 . When the sample was larger than the cross section of the beam (red curve), σ of the principle correlation peak was 0.022 nm −1 . S-5

Calculation of the Ion Surface Excess
The excess of molar species i was computed using where c 0 is the concentration in the mid-plane between the two bilayers, and z is the coordinate normal to the membrane surface. z GDS was determined as the z coordinate where the water excess Γ H 2 O was zero. The integration is taken from the center of the bilayer at z = 0 to the center of the solution at Z. We present the normalized excess, Γ i /c 0 , to account for the concentration differences along the z direction. The concentration of ions in the center of the bilayer was assumed to be zero. The results were averaged over both membrane interfaces. Figure 3B presents the normalized excess ions near the bilayers on the basis of our simulations. The strong synergistic binding of HCO -3 and Ca 2+ is reflected by their drastically higher normalized excess when simultaneously present in the solution. Note that in Figure 3A we have set the origin of the z-axis to z GDS , i.e. the Gibbs-Luzzati plane.

1
We define an interfacial radial distribution function (RDF), with the aim of sampling interfacial structure. Ordinary RDFs are widely used and can be easily interpreted. To sample the structure within an interface it is practical to restrict the sampling to the particles inside the interface (see Figure S5). In particular, 2D RDFs (in other words, with a radius containing only x and y components) are contaminated with information from the bulk and 1 Taken from the SI of the preprint "Calcium ions promote membrane fusion by forming negative-curvature inducing clusters on specific anionic lipids" Christoph Allolio, Daniel Harries https://www.biorxiv.org/content/10.1101/2020.04.29.068221v1.full S-6 from the opposing leaflet, as are ordinary RDFs of bound ions.
Following Ben-Naim S1 we begin with the generic pair distribution function which describes the probability to find one particle of type A in the section of space between X and X + dX and a particle of type B between X and X + dX . In particular we define the pair-correlation function g AB by Assuming the sytem is fluid, we expect that for long distances g AB (X, X ) → 1, which means the density of particle B at X is uncorrelated to the density of particle A at X. If the fluid (interface) is also homogeneous, then it is possible to set an arbitrary origin at X and then write g AB (R(X, X )), as a function of only the distance R between the particles. Otherwise this is not generally possible, but one loses some information when transforming to ρ (2) In the next step, we assume that our system only consists of a particular region X, X ∈ S in which the interface is located. This region corresponds to an open system, so that the probability is replaced by an average number of particles. In particular, if the bulk is large enough, the interfacial system is in a grand canonical ensemble, with a chemical potential of particles inside the interface equal to that of the corresponding particles in the bulk. We can still define a g AB (X, X ), which is arbitrary outside the region, but inside will be equivalent to the definition for the full system. g AB (X, X ) is still reducible to g AB (R(X, X )), for the points in the region, provided that the matter inside the region is a homogeneous fluid.
This condition is not generally fulfilled for fluid interfaces. Nevertheless, for example in a flat interface along the z direction, translational symmetries exist in the x, y directions and S-7 rotational symmetry exists in the azimuthal angle, φ, so the expression should become exact for small enough intervals dz and dθ, respectively. In particular, in an open bulk system g AB (R) is the same as in the bulk over the entire domain. The number of particles of type B in a spherical volume dR around particle A at X = 0 is which is known as the coordination number. Note, that the integration order can be chosen at will, provided the limits are chosen correctly, however only the innermost integral limit can depend on a function, so that the outermost limit has to go over all possible radii. For the region S, we chose to carry out the integration in R last: The coordination number N CN includes only atoms in the interfacial region S. As g AB (R) depends only on R, regardless of the actual form of the limits, it will be integrated as a constant over these two coordinates. Before the final step this part will therefore accumulate a prefactor dV S (R) dR , with V being the volume of region S. We include the Jacobian contribution R 2 into this volume, because this removes the coordinate system dependency. The final integration will therefore be: The coordination number can be extracted from simulation by simple counting of adjacent particles in slices, up to to the cutoff R . Up to this point the discussion has been completely general to any region, provided it is sufficiently homogenous. As we only have to define V (R ), we chose the geometry in Figure S5, where the interface consists of a flat region, limited by water on the outside (high z) and lipid on the inside (low z). The volume of a S-8 spherical cap is known S2 to be with h as the height of the spherical cap. This volume has to be subtracted from the full volume of the sphere, so that where the subtraction only occurs if r > h or r > h respectively.
with d, d defined as distances along the z axis from the box limits as shown in Figure S5. p z is the z position of the particle and Z is the extension of the box in the z direction. The origin of the z axis was set to the membrane center. This center was computed using the z positions of the membrane lipids as computed from the lipid definitions used for the ReSIS surfaces. S3 In this way sampling was reduced to one half space along the local membrane normal, taking the membrane center as the origin. The interval defining the interface was d, d =1 nm. In other words, the pair distribution function includes roughly only the interfacial region, as there is ≈2 nm of bulk water density, and a membrane thickness of ≈3 nm. We use these cutoffs to compute g AB (R) from the coordination numbers in the simulation. The half space below the membrane was sampled in the same way, so that coordination numbers include sampling over both.
S-9 Figure S5: Schematic illustration of the calculation of the interfacial RDF. At small r, the calculation is identical to normal spatial distribution functions. At large r, parts of the bulk and the lipid interior are not included, restricting the function to the interface.

Water Orientation
We use the water orientation along the membrane normal as a proxy for the electric charge on the membrane, see Figure S7. In the neutral system the distribution is flat in the space between membranes, orientation is exclusively through the membrane dipoles. The center of the water layer is set to z = 0, so that the distribution is flat. The degree of preferential adsorption of one species leads to a field orienting the water molecules, as shown by the slope of the orientation. This field is partially neutralized when adding HCO -3 . For neutral membranes, the water orientation has been used as an order parameter to estimate the decay of hydration repulsion. S4,S5 S-11  Figure S7: Orientation of normalized water dipoles along the z axis.