Scientific paper
The Application of the Replica Ornstein-Zernike Methodology for Studying Ionic Membrane Equilibria
Barbara Hribar-Lee
Faculty of Chemistry and Chemical Technology, University of Ljubljana, Aškerčeva c. 5,
SI-1000 Ljubljana, Slovenia
* Corresponding author: E-mail: barbara.hribar@fkkt.uni-lj.si
Received: 30-01-2012
Dedicated to Prof. Dr. Gorazd Vesnaver on the occasion of his 70h birthday
Abstract
The applicability of the replica Ornstein-Zernike (ROZ) method, often used to explore the partly quenched systems, for studying ionic equilibria in membrane-like systems was investigated. A simple +1:-1 electrolyte solution was modelled as a primitive model electrolyte on the McMillan-Mayer level of approximation. The particles mimicking membrane were modelled as permeable hard spheres frozen in space in their positions. The ion-membrane potential was chosen to conform as closely as possible to the experimental data. The density profiles of ions near the membrane, as well as adsorption isotherms, were calculated using the ROZ equations complemented by the hypernetted chain (HNC) approximation. The method provides reasonable results in qualitative agreement with experimental observations.
Keywords: Membrane equilibira, primitive model electrolyte, replica Ornstein-Zernike integral equations, hypernet-ted-chain approximation, adsorption isotherms
1. Introduction
One of the essential bio-membrane functions is being a barrier to solute flux between aqueous compartments in the living organisms. The membranes forming the cell walls are usually formed as phospholipid bilayers and they interact with the small solute molecules, ligands, both ions and nonelectrolytes. The small molecules bind to the membrane, they can be either adsorbed at the surface or partitioned into the interior of the bilayer. In the second case, the membrane is characterized as being permeable to the molecules, and the rate of permeability is given by partition coefficients.1 When ions are present in a system, the total electric potential energy profile of the membrane must be considered. The electrical component of the free energy of an ion near or within the membrane is important, not only for the understanding how divalent and monovalent metal ions bind to the membrane surface, but also for estimating the local pH at the membrane surface, and for modelling how some enzymes and ion channels might be regulated by the voltage across the bilayer. Considering the structure of the lipid membrane the internal dipol potential is formed inside the bilayer. The dipo-
les oriented at the membrane surface result in a positive potential in the center of the phosphotidylcholine in the membranes surrounding the biologicl cells. This internal dipol arises due to the orientation of the fatty ester car-bonyls of each phospholipid molecule. The resulting electric potential profile shows that in the case of small ions the permeability barrier is very large for both ionic species and is slightly lower for the anions than for the ca-tions.1 Properties of solution-membrane interface play a large role in medical and pharmacological phenomena,2-5 as well as in technological processes.6
Membrane equilibria can be treated as partly-quenched systems with directional dependent potential.7 Bryk et al.8 used a density functional approach to study phase behaviour of a Lennard-Jones fluid in a system of slit-like pores separated by semipermeable walls. Same approach, in a combination with the Monte Carlo computer simulations, was used by Boda et al. to study selective partitioning of two restricted primitive model electrolytes across the membrane, permeable to only one electrolyte.9 Further, the electrostatic potential near lipid membrane in an ionic solution was calculated using the Poisson-Boltz-mann (PB) approach,10 the usage of this mean-field ap-proache for studying the electrostatic properties of the
membranes is in detailes discussed by Andelman.11 Braca-montes et al.12 used ROZ integral equation theory to study a hard sphere fluid in an array of permeable obstacles. The fluid-matrix potential they applied12 qualitatively correctly describes the real ligand-membrane potential and can be, with an appropriate choice of adjustable parameters, fitted to describe the behaviour of ionic ligands. The authors tested the theory against grand canonical Monte Carlo simulation results and obtained excellent agreement for both, structural and thermodynamic properties.
In this work, the ligand-membrane potential applied by Bracamontes et al.12 was used and appropriate parameters were chosen to study the partitioning of simple monovalent ions in this model membrane using the ROZ/HNC approximation. The applicability of the method for studying ionic membrane equilibria was established. The paper is organised as follows: after a brief introduction, the model and method are described, following by the results and discussion. The conclusions are given at the end.
2. Model and Method
2.1. Primitive Model Electrolyte Distributed in a Model Membrane
The model under investigation consists of two components: the membrane subsystem is modelled as a quenched one-component neutral fluid, which is called the matrix, and the penetrating electrolyte is an annealed ionic fluid which thermally equilibrates in the presence of the matrix species. The notation used in this paper is the same as used before: the superscripts 0 and 1 correspond to the matrix and the annealed fluid species, respectively.7 The matrix is representing the biological cells surrounded by cell membranes and is formed by neutral hard spheres with the diameter o0 = 26 A. Note that the size of the matrix particles does not mimick any system in particular; the size of the lipid bilayer membranes, however, is usually considered to be between 20 and 30 A.13 The concentration of matrix particles was the same in all cases, and namely c0 = 0.02 mol dm-3. The interaction potential between matrix particles is the hard sphere potential:
u
oo
foo.r < = 10, r >
< cr0
(1)
The structure of the matrix corresponds to an equilibrium state of a hard sphere fluid of concentration c0. The matrix structure was obtained by solving Ornstein-Zerni-ke equation in the hypernetted chain (HNC) approximation14 without any loss of generality. Other closures can be used for this purpose as well.
The model of the annealed fluid is similar to that used in previous works15-19; it corresponds to an electro-neutral system of charged hard spheres with charges
z+ = +1 e0 and z-= -1 e0, and diameters o+ = 3.87 A, and = 3.62 A, mimicking sodium (Na+) and chloride (Cl) ions, respectively.20 The structureless solvent in this model is considered as a dielectric continuum with a dielectric constant er = 78.5, corresponding to water at TL = 298.15 K. The fluid-fluid interaction potential reads:
(2)
where LB = $e0j4neoer is the Bjerrum length, and equals to 7.14 A for aqueous solutions at 298.15 K, studied here. /3 = 1/kBT1 where kB is the Boltzmann constant and T1 the absolute temperature of observation.
The fluid-matrix interaction is determined by electrostatic lipid bilayer potential and is taken from the paper of Bracamontes et al.12:
(3)
where ey is the depth of the attractive part of the potential, U0 is the height of the barrier, and wy is the half-width of the barrier. The potential is centred at o0/2. By adjusting the parameters one can go from almost entirely repulsive potential (small values of eij) to attractive one (larger values of eij). The purpose of this paper was not to exactly reproduce the ion-membrane potential, but to capture its main features, that is a large barrier for both ionic species, which is slightly lower for the anions than for the cations, and, on the other hand, a narrow and shallow attraction.13 The potential for the parameters used is shown in Figure 1: Pe+0 = Pe_0 = 0.5, w+0 = w_0 = 2 A, and pu;°= 5, PU-0 = 2.
6
+ ion	
	- ion
	1
	' 1 1
	f 1
	
