Influence of Mixed and Multivalent Counterions in Overcharging of

It is seen that the depth of the minimum in the excess ground-state energy (over the neutral reference state) versus the number of overcharging counte...
0 downloads 0 Views 155KB Size
Download PDF
11802

Langmuir 2004, 20, 11802-11810

Influence of Mixed and Multivalent Counterions in Overcharging of DNA-like Spherocylindrical Macroions A. K. Mukherjee Department of Physics, University of Puerto Rico, San Juan, Puerto Rico 00931-3343

K. S. Schmitz Department of Chemistry, University of MissourisKansas City, Kansas City, Missouri 64110

L. B. Bhuiyan* Department of Physics, University of Puerto Rico, San Juan, Puerto Rico 00931-3343 Received March 1, 2004. In Final Form: October 1, 2004 The influence of multivalent and mixed valency counterions on the ground-state energetics of overscreening of a core DNA-like model (sphero)cylindrical macroion is investigated using an earlier developed energy minimization numerical simulation algorithm. The effects of mono-, di-, tri-, and tetravalent counterions, and mixed valency (mono- plus di-) counterions are compared and contrasted. It is seen that the depth of the minimum in the excess ground-state energy (over the neutral reference state) versus the number of overcharging counterions increases as counterion valency changes from mono- to tetra- testifying to the efficiency of the overcharging process due to multivalent counterions. The influence of (i) the presence of mixed valency counterions and (ii) counterion size on the energetics is also investigated.

1. Introduction Short pieces of deoxyribonucleic acid (DNA) physically behave as rigid cylinders and, therefore, may serve as a model experimental system against which theories and simulations for uniformly charged cylinders of finite length and thickness may be compared. This includes not only local properties such as counterion distribution about the molecule but also long range phenomena such as bundling and gel formation. Highly monodispersed preparations of cylindrical DNA can be obtained by nuclease digestion of chromatin. DNA wraps around an octameric complex of histones to form a “bead-and-string” chromatin structure as the first order of folding the genetic information in eukaryotic cells. Nuclease digestion of chromatin material and subsequent extraction of the histones result in a highly monodisperse sample of core DNA, which for our purpose of simplifying the notation we denote as “cDNA”. The physical characteristics of the cDNA is a rigid cylinder of diameter of about 20 Å and a length of 150 × 3.4 Å ) 510 Å, or about one persistence length. One of the interesting properties of cDNA is that of an ionic strength induced “phase separation” into coexisting “free” and “gel” domains. We note the dynamic light scattering (DLS) study of Fulmer, Benbasat, and Bloomfield1 on cDNA. The DLS method monitors the fluctuations in the intensity of scattered light and represents this information as an autocorrelation function denoted as C(t). In the case of cDNA, the data of Fulmer et al. indicate that C(t) is a single-exponential decay at ionic strengths above 10-2 M but exhibits two decay modes at lower ionic strengths. They interpret the latter result as a manifestation of two environments for their cDNA: those free to diffuse in the solution to account for the fast decay and those in a “gel” to account for the slow mode. When asked if there was any other “evidence” that gels were formed, (1) Fulmer, A. W.; Benbasat, J. A.; Bloomfield, V. A. Biopolymers 1981, 20, 1147.

Bloomfield gave the response that when the sample cell was inverted none of the solution flowed from the cell (V. A. Bloomfield, personal communication). It was also reported that DNA extracted from dinucleosome fractions also exhibits two decay rates, where the fast mode was consistent with that of a rigid cylinder.2 The two phases inferred from the DLS experiments are not macroscopic phases with well-defined boundaries but are what might be called “microphases” of freely diffusing cDNA among “bundled” DNA regions. Such a “microphase separation” in the decay rates of the cDNA in low ionic strength likely suggests some form of “attraction” between the particles. There are several studies in the literature in which an apparent “attraction” occurs between macroions of like charge. The roots of many of these studies are found in the pioneering papers of Kirkwood and Shumaker3,4 for globular proteins (spherical symmetry) and of Oosawa5 for rodlike particles. Kirkwood and Shumaker showed that fluctuations in the net charge of the charged spheres led to a negative chemical potential for interacting spheres if the net charge on the proteins is small and, hence, the inference of possible charge reversal as a contributor to the attraction. Oosawa likewise examined the effect of correlations in the charge density of two interacting rods, 〈δF1δF2〉. If this term is negative, then an attraction between two rods was found if the separation distance was of the order of the Bjerrum length. Such “anticorrelations” in charge density arise from redistribution of counterions on charged surfaces, as was shown later by Rouzina and Bloomfield6 for two flat surfaces, or from the relative orientation of rods with discrete charge distribu(2) Schmitz, K. S.; Lu, M. Biopolymers 1984, 23, 797. (3) Kirkwood, J. G.; Shumaker, J. B. Proc. Natl. Acad. Sci. U.S.A. 1952, 38, 855. (4) Kirkwood, J. G.; Shumaker, J. B. Proc. Natl. Acad. Sci. U.S.A. 1952, 38, 863. (5) Oosawa, F. Biopolymers 1968, 6, 1633. (6) Rouzina, I.; Bloomfield, V. A. J. Phys. Chem. 1996, 100, 9977.

