Charged Biological Membranes Repel Large Neutral Molecules by Surface Dielectrophoresis and Counterion Pressure

Macromolecular crowding is the usual condition of cells. The implications of the crowded cellular environment for protein stability and folding, protein–protein interactions, and intracellular transport drive a growing interest in quantifying the effects of crowding. While the properties of crowded solutions have been extensively studied, less attention has been paid to the interaction of crowders with the cellular boundaries, i.e., membranes. However, membranes are key components of cells and most subcellular organelles, playing a central role in regulating protein channel and receptor functions by recruiting and binding charged and neutral solutes. While membrane interactions with charged solutes are dominated by electrostatic forces, here we show that significant charge-induced forces also exist between membranes and neutral solutes. Using neutron reflectometry measurements and molecular dynamics simulations of poly(ethylene glycol) (PEG) polymers of different molecular weights near charged and neutral membranes, we demonstrate the roles of surface dielectrophoresis and counterion pressure in repelling PEG from charged membrane surfaces. The resulting depletion zone is expected to have consequences for drug design and delivery, the activity of proteins near membrane surfaces, and the transport of small molecules along the membrane surface.


Charged membrane induced dielectrophoretic force on a neutral spherical object immersed in an ionic solution (derivation of the expression for FD).
The classical expression for the DEP force acting on a neutral, homogeneous, spherical particle of radius a is 1 : where p and w are the absolute permittivities of the particle and the surrounding medium, respectively.E is the modulus of the electric field (or the rms in oscillating fields).The fraction enclosed in brackets is the Clausius-Mossotti factor, related to the effective polarizability of the particle.A first correction to this expression involves solving exactly the electrostatic potential both inside and around the particle.The general solution of Laplace equation for the electric potential within the particle under azimuthal symmetry can be expressed in terms of Legendre polynomials (using spherical coordinates r and  with symmetry around the z axis): ( ) ( ) The external potential is obtained after solving Poisson-Boltzmann equation in its linear approximation in the region outside the spherical particle.The general solution can be expressed as: where In and Kn are modified Bessel functions.For the case of an infinite charged plane with surface charge density  in a solution of Debye length  -1 , the external potential at distance z is given by: and, after expansion in series, it reads: This expression accurately captures the external potential with only a few terms (typically up to n = 3) in the series expansion.Once the electric potential is mapped, we calculate the force on the boundary of the particle by integrating the Maxwell stress tensor, as done in other biomolecular applications as 2,3 where Eo and Ei are the electric field outside (o) and inside (i) the particle boundary surface and n(r) is the unit normal vector pointing away from that boundary.This expression should be further integrated to obtain the DEP force modulus as a series expansion where ā   a is the dimensionless particle radius and leads to the following expression for FD, the z-component of the DEP force, normal to the charged plane: where Kn are modified Bessel functions and the coefficients An are functions of  and  given by ( ) The series in Eq. 7 converges rapidly and a few terms (typically up to n = 5) are enough to get the desired accuracy.The dielectrophoretic free energy of the particle as a function of the distance d to the charged interface can be expressed as the work needed to bring the particle from infinity to z = d.By grouping into prefactor KD all terms non-dependent on z we get: ( ) where KD includes the dependence on the other parameters (particle radius, ionic strength, surface charge density and permittivities).

Pressure exerted by the Electric Double Layer (derivation of the expression for FH)
In addition to dielectrophoretic force, the particles in solution are subject to hydrodynamic forces.
While the former are generally obtained from the Maxwell tensor (or Maxwell stress tensor), hydrodynamic forces are analogously obtained from the integration of the so-called hydrodynamic stress tensor (or pressure tensor): where  is the viscosity of the fluid, u(r) is the solvent velocity field and p(r) is the local pressure.
For steady-state conditions, mass conservation leads to zero divergence of u(r) and Eq. 10 simplifies to: The pressure p(r) follows from solving the Navier-Stokes equation, which in our case simplifies to the form: where (r) is the local charge density due to the mobile ions in solution.By using the Poisson Boltzmann equation for (r) in Eq. 12 the local pressure depends exclusively on the electrostatic field through the relationship 4 : The expression on the right, valid for low potentials, allows using of the expression for the external potential outside the particle (Eq. 3) and, after integrating over the external particle surface, one obtains: for the z component of the hydrodynamic force, where n represents the radial part of the potential outside the particle: ( ) ( ) The coefficients   are given by Eq. 8 and the coefficients   are obtained from the solution of the linear Poisson-Boltzmann equation.
Alternatively, it is possible to use the expression of the (simpler) potential inside the particle, since both coincide on the surface of the particle (r = a), thanks to the boundary condition on the potential.The corresponding hydrostatic free energy of the particle as a function of the distance d to the charged interface can be expressed as the work needed to bring the particle from infinity to z = d.By grouping into prefactor KH all terms non-dependent on z we get: where KH includes the dependence on the other parameters (particle radius, ionic strength, surface charge density and permittivities).

