Potential of mean force and transient states in polyelectrolyte pair complexation

Polyelectrolyte complexes (PECs) consisting of oppositely charged polymer chains have been studied for decades owing to their fundamental importance in biophysics and technological applications.^1,2 Common examples include delivery vehicles for gene therapy and oral vaccination,^3–6 neurofilament association,⁷ membrane support for filtration processes,^8,9 functional coatings,¹⁰ and highly ordered macromolecular structures for water treatment and mining.^11–13 A major focus of previous investigations has been on the physicochemical nature of the complex-forming polyelectrolytes (PEs) and environmental conditions that control the structural and topological properties. The asymmetry between polyions of opposite charges leads to PECs with a broad range of conformations and size distributions, suggesting that the PEC size and topology^1,2 as well as assembly kinetics^14,15 and macroscopic phase behavior^16–21 can be controlled by alternating the polymer chemistry as well as solution conditions.

In particular, polycations for gene delivery are designed to associate with nucleic acids so that the gene transfer vectors are able to overcome the intracellular barriers such as the plasma membrane, the endosome, and the nuclear membrane.³ Successful gene delivery depends upon the ability of the vector to adopt different metastable and transient structures or possess certain properties at different intracellular environments. It has been shown, for instance, that the complex of small interfering RNA (siRNA) with poly(amidoamine) (PAMAM) undergoes aggregation and precipitation at different stages of delivery, and the dynamic behavior is reported to possibly facilitate the transaction efficiency.⁶ Different binding characteristics were reported when a DNA helix interacts with synthetic polymers of different charge density (e.g., polyethylenimine (PEI) and poly-l-lysine (PLL)).²² Because the effective polymer charge density is closely linked to counterion-condensation effects,^23–27 the appearance of different complexation and transient structures depends on the ionic strength and counterion valence.

In order to understand complex formation between PEs on a detailed molecular level (in both pair or many-body PE systems), monomer-resolved computer simulations are a powerful tool. Examples include atomistically resolved systems^28–33 as well as coarse-grained models.^34–42 In addition to detailed structural and topological insights into the final complexed state, those studies indicate that the primary driving forces of PE complexation for weakly charged PEs arise as expected from electrostatic attraction. By contrast, association between highly charged PEs is dominated by the release of condensed counterions if the polymer charge density is above the threshold value for the phenomenon of counterion-condensation to occur.^23–27 As systematically shown by Ou and Muthukumar using coarse-grained simulations and a mean-field lattice model, counterion-release effects have remarkable consequences on the thermodynamics of complexation, which is of mostly entropic origin for highly charged PEs, while mostly enthalpic for the weakly charged ones.⁴⁰ The entropic effect on the binding free energy appears to be consistent with the classical Record-Lohman picture for PE pair complexation^43,44 and has implications for statistical theories describing macroscopic phase behavior.^16–20 

Despite the extensive simulation work on PEC formation on a pair level, to the best of our knowledge, little attention has been given to the variation of the system’s free energy and structures along the PE pair association pathway. The previous simulation studies indicate very fast and cooperative association processes once the isolated PE coils come into first contact.^34,37,40 Here, apparently, transitional but very distinctly stretched states of the PEs appear to govern the intermediate states of kinetic association towards the final complex. Very recently, Peng and Muthukumar⁴² calculated the potential of mean force (PMF) along the center-of-mass (COM) distance reaction coordinate and reported a constant force of attraction (linear PMF) for a large distance range comparable to the PE contour length. A systematic exploration of the pair PMFs for different PE charge densities, however, is still absent, as well as their quantitative description by a tractable theory. Also the nature of the transient states right at the onset of attraction and during complexation is not fully characterized, despite the need for a better understanding of PEC metastable states and preceding association kinetics, e.g., in gene delivery and therapy.^3,6,29

In the present work, we are concerned with the process of complex formation between two oppositely charged PEs right at the onset of the association and the intermediate range before collapsing into the final state. For this, we investigate a pair of two oppositely charged PEs of the same size and charge density using implicit-water, explicit-salt Langevin simulations of a generic (bead-spring) PE model for various charge densities. Our main focus is placed on the PE configurations and the association free energy as well as the number of released ions along the center-of-mass distance reaction coordinate. Above the counterion-release threshold charge density, we find very distinct sliding-rod association pathways accompanied by an essentially constant mean force. We show that the latter can be described by a combination of a simple Debye-Hückel (DH)-like theory for associating rods and the counterion-release entropy, which dominates the PE chains with large charge densities. Importantly, we observe that the mean force is discontinuous at the onset of complexation. We explain this observation by a bimodal structural behavior that arises from the coexistence of interconnected-rod and isolated-coil states. The two intermediate states of the PE complex are balanced by subtle counterion-release mechanisms and separated by a free-energy barrier due to unfavorable stretching entropy of the PE chains. We discuss possible consequences of our findings briefly in Sec. IV of this work.

