Study of Hydrophobic Clustering in Partially Sulfonated Polystyrene

Sep 23, 2016 - Illustration of the solvation shell around monomers/clusters of a coil/globule polymer in solution. The yellow spheres marked with numb...
0 downloads 18 Views 3MB Size
Article pubs.acs.org/Macromolecules

Study of Hydrophobic Clustering in Partially Sulfonated Polystyrene Solutions with a Systematic Coarse-Grained Model Ran Zhang* and Nico F. A. van der Vegt Eduard-Zintl-Institut für Anorganische und Physikalische Chemie and Center of Smart Interfaces, Technische Universität Darmstadt, Alarich-Weiss-Straße 10, D-64287 Darmstadt, Germany S Supporting Information *

ABSTRACT: Systematically coarse-grained (CG) force field models are state point dependent and therefore are usually not transferable to different chemistries and thermodynamic conditions. Fragment-based approaches for deriving CG polymer force fields increase the transferability of the models, but several challenges remain, in particular to describe conformational transitions of polymers driven by changes in the solvent environment or chemical composition of the chain. Herein, we describe a coarse-graining approach for deriving implicit-solvent CG polymer models, which can be used to study hydrophobic collapse transitions in aqueous solutions. On the basis of an earlier reported CG model for sodium poly(styrenesulfonate) (PSSNa) [Macromolecules 2012, 45, 2551−2561], we describe an approach to rescale nonpolar pair potentials of mean force based on a solvent-accessible-surface-area argument, thus allowing to reduce the fraction of charged monomers in the chain and to account for partial chain collapse. Based on a comparison with extensive detailedatomistic, explicit-solvent, simulations, it is demonstrated that the rescaling function for nonpolar interactions is chain length dependent. Extrapolation of the rescaling factor to long chains provides a new way to perform CG simulations of polymer solutions with moderately hydrophobic polyelectrolyte chains. Upon application to a 50%-charged PSSNa chain we observe pearl-necklace conformations in a cascade evolving manner with increasing chain length, in agreement with theoretical predictions for polyelectrolytes in poor solvent.

1. INTRODUCTION The properties of polyelectrolytes and other water-soluble macromolecules depend on their chemical structure and architecture, the concentration of (specific) salt, counterions, and cosolutes present in the solution. Properties, such as chain conformations and lower critical solution temperatures, depend on direct ion−macromolecule interactions as well as on interactions of ions or cosolutes with the first hydration shell of the macromolecule.1−3 While computer simulations provide useful information at the atomic level, detailed-atomistic simulations of these systems remain a huge challenge due to their high computational cost, which results, to major extent, from modeling the solvent molecules explicitly.4−14 Generic bead−spring polyelectrolyte models15−17 that replace the explicit solvent by a dielectric continuum do not account for specific, short-range interactions; therefore, alternative models are required which are computationally efficient but keep sufficient chemical information. These alternative, coarsegrained (CG) models can be obtained by systematic (bottom-up) coarse-graining methods in which the parameters of the CG model are derived from a fine-grained, e.g., atomistic, model of the same system.18−27 An implicit-solvent, CG model for sodium poly(styrenesulfonate) (PSSNa) chains with 100% sulfonation degree has been derived in a previous work by one of the present authors.26 This model successfully reproduces the local and © XXXX American Chemical Society

global chain conformations as compared with a short atomistic parent chain in water. An important feature of this CG model is that the bonded and nonbonded interaction terms have been developed independently with a fragment-based procedure.26,27 This procedure derives the nonbonded parameters of the CG interaction sites from fine-grained atomistic models of small molecules representing those sites. The bonded parameters for bond stretching, angle bending, and torsional rotations in the CG model are derived independently from the distributions of these degrees of freedom sampled by an isolated atomistic chain. Combining these bonded and nonbonded contributions provides a CG model which can potentially be used for transferability studies in which thermodynamic (e.g., salt28,29) or chemical (monomer properties) conditions are varied. Herein, we describe a way to improve the transferability of the model of Li et al., aiming to simulate PSSNa chains with lower sulfonation degrees. Partially sulfonated (sulfonation degrees smaller than 0.9) PSSNa chains belong to the category of “hydrophobic polyelectrolytes”.30,31 While still bearing the nature of polyelectrolytes, these PSS chains collapse in water due to the hydrophobic interactions between uncharged styrene monomers. PSSNa chains with a sulfonation degree f between Received: May 27, 2016 Revised: September 15, 2016

A

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

experimental results based on large-scale simulations with system sizes up to 40 nm. This comparison allows to validate not only the CG model of Li et al.26 but also the atomistic model on which it is based.

0.3 and 0.9 adopt cylindrical-shaped conformations in aqueous solvent or, more interestingly, a “pearl-necklace” structure, i.e., dense monomer clusters connected by extended chain strings.31 With the implicit-solvent CG model of Li et al.,26 the simulations however do not show evidence for these types of structures; instead, extended coil conformations are observed. This is most likely caused by the fact that hydrophobic clustering of uncharged beads cannot be described with the implicit-solvent nonbonded interactions of Li et al. These interactions correspond to pair potentials of mean force (PMF) calculated from fine-grained simulations with two atomistic fragment molecules in a pure-water box and, due to the lack of incorporating effective multibody terms, do not correctly describe effective bead−bead interactions in fluctuating inhomogeneous environments. This problem is common for implicit-solvent CG models, which inevitably fail to reproduce the global aggregation structures or assembly dynamics.32−36 Therefore, the predictive modeling of hydrophobic clustering effects needs to be addressed in developing implicit-solvent nonbonded interaction potentials by systematic coarse-graining methods. However, there exist several challenges in deriving an implicit-solvent CG potential, which realistically describes hydrophobic aggregation and collapse of polymer chains or chain fragments in aqueous media. On the one hand, the computational effort of introducing multibody potential terms may become unaffordable while, on the other hand, the strength of hydrophobic bead−bead pair interactions required to observe collapse depends on the degree of polymerization (DP) of the chain as will be illustrated later. This dependence on chain length has not been considered in recent simulation work of partially sulfonated PSS chains in water26,37 yet is essential in order to observe pearl-necklace-type structures in implicit-solvent simulations of long chains with systematically coarse-grained models. Moreover, the DP dependence of the hydrophobic potential makes it difficult to parametrize the nonbonded potentials between nonpolar beads from PMF calculations.26 We will herein describe a method to readjust the available pairwise nonbonded potentials,26 accounting for the average multibody nature of the problem in an approximate, analytical manner. The method relies on the idea that solvation of a nonpolar group is energetically penalized in proportion to its solvent-exposed surface area.38 The pair PMFs of ref 26 quantify the effective strength of hydrophobic interaction between two nonpolar solutes in an aqueous medium with no further solutes present. Therefore, they do not realistically describe the strength of interaction between two solutes in a larger hydrophobic aggregate in which their solvent exposure is smaller. Instead of developing a new coarse-graining model that employs multibody terms for nonbonded interactions, we keep a fragment-based pairwise parametrization approach and rescale the nonbonded (hydrophobic) interactions with a suitable factor that accounts for changes in solvent-exposed surface area. The rest of this paper is organized as the following: In section 2, we outline the theoretical foundation of potential rescaling. In section 3, we present the atomistic and coarsegrained models and illustrate how to use potential rescaling in hierarchical, atomistically detailed, and coarse-grained simulations. The results are presented and discussed in section 4. The conclusions of this work are summarized in section 5. In the Supporting Information, we present details of the CG model and discuss the accuracy of this model, originally developed by Li et al.,26 by comparing the conformational properties of long, dilute PSSNa chains with corresponding

