Multi-valent Ion Mediated Polyelectrolyte Association and Structure

Polyelectrolytes are commonly used to chelate multi-valent ions in aqueous solutions, playing a critical role in water softening and the prevention of mineralization. At sufficient ionic strength, ion-mediated polyelectrolyte--polyelectrolyte interactions can precipitate polyelectrolyte--ion complexes, a phenomenon known as"like-charge attraction". While the significant influence of small ions on polyelectrolyte solution phase behavior is recognized, the precise molecular mechanisms driving the counterintuitive phenomenon remain largely elusive. In this study, we employ all-atom molecular dynamics simulations to investigate the molecular mechanism of like-charge attraction between two poly(acrylic acid) (PAA) chains in solution. We find that moderate quantities of Ca$^{2+}$ ions induce attraction between PAA chains, facilitated by the formation of PAA--Ca$^{2+}$--PAA bridges and a significant increase in the coordination of Ca$^{2+}$ ions by the PAA chains. At high Ca$^{2+}$ number densities, ion bridges are disfavored due to electrostatic screening, yet the chains are still attracted to each other due to solvent-mediated interactions between the chains and their chelated ions. The insights gleaned from this study not only enrich our understanding of the intricate mechanism of like-charge attraction between polyanions in solution but also illuminate the influence of multi-valent ions on polyelectrolyte interactions.


Introduction
Aqueous polyelectrolyte (PE) solutions find widespread utility in diverse fields, including water treatment, 1 drug delivery, 2 and scale (CaCO 3 ) inhibition, 3 among others.For many of these applications, it is crucial to control the behavior of the polyelectrolyte in the solution phase.5][6][7][8] Therefore, a detailed understanding of the relationship between polyelectrolyte structure and PE-ion interactions is essential for the systematic design of advanced polyelectrolyte materials and additives.
][13] The effectiveness of PAA, and other scale inhibitors, relies heavily on the polymer's ability to chelate many Ca 2+ while preventing precipitation of the polyelectrolyte and chelated ion complex.Past experimental studies, such as those by Huber, 14,15 have demonstrated that the addition of Ca 2+ induces a transition in polymer conformation from an extended coil to a collapsed chain, which results in the precipitation of PE-ion complexes once the polymer's binding capacity is exceeded.Subsequent work conducted by Sinn et al. 9 found that a screened-Coulomb ion-polyelectrolyte interaction model could not explain the observed precipitation phenomenon, suggesting a more detailed understanding of Ca 2+ mediated interactions is necessary.
][18][19][20][21][22] Coarse-grained theoretical models have captured the addition of multivalent ions leading to chain collapse.4][25][26] While theory and experiments predict precipitation of polyelectrolyte-Ca 2+ complexes at certain conditions, to design novel polyelectrolytes that stay soluble in aqueous solution, we need to understand the molecular principles that govern the Ca 2+ induced association between likecharged polyelectrolytes.Molecular dynamics (MD) simulations provide a framework for such an investigation.
Several seminal MD studies have investigated the behavior of polyelectrolytes in aqueous CaCl 2 solutions.The work of Molnar and Rieger 27 provided evidence for the attraction of polyanions in solution by showing that two PAA chains in a solution were more prone to association as the number of Ca 2+ ions increased.However, the limited simulation box size (6 nm) and integration time (10 ns) precluded observation of different ion binding en-vironments, polymer conformations, and transitions between states.Subsequent single-chain PAA simulations by the Parrinello group 28,29 corroborated the stability of ion bridging and hypothesized that interchain ion bridges were responsible for the observed attraction between chains.
While these prior studies have yielded valuable insights, our preceding work 30 underscored a critical limitation: using classical force fields to model electrostatic interactions involving Ca 2+ without accounting for polarization effects results in exceedingly long Ca 2+ -PAA relaxation times.Such models can erroneously predict charge inversion of PAA-Ca 2+ complexes, a phenomena that is not supported by experimental data. 9,10This discrepancy implies that earlier MD investigations concerning Ca 2+ mediated like-charge interactions could suffer from insufficient sampling of PAA-Ca 2+ populations and polyelectrolyte conformations.
To overcome these challenges, our previous work 30 utilized the Electronic Continuum Correction (ECC) method, modeling the electronic polarization effects in a mean-field manner.This approach enabled us to mirror the polyelectrolyte-ion binding capacities observed in experiments. 9,10We also calculated the adsorption isotherm for Ca 2+ ions on a single PAA chain.The findings showed a transition from an extended coil to a collapsed chain conformation as the number of Ca 2+ ions increased, a result consistent with experimental observations. 10,14e attributed this transition to the formation of intrachain ion bridges which caused a chain collapse.Despite these advancements, there remain unresolved questions.For example, it is unclear how Ca 2+ ions mediate interchain interactions, and whether the formation of interchain ion bridges constitutes the dominant mechanism driving like-charge attraction between PAA chains in solution.
In the present study, we use our previously developed all-atom MD model 30 to investigate the molecular mechanism of like-charge attraction between two PAA chains in solution.Our enhanced sampling protocol provides an efficient exploration of the polyelectrolyte conformational space and ion binding environments, enabling the identification of the dominant conformations of PAA chains in solution.We find that at moderate Ca 2+ numbers, the PAA chains are attracted to each other due to the formation of direct PAA-Ca 2+ -PAA bridges between the chains as well as a significant increase in the number of ion and polymer contacts.However, at high Ca 2+ numbers, direct ion bridges between chains are disfavored due to electrostatic screening, yet the chains are still attracted to each other mainly due to solventmediated interactions between chelated ions on one chain and the carboxylate groups on the other chain.We then compute the precipitation concentration of PAA in solution using the second osmotic virial coefficient to quantify the net attraction between PAA chains.To analyze the dominant conformations of PAA chains interacting in solution and to investigate potential transition pathways between these conformations, we employ machine learning techniques.
The rest of the manuscript is organized as follows.We first describe our models and specific parameters used in the enhanced sampling algorithms for the molecular dynamics simulations.Next, we report the numerical data obtained from our calculations and discuss their implications.We then conclude the article with an outlook on the path forward for molecular design of antiscalant polyelectrolytes.