Consider an aqueous solution containing two polyelectrolyte (PE) chains of the same size but opposite charge in the presence of salt ions at a finite concentration. We treat the solvent implicitly via a uniform dielectric background with a dielectric constant of water at room temperature, ϵ_r = 78.44. The PE chains are represented by a coarse-grained beads-on-the-string model.⁴⁵ Approximately, each bead represents a monomer for a realistic PE chain. In addition to electrostatic interactions, the polymer beads interact with each other in terms of the Lennard-Jones (LJ) potential,

with a diameter σ_LJ = 0.3 nm and an energy ϵ_LJ = 0.1 k_BT. Monovalent salt ions are modeled explicitly as charged beads with the same LJ potential as that for the PE beads. The valency of small ions is z_± = ± 1.

The PE chain connectivity is imposed by a harmonic potential,

where l represents the distance between consecutive beads and l₀ = 0.4 nm is the equilibrium bond length. The spring constant is k_b = 4100 kJ mol⁻¹ nm⁻². To account for the polymer backbone flexibility, we restrain the bond angle also by a harmonic potential,

where ω is the angle between a triplet of adjacent beads and ω₀ = 120° its prescribed equilibrium value. The potential constant is k_a = 418 kJ mol⁻¹ rad⁻².⁴⁶ 

In this work, we consider relatively short PE chains with the total number of beads N_b = 25, close to the degree of polymerization for PEs used in some related experimental studies.^6,47–50 The contour length of the PE chains is L_c ≃ (N_b − 1) l₀sin(ω₀/2) ≃ 8.3 nm, where the sine-function takes care of the bond angular constraints in a moderately stretching regime. Each bead carries a bare Coulomb charge $|qb|=e|zb±|$⁠, with e being the unit charge and $|zb+|=|zb−|≡zb$ the valency. The total charge of each PE chain is |Q_b| = N_b|q_b|. Since we consider complex formation between two anti-symmetric PE chains, the final complex is electroneutral. The two polymers have the same properties except with opposite charges, that is, Q_b > 0 for one PE chain, and Q_b < 0 for the other.

Our computer simulations are based on the Langevin dynamics (LD). The equation of motion for each bead is given by⁵¹ 

where m_i and ξ_i are the mass and the friction coefficient of the ith bead (or small ion), respectively. All beads have the same unit mass, which is set to minimize the inertia effects and enhance sampling. The choice of particle mass does not affect the equilibrium properties of polymers studied in this work. The total potential energy of the entire system is the sum of all the contributions U = U_Coul + U_LJ + U_bond + U_angle. The random force R_i(t) in Eq. (4) has zero mean and its autocorrelation function satisfies the fluctuation-dissipation theorem,⁵² 

The Langevin friction is chosen as ξ_i = 1.0 ps⁻¹ such that it dissipates energy at constant temperature T = 300 K on the time scale much faster than those governing the dynamics of the polymer system. To integrate the equations of motion, we employ the leap-frog algorithm with a time step of 2 fs. The simulations are carried out by the GROMACS 4.5.4 software package.⁵³ 

In the production runs, we use a cubic box with a side length of L = 30 nm with periodic boundary conditions in all three directions. The center of mass translation of the system is removed in every 10th integration time step. The electrostatic interactions are treated with the Particle-Mesh-Ewald method (PME),⁵⁴ where the long-range potential is evaluated in the reciprocal space using the Fast Fourier Transform (FFT) with a grid spacing 0.12 nm and the cubic interpolation of the fourth order. We use a cutoff radius r_cut = 3.0 nm for the short-range electrostatics and the LJ interactions. The choice of the cutoff is verified by reference simulations with increased cutoff value r_cut = 5.0 nm. For PE chains with the highest charge density simulated in this work (Manning parameter ξ = 2.31, defined below), where the electrostatic interactions are most significant, the potential of mean force (PMF) curves deviated less than 6% from the results treated with r_cut = 3.0 nm. In all simulations, the PEs are immersed in a salt solution with N_i = 325 pairs of cations and anions, resulting into a salt bulk concentration of c₀ = N_i/L³ ≃ 20 mM.

