Field‐theoretic simulations of polyelectrolyte complexation

Abstract

This viewpoint article is intended as a brief introduction to the emerging subject of field-theoretic simulations (FTS) of charged polymers. While the direct numerical simulation of field theory models has begun to impact several traditional areas of polymer science, including blends and block copolymers, polyelectrolytes have hitherto not been the subject of field-theoretic simulations. Here we report on a preliminary FTS study of polyelectrolyte complexation that demonstrates the potential of this novel numerical approach. Polyelectrolytes are ubiquitous in nature and in applications ranging from personal care products to paints, coatings, and processed foods. Indeed, practically all biopolymers are polyelectrolytes. In the application context, the introduction of dissociable groups is one of the most powerful ways to confer water solubility on a polymeric material. Scientifically, the polymer bound charges, which are compensated by a sea of oppositely charged counterions, produce a coupling between chain conformations and electrostatics that leads to an incredible richness of polyelectrolyte phenomena. However, this richness comes with a price: charged polymers are among the most difficult polymer systems to study theoretically or to simulate on the computer.1-4 The explanation lies in the long range nature of Coulomb interactions—charged segments feel each other at much larger distances than segments in neutral polymer systems. The situation becomes particularly challenging for the dense polyelectrolyte complexes that are the subject of the present work. Analytical theories based on assumptions of low concentration or weak interactions break down, and equilibration times in numerical simulations become prohibitively long. To conduct such simulations, polyelectrolytes are often modeled as coarse-grained chains of charged beads and the counterions are taken to be point particles, usually embedded in an implicit solvent with uniform dielectric properties. The long-ranged character of Coulomb interactions is problematic for such “particle-based” modeling approaches, however, because Ewald sums and other expensive computational techniques are required to evaluate contributions to electrostatic energies and forces that extend beyond the computational cell.5-7 The confluence of this long-range effect with the intrinsically slow kinetics of particle-based simulations at high density and molecular weight produces significant challenges for these methods. A different strategy to tackle such problems was proposed long ago by Edwards.8 The idea is to replace the coordinates and momenta of particles (polymer segments) with collective variables, or fields. One can introduce, for example, fields that describe the density of polymer segments and the density of charge. Moreover, by augmenting these fields with certain conjugate fields, such as a chemical potential (conjugate to segment density) and an electrostatic potential (conjugate to charge density), it is possible to exactly transform any classical (equilibrium) particle-based model into a statistical field theory.9 This field-based description is particularly useful for dense polymer systems, such as concentrated solutions and melts, where there is strong overlap among polymers. In such cases, the mean-field approximation, also known as self-consistent field theory (SCFT), can be applied with confidence. Within this approximation, the saddle point of the Hamiltonian dominates the partition function, and the field fluctuations around the saddle point are ignored.10 The SCFT method has yielded some remarkable analytical and numerical results for a wide variety of dense equilibrium polymer systems, perhaps most notably in the field of block copolymers.11, 12 Besides its restriction to systems at equilibrium, a major limitation of SCFT is that the assumption of negligible field fluctuations breaks down rapidly as polymers are diluted in a solvent. This is especially problematic for the study of polyelectrolytes, since they are highly solvated in most applications. Moreover, polyelectrolytes are characterized by very strong charge correlations in addition to the density correlations of neutral polymer solutions. Evidently the well-tuned machinery of SCFT is not the appropriate tool for investigating solutions of charged macromolecules. One way to address the limitations of SCFT, while still preserving the field-based collective variable approach, is to return to the exact statistical field theory and to account for the field fluctuation effects. Analytically, this procedure leads to a so-called “loop” expansion, in which systematic corrections to the free energy or other thermodynamic quantities (one-loop, two-loop, etc.) are developed in terms of integrals over field fluctuations about the saddle point of the theory. Such loop expansions can be further augmented by renormalization techniques in cases of strong fluctuations.1, 13 While powerful in the context of homogeneous phases of charged and uncharged polymers, these analytical methods can be very difficult to implement in more general inhomogeneous situations where interfaces or mesophases are present. It has recently been demonstrated that statistical field theory models of polymers can also serve as the basis for computer simulations—so-called field-theoretic simulations (FTS).9, 14 While similar conceptually to numerical SCFT, field-theoretic simulations aim to numerically sample the statistically important field configurations of the full theory, rather than just the saddle point configuration. This statistical sampling is problematic because the Hamiltonians of the relevant field theories are complex, rather than strictly real; a manifestation of the famous “sign problem”. We have found that a powerful way to circumvent this difficulty is to adopt the complex Langevin stochastic procedure,15, 16 which adaptively samples field configurations along nearly constant phase trajectories. Although FTS is computationally more expensive than SCFT, recent numerical advances have made high resolution FTS feasible for a wide variety of polymer systems.17 Polyelectrolytes have been the subject of extensive theoretical and computational research for decades.1-4, 18-33 Statistical field theory models have played a significant role in these theoretical investigations, and both mean-field and non-mean-field approaches have been employed to gain insights into the structure and thermodynamics of a wide variety of polyelectrolyte systems.21-28 Prior to this work, however, there has been no general numerical tool for simulating a field theory model of polyelectrolytes without the use of simplifying approximations. This article describes the first application of FTS to polyelectrolytes, and specifically to a phenomenon that is known to occur in aqueous mixtures of two oppositely charged polymers. Under appropriate conditions of charge density, solvent quality, and molecular weight, a phase transition can occur in such a system, with two phases being formed: one rich in both polymers and the other consisting of nearly pure solvent.18, 19, 34-40 This process is usually referred to as polyelectrolyte complexation in the physics community or as complex coacervation in the physical chemistry, colloid science, and biological communities. The resulting polymer-rich coacervate phase has two important properties: it is dense yet liquid, and it is charge neutral. Coacervates and related polyelectrolyte complexes have a variety of important applications. For example, complexes can be employed as carrier systems for charged macromolecules including protein drugs, enzymes, and DNA. In such systems, one of the polyelectrolytes serves as a chaperon with properties that can be tailored to assist targeted delivery, while the other represents the macromolecular payload.18 Another industrial use of complex coacervation is in cements, glues, or adhesives, where the coacervate phase can serve as a precursor to the formation of a solidified structural material. Examples of such coacervate precursors are also encountered in nature; for example, sand-castle worms produce a strong protein-based cement that sets in sea water and is used to construct meter-scale reefs out of sand and shells.41 Other existing and potential applications of polyelectrolyte complexes include water purification, DNA sensors,42 and encapsulation of food and pharmaceuticals. Polyelectrolyte complexation is driven by several competing factors.18, 19, 34-39 In addition to the direct electrostatic attraction between oppositely charged polyions, there are other factors playing an important role: electrostatic screening by small ions (salt or counterions), excluded volume interactions of polymer backbones mediated by the solvent, small ion translational entropy, and polymer conformational entropy. In many cases these factors are competitive rather than reinforcing, so it can be difficult to anticipate the conditions under which coacervate phases will form. For example, excluded volume interactions oppose complexation, while the translational entropy of counterions tends to favor it. To investigate the physics of complex coacervation we have chosen a simple yet fundamental model system [Fig. 1(a)]. Specifically, we consider a symmetric polycation-polyanion mixture in an implicit solvent without salt. Such a system could be realized by mixing a polyacid with a polybase in water. The symmetric assumption implies that the molecular weights and charge densities of the two types of polymers are the same, and they are combined in equal amounts. The polymer backbones are modeled as flexible (Gaussian) chains of length (number of statistical segments) N, differing only by the sign of the charge per statistical segment (charge density), ±σ. We employ a canonical ensemble in which n polyanions and n polycations are mixed to form a solution of volume V. The chains are assumed to interact by means of electrostatic and excluded volume interactions, characterized by the Bjerrum length lB = e2/ϵkBT and the excluded volume parameter u0, respectively. The solvent dielectric constant is denoted by ϵ, and e is the fundamental unit of charge. Note that the addition of salt, consideration of asymmetric polyions, or explicit inclusion of the solvent constitute only minor modifications of the formalism, although these variations will not be pursued here. Two model polyelectrolyte mixtures. Both are symmetric, with half of all chains being positively charged (shown in black) and the other half negatively charged (in gray). (a) The total charge ±σ N is evenly distributed over the entire length of each chain. (b) The same total charge is concentrated in blocks of N/4 segments near one end of each chain. Uncharged portions are shown in white. It is important to note that this field theory is an exact reformulation of the original particle-based model. While the fundamental model is not new, previous field-theoretic approaches to polyelectrolyte complexation22, 23, 27 have simplified the field theory by use of weak inhomogeneity expansions. This simplification is not exact, nor is it desirable or necessary. We further note that the long-ranged Coulomb interaction in the particle description is replaced by a short-ranged “square gradient” interaction |∇ϕ|2 in the field theory. From a computational standpoint, this is very attractive. The field-theoretic representation is also convenient for analytical studies of field fluctuation effects, such as loop expansions. A convenient place to begin our discussion of the field theory model is with the mean-field approximation, or equivalently, SCFT. The mean-field solution neglects all field fluctuations and replaces the sum over field configurations in eq 1 by a single most probable configuration–a saddle point located off the real axis in the complex plane. For the present model (subject to periodic boundary conditions) the saddle point configuration, obtained from the simultaneous equations δH / δw = 0 and δH / δϕ = 0, is homogeneous (a constant pure imaginary number) in both fields. Furthermore, because of global charge neutrality, the saddle point solution has a distinctive feature: the mean electrostatic potential is constant so that every positive charge is compensated by an equal negative charge. It follows that the Coulomb interactions are irrelevant in the mean-field limit and the polyelectrolyte mixture should behave exactly as an analogous mixture of neutral polymers. This leaves open no possibility for complexation, and indeed we shall see that it is necessary to include fluctuations and account for charge correlations to obtain a coacervate phase. Thus, the mean-field approximation breaks down completely in the case of a simple mixture of cationic and anionic polyelectrolytes. Analytically, the next level of sophistication is to account for quadratic field fluctuations about the saddle point solution in the evaluation of eq 1. This is the leading term in a systematic loop expansion13 known as the Gaussian or one-loop approximation. Calculations of this type have been reported previously for models of polyelectrolyte complexation.22, 23, 27 For our symmetric polycation-polyanion mixture, fluctuations of the w and ϕ fields decouple at the one-loop level. Moreover, it turns out that the system is parameterized by only three dimensionless combinations of variables: reduced polymer concentration C = 2nR/V, reduced excluded volume parameter B = u0N2 / R, and reduced Bjerrum length E = KdlBσ2N2 / R. Here Rg = (Nb2 / 2d)1/2 is the size of an ideal non-interacting polymer (which coincides with its radius of gyration in 3D), b is the statistical segment length, and d is the number of dimensions. Thus, the one-loop approximation provides relevant combinations of parameters responsible for the thermodynamic state of the system.43 It should be noted that the electrostatic interactions are manifested only in the parameter E. Beyond correlation lengths, the one-loop approximation provides an analytical expression for the Helmholtz free energy that can be used to derive expressions for standard thermodynamic quantities. For example, the fluctuational electrostatic contribution to the osmotic pressure scales as −Ad(EC)d/4kBTR ∼ −kBT/ξPEd, where Ad > 0 is a known numerical coefficient dependent on dimensionality d. This expression is again qualitatively different from the analogous electrostatic term in Debye–Hückel theory for simple electrolytes. Because the net contribution from charge correlations is negative, the one-loop theory can predict polyelectrolyte complexation and the coexistence of dilute and coacervate phases. The spinodals and binodals for this two-phase region follow directly from the one-loop osmotic pressure expression. While some of these predictions were obtained previously by slightly different methods,22, 23, 27 their validity has been difficult to assess because the next term in the loop expansion (two-loop order) is very tedious to evaluate. With the advent of the FTS method, we now have a numerical technique that can be used to simulate the full field theory without the assumption of weak charge and density correlations that underpins the loop expansion. Computer simulations of the full field theory can provide test beds for the analytical loop expansions, just as particle simulation methods complement analytical virial expansions in the particle-based theory. The advantage (or disadvantage) of FTS over conventional particle-based simulation techniques can be assessed by comparing computational costs. The cost of particle-based simulations depends on the total number of atomistic or coarse-grained particles in the system. State-of-the-art particle-based methods require a number of operations per MD step or MC cycle on the order of nN ln nN for n polyelectrolytes, each with N beads or segments.6 On the other hand, the cost of an FTS field update depends on the spatial discretization of the system, rather than the number of particles. When the computational cell is divided into a lattice of size M, each update requires of order NM ln M operations.9 In a melt, M < n if the lattice spacing Δx in the FTS approach can be taken larger than Vp1/3, where Vp ∼ N is the volume of a polymer. This is met, for example, by choosing Δx to be a fraction of Rg ∼ N1/2. Thus, FTS has an obvious computational advantage for dense systems of long polymers in which a coarse computational grid suffices to capture mesoscopic structure at and beyond the Rg scale. On the other hand, particle-based methods are advantaged under more dilute conditions, when fewer molecules need to be described on microscopic scales, or when full chemical details are required. These considerations are rough guides, and further studies are needed to fully elucidate the conditions under which FTS is competitive with existing simulation techniques. We employ complex Langevin (CL) sampling in our FTS to avoid the sign problem associated with the complex Hamiltonian H and the non-positive definite character of the statistical weight exp(−H) in the field theory. This method was originally devised as a strategy for sampling general types of quantum field theories with complex actions15, 44 and has been more recently applied in polymer physics.16, 45, 46 The idea behind this method is to extend the real fields into the complex plane and to compute ensemble averages of observable quantities by sampling fields along a stationary stochastic trajectory in the complex function space. While this extension to the complex plane doubles the number of field degrees of freedom, the relevant statistical weight becomes a real non-negative distribution of fields P[W,Φ], where W ≡ wR + iwI and Φ ≡ ϕR + i ϕI. A stochastic CL dynamics in the complex function space is employed to generate a Markov sequence of complex fields with stationary distribution P[W,Φ]. The complex Langevin equations are ∂W/∂t = −λ (δH/δW) + η and a similar expression for the Φ field. Here W and δH/δW (or Φ and δH/δΦ) are complex, but the thermal noise η(r,t) is a real Gaussian white noise, satisfying the usual fluctuation-dissipation theorem with dissipative coefficient λ. Since the thermal noise is placed asymmetrically only on the real part of the CL dynamics, the imaginary parts of the equations ensure the sampled fields have nearly constant phases HI along the Langevin trajectory. This removes rapid oscillations and improves convergence. Recent improvements in stochastic integration algorithms for the CL equations have enabled high-resolution, three-dimensional (3D) field-theoretic simulations.17 Here we report on the application of CL-FTS to polyelectrolyte complexation phenomena, specifically to the same symmetric binary polycation-polyanion model that was used for the analytical one-loop calculations. A particular focus of our CL-FTS study is the location and size of the two-phase region where coacervation takes place. An example of such a phase diagram (in 3D) is provided in Figure 2. Spinodals and binodals for the field theory of eq 2 are surfaces in the three-parameter space of the reduced variables C, B, and E. The figure represents a cross-section of this three-dimensional space by a plane E = 14,400; hence the diagram involves only the C and B variables. The diagram features both the CL-FTS results (symbols and dotted-line fit) and the one-loop analytical approximation to the binodal (solid line) and the spinodal (dashed line). The one-phase region (disordered homogeneous phase) is above and to the right of the lines, and the two-phase region is below and to the left. The tie lines in the two-phase region are horizontal (constant B) and connect nearly pure solvent (C ≈ 0) with a coacervate phase at the binodal concentration. The two symbols for the CL simulation data at each concentration C represent the hysteresis upon varying the solvent quality B, with the upper symbol corresponding to superheating (increasing B) and the lower one corresponding to supercooling (decreasing B). Phase diagram for the symmetric polyelectrolyte mixture in 3D plotted in reduced polymer concentration C and reduced excluded volume B at fixed reduced Bjerrum length E. Solid and dashed lines are the analytical one-loop binodal and spinodal, respectively. Symbols are the result of complex Langevin simulations; the dotted line is a power-law fit. Numerical simulations were conducted in a cubic cell of size 4Rg × 4Rg × 4Rg with periodic boundary conditions. Remarkably, the analytical and numerical results nearly coincide in the high concentration region of the figure, despite the limitations of each method. The numerical results are subject to finite cell size and chain discretization limitations, while the analytical predictions neglect two-loop and higher order terms in fluctuations. Nonetheless, our numerical supercooling result practically follows the analytical spinodal, and the analytical binodal yields nearly the same exponent (−1.31) as is obtained from a power-law fit to the numerical superheating points (−1.40). The cause for the discrepancy between theory and simulation for the location of the binodal is unclear at present, although the overall semi-quantitative agreement indicates that both approaches have utility for this class of problems. Our future work in polyelectrolyte complexation will extend beyond the symmetric model presented here to include unequal chain lengths, counterions, salt, and explicit solvent. The latter will allow for the treatment of polyelectrolytes with hydrophobic backbones. It is important to emphasize that these are straightforward extensions that do not complicate the particle-to-field transformation of the model, nor do they make analytical loop expansions or CL simulations any more difficult. Another extension is to investigate the effect of charge distribution and polymer architecture on polyelectrolyte complexation. This brings us into the realm of block copolyelectrolytes and branched and dendritic polyelectrolyte systems—all of considerable experimental interest. As a preliminary example, we have explored a simple variation of our binary polycation–polyanion model where all the charge of each species is concentrated into a small “charged block” at the chain end; compare Figure 1(b). The remaining longer portion of each chain constitutes a “neutral block”. In our simple model, the charged blocks should drive complexation because of the correlation-induced electrostatic attraction, while the neutral blocks repel each other due to the excluded volume interactions in the (assumed) good solvent. These competing tendencies can drive aggregation to form micellar structures with charged blocks in the micelle cores and the neutral blocks in the coronas. At relatively low concentrations, the micelles are expected to form a disordered micellar phase, while at higher concentrations they can pack into periodic lattices to form mesophases of various symmetry. Such structured coacervate phases would seem to be of considerable interest for a variety of applications. In accord with expectations, complex Langevin simulations of the symmetric binary block copolymer model show that the redistribution of charges along the polyelectrolyte chains leads to the formation of aggregates qualitatively similar to micelles. Figure 3 illustrates the dramatic difference between the phases formed by uniformly charged chains and unevenly charged chains with the same total charge concentrated into a block of length N/4. The uniformly charged chains [Fig. 3(a)] form macrophases, that is, large homogeneous phases (a dilute phase and a coacervate phase) constrained only by the size of the computational cell. On the other hand, the oppositely charged block copolymers [Fig. 3(b)] form a mesophase with structure on the scale of Rg. The core of each micelle predominantly consists of charged blocks, while the neutral blocks form the repulsive corona. We are aware of at least one body of experimental work that confirms these (as yet) qualitative predictions.47, 48 Snapshots of equilibrated field-theoretic simulations. Normalized total polymer density is presented for the cases of (a) uniformly charged chains and (b) unevenly charged chains, with the same amount of total charge concentrated on 25% of the chain length. Macrophases are formed in the uniformly charged case, and a mesoscopic structure develops in the case of block polyelectrolytes. Both systems are in 2D with periodic boundary conditions. Cell size is 16Rg × 16Rg, and the same parameters C = 6.0, B = 0.3, and E = 64,000 were used in both cases. In summary, field-theoretic simulation methods based on the introduction of auxiliary fields have considerable promise for the investigation of polyelectrolyte complexation phenomena. These methods possess several distinctive and attractive features. First, FTS involves exact Hamiltonians and accounts for arbitrarily large fluctuations and strong inhomogeneities. Thus the technique is not limited in the same ways as the available analytical tools. Second, field-theoretic methods are especially convenient for treating the long-range Coulomb interaction, which is replaced by a short-ranged square-gradient operator in the auxiliary field representation. Third, FTS methods are computationally advantaged over particle-based simulations when the systems are dense, polymer chains are long, and the length scale of interest is large (mesoscopic). Finally, inclusion of additional species (either long chains or small ions) neither does lead to a substantial elaboration of the field-theoretic models, nor does it complicate the simulations. With continued advances in algorithms, we anticipate that FTS methods will enable the study of broad classes of polyelectrolyte systems and phenomena that are beyond the reach of current analytical methods and particle-based computer simulations. The importance of analytical theory and particle-based simulations will not be diminished, however, by the emergence of FTS techniques. We see these approaches as complementary, with each contributing valuable insights into the rich structure and thermodynamics of charged polymeric fluids. The authors are grateful to Fyl Pincus and Kirill Katsov for many valuable discussions and advice. Acknowledgement is made to the Donors of the American Chemical Society Petroleum Research Fund, the Institute for Collaborative Biotechnology, Rhodia Corporation, and the Mitsubishi Chemical Corporation for the support of this research.