10.1021/la049468w CCC: $27.50 © 2004 American Chemical Society Published on Web 11/24/2004

Overcharging of Spherocylindrical Macroions

tions, as shown by Allahyarov and Lo¨wen7 for two rods modeled with DNA discrete charge distributions. Ha and Liu8 showed that nonlinear effects in the interaction potential may likewise stabilize bundles of DNA rods. The idea of attraction via charge correlations is encountered in recent computer simulation data in which an attraction induced by such correlations can overcome the repulsion between two particles of like charge. Calculations show this occurs for spheres,9-11 parallel rods,12,13 planar surfaces,12 and ordered arrays of cylinders.14-17 We return briefly to the model of Kirkwood and Shumaker3,4 in which the pair potential is represented as a series expansion. The term that gives rise to the attraction has the square of the charge fluctuation. This leaves open the possibility that one, or both, of the macroions may exhibit a reversal in net charge. That is, more counterions become attached to the particle’s surface (or its vicinity) than is needed to compensate for the surface charge of the particle. Indeed, the possibility of attraction between a pair of like charged macroions through surface charge fluctuations leading to a charge reversal of one of the pair has also been suggested in the numerical studies by Messina et al.18,19 We ought also to cite here two recent, relevant experiments, which speak directly to the related concept of an effective macroion charge in a charge stabilized solution. In one, Ohshima et al.20 performed small-angle X-ray scattering on solutions of dendrimers with mono- and divalent counterions. They found that the counterions to be associated with the dendrimer ions (divalent ions more strongly so than monovalent ions) leading to counterion mediated attraction whose intensity depended on effective charge density and effective charge number. In another, Zhang et al.21 did light and neutron scattering experiments on polystyrene sulfonates again in the presence of mono- and divalent counterions. Two important observations were (i) an enhancement of the interchain correlation length within the polyelectrolyte domains for divalent counterions relative to monovalelent ones and (ii) a reduction in the size of these domains also for divalent ions. Both of these findings were attributed to a reduction in effective polymer charge density for divalent counterions. We should emphasize also that charge reversal itself is an experimental phenomenon having been observed in electrophoeretic mobility experiments over 5 decades ago.22,23 It is the phenomenon of overcharging and its properties that is the object of the

Langmuir, Vol. 20, No. 26, 2004 11803

present study. The many faceted significance of overcharging or charge inversion/reversal, some of which have been outlined above, has spawned a substantial interest in the subject.18,19,24-26 We also refer the reader to a couple of excellent recent reviews.27,28 Even though there are ample studies on the phenomenon of “attraction between particles of like charge”, many of their interpretations may be related to charge fluctuations, whether it is the local counterion distributions at fixed “effective charge” or a fluctuation in the net charge itself. It is thus necessary to understand how the counterions are distributed in an equilibrium configuration on the isolated surface in order to gain knowledge as to how fluctuations in these distributions might evolve to a dominant attraction force between similar macroions. Furthermore, it is of relevance to characterize such equilibrium distributions that include stable overcharged states in terms of physical parameters such as macroion geometry, counterion valency, and relative macroion/ counterion sizes. The present paper reports such a study tailored to the cDNA system. We reiterate that in this work we will confine ourselves to these aspects of overcharging in relation to cDNA and not delve into likecharge attraction at this time. In an earlier paper,29 henceforth to be denoted by I, we considered overcharging in isolated macroions of different shapes and charge distributions using simulation techniques. The ground state energetics was found to be influenced substantially by the macroion geometry and/ or charge distribution. One of the geometries studied was the spherocylinder with a uniform line charge distribution along the cylinder axis. This spherocylindrical macroion model is seemingly a natural and an interesting candidate to mimic a cDNA cylinder. In this paper we study further aspects of overcharging such as overcharging by mixtures of counterions with different valencies by treating such a model cDNA cylinder by (i) mono-, di-, tri-, and tetravalent counterions, and (ii) mixed valency counterions. The effect of counterion size on overcharging is also investigated. The paper is organized as follows. In the next section details of the model and the simulation techniques are given. In section 3 the results for a cylinder of cDNA dimensions are discussed, and in section 4 some concluding remarks are made. 2. Model and Methods