According to the Onsager–Manning–Oosawa theory,^23–27 the net charge of a highly charged PE is renormalized by partial condensation of the surrounding counterions. The condensed and free ions are considered separately as two different liquids. While the counterions condensed at the surface of the PE backbone form a strongly correlated liquid, the free ions in the diffusive double layer^55,56 can be approximately described by a Debye–Hückel (DH)-like mean-field theory. The analytic approach will be utilized later in this work to develop a simple expression for the PMF between two anti-symmetric PEs.

The propensity of counterion condensation is characterized by the Manning parameter,^23–27 

Here, |z_±| stands for the counterion valency (i.e., |z_±| = 1 in our case), λ = |q_b|/[l₀sin(ω₀/2)] is the bare line charge density, and l_B = e²/(4πϵ₀ϵ_rk_BT) stands for the Bjerrum length, which is l_B = 0.71 nm for water at room temperature as used in our simulations. Note that due to thermal fluctuations of the flexible chain, the bare charge density is defined as above strictly only in the low temperature limit.

The Onsager–Manning–Oosawa theory^23–27 predicts that counterion condensation occurs if ξ > 1, i.e., if the Manning parameter exceeds unity. In this work, we consider a systematic scan of the Manning parameter in the range from ξ = 0.29 to ξ = 2.31, corresponding to the bare charge density from λ = 0.412 to 3.3 e nm⁻¹ tuned by the monomeric (partial) charge between q_b = 0.142 86e and q_b = 1.142 86e, respectively. When ξ > 1, the two-state model predicts a fraction θ = 1 − 1/ξ of counterions in the condensed state. The conventional model for ion distribution in the diffusive double layer is based on a cylindrical cell model with the PE chain placed at the center.^57,58 If the polymer backbone is treated as a cylinder, the effective, renormalized charge density of such a cylinder with ξ > 1 is then

resulting in a Manning parameter of unity. The cell model provides simple expressions for the condensation threshold radius r₀, the Manning parameter ξ, and the fraction θ of condensed counterions.⁵⁹ It predicts that the fraction θ of condensed ions is independent of the cell size.^25–27,60 When ξ < 1, counterion condensation does not occur and the PE bare charge and the effective charge are the same.

In order to estimate the number of condensed counterions, N, in our simulations, we count an ion as condensed if it is located within a radial distance r₀ from any bead of a PE chain with the opposite electrostatic charge. To avoid double counting, each ion is counted only once according to the nearest distance to the polymer beads. Understandably, the radius r₀ lacks a rigorous definition (see, e.g., the discussion in Refs. 25 and 61), in particular for flexible PEs in an electrolyte solution. In this work, we propose an unambiguous procedure to account for counterion condensation: the radius r₀ is fixed such that the Manning prediction of N = λ/e(1 − 1/ξ) L_c is obeyed for each PE chain. This procedure yields r₀ = 0.5, 0.56, 0.63, 0.68, and 0.74 nm for ξ = 1.15, 1.44, 1.73, 2.02, and 2.31, respectively. These values are reasonable because it is known that r₀ increases monotonically with increasing ξ and that it should be identical to the effective PE (modeled as rod) radius for ξ = 1,²⁵ which is about 0.3–0.4 nm for our PE model. Similar values of r₀ can be estimated from an inspection of the radial distribution of counterions (not shown) around the PE beads (see, e.g., Refs. 46, 47, and 61).

To obtain the PMF between two oppositely charged PEs, we use steered Langevin Dynamics (SLD) simulations with the “pull-code” as implemented in GROMACS.⁵³ Here, the distance r between the centers-of-mass (COM) of the two PE chains serves as the reaction coordinate, which is constrained by an external harmonic potential. The fictitious potential is time-dependent such that it exerts a pulling force to steer the PE chains moving in a prescribed direction with a prescribed velocity v_p. An exemplifying snapshot of two PEs before association is shown in Fig. 1(a).

We have empirically tested a variety of steering velocities to make sure that the PE drift is slow enough to sample the equilibrium state, i.e., the simulation results are independent of the choice of v_p. For all results reported in this work, the steering velocity v_p = 0.2 nm/ns is used along with the harmonic force constant K = 2500 kJ mol⁻¹ nm⁻². During the simulation, we steer the constrained PE chains approaching each other from a well separated state (r ∼ 12 nm) to the final state (r ∼ 0.5 nm). The friction force f_f = − mξ_iv_p is subtracted from the constraining force and the result is averaged within a specific interval of the discrete spacing Δr to obtain the mean force of the interaction potential. After that the PMF profile is acquired with a backward integration. Since one of the PE chains is radially constrained in a three-dimensional space, we need to subtract the center-of-mass translational entropy,^25,62–64

where W^I(r) is the integrated mean force and D = 3 is the dimensionality of the external constraint.