2. THEORY A schematic picture explaining the multibody nature of the problem is given in Figure 1: for a coil-like “hydrophilic”

Figure 1. Illustration of the solvation shell around monomers/clusters of a coil/globule polymer in solution. The yellow spheres marked with numbers represent certain nonbonded monomers along the chain (with other parts of the chain omitted for simplicity), while the smaller gray spheres represent the solvent molecules. Blue dashed circles schematically show how the solvated surface area should be considered in an implicit-solvent CG model, while the red ones show how existing nonbonded potentials (pair PMFs) account for the solvated surface area in a hydrophobic clustering case.

polymer, nonbonded pair PMFs can be used in CG simulations because the process of monomer association frequently involves only two monomers. Following a bottom-up coarsegraining scheme, these pair PMFs are derived from atomistic explicit-solvent simulations. These simulations provide the free energy of interaction between a solvated nonbonded solute pair (projected on the intermolecular distance), i.e., a two-body, or pair, PMF, by means of suitable averaging over solvent degrees of freedom and over internal degrees of freedom of the nonbonded solute pair selected to represent the corresponding chemical groups on the PSS chain.26 However, for the hydrophobic case where clusters of more than two nonbonded monomers associate (globule in Figure 1), the above pair PMFs do not correctly describe the stability (free energy) of the resulting cluster. In this case, the application of pair PMFs incorrectly assumes that the free energy of a hydrophobic cluster (of any size) is proportional to a pairwise sum of twobody PMFs (each PMF corresponds to the free energy of a “two-bead cluster” in bulk aqueous environment). The free energy of the hydrophobic cluster, assumed by means of this approach, is too high because the solvent-accessible surface area of a monomer inside a hydrophobic cluster will in reality be smaller than the solvent-accessible surface area, which is effectively assumed when a pairwise sum of two-body PMFs is used as a measure of the cluster total free energy. Clearly, the use of pair PMFs in coarse-grained simulations of polyelectrolyte solutions will not provide a correct physical description of hydrophobic polyelectrolytes. In these systems, the balance of surface energy and electrostatic repulsion determines the stability of globular, cylindrical, or pearlnecklace-type structures.39−41 We here derive a rescaling procedure of the nonbonded pair potentials for the CG chain beads using a solvent-accessible surface area argument. We assume that the free energy of the polyelectrolyte chain contains an electrostatic component (selfB

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

where U0 can be considered as a constant representing the surface energy of the corresponding chain blocks in a coil state. However, with the pairwise nonbonded potentials of Li et al.,26 the surface energy of such a globule is proportional to U2 and contains contributions from separate solvation shells around any possible clusters of two nearest beads (Figure 2c):

repulsion between the charged ionizable groups) and a nonelectrostatic component, which is proportional to the solvated molecular surface area of the uncharged groups. Because the charged groups are hydrated independent of the conformation adopted by the chain (and electrostatic interactions are explicitly accounted for in the CG model), our attention focuses on the latter free energy component. For a globular polymer, this component is proportional to the number of uncharged beads at the surface38,39 (we assume that the CG polyelectrolyte chain can be described as a connected sequence of beads). A bead corresponds to (part of) the chemical repeat unit of the real chain (cf. section 3) and can thus be expressed as U ≈ u0R glob 2/ξ 2

U′ ≈

1 nCnU2 2

(4)

and the difference would be

ΔU ′ =

1 nCnU2 − U0 2

(5)

where n is the number of beads in such a cluster and Cn the average coordination number of a bead in a cluster of size n. As illustrated in Figure 2c, the surface energy now equals the number of such shells (1/2nCn) multiplied by the surface energy of each shell (U2). This illustrates how the original nonbonded potentials will function in the CG simulations in the case of hydrophobic clustering. In order to reach the actual final state described in eq 1, extra energy is needed:

(1)

where Rglob is the size of the globule and ξ the size of the bead (cf. Figure 2a) and u0 corresponds to the energy required to

⎛R 2 nCnR 2 2 ⎞ glob ⎟ ΔUe = U − U ′ ≈ u0⎜⎜ 2 − 2ξ 2 ⎟⎠ ⎝ ξ

(6)

Since the size Rglob of a spherical globule is proportional to ξn1/3, the expression of ΔUe can be simplified as ⎛ 22/3nCn ⎞ ⎟⎟ ΔUe ≈ u0⎜⎜n2/3 − 2 ⎠ ⎝

Figure 2. Thermodynamic cycle used to derive the compensation of surface energy on a pairwise interaction level. (a) Hydrophobic cluster of five nonbonded electrostatic beads, with Rglob and ξ being the cluster size and bead size, respectively. The dashed circle represents the schematic solvent shell, (b) separated nonbonded beads before hydrophobic clustering, each with a solvation shell, and (c) schematic solvation shell formation on the pairwise interaction level. Two close beads share one solvation shell, and there are five such shells; (d) the spherical approach of such a solvent shell in (c) with its size R2. ΔU, ΔU′, and ΔUe represent the differences of surface energies between states connected by arrows.