(7) Allahyarov, E.; Lo¨wen, H. Phys. Rev. E 2000, 62, 5542. (8) Ha, B. J.; Liu, A. J. Phys. Rev. Lett. 1998, 81, 1011. (9) Grøenbech-Jensen, N.; Beardmore, K. M.; Pincus, P. Physica A 1998, 261, 74. (10) Linse, P. J. Chem. Phys. 2000, 113, 4359. (11) Linse, P.; Lobaskin, V. Phys. Rev. Lett. 1999, 83, 4208. (12) Guldbrand, L.; Jo¨nsson, B.; Wennerstro¨m, H.; Linse, P. J. Chem. Phys. 1984, 80, 2221. (13) Deserno, M.; Arnold, A.; Holm, C. Macromolecules 2003, 36, 249. (14) Guldbrand, L.; Nilsson, L. G.; Nordenskio¨ld, L. J. Chem. Phys. 1986, 85, 6686. (15) Nilsson, L. G.; Guldbrand, L.; Nordenskio¨ld, L. Mol. Phys. 1991, 72, 177. (16) Lyubarstev, A. P.; Tang, J. X.; Janmey, P. A.; Nordenskio¨ld, L. Phys. Rev. Lett. 1998, 81, 5465. (17) Stevens, M. J. Phys. Rev. Lett. 1999, 82, 101. (18) Messina, R.; Holm, C.; Kremer, K. Phys. Rev. Lett. 2000, 85, 872. (19) Messina, R.; Holm, C.; Kremer, K. Eur. Phys. J. E 2001, 4, 363. (20) Ohshima, A.; Konishi, T.; Yamanaka, J.; Ise, N. Phys. Rev. E 2001, 64, 051808. (21) Zhang, Y.; Douglas, J. F.; Ermi, B. D.; Amis, E. J. Chem. Phys. 2001, 114, 3299. (22) Bungenberg de Jong, M. G. In Colloid Science; Kruyt, H. R., Ed.; Elsevier: New York 1949; p 335. (23) Strauss, U. P.; Gershfeld, N. L.; Spiera, H. J. Am. Chem. Soc. 1954, 76, 5911.

The spherocylindrical macroion model and the simulation methods employed in this work have been chronicled at length in I. We will therefore be brief in our discussion here and will restrict ourselves to outlining the principal features of the model and simulations apropos of cDNA. For convenience we will also follow closely the notations used in I. The spherocylinder consists of a rigid cylinder of length L and radius rm with two hemispherical end-caps of the same radius. As commonly used for most DNA models, the total macroion charge -Zm|e| (Zm being the macroion valency and |e| the magnitude of the electron charge) is distributed uniformly along the axis length L. The (24) Shlovskii, B. I. Phys. Rev. Lett. 1999, 82, 3268. (25) Shlovskii, B. I. Phys. Rev. E 1999, 60, 5802. (26) Messina, R.; Holm, C.; Kremer, K. Phys. Rev. E 2001, 64, 021405. (27) Grosberg, A. Y.; Nguyen, T. T.; Shlovskii, B. I. Rev. Mod. Phys. 2002, 74, 329. (28) Quesada-Perez, M.; Gonza´lez-Tovar, E.; Martin-Molina, A.; Losada-Cassou, M.; Hidalgo-A Ä lvarez, R. Chem Phys Chem 2003, 4, 234. (29) Mukherjee, A. K.; Schmitz, K. S.; Bhuiyan, L. B. Langmuir 2003, 19, 9600.

11804

Langmuir, Vol. 20, No. 26, 2004

Mukherjee et al.

(

)

(((L/2) - zi)

ψ1i ) tan-1

r cos(θi)

(

)

(L/2) + zi

ψ2i ) tan-1

r cos(θi)

(2)

(3)

and

(

θi ) sin-1

Figure 1. Geometry of a spherocylinder.

spherocylinder is surrounded by Nc + n small, simple, and rigid counterions of common radius rc residing on its surface, where Nc is the number of “neutralizing” counterions and n is the number of overcharging counterions. Thus the neutral state of our model system is specified by Zm|e| ) -Ncq, where q ) Zc|e| is the charge of a counterion of valency Zc. In this work we have allowed mixed valency counterions to be present such that the valencies of counterions comprising Nc and n may be different. The geometry of a spherocylinder can be seen in Figure 1. As mentioned in the Introduction, we have chosen the central cylinder of cDNA to be of the length L ) 510 × 10-10 m and consisting of 150 base pairs to give the length per charge unit b ) 1.7033 × 10-10 m. The radius of the spherocylinder is rm ) 10 × 10-10 m, while the counterion radius is taken to be rc ) 1 × 10-10, 1.8 × 10-10, or 2.5 × 10-10 m. As in I, we have taken the relative permittivities of the solvent outside the cylinder and of the interiors of both the spherocylinder and the counterions to be the same so that no dielectric discontinuity occurs. The value of r ) 16 used in the present calculations is the same as that used in I and by Messina et al.26 in their molecular dynamics (MD) simulation study of an overcharged spherical macroion, although admittedly, the exact value of r is not relevant for our purposes. We would like to remark here that the hemispherical end-caps of the spherocylinder have probably little to do with the cylinder of cDNA. However, these are artifacts of the model (cf. ref I) and from past experience have little bearing on the energetics of a long spherocylinder (L . rm). Following Deserno et al.,13,30 the model used here might also be dubbed a DNA-like system. The analytic expressions for the Coulomb interaction energy of the system of counterions and a spherocylindrical macroion with a line charge have been given in eq 7 of I, which we reproduce here for completeness