0 2 4 6 8 10 12 14 16 18 20 Distance from the center of the cell in A
Figure 1. The ion-matrix potential for the following parameters: Pe+0 = Pe-0 = 0.5, w+0 = w-0 = 2 A, and pU+00 = 5, PU-0 = 2. c0 = 26 A. The solid line represents the potential for the cation, and the dashed line for the anion. The membrane is centered at c„/2 = 13 A.
2. 2. Replica Ornstein-Zernike Equations in the HNC Approximation
The structure and thermodynamics for the model described were calculated using the replica Ornstein-Zer-nike equations with the HNC closure. The method was described in details elsewhere7' 15-19' 21-28 so only the main features will be discussed here.
First the matrix structure was obtained in terms of total pair correlation functions h00, by solving the Orn-stein-Zernike (OZ) equation for a one component fluid:
h00-c00 = c00*p°h00
(4)
where c00 is the direct correlation function' and the symbol * denotes convolution in r space. p0 = c0NA is the number density of the matrix particles, and NA is the Avogadro number. The OZ equation was solved in the hypernetted-chain (HNC) approximation: ln[h00 + 1] = -fa00 + h00 - c00.
The ROZ equations for the fluid-fluid and fluid-matrix correlations read7'17:
where c(s) denotes the short range part of the direct correlation function coming from the renormalization procedure. The mean activity coefficient of the +1:-1 annealed electrolyte was calculated as: P/Jex = ln/± = -Jr(ln/+ + lnj-).
Further, we calculated a coordination number of ions i inside the matrix as12:
Nci(R) = 4n pî fRgi0°r2 dr
(8)
where g10 = h10 + 1, is the ion-matrix pair distribution function. The excess charge distribution was calculated as
NJr) = Nr+(r) - N(r).
3. Results and Discussion
First we present the ion-matrix profiles, i. e. the fluid-matrix pair distribution functions (pdfs), gfm(r), for two different electrolyte concentrations in two different matrices that differ in the width of the barrier, w. In Figure 2 the ion-matrix profile, gfm(r), is shown for two concentrations
(5)
The h12 and c12 are the blocking contributions to the total and direct correlation function, respectively.7,17
In the present contribution the ROZ equations were supplemented by the HNC closure in the form:1519
hmn _ exp	+ ftmn _ £mn] _ ±
h12 = exp[/l12 - C12] - 1
(6)
where the superscripts m and n take values 0 and 1.
Before solving numerically, the set of equations 5 and 6 was renormalized; this procedure, for a non-charged matrix, is in details described in17. The ROZ equations were solved by direct iteration with 1015 integration points, and integration step 0.005 A.
The excess chemical potential of the annealed ions was calculated via:18
of NaCl model solution (a: c1 = 0.1 M, and b: c1 = 1.0 M) within the matrix with a relatively narrow barrier, characterized by the parameters: fae+0 = fae-0 = 0.5, w+0 = w-0 = 2 A, and fHJ+0° = 5, fiU^ = 2. The full line represents the profile for the cation, and the dashed line for the anion.
There are practically no cations present within the barrier representing the membrane. Due to the lower barrier for the anions, there are some of them present within the barrier, causing the concentration of cations to be slightly higher than that of the anions within the matrix particles (cells). The effect is more pronounced for a higher electrolyte concentration (e. g., Figure 2. b). This is in qualitative agreement with the experimental observations: the small ions are virtually excluded from the membrane walls where their energy barrier is at least 40 kcal/mol.13 Also, the lipid membrane barriers are at least three times more permeable to chloride compared to sodium or potassium.2930
a)
1.6
b)
1.6 _ 1.2
e
0.8 0.4 0
Jl		
	i Í » t v..,./	