Methods
Molecular Dynamics Simulations.Molecular dynamics simulations were performed using GROMACS (version 2022.3)5][36] The PAA chains were modeled as atactic and fully deprotonated (pH ≥ 11) with Na + counter-ions, and constructed using CHARMM-GUI 37,38 with 16 monomers per chain.The PAA chains were then solvated in a 12 nm cubic box of SPC/E water 39 with Packmol. 40he general AMBER force field (GAFF) 41,42 was employed to model the PAA chains, following the protocol used in our single-chain PAA studies 30 and originally validated by Mintis and Mavrantzas. 43As electrostatic interactions are known to be overly strong in non-polarizable force fields, 44 we utilized the electronic continuum correction (ECC) method 45,46 to more accurately describe the electrostatic interactions between PAA and Ca 2+ ions.The ion parameters for Ca 2+ and Cl − were taken from the electronic continuum correction with rescaling (ECCR) parameters of Martinek et al., 44 while those for Na + were taken from the ECCR optimized parameters of Kohagen et al. 47 It has been shown that scaling polyelectrolyte and small ion charges is required to reproduce experimental binding results accurately. 30,48As a result, we applied the ECC method to all PAA partial charges.
For the van der Waals interactions, a cutoff of 1.2 nm was chosen.The long-range electrostatic interactions were calculated via the PME method 49,50 with a real space cutoff of 1.2 nm.The LINCS algorithm 51 was utilized to constrain all bonds with hydrogen atoms, and the system equations of motion were integrated using the leap-frog algorithm.
After solvation and ion addition, the system was energy-minimized using the steepest descent algorithm for approximately 100,000 steps.The system was then equilibrated in the NVT ensemble at 300 K for 10 ns with 1 fs time steps, using the Nosé-Hoover thermostat 52,53 with a 0.25 ps relaxation time constant.Further equilibration was carried out for another 10 ns in the NPT ensemble at 300 K and 1 bar using the Parrinello-Rahman barostat 54 with a 5 ps relaxation time constant.
The production phase simulations were conducted in the NVT ensemble at 300 K using a time step of 2 fs.The sampling of the equilibrium ensemble was enhanced by applying a combined well-tempered metadynamics 55-57 and HREMD 58,59 protocol.The well-tempered metadynamics protocol was used for improved sampling of the two-chain distance free energy surface (FES). 56,57Though metadynamics sampled the two-chain distance effectively, the slow relaxation of the polymer conformations and ion binding environments prevented adequate convergence of the FES.Consequently, we em-ployed the HREMD protocol to efficiently sample the polymer conformational space 60 and ion binding environments. 61Concurrently, the welltempered metadynamics protocol was used for improved sampling of the two-chain distance free energy surface (FES). 56,57oordinate exchanges between replicas were attempted every 100 steps, and the number of replicas was set to ensure an average exchange acceptance rate of approximately 25%.The systems with 0, 8, 16, and 32 Ca 2+ had 16 replicas, while those with 64 and 128 Ca 2+ had 24 replicas.In each replica, the electrostatic, Lennard-Jones, and dihedral interactions for all solute atoms were scaled by a parameter λ, ranging from 1 (unbiased) and 0.68 (most biased) geometrically distributed to improve exchange acceptance rates. 59ach replica contained an independent welltempered metadynamics bias potential with an initial Gaussian height of 1.2 kJ/mol, Gaussian width of 0.025 nm, deposition stride of 500 steps, and bias factor of 10.The metadynamics biases of each replica were evaluated in the Metropolis coordinate exchange probability.The collective variable (CV) in metadynamics was the distance between the centers of mass of the two PAA chains, facilitating the sampling of various associated and dissociated states.An additional harmonic wall was placed at 6 nm to prevent sampling beyond the minimum image.Replicas were equilibrated for 25 ns at 300 K before production sampling, and the metadynamics bias was added.Each replica simulation was run for at least 250 ns, yielding a total production simulation time of 29.4 µs.
The reported data in our study were collected from the λ = 1 replica, and statistics were calculated by reweighting this data based on the Boltzmann factor of the metadynamics and harmonic wall biases.Uncertainties, when reported, were estimated using block averaging on correlated MD data to obtain the standard error of the sample mean, which was then converted to 95% confidence intervals.Analysis of the simulation data involved the use of MD-Analysis 62,63 and custom Python scripts, while VMD 64,65 was employed for visualizing the generated trajectories.
Stability of Polyelectrolyte Solution.Determination of the solution phase boundary from molecular simulation is non-trivial, as phase transitions are macroscopic events, and the range of system configurational space is large.
This boundary typically comprises a polymer-dilute branch and a polymer-rich branch interconnected at the critical concentration, and as shown below, the spinodal of the dilute branch can be estimated with a leadingorder expansion of the osmotic virial equation of state.In contrast, the branch corresponding to higher polymer concentrations would require many chain simulations and the enhanced sampling approaches detailed above would become computationally prohibitive.Despite the difference in polymer concentrations and the complexities involved in their prediction, one can anticipate similar key physics and interactions on both branches of the spinodal.This is due to the fact that precipitation is primarily driven by local interactions between the polymer chains and the ions.Thus, regardless of the concentration, the essential characteristics of these interactions remain the same.
To estimate the spinodal from the simulation data, the primary property of interest from the simulations is the potential of mean force (PMF) between the two PAA chains as a function of the center of mass separation distance.
The PMF provides a local measure of the relative stability of associated and dissociated chain configurations.If the PMF at narrow separations is negative, the chains favor association, which may lead to precipitation of the polyelectrolyte at sufficiently high concentrations.
Following previous studies, [66][67][68] we can then calculate the second osmotic virial coefficient as where r denotes the center of mass distance between the two PAA chains, β = 1/k B T , and U PMF (r) is the PMF obtained from the simulation.Rigorously, the PMF should be calculated at fixed chemical potential of the ions, but constant chemical potential simulations are computationally prohibitive for the systems stud-ied here.However, the PMF calculated at fixed number density of ions (canonical ensemble) has been shown to reproduce experimental trends in phase behavior for similar macromolecules with added salt. 67 positive B 2 (T ) indicates a repulsive interaction, while a negative B 2 (T ) indicates an attractive interaction.To calculate the maximum polyelectrolyte density of the suspension before the onset of polyelectrolyte-ion precipitation, we utilized the leading-order term in the virial expansion of the osmotic pressure 69,70 to find the spinodal of the polyelectrolyte solution as The spinodal polyelectrolyte chain concentration allows us to quantitatively determine the effects of Ca 2+ ions on the polyelectrolyte solution phase behavior and evaluate the emergence of like-charge attraction.
Autoencoder Neural Network.We employed machine learning models to analyze the large amount of simulation data generated and to map the polyelectrolyte conformational space in a low-dimensional representation.Autoencoder (AE) networks, a type of neural network, have emerged as powerful tools in understanding and predicting complex systems.
2][73][74] An AE network consists of an encoder that maps the input to a lower-dimensional representation and a decoder that reconstructs the original input from the lower-dimensional representation.Both networks are symmetric, containing the same number of layers and neurons.The encoded lower-dimensional representation, known as the latent space, provides a compressed representation of the input data and aids in visualizing the phase space.
In line with the approach proposed by Bandyopadhyay and Mondal, 73 we trained the AE network using the pairwise distances between the C α atoms of the PAA backbone.The encoder was built with three hidden layers com-prising 496, 128, and 32 neurons, fully connected with a latent space output of 2 dimensions.The training was conducted over 150 epochs using the Adam optimizer with a learning rate of 0.001 and batches of 256 randomly chosen conformations.We used the mean squared error between the input and output as the loss function and applied an L 2 penalty of 0.00001 for weight regularization.Furthermore, to test the model, 20% of the data was withheld.All aspects of model development and training were facilitated by the PyTorch library. 75,76