E ) λe where

ZmZc L

N

(

ln ∑ i)1

)

{1 - sin (ψ2i)} cos(ψ1i) {1 + sin(ψ1i)} cos(ψ2i)

(1)

)

zi - (L/2) r

(4)

where λe ) e2/(4π0r) is the Bjerrum length, o is the vacuum permittivity, and the angles ψ1i and ψ2i describe the angles of the ith bound counterion with z-coordinate zi as shown in Figure 1. Equation 1, the positive sign in eq 2, and θi ) 0 apply to the situation when the ion is on the cylindrical part of the spherocylinder, while the negative sign in eq 2 is taken when the ion is on a hemispherical end-cap. The simulation strategy is based upon ideas of energy minimization (EM). Since we are interested in the ground state (T ) 0) energy configurations, the standard, traditional Metropolis algorithm is precluded. Instead we rely on an EM procedure whereby a random “move” of a randomly chosen particle is accepted if and only if it leads to a lowering of the total energy. The ground state energetics is achieved in a two-stage process. In the first stage the ground state energy of a neutral cDNA cylinder plus counterions is calculated. The requisite number of neutralizing counterions is distributed randomly on the surface of the spherocylinder and the total electrostatic energy calculated. A counterion is then chosen at random, which is then moved randomly along the surface through a preset “displacement parameter” with the “move” being accepted if the recalculated total energy is lower, and so on. On the average, passes (number of “moves” per particle”) of the order 106 were made before the ground state configuration was achieved. An important feature of the algorithm is the displacement parameter, which adjusts itself continuously as the run proceeds. This helps in avoiding the system being trapped in metastable states (local minima). The second stage consists of adding “overcharging” counterions, one at a time, to the system and repeating the whole process to determine the ground state following each addition. For further details of the energy minimization process we refer the reader to ref I. 3. Results and Discussion We first show in Figure 2 the ground state configuration of a neutral cylinder for mono-, di-, tri-, and tetravalent counterions. The mutual repulsive interactions between the “bound” counterions result in their being symmetrically distributed over the cylindrical surface of the cylinder the symmetry being about the waist of the cylinder. Note that the hemispherical end-caps are devoid of any ion. This is due to these end-caps being energetically unfavorable locations. The configurations in Figure 2 are our reference states for the overcharging simulations and, thus, define our “zero of energy”. Starting from the neutral cases in Figure 2 the “overscreened” states are built up by the addition of like counterions onto the macroion surface. Shown in Figure 3 is the reduced “minimized excess energy”, ∆En/kBT, corresponding to each of the four cases of Figure 2 as a function of the number n of “overcharging” counterions. (30) Deserno, C.; Holm, C.; May, S. Macromolecules 2000, 33, 199.

Overcharging of Spherocylindrical Macroions

Figure 2. Ground (neutral) state energy configurations of (a) 300 monovalent, (b) 150 divalent, (c) 100 trivalent, and (d) 75 tetravalent counterions on a model spherocylindrical cDNAlike macroion (Zm ) 300). The filled circles represent counterions on the front side of the cylinder, while the open circles represent counterions on the backside of the cylinder. Only half of a cylinder is shown as the other half has the same configuration by symmetry.

Figure 3. Excess ground-state electrostatic energy (over the neutral reference state) of a model spherocylindrical cDNAlike macroion (Zm ) 300) as a function of overcharging counterions. The curves (from top to bottom) correspond to monodi-, tri-, and tetravalent counterions. The solid lines are from the analytical modified Scatchard model.

These data indicate that overcharging is favorable over the range examined. The characteristics of these plots are similar to those seen earlier in the overcharging plots in I and in the work by Messina et al.26 However, the

Langmuir, Vol. 20, No. 26, 2004 11805

Figure 4. Ground (neutral) state energy configurations of mixed valency counterions on a model spherocylindrical cDNAlike macroion (Zm ) 300): (a) 150 mono- and 75 divalent counterions; (b) 60 mono- and 120 divalent counterions. Circles represent monovalent counterions, while squares represent divalent counterions. The meaning of filled and open symbols is as in Figure 2.

