Field‐theoretic simulations of polyelectrolyte complexation

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 In a M n the lattice in the FTS can be taken larger than where ∼ N is the volume of a This is for example, by to be a of Rg ∼ Thus, FTS has an computational advantage for dense systems of long polymers in which a computational to structure at and beyond the Rg On the other hand, particle-based methods are under more dilute to be on or full chemical are are and further studies are to the conditions under which FTS is competitive with existing simulation We employ complex Langevin sampling in our FTS to the sign with the complex Hamiltonian and the character of the statistical weight in the field theory. This method was as a strategy for sampling general types of field theories with complex and has been more recently applied in polymer The idea this method is to extend the real fields into the complex and to ensemble of quantities by sampling fields along a stochastic in the complex While this to the complex the number of field of the relevant statistical weight becomes a real of fields where and A stochastic in the complex is employed to a of complex fields with The complex Langevin equations are = and a similar expression for the Here and (or and are complex, the is a real Gaussian the with coefficient the is only on the real of the the imaginary of the equations the fields have nearly constant phases along the Langevin This and in stochastic for the equations have field-theoretic Here we report on the application of to polyelectrolyte complexation specifically to the same symmetric polycation-polyanion model that was used for the analytical one-loop A of our study is the and size of the two-phase region where coacervation of such a phase (in is in and binodals for the field theory of eq are in the of the reduced and E. The represents a of this by a E = the only the C and B The both the results and and the one-loop analytical approximation to the and the The region homogeneous is and to the of the and the two-phase region is and to the The in the two-phase region are and nearly pure solvent with a coacervate phase at the The two for the simulation at each concentration C the the solvent with the to and the one to for the symmetric polyelectrolyte mixture in in reduced polymer concentration C and reduced excluded volume B at reduced Bjerrum length E. and are the analytical one-loop and respectively. are the of complex Langevin the is a simulations were in a cell of size with periodic boundary the analytical and numerical results nearly in the high concentration region of the the limitations of each The numerical results are subject to cell size and chain discretization while the analytical predictions and order terms in our numerical practically follows the analytical and the analytical nearly the same as is obtained from a to the numerical The for the between theory and simulation for the of the is at although the that both approaches have for this of in polyelectrolyte complexation will extend beyond the symmetric model to include chain lengths, counterions, salt, and explicit solvent. The will for the of polyelectrolytes with It is important to that these are that not the of the nor they analytical loop expansions or simulations any more Another is to investigate the effect of charge and polymer on polyelectrolyte This into the of block and and polyelectrolyte of a preliminary example, we have a simple of our model where all the charge of each is concentrated into a small at the chain The of each chain a In our simple the charged blocks should complexation because of the electrostatic while the neutral blocks each other to the excluded volume interactions in the solvent. competing can to form with charged blocks in the and the neutral blocks in the low the are to form a while at they can into periodic to form mesophases of Such coacervate phases to be of for a variety of applications. In with complex Langevin simulations of the symmetric block model that the of along the polyelectrolyte chains leads to the formation of qualitatively similar to the between the phases by charged chains and charged chains with the same total charge concentrated into a block of length The charged chains [Fig. form that homogeneous phases (a dilute phase and a coacervate only by the size of the computational On the other hand, the oppositely charged block [Fig. form a with structure on the of The of each of charged while the neutral blocks form the We are of at one of that these of field-theoretic simulations. total polymer density is for the cases of (a) charged chains and (b) charged with the same of total charge concentrated on of the chain are in the charged and a structure in the case of block polyelectrolytes. Both systems are in with periodic boundary size is and the same parameters C = B = and E = were used in both In field-theoretic simulation methods based on the introduction of fields have for the of polyelectrolyte complexation phenomena. methods several distinctive and FTS exact Hamiltonians and for fluctuations and strong the technique is not in the same ways as the analytical field-theoretic methods are especially convenient for the long-range Coulomb which is replaced by a short-ranged in the field FTS methods are computationally over particle-based simulations the systems are polymer chains are and the length of is inclusion of long chains or small to a of the field-theoretic nor it the simulations. With advances in we anticipate that FTS methods will the study of of polyelectrolyte systems and that are beyond the of analytical methods and particle-based computer simulations. The of analytical theory and particle-based simulations will not be however, by the of FTS We see these approaches as with each insights into the rich structure and thermodynamics of charged polymeric The are to and for many and is made to the of the the for and the for the of this