We use an analogous procedure to calculate the free energy of stretching a single PE chain. In that case, the distance between the head and tail monomers from the same PE chain serves as the reaction coordinate.

In order to characterize and demonstrate the mutual alignment of two highly charged PEs during complex formation, we introduce the order parameter m_pp defined as

Here, $t→ic$ is the intra-PE direction vector connecting one terminal bead and the central bead of the same polymer i, with i = 1, 2, whereas $t→cc$ is the inter-PE direction vector linking the central beads of both PEs. The notation ‖…‖ represents the norm of a vector, and 〈…〉 denotes the ensemble average. For stretched, rod-like polymer configurations, we expect parallel alignment of both vectors $t→ic$ with $t→cc$⁠, i.e., the PE chains are aligned with themselves and with the connection axis. In that case, the order parameter approaches unity for both PE chains, $mppi≲1$⁠. If two parallel rods are perpendicularly approaching, the order parameter stays as zero for both PE chains, $mppi≃0$⁠. If rods are perpendicularly approaching, the order parameter is either $mpp1≃0,mpp2≃1$ or $mpp1≃1,mpp2≃0$⁠. If the PEs are uncorrelated, we would expect an average value of the order parameters much closer to $mppi≃1/2$⁠. We will plot and discuss the distance-resolved value for $mpp(r)=[mpp1(r)+mpp2(r)]/2$⁠, which is averaged over the two PE chains.

Moreover, we monitor the distance-resolved end-to-end distance R_ee(r) of the PEs as they approach each other. A stretching of the PE chains can then be identified simply from the increase of R_ee(r) at a certain COM distance.

To attain a better understanding of the thermodynamic driving forces and the mechanisms of PE–PE association, we now introduce a simple electrostatic model for the PMF between highly charged PE chains (ξ > 1). We assume that the PEs adopt an isolated, coil-like state when they are far apart, as sketched in Fig. 1(a). In this case, the PEs interact roughly with a DH-like potential $VDH(r)∝Qr2lBexp(−κr)/r$⁠, where Q_r ≃ 11.7e is the renormalized net charge of each PE and $κ=8πlBc0$ the usual DH screening parameter. For distances larger than the contour length r ≳ L_c, this accounts only for weak attraction (≲k_BT) and will be discarded in the following discussion. We further assume that during their association at smaller distances r ≲ L_c, the two PEs are in a parallel sliding conformation, as depicted in Fig. 1(b∗), until they reach a final state with r close to zero. In the latter, the PE chains have collapsed into a globular complex, cf. Fig. 1(c). Consequently, we assume that the PMF is dominated by the parallel sliding process for a wide range of separations.

We model the two approaching PE chains in their fully stretched association configuration with the simplest analytically tractable model, that is, two parallel, infinitely thin, and oppositely charged rods of length L_r. We define the effective rod length L_r as the COM distance at the onset of the elongated PE configuration, which can be read off in Fig. 2(a) at the intersections of the diagonal dashed line, R_ee = r, and the profile of the end-to-end distance R_ee(r). The effective length L_r increases with increasing ξ, since larger electrostatic attraction promotes the jump into the elongated state at larger separations. For small values of ξ, the effective length L_r is significantly smaller than the contour length L_c = 8.3 nm of the chains, since the attraction in that case is not strong enough to extend the chains considerably. However, for higher charges, that is, for higher values of ξ, the effective length approaches L_c.

The COM distance of the two rods is again denoted by r, with the parallel component r_∥ and the perpendicular component r_⊥, as shown in Fig. 1(b∗). The latter is kept fixed and reflects the closest distance of the two PEs, which is related to the LJ diameter σ_LJ. We now assume that the major interaction contribution to the total free energy arises only from the neighboring parallel segments of length L_r − r from each rod. In the following, we will evaluate in detail the electrostatic contribution based on a DH approach and additionally account for the purely entropic contribution of the counterion release effect in the scenarios with ξ > 1.

We first evaluate the DH energy corresponding to the electrostatic interaction between two parallel, partially neighboring rods. Assuming r_⊥ ≪ L_r, we may neglect edge effects and approximate the pair potential between the neighboring segments as those from infinitely long rods. With the pairwise additive approximation, we can derive the overall interaction energy then as⁶⁰ 

Here, β = 1/(k_BT), λ is the bare line charge density of a PE rod, which can be renormalized by condensed counterions. The fraction of neutralization is θ = 0 for ξ < 1 and θ = 1 − 1/ξ if ξ > 1. Note that in the latter case, the effective charge density stays as e/l_B, irrespective of the COM distance r. Performing the above integral yields^60,65