depth and location of the minimum in the profile are dependent upon the valency of the counterions. The fact that the depth of the minimum decreases with higher valency counterions indicates that overcharging becomes more efficient as counterion valency increases. For some time now there has been active research, on both theoretical and experimental fronts, on the relative adsorption by DNA of counterions which differ in their valencies. The interest lies in the observation that multivalent counterions can cause dramatic conformational changes in the structure of DNA. It is of interest to investigate the competing nature of multivalent counterions with their monovalent counterparts. We therefore examined the overcharging of our model cylinder in the simultaneous presence of mono- and divalent counterions. Shown in Figure 4a is the neutral reference case in which 150 mono- and 75 divalent counterions are employed, while Figure 4b shows the case with 60 mono- and 120 divalent counterions. Generally more divalent ions are located around the vicinity of the waist of the cylinder, whereas the monovalent ions tend to be more spread out. What is more interesting is the fact that as the number of monovalent ions decreases in going from Figure 4a to Figure 4b, some of their places near the center are occupied by the divalent ions. Note, for example, that near the bottom of Figure 4b there are no monovalent ions. The explanation may be linked to the differences in the energy gain due to a monovalent counterion versus that due to a divalent counterion. It is certainly more favorable energetically for a divalent ion rather than a monovalent ion to be located near the waist of the cylinder. We next examined the relative stabilities of counterions of mixed valencies. This is shown in the excess energy plots of Figure 5. For each neutral situation depicted in Figure 4, overcharged ground state is achieved by adding

11806

Langmuir, Vol. 20, No. 26, 2004

Figure 5. Execess ground-state electrostatic energy (over the neutral reference state) of a model spherocylindrical cDNAlike macroion (Zm ) 300) as a function of overcharging counterions (in the presence of mixed valency counterions): filled diamond, 150 mono- and 75 divalent counterions with overcharging by monovalent counterions; open diamond, 150 mono- and 75 divalent counterions with overcharging by divalent counterions; filled circle, 60 mono- and 120 divalent counterions with overcharging by monovalent counterions; open circle, 60 mono- and 120 divalent counterions with overcharging by divalent counterions. The solid lines are from the analytical modified Scatchard model.

Figure 6. Ground (neutral) state energy configurations of 150 divalent counterions of different radii on a model spherocylindrical cDNA-like macroion (Zm ) 300): (a) rc ) 1 × 10-10 m; (b) rc ) 1.8 × 10-10 m; (c) rc ) 2.5 × 10-10 m. The meaning of the symbols is as in Figure 2.

either mono- or divalent counterions. Once again consistent with earlier observations we notice that the curves corresponding to the addition of divalent counterions lie below those that correspond to the addition of monovalent ions. Figures 6 and 7 pertain to the influence of counterion size on energetics of overcharging. In Figure 6 three neutral state distributions of 150 divalent counterions are shown for rc ) 1 × 10-10, 1.8 × 10-10, and 2.5 × 10-10 m, respectively, while in Figure 7 the corresponding ground-state energy versus the number of overcharging

Mukherjee et al.

Figure 7. Excess ground-state electrostatic energy (over the neutral reference state) of a model spherocylindrical cDNAlike macroion (Zm ) 300) as a function of overcharging counterions of different radii. The curves correspond to (a) rc ) 1 × 10-10, (b) 1.8 × 10-10, and (c) 2.5 × 10-10 m, respectively. The solid lines are from the analytical modified Scatchard model.

counterions plots are given. It is worth noting here that the relative distance between the axial macroion charge and a counterion center is smallest for rc ) 1 × 10-10 m and largest for rc ) 2.5 × 10-10 m. This leads to the curve for rc ) 1 × 10-10 m lying below that for rc ) 1.8 × 10-10 m and so on. The physical features of “overcharging” are quite simple and straightforward (see for example, Messina et al.26 and Grosberg et al.27). Essentially, the finite dimensions of the macroion provide an “envelope” about a region in space that counterions cannot penetrate. The shape and dimensions of this envelope define the macroion surface. The counterions, of finite size, will seek the lowest energy configurations, which in turn shape the distribution pattern of surface-associated counterions. The first counterion will be found at the location of minimal separation distance between the counterion and the macroion consistent with the maximum in the potential field set up by the macroion charge distribution. Of course the second counterion will seek to do the same, but now the counterion-counterion interaction comes into play and the two counterions adjust their positions accordingly. When a neutralizing number of counterions are placed on the macroion surface, there still remains “energy gaps” for the placement of additional counterions that leads to a negative energy contribution. The total energy continues to decrease until upon addition of a counterion the repulsive interaction energy associated with the counterions becomes greater than the net attractive interaction of the macroion-counterion charge complex. Although the total energy still remains negative, the energy curve ∆E/kBT begins an upturn with further addition of counterions with the relative size and charge of the counterions determining the shape, depth, and minimum location in the energy ∆E/kBT verses n profile. The sphere and spherocylinder are topologically equivalent structures. That is, one may stretch a sphere along its north-south axis to obtain the shape of a cylinder with hemispherical caps, viz., a spherocylinder. The stretching process also changes the point charge at the center of the sphere to a line charge along the coordinate axis of the cylinder. The energy minimization simulations

Overcharging of Spherocylindrical Macroions