Results and Discussion
The 16-mer PAA two-chain PMF curves for each number of Ca 2+ ions (N Ca 2+ ) are shown in Figure 1.The reference state for each system is the free energy at r = 4 nm, where the average two-chain interactions have plateaued for all calcium-containing systems.In the absence of Ca 2+ ions, the system has repulsive (∆F > 0) interactions at all distances, as expected.At small numbers of Ca 2+ ions (N Ca 2+ ≈ 8), Ca 2+ adsorption reduces the long-ranged chain repulsion and a metastable association well emerges at ∼1.2 nm.As N Ca 2+ increases, the PMF develops a globally stable well at short distances (r ≤ 1.8 nm) with a large barrier separating the associated and dissociated states.Higher N Ca 2+ shifts the PMF well to larger distances and the barrier height vanishes.Increased electrostatic screening at higher N Ca 2+ reduces the long-ranged repulsion between the chains, and likely contributes to the observed reduction in the barrier height.We next estimated the precipitation conditions of the polymer-ion complex using the osmotic virial equation of state.The second osmotic virial coefficient (B 2 ) physically represents the two-body interaction between polymer chains such that a positive B 2 value indicates repulsion, and a negative B 2 value indicates attraction.We calculate B 2 using Eq. 1 with the approximation that the infinite integral domain can be replaced with a finite domain of r ∈ [0, 6] nm, which is appropriate given the short-ranged interactions and the PMF curves plateau, which causes the integrand to vanish.We do not include the N Ca 2+ = 0 sys- tem in the calculation of B 2 because the PMF is repulsive at all distances.For systems with net-attraction (B 2 < 0), the corresponding maximum polymer concentration at which the polymer-ion complex is predicted to precipitate is calculated using Eq. 2 and shown in the inset of Figure 2. Interestingly, the system with 16 Ca 2+ ions exhibits a positive B 2 value, suggesting no precipitation of the polymer-ion complex (to leading order), despite possessing an attractive PMF well of ∼4k B T .We additionally observe a non-monotonic trend in the second osmotic virial coefficient.This trend results in the spinodal polymer concentration be-ing higher for 128 Ca 2+ than 64 Ca 2+ , which implies a salting-in effect.The non-monotonic behavior of the second osmotic virial coefficient with increasing N Ca 2+ qualitatively agrees with theoretical predictions of salting-in effects reported by Wittmer et al. 25 at higher ionic strengths.However, unlike the theoretical prediction of polymer charge inversion via ion chelation, our simulations did not show such behavior.Instead, the number of calcium ions bound to the polymer increased and saturated below the 16 Ca 2+ count needed for charge neutrality (see Figure 3).The number of calcium ions bound to the two 16-mer PAA chains saturates around 9 ions, which is not enough to neutralize, much less invert, the charge of the polymer.These results are consistent with our single-chain 32-mer PAA simulations, 30 which showed that the number of calcium ions bound to the polymer saturates at approximately 0.4 ions per monomer.Experimental studies on longer chains have similarly measured a binding capacity of about 0.3 Ca 2+ per monomer. 9,10t high N Ca 2+ , the electrostatic screening reduces the formation of ion bridges between chains, leading to a slightly weaker and longer-ranged attraction, which is consistent with our single-chain PAA simulations 30 and de la Cruz et al . 16This phenomenon is the source of the salting-in behavior.The average number of bridging Ca 2+ ions approaches zero at 128 Ca 2+ ions, yet the chains are still attracted to each other, as seen in the PMF curves (Figure 1).However, the chain-associated state remains favorable due to the presence of 9 single-chain adsorbed Ca 2+ ions.These single-chain adsorbed Ca 2+ are often chelated by a single carboxylate group, which locally inverts the effective monomer charge from −1 to +1.The positive charge of the effectively monovalent cation monomer screens the long-ranged chain repulsive interactions and promotes chain association.
The depiction of Ca 2+ adsorption in Figure 3 elucidates how Ca 2+ facilitates the association of PAA chains.However, this representation only shows the structure of associated states and not the overall driving forces propelling the chain association.Our single-chain studies 30 have shown that isolated chains similarly adsorb Ca 2+ and exhibit intrachain ion bridging.We observe that as the two chains approach one another, the calcium-mediated interactions become more pronounced through both an increase in ion bridging and calcium adsorption, as shown in Figure . 4.
The above trend can be quantified by evaluating the number of contacts between carboxylate groups and Ca 2+ ions.We establish a 'contact' when the distance between a carboxylate group and a Ca 2+ ion falls below 0.35 nm.This cut-off was determined by the location of the first minima in the radial distribution function for Ca 2+ and carboxylate oxygen atoms.The proximity denotes direct PE-ion interactions, which eliminate solvent-mediated interactions.
Figure 4 shows the free energy landscape of PE-ion contacts as a function of the interchain distance for systems with 32 and 128 Ca 2+ ions, respectively.We focus on these two systems, as 32 Ca 2+ ions is the minimum number of ions in this study that leads to B 2 < 0 and 128 Ca 2+ ions exhibits non-monotonicity in the B 2 curve (Figure 2).For 32 Ca 2+ ions, the stable associated states have a larger number of PE-ion contacts than the dissociated states.As the chains approach each other, more Ca 2+ ions adsorb onto the chains due to the higher density of carboxylate groups.This not only decreases the electrostatic repulsion between the chains but also aids in forming ion bridges.However, for 128 Ca 2+ ions, the chains are already saturated with Ca 2+ ions in the dissociated states, and so the relative increase in the number of PE-ion contacts is less with further decrease in the interchain distance.In addition, the overall free energy surface (FES) valley is shallower and shifted to larger chain center of mass distances, which is consistent with the decreased number of bridging Ca 2+ ions (Figure 3).
We sought to identify the dominant polymer conformations of the most stable associated states in the 32 Ca 2+ system that contribute to the two-chain attractive interactions and to explore the relative importance of bridging and single-chain adsorbed Ca 2+ ions.Traditional collective variables, such as the radius of gyration or the angle between principal radius of gyration vectors, proved insufficient to discern the dominant conformations due to the multitude of metastable states and the broad distributions of these variables.To address these challenges, we employed machine learning techniques, specifically an autoencoder, to learn a low-dimensional representation of the conformational space.
Figure 5 illustrates the two-dimensional latent space projection (z 1 , z 2 ) of the input data, which is reweighted by the Boltzmann factor of the metadynamics bias potential and subsequently transformed into a free energy surface.The reference state for the free energy surface was set as the most stable minima within the latent space and conformations near relevant minima are rendered.The chains are colored to indicate the conformations more clearly, and ions within 0.35 nm of the chains are shown.Water molecules are omitted for clarity.Remarkably, the autoencoder identifies different ion bridging environments without explicit training on the number of bridging ions or their coordinates.
The most stable conformations reside in the upper energy basin and contain 4 bridging Ca 2+  with 6 single-chain adsorbed Ca 2+ .Within this arrangement, one of the chains adopts a collapsed conformation (colored red) and is partially wrapped around the other chain, which is in an extended conformation (colored blue).
The collapsed chain facilitates the formation of ion bridges with the extended chain, and each chain adsorbs three additional Ca 2+ ions, which serve to neutralize the polymer charge and mitigate unfavorable carboxylate-carboxylate interactions.
The central energy basin hosts a more diverse set of conformations, stabilized by 3 bridging Ca 2+ ions.Chains in this region can assume either extended or partially collapsed conformations.The most stable conformation maximizes the number of ion bridges, but interestingly, the autoencoder also identifies conformations with fewer ion bridges that are stabilized by extra single-chain adsorbed Ca 2+ ions.Both metastable basins at 0.18 and 0.31 k B T contain 8 single-chain adsorbed Ca 2+ .Additional frames depict local minima along the transition path between metastable conformations, as well as a possible transition between the basins of 3 and 4 bridging calcium ions.
The basin with 3 bridging calcium ions at 0.74 k B T contains 7 single-chain adsorbed Ca 2+ and highlights the relative importance of ion adsorption and bridging at moderate N Ca 2+ .Com-pared to the minimum free energy conformation observed, the total number of Ca 2+ on the polyelectrolyte complex is the same in both conformations, but a bridging calcium ion has been replaced by a single-chain adsorbed calcium ion.The 0.74 k B T basin also has a relatively extended chain conformations, which likely yields an increase in the chain conformational entropy.The loss of the ion bridge has a modest impact on free energy (< 1 k B T ), indicating that the ion bridging does not dominate the free energy landscape.The cost of losing an ion bridge can be further reduced to 0.18 k B T by the addition of 1 single-chain adsorbed calcium ion.
As seen in Figure 3, ion bridges are not dominant at N Ca 2+ = 128.However, the presence of ion bridges does not vanish, as seen in Figure 6.The left panel depicts a conformation with 4 ion bridges and the right panel depicts a conformation with no ion bridges that largely interacts through via solvent mediation.Both conformations are within 0.1 k B T of the minimum in the interchain potential of mean force and are highly populated.However, the majority of conformations in the basin contain no ion bridges, as the increased electrostatic screening reduces the cost of separating the chains and leads to a slightly longer-ranged interaction well (Figure 1).We observe that the increased flexibility of polymer conformations with more Ca 2+ ions allows for a broader range of conformations, where ion bridges can be destroyed and additional ions adsorbed with almost no free energy difference.The conformational flexibility may play a role in the salting-in behavior, as the chains can relax into conformations that are more favorable for solvation.