where K₀(κr_⊥) is the Bessel function of the second kind. According to Eq. (11), the electrostatic DH part of the PMF is linear in the center-of-mass distance r, which is a logical outcome of our assumption, since only neighboring segments of the parallel aligned rods contribute. Note again that because of the charge renormalization effect for ξ ≥ 1, this contribution does not explicitly depend on the charge density.

For ξ ≥ 1, in addition to the direct electrostatic interactions between polyions, we may estimate the contribution to the free energy of PE complexation due to the release of counterions. For association between two oppositely charged PEs with the Manning parameter ξ ≫ 1, the free energy is dominated by the entropically favored release of condensed counterions. The idea was first proposed by Lohman and co-workers^43,44 and supported later by other theoretical investigations including coarse-grained (but explicit-salt) computer simulations.^40,46 The entropy gain upon releasing n bound counterions is

where c_b corresponds to the effective density of the condensed ions in the vicinity of the PE and c₀ is the ion bulk density. The corresponding free energy gain due to counterion release is then given by

For complex formation between two PEs of opposite electric charge, each PE chain carries condensed counterions of a line density λθ, where θ is the fraction of the PE charge neutralization by counterions. When the PE segments approach each other, all counterions in the overlapping region are liberated and released into the bulk solution. The number of released ions is n = 2 × λθ(L_r − r)/e. In terms of the Manning parameter, n can be expressed as

The density of the bound ions^46,47 is defined by the number of condensed counterions N and their occupied volume V, i.e., c_b = N/V. The condensed counterions reside around each rod at a radial distance between the inner radius σ_LJ and the outer radius r₀. The occupied volume thus corresponds to that of two hollow cylinders, V = 2s_c(L_r − r), where $sc=π(r02−σLJ2)$⁠. Accordingly, the density of the bound ions for ξ > 1 can be estimated as

The free energy contribution from the released counterions is then given by

The above expression is valid only for ξ > 1; otherwise, W_ion = 0 since there are no counterions condensed on the PE surface. For a salt concentration c₀ = 20 mM and, for instance, ξ between 1.15 and 2.31, Eq. (15) predicts c_b ≃ 0.7 M and 2.12 M, respectively, which corresponds to 3.6−4.7 k_BT of dissociation free energy per single released ion, respectively.

Finally, we sum up the DH contribution W_DH and the counterion release W_ion, viz.,

Equation (17) predicts a PMF between PE chains to be linearly dependent on the COM distance r. As previously concluded in related simulations of PE complexation,⁴⁰ the PMF is dominated by electrostatic enthalpy for ξ < 1 and counterion-release entropy for ξ ≫ 1.

It should be noted that Eq. (17) is valid only for a finite length of overlapping rod segments, r < L_r. It neglects smaller contributions for larger distances r > L_r. The separation distance between the rods r_⊥ is approximately equal to the LJ diameter σ_LJ. However, r_⊥ could be smaller than that since a bead of one PE chain can sit in the region between the two neighboring beads of the other PE chain. In this work, we obtain the value for r_⊥ by the best fit of Eq. (17) to the simulation data for ξ = 1.0, yielding r_⊥ = 0.28 nm. We use this value for all other ξ-values.

Finally note that the two PEs in the overlapping parts are almost free of the surrounding ions and therefore interact as (almost) bare charges via unscreened Coulomb force. However, the total free energy of the PEs, which is the work of bringing the two PEs from large separations to close distance, cannot be computed in terms of the unscreened electrostatics, as at large distances the screening always becomes significant. Furthermore, the electrostatic attraction should comprise the portion of the energy penalty to deprive the condensed counterions on the overlapping part. Thus, attraction only in terms of bead bare charge will exclude that energy and eventually lead to an underestimate of the PMF. As an approximation, we adopt the effective charge to implicitly account that energy penalty. Our intention with this simple model is to qualitatively describe the main characteristics of the rather complex association process obtained by the simulations. As seen by the nice agreement between the theory and the simulations later in Fig. 3(b), the used model is able to capture the basic features of the association.

A first insight into the configurations of the PE chains along their approach to association can be attained by analyzing the distance-resolved end-to-end distance R_ee(r) of the individual PEs. Figure 2(a) shows the average R_ee(r) of the antisymmetric PEs with different Manning parameters. At large separations, where the PE chains are in a coil-like state and independent of each other, cf. Fig. 2(a), R_ee(r) is constant. The effective size of the coil increases with increasing ξ due to self-electrostatic repulsion within the PE chain. For highly charged chains, ξ ≫ 1, significant and steep, almost discontinuous increase of the end-to-end distance can be observed at COM distances close to the contour length r^∗ ≃ 5–7 nm. The exact value of the critical distance r^∗ depends on ξ. For larger ξ, the critical distance increases towards the maximum contour length and the magnitude as well as steepness of the jump grow. This behavior strongly indicates stretching of the PE chains towards each other and “handshake”⁷ at a critical distance r^∗ to assume a maximum end-to-end length $Ree∗$ as shown in the snapshots in Fig. 1(b). An observed linear relation $Ree∗∝r∗$⁠, cf. Fig. 2(a), further consolidates that proposition. At small distances, r ≲ 2 nm, the PE chains seem to collapse again, with R_ee(r) approaching the values that correspond to the isolated coil states.