reveal some interesting features of the distribution of counterions as one deforms the sphere to the topological equivalent spherocylinder. As shown in Figure 1 of I, the counterions on the surface of a sphere are arranged in regions of hexagon and pentagon shapes with counterions at the vertexes and at the center of each structure. The hexagon and pentagon shapes are thus composed, respectively, of two types of triangles as discerned by connecting the perimeter counterions with the central counterion of each structure. When the sphere is deformed to the spherocylinder, the counterions appear to arrange in a regular helical array, as shown in Figures 8 and 9. These figures show, respectively, the magnified view of sections of cylinders corresponding to the neutral state configurations in Figures 2 and 6. Since the counterion density is greater at the waist of a cylinder and decreases as one moves toward the ends, the apparent pitch of the helical array varies with distance. We now slice the spherocylinder surface along the cylinder axis (the northsouth line of the sphere) and open it to form a planar surface without disturbing the locations of the counterions. We will again see a distribution of triangles on this new planar geometry, but now the triangles are of variable size. Because of the original cylindrical symmetry, the triangles are of comparable size in the east-west direction but vary in the north-south direction. The minimum triangle size is along the equatorial line and becomes progressively larger as one travels north or south. Compare this distribution with that anticipated for counterions distributed on a uniformly charged planar surface of finite dimensions, where the counterion density will increase as one moves from the center of the plane to the perimeter due to the cumulative repulsion between the counterions. Again the distribution of counterions can be described as triangles, with the larger triangles now found in the central region of the plane. This inversion of the location of triangle size is due to the fact that the spherocylinder and planar surfaces are not topologically equivalent and the charge distribution of the surfaces differs. In contrast, in an infinite planar structure the array of counterions is infinite from the vantage point of any counterion, so that only a single triangle dimensions is required to describe the counterion distribution. The essential point that emerges from this discussion is that energy minimization simulations will distribute the counterions in triangular structures regardless of the shape of the surface. The minimum energy simulations in this paper are for the ground state, which means that T ) 0. Since the experiments of Fulmer et al.1 were performed under room temperature conditions, one might be curious to ask if “overcharging” might conceivably play a role in the interpretation of the experiments. What is necessary is to have the “effective charge” of these cylinders be greatly reduced to near neutrality while within the cluster. Without pressing this point further we would like to remark that there is some evidence in the literature31,32 that for a cluster of spherical macroions the central sphere can essentially be neutralized by the excess of counterions drawn to the interior of the cluster by the combined electric fields of the component macroions. The possible role of “overcharging” on the stability of “bundles” of DNA will be discussed in another paper. 3.1. Analytical Fit to the Simulation Results. It is of value to quantify the EM energetics in Figures 3, 5, and (31) Schmitz, K. S. Langmuir 1999, 15, 4093. (32) Schmitz, K. S. In Handbook of Polyelectrolytes and Their Applications; Tripathy, S. K., Kumar, J., Nalwa, H. S., Eds.; American Scientific Publishers: Stevenson Ranch, CA, 2002; Vol. 3, p 195.

Langmuir, Vol. 20, No. 26, 2004 11807

Figure 8. Magnified view of sections of cylinders corresponding to the ground (neutral) state configurations in Figure 2 showing helical patterns formed by the counterions: (a) 300 monovalent counterions, (b) 150 divalent counterions, (c) 100 trivalent counterions, and (d) 75 tetravalent counterions. The spherocylindrical cDNA-like macroion has the charge -300|e|. The meaning of the symbols is as in Figure 2.

7 on the basis of an equation for the purposes of comparison with other models and data on the phenomenon of overcharging. We therefore adopt the modified Scatchard model of I. The Scatchard model for binding ligands to biopolymers33 is a combinatorial model which takes into consideration the interaction of the ligand with the binding

11808

Langmuir, Vol. 20, No. 26, 2004

Mukherjee et al.

distance dependence of the respective electrostatic intereactions. To emphasize what these parameters mean in terms of their approximate methods, we will point out how their values can be obtained through simulations and how these values may differ for spherical and cylindrical geometries. In calculating 〈A〉, one simply moves a counterion to all locations on the surface and then calculates the counterion-surface interaction for each location and then an average is taken over all locations. In the case of a spherical macroion all surface locations are identical; hence 〈A〉 is the unscreened Coulomb interaction between two points separated by the sum total of their radii. In the case of a cylinder with a central line charge, the counterion-surface interaction is greatest at the waist of the cylinder. Hence the value of 〈A〉 for the cylinder depends not only on the radii of the cylinder and the counterion but also on the length of the cylinder. Similarly 〈B〉 is averaged over all possible locations of a pair of particles in the absence of all other particles. As only the bare interaction between two particles is involved in evaluation of 〈A〉 and 〈B〉, the interionic correlations are neglected within the context of the model thus imparting a mean-field character to the model. As shown in I the average energy 〈EN〉 for N ) Nc + n bound counterions for the combinatorial model can be expressed, upon separation of the Nc and n contributions, as a quadratic equation in the overcharging counterions n