The limiting case where only two beads are interacting with each other, i.e., the value of ΔUe for n = 2 and C2 = 1, provides lim n → 2 ΔUe = 0

ΔUe|pair ≈

ξ2

− U0

22/3nCn 2

1 nCn 2

) = u ⎛⎜

⎞ 2 2/3 ⎟⎟ 2 − 0⎜ 1/3 ⎝ Cnn ⎠

(9)

For large clusters, we obtain ΔUe|pair ≈ −22/3u0 ,

where R2 is the size of such a cluster. In implicit-solvent CG simulations, U2 corresponds to the bead−bead pair PMF. A conformational transition of the hydrophobic PSSNa chain from the coil to the globular state actually involves a complicated process. In such a coil−globule transition, the collapsing part of the chain is desolvated while a new solvation shell is formed around the aggregate. The surface energy of the collapsed chain is smaller than the one before the transition. The difference of the surface energies is ΔU = U − U0 ≈ u0

(

u0 n2/3 −

(2)

R glob 2

(8)

In the case of hydrophobic clustering where more than two beads interact simultaneously, ΔUe is not zero and has to be considered when using pair PMFs. Ignoring this factor will lead to different predictions of polymer conformational properties compared with those from the parental atomistic simulations. The energy correction for each interacting bead pair in the globule is

bring one uncharged bead from inside the hydrophobic cluster to the surface, while Rglob2/ξ2 is proportional to the number of such CG beads on the surface of the hydrophobic cluster. The surface energy of a cluster containing only two beads is here defined as U2 and corresponds to the surface energy of two aggregating beads with an effective spherical shape (see Figure 2d):

U2 ≈ u0R 2 2/ξ 2

(7)

n≫1

(10)

The expression of ΔUe|pair in eq 9 shows a dependence on n to the power of −1/3, while the dependence of Cn on n is unknown yet and treated as a constant. Theoretically, this −1/3 dependence can be applied to systematic CG simulations of hydrophobic PSS in implicit water. Note that this theoretical derivation does not grant quantitative rescaling of the nonbonded CG potentials, since the calculation above assumes spherical clusters in which only the surface beads are solvent exposed. Nevertheless, it draws the skeleton of the rescaling function for hydrophobic PSS bead interactions and theoretically predicts the dependence of the rescaling function on the length of the polymer.

(3) C

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

3. MODEL 3.1. Mapping Schemes. The mapping scheme of ref 26 has been used. The PSS monomer, which consists of the backbone (−CH2−CH−), a phenyl ring, and a sulfonate group, is mapped on three CG beads (see Figure 3). Nonsulfonated

syndiotactic single PSSNa chains in solution with a sulfonation degree of 50% (with alternating charged and uncharged monomers) and different degrees of polymerization: 8, 16, 25, 32, 48, and 64. A total simulation time of 70 ns was used for each DP with the last 20 ns for data extraction. Averaging was performed over 2 (DP = 8, 16, 25, and 32) or three (DP = 48 and 64) independent runs. The simulation box was constructed consisting of a single PSS chain together with its counterions in water. The number of water molecules for DP = 8, 16, 25, 32, 48, and 64 were 3607, 7228, 11 275, 13 689, 21 652, and 29 329, respectively. Periodic boundary conditions were implemented in 3 dimensions, and the monomer concentration of all systems was ≈0.12 mol/L. 3.2.2. Coarse-Grained Simulations. All CG simulations have been performed in implicit water solutions under constant NVT conditions at T = 298 K. The box volumes used for the CG single-chain simulations were chosen according to the average box volumes of the corresponding atomistic NPT simulations. Counterions (Na+) were treated explicitly in the CG model. The stochastic dynamics integrator of GROMACS was employed with time coupling constant of 2 ps and an integration time step 4 fs. A simulation time of 800 ns was used with the last 50% simulation time for data sampling. Bonded and nonbonded (V short ) potentials were presented in GROMACS tabulated forms. Long-range electrostatic interactions were calculated with PME. Both the real-space contribution of PME and the short-range tabulated interactions shared the same cutoff of 1.2 nm. Fully sulfonated PSS chains (DP = 126) were used for model validation (please see Figure S8 in the Supporting Information) according to the reference experimental research on PSSNa in semidilute solutions.54 The radius of gyration (Rg) of the simulated chains was analyzed with a WLC (worm-like-chain) approach54,55 to mimic the PSSNa chains with a polydispersity around 1.15 reported in the experiments.54 In Table SI we presented the details such as concentration, simulation time, and tacticity of the PSSNa chains. 3.3. Coarse-Graining Methods. 3.3.1. Bonded Potentials. The CG bonded potentials of the PSS chain, including the bond, angle, and dihedral torsion potentials, ensure that the CG model inherits the chemical information from its parent atomistic model (such as polymer tacticity and local conformation). These potentials are derived by performing atomistic MD simulations of isolated random walks in vacuum at a selected temperature T = 298 K, where the probability distributions of the CG bond lengths, bending angles, and dihedral angles are extracted and transformed into potentials of mean force through Boltzmann inversion. Bonded potentials for isotactic and syndiotactic CG PSS chains have been previously derived by separately sampling atomistic random walks with the corresponding tacticity.26,43 Because these random walks take only local intrachain interactions into account within a finite monomer cutoff range along the backbone, the CG bonded potentials include only the effects of chain stiffness as determined by local chemical details.26 The bonded potentials were taken from ref 26. The same bonded potentials have been used for fully sulfonated and partially (50%) sulfonated chains. 3.3.2. Nonbonded Potentials. As a complementary factor to the bonded potentials for CG models, nonbonded potentials describe all the interactions beyond the reach of bonded ones in the model. In the system of PSSNa solution, the nonbonded interactions involve short-ranged interactions and long-ranged

Figure 3. Chemical structure and mapping scheme of coarse-grained beads for (a) fully sulfonated PSS and (b) 50% sulfonated PSS chains with alternating charged and uncharged monomers. (c) Schematic CG chain with black lines representing the CG bonds for 50% sulfonated PSS with alternating, sulfonated and non-sulfonated, repeat units. Coarse-grained beads A, B, and C are colored in blue, yellow, and red, respectively.

monomers are mapped on two CG beads. Following the earlier approach,26,42,43 bead A is placed at the mass center of the −CH−CH2−CH− group between two adjacent monomers along the PSS backbone, with the −CH− groups weighted with half their masses; bead B represents all the atoms of the phenyl ring and bead C the sulfonate group. The backbone of the CG chain is established with an −A−B−A−B− connecting style (see Figure 3c), while bead C is connected as a “side group” to the bead B. This mapping scheme allows to distinguish between specific tacticities (such as isotactic, syndiotactic, and atactic) of the PSS chain.43 Charged beads carry the net charge of the corresponding groups, i.e., qA = 0, qB = −0.505, and qC = −0.495. The sodium counterion is described in the same way as in former work with qNa = 1.26 3.2. Simulation Details. The molecular dynamics package GROMACS 4.644,45 was used to perform all the simulations reported in this paper. 3.2.1. Atomistic Simulations. For the atomistic modeling of NaPSS, we used a united atom model based on the TraPPE force field.46−48 We employed the partial charges assigned on PSS monomer as in previous work,26 with same the force field parameters for the counterions (Na+).49 The SPC/E water model was used.50 The LINCS method51 was used to constrain all the bond lengths of PSS, while the SETTLE algorithm was used for water molecules.52 For nonbonded short-range interactions, a cutoff distance of 1.0 nm was used for the van der Waals interactions. Long-range electrostatic interactions were calculated with the particle-mesh-Ewald (PME) method using a 1.0 nm real space cutoff. All atomistic simulations were performed under isothermal−isobaric (NPT) conditions using the stochastic dynamics integrator of GROMACS. An inverse friction constant of 2 ps and an integration time step of 1 fs were used. The Berendsen barostat was used to couple the system to an external bath at 1 bar with time coupling constant at 1 ps.53 The temperature was set at 298 K. We simulated D

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Coulomb interactions as VPMF = Vshort,ij(r) + qiqj/(4πε0εrr), where qi and qj are the charges of bead types i and j (B, C, or Na+), r is the distance between two beads, ε0 is the dielectric permittivity of vacuum, and εr is the relative dielectric permittivity of water. The CG bead charge q is taken as the sum of the partial charges on the corresponding atom groups. In the development of nonbonded potentials, VPMF is usually obtained by calculating pair PMFs between two small molecules chemically relevant to the PSS CG beads in allatom explicit-solvent simulations.26 In this way, both the internal degrees of freedom of such molecules and the solvent degrees of freedom are considered. Note that for PSSNa solution the term VPMF contains contributions from both the van der Waals and the electrostatic interactions. Vshort(r) is obtained by subtracting the Coulomb term from VPMF. Thus, Vshort(r) depends explicitly on distance r and implicitly depends on state (temperature, pressure). In the CG simulation, Vshort(r) is implemented in a tabulated form, while the remaining long-ranged Coulomb interactions are computed using the particle mesh Ewald (PME) algorithm.56 The short-range part of the potential requires a correction as explained in section 2. While in biomolecular implicit-solvent force fields atom-specific terms are used to account for implicit solvation,57−59 we here use a rescaling approach which guarantees agreement of the radius of gyration with results from fine-grained simulations. A multiplicative rescaling function f(DP) is used for nonbonded pair potentials between chosen (A−A, A−B, and B−B, while B only refers to the unsulfonated phenyl ring) pairs. The potential Vshort(r) is multiplied with f in regions where it is negative, while being multiplied with 1/f in regions where the potential is positive. This is demonstrated in Figure 4 for the A−B nonbonded interaction.

Figure 5. (a) Rg of alternatingly sulfonated PSS (CG) plotted against the rescaling factor f for different DP and (b) Values of Δf, which provide a match between the atomistic (atm) and CG Rg values, as a function of polymer chain length, DP. The dotted line is fitted to the data with a functional shape y = a/xb + c. The colored arrows target the corresponding scaling factors where atomistic and CG Rg values match. (c) Double-logarithmic plot of (Δf − c) and DP with data at DP = 3 excluded. The fitted linear curve yields a scaling exponent b at 0.381 ± 0.045, a = 3.106 ± 0.175, and c = −2.045 ± 0.337. A dashed line with the predicted b at 1/3 is also presented as a comparison. Figure 4. Comparison of the original A−B nonbonded pair potential (original AB) and rescaled ones (starting from the first attracting well) with different f values. For the AA pair we only implemented the rescaling when Vshort(r) is positive.

value of f (for DP = 8, the Rg from atomistic simulations is 0.66 ± 0.03 nm, which can be reproduced with f values between 1.2 and 2). Clearly, however, f should theoretically approach the value 1 in the limiting case of small DP where hydrophobic clusters contain at most two beads. A theoretical limit where f = 1 (no modifications made) is therefore proposed and set at DP = 3. This choice was made based on the observation that at DP = 3 the bonded interactions within the molecule prevent clustering of the hydrophobic beads, while for DP > 3 this is theoretically possible. Figure 5b shows the quantity Δf (Δf = 1 − f) versus DP. Because the rescaling factor f can be viewed as the ratio of the modified attractive pair potential (well depth of the first minimum) over its original one in Figure 4, Δf is proportional

Following this rescaling approach, CG simulations of a 50% alternatingly sulfonated PSS chain have been carried out with a series of f values, and the averaged radius of gyration Rg is plotted against f for different DP (from 8 to 64 as shown in Figure 5a). For each chain length, a nonlinear fitting is done so that Rg values obtained from a series of atomistic simulations with different chains lengths can be used to track the corresponding f value. However, this does not work properly with DP < 10, as the fitting curve of Rg versus f is almost parallel to the f-axis, and it is impossible to track the corresponding E

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules to ΔUe. Nonlinear fitting of the data yields the scaling exponent b ≈ 0.38 ± 0.1, which is close to our theoretically predicted value 1/3 (eq 9). Possibly, the deviation of the fitted scaling exponent from the theoretical value originates from the dependence of Cn on DP, which is treated as a constant in the theoretical derivation. Nevertheless, the fitted relationship basically agrees with the theoretical approach (section 2) and confirms the dependence of excess energy ΔUe on DP.

for DP > 40. Without rescaling, Rg is overestimated for every DP. The atomistic simulation data shown in Figure 6b indicate a transition in the chain conformation. While for low DP we observe Rg ∼ DP1/3, the scaling changes to Rg ∼ DP at high DP, indicating a transition from a globular to a rod-like-shaped structure. Similar transitions have been reported in the literature with atomistic simulations and basic agreement could be found.37 This transition could be identified by a combination of simulation snapshots as well as by examining the nonsphericity parameter of the chain. In Figure 7 we show the representative conformations of 50% sulfonated PSS chains with DP varying between 16 and 64. The atomistic snapshots clearly show a compact structure, but this structure slowly extends with increasing DP. At larger DP, such as DP = 48 and 64, the structure has gradually shifted from spherical to cylindrical, and the profile of averaged nonsphericity parameters at the corresponding DP values echoed this shift (Figure 8). The spherical−cylindrical transition found here, and the globule−coil transition reported in Carrillo’s paper,37 suggest that the surface energy that maintained the spherical structure has become inferior to the electrostatic repulsion with increasing DP. At this transition point, which could be roughly considered as DP = 48 in our case, Rg approximately equals 1.2 nm, which is very close to the experimental Rg value (≈1.24 nm) of densely packed beads along the PSS chain at a sulfonation degree of 64%.31 When DP < 48, the partially sulfonated PSS chain maintained a compact structure. Inside such compact globules, the distributions of hydrophobic monomers (unsulfonated phenyl rings) and sulfonate groups with respect to the center of chain mass present different profiles at DP = 32 as shown in Figure 9, with charged sulfonate group pointing outside and unsulfonated phenol rings acting as the hydrophobic clustering “central core”, in agreement with Carrillo’s findings.37 Figure 9 moreover shows that the CG simulations with “functional” rescaling reproduce the atomistic RDF profiles for both the unsulfonated phenyl rings and the sulfonate groups. Simulations of the CG model without rescaling fail to reproduce the changes in chain conformation with chain length and show extended structures for the whole DP range (Figure 7 marked with “no rescaling”). With functional rescaling, CG simulations exhibited an excellent agreement with their parental atomistic simulations by reproducing both the global conformations (Figure 7 marked with “functional rescaling”) and shape parameters (Figure 8 labeled with “functional rescaling”). It has been reported that an ambiguously designed generic model would lead to obvious structural mistakes and artificial properties compared with experimental data, due to the lack of consideration of local chemical structures and specificities.32−34 Similar conclusions may be drawn as well, even if we used systematic CG models to maintain important chemical specificities, but with improper nonbonded potentials such as “universal” rescaling in Figure 6. However, questions still remained after the agreement was achieved for atomistic simulation and rescaled systematic CG simulation. For instance, by increasing DP, there would be a global conformational transition of the hydrophobic polyelectrolyte in solution, but the atomistic simulation so far cannot produce reliable results beyond DP = 64, since system equilibration is difficult to reach as the computational time becomes unaffordable. This means that further atomistic-CG validation of the rescaling function would be unfeasible.

4. RESULTS We performed simulations of 50% sulfonated (alternatingly) syndiotactic PSSNa single-chain solutions with both atomistic and systematic CG models in order to validate the rescaling procedure applied to nonbonded interactions between uncharged beads and to study the properties of hydrophobic polyelectrolytes. Again, at first we used the global property Rg as the standard for result comparison. Then, detailed local conformations were compared, such as the average angles and torsions along the chain (see Supporting Information). In Figure 6a we present the simulation results for Rg. With the “functional” rescaling, the atomistic simulation results are reproduced by construction: the influence of chain length on the rescaling factor f is determined at every DP. If, instead, the rescaling factor determined for DP = 25 is used for all other chain lengths as well (“universal” rescaling), Rg is overestimated

Figure 6. (a) Comparison of atomistic (atm) and CG simulation results (with different rescaling methods of nonbonded potentials) for the radius of gyration (Rg) of 50% sulfonated PSSNa chains at different degrees of polymerization (DP). In “universal” rescaling, the rescaled nonbonded potentials derived from atomistic simulations at DP = 25 were used for all other chain lengths. The “no rescaling” results were obtained with the original nonbonded potentials. The “functional” rescaling results were obtained using a rescaling factor f determined at each chain length from the fitted curve in Figure 5b. (b) Doublelogarithmic plot of atomistic simulation results. The solid lines (with slopes of 1/3 and 1, respectively) serve to guide the eye. F

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 7. Snapshots showing global conformation of 50% sulfonated PSSNa. “atm” stands for atomistic simulations; “no rescaling” and “functional rescaling” refer to CG simulations without and with “functional” rescaling of the pair potentials, respectively.

enlightened by Lord Rayleigh’s charged droplet theory,60 a commonly accepted explanation brought forward by Dobrynin et al. through theoretical derivation combined with simulations.61 Dobrynin et al. proposed a theoretical solution of the charged chain conformation in poor solvent and reached to the conclusion that with increasing chain length a “pearl-necklace” cascade increasing pattern may occur. The pearl-necklace chain is constructed by a growing number of densely packed spheres separated by stretched chain blocks, and the argument has been made that this conformation has a lower free energy than the cylindrical globule.61 This special conformation can hardly be observed in atomistic simulations due to long time and large length scale requirements37,62 but has been justified by generic model simulations, theoretical approaches,63−66 and experimental studies of partially sulfonated PSSNa solutions,31 which act as a solid reference for systematic CG simulations for PSSNa solutions. Therefore, a possible validation of our rescaling algorithm at increased DP is to check if it follows the principles of such a cascade evolving pattern. A series of larger DPs ranging from 80 to 160 were considered with the corresponding functional rescaling factor (read from Figure 5b) applied to the 50% sulfonated CG PSS

Figure 8. Nonsphericity parameter δ of 50% a sulfonated PSSNa chain in water, where δ = [(L12 − L22)2 + (L22 − L32)2 + (L32 − L12)2]/ [2(L12 + L22 + L32)2], and L1, L2, and L3 are eigenvalues of the Rg tensor. The nonsphericity parameter varies from δ = 0 to δ = 1; 0 indicates a perfect spherical shape and 1 a rodlike structure.

Another question corresponds to the situation after the conformational transition with larger DP. Is there any stable global structure that the hydrophobic polyelectrolyte in solution would have? The answer to the second question is

Figure 9. Comparison of radial distribution function (RDF) of (a) unsulfonated phenyl ring and (b) sulfonate group around the center of mass of the PSS chain at DP = 32. The atomistic coordinates were mapped on the CG beads for the RDF calculation. G

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules chain. The Rg of PSS with DP = 48−160 is plotted in Figure 10. As guided by the fitted linear curve, Rg of PSS chains at DP ≥

increased monomer numbers in both the compact domains and the string connecting them. This conformational transition from globule to dumbbell indicates that the electrostatic repulsion exceeds the surface energy of the globule. This results in the formation of two compact domains in which the balance between surface energy and electrostatic repulsion is restored. A third domain appears when DP = 128, and this structure remains unchanged after reaching DP = 144. Afterward, a fourth domain appears when DP = 160. Fluctuations of these structures were observed in the simulations, corresponding mainly to transitions between n condensed beads with n + 1 or n − 1 beads along the chain. This happened when the pearlnecklace conformation evolved with the increase of chain length. Yet, with rescaling, fully extended or globular conformations for those larger DPs were never observed. The appearance of these domains is in agreement with the theoretical description, i.e., eq 2.12 in Dobrynin’s paper.61 The above observation of domain formation along the chain corresponds to a pearl-necklace cascade evolving pattern. Thus, we achieved the consistency with theoretical prediction on both global conformation and scaling behavior of PSS chains in water solution, while the important local chemical information on the polymer was integrated in the systematic CG model. A parallel experimental report could not be found since the observations are mostly made in semidilute solutions and the DP of related PSS is usually larger and of at the order of 103.

Figure 10. Rg of 50% sulfonated PSSNa chains in water plotted against DP. Data from CG simulations with “functional” rescaling of the nonbonded potentials.

48 is proportional to DP. Although this linear dependence on DP has been shared by other representative conformations that hydrophobic polyelectrolytes might take in solution, such as the cylindrical61 or simply extended,39 the following observation of PSSNa in solution (Figure 11) clearly demonstrates that from

5. SUMMARY AND CONCLUSIONS Hydrophobic clustering effects play an important role in the self-assembly and aggregation behavior of biological and synthetic systems. Examples include physical gelation and the formation of lipid bilayers, capsules, and micelles. Implicitsolvent coarse-grained simulations of specific systems can be performed with systematically coarse-grained models and rely on the inclusion of potential energy functions, which properly balance the hydrophilic (electrostatic) and hydrophobic interactions operative in these systems. The omission of explicit solvent molecules, however, complicates the coarsegrained simulations because in aggregating systems hydrophobic interactions contain multibody contributions. This work provides a simple solution to overcome this issue for a hydrophobic (partially collapsed) polyelectrolyte composed of sodium poly(styrenesulfonate) with a sulfonation degree of 50%. The idea described in this paper is to rescale pair potentials of mean force between uncharged beads on the polymer (derived from atomistic simulations) with a chainlength-dependent scaling factor. This factor effectively corrects for changes in solvent-accessible molecular surface area if the hydrophobic cluster size increases from two nonpolar beads surrounded by water to a larger globule in which only the surface beads are solvent-exposed. Although an appropriately rescaled pairwise interaction potential does not capture the full multibody nature of hydrophobic clustering in aqueous solution, and needs to be tested for other systems, we herein demonstrate that it is accurate enough to reproduce the properties of an atomistic system composed of partially sulfonated polystyrene chains. We note that the validity of a pairwise decomposition of the n-body solvation potential has previously been successfully tested for flexible chains in a model solvent, showing that the effects of solvation on chain conformations can be captured with effective two-body interactions.67,68 We further show that the rescaling approach proposed in the present work provides a way to realistically

Figure 11. Snapshots of PSS chains (50% sulfonation degree) at DP = 64, 80, 96, 128, 144, and 160 showing the cascade evolving of pearlnecklace conformation.

DP = 64 to 80 a transition from a cylindrical globule to a dumbbell conformation occurs. This transition is difficult to observe in all-atom simulations.37 Mantha et al. reported an observation of pearl-necklace conformation of partially sulfonated PSS of DP = 48 by varying the sulfonation and charge distribution.62 The authors used a bottom-up CG method with explicit CG water; therefore, the dependence of hydrophobic clustering on DP could theoretically be ignored in the development of CG potentials. However, the dependence of the chain conformation on DP was not investigated.62 Next, by further increasing DP from 80 to 96, the chain conformation (Figure 11) maintains the dumbbell shape with slightly H

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

(4) Ennari, J.; Elomaa, M.; Sundholm, F. Modeling a polyelectrolyte system in water to estimate the ion-conductivity. Polymer 1999, 40, 5035−5041. (5) Lyubartsev, A. P.; Laaksonen, A. Molecular dynamics simulations of DNA in solution with different counter-ions. J. Biomol. Struct. Dyn. 1998, 16, 579−592. (6) Duca, J.; Hopfinger, A. J. Molecular modeling of polymers. Comput. Theor. Polym. Sci. 1999, 9, 227−244. (7) Holm, C.; Joanny, J. F.; Kremer, K.; Netz, R. R.; Reineker, P.; Siedel, C.; Vilgis, T. A.; Winkler, R. G. In Polyelectrolytes with Defined Molecular Architecture II; Springer: Heidelberg, 2004; Vol. 166, pp 3− 17. (8) Hoda, N.; Larson, R. G. Explicit- and implicit-solvent molecular dynamics simulations of complex formation between polycations and polyanions. Macromolecules 2009, 42, 8851−8863. (9) Qiao, B.; Cerda, J. J.; Holm, C. Atomistic study of surface effects on polyelectrolyte adsorption: Case study of a poly(styrenesulfonate) monolayer. Macromolecules 2011, 44, 1707−1718. (10) Algaer, E. A.; van der Vegt, N. F. A. Hofmeister ion interactions with model amide compounds. J. Phys. Chem. B 2011, 115, 13781− 13787. (11) Park, S.; Zhu, X.; Yethiraj, A. Atomistic simulations of dilute polyelectrolyte solutions. J. Phys. Chem. B 2012, 116, 4319−4327. (12) Delle Site, L.; Holm, C.; van der Vegt, N. F. A. Multiscale approaches and perspectives to modeling aqueous electrolytes and polyelectrolytes. Top. Curr. Chem. 2011, 307, 251−294. (13) Rodriguez-Ropero, F.; van der Vegt, N. F. A. Direct osmolyte− macromolecule interactions confer entropic stability to folded states. J. Phys. Chem. B 2014, 118, 7327−7334. (14) Rodriguez-Ropero, F.; Hajari, T.; van der Vegt, N. F. A. Mechanism of polymer collapse in miscible good solvents. J. Phys. Chem. B 2015, 119, 15780−15788. (15) Stevens, M. J.; Kremer, K. The nature of flexible linear polyelectrolytes in salt free solution: A molecular dynamics study. J. Chem. Phys. 1995, 103, 1669−1690. (16) Limbach, H. J.; Holm, C. End effects of strongly charged polyelectrolytes: A molecular dynamics study. J. Chem. Phys. 2001, 114, 9674−9682. (17) Limbach, H. J.; Holm, C. Single-chain properties of polyelectrolytes in poor solvent. J. Phys. Chem. B 2003, 107, 8041− 8055. (18) Lyubartsev, A. P.; Laaksonen, A. Calculation of effective interaction potentials from radial distribution functions: A reverse monte carlo approach. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 52, 3730−3737. (19) Tschöp, W.; Kremer, K.; Batoulis, J.; Bürger, T.; Hahn, O. Simulation of polymer melts. I. Coarse-graining procedure for polycarbonates. Acta Polym. 1998, 49, 61−74. (20) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving effective mesoscale potentials from atomistic simulations. J. Comput. Chem. 2003, 24, 1624−1636. (21) Izvekov, S.; Voth, G. A. A multiscale coarse-graining method for biomolecular systems. J. Phys. Chem. B 2005, 109, 2469−2473. (22) Noid, W. G.; Chu, J. W.; Ayton, G. S.; Krishna, V.; Izvekov, S.; Voth, G. A.; Das, A.; Andersen, H. C. The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models. J. Chem. Phys. 2008, 128, 244114. (23) Shell, M. S. The relative entropy is fundamental to multiscale and inverse thermodynamic problems. J. Chem. Phys. 2008, 129, 144108. (24) Brini, E.; Marcon, V.; van der Vegt, N. F. A. Conditional reversible work method for molecular coarse graining applications. Phys. Chem. Chem. Phys. 2011, 13, 10468−10474. (25) Ganguly, P.; Mukherji, D.; Junghans, C.; van der Vegt, N. F. A. Kirkwood−buff coarse-grained force fields for aqueous solutions. J. Chem. Theory Comput. 2012, 8, 1802−1807. (26) Li, C.; Shen, J.; Peter, C.; van der Vegt, N. F. A. A chemically accurate implicit-solvent coarse-grained model for polystyrenesulfonate solutions. Macromolecules 2012, 45 (5), 2551−2561.

describe the conformations of partially collapsed polyelectrolyte chains. The dependence of the rescaling factor on the chain length has been theoretically derived and confirmed by results of coarse-grained and extensive, detailed-atomistic, simulations. Extrapolation of the rescaling factor of short chains to large chain lengths, using the theoretical relationship derived in this work, is shown to yield models that verify classic hydrophobic polyelectrolyte behavior in aqueous solution, i.e., a pearlnecklace conformation with a cascade grow pattern. It should be noted that while the same functional form of the rescaling factor can in principle be applied to different polymer backbones (with, e.g., different sulfonation degrees), a reparametrization of the rescaling function is required in those cases. With the modified systematic coarse-graining method described in this work, we extend the transferability of fragment-based coarse-grained models and expect that other types of hydrophobic aggregation behavior may be studied as well.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.6b01132. Graphs of the bonded and nonbonded CG potentials developed from atomistic simulations, including the rescaled nonbonded potentials for different DP; comparison of detailed local conformations (graphs of probability distributions of CG bonds, angles, and torsions) between atomistic and CG models; validation of the fully sulfonated PSSNa CG model against experimental results by means of matching the radius of gyration (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected]; Tel +86-431-85262516; Fax +86-431-85262969 (R.Z.). Present Address

R.Z.: State Key Laboratory of Polymer Physics and Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun 130022, P. R. China. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We are grateful to Samira Hezaveh for setting up the coarsegrained simulations discussed in the Supporting Information. This project has been carried out as part of the Collaborative Research Center Transregio TRR146 “Multiscale Simulation Methods for Soft Matter Systems” funded by the German Science Foundation.



REFERENCES

(1) Zhang, Y.; Cremer, P. S. Interactions between macromolecules and ions: The hofmeister series. Curr. Opin. Chem. Biol. 2006, 10, 658−663. (2) Ball, P. Water as an active constituent in cell biology. Chem. Rev. 2008, 108, 74−108. (3) Zhang, Y.; Cremer, P. S. Chemistry of hofmeister anions and osmolytes. Annu. Rev. Phys. Chem. 2010, 61, 63−83. I

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules (27) Villa, A.; Peter, C.; van der Vegt, N. F. A. Self-assembling dipeptides: Conformational sampling in solvent-free coarse-grained simulation. Phys. Chem. Chem. Phys. 2009, 11, 2077−2086. (28) Hess, B.; Holm, C.; van der Vegt, N. F. A. Modeling multi-body effects in ionic solutions with a concentration dependent dielectric permittivity. Phys. Rev. Lett. 2006, 96, 147801. (29) Shen, J.; Li, C.; van der Vegt, N. F. A.; Peter, C. Transferability of coarse grained potentials: Implicit solvent models for hydrated ions. J. Chem. Theory Comput. 2011, 7, 1916−1927. (30) Qu, D.; Baigl, D.; Williams, C. E.; Möhwald, H.; Fery, A. Dependence of structural forces in polyelectrolyte solutions on charge density: A combined afm/saxs study. Macromolecules 2003, 36, 6878− 6883. (31) Spiteri, M. N.; Williams, C. E.; Boue, F. Pearl-necklace-like chain conformation of hydrophobic polyelectrolyte: A sans study of partially sulfonated polystyrene in water. Macromolecules 2007, 40, 6679−6691. (32) Chang, R.; Yethiraj, A. Dilute solutions of strongly charged flexible polyelectrolytes in poor solvents: Molecular dynamics simulations with explicit solvent. Macromolecules 2006, 39, 821−828. (33) Park, S.; Zhu, X.; Yethiraj, A. Atomistic simulations of dilute polyelectrolyte solutions. J. Phys. Chem. B 2012, 116, 4319−4327. (34) Reddy, G.; Yethiraj, A. Implicit and explicit solvent models for the simulation of dilute polymer solutions. Macromolecules 2006, 39, 8536−8542. (35) Foley, T. T.; Shell, M. S.; Noid, W. G. The impact of resolution upon entropy and information in coarse-grained models. J. Chem. Phys. 2015, 143 (24), 243104. (36) Villa, A.; Peter, C.; van der Vegt, N. F. A. Transferability of nonbonded interaction potentials for coarse-grained simulations: Benzene in water. J. Chem. Theory Comput. 2010, 6, 2434−2444. (37) Carrillo, J.-M. Y.; Dobrynin, A. V. Detailed molecular dynamics simulations of a model napss in water. J. Phys. Chem. B 2010, 114, 9391−9399. (38) Chandler, D. Interfaces and the driving force of hydrophobic assembly. Nature 2005, 437, 640−647. (39) Dobrynin, A. V.; Rubinstein, M. Theory of polyelectrolytes in solutions and at surfaces. Prog. Polym. Sci. 2005, 30, 1049−1118. (40) de Gennes, P. G.; Pincus, P.; Brochard, F.; Velasco, R. M. Remarks on polyelectrolyte conformation. J. Phys. (Paris) 1976, 37 (12), 1461−1473. (41) Dobrynin, A. V.; Colby, R. H.; Rubinstein, M. Scaling theory of polyelectrolyte solutions. Macromolecules 1995, 28, 1859−1871. (42) Harmandaris, V. A.; Reith, D.; van der Vegt, N. F. A.; Kremer, K. Comparison between coarse-graining models for polymer systems: Two mapping schemes for polystyrene. Macromol. Chem. Phys. 2007, 208 (2), 2109−2120. (43) Fritz, D.; Harmandaris, V. A.; Kremer, K.; van der Vegt, N. F. A. Coarse-grained polymer melts based on isolated atomistic chains: Simulation of polystyrene of different tacticities. Macromolecules 2009, 42, 7579−7588. (44) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Gromacs: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (45) Lindahl, E.; Hess, B.; van der Spoel, D. Gromacs 3.0: A package for molecular simulation and trajectory analysis. J. Mol. Model. 2001, 7, 306−317. (46) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable potentials for phase equilibria. 6. United-atom description for ethers, glycols, ketones, and aldehydes. J. Phys. Chem. B 2004, 108, 17596−17605. (47) Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (48) Vishnyakov, A.; Neimark, A. V. Specifics of solvation of sulfonated polyelectrolytes in water, dimethylmethylphosphonate, and their mixture: A molecular simulation study. J. Chem. Phys. 2008, 128, 164902. (49) Hess, B.; van der Vegt, N. F. A. Cation specific binding with protein surface charges. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 13296−13300.

(50) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269−6271. (51) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. Lincs: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (52) Miyamoto, S.; Kollman, P. A. Settle: An analytical version of the shake and rattle algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952−962. (53) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (54) Nierlich, M.; Boué, F.; Lapp, A.; Oberthur, R. Radius of gyration of a polyion in salt free polyelectrolyte solutions measured by s. A. N. S. J. Phys. (Paris) 1985, 46, 649−655. (55) Oberthuer, R. C. Radius of gyration versus molar mass relations from light scattering for polydisperse non-gaussian chain molecules. Makromol. Chem. 1978, 179, 2693−2706. (56) Darden, T.; York, D.; Pedersen, L. Particle mesh ewald: An n· log(n) method for ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089. (57) Kleinjung, J.; Fraternali, F. Design and application of implicit solvent models in biomolecular simulations. Curr. Opin. Struct. Biol. 2014, 25, 126−134. (58) Kleinjung, J.; Scott, W. R. P.; Allison, J. R.; van Gunsteren, W. F.; Fraternali, F. Implicit solvation parameters derived from explicit water forces in large-scale molecular dynamics simulations. J. Chem. Theory Comput. 2012, 8, 2391−2403. (59) De Simone, A.; Dodson, G. G.; Verma, C. S.; Zagari, A.; Fraternali, F. Prion and water: Tight and dynamical hydration sites have a key role in structural stability. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 7535−7540. (60) Rayleigh, L. On the equilibrium of liquid conducting masses charged with electricity. Philos. Mag. 1882, 14, 184−186. (61) Dobrynin, A. V.; Rubinstein, M.; Obukhov, S. P. Cascade of transitions of polyelectrolytes in poor solvents. Macromolecules 1996, 29, 2974−2979. (62) Mantha, S.; Yethiraj, A. Conformational properties of sodium polystyrenesulfonate in water: Insights from a coarse-grained model with explicit solvent. J. Phys. Chem. B 2015, 119, 11010−11018. (63) Chodanowski, P.; Stoll, S. Monte carlo simulations of hydrophobic polyelectrolytes: Evidence of complex configurational transitions. J. Chem. Phys. 1999, 111, 6069−6081. (64) Lyulin, A. V.; Dünweg, B. D.; Borisov, O. V.; Darinskii, A. A. Computer simulation studies of a single polyelectrolyte chain in poor solvent. Macromolecules 1999, 32, 3264−3278. (65) Migliorini, G.; Lee, N.; Rostiashvili, V.; Vilgis, T. A. Polyelectrolyte chains in poor solvent. A variational description of necklace formation. Eur. Phys. J. E: Soft Matter Biol. Phys. 2001, 6, 259−270. (66) Solis, F. J.; Olvera de la Cruz, M. Variational approach to necklace formation in polyelectrolytes. Macromolecules 1998, 31, 5502−5506. (67) Taylor, M. P.; Petersen, G. M. Solvation potentials for flexible chain molecules in solution: On the validity of a pairwise decomposition. J. Chem. Phys. 2007, 127, 184901. (68) Taylor, M. P.; Ye, Y.; Adhikari, S. R. Conformation of a flexible polymer in explicit solvent: Accurate solvation potentials for lennardjones chains. J. Chem. Phys. 2015, 143, 204901.

J

DOI: 10.1021/acs.macromol.6b01132 Macromolecules XXXX, XXX, XXX−XXX