The behavior of R_ee(r) for the associated states strongly correlates with the alignment order parameter m_pp, shown in Fig. 2(b). At large distances r ≳ 6.5 nm, the PE chains are independent of each other and uncorrelated in any alignment, as indicated by Fig. 1(a). In that case, the order parameter tends to the value m_pp ≃ 1/2. At the intermediate distances r ≃ 5–7 nm, again cf. Fig. 1(b), the ends of different PEs jump together and the polymer chains become elongated and m_pp rises to high values, indicating parallel alignment of stretched PEs. The order parameter in this regime increases with increasing ξ. Particularly for ξ = 2.31, the order parameter reaches m_pp = 0.95, implying that the two PEs get almost completely stretched and perfectly aligned. At smaller distances, r ≲ 2 nm, the order parameter tends to the value m_pp ≃ 1/2, and the previously extended PEs collapse into a compact globule. In a relatively large spatial regime roughly appearing at distances between 1 − 2 nm < r < 5 − 7 nm (with the exact values depending on ξ), both PEs slide along each other in a stretched and parallel configuration, with their orientations pointing along the connection axis.

An interesting question arises as to why the PEs tend to associate via the elongated sliding configuration rather than via other possibilities, such as perpendicular or parallel sideward approach. Fixing the COM separation between the two PEs, the elongated configuration allows the terminal segments on the respective PEs coming into a close contact. In presence of the strong electrostatic attraction and counterions release, this configuration has the smallest free energy compared with all other realizations. Similar elongated sliding approaches of PEs have been reported (but remained undiscussed) in a number of previous publications in coarse-grained as well as atomistic-level simulations.^40,42,66 This sliding mechanism has implications on the interpretation and theoretical description of the PMFs along their COM distance reaction coordinate, as discussed in the following.

In Fig. 3, we present the PMF, W(r), between the two PE chains with different Manning parameters ξ. Here the salt concentration is again c₀ ≃ 20 mM. Figure 3(a) shows the PMFs for ξ ≤ 1 and (b) for ξ ≥ 1. The three generic stages of the PE–PE association discussed before are reflected also in the behavior of the PMF: When the PE chains are far apart, r ≳ 6.5 nm, they do not significantly interact with each other (on the shown scale) regardless of the Manning parameter. Once the PE chains begin to associate at r^∗ = 5–7 nm and stretch at intermediate separations 1–2 nm < r < r^∗, W(r) becomes nearly a linear function of r, with the slopes increasing with increasing ξ. Consequently, the mean force f = − dW(r)/dr, shown in the insets, can be regarded as nearly constant in that r-range, as reported already in previous simulations.⁴² However, some slope in the mean force is clearly visible, especially for the smaller ξ-values, so the assumption of a strictly constant mean force is not generally true. Strikingly, the mean force f = − dW(r)/dr exhibits a discontinuity at r^∗ for ξ ≳ 1, see the inset to Fig. 3(b), a fact that will be discussed later in more detail. At smaller separations r < 1–2 nm, the PE chains tend to intertwine into a collapsed globule and by that further increasing the association free energy. In this regime, the corresponding attraction grows even stronger (superlinear) with the distance. The value of W(r ≃ 0) at the associated state represents the free energy of PE–PE complexation. Its value increases with the Manning parameter ξ, as can be expected from increased Coulomb attraction between both PEs. The thermodynamics in terms of enthalpy and entropy of the final complex was investigated in detail by Ou and Muthukumar.⁴⁰ 