10 ,» 15 r/A
20
25
Figure 2. The ion-matrix pair distribution function, gfm(r) for the following parameters: Pe+0 = Pe-0 = 0.5, w+0 = w-0 = 2 A, and = 5, P^0° = 2; C0 = 26 A, c0 = 0.02 mol dm-3. The solid line represents the gfm for the cations, and the dashed line for the anions. a) Cj = 0.1 mol dm-3, b) Cj = 1.0
J0 - ¿.u n, ^
mol dm-3.
As the concentration of the electrolyte increases, the cations get unequally distributed on both sides of the barrier (membrane); the "contact" value of gf+(r), i. e. the value of gf+(r) at a distance that corresponds to vanishing repulsion, is higher on the inner side of the barrier than on the outside, despite the symmetry of the fluid-matrix potential.
The results for the fluid-matrix pair distribution functions for a case of a matrix with somehow broader barrier (w+0 = w-0 = 4 A) are shown in Figure 3. The shapes of the functions are very similar as observed for a narrower barrier above. Since the ions are here excluded from a larger space, the ions within the matrix particles start forming some kind of layers (the oscila-tion in the gfm function on the inner side of the membrane).
The distribution of the ions within the matrix can also be represented as a coordination number (defined by eq. 8) of different ion species. Figure 4 shows the coordi-
nation number for cation and anions, as a function of a distance from the membrane center, and namely for the two cases presented in Figure 2.
As expected, the coordination numbers of both ionic species are lower in the case of lower concentration of the NaCl (0.1 mol dm-3). Similar trends were observed by studying more detailed lipid membrane models via molecular dynamics simulations.31
The coordination numbers for the model electrolyte within a matrix with somehow broader barrier (w+0 = w-0 = 4 A) are shown in Figure 5.
Similarly, as in the case of the fluid-matrix pdfs, the width of the membrane barrier does not change the observed trends in behavior qualitatively. While there are slightly more positive than negative ions on the inner side of the membrane, the opposite is true in the close proximity of the outside edge of the barrier. As expected, there are less ions within the matrix in the case of broader membrane barrier (wi0).
0	5
Figure 3. Same as Figure 2, but for w+0 =
20
--4 k.
w
a) 0.5		b) 5		
				
	/f'J		i/—"	
0.4	/	4	h	
zFO.3	J	"□3	/	
	J	•z.		
