J. Phys. Chem. B 2000, 104, 5161-5170
5161
Complete Asymptotic Solution of Cylindrical and Spherical Poisson-Boltzmann Equations at Experimental Salt Concentrations Irina A. Shkel, Oleg V. Tsodikov, and M. Thomas Record, Jr.* Departments of Chemistry and Biochemistry, UniVersity of Wisconsin-Madison, Madison, Wisconsin 53706 ReceiVed: January 7, 2000; In Final Form: March 9, 2000
We report an exact analytic representation of the nonlinear Poisson-Boltzmann (PB) potential as a function of radial distance from a cylindrical or spherical polyion in solutions containing a symmetrical electrolyte, in the form of an asymptotic series in elementary functions, generally valid at radial coordinates larger than Debye length. At sufficiently high salt concentrations, where the ratio of Debye length (κ-1) to the polyion radius (a) is sufficiently small ((κa)-1 e 1), the asymptotic series is valid at any distance from the polyion surface. This analytic representation satisfies exactly the complete nonlinear Poisson-Boltzmann equation, subject to the boundary condition on the derivative of potential at infinity, and therefore contains one integration constant, which in this salt range we determine to an accuracy of order (κa)-2. Because it explicitly introduces for the first time all the terms which arise due to nonlinearity of the PB equation, this analytic representation clarifies the connection between the exact solution of the PB equation and various approximations including the Debye-Hu¨ckel approximation (the solution of the linearized PB equation). From these considerations we obtain a new approximate solution designated “quasi-planar” and expressed in elementary functions, which we show to be accurate at any distance from the polyion surface at typical experimental salt concentrations (e.g., 0.1 M 1:1 salt concentration for double-stranded DNA, where the PB equation retains its accuracy by comparison to Monte Carlo simulations). We apply our analysis to the calculation of the electrostatic free energy and the salt-polyelectrolyte preferential interaction (Donnan) coefficient (Γ).
Introduction The Poisson-Boltzmann (PB) equation is widely used to calculate the electric potential around a macromolecular or colloidal charged particle in electrolyte solutions, beginning with Gouy-Chapman (1910) double layer theory.1 Exact analytic solutions are known for symmetric electrolyte in solution with an infinitely dilute planar polyion,1 and subsequently, for a cylindrical polyion in salt-free solution.2 However, the analytic solution of the PB equation for highly charged spherical or cylindrical surfaces (applicable to micelles, nucleic acids and other polyelectrolytes) in added salt, has remained a challenge. Benham3 showed that the exact solution of the cylindrical PB equation must be a new transcendental function. The DebyeHu¨ckel (DH) linearization4 was the first approximate analytic treatment of the equation, applicable to spherical or cylindrical polyions only at high salt concentrations and/or low surface charge density. Although numerical solutions of three-dimensional finite-difference PB equations for any specified biopolymer structure are now generally possible,5 analytical solutions for simple spherical or cylindrical particles remain of great importance. For example the counterion condensation (CC) theory6,7 has been uniquely valuable because it provides simple analytic expressions for ion-polyion thermodynamic quantities, polyion surface potential, and surface counterion concentration in terms of fundamental polyion structural parameters and salt concentration, applicable at sufficiently low salt concentrations for use in analytic theory to predict the effects of these variables on processes involving rodlike polyions. Thermodynamic * Corresponding author: Department of Biochemistry, University of Wisconsin-Madison, 433 Babcock Dr., Madison, WI, 53706-1544. Phone: (608) 262-5332. Fax: (608) 262-3453. E-mail:
[email protected].
expressions derived as low-salt (CC) limiting laws are consistent with corresponding PB limiting law results,8 but an analytical representation of the PB potential is needed to extend them to high salt concentration. With this motivation, the present work develops an analytic representation of the solution of the PB equation in terms of an asymptotic series in elementary functions. Such asymptotic series solutions often are more useful than convergent series because fewer terms are necessary for the accurate approximation of the exact solution.9 In this paper, we present two asymptotic forms of the exact solution to the cylindrical and spherical nonlinear PB equations at infinite dilution of the polyion and at sufficiently high salt concentration. In the Results Section, we present complete asymptotic expansion of the exact solution. The expansion is an infinite asymptotic series in elementary functions, which satisfies exactly (without approximation) the nonlinear PB equation and the boundary condition on the decay of the potential at infinity. Thus, all coefficients in this series are expressed through only one constant, which must be determined from the boundary condition on the polyion surface. The Results Section presents a procedure for calculation of this constant to an accuracy of any power of the small parameter (κa)-1. We also show that at sufficiently high salt concentration ((κa)-1 < 1) the solution is valid at any distance from the polyion surface. An alternative asymptotic form of the solution, designated the “quasi-planar” approximation, takes into account terms that are most significant at the polyion surface. This approximation is valid in all space at (κa)-1 < 1, which includes the experimental range of salt concentration Cb > 0.1 M in double-stranded DNA solutions. In the Discussion Section, the complete asymptotic solution is used to compare various known
10.1021/jp0001182 CCC: $19.00 © 2000 American Chemical Society Published on Web 05/03/2000
5162 J. Phys. Chem. B, Vol. 104, No. 21, 2000
Shkel et al.
approximations with the quasi-planar approximation. We show that at salt concentrations where (κa)-1 < 1, the quasi-planar approximation is always better for calculations of surface potential and relevant thermodynamic properties (such as the electrostatic free energy and the preferential interaction coefficient Γ) than widely used approximations including solutions of linearized or planar PB equations. The solution presented in this paper is obtained for symmetrical (z:z) salt of any charge. However quantative applications of the PB equation are generally restricted to the case of 1:1 salt (z ) 1) or mixtures containing excess 1:1 salt because of strong correlations between multivalent ions neglected in the PB equation.10-12 Poisson-Boltzmann Equation. The form of the PoissonBoltzmann equation describing the electric potential around a uniformly charged cylindrical (j ) 1) or spherical (j ) 2) particle immersed in a symmetrical (z:z) electrolyte at infinite dilution of the polyion is
1 d jdy r ) κ2 sinh(y) rj dr dr
( )
(1)
dy σ* )dr r)a a
|
(2)
dy )0 rj dr r)∞
(3)
|
where y t zeψ/kBT, κ2 t 2nbz2e2/DD0kBT, σ* t σeza/DD0kBT. Here, ψ and y are actual and reduced potentials respectively, r is the radial coordinate, nb is the bulk (zero potential) concentration of either cation or anion (which are equal to one another because the polyion is at infinite dilution), a, σ*, and σ are the radius, reduced, and actual surface charge densities of the polyion, κ is the inverse of the Debye length, e is the proton charge, z is the cation valency, D and D0 are the dielectric constants of the bulk solution and vacuum respectively, kB is the Boltzmann constant, and T is the temperature. Surface charge density, σ, is the ratio of the total polyion charge to its surface area. For a cylindrical polyion σcyl ) e/(2πab), where b is the axial charge separation (length per elementary charge). For a sphere σsph ) Q/(4πa2), where Q is the charge of the sphere. Then reduced charge density is
σ/cyl )
2
ez ) 2ξz, 2πDD0kBTb
σ/sph )
Qez 4πDD0kBTa
(4)
where ξ is the reduced linear (axial) charge density of the cylinder. Charge can be distributed on the surface or in the volume of the cylinder or the sphere as long as the symmetry of the distribution is preserved. In the PB approximation, concentrations of positively (n+), and negatively (n-) charged ions are related to potential, n( ) nb exp(-y), and to each other through the bulk concentration nb (where y ) 0): n+n- ) n2b. The concentration n is the number of ions per unit volume (1/m3). Molar concentration (C, mol/L) is related to n by n ) CNA × 103, where NA ) 6.02 × 1023 is Avogadro’s number. The planar PB equation corresponds to the case j ) 0 in eq 1 when the inner boundary condition (eq 2) is applied at r ) 0 and the surface charge density is σ/plane ) e2z/(DD0kBTa) where in the planar case a is redefined as a ) A1/2, where A is the surface area per elementary charge. We consider conditions where the ratio of Debye length (κ-1), to the particle radius (a) is sufficiently small, and introduce it as a new parameter:
≡ (κa)-1 < 1
(5)
It is useful to rewrite the equation and boundary conditions eqs 1-3 in terms of the radial coordinate reduced by the product of the particle radius and parameter , X ) r/(a).
1 d j dy X ) sinh(y) Xj dX dX
(6)
dy dy ) -σ*, Xj )0 dX X)1/ dX X)∞
(7)
( )
|
|
The equation and boundary conditions (eqs 6-7) are formulated through two unitless parameters, and σ*. Although the initial set of parameters describing physical characteristics of the system under consideration is more extensive (nb, a, b, DD0, kBT, ez, Q), the two unitless parameters, and σ*, are the minimal possible set of independent parameters for this problem. The regime of < 1 for a polyion with fixed charge and geometry corresponds to the range of salt concentrations nb >DD0kBT /(2z2e2a2). For double-stranded DNA in aqueous 1:1 salt solution, < 1 corresponds to Cb > 0.1 M. Asymptotic expansions of a function at X f X0 (cf. Appendix A, refs 9, 13, 14, and 15) are very useful for function calculations since a few terms are usually enough for good accuracy if X is sufficiently close to X0.9 The linearized cylindrical PB equation
|
d2y 1 dy ) y, + dX2 X dX
dy )0 dX X)∞
is a good example of the practical usefulness of asymptotic series. The solution of this equation and boundary condition is the modified Bessel function of zero order y ) MK0(X) where M is a constant.16 This function cannot be expressed in a closed form in terms of elementary functions. K0(X) has the following complete asymptotic expansion at large X (Xf∞):9
K0(X) )
x
π
e-X
2X
i)∞
∑ i)0
((i - 1/2)!)2 i!(-2X)i
(8)
This series is not convergent, but it is an exact solution because it satisfies the linearized PB equation exactly. The series eq 8 truncated to a few terms can be used for the calculation of K0(X) with a high degree of accuracy. The larger the value of X, the better the accuracy of the truncated series.9 By analogy with the linearized PB equation, it is reasonable to expect that the exact solution of the nonlinear PB equation has similar complete asymptotic expansion, containing the complete asymptotic expansion of K0(X). Results Complete Asymptotic Solution of PB Equation. Here we develop the complete asymptotic solution of the nonlinear PB equation at large values of radial coordinate (X f ∞). It can be constructed by starting from complete asymptotic expansion for the Bessel function of zero order, K0(X), eq 8. At large X, the potential y is small, and the righthand side of the PB equation can be expanded in a power series of y:
( )
n)∞ dy Xj ) sinh(y) ) y2n+1/(2n + 1)! j dX dX n)0 X
1 d
We seek a solution in the form
∑
(9)
Complete Asymptotic Cylindrical and Spherical P-B Solutions i)∞
y)
yi ∑ i)0
where each next term is asymptotically small compared to previous one, i.e.,
yi+1(X)
lim
yi(X)
Xf∞
f0
for any i. Thus the first term, y0, can be found from the equation as its approximate solution subjected to condition that the remainder decreases to zero faster than y0 itself, y(X) ) y0(X) + R0(X), and
R0(X)
lim
y0(X)
Xf∞
f0
J. Phys. Chem. B, Vol. 104, No. 21, 2000 5163 in the solution (k g 1) appear because of the nonlinearity of the PB equation. (In the planar case, all coefficients Ai,2k are zero except A0,2k, i.e., eq 10 contains terms that correspond to i ) 0 only). Series eq 10a is an asymptotic series at infinity since each function decays faster than the previous one as X f ∞. It obviously satisfies the boundary condition at infinity. To show that eq 10a is the complete asymptotic expansion of the exact solution of the PB equation, we need to show that it satisfies the PB equation. We substitute the series given by eq 10a into eq 9 and equate coefficients of terms e-(2k+1)XX-(2k+1)j/2-i. To prove that the series eq 10a is a solution, one can show that no terms with functional forms differing from e-(2k+1)XX-(2k+1)j/2-i are present and that all coefficients Ai,2k are determined uniquely from the substitution. After equating coefficients of terms e-(2k+1)XX-(2k+1)j/2-i, all coefficients Ai,2k are expressed in terms of A00 by induction:
Ai,0 ) -Ai-1,0(i - 1 + j/2)(i - j/2)(2i)-1,
After the first N terms in the expansion are known, the procedure of obtaining the next one is similar: one substitutes i)N-1
y(X) )
Ai,2k )
yi(X) + yN(X) + RN(X) ∑ i)0
into the original equation, and finds an approximate solution of the resulting equation for yN(X) subject to the condition
lim
RN(X)
Xf∞
yN(X)
f0
It is straightforward to verify that starting from the function e-XX-1/2 for the linearized cylindrical PB equation and following the algorithm above, one obtains the complete asymptotic expansion of the Bessel function K0(X), eq 8. By following the same procedure for the nonlinear PB equation, we find that the complete asymptotic expansion of K0(X) is a part of the complete asymptotic expansion of the solution of the nonlinear PB equation, and the subsequent terms have the form of e-mXX-n/2, where m and n are integer numbers such that m > 1, and n g m. Thus, the complete asymptotic expansion of the exact solution of eq 6 subject to one boundary condition at infinity (eq 7) is found as a series in the inverse powers of coordinate X multiplied by exponentially decaying functions:
y)
k)∞ -(2k+1)X i)∞A
∑
e
∑
k)0 X(2k+1)j/2 i)0
( (
) )
nek
1
q
1
1
In eq 10b the first row of terms (k ) 0) is a well known solution of the linearized PB equation, the Debye-Hu¨ckel approximation.4,16 In the cylindrical case, the first row of eq 10b is the complete asymptotic expansion of the modified Bessel function of zero order, A00(2/π)1/2K0(X). In the spherical case, the coefficients Ai0 are equal to zero for all i g 1, so that only the first term, A00e-x/X, remains in the first row. All other rows
q
1
q
-Ai-1,2k(2k + 1)(2i - 2 + 2kj)
(11)
-Ai-2,2k(i - 2 + (2k + 1)j/2)(i - 1 - j + (2k + 1)j/2)] Here the sum ∑′sp is the sum over all sequences of integer sp, 1 e p e q subject to the following conditions: ∑psp ) 2n + 1; ∑psplp ) k - n; ∑pspmp ) i. It is important to note that Ai,2k ) A2k+1 00 Bi,2k, where the terms Bi,2k depend only on i and k and not on A00. Equation 11 is used to calculate exact coefficients Bi,2k by a computer program written in this laboratory (available upon request). The first several coefficients Bi,2k are presented in Table 1. The single unknown coefficient A00 is determined from the boundary condition on the surface (see the following section). The range of validity of the solution (eq 10a) is the range of the radial coordinate where each successive term in eq 10a decays faster than the previous one. Two comparisons are required: between terms with different powers of X at the same exponent, and between terms with different exponents.
Ai+1,2ke-(2k+1)X X(2k+1)j/2+i (2k+1)j/2+i+1
X(2k+3)j/2+i
(10b)
q
p
-(2k+1)X
Ai,2ke
Ai,2k+2e-(2k+3)X X(2k+1)j/2+i
)
1
′Ams ,2l ...Ams ,2l ∑ ∑ s !...s ! n)1 s
(10a)
Xi
A10 A20 e-X + 2 + ... + A00 + j/2 X X X A12 A22 e-3X A02 + + 2 + ... + ... X X3j/2 X -(2k+1)X A A2,2k e 1,2k A0,2k + + 2 + ... + ... (2k+1)j/2 X X X
(
[
(2k + 1)2 - 1
X
i,2k
or equivalently
y)
1
∼
1 ; inequality 13 is also satisfied at these r, therefore r/a > represents the range of validity of the complete asymptotic solution eq 10a. At high salt concentrations ( < 1), the asymptotic solution is valid over the entire space (a e r < ∞). Note that these two are necessary conditions. For practical calculations, series eq 10a is truncated and its range of applicability is determined by estimation of the remainder. Quasi-Planar Solution. For the series (eq 10b) at the polyion surface, where X ) 1/, we observe that all first column terms (A0,2ke-(2k+1)XX-(2k+1)j/2) are of the same order of magnitude at the surface because A00j/2e-1/ ∼ 1 (as shown below). All second column terms (A1,2ke-(2k+1)XX-1-(2k+1)j/2) are of the order of at the surface; third column terms are of order 2 etc.
5164 J. Phys. Chem. B, Vol. 104, No. 21, 2000
Shkel et al.
TABLE 1: Relative Expansion Coefficients Bi,2k ) Ai,2k/A2k+1 in the Cylindrical and Spherical Poisson-Boltzmann Solutions 00 (Eqs 10-11) cylindrical (j ) 1)
spherical (j ) 2)
k
i)0
i)1
i)2
i)3
i)4
i)0
i)1
i)2
i)3
0 1 2 3 4
1 1/ (3×24) 1 /(5×28) 1 /(7×212) 1/ (9×216)
-1/23 -3/27 -3/211 -3/215 -3/219
9 7 /2 71/ 11 2 89 15 /2 107 19 /2 125/ 23 2
-75/210 -3215/3×214 -1485/218 -1965/222 -2517/226
3675 15 /2 79 521/ 19 2 114 071 23 /2 157 237 27 /2 211 275/ 31 2
1 1/ (3×24) 1 /(5×28) 1 /(7×212) 1/ (9×216)
0 -1/25 -1/29 -1/213 -1/217
0 7/ 7 2 7 /(3×29) 205 2 16 /(3 ×2 ) 2413/ (5×32×221)
0 -15/27 -35/(3×210) -863/(33×215) -90271/(5×33×223)
Therefore, close to the surface, we should sum only the firstterms in each parentheses (A0,2ke-(2k+1)XX-(2k+1)j/2) and neglect the rest. In contrast, at r f ∞, it is better to sum only the first row (Debye-Hu¨ckel approximation) and neglect the other rows. Since each column has an infinite number of terms, evaluation would be greatly facilitated by determining the functional form of each sum. We use the method of multiple scales,13 which yields the solution of the nonlinear PB equation (eq 6) in the form of a uniformly valid asymptotic series in powers of , i)∞
y)
iyi ∑ i)0
where all yi are of the same order with respect to . The details of the derivation of this solution are presented in Appendix B. If the potential y is approximated by only the zero order term in this expansion, then y ) y0 + O(), where
y0 ) 2 ln
(
)
1 + ˜t x-j/2 e-(x-1)/ 1 - ˜t x-j/2e-(x-1)/
(14)
where ˜t ) -p˜ -1 + (1 + p˜ -2)1/2, p˜ ) σ*/(2 + j), x is the reduced radial coordinate (x ) r/a ) X) and O() is the symbol for the remainder meaning that the remainder is on the order of as f 0. Constant ˜t is calculated by applying the boundary condition at the surface to the zero order term (y0) only. We call the approximate solution obtained by this method “quasiplanar” (QP) asymptotic expansion because its leading term, eq 14, is similar to the solution of the planar PB equation.1 Indeed, the case j ) 0 of eq 14 is the exact planar solution. The differences between eq 14 and the planar solution are in the power of radial coordinate multiplied by the exponential functions and in expressions for constants ˜t, p˜ , and q˜ , symbols chosen to distinguish these quantities from analogous constants t, p, and q in the planar solution.17,18 If y0 is expanded at large radial coordinate x and expressed in terms of the radial coordinate X ) x-1, it has only terms of various powers of e-XX-j/2, and this expansion yields exactly the first column of the complete asymptotic solution eq 10b. Thus, y0 can be considered a truncated form of the complete asymptotic solution where terms from only the first column (i ) 0) are taken into account. The expansion of y0 and the complete solution should have equal coefficients at each power of the function e-XX-j/2, which yields the following expression for the coefficient A00:
A00j/2e-1/ ) 4t˜ + O()
(15)
If two terms are retained in the asymptotic expansion then y ) y0 + y1 + O(2) and the expression for y1 is
(
)
jt˜x-j/2 e-(x-1)/ ˜t 2e-2(x-1)/ 2 - j y1 ) c + + 1 2x xj+1 1 - ˜t 2x-je-2(x-1)/
(16)
i)4 0 155/ 9 2 2327/ 2 13 (3 ×2 ) 40 621/ 3 19 (3 ×2 ) 11 313 029/ 3 2 26 (3 ×5 ×2 )
where c1 ) (t˜4 - ˜t2(4 - j/2) - 1 + j/2)/(1 + ˜t2) is the constant found from the boundary condition (∂y1/∂x)|x)1 ) 0. Again, y1 expanded at small e-XX-j/2 yields terms from the column i ) 1 and small (on the order of ) corrections to the terms from the column i ) 0 of the complete asymptotic solution. By comparing coefficients at corresponding terms of column i ) 0, we obtain the next approximation for A00
A00j/2e-1/ ) 4t˜(1 - jc1/4) + O(2)
(17)
In principle, consecutive terms in this approximation can be found and the ith term, yi, yields the ith column of the complete asymptotic solution and corrections of higher order (∼i) with respect of to the previous columns. This allows one to calculate A00 with any accuracy. Quasi-planar expansion, y ) y0 + y1 + O(2), is an expansion in small parameter valid at < 1. For typical parameters of double-stranded DNA (j ) 1, a ) 10 Å, b ) 1.7 Å, σ* ) 2ξ ) 8.4) we find that the two term expansion y = y0 + y1 is accurate to within 1% for potential near the polyion surface at < 1 (Cb > 0.1 M, which includes experimentally relevant conditions) and with the same accuracy for surface concentration at < 0.67 (Cb > 0.2 M). Thermodynamic Applications Examples of applications of the complete solution in the form of the quasi-planar expansion are the calculations of the electrostatic free energy and the salt-polyelectrolyte preferential interaction (Donnan) coefficients, Γ, for a model relevant to DNA solutions: a charged cylinder in symmetrical 1:1 electrolyte. The application of Γ to thermodynamic analysis of saltDNA interactions and the effect of salt concentration on DNA processes have been recently reviewed.19 Within the PB approximation, the Coulombic part of Γ is calculated either by integrating ion distributions over the volume per elementary charge of the polyion20-22 or by differentiating the electrostatic free energy, fel, which is calculated as the work of charging the polyion surface in the presence of the salt ions21,22
) Γcoul i coul Γcoul + + Γ- ) -
r)∞ (ni - ni,b) dV ∫r)a
∂fel , ∂µ
fel ) NAA
∫0σψ(σ′) dσ′
(18) (19)
where i is positive or negative, µ ) RT ln Cb is the ion chemical potential, R is the gas constant, and A is the surface area per elementary charge. More detailed formulae and explanation of different notations are presented in Appendix C. Further discussion deals with Coulombic part of the preferential interaction coefficient as it is defined by eq 18 or eq 19. The superscript “coul” and subscript “-” will be omitted for simplicity, Γcoul - ) Γ.
Complete Asymptotic Cylindrical and Spherical P-B Solutions TABLE 2: Comparison of Preferential Interaction Coefficients Γ for Different Cylindrical Polyions Including Double-Stranded (ds) DNAa and Single-Stranded (ss) DNAb σ*/2 2 4 6 8 10
κa 1 2 3 4 5
Cb (M) s s s s s
-Γnumc -ΓQP,1c
δΓQP,1 δΓPBA,2 (%)d -ΓPBA,2c (%)d
2.411 2.631 2.712 2.755 2.780
0.288 0.248 0.231 0.222 0.216
0.300 0.252 0.233 0.224 0.217
-4.0 -1.8 -1.0 -0.63 -0.47
0.216 0.243 0.229 0.221 0.216
ynum
9.7 1.9 0.73 0.39 0.19
4.2 0.87 1.04 1.47 1.8 2.32 2.75 3.28
0.07 0.1 0.2 0.3 0.5 0.7 1
4.189 3.862 3.245 2.898 2.482 2.223 1.964
0.170 0.182 0.212 0.230 0.257 0.276 0.296
0.178 0.190 0.217 0.235 0.260 0.278 0.297
-5.4 -4.4 -2.8 -2.1 -1.4 -1.1 -0.8
0.152 0.169 0.203 0.224 0.253 0.272 0.293
10.2 7.3 3.7 2.4 1.4 1.0 0.6
2.1 0.73 1.03 1.26 1.63 1.92 2.3
0.1 0.2 0.3 0.5 0.7 1.0
2.990 2.460 2.172 1.833 1.628 1.426
0.257 0.284 0.300 0.323 0.338 0.353
0.272 0.294 0.309 0.329 0.343 0.357
-5.7 -3.9 -3.0 -2.1 -1.6 -1.1
0.215 0.258 0.282 0.310 0.329 0.347
16.5 9.0 6.2 3.8 2.6 1.7
a ds DNA: a ) 10 Å, b ) 1.7 Å, σ*/2 ) ξ ) 4.2. b ss DNA: a ) 7 Å, b ) 3.4 Å, σ*/2 ) ξ ) 2.1. c ΓQP,1 is calculated from eq 21, ΓPBA,2 is from three term expression from Van Aken et al. (ref 21) (cf. eq 47 Appendix C), and Γnum is from the numerical solution. d δΓ ) (1 - Γ/Γnum) 100%.
We derived expressions for fel and Γ using as a starting point the two term quasi-planar (QP) expansion for the surface potential (eq 43, Appendix C).
[
fel,QP,1 ) 2RTp˜ -1 2p˜ ln(p˜ + q˜ ) - 2q˜ + 2 + j q˜ - 1 - 2 ln
(q˜ +2 1))] (20)
(q˜ - 1) 1 1 - q˜ -ΓQP,1 ) + - j 2 2p 4pq˜
(21)
(
2
where p ) σ*/2, p˜ ) p(1 + j/2)-1 and q˜ ) (1 + p˜ 2)1/2. Subscripts QP and 1 indicate that fel is calculated with two term QP expansion for surface potential and that a consistent level of approximation is maintained in the expressions for fel and Γ as described in Appendix C. In the limiting case of weak Coulombic interaction, σ*/2 , 1 (low polyion charge, σ*/2 , 1, or high salt concentration, , 1), eq 21 can be expanded in powers of σ*.
σ* 1 -ΓQP,1 ≈ 2 2(2 + j)2
(22)
The first term of this expansion is the ideal value of Γ in the absence of interactions. The second term provides the functional form of salt dependence of Γ at these conditions (σ*/2 , 1). The accuracy of ΓQP,1 calculated using eq 21 is summarized in Table 2 for the cylindrical case (j ) 1, σ* ) 2ξ). In the range 0.8 < κa < 2, which includes typical experimental salt concentrations Cb ) 0.1 - 0.2 M for single-stranded and doublestranded DNA, the error is less than 5% which is sufficient for comparison with experimental measurements of salt dependence of binding constant.19 Equation 21 is even more accurate at higher values of κa (higher salt concentration). Equation 21 for Γ, derived here from the quasi-planar expansion can be compared with that obtained from the planar
J. Phys. Chem. B, Vol. 104, No. 21, 2000 5165 expansion by Van Aken et al., 1990,21 and given here as eq 47 of Appendix C. The two term expression for ΓQP,1 (eq 21) is simpler in functional form than the three term expression ΓPBA,2. Numerical comparison of ΓQP,1 and ΓPBA,2 is shown in Table 2 for three cases: (i) p ) 2, 1 e κa e 5 (values used in Van Aken et al.21). (ii) Double-stranded DNA, ξ ) 4.2 (a ) 10 Å, b ) 1.7 Å) in 1:1 salt solution at concentration from Cb ) 0.07 M (κa ) 0.87) to Cb ) 1 M (κa ) 3.28). (iii) Single-stranded DNA, ξ ) 2.1 (a ) 7 Å, b ) 3.4 Å) in 1:1 salt solution at a concentration from Cb ) 0.1 M (κa ) 0.73) to Cb ) 1 M (κa ) 2.3). In all three cases the accuracy of ΓQP,1 is better than ΓPBA,2 at 1 < κa < 2. Particularly, in the typical experimental range for both double- and single-stranded DNA (0.1 M < Cb < 0.2 M) ΓQP,1 is twice as accurate as ΓPBA,2, and more accurate than a typical experimental measurement of a Donnan coefficient.23 For both double- and single-stranded DNA both expressions are accurate approximations to the numerical PB solution (less than 2% error) at sufficiently high salt concentration (above 0.3 M for double-stranded DNA, and above 0.5 M for singlestranded DNA). Discussion Complete Asymptotic PB Solution and Comparison of Alternative Linearized (DH) Solutions for PB Potential. Analytical solutions for the dependence of potential around a cylindrical or spherical polyion on distance, salt concentration, and polyion structural quantities have been sought for over 80 years. A very large number of publications (more than a thousand in the past decade) have addressed the PB equation for various symmetries from both analytical and computational standpoints. However an explicit analytical expression in elementary functions for the PB potential for either the cylindrical or spherical polyion in electrolyte solution has not been obtained. The high salt concentration solution is of particular interest for comparison with the linearized (DH) approximation, still widely used for potential calculations.24 The solution obtained in this paper (eq 10a) is the first complete asymptotic solution of the nonlinear PB equation for cylindrical or spherical symmetry at large radial coordinates expressed explicitly in elementary functions. By contrast other solutions17,25-31 have only a few terms in explicit form, and the only previously published complete solution for spherical symmetry32,33 involves integral relationships between successive terms. The complete asymptotic expansion eq 10a can be generalized for arbitrary electrolyte composition. The complete asymptotic solution eq 10a provides the behavior of the exact solution of the PB equation at large enough radial distance (r/a > ) under any physical conditions. Thus, any other asymptotic or approximate solution, valid at large radial coordinates, should have matching asymptotic behavior at infinity. Therefore comparing any other approximation with the complete solution at infinity will provide information about the accuracy of this approximation and any terms that are neglected. The classical approximate solution at large radial coordinate is the Debye-Hu¨ckel (DH) solution of the linearized PB equation, which is only a part of the complete solution, as shown above. The complete solution yields an analytic expression for the constant M which multiplies the function K0(X) or e-XX-1 in the DH cylindrical or spherical solutions correspondingly. The linear part of the complete solution (first row, eq 10b) is ylnr,cyl ) A00(2/π)1/2K0(x-1) for a cylinder and ylnr,sph )
5166 J. Phys. Chem. B, Vol. 104, No. 21, 2000
Shkel et al.
Figure 1. Approximate solutions for the reduced Poisson-Boltzmann potential (absolute value), y, as a function of reduced radial coordinate, r/a, for the cylinder (A) and the sphere (B): numerical solution ynum (solid line), the DH approximation yDH (upper dashed line), and the linear part of the complete solution, ylnr, with the constant from eq 17 (lower dashed line), at ) 0.43 and σ* ) 8.4. (Corresponding parameters for the DNA model: cylinder radius a ) 10 Å, charge separation b ) 1.7 Å, salt concentration Cb ) 0.5 M, corresponding parameters for a sphere: radius a ) 10 Å, charge Q ) 11.88 e, salt concentration Cb ) 0.5 M.)
Figure 2. Comparison of approximate solutions with the numerical solution. The magnitude of deviation of the reduced potential from the numerical result |∆y| ) |y - ynum| is plotted as a function of reduced radial coordinate, r/a, for the cylinder (A) and the sphere (B) at ) 0.43 and σ* ) 8.4. The error of the Debye-Huckel approximation (yDH) is represented by the upper dashed curve, the linear part of the complete solution (ylnr) is the lower dashed curve, the planar solution is the upper solid curve and the quasi-planar approximation (yQP,0) is the lower solid curve.
A00e-x//x for a sphere, where A00 is determined by eq 15 or eq 17. The application of the inner boundary condition (eq 2) to these linearized solutions yields the following expressions for potential, which are traditionally called DH solutions: for a cylinder, yDH,cyl ) σ*K0(x/)(K1(1/))-1, and for a sphere yDH,sph ) σ*e-(x-1)/(x(1 + ))-1. Stigter34 observed that yDH,cyl overestimates the PB potential at larger distances (2-fold for the conditions in ref 34). By fitting the Bessel function to the numerical solution at large radial coordinates, Stigter found that the multiplicative constant in front of the Bessel function should be smaller in order to accurately describe the potential at larger distances from the surface. Here we obtain the analytic form of this constant (eqs 15, 17) from the complete solution. Figure 1A-B compares the DH linearized solutions (yDH defined above) with the linear part of the complete solutions (ylnr) and numerical solutions (ynum) for cylindrical (1A) and spherical cases (1B). Conditions for the figures are selected for which ynum(a) > 1: a ) 10 Å, b ) 1.7 Å, Cb ) 0.5 M, σ* ) 8.4, ) 0.43, y(r ) a) = 2 - 3, the zero order approximation for constant A00 (eq 15) is used to plot ylnr. As expected, neither linearized solution is accurate at the surface. The solution ylnr is much closer to the numerical solution than yDH for r/a > 1.5. The use of eq 15 (eq 17) for A00 provides an accuracy of ylnr,cyl better than 10% (4%) at r/a > 1.5 and ) 0.43, while the DH approximation yDH,cyl has ∼34% error at the same conditions. For the spherical case, the corresponding errors are 7% (eq 15), 2.5% (eq 17), and 21% (DH). This comparison shows that the complete solution provides a much more accurate expression for A00 than the DH approximation (which yields σ*(π/2)1/2K1(1/))-1 for the cylinder and σ*e1/(1 + )-1 for the sphere). The two linearized solutions, ylnr and yDH differ less with increasing salt concentration (y(a) < 1) and have equal limits as f 0 (very high salt concentration). For low salt, > 1, the determination of A00 is a separate challenging problem since the complete asymptotic solution eq 10 is not valid at the surface of the polyion and another representation of the exact solution is needed.
Comparison of the Quasi-Planar (QP) Approximation with Linearized (DH) Solution: Importance of Nonlinear Terms. The quasi-planar approximation for y, eqs 14 and 16, is more accurate than the linearized potential (ylnr) for calculations near the polyion surface. Since ylnr is the solution of the linearized equation, all the other terms in the complete solution result from the nonlinearity of the PB equation. At infinity, the complete solution (eq 10a) is an asymptotic series since each successive function decays faster. Since exponentials decay much faster than any power of the same variable, the entire row k in the complete solution is much smaller than the first term in row (k - 1) at large enough radial coordinates. Consequently ylnr is the most significant part of the solution. However, at small values of the coordinate, terms with different exponents become comparable and the quasi-planar approximation is the better one to use. This result is particularly important for highly charged surfaces, σ* > 1, at salt concentrations where < 1 and y(a) > 1. Figure 2A-B shows a comparison between the numerical solution of the PB equation (ynum) and various approximations for the cylindrical (2A) and spherical (2B) cases and parameters σ* ) 8.4, ) 0.43. One can see that at these conditions, the zero-order term of the quasi-planar approximation (yQP,0) is as good as the linearized approximation (ylnr) at large distances and significantly better at the surface where the relative error of the quasi-planar approximation is ∼4% and the error of ylnr is ∼18%. While the quasi-planar approximation can be improved further by using the next term, y1 (see Table 3), which is expressed in elementary functions, the approximation based on the DH solution32,33 cannot. Since the quasi-planar approximation is an asymptotic series in powers of , its accuracy becomes better as decreases. The accuracy also depends on the parameter σ*. Table 3 shows the relative error for the surface potential of linearized (ylnr) and the quasi-planar (yQP) approximation as a function of at two values of σ*. In both cases, the quasi-planar approximation is much better at the surface than ylnr for any < 1. For applications to doublestranded DNA (σ* ) 2ξ ) 8.4), the quasi-planar approximation
Complete Asymptotic Cylindrical and Spherical P-B Solutions TABLE 3: Comparison of Accuracy of Quasi-Planar (QP) and Linear Approximations to the Surface Potential for Cylindrical and Spherical Polyionsa Cb (M)b
σ* c 8.4 ycyl num δyQP,0 (%)d,e δyQP,1 (%)d,e δylnr (%)d,f ysph num δyQP,0 (%) δyQP,1 (%) δylnr (%)
2
ycyl num δyQP,0 (%) δyQP,1 (%) δylnr (%) ysph num δyQP,0 (%) δyQP,1 (%) δylnr (%)
0.1
0.3
0.5
0.7
1
1.5
2
3.861 10.4 -1 33.7 3.41 13.3 -6.7 26.3
2.89 5.9 g 22.6 2.60 7.6 -2.7 17.2
2.479 4.2 s 17.9 2.25 5.4 -1.7 13.5
2.22 3.3 s 15.0 2.03 4.3 -1.2 11.1
1.961 2.4 s 12.2 1.81 3.2 -0.8 8.9
1.687 1.6 s 9.4 1.57 2.2 -0.5 6.7
1.507 1.2 s 7.8 1.41 1.7 -0.3 5.5
1.302 6.1 2.6 16.2 0.977 2.3 -2.5 4.1
0.867 2.7 1.6 8.8 0.714 0.9 -0.6 2.0
0.706 1.8 1.2 6.9 0.603 0.6 -0.3 1.3
0.613 1.3 0.9 5.9 0.535 0.4 1.0
0.526 1.0 0.7 3.8 0.469 0.3 s 0.8
0.440 0.7 0.5 2.6 0.400 s s 0.5
0.387 0.5 0.4 1.9 0.356 s s 0.4
a Superscripts cyl and sph indicate cylindrical and spherical cases. Relative errors of the surface potential (in %) is calculated as δylnr,QP ) (1 - ylnr,QP(a)/ynum(a))100% between the numerical solution (ynum) and corresponding approximation. c The quasi-planar approximation is calculated with one term yQP,0 ) y0 (eq 14) or two terms yQP,1 ) y0 + y1 (eqs 14, 16), respectively. d The linear part of the complete solution ylnr is the first row of the complete solution eq 10b with constant A00 from eq 17. e The accuracy of numerical surface potential is 0.2%; cases where the approximate solution deviates by less than 0.2% are indicated by s. b
is better over the range 0.1 < < 1, which corresponds to 7 > Cb > 0.1 M and includes experimentally relevant salt concentrations. Other Solutions. Many approximate solutions of the PB equation have been previously reported,17,25-33 which generally are leading terms in an asymptotic expansion of the exact solution. For example, the DH solution of the linearized spherical PB equation4 was shown to be the leading term of an asymptotic expansion of the solution at large radial coordinate (at infinity),32 and the solution of the planar PB equation17,18,21,25,26 is the leading term in an asymptotic expansion in a small parameter in the vicinity of the surface of a polyion, as was very accurately demonstrated by Chew and Sen17 using the matched asymptotic expansion method for the spherical PB equation. Even an approximate solution obtained by perturbation methods may be a part of an asymptotic expansion, as we demonstrate for solutions by Ohshima et al.,27 White,28 and Ohshima29 for spherical and cylindrical geometries. For spherical symmetry, the original solution of Gronwall and La Mer32,33 is a convergent series with the DH approximation as the leading term. Each successive term is expressed through previous terms by integral relations and is not expressed in elementary functions. Nevertheless, from the consideration of the power of the exponent in each row, one may conclude that each successive term in this solution corresponds to a row in the complete solution. Thus, each nonlinear term from the complete solution (terms with higher order exponents, e-(2k+1)XX-(2k+1)j/2, where k > 0) is taken into account in the (2k + 1)th term of the Gronwall-La Mer solution, which requires a large number of terms in order to achieve a good accuracy near the surface. In contrast, the leading term in the quasi-planar approximation takes into account all the most significant nonlinear terms at once, which makes it much more accurate on the surface than the leading terms of the other
J. Phys. Chem. B, Vol. 104, No. 21, 2000 5167 approximation expressed in elementary functions. This is obvious from Figure 2A and B, where the leading terms of the two solutions are compared. This comparison illustrates that the asymptotic solution (e.g., the quasi-planar approximation) is much better for practical purposes than the convergent series.32,33 As compared with the quasi-planar (QP) approximation, the perturbation method developed by Ohshima et al.,27 White,28 and Ohshima29 produces the same zero order term (y0) for the spherical case and a similar term for the cylindrical case. The difference is that the Bessel function K0(X) replaces e-XX-1/2. This approach requires numerical calculation and does not yield further terms in explicit form. Our method yields subsequent terms in elementary functions, and shows a clear connection between the QP approximation and the complete asymptotic solution. The solution by Ohshima29 for cylindrical geometry is expressed through Bessel function (K0(X)) which makes it inapplicable to calculations of the type presented in the Thermodynamic Applications Section. Weisbuch and Gueron,25 Chew and Sen,17 and recently Rouzina and Bloomfield26 proposed planar-based asymptotic (PBA) approximations for radial potential distribution. The same approach was used by Lekkerkerker18 to obtain an approximate solution for the surface potential. Both PBA and QP approximations are asymptotic expansions in small parameter , and therefore approach each other and the exact solution as f 0 (high salt) in the region where both are valid (near the polyion surface). The QP approximation is valid at any radial coordinate (uniformly valid), while the PBA approximation is valid only in the vicinity of the surface. To construct a uniformly valid approximation from the solution of the planar PB equation, one needs to match the PBA solution to one in the region far enough from the surface and thereby obtain a composite solution,17 which is in a complicated functional form. Although both PBA and QP solutions are asymptotic expansions in powers of small parameter , they are obtained at different conditions on the second parameter; PBA is the expansion at p ) σ* ) const, QP is the expansion at p˜ ) p/(1 + j/2) ) const. The difference it produces for the accuracy of the surface potential is illustrated for leading terms of PBA18,25,26 and QP expansions in Figure 2 under the same physical conditions as for Figure 1. A comparison of the accuracy of calculating surface potential by PBA and QP expansions (eqs 43, 44 Appendix C) shows that the PBA solution requires three terms to achieve the same accuracy for surface potential (better than 1% for both double-stranded and single-stranded DNA at Cb g 0.1 M) as the two-term QP approximation. The previous section provides examples of calculation of thermodynamic properties (the electrostatic free energy, fel, and the preferential interaction coefficient, Γ) of polyelectrolyte solution based on surface potential of a polyion. The expression for the preferential interaction coefficient calculated from the two term QP approximation (ΓQP,1) is simpler in functional form and more accurate at < 1 than the corresponding expression (ΓPBA,2) calculated from the three term PBA approximation. The quasi-planar solution also clarifies the similarities and differences between PB solutions for three different geometries: planar, cylindrical, and spherical. The planar solution is considered to be a good approximation for the potential near highly charged polyion surface which is usually referred to as the planar screening regime.25,26 The zero-order term of the QP solution, eq 14, differs from the planar solution in two ways: first, there is a multiplier x-j/2 in front of the exponent function; second, the constant ˜t (or equivalently, p˜ and q˜ ) depends on the
5168 J. Phys. Chem. B, Vol. 104, No. 21, 2000 geometry. As a result, the QP solution predicts that near the polyion surface (where x ≈ 1 and x-j/2 is insignificant) the potential y indeed behaves like the planar solution but with different effective surface charge density, σ/eff ) σ*/(1 + j/2), which is lower than the actual σ*. The difference between actual and effective charge densities is larger for more curved surfaces, i.e., it is larger for a sphere (j ) 2) than for a cylinder (j ) 1). Consequently, the leading term of QP approximation, eq 14, is a better approximation than the planar solution for both cylindrical and spherical polyions. Moreover, it has comparable accuracy for the cylinder and the sphere (Figure 2) whereas the planar solution is known to have larger error for the sphere than for the cylinder. As the salt concentration increases ( f 0), the effective charge density approaches the actual surface charge density (σ*eff f σ*), the zero-order term of the QP approximation (yQP,0) behaves near the polyion surface more and more like the planar solution, and the difference between the three geometries vanishes. Taking into account the first-order term in the quasi-planar approximation (yQP,1 ) y0 + y1) significantly improves the accuracy of the approximation (see Table 3). Table 3 shows that yQP,1 has better accuracy for cylinder than for sphere, and for highly charged surfaces (σ* ) 8.4). Thus for parameters corresponding to double-stranded DNA, the surface potential is accurate to within 1% for bulk concentrations Cb > 0.1 M. Therefore our solution is accurate in the typical experimental range of salt concentrations (0.1-0.2 M). In this range, the PB potential is theoretically shown to be a good approximation of the potential of mean force35 and PB calculations of radial ion distributions and thermodynamic coefficients agree well with Monte Carlo calculations for the same model performed with a reasonable choice of ionic radii.10,36 Conclusions From the novel, complete asymptotic solution of the nonlinear cylindrical and spherical PB equation (eq 10a), we have developed a new quasi-planar approximate solution, and have compared it with previously published approximate solutions by considering terms neglected in them. Both the complete asymptotic solution and the quasi-planar approximation are expressed explicitly and in elementary functions. The quasiplanar solution is an expansion in the small parameter which is uniformly valid with respect to radial coordinate; it exhibits clear advantages compared to the DH solution (which is less accurate at the polyion surface) and to the planar-based expansion (which is valid in the vicinity of the polyion surface only and breaks down at larger r). Analysis based on the complete solution shows that for a given number of terms the quasi-planar approximation must be the most accurate approximation in the vicinity of the charged surface, because it takes into account the most significant nonlinear terms. Acknowledgment. We thank Dr. Charles Anderson for helpful discussions, Dr. Charles Woodbury and an anonymous referee for critical reading of the manuscript and comments. This work was supported by NIH grants GM34351 and GM47022. Appendix A. Asymptotic Expansions Here we outline the most relevant features of asymptotic expansions (cf. refs 9, 13-15 for more complete discussions). An asymptotic expansion of function y(X) is a series of functions yi(X) arranged in such a way that each successive function decays faster than the previous one in the vicinity of a particular
Shkel et al. point, X0, i)∞
y(X) )
yi(X) ∑ i)0
where yi+1(X)/yi(X) f 0 as X fX0, where X0 can be finite or infinite. If function y(X) is a solution of a differential equation, then its asymptotic expansion is a form of the exact solution of the equation if all terms in the expansion are known. A series that exactly satisfies the equation and boundary conditions is called a complete asymptotic expansion of the solution.9 Any other asymptotic or approximate solution, valid in the vicinity of X0, should have matching asymptotic behavior. (That is, the solution itself or its expansion should match exactly a part of the complete solution or the entire complete solution as X f X0; the phrase “to match” means here “to have the same functional form with the same constants”). Although asymptotic series may diverge, an optimal number of terms (N) can be found at a given value of independent variable (X) in an asymptotic series, where the accuracy of the first N terms of the series is the best. The accuracy of a truncated asymptotic series is determined by the remainder RN, defined as the difference between the function and its truncated asymptotic expansion i)N
RN(X) ) y(X) -
yi(X) ∑ i)0
According to asymptotic theory,15 the remainder RN(X) is on the order of the first term neglected in the truncated series. In practice, one terminates the series as soon as the desired accuracy is achieved or when yi+1 is larger that the last term (yi(X)) in the truncated series for any X of interest.15 The accuracy of the truncated asymptotic series also depends on the value of X and improves as X becomes closer to X0. Thus, the asymptotic series approximates the function as X f X0 (in contrast to convergent series, which becomes closer to the function as N f ∞). Since asymptotic theory does not impose any restrictions on the functional form of yi(X) (only on the rate of its decrease, at X f X0), a function y(X) can have more than one asymptotic expansion in the vicinity of X0. It appears to be very practical for the case considered here since one of the asymptotic expansions presented here (quasi-planar) is composed of functions which are sums of a large number of terms from another asymptotic expansion (the complete asymptotic expansion). Thus, a small loss in functional simplicity is compensated for by a very large improvement in accuracy, where only one term from the quasi-planar expansion provides the same accuracy as a large number of terms from the complete asymptotic expansion. Appendix B. The Method of Multiple Scales To obtain a uniformly valid asymptotic expansion of the exact solution of the PB equation, we apply the method of multiple scales13 to eq 1 and the boundary conditions (eqs 2-3)
r-jd(rjdy/dr)/dr ) κ2 sinh(y)
|
|
dy σ* dy ) - , rj )0 dr r)a a dr r)∞ We introduce two independent variables
(23) (24)
Complete Asymptotic Cylindrical and Spherical P-B Solutions
x1 ) (r/a - 1)/, x2 ) r/a
(25)
We seek the asymptotic solution in the form of an expansion in a series in powers of small parameter ∞
y(x1, x2, , σ*) )
jyj(x1, x2, σ*) ∑ j)0
(26)
By replacing derivatives in eq 23 using ad/dr ) -1∂/∂x1 + ∂/∂x2, and a2d2/dr2 ) -2∂2/∂x21 + 2-1∂2/∂x1∂x2 + ∂2/∂x22, we obtain the equation and the boundary conditions in terms of new variables x1 and x2, 2 ∂2y j ∂y 2j ∂y ∂ 2y 2∂ y + + 2 + + ) sinh(y) ∂x1∂x2 x1 ∂x2 x2 ∂x2 ∂x21 ∂x22
(
)|
1 ∂y ∂y + ∂x1 ∂x2
xj2
(
) -σ*
)|
(
∂x21
)|
1 ∂y0 ∂y0 + ∂x1 ∂x2
x1)0,x2)1
)0
(29)
) sinh(y0)
(
)|
1 ∂y0 ∂y0 + ∂x1 ∂x2
(
)
1 + t(x2)e-x1
1 - t(x2)e-x1
(32)
∂2y1
∂2y0 j ∂y0 y cosh(y ) ) -2 1 0 2 ∂x1∂x2 x2 ∂x1 ∂x1
(33)
At infinity, y1 should decay as fast as y0 or faster. This is possible if the nonhomogeneous part of equation
)
∂2y0 j ∂y0 -2 ∂x1∂x2 x2 ∂x1
)
(34)
1 - ˜t x2-j/2e-x1
∂x21
-2x1 -4x1 1 + 6t˜2x-j + ˜t 4x-2j 2 e 2 e
- y1
-2x1 2 (1 - ˜t 2x-j ) 2 e
)
8jt˜3x-1-3j/2 e-3x1 2 -2x1 2 (1 - ˜t 2x-j ) 2 e
(35)
( (
)| )|
1 ∂y1 ∂y1 + ∂x1 ∂x2
xj2
(31)
Equation 30 for y0 does not contain derivatives with respect to x2. It leads to a function t(x2) in place of an integration constant. Since we artificially introduced an extra variable, we cannot determine this function from only the boundary condition on the surface. An additional condition is needed. In the method of multiple scales, an additional condition is obtained from the requirement that the next term in the expansion is small everywhere compared to the previous one.13 This produces an equation for the function t(x2) and, at the same time, guarantees that the obtained expansion is uniformly valid. This condition is derived from the equation for the next term, y1:
(
∂2y1
)0 x1)∞,x2)∞
(
1 + ˜t x2-j/2e-x1
The procedure can be carried on for any number of terms. We briefly outline it here for the first-order term, y1. The equation and boundary conditions for y1 after substituting y0 are:
(30)
The equation for the zero order term (y0) is the planar PB equation with respect to x1, its solution subject to the boundary condition at infinity is
y0 ) 2 ln
which allows us to formulate the boundary condition eq 31 in terms of ∂y0/∂x1 only. Finally, the zero order term in the solution is completely determined:
y0 ) 2 ln
x1),∞,x2)∞
) -σ*, xj2
(∂y0/∂x2)|x1)0,x2)1 ) (j/2)(∂y0/∂x1)|x1)0,x2)1
(28)
After expanding sinh(y), sinh(y) ) sinh(y0) + (y1 + 2y2 + ...) cosh(y0) + (y1 + 2y2 + ...)2sinh(y0)/2 + ..., and substituting this expansion into eq 23 and boundary conditions eqs 24, we equate terms at the same powers of and obtain equations for each term in the expansion. For the zero term:
∂ 2y 0
decays faster than e-x1. The two nonhomogeneous terms behave at infinity as 2∂2y0/∂x1∂x2 ∼ -8t′(x2)e-x1 and jx-1 2 ∂y0/∂x1 ∼ -x1, where t′(x ) ) dt/dx . By requiring terms of e -4jt(x2)x-1 2 2 2 order e-x1 to cancel each other, we obtain the equation for t(x2), jt(x2) + 2x2t′(x2) ) 0, and its solution t(x2) ) ˜tx-j/2 2 . The constant ˜t now can be obtained from the boundary condition at the surface ˜t ) -p˜ -1 + (1 + p˜ 2)1/2, where p˜ ) σ*/(2 + j). In determining constant ˜t we used the identity
(27)
x1)0,x2)1
1 ∂y ∂y + ∂x1 ∂x2
J. Phys. Chem. B, Vol. 104, No. 21, 2000 5169
1 ∂y1 ∂y1 + ∂x1 ∂x2
)0
(36)
x1)0,x2)1
)0
(37)
x1)∞,x2)∞
Solution y1 is
(
-x1 jt˜x-j/2 2 e ˜t 2e-2x1 + k(x2) y1 ) -2x1 1 - ˜t 2x-j x1+j 2 e 2
)
(38)
The function k(x2) is obtained from the similar requirement that the nonhomogeneous part of the equation for y2 does not contain terms of order e-x1 at infinity, k(x2) ) c1 + (2 - j)/(2x2). The constant c1 is determined from the boundary condition at the surface eq 36 in the leading order of , c1 ) ˜t4 - ˜t2(4 - j/2) 1 + j/2)/(1 + ˜t2). Finally, the first-order solution is y ) y0 + y1 + O(2), where
(
)
-x1 jt˜x-j/2 2 e ˜t 2e-2x1 (2 - j) c + 1+j + y1 ) 2 -j -2x1 1 2x2 1 - ˜t x2 e x2
(39)
Equations 34 and 39 yield the expression for the surface potential in two term quasi-planar approximation:
yQP,1|r)a ) (y0 + y1)|r)a ) 2 ln
(11 -+ ˜˜tt) + j1 2t+˜ ˜t 3
2
(40)
Appendix C. Thermodynamic Expressions for Electrostatic Free Energy (fel) and Preferential Interaction Coefficient (Γ) from Quasi-Planar and Planar-Based Approximations Two types of notation are used for planar based asymptotic (PBA) approximation: p and q are used by Lekkerkerker18 and
5170 J. Phys. Chem. B, Vol. 104, No. 21, 2000
Shkel et al.
van Aken et al.,21 t is used by Chew and Sen17 for planar-based expansion; p˜ , q˜ , and ˜t are introduced here for quasi-planar expansion. Note that p and q (as well as p˜ and q˜ ) are not independent constants, the two are introduced for convenience (to avoid a square root in expressions).
p ) σ*/2, q ) x1 + p2, t ) p˜ )
q-1 p
(41)
p q˜ - 1 σ* ) , q˜ ) x1 + p˜ 2, ˜t ) (42) 2 + j 1 + j/2 p˜
The expression for surface potential in the quasi-planar approximation given by eq 40 in terms of ˜t, in terms of p˜ and q˜ has the form:
(q˜ - 1)2 yQP,1|r)a ) 2 ln(p˜ + q˜ ) + j p˜ q˜
(43)
For comparison, the expression for surface potential in the planar-based expansion18 in two notations is:
( )
2t(1 - t2) 1+t - j + 1-t 1 + t2 (1 - t2)2ln(1 - t2) (3 + t2)(1 - t2)3t 2 j2 + j(j -1) ) 2(1 + t2)3 t(1 + t2) (2q + 1)(q - 1)2 2(q - 1) 2 ln(p + q) - j + 2 j2 pq p3q3
yPBA,2|r)a ) 2 ln
(
)
(
)
2 ln((q + 1)/2) (44) j(j - 1) pq The electrostatic free energy is calculated from the two term QP planar potential, eq 43, by integrating it with respect to surface charge density, and consequently, has two terms in powers of
[
fel,QP,1 ) 2RTp˜ -1 2p˜ ln(p˜ + q˜ ) - 2q˜ + 2 +
(
j q˜ - 1 - 2 ln
(q˜ +2 1))] ) 2RT[2 ln11 -+ ˜˜tt - 2t˜ +
(
j ˜t +
1 - ˜t 2 ln(1 - ˜t 2) ˜t
)]
(45)
In derivation of the expression for Γ, we used the electroneutrality condition, Γ+coul - Γ-coul ) 1, which follows from eq 18 and yields -2Γ- ) 1 + ∂fel/∂µ (eq 19). The last term in this expression involves the derivative of fel with respect to salt concentration. Since salt concentration enters in both p˜ and , the resulting expression for ∂fel/∂µ does not appear as a series in powers of . Simply expanding ∂fel/∂µ in powers of and truncating it at two terms introduces an error which can be partly avoided. Instead, we expand the product p(∂fel/∂µ) in powers of , which allows us to preserve the zero order term in ∂fel/∂µ and neglect only the 2 order contribution to the first-order term. Finally, the expression for Γ is
(q˜ - 1) 1 1 - q˜ -ΓQP,1 ) + - j . 2 2p 4pq˜ 2
(46)
The three term expression for the preferential interaction coefficient ΓPBA,2 is obtained in ref 21 from fel calculated from the three term planar-based expansion of the surface potential (eq 44).
((
1 - q 2 j2 1 1-q 2 + - j + -1 -ΓPBA,2 ) + 2 2p 2pq 2p 4 q(q + 1) j(j - 1) 2 1 + q z)1 ln z 2 ln dz (47) + + z)2/(q+1) 3 2 q 2 z -1 q
)
( (
)
∫
))
References and Notes (1) Verwey, E. J. W.; Overbeek, J. T. G. Theory of Stability of Lyophobic Colloids; Elsevier: Armsterdam, 1948. (2) Fuoss, R. M.; Katchalsky, A.; Lifson, S. Proc. Natl. Acad. Sci. U.S.A. 1951, 37, 579-589. (3) Benham, C. J. J. Chem. Phys. 1983, 79(4), 1969-1973. (4) Debye, P.; Huckel, E. Phys. Zeit. 1923, 24, 185-206. (5) Hecht, J. L.; Honig, B.; Shin, Y.-K.; Hubbell, W. L. J. Phys. Chem. 1995, 99, 7782-7786. (6) Manning, G. S. J. Chem. Phys. 1969, 51, 924-933. (7) Manning, G. S. Q. ReV. Biophys. 1978, 11, 179-246. (8) Anderson, C. F.; Record, M. T., Jr. In Structure and Dynamics: Nucleic Acids and Proteins; Clementi, E., Sarma, R., Eds.; Adenine: New York, 1983; pp 301-319. (9) Dingle, R. B. Asymptotic Expansions: Their DeriVation and Interpretation; Academic Press: London, 1973. (10) Ni, H.; Anderson, C. F.; Record, M. T. Jr. J. Phys. Chem. 1999, 103, 3489-3504. (11) Ennis, J.; Marcelja, S.; Kjellander, R. Electrochem. Acta 1996, 41(14), 2115-2124. (12) Shklovskii, B. I. Phys. ReV. E 1999, 60(5), 5802-5811. (13) Nayfeh, A. Perturbation Methods; John Wiley & Sons: New York, 1973. (14) Van Dyke, M. Perturbation Methods in Fluid Dynamics; Academic Press: New York, 1964. (15) Murray, J. D. Asymptotic Analysis; Clarendon Press: Oxford, U.K., 1974. (16) Hill, T. L. Arch. Biochem. Biophys. 1955, 57, 229-239. (17) Chew, W. C.; Sen, P. N. J. Chem. Phys. 1982, 77(4), 2042-2044. (18) Lekkerkerker, H. N. W. Physica A 1989, 159, 319-328. (19) Record, M. T., Jr.; Zhang, W.; Anderson, C. F. AdVances in Protein Chemistry 1998, 51, 281-353. (20) Anderson, C. F.; Record, M. T. Jr. Biophys. Chem. 1980, 11, 353360. (21) van Aken, G. A.; Lekkerkerker, H. N. W.; Overbeek, J. T. G.; de Bruyn, P. L. J. Phys. Chem. 1990, 94, 8468-8472. (22) Sharp, K. A. Biopolymers 1995, 36, 227-243. (23) Strauss, U. P.; Helfgott, C.; Pink, H. J. Chem. Phys. 1967, 71, 2550-2556. (24) Fogolari, F.; Zuccato, P.; Esposito, G.; Viglino, P. Biophys. J. 1999, 76, 1-16. (25) Weisbuch, G.; Gueron, M. J. Phys. Chem. 1981, 85, 517-525. (26) Rouzina, I.; Bloomfield, V. A. J. Phys. Chem. 1996, 100, 42924304. (27) Ohshima, H.; Healy, T. W.; White, L. R. J. Colloid Interface Sci. 1982, 90(1) 17-26. (28) White, L. R. J. Chem. Soc., Faraday Trans. 2 1977, 73, 577-596. (29) Ohshima, H. J. Colloid Interface Sci. 1998, 200, 291-297. (30) Dukhin, S. S.; Semenichin, N. M.; Shapinskaya, L. M. Dokl. Akad. Nauk USSR 1970, 193, 385-388. (31) Parlange, J.-Y. J. Chem. Phys. 1972, 57, 376-377. (32) Gronwall, T. H.; LaMer, V. K.; Sandved, R. Physik. Z. 1928, 29, 358-393. (33) La Mer, V. K.; Gronwall, T. H.; Greiff, L. J. J. Phys. Chem. 1931, 35, 2245-2288. (34) Stigter, D. J. Colloid Interface Sci. 1975, 53(2), 296-306. (35) Fixman, M. J. Chem. Phys. 1979, 70, 4995-5005. (36) Le Bret, M.; Zimm, B. H. Biopolymers 1984, 23, 271-285.