〈EN〉 ) g0 + g1n + g2n2

(5)

where

{

g0 ) -ZmZc〈A〉Nc + Zc2〈B〉

}

Nc(Nc - 1) 2

(6)

represents the total interaction of only the neutralizing counterions Nc and the other two coefficients are

g1 )

1 2 Z {2Nc(〈B〉 - 〈A〉) - 〈B〉} 2 c

(7)

1 2 Z 〈B〉 2 c

(8)

and

Figure 9. Magnified view of sections of cylinders corresponding to the ground (neutral) state configurations in Figure 6 showing helical patterns formed by the counterions (a) rc ) 1 × 10-10, (b) rc ) 1.8 × 10-10, and (c) rc ) 2.5 × 10-10 m. The spherocylindrical cDNA-like macroion has the charge -300|e|. The meaning of the symbols is as in Figure 2.

surface and the average pairwise interaction between the bound ligands. It is a combinatorial model because the average pairwise interaction is assumed to be independent of the relative locations of the two interacting ligands. It has been shown34 for the exact linear model that the combinatorial approach is valid if the pairwise interaction is small relative to the thermal energy, but fails for strong interactions. In the present case the average interaction of the isolated univalent counterions (ligand) with the univalent surface is denoted by 〈A〉, and that for the average pairwise interaction between two univalent charged counterions is denoted by 〈B〉. That is, the parameters 〈A〉 and 〈B〉 represent only the reciprocal (33) Scatchard, G. Ann. N. Y. Acad. Sci. 1949, 660. (34) Schmitz, K. S. Biopolymers 1977, 16, 143.

g2 )

Note that if n ) 0 we have

g0 ) 〈ENc〉

(9)

The parameter g0 therefore represents the energy of the neutralized macroion-counterion complex relative to an overcharged such complex and is thus the underlying “reference state” for our model. Therefore the difference 〈EN〉 - 〈ENc〉 gives the energy of the system with the neutralized surface as the reference state, viz.

∆En ) 〈EN〉 - 〈ENc〉 ) g1n + g2n2

{

) Zc2Nc〈B〉 1 -

}

1 n 1 n2 〈A〉 + 2 Nc 2 Nc 〈B〉

(10a) (10b)

The shapes of the plots in Figures 3, 5, and 7 may now

Overcharging of Spherocylindrical Macroions

Langmuir, Vol. 20, No. 26, 2004 11809

Table 1. Comparison of the Parameter C (eqs 13 and 14) Obtained from the EM Simulations and from Modified Scatchard Model Fit to the Simulations for Counterions with Differing Valenciesa simulation monovalent divalent trivalent tetravalent

Cfit

Cb

nmin

〈A〉/〈B〉c

-0.123 -0.177 -0.222 -0.257

-0.130 -0.180 -0.230 -0.253

20 13 12 10

1.065 1.090 1.115 1.126

a The entries monvalent-, divalent, trivalent, and tetravalent correspond to the similarly labeled curves in Figure 3. b Calculated from eq 14. c Calculated from eq 12.

Table 2. Comparison of the Parameter C (eqs 13 and 14) Obtained from the EM Simulations and from Modified Scatchard Model Fit to the Simulations for Mixed Valency Counterionsa simulation a b c d

Cfit

Cb

nmin

〈A〉/〈B〉c

-0.181 -0.118 -0.219 -0.143

-0.228 -0.150 -0.182 -0.111

21 14 21 13

1.114 1.075 1.091 1.056

a The entries a, b, c, and d correspond to the similarly labeled curves in Figure 5. b Calculated from eq 14. c Calculated from eq 12.

Table 3. Comparison of the Parameter C (eqs 13 and 14) Obtained from the EM Simulations and from Modified Scatchard Model Fit to the Simulations for Counterions with Different Radiia simulation a b c

rc

Cfit

Cb

nmin

〈A〉/〈B〉c

1.0 × 10-10 m 1.8 × 10-10 m 2.5 × 10-10 m

-0.181 -0.177 -0.173

-0.193 -0.180 -0.180

15 14 14

1.096 1.090 1.090

a The entries a, b, and c correspond to the similarly labeled curves in Figure 7. b Calculated from eq 14. c Calculated from eq 12.

be interpreted. From eq 10 it is easily seen that for ∆En < 0, we require

{

}

〈A〉 1 (1 - n) > 12Nc 〈B〉

(11)

For n, Nc positive integers (n/2Nc)(1 - n) < 0, so that 〈A〉/〈B〉 > 1. (Note also that g1 < 0.) Indeed in I (Table 7) this relation was satisfied and found to be remarkably consistent for macroions of different geometries, and as we shall see shortly, the relation also holds true in the present case. In this model eq 10a implies that ∆En will be approximately linear for low n and quadratic at higher n passing through a minimum since g1 < 0, g2 > 0. It was seen in I that at the minimum of a ∆En versus n curve, one obtains