0.2	y	2	y	
0,1		1		
0		0		
6
r/A
10
12
6
r/A
10
12
Figure 4. The coordination number for cations (NC+(r)), and anions (Nc_(r)) gfm(r), for two different electrolyte concentrations (a: Cj = 0.1 mol dm-3, and b: c1 = 1.0 mol dm-3) for the matrix defined with the parameters: Pe+0 = Pe_0 = 0.5, w+0 = w-0 = 2 A, and = 5, P^0 = 2; c0 = 26 A, c0 = 0.02 mol dm-3. The solid line represents the function for the cation, and the dashed line for the anion.
a)
b)
0 2 4 6 8 r/A
Figure 5. Same as in Figure 4, but for w+0 = w-0 = 4 Â.
In many cases, the membranes do not directly perceive the ionic conditions of the surrounding bulk solution but are in contact with a modified medium, the composition of which is determined by the interaction between the solution and the outer wall of the membrane.32 We therefore analysed the net charge in this model membrane. Due to unequal degree of penetration of positive and negative ions into the model membrane, a net charge inside the membrane is observed, Nch(r) = NC+(r) - NC (r), see Figures 6 (wi0 = 2 À), and 7 Cwi0 = 4 À+).
The inside of the cell is slightly positively charged, while within the membrane (barrier), which is slightly lower for negative ions, the charge becomes negative, causing the charge oscillation in the outside solution. Note, however, that, although it depends on the electrolyte concentration, as well as on the width of the membrane barrier, the net charge on both sides of the barrier is very small even in the case of the highest electrolyte concen-
tration studied here (i. e., = 1.0 mol dm 3, Figures 6 and 7 b).
Finally, in Figure 8 we present the mean activity coefficient, y± (a), and the chemical potential, fy = P/jex + ln(c1/mol dm-3) (b), as a function of the invading electrolyte concentration, in the two different matrices studied. Both properties are also shown for the bulk model electrolyte, for comparison. The presence of the matrix particles does not change the qualitative behavior of the activity coefficient of the invading electrolyte as a function of concentration. At low concentrations of the electrolyte the mean activity coefficient decreases with the electrolyte concentration due to strong attractive forces, and starts increasing at larger concentrations. The presence of the membrane particles increases the activity coefficients, the effect is stronger for a broader membrane barrier. As a results, the chemical potential of the invading electrolyte increases, in comparison with the bulk electrolyte of the same concentration.
a)
b) ■
Figure 6. The net charge, Nch(r), within the matrix characterized by Pe+0 = Pe-0 = 0.5, w+0 = w-0 = 2 A, and = 5, P^0° = 2; = 26 A, c0 = 0.02 mol dm-3, for two different electrolyte concentrations (a: c = 0.1 mol dm-3, and b: cj = 1.0 mol dm-3).
a)	b)
Figure 7. Same as in Figure 6, but for w+0 = w-0 = 4 A.
a)
b)
Figure 8. The mean activity coefficient, y± (a), and the chemical potential, PjU (b), as a function of the invading electrolyte concentration, in the two different types of membranes, wi0 = 2 A (solid lines), and wi0 = 4 A (dashed lines), shown here. The dotted lines show the results for the bulk (unperturbed) electrolyte.
4. Conclusions
A simple model of the ionic membrane potential was used to investigate the applicability of the Replica Ornstein-Zernike methodology for studying ionic equilibria. The parameters for the primitive model electrolyte were chosen to mimic the NaCl aqueous solutions, while the matrix was represented as a subsystem of quenched hard spheres with the membrane-like surface. Despite of the simplicity of the ion-membrane potential, it has the flexibility to fine-tune the distribution of positive and negative ions within the membrane by changing the parameters of the model. The results are qualitatively similar to those obtained by more sophisticated models, as well as to experimental observations. The model in the present form is rather simple. However, it permits several extensions that can put it closer to the systems of experimental research. In particular, one can exploit different sets of parameters for fluid-membrane interactions to tune partitioning of ions in the membrane with respect to the bulk solution. On the other hand, it is of interest to employ the model with the purpose of studying the additional uncharged species delivery into membranes. Some of these issues are now under study in our laboratory.
A more systematic investigation, where the limitations of the ROZ/HNC theory (its thermodynamic consistency, the convergence depending on the diameter of the quenched particles and their packing fraction, etc.) will be thoroughly explored, is already in progress.
5. Acknowledgements
I thank V. Vlachy and O. Pizio for helpful comments and critical reading of the manuscript, and M. Luksic for usefull discussions. The useful of the Ministry of Science and Education of Slovenia under grants P1-0201, and J1-4148 is gratefully acknowledged.
6. References
1.R.	B. Gennis, Biomembranes. Molecular Structure and Function, Springer: London, 1989.
2.	R. H. Guy, D. H. Honda, Int. J. Pharm., 1984, 19, 129-137.
3.	V. M. Knepp, R. H. Guy, J. Phys. Chem., 1989, 93, 68176823.
4.	B. J. Forrest, J. Mattai, Biochemistry, 1985, 24, 7148-7153.
5.	A. Enders, Biochim. Biophys. Acta, 1990, 1029, 43-50.
6.	M. Ulbricht, Polymers, 2006, 47, 2217-2262.
7.B.	Hribar-Lee, M. Luksic, V. Vlachy, Annu. Rep. Prog. Chem., Sect. C, 2011, 107, 14-46.
8.	P. Bryk, A. Patrykiejew, J. Reszko-Zygmunt, S. Sokolowski, Mol. Phys., 1999, 96, 1509-1516.
9.	D. Boda, D. Henderson, R. Rowley, S. Sokolowski, J. Chem. Phys., 1999, 111, 9382-9388.
10.	K. E. Forsten, R. E. Kozack, D. A. Lauffenburger, S. Subra-maniam, J. Phys. Chem, 1994, 98, 5580-5586.
11.	D. Andelman, Electrostatic Properties of Membranes: The Poisson-Boltzmann Theory, in R. Lipowsky, E. Sackmann, ed.: Handbook of Biological Physics, vol. 1, Elsevier, London, 1995.
12.	L. I. Bracamontes, O. Pizio, S. Sokolowski, A. Trokhymc-huk, Mol. Phys., 1999, 96, 1341-1348.
13.	B. H. Honig, W. L. Hubbell, R. F. Flewelling, Ann. Rev. Biophys. Biophys. Chem., 1986, 15, 163-193.
14.	J.-P. Hansen, I. R. McDonald, Theory of Simple Liquids, Elsevier, London, 2006.
15.	B. Hribar, O. Pizio, A. Trokhymchuk, V.Vlachy, J. Chem. Phys., 1998, 109, 2480-2489.
16.	B. Hribar, V. Vlachy, A. Trokhymchuk, O. Pizio, J. Phys. Chem. B, 1999, 103, 5361-5369.
17.	B. Hribar, V. Vlachy, O. Pizio, J. Phys. Chem. B, 2000, 104, 4479-4488.
18.	B. Hribar, V. Vlachy, O. Pizio, J. Phys. Chem. B, 2001, 105, 4727-4734.
19.	B. Hribar, V. Vlachy, O. Pizio, Mol. Phys., 2002, 100, 3093-3103.
20.	J. P. Simonin, O. Bernard, L. Blum, J. Phys. Chem. B, 1998, 102, 4411-4417.
21.	O. Pizio, S. Sokolowski, J. Phys. Studies, 1998, 2, 296-321.
22.	O. Pizio, in: Computational Methods in Surface and Colloid Science (Surfactant science series, vol. 89), M. Borowko, ed., Kluwer, Marcel Dekker, New York, 2000, 293-345.
23.	A. D. Trokhymchuk, O. Pizio, M. F. Holovko, S. Sokolowski, J. Phys. Chem. 1996, 100, 17004-17010.
24.	V. Vlachy, H. Dominguez, O. Pizio, J. Phys. Chem. B, 2004, 108, 1046-1055.
25.	T. Urbič, V. Vlachy, O. Pizio, K. A. Dill, J. Molec. Liquids, 2004, 112, 71-80.
26.	M. Lukšič, V. Vlachy, O. Pizio, Acta Chim. Slov., 2006, 53, 292-305.
27.	O. Pizio, S. Sokolowski, Phys. Rev. E, 1997, 56, R63.
28.	A. Trokhymchuk, O. Pizio, M. Holovko, S. Sokolowski, J. Chem. Phys., 1997, 106, 200-209.
29.	D. W. Deamer, J. Bramhall, Chemistry and Physics of Lipids, 1986, 40, 167-188.
30.	G. Eisenman, R. Horn, J. Membrane Biol., 1983, 76, 197225.
31.	A. A. Gurtovenko, M. Miettinen, M. Karttunen, I. Vattulai-nen, J. Phys. Chem. B, 2005, 109, 21126-21134.
32.	H. Sentenac, C. Grignon, Plant Physiol., 1981, 68, 415-419.
Povzetek
Namen dela je bil proučiti možnost aplikacije »Replika« Ornstein-Zernikovih enačb (ROZ), ki se običajno uporabljajo za študij sistemov z delno zamrznjenimi prostostnimi stopnjami, za študij membranskih ravnotežij. Za študij sem uporabila enostaven McMillan Mayerjev model +1:-1 elektrolita, membrane pa je predstavljal skupek »zamrznjenih« togih kroglic. Potencialna funkcija med raztopino elektrolita in membrano je bila izbrana tako, da je čim bolje opisala eksperimentalna opažanja. S pomočjo ROZ/HNC približka sem izračunala razporeditev ionov ob membrani, kot tudi adsorp-cijske izoterme. Rezultati se kvalitativno dobro ujemajo z eksperimentalnimi podatki.