Relative strength of the two contributions (FD and FH).
When using the analytical expressions of the dielectrophoretic and hydrodynamic force acting on PEG molecules in ionic solution near a charged membrane, a question arises whether the dielectrophoretic and the counterion pressure contributions to the overall repulsive force are similar or one of them prevails over the other.
Figure S1 shows the ratio FH /FD at the distance to the membrane where total energy becomes 1 kBT, for ionic strengths ranging from dilute (10 mM) up to concentrated (1 M) solutions and PEG size typical of low-medium MW (roughly 600-6000 Da) near a charged planar surface with  = 0.34 C/m 2 (this is the typical surface charge density of a negatively charged membrane with serine polar groups).At physiological concentrations (100 mM), FH is several times larger than FD.

MD simulations
System preparation.PEG molecule models of various molecular weights (600, 1000, 1540, 2000, 3400 and 6000) have been prepared using the Charmm-gui 5 service (Polymer builder-Polyethylene oxide) by selecting the corresponding number of monomers with the best approximation to the target molecular weight (respectively 13, 22, 35, 45, 77 and 136 monomers for the above molecular weights).These already solvated PEG models have been downloaded in GROMACS 6 and/or NAMD 7 format to prepare the systems considered in this work.Similarly, lipid membrane models have been prepared with Charmm-gui with compositions of 100% DPPC, DPPS and DPTAP and 8:2 DPPC-DPPS composition.In all cases the resulting membrane includes around 180 lipid molecules distributed over the two bilayers.These model bilayer membranes have been downloaded from the Charmm-gui service both in GROMACS and/or NAMD native format.VMD 8 was used to extract and combine the PEG molecules and membranes obtained in the previous steps into a single system.The VMD's "add solvation box" and "add ion" extensions were then used to add water and K + and Cl -ions to a final ionic concentration of 0.1 M (or as needed) and to achieve an electrically neutral system.The final system box obtained, before minimizing and equilibration in all cases studied, had approximate dimensions of 85×85×265 Å and consisted of approximately 150,000 atoms, including ~42,000 water molecules, with minor variations depending on the specific system considered.The TopoGromacs 9 was used to convert the combined system CHARMM topology to GROMACS from within VMD.The final systems obtained following the above procedure have been relaxed following a series of steps which are detailed below (Obtained with minor variations of the protocol suggested by the Charmm-gui 5 ).In all cases, the 2021 version of GROMACS 6 running on a GPU-CUDA system has been used.
Production.Production was done with not less than 5,000,000 steps with a step size of 2 fs in NPT (10 ns), with PME electrostatics and a cutoff of 1.2 nm and Van Der Waals interaction with a Verlet cut-off of 1.2 nm and a force-switch modifier of 1 nm.All H links have been constrained with LINCS.Nose-Hoover thermostat with coupling constant 1 ps and reference temperature of 298.15 K. Parrinello-Rahman isotropic barostat with coupling constant 5 ps, reference pressure 1 bar and compressibility 4.5•10 -5 bar -1 .