In the same figures, we plot the PMF predictions of our simple theory given by Eq. (17) by dashed lines. We note again that for each ξ, we determine the effective length L_r at the onset of attraction of stretched PEs directly from Fig. 2(a) and use it as an input to the theory, which is applicable for r < L_r. Another parameter in our analytical model, the rod distance r_⊥ = 0.28 nm, is fixed by fitting the theory to the linear part of the PMF for ξ = 1, where PE chains are mostly stretched and the theory is expected to be most reliable. For all other ξ-values, the model now delivers a prediction. We see in Fig. 3(a) that for ξ ≤ 1, the discrepancy between the theory and the simulations grows in relative terms with decreasing ξ at the whole attractive range of r. This is expected as the assumption of a parallel rod-like sliding mechanism becomes less accurate for decreasing ξ. In contrast, for ξ > 1 the linear trend of the analytical prediction agrees very well, even quantitatively, with the simulated PMFs for intermediate distances 2.5 nm ≲r ≲ 6.5 nm. Note that in this regime, the mean force is relatively constant, as assumed in the theory. At very small distances r ≲ 2 nm, the PE chains collapse into a globule and therefore the rod model clearly breaks down. This shows that the first steps of the PE association for highly charged PEs can be very well captured by the analytical model based on two rigid sliding rods together with counterion release entropy. Note that the differences between the PMF for ξ = 1 and the PMFs for ξ > 1 are solely provided by counterion-release entropy. The enthalpic Coulomb part (attraction of rods with renormalized charge) becomes less important with increasing ξ.

As a minor but interesting note, we now estimate the relative time period of the association process into the final complex with respect to ordinary diffusion. The typical diffusion time over a length L in a simple Rouse picture⁶⁷ would be τ_D ≃ N_bξ_mL²/k_BT, where ξ_m is the friction constant of a monomer. This has to be compared to a macromolecule associating a length L with a speed L/τ_assoc = f/(N_bξ_m) under the influence of a driving force f. Comparing the time scales, we obtain simply τ_D/τ_assoc ≃ βfL. Plugging in our calculated values of the mean force (insets to Fig. 3) for L on a nanometer scale we see that, for not too small Manning parameters ξ, the association for electrostatic and counterion-release driven PE complexation can be easily 1–2 orders of magnitude faster than a simple diffusion-dominated association process. These short time scales have indeed been observed in coarse-grained computer simulations of unrestrained PE complexation.^34,37,40

In order to further corroborate our theoretical assumptions on the counterion release effect, we now examine the number of released ions during the association process. Figure 4 presents the number n(r) of released counterions per PE (averaged over both chains) for various ξ values resolved in COM distance r. In the coil phase, where both PEs are independent of each other, the amount of condensed counterions reaches its maximal value. As mentioned before, we define the Manning radius r₀ for each value of ξ such that the theoretical prediction for N = N_bz_b(1 − 1/ξ) can be reproduced for an isolated rod. At the critical distance r^∗, the number of condensed counterions experiences a discontinuous jump ΔN. Here, the coil phase goes over to the elongated phase where the end parts of the PEs stick together and by that trigger the release of ΔN counterions. The latter increases with ξ, where in the case of ξ = 2.31, more than 3 counterions are liberated per PE. After the jump, a linear decrease in the number of condensed counterions demonstrates the progressive release of counterions from the overlapping segments of PEs, as predicted by our theory, Eq. (14). Finally, all counterions are released in the final complex at r = 0, where both PEs completely neutralize each other. Here, it can be seen that between 2 and 15.5 ions are released per PE into the bulk for the span of ξ between 1.15 and 2.31. Assuming approximately 5 k_BT per ion (see methods after Eq. (16)) in the case of ξ = 2.31, we end up with a total free energy of complexation about 155 k_BT, in good agreement with the PMF data in Fig. 3(b).

The discontinuities in the mean force, condensed counterion number N(r), as well as the end-to-end distance R_ee clearly suggest that the system undergoes an abrupt change at the critical distance r^∗ of the approaching PE chains. In an attempt to characterize this transition in more detail, we have calculated the probability distribution P(R_ee) of the PE end-to-end distance for various fixed (harmonically constrained) COM distances r in the range from 6.8 nm to 7.2 nm for the case of ξ = 2.31, thereby crossing its critical distance r^∗, shown in Fig. 5(a). For the two largest distances as well as for the smallest distance r, we find single-peaked distributions, corresponding to the well defined single states, coil and stretched PE chains, respectively. At around r^∗ = 6.9–7.0 nm, however, a bimodal distribution between the two states appears, indicating a structural coexistence, separated by unlikely states. The transition free energy profiles F_trans(R_ee) = − k_BTln P(R_ee) are plotted in Fig. 5(b). In the bimodal states, the potential barrier separating the two states has a height of the order of ΔF_trans ∼ k_BT. Note also that the peak position of the extended state is located at values of about L_r ≃ 7.2 nm, somewhat larger than r^∗. The reason is that the stable state apparently needs some finite overlap, i.e., a critical attraction to warrant a stable state.