Conclusions
This investigation has uncovered the fundamental factors underpinning the like-charge attraction observed within polyanion solutions mediated by multi-valent cations.We have found that two primary mechanisms contribute to this phenomenon.
The initial driving force is the increase in the number of favorable polyelectrolyte-ion contacts in chain-associated states when compared to their respective dissociated states (Figure 4).At lower Ca 2+ numbers, these contacts lead to ion bridges between the chains, which are important in setting the length-scale of the attraction.This is manifested in the position of the PMF minima in Figure 1.
Our autoencoder analysis revealed that once a sufficient number of ion bridges are formed (3 for the 32 Ca 2+ system), the free energy difference between forming and breaking additional ion bridges is small compared to thermal energy.Notably, the ion bridge stability is contingent on electrostatic screening, with higher Ca 2+ numbers rendering these bridges unnec- essary for attraction.At high ionic strengths, chains are saturated with Ca 2+ ions and increases in electrostatic screening decrease the favorability of interchain ion bridges, yet the chains still experience a net attraction.This attraction remains similar in magnitude even when the relative increase of PE-ion contacts in chain-associated states is less pronounced.
The observed attraction at high Ca 2+ numbers hints at a second factor driving the attraction: favorable chain-chain interactions that are enabled by the PE-ion contacts on individual chains.We observed that as the two chains approach one another, they increase the number of PE-ion contacts to stabilize the associated chain states.
Our study is consistent with earlier theoretical findings, 16,[18][19][20] which propose that divalent ions can induce chain association.However, we propose that precipitation might not be dominated by ion bridging induced chain collapse or chain neutralization, as generic polyelectrolyte models may suggest.As the number of Ca 2+ increases, chains are saturated with ions, but the saturated PE-ion complex still carries a net-negative charge (Figure 2).At high Ca 2+ numbers, electrostatic screening decreases the favorability of direct ion bridging, yet the chains still experience a net attraction due to solvent-mediated interactions between the chains and chelated ions.
Our investigation into chain association builds on the previous research of the Parrinello group 28,29 , focusing on longer chains that can accommodate a wider range of Ca 2+ adsorption environments, such as ions chelated by multiple, non-neighboring monomers.Figure 5 illustrates the relative favorability of ion bridging, and our observations show that the most stable associated state arises when at least one of the chains adopts an extended conformation, with ion bridges forming exclusively between chains.Even minor changes in the relative polymer conformations can significantly alter the ion binding environments and the favorability of chain association, as shown in the upper basin of Figure .5.
Looking forward, potential strategies for maintaining aqueous polyelectrolyte solution stability could involve the introduction of bulky or non-polar side-chains at lower ionic strengths to curtail the short-ranged ion bridging.At higher ionic strengths, solution stability might be achieved by decreasing the pH to protonate carboxylate groups and subsequently the number of adsorbed ions.The competition between interchain ion bridging and single-chain ion chelation elucidates a more comprehensive understanding of like-charge attractions in polyanion solutions.
(   the ion binding/bridging, they provide a local measure of the interchain interaction distance.
The carboxylate carbon RDFs show a well-defined peak at 0.5 nm that is representative of the distance between carboxylate groups bridged by a single Ca 2+ ion.Notably, the height of this peak shows the same non-monotonic trend as the number of bridging Ca 2+ ions seen in the main text with increasing Ca 2+ ions.Ions bound to a single chain are green, and bridging ions are blue.The 0 Ca 2+ data comes from associated states where the interchain distance was less than 2 nm, as the PMF was repulsive at all distances.The dotted lines are guides to the eye.

Additional Free Energy Surfaces
We show the free energy surfaces for the number of calcium and carboxylate oxygen contacts as a function of the interchain center of mass distance for systems not shown in the main text in Figure 5.The 8 Ca 2+ system favors dissociated chain states, as there are not enough Ca 2+ ions to form sufficient numbers of polymer-ion contacts to stabilize the two-chain complex.
The 16 and 64 Ca 2+ systems exhibit a qualitatively similar free energy surface to the 32 Ca 2+ system in the main text.
8 Ca 2+     To gain insight into the relationships between the latent space features and relevant physical quantities influencing the conformational space, we analyzed the linear correlation between a few selected variables in Figure 8.The latent space features exhibit weak correlations with each other, implying that they capture distinct and complementary aspects of the system.This observation suggests that the autoencoder effectively performs dimensionality reduction and representation learning, allowing it to extract relevant features while preserving important information about the system.The moderate correlation with polymer radii of gyration is consistent with our previous single-chain studies, highlighting the significance of polymer size in determining calcium adsorption behavior.The relatively weak correla-within 1 k B T of the minimum in the interchain PMF (0.70-1.07 nm), focusing on specific conformational states near the minimum.

Figure 1 .
Figure 1.Potential of mean force for two 16-mer PAA chains with varying number of Ca 2+ ions.

Figure 2 .
Figure2.PAA second osmotic virial coefficient as a function of Ca 2+ ions in the system.Inset shows the corresponding polymer concentration at which the polymer-ion complex is predicted to precipitate.The dotted lines are guides to the eye.

Figure 3 .
Figure3.Number of calcium ions bound to PAA carboxylate groups as a function of the number of calcium ions in the system.Data are plotted within 1 k B T of the minimum in the interchain potential of mean force, ensuring the relevance to the system's most stable configurations.Calcium ions bound to a single chain are depicted in green, while bridging ions between chains are shown in blue.The dotted lines are guides to the eye.

Figure 4 .
Figure 4. Two-dimensional free energy surfaces for two 16-mer PAA chains with 32 (left panel) and 128 (right panel) Ca 2+ ions.The horizontal axis represents the center of mass distance between the PAA chains, while the vertical axis denotes the summed coordination number of carboxylate oxygen atoms about Ca 2+ ions.Isolines are drawn at 1, 2, 4, and 6 k B T .

Figure 5 .
Figure 5. Free energy landscape of the autoencoder latent space for two 16-mer PAA chains with 32 Ca 2+ ions within 1 k B T of the minimum in the interchain potential of mean force.Chain conformations near relevant minima are visualized, along with ions within 0.35 nm of the chains.The tuple below each conformation displays the corresponding free energy, the number of bridging calcium ions, and the number of single-chain adsorbed calcium ions.Isolines are drawn at 1, 2, and 4 k B T .

Figure 6 .
Figure 6.Representative conformations of two 16mer PAA chains with 128 Ca 2+ ions (within 0.1 k B T of the minimum in the interchain potential of mean force).The visualization depicts the two PAA chains as dark blue and red, Na + ions as magenta, Ca 2+ ions as yellow, and interchain bridging Ca 2+ ions as light blue.Left panel: The chains are in direct contact, bridged by 4 Ca 2+ ions, and have 5 single-chain adsorbed Ca 2+ ions.Right panel: The chains have no bridging Ca 2+ ions but feature 12 single-chain adsorbed Ca 2+ ions.Note that the red chain's principal axis is aligned into the page for the right panel.
Alec Glisman, †, § Sriteja Mantha, †, § Decai Yu, ‡ Eric Wasserman, ¶ Scott Backer, ¶ and Zhen-Gang Wang * , † †Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA ‡The Dow Chemical Company, Core R&D, 633 Washington St., Midland, Michigan 48674, USA ¶The Dow Chemical Company, Consumer Solutions R&D, 400 Arcola Road, Collegeville, Pennsylvania 19426, USA §Contributed equally to this work E-mail: zgw@caltech.eduPolyelectrolyte-IonStructureIn this section, we present additional details on the structure of the polyelectrolyte-ion complexes for the 16-mer poly(acrylic acid) (PAA) systems.The potential of mean force (PMF) curves as a function of the interchain distance shown in the main text are coarsegrained metrics of the interchain interactions and do not provide information about the local structure of the complexes.In Figure1, we present the radial distribution function (RDF) of side-chain carboxylate carbon atoms, denoted C cb , on one chain (A) with respect to C cb 1 arXiv:2311.10914v1[cond-mat.soft]17 Nov 2023 atoms on the other chain (B).As the C cb atoms are on the side-chains and directly involved in

Figure 2 .
Figure 2. Number of sodium ions bound to PAA carboxylate groups as a function of the number of calcium ions in the system.Data plotted within 1 k B T of the minimum in the interchain potential of mean force.Ions bound to a single chain are green, and bridging ions are blue.The 0 Ca 2+ data comes from associated states where the interchain distance was less than 2 nm, as the PMF was repulsive at all distances.The dotted lines are guides to the eye.

Figure 5 .
Figure 5. Two-dimensional free energy surfaces for two 16-mer PAA chains with 8 (upper left panel), 16 (upper right panel), and 64 (bottom panel) Ca 2+ ions.The horizontal axis denotes the chain center of mass distance, and the vertical axis denotes the summed coordination number of carboxylate oxygen atoms around Ca 2+ ions.Isolines are drawn at 1, 2, 4, and 6 k B T .

Figure 7 .
Figure 7. Top panel: Difference between the input and output of the AE for the pairwise distance between all backbone alpha carbons.Middle panel: Scaled input features to the AE.Bottom panel: Reconstructed features from the latent space.

Figure 7
Figure 7 provides a visual representation of the AE's performance on a set of 6 randomly selected conformations.The top panel illustrates the difference between the scaled input data and the reconstructed output of the AE.The middle panel displays the scaled input features to the AE, while the bottom panel presents the reconstructed features from the latent space, showcasing the capability of the AE to capture important features of the input data.

Figure 8 .
Figure 8.The correlation matrix of the latent space features z i with physical features of interest for the 16-mer PAA chains with 32 calcium ions within 1 k B T of the minimum in the interchain PMF.The physical features are the interchain distance r 12 , the chain radius of gyration R g , the total number of contacts between C α atoms, and the total number of contacts between Ca 2+ ions and carboxylate carbons, respectively.
Figure 1.RDF of carboxylate carbon (C cb ) atoms on one chain with respect to carboxylate carbon atoms on the other chain.
Ca 2+ outcompete Na + for binding to the carboxylate groups, and the number of Na + bound to the carboxylate groups decreases with increasing Ca 2+ , as shown in Figure2.Even at high Ca 2+ numbers, one Na + remains bound to the polyelectrolyte-ion complex on average and does not contribute significantly to the ion bridging behavior.
Figure 4. Left panel: RDF of sodium ions around carboxylate oxygen atoms.Right panel: RDF of sodium ions around carboxylate carbon atoms.