PMF computation: Accelerated Weight Histogram (AWH).
The equilibrated systems as described in the previous procedure were used as input for computing the Potential of Mean Force (PMF) using the AWH method as implemented in GROMACS.We selected as a reaction coordinate the distance z between the lipid membrane and PEG mass center, with z ranging from 3.5 to 7 nm, depending on the PEG size.We defined 8 replicas for the AWH with a per replica simulation time of 70 ns (giving a total simulation of 630 ns per simulated system) and with all the remaining simulation parameters as in the production step.Additional systems were prepared containing 120 PEG molecules and a lipid membrane (DPPC or DPPS) in a simulation box with enough water molecules and ions to simulate a 10% w/w PEG / 0.1 M KCl water solution.These systems were run for sufficiently long simulation time to extract the concentration profile of PEG molecules in the vicinity of the membrane.The PMF was obtained by adjusting the mass fraction profiles to an exponential decay from the membrane surface.These simulations allowed to explore any potential crowding effect into the PMF profile by comparing PMF for a single PEG molecule with that for a 10% PEG solution.
Analysis.The output was analysed using the gmx awh command in GROMACS.The PMF obtained in presence of a neutral DPPC membrane was subtracted from the corresponding to a charged membrane (DPPS or DPTAP) to eliminate any entropic contribution from the final PMF profile.Since the lipid bilayer does not have a well-defined region to ascribe the charge distribution (especially in the case of the charged bilayer: DPPS) it has been considered that it approximately corresponds to the positions of the O2L oxygen of the DPPS/DPPC molecule.The final position considered corresponds to the average position of all O2L oxygens (position averaging) with respect to the center of the bilayer, and along the trajectory considered (time averaging).In the case of DPTAP membranes, the selected atom was the CL from the headgroup.All these analyses have been implemented in a program written in python 10 using the vmd-python module (https://vmd.robinbetz.com/) and several numpy module routines from python that write the results in *.xlsx files (pandas python module) for further analysis and representation.

Figure S1 .
Figure S1.Ratio FH/FD between the hydrostatic and dielectrophoretic forces for neutral dielectric particles (p = 10o) of different size in solutions of varying concentrations near a planar surface with charge density  = 0.34 C/m 2 .The ratio corresponds to the position of the particle where total energy becomes 1 kBT.

Figure S2 .
Figure S2.Scaling of the PEG self energy (PMF) obtained from MD simulations with the hydrodynamic radius of the PEG molecule.Values for Rh are taken from 11 .Labels indicate PEG MW.

Figure S4 .
Figure S4.Composition space analysis of NR data for 20:1 D2O:PEG 600 by weight (100 mM KCl, 10 mM tris, pH 7.4).(A) Volume occupancy profiles of each molecular component.For the bilayer + substrate sum and the PEG 600 profile, shading shows the 68% and 95% confidence intervals, while the dark curve shows the best fit.(B-D) Detail of depletion region near the bilayer surface for three lipid compositions: 1:1 DOPC:DOPS (negatively charged), DOPC (neutral), and 1:1 DOPC:DOTAP (positively charged).Shading represents the 68% and 95% confidence intervals.The vertical solid line represents the median distance from the hydrophobic interface at which the PEG 600 density drops to half its solution density.The dashed and dotted lines represent the 68% and 95% confidence intervals of this value, respectively.

Figure S5 .
Figure S5.Neutron reflectivity data for a 1:1 DOPC:DOPS lipid bilayer supported on a silicon substrate and exposed to 5:1 D2O:PEG 600 by weight (100 mM KCl, 10 mM tris, pH 7.4).(A) Neutron reflectivity scaled to the Fresnel reflectivity, i.e. the reflectivity expected between that solution and silicon in the absence of any interfacial structure.Each curve represents the reflectivity of the surface when the reservoir is filled with the solution composition indicated.Error bars are 68% confidence intervals estimated from Poisson counting statistics.(B) Neutron scattering length density (nSLD) profiles of the interface.Shading shows 68% and 95% confidence intervals.

Figure S6 .
Figure S6.Composition space analysis of NR data for 5:1 D2O:PEG 600 by weight (100 mM KCl, 10 mM tris, pH 7.4).(A) Volume occupancy profiles of each molecular component.For the bilayer + substrate sum and the PEG 600 profile, shading shows the 68% and 95% confidence intervals, while the dark curve shows the best fit.(B-D) Detail of depletion region near the bilayer surface for three lipid compositions: 1:1 DOPC:DOPS (negatively charged), DOPC (neutral), and 1:1 DOPC:DOTAP (positively charged).Shading represents the 68% and 95% confidence intervals.The vertical solid line represents the median distance from the hydrophobic interface at which the PEG 600 density drops to half its solution density.The dashed and dotted lines represent the 68% and 95% confidence intervals of this value, respectively.

Figure S7 .
Figure S7.Neutron reflectivity data for a 1:1 DOPC:DOPS lipid bilayer supported on a silicon substrate and exposed to 10:1 D2O:PEG 600 by weight in low ionic strength conditions (5 mM KCl, 5 mM tris, pH 7.4).(A) Neutron reflectivity scaled to the Fresnel reflectivity, i.e. the reflectivity expected between that solution and silicon in the absence of any interfacial structure.Each curve represents the reflectivity of the surface when the reservoir is filled with the solution composition indicated.Error bars are 68% confidence intervals estimated from Poisson counting statistics.(B) Neutron scattering length density (nSLD) profiles of the interface.Shading shows 68% and 95% confidence intervals.

Figure S8 .
Figure S8.Composition space analysis of NR data for 10:1 D2O:PEG 600 by weight in low ionic strength buffer (5 mM KCl, 5 mM tris, pH 7.4).(A) Volume occupancy profiles of each molecular component.For the bilayer + substrate sum and the PEG 600 profile, shading shows the 68% and 95% confidence intervals, while the dark curve shows the best fit.(B-D) Detail of depletion region near the bilayer surface for three lipid compositions: 1:1 DOPC:DOPS (negatively charged), DOPC (neutral), and 1:1 DOPC:DOTAP (positively charged).Shading represents the 68% and 95% confidence intervals.The vertical solid line represents the median distance from the hydrophobic interface at which the PEG 600 density drops to half its solution density.The dashed and dotted lines represent the 68% and 95% confidence intervals of this value, respectively.

Figure S9 .
Figure S9.Neutron reflectivity data for a SiO2 film on a silicon substrate and exposed to PEG 6000 solutions in low ionic strength conditions (5 mM KCl, 5 mM tris, pH 7.4), taken on the LIQREF reflectometer.(A) Neutron reflectivity (unscaled).Each curve represents the reflectivity of the surface when the reservoir is filled with the solution composition indicated.Error bars are 68% confidence intervals estimated from Poisson counting statistics.Solid lines are optimized solutions to the structural model.(B) Neutron scattering length density (nSLD) profiles of the interface.Shading shows 68% and 95% confidence intervals.(C) Volume occupancy profiles of each molecular component.Shading shows the 68% and 95% confidence intervals, while the dark curve shows the best fit.(D) Derived energy of interaction at the SiO2 surface.Shading shows the 68% and 95% confidence intervals; error bars for CANDOR data on an identical preparation represent 68% confidence intervals.

Table S1 .
Analytical prediction of the size of the exclusion zone for PEG molecules of varying molecular weight