So, what is the reason for this bimodal distribution? Consider first the PE chains in the coil state, i.e., when r > r^∗. In a rare fluctuation, the chains stretch out, accompanied by a significant loss in conformational entropy, and may achieve their handshake by overlapping with one or a few more monomers at the ends. For ξ > 1, we have shown that this “first touch” will be accompanied by a significant release of counterions, contributing a large favorable entropy of several k_BT per released ion. Apparently, this gain in counterion entropy is large enough to compensate for the loss in the conformational entropy of PE chains such that a coexistence can be established in this restrained equilibrium. Quantitatively, we estimate the stretching entropy of a single PE from our steered Langevin simulations. Results are shown in Fig. 6(a) for various ξ values, including the neutral reference ξ = 0. For ξ > 1, the PE stretching is a bit easier than for a neutral polymer owing to the internal electrostatic repulsion. However, for example, for ξ = 2.31, we measure an appreciable entropy loss of about 10 k_BT per PE chain for the extension L_r ≃ 7.5 nm. This loss in the chain conformation entropy needs to be compensated by the handshake and forthcoming release of counterions in order to establish the coexistence between associated and free states. For this particular example, thus 20 k_BT are needed to be compensated in total for both chains. For ξ = 2.31, we estimated 5 k_BT per released ion, so that about four counterions were required. Glancing back at Fig. 4, however, we see that for ξ = 2.31 in total seven ions are released at the handshake. The apparent discrepancy can be reconciled by investigating the number of condensed and released counterions during stretching, see Fig. 6(b). Apparently, ions are released during stretching, very likely due to the increasing effective PE length (and decreasing charge density) during stretching. At the relevant stretched state at R_ee = L_r ≃ 7.5 nm for ξ = 2.31, a mean of 1.5 ions per chain (thus, three in total) is released. In other words, only four counterions are released in the actual handshake, in agreement with the needed compensation stretching penalty mentioned above. These subtle structural effects on a single-ion level are important for a quantitative interpretation of the transient states in complexation. The free energy barrier between the coil and extended state must then be clearly attributed to the entropy of intermediate stretching.

In summary, we have studied PE structure variations and the resultant PMF profiles along the PE–PE center-of-mass reaction coordinate for PE pair complexation with a focus on intermediate association ranges for various PE charge densities. For charge densities above the condensation threshold, we observed and analyzed in detail a (fast) sliding-rod-like process preceding the PE complexation. We introduced an abstract model leading to an analytical expression for the PMF. The latter predicts a PMF virtually linear in center-of-mass distance at the onset of complexation of the first “handshake” until collapsing into the final complex, in good agreement with the computer simulations. Furthermore, a detailed inspection of the mean force profile uncovered a discontinuity at the onset of complex formation, which is also embodied as a jump in numerous other simulation measures. We demonstrated that the discontinuity can be attributed to the presence of a free-energy barrier stemming from cooperative counterion-release effects and single PE stretching entropy. Perhaps even more complex transient behavior can be expected when tuning chain stiffness, chain length, counterion valency, and the strength of the electrostatic interaction,^35,39 or including the action of explicit solvent.^68–70 

Because of the drastic changes in the topology of polymer chains, we suspect that the metastable states that emerge preceding to PE complexation may be relevant for realistic biological systems. For example, polycations for gene delivery are designed to complex with nucleic acids such that the gene transfer vectors are able to overcome the intracellular barriers such as the plasma membrane, the endosome, and the nuclear membrane. Successful gene delivery depends upon the ability of the vector to adopt different metastable structures or possess certain properties at different intracellular environments.^3,6,22 Here, one could suspect that transiently stretched or coiled states could provide some function, possibly some disorder-order based signaling, in a specific environment. In that respect, it is also interesting to note that transient disorder-to-order transitions, coupled with the adoption of different structures with different partners, are a common feature of intrinsically disordered proteins (IDPs)⁷¹ whose capacity for binding diversity plays important roles in both protein-protein interaction networks and likely also in gene regulation networks.

We finally note that the discontinuity observed in the PE–PE association process may also emerge during the complexation of PEs and globular proteins as indicated in recent coarse-grained computer simulations of more⁴⁷ or less⁴⁶ resolved protein models with heterogeneous charge distributions. In these examples, counterion-release effects after PE binding to protein surface charge patches have been identified as the main interaction responsible for the PE–protein complexation.

X.X. acknowledges funding from the Chinese Scholarship Council (CSC). Research at UCR is financially supported by the U.S. National Science Foundation (Grant No. NSF-CBET-0852353). M.K. and J.D. acknowledge funding from the ERC (European Research Council) within the Consolidator Grant with Project No. 646659–NANOREACTOR.