nmin )

(

)

〈A〉 1 + Nc -1 2 〈B〉

(12)

Thus the value of 〈A〉/〈B〉 can be estimated from nmin. Tables 1, 2, and 3 give such estimated 〈A〉/〈B〉 from nmin for the various curves in Figures 3, 5, and 7, respectively. We see again that these values are generally consistent and have very similar magnitude as that in Table 7 of I.

If 〈A〉 and 〈B〉 were truly independent of N, these could be expressed in terms of the first and second overcharging energies ∆E1 and ∆E2 in a manner similar to that used by Messina et al.,26 who eliminated the geometric factor R in the Wigner crystal (WC) model in favor of simulated ∆E1. We rewrite eq 10b in the form

∆En )

∆E1 {N C + n(n - 1)} NcC c

(13)

where

(

C)2 1-

)

〈A〉 〈B〉

(14)

Using now ∆E1 from the simulations and treating C as a fitting parameter, now identified as Cfit, we obtained the solid lines in Figures 3, 5, and 7. C was adjusted by trial and error to yield an optimum fit of eq 13 to the simulated value of ∆En. It is worth noting here the reasonable consistency between the Cfit and the corresponding value from the modified Scatchard approach, C ()2(1 - 〈A〉/〈B〉|min)), is shown in Tables 1, 2, and 3. It is of interest to compare our results for the cylinder with line charge with that of the spherical model with a point charge with regard to the counterion valency dependency of ∆En. Messina et al.26 used a Wigner crystal (WC) expression for ∆En (cf. eq 1 of I)

∆En ) -

{ (

RλeZc2 A1/2

Nc3/2 1 +

n Nc

)

3/2

}

- Nc3/2 +

n2Zc2λe 2R (15)

where R is the macroion radius, R is the aforementioned WC geometric factor, and A is the surface area of the macroion. Expanding (1 + n/Nc)3/2 binomially up to linear terms in (n/Nc) and using the neutrality condition Nc ) Zm/Zc yields

∆En ) aZc3/2 + bZc2

(16)

where a and b are constants. In particular for ∆En < 0 implies

(

)

3RZm1/2 a ) >0 b (4π)1/2 An analogous analysis for the modified Scatchard model leads to the approximate relation

∆En ) cZc + dZc2

(17)

where again c and d are constants and c/d ()(2Zm(〈A〉/〈B〉 - 1)/(n - 1))) > 0 for ∆En < 0. Equation 17 is qualitative with the curves in Figures 3 and 5. Finally, eq 10b is also useful in interpreting the behavior of ∆En with respect to counterion size shown in Figure 7. Since the counterions lie on the surface of the macroion and they do not overlap at these densities, the counterion diameter enters only through the distance of closest approach between a counterion and the (macroion) line charge. For a smaller counterion 〈A〉 will therefore be greater than the corresponding 〈A〉 for a larger counterion. The available macroion surface will also be less and as a result 〈B〉 will also increase. Consequently the quantity

11810

Langmuir, Vol. 20, No. 26, 2004

〈A〉/〈B〉 will be approximately constant, which is indeed what is observed in Table 3 for the three counterion diameters. In the right-hand side of eq 10b, we now have the quantity Zc2Nc(1 - n/2Nc + n2/2Nc - 〈A〉/〈B〉) remaining roughly constant so that by the same equation ∆En ∝ 〈B〉. Since ∆En is negative, the curve corresponding to rc ) 1 × 10-10 m lies below that for rc ) 1.8 × 10-10 m and so on and so forth. 4. Concluding Remarks The main achievement of this paper has been a further exploration of overcharging of nonspherical macroions by counterions of different valencies and mixtures of such counterions. The influence of counterion size on the ground state energetics of the macroion was also examined. Energy minimization methods developed in an earlier paper were used to determine both the energy and the distribution of counterions of finite size confined to the spherocylindrical surface in which the macroion charge is represented as an interior line charge along the cylinder axis. Relative to the neutral state, the energy of the

Mukherjee et al.

overcharged system is negative, indicating a stable system. As a function of excess surface bound counterions, the energy at first decreases in value, attains a minimum, and then increases for the systems examined in this study. A mean field type, analytical modified Scatchard model was found to predict qualitatively or better the energy gain characteristics of overcharged spherocylinders. A novel feature is that, for the cylindrical surface, the neutral state counterion distribution appears to be that of parallel helical structures wrapped about the cylindrical surface. Although such overcharged structures were obtained for isolated cylinders at T ) 0 (ground state), it is suggested that for clusters of charged cylinders some form of overcharging may still occur for room-temperature conditions. Acknowledgment. Support of an NSF Grant 0137273 is acknowledged. L.B.B. also acknowledges an institutional Grant from FIPI, University of Puerto Rico. LA049468W