Association Equilibrium for Cross-Associating Chains in a Good

The difference between terminal and nonterminal sites becomes less .... vs the fraction of p-chains for constant content of solvent at 365.4 K. Symbol...
2 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Association Equilibrium for Cross-Associating Chains in a Good Solvent: Crowding and Other Nonideality Effects Igor Yu. Gotlib,* Ivan K. Malov, Alexey I. Victorov, and Mikhail A. Voznesenskiy St. Petersburg State University, 7/9 Universitetskaya nab., St. Petersburg 199034, Russia S Supporting Information *

ABSTRACT: Association equilibrium has been studied by molecular dynamics (MD) for mixtures of cross-associating molecules (n-decamer+p-dimer and n-decamer+p-decamer) in a good solvent. Each monomer of n-decamers carries an associative site (n-sticker); each molecule of the second component contains two terminal associative sites (p-stickers). To model the univalent association between the n-sticker and the p-sticker, a technique based on introduction of dummy atoms has been used. We report MD data on the effects of temperature, chain flexibility, and location of the sticker within the chain on the association equilibrium. We find that the presence of nonassociating monomer units of p-chain has a substantial effect on the association equilibrium. This effect is similar to “crowding” in reactive mixtures known to be caused by the presence of inert molecules. Widely used mean field theories of associating chains (e.g., SAFT or Semenov− Rubinstein theory) do not account for the effect of crowding caused by the inert fragments of reactive chains. We introduce simple empirical corrections for crowding that describe association equilibrium in the presence of nonassociating fragment in a chain-like molecule.

1. INTRODUCTION Mixtures of reversibly associating polymer chains with tunable association strength may form responsive gels and find numerous applications in engineering artificial tissues; controlling release of drugs; and as self-healing coatings, adhesives, etc.1−3 Extensive computer simulations of systems containing associative (sticky) chains have been devoted to viscoelastic and transport properties, kinetics of clustering and network formation, sol−gel transitions, swelling of gels, and release of chemicals, see e.g., refs 4 and 5 and references therein. For a solution of reversibly associating chains, a key issue is the equilibrium balance between the bonded and free associative groups (stickers) because this balance controls formation and properties of a transient network. For describing association equilibrium between specifically interacting groups, the statistical associating fluid theory (SAFT) has been used most widely.6−12 Other analytical approximations for association equilibrium between the sticky chains have also been developed for self-associating chains,13,14 cross-associating chains,15,16 as well as for systems where self-association may compete with cross-association.17−20 Most of the simulation studies of associative chains are devoted to telechelic polymers with multiple participation of stickers in the cross-links.4,21−24 Previous work on chains with uni-univalent association includes Monte Carlo (MC) simulations where association is treated as formation of a new (“united”) “atom”,25 MC simulation of chain association to stickers placed regularly on a lattice,26 and MC simulation of telechelic chains with univalent terminal ligands forming complexes with trivalent metal ions.27 A hybrid molecular dynamics/MC technique has been applied © 2016 American Chemical Society

to study thermodynamic and kinetic properties of reversibly associating polymers that contain univalent sticky monomer units.23,28 MC simulations have also been performed for chains formed by reversibly associating spherical particles with sticky sites on the surface (“patchy colloids”).29,30 To the best of our knowledge, there has been no simulation studies of association equilibrium between dissimilar chains with a realistic interaction potential for the stickers that may form a univalent specific bond. Important examples of such systems are mixtures of polymers containing either proton donor or proton acceptor groups,31,32 polyanion−polycation gels,33 and mixtures of a polymer containing ligand groups and a complex-forming organometallic compound.34 In this work we perform molecular dynamic (MD) simulation for mixtures of cross-associating molecules, n-decamer+p-dimer and n-decamer+p-decamer, in a good solvent. Each monomer of n-decamers carries an associative site (n-sticker); each molecule of the second component contains two terminal associative sites (p-stickers). We examine effects causing nonideality of the association equilibrium, particularly the effect of crowding in the system containing n-decamer and p-decamer. This effect, i.e., the shift of association equilibrium caused by adding an inert component (the crowder), has been widely discussed for biomolecular systems.35−40 More generally, crowding is defined as “the nonspecific influence of steric repulsions on specific Received: December 22, 2015 Revised: May 27, 2016 Published: June 30, 2016 7234

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B

backbone atoms and solvent is zero while short-range repulsion acts between all nonassociating dummy atoms. This model allows to simulate uni-univalent cross-association (that may occur, e.g., due to hydrogen bonding or donor−acceptor interactions) by means of standard short-ranged model potentials that are simple to use in atomistic MD simulations. For all possible combinations of sites (atoms), the pair potentials for nonbonded interactions (describing dispersion forces and reversible specific association, not permanent chemical bonds or valence angles) are presented in Table 1.

reactions and processes that occur in highly volume-occupied media”.40 Details of simulation are explained in Section 2. The results are presented in Section 3. We start with mixtures containing n-decamer and p-dimer (where associating molecules do not carry inert fragments) and discuss the effect of different factors on association equilibrium: chain flexibility, location of stickers, and preferential solvation. Next we consider the n-decamer−p-decamer system where each p-decamer carries eight inert monomers, and we find that the presence of the inert parts of p-chains makes an essential additional contribution to nonideality of the association equilibrium. We discuss limitations of existing association theories in describing the crowding effect of this type and propose empirical corrections that help to take this effect into account.

Table 1. Nonbonded Interaction Potentials Uα−β for Different Combinations of Atom Types α and β in the Model Systems

2. MODEL AND SIMULATION PROCEDURE We model Lennard-Jones chains (decamers) or dumbbells (dimers) decorated with dummy atoms that serve to model strong short-range attraction between specific sites resulting in uni-univalent association. The solvent consists of single LennardJones atoms (AS). Other molecules include three different types of monomers: n-decamer consists of ten n-monomers, p-dimer consists of two p-monomers, p-decamer consists of two p-monomers at the ends separated by eight o-monomers. Association is only possible between n-monomers and p-monomers while interactions between nonassociating units (n−n, p−p, o−o, n−o, p−o) are (approximately) energetically equivalent to solvent−solvent and solvent−monomer interactions. Each monomer unit is a combination of two sites: the backbone atom (AN, AP, AO) and the dummy atom (DN, DP, DO)see Figure 1. The interaction energy of a dummy atom with the

αβ

AS

AN

AP

AO

DN

DP

DO

AS AN AP AO DN DP DO

LJ′a LJ′ LJ′ LJ′ 0 0 0

LJ′ LJb LJ LJ 0 0 0

LJ′ LJ LJ LJ 0 0 0

LJ′ LJ LJ LJ 0 0 0

0 0 0 0 LJRc LJA LJR

0 0 0 0 LJAd LJR LJR

0 0 0 0 LJR LJR LJR

a

LJ′: Lennard-Jones potential for interactions of solvent molecules with other solvent molecules or with monomers in chains, ⎛⎛ σ ′ ⎞12 ⎛ σ ′ ⎞6 ⎞ Uα − β(r ) = 4ε′⎜⎜⎜ ⎟ − ⎜ ⎟ ⎟⎟, NAε′ = 3.0735 kJ/mol, ⎝r ⎠⎠ ⎝⎝ r ⎠ σ ′ = 0.436759 nm

b LJ: Lennard-Jones potential for nonspecific interactions between monomers in chains,

⎛⎛ σ ⎞12 ⎛ σ ⎞6 ⎞ Uα − β(r ) = 4ε⎜⎜ ⎟ − ⎜ ⎟ ⎟, NAε = 3.375 kJ/mol, ⎝r⎠ ⎠ ⎝⎝ r ⎠ σ = 0.43 nm c

LJR: Lennard-Jones repulsion between nonassociating dummy atoms, ⎛σ ⎞12 Uα − β(r ) = 4εNPR ⎜ NPR ⎟ , NAεNPR = 24.894 kJ/mol, ⎝ r ⎠ σNPR = 0.3 nm d

LJA: associative interaction, ⎛⎛ σ ⎞12 ⎛ σ ⎞6 ⎞ Uα − β(r ) = 4εNP⎜⎜ NP ⎟ − ⎜ NP ⎟ ⎟ , ⎝ r ⎠⎠ ⎝⎝ r ⎠ NA εNP = 24.894 kJ/mol, σNP = 0.1 nm

where NA is the Avogadro number, r is the distance between interacting atoms.

Dummy atoms DN and DP act as stickers for associating chains; the parameter εNP controls the strength of association of an n−p pair. A system without association, where UDN−DP(r) is identical to the interaction potential between nonassociating dummy atoms, has also been simulated as a reference. The solvent is nearly athermal: at large distances, a monomer unit in the chain or dumbbell is “seen” as a LennardJones interaction center that is approximately equivalent to a solvent molecule. This is achieved by choosing the parameters ε,σ (interactions between backbone atoms), ε′,σ′ (solvent− backbone and solvent−solvent interactions), and εNPR,σNPR

Figure 1. A schematic view of n- and p-molecules and the vector equations for the location of a dummy atom ((DN)i, (DP)i, (DO)i) relative to the backbone atom ((AN)i, (AP)i, (AO)i); i is the number of the monomer unit in the chain, from 1 to 10. 7235

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B

calculated for the system with the same percentage composition and N = 50 000 show too strong fluctuations making it difficult to assess how the average values of those quantities change with the system composition. Standard periodic boundary conditions are used. The equilibration time varies from 2 to 6 ns and the length of a productive run is between 1 and 4.5 ns, depending on the system composition (longer simulation runs have been performed for dilute solutions; getting stable fractions of associated/ nonassociated stickers can take a long time, so the equilibration time has been higher than the productive run time in some cases). The numbers of solvent molecules, n-, p-, and o-monomers in the cell are denoted by Nsolv, Nn, Np, and No, respectively. For n-decamer−p-dimer solutions, No = 0; for n-decamer−p-decamer solutions, No = 4Np. The number of n-molecules in the cell is Nn/10 and the number of p-molecules is Np/2. The mole fractions of the monomer units are given by xn = Nn/N and xp = (Np + N0)/N. The mole fraction of solvent (in terms of nondummy atoms) is given by xsolv = 1 − xn − xp. The compositions of simulated mixtures fall on lines shown in Figure 2. The composition lines 1,

(repulsion between nonassociating dummy atoms) so that 6 6 4ε′σ′12 = 4εσ12 + 4εNPRσ12 NPR and 4ε′σ′ = 4εσ . The athermal solvent is a typical example of a good solvent where mixing is thermodynamically favored.41 The pair potentials are modified by the “force switch” cutoff procedure when the force is multiplied by a “switch function” W(r), ⎧1 if r < rsw ⎪ W (r ) = ⎨ S(r ) if rsw ≤ r < rc , ⎪ ⎪0 if r ≥ rc ⎩

(1)

where S(r) is defined in such a way that the force, and hence the potential, are continuous and continuously differentiable;42,43 in the present work, S(r) is constructed as a cubic polynomial. For DN−DP interactions, rsw = 0.16 nm, rc = 0.2 nm (the attraction that causes association is short-ranged, the energy minimum is at 0.112 nm, and the width of the potential’s well is somewhere between 0.16 and 0.2 nm). For all other interactions, rsw = 1.6 nm and rc = 2 nm. For chemical bonds in chains and dumbbells, a harmonic 1 potential is used, Ubond = 2 k(r − r0)2 , where r is the distance between bonded atoms, r0 = 0.33 nm, k = 17 000 kJ mol−1 nm−2. Valence angle potential in chains is calculated as 1 Uangle = 2 kθ(cosθ − cosθ0)2 , where θ is the angle between adjacent bonds, θ0 = 130°, kθ = 85 kJ mol−1. These bond and angle parameters, as well as the values of ε and σ for the backbone atoms, are taken from a coarse-grained model of poly(ethylene glycol)44 that can be regarded as a simple example of a polymer with associative sites (oxygen atoms with lone electron pairs). This potential describes a flexible chain, with a wide angle distribution at temperatures of a few hundreds of Kelvin. The angle distribution may vary depending on the molecular environment, affecting the chains’ associative behavior. When the valence angle becomes close to 180°, the sticker moves toward the backbone atom and repulsion between AN and AP sites overrides DN−DP attraction: the monomer in the vertex of the angle loses its association capability. To check for a specific contribution from this effect, the system with more rigid chains (kθ = 850 kJ mol−1) has also been simulated. The association criterion is the proximity of stickers: two monomer units are considered as “bonded” when the DN site is closer than 0.16 nm to the DP site. The whole set of energetic and geometric characteristics of association in the model can be considered a rough description of hydrogen bonding or a donor−acceptor complex-forming interaction between two molecular speciesa chain consisting of several associative units and a dumbbell or a telechelic chain acting as a cross-linker. The MD simulation is performed using the GROMACS package43,45 on computers in Saint Petersburg State University computational resource center and on “Lomonosov” supercomputer in the Supercomputing Center of Lomonosov Moscow State University.46 A mixture containing given amounts of the solvent, n-component and p-component, is equilibrated using the Nosé−Hoover thermostat (T = 365.4 K = 0.9ε/kB and T = 406 K = ε/kB) and Martyna−Tuckerman−Tobias−Klein barostat (p = 10 bar). The total number N of nondummy atoms (solvent particles and backbone atoms in monomer units) in the MD cell is 50 000 except for solutions with low concentration of n- or p-component where N = 400 000. The model system with N = 400 000 has been simulated when the quantities

Figure 2. Simulated mixtures’ compositions.

2, 3, and 7 correspond to fixed values of xn + xp (0.2, 0.5, 0.8, and 1, respectively). Upon removing all nonterminal monomers (o-monomers) from a p-decamer molecule, it transforms into a p-dimer. Replacing all the o-monomers by the solvent molecules while keeping Nn and Np unchanged in n-decamer+p-decamer +solvent mixtures along the lines 1, 2, and 3, we obtain composition lines 4, 5, and 6 (xn + 5xp = 0.2, xn + 5xp = 0.5 and xn + 5xp = 0.8, respectively) for the n-decamer−p-dimer solution.

3. RESULTS AND DISCUSSION 3.1. The n-Decamer−p-Dimer System. Figure 3 shows the “apparent” constants of association equilibrium between the stickers K ass = =

7236

(Npairs/V ) ((Nn − Npairs)/V )((Np − Npairs)/V ) pn pp (1 − pn )(1 − pp )(Npairs/V )

(2)

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B

Figure 4. “Apparent” association equilibrium constants for different types of n-monomers (associating sites in n-chains) in the n-decamer− p-dimer systems at 365.4 K (composition line 3) from MD simulation. The error bars show the mean square fluctuation.

Figure 3. “Apparent” association equilibrium constant vs total density of n- and p-stickers for the n-decamer−p-dimer system at two temperatures. Results are shown for systems containing flexible (kθ = 85 kJ mol−1) and rigid (kθ = 850 kJ mol−1) decamers. The error bars show the mean square fluctuation.

for the chain ends than for monomer units in the middle of a chain, and for monomers in the middle, it decreases with increasing chain rigidity much more significantly than for terminal monomers. Figure 5 shows the system volume per nondummy atom (“per monomer”), the total configuration energy per monomer, and

for the model n-decamer−p-dimer system. Here, Npairs is the number of associated n−p pairs of stickers in the MD cell, V is the cell volume, and pn and pp are the fractions of n- and p-stickers, respectively, that are part of the associated pairs. The total density of n- and p-stickers is given by (Nn + Np)/V = (xn + xp)/ v*, where v* = V/N is the volume per nondummy atom. Figure 3 demonstrates a systematic deviation from the equilibrium in an ideal mixture of reacting sites where eq 2 gives the true equilibrium constant that does not depend on the mixture composition. Here, on the contrary, the “apparent” equilibrium constant depends on composition, clearly showing a correlation with the total density of stickers, (Nn + Np)/V. Except for dilute solutions, Kassoc tends to increase with (Nn + Np)/V, especially in concentrated mixtures (see the upturns for lines 3 and 7 in Figure 3) where associated n−p pairs are present in substantial amounts. There is no qualitative change in these graphs when flexible n-chains are replaced by rigid ones but the absolute value of Kassoc decreases 1.4−2 times. Strictly speaking, the associating sites in n-chains are not equivalent but can be divided into five types, according to their location in the chain, from 1 (terminal sites) to 5 (the sites farthest from the chain end). In the model system with flexible n-molecules, Kass is much lower for terminal sites than for those in the middle of the chain (Figure 4; Supporting Information, Figure SI1). For the nonterminal positions, the differences in Kass are small; only the sites of type 2 are somewhat less associative than the other middle-chain sites in mixtures with high concentration of n-component. In the system with rigid chains, Kass for the nonterminal sites becomes much lower and is indeed close to Kass for the terminal sites. These results can be explained by the fact that associative bonding is accompanied by entropic loss due to lost conformational freedom. A nonassociated terminal site has more conformational freedom than a nonterminal site, and a nonterminal site has more conformational freedom in a flexible chain than in a rigid one (while for chain ends, the restricting effect of chain rigidity is not significant). Hence, the entropic loss caused by bonding is larger for terminal sites than for nonterminal sites, and for nonterminal sites, it is larger in flexible chains than in rigid ones. Therefore, Kass is lower

Figure 5. MD simulation results for the n-decamer−p-dimer systems at 365.4 K (composition line 3). Left: The volume per monomer, v* = V/N, * /ε = and the total configuration energy per monomer in ε units,Econf Econf/Nε. Right: The average valence angle in n-chain, θav.

the average valence angle in n-chain as functions of xn along the composition line 3. These characteristics for all the compositions simulated, along with the differences in energetic characteristics between the rigid-chain and the flexible-chain systems of the same composition (the difference in the total configuration energy per monomer in ε units, E*conf/ε, in contributions to E*conf/ε from DN−DP interactions, and in contribution to Econf * /ε from deformation of valence angles) are presented in the Supporting Information (Figure SI2). When the chain rigidity of the n-decamer is increased from kθ = 85 kJ mol−1 to kθ = 850 kJ mol−1, the volume of the system almost does not change while the average configurational energy Econf (that is, the total potential energy of the system) becomes 7237

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B less negative, particularly in mixtures with high concentration of n-chains. An analysis of different energy contributions (Figure SI2) shows that in the systems containing significant amounts of both n- and p-molecules, this energy change is mostly due to lower energy gain from DN−DP association. For solutions where xn is significantly higher than xp, the increase of valence angle deformation energy and a less negative contribution of Lennard-Jones interactions between nonassociating sites play the decisive role. The difference in DN−DP interactions’ contribution to the energy is partly due to a change in the number of associated n−p pairs as Kass is lower in the rigid-chain systems. However, the average DN−DP interaction energy in a pair of associated stickers is also less negative for the rigid chains: at 365.4 K, for composition lines 3 and 7 (concentrated solutions and mixtures containing no solvent), this value is between −5.77ε and −5.80ε for the flexible chains and between −5.66ε and −5.68ε for the rigid chains. Even more affected by the chain rigidity is the average energy gain from association, Eass, i.e., the difference between the average interaction energies (including DN−DP as well as AN−AP contributions) for a pair of associated monomer units and a pair of nonassociated ones calculated for all pairs with AN−AP distance less than a certain value (here taken to be 1 nm). This gain becomes less negative by 0.4ε−0.6ε when the flexible chains are replaced by the rigid chains. Although the monomer’s propensity to association, in principle, depends on the valence angle, it is noteworthy that there is no significant correlation between the angle distribution and the difference in association equilibria for the systems having different kθ and same composition. For instance, the flexible chains are stretched less than the rigid chains at low xn while at high xn, the reverse is true (see θav in Figure 5, right; see also Figure SI2, bottom, composition lines 2, 3, 5, 6, 7). This may be easily explained: when the fraction pn of associated n-stickers increases with increasing xp/xn ratio, the average valence angle diminishes as association is more thermodynamically advantageous for lower angles (when the dummy atom moves away from the backbone). Nevertheless, the ratio Kass(flexible chains)/ Kass(rigid chains) remains the same, about 1.5. Hence it is reasonable to suggest that, while the chain flexibility is an important factor affecting Kass (and also influencing the system’s average energetic characteristics), the effect is caused mostly by the fact that with more flexible chains, a more thermodynamically favorable local molecular arrangement can be attained. These arrangements result from concerted associative (DN−DP) and nonassociative interactions when formation of an associated pair leads to a regrouping of nonassociating monomer units and solvent molecules around it. For concentrated solutions (composition line 3), Figure 6 presents the average energy gain from association, Eass, and the “apparent average association enthalpy”, ΔHass = −kB

ln(K ass(T1)/K ass(T2)) (1/T1) − (1/T2)

Figure 6. Average energy gain from association for a monomer unit and the average association enthalpy in 365.4 K−406 K temperature range for the n-decamer−p-dimer systems along the composition line 3. The errors are estimated from Kass fluctuations.

content of solvent, the rearrangement involves primarily molecules of the solutes, and forming of an associated pair favors formation of another one by nearby monomer units, so that for mixtures where both xn and xp are high, (−ΔHass) > (−Eass). This cooperative behavior of aggregating molecules probably affects both nonterminal and terminal n-monomers and is likely to be responsible for a significant increase of Kass and (−ΔHass) in such concentrated mixtures. In principle, the molecular rearrangement effect may be treated as preferential solvation in a broader sense including not only interactions with solvent molecules but also nonassociative interactions between monomer units. Density and concentration fluctuations are an additional factor affecting the equilibrium characteristics. It may play a particularly important role at conditions close to phase separation. For the simulated system, such conditions are likely to be expected at high solvent content. Due to the fluctuations, small regions with higher local concentration of associating sites can be formed giving higher total number of associated n−p pairs. This may explain a minor increase of Kass upon dilution in dilute solutions (Figure 3, 406 K, lines 4 and 5). To the best of our knowledge, there is no analytical theory of associating fluids that describes complex cooperative effects observed in our simulations. In the Semenov−Rubinstein (SR) mean field theory, the contribution of bonds between stickers to the system’s partition function is estimated from combinatorial arguments with a correction for spatial restrictions imposed by the associative bonding. Kass is equal to the association parameter λ = vb exp (E̅*/kBT), where vb is the “bond volume” that reflects the steric requirements for bond formation and E̅* is the mean field energy of an associative bond.13,16 Formally, the nonideality factors may be included in vb and E̅ * but at the cost of adjusting these parameters. In contrast to the SR theory that is suited for description of an ideal association equilibrium with a constant value of Kass, SAFT-based approaches6−12 in principle reflect the density and

(3)

calculated from the MD data for T1 = 365.4 K and T2 = 406 K. The difference between ΔHass and Eass and the dependence of ΔHass on the system composition can be explained by the effect of rearrangement of molecular environment on thermodynamics of association. In dilute solutions, the rearrangement of surrounding solvent molecules (partial desolvation of associating sites/groups) makes association less favorable in comparison with a hypothetical case of an environment-independent interaction of two monomer units. In mixtures of low (or zero) 7238

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B concentration dependence of Kass identifying it with the association strength Δ defined as Δ = 4π



⎞ ⎛ ⎛ Φ (12) ⎞ gR (12)⎜⎜exp⎜ − ass ⎟ − 1⎟⎟ kBT ⎠ ⎠ ⎝ ⎝

equilibrium in a wide composition range, it needs corrections owing to the effects of local structure rearrangement and also because the associativity of a site in a flexible chain depends on the local conformation while the conformational distribution depends on composition. 3.2. The n-Decamer−p-Decamer System. Figure 7 shows Kass values for the n-decamer−p-decamer model system together

r 2dr ω1ω2

(4)

where gR(12) is the pair correlation function in the reference (nonassociating) fluid, Φass(12) is the association potential, and ⟨⟩ω1ω2 denotes an unweighted average over all relative orientations of the associating sites. When a hard-sphere system is used as the reference fluid, the association strength can be approximated as9 Δ = κassd3g (d)(exp(εass /kBT ) − 1)

(5)

where g(d) is the contact value of the pair correlation function for the spheres (d is the contact distance), εass is the association energy, and κass is the “dimensionless association volume” that incorporates local orientational/steric effects.9,47 A model system with soft repulsion may also be chosen as a reference. Then, the simplest approximation for Δ is given by Δ = 4πFσ3I where F = exp(εass/kBT) − 1, σ is the geometric parameter of the “reference” Lennard-Jones interaction, and I is the so-called “association kernel”. The latter is an integral that depends on g(r) in the reference fluid and on local anisotropy of the association potential.10,48 We simulated the nonassociating n-decamer+p-dimer+solvent reference fluid and computed the integral (6)

Figure 7. “Apparent” association equilibrium constants for the n-decamer− p-dimer and n-decamer−p-decamer systems with flexible chains from MD simulation at 365.4 and 406 K.

where g(r) is the radial AN−AP correlation function in the nonassociating mixture, and r0 is the distance corresponding to the first peak of g(r). This integral characterizes local correlations between segments in the reference fluid and is not related to the local structure of the associating system. The ratio Kass/G measures the product of (exp(εass/kBT) − 1) by the “local association volume”, the quantity denoted as κass in eq 5 or incorporated as a local anisotropic part in the association kernel I (although generally this part can not be separated analytically from the rest of the kernel). The association energy (that may be approximated as Eass as defined above) is almost independent of the mixture composition (Figure 6). Hence the composition dependence of Kass/G shows the composition dependence of the “local association volume” that is influenced by local association anisotropy and other entropic factors including the restructuring of molecular environment owing to association of monomers. Figure SI3 in the Supporting Information presents Kass/G for different types of n-sites in flexible and rigid chains at 365.4 K (r0 = 0.478 nm) from our MD-data. For terminal sites, Kass/G varies with the solutes’ concentration less strongly than Kass but in relatively concentrated solutions (xn + xp ≥ 0.5), displays a composition dependence pattern similar to Kass, increasing in mixtures that have low concentration of solvent and high concentration of associated pairs, probably due to the aggregation cooperativity described above. For nonterminal sites in mixtures with high content of solvent, Kass/G increases with dilution. This correlates with a decrease of θav that makes association more sterically favorable and is stronger for flexible chains than for rigid chains (see Figure 5 and Figure SI2). Basically, as long as Kass/G remains approximately constant, the SAFT approach is expected to work satisfactorily for the n-decamer−p-dimer system. However, to describe the association

with those for the n-decamer−p-dimer system (in both cases, n-chains are flexible). A strong crowding effect is observed: an increase of the content of nonassociating o-monomers (intermediate inert units in the p-chain) leads to a significant increase of Kass, particularly in concentrated solutions. In mixtures of low concentration of p-component (where the fraction of o-monomers approaches zero), Kass becomes approximately equal to that in the solution containing n-decamers and p-dimers with the same total density of associating monomer units. This association-enhancing effect of inert fragments of p-chains contrasts with the slight decrease of Kass that is observed when small free inert molecules of a solvent are added to n-decamer−p-dimer solutions at moderate concentrations of the solutes. The “apparent” association equilibrium constants for n-sites having different location along the chains’ backbone are presented in Figure SI4. Similar to the n-decamer−p-dimer system, Kass is lower for terminal monomer units than for the nonterminal ones. The difference between terminal and nonterminal sites becomes less significant in mixtures with low solvent content. Figure 8 presents the volume per monomer, v*. For the solvent-rich mixtures, the dependence of v* on the n-chains/ p-chains concentration ratio is weak. It becomes more pronounced in concentrated mixtures. At some low (but nonzero) concentration of solvent, v* has a minimum. This is another indication of a collective molecular rearrangement effect: when a small amount of solvent is added to a mixture of n- and p-decamers, slightly more compact molecular packing becomes possible for the solvated chains, even if the solvent particles are of approximately the same size as the chains’ monomer units. 3.2.1. Crowding Corrections. Crowding effects are usually described by expressing the law of mass action in terms of the activities of reacting species.36,38,49−51 This approach is difficult

G=

∫0

r0

4πr 2g (r )dr

7239

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B

The latter depends on the density and composition of the mixture and hence may reflect nonideality effects on association equilibrium. Nevertheless, because of the symmetry of n-decamer− p-decamer mixtures, here the reference nonassociating system would be the same for all different compositions with the same (xn + xp). Hence the SAFT expression for Δnp may not give the substantial crowding effects observed in our MD-simulation. To describe the crowding effects, we assume that a part of the volume of the system becomes unavailable for the reacting particles due to the presence of the crowder. We introduce an average free volume Vf(i) for each reacting species i that is related to the configurational partition function of the mixture, Zconf = ∏i(Vf(i))Ni. Then, Vexcl(i) = V − Vf(i) is the excluded volume for the species i. In our system, we have three values: Vf(pairs) for bonded pairs, Vf(nfree) for nonassociated n-stickers, andVf(pfree) for nonassociated p-stickers. The free energy of association in this nonideal mixture of bonded pairs and free stickers is given by

Figure 8. Volume per monomer for the n-decamer−p-decamer system along (xn + xp) = const composition lines from MD simulation at 365.4 and 406 K.

to apply to our system that may contain polydisperse clusters of chains (infinite clusters in case of gelation) and there is no definite scheme of reaction between molecular species. On the other hand, the activity of stickers located on the chains is not a well-defined quantity. In this work we revert to an heuristic approach that corrects the SR theory for the effect of excluded volume as this effect plays the leading role in crowding phenomena.37,51 In the SR theory the free energy of association is calculated by comparing an ideal gas mixture of bonded pairs and free stickers with an ideal gas mixture of totally nonreactive stickers. This is readily seen by writing eq 8 of ref 16 in the following form: fst ≡

fst̃ = fst − ρpairs ln

ρ

pairs ρnfree ρpfree

+ ρpfree [ln ρpfree − 1] − ρn [ln ρn − 1] − ρp [ln ρp − 1]

(

v

free

free

)

) = Vf (p

(

v

)

) = V * 1 − 0.8xp v 0* .

v 0.8xp 1 1 = − 0 K ass K0 K 0 v*

where ρn = Nn/V, ρp = Np/V, ρpairs = N pairs /V; = ρn − ρpairs; and ρfree n = ρp − ρpairs are the densities of free stickers in the reactive mixture, E̅* is the bond energy, and vb is the bond volume. Minimization of fst with respect to the number of pairs

(11)

where K0 is the “apparent” association constant in a noncrowded system (at xo → 0, implying xp → 0) that incorporates the nonideality contributions not related to crowding. Assuming that all nonideality factors other than crowding do not change significantly along the composition lines xn + xp = 1 − xsolv = const, we take the mixture with the same xsolv and xp → 0 as the reference noncrowded system. 2) Assume that free stickers of both type “feel” similarly the presence of the crowder (o-monomers) and hence have equal free volume while the bonded pairs are substantially screened within the clusters; their free volume is not sensitive to the crowder and does not depend on its content (at a given fraction of solvent it may be then included in K0). We have Vf(pairs) = const,

= K ass where

Kass = λ = vb exp(E̅ */kBT). Within this model, all changes in the thermal and internal degrees of freedom that occur upon formation of a bond between the stickers have been lumped in parameters E̅ */kBT and vb. As long as the composition dependence of Kass is not included explicitly in the equation through these parameters, the effects of crowding may not be reflected by this model. From eq 7, the mass action law may also be expressed in a different form:16 1 + K ass(ρn + ρp ) 2K ass

v

(

)

Vf (n free) = Vf (p free ) = V 1 − 0.8xp v 0* . Within this ap-

2 1 + K ass (ρn − ρp )2 + 2K assρp + 2K assρn

2K ass

(10)

From eq 10 we have

(7)

ρpairs free free ρn ρp

V · Vf (pairs) Vf (n free)· Vf (p free )

of p-decamers is 1 − 0.8xp v 0* , and the free volume is

ρfree n



= K ass = K 0

Vf (pairs) = Vf (n

(∂fst/∂ρpairs = 0) gives the mass action law,

(9)

where K0 = λ = vbexp(E̅*/kBT). Below we test three empirical expressions for Vf(i) that take into account the effect of excluded volume. 1) Assume that each of the eight nonassociative units (o-monomers) in a p-decamer excludes a constant volume v0 for the bonded pairs and free monomers. The fraction of total volume excluded by the inert part

= ρpairs [ln ρpairs − 1] + ρnfree [ln ρnfree − 1]

ρpair =

Vf (n free)·Vf (p free )

where fst is calculated from eq 6. The free volumes depend on the volume and composition of the mixture. Assuming that they are independent of ρpairs and calculating ∂fst̃ /∂ρpairs = 0, we obtain

Fst kTV

⎡ E̅ * ⎤ − ρpairs ⎢ln v b + ⎥ kBT ⎦ ⎣

V ·Vf (pairs)

proximation, eq 11 is replaced by v0 0.8xp 1 1 = − 1/2 K ass K 01/2 K 01/2 v*

(8)

It has been shown16 that SAFT gives an equation of the same form with Kass replaced by Δnp, the association strength parameter. 7240

(12)

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B Eqs 11 or 12 make possible to calculate Kass along a composition line xn + xp = 1 − xsolv = const if one knows K0 in the noncrowded system and the composition dependence of the ratio v0/v* along this line. Generally, one can expect v0 to be close to the effective size of inert monomer unit that, in turn, is likely to be close to v*. A preliminary estimate of the concentration dependence of Kass in a crowded system can be made using v0 ≈ v* as a first approximation. However, more accurate predictions should take into account that different temperature dependence for v* and v0 is to be expected (v* increases with temperature while v0 may decrease). Knowledge of the composition and temperature dependence of v* and v0 makes possible a quantitative description with the aid of eqs 11 or 12. 3) For a fluid of chains formed by tangent hard spheres, the average free volume may be estimated from SAFT. Assuming that Vf = Vf(pairs) = Vf(nfree) = Vf(pfree) (the same assumption as in the case 1), we obtain kBT ln(V /Vf ) = achain + a seg chain

Figure 9. Reciprocal (left) and the reciprocal square root (right) of the “apparent” association constant vs the density of the crowding (inert) segments of p-chains (0.8xp/v*) in the n-decamer−p-decamer system from MD simulation at 365.4 K. Composition lines are shown in Figure 2.

for low-solvent-content or solvent-free mixtures. The latter seems to be more physically meaningful. These results indicate that eq 12 is somewhat superior at high concentrations of the solutes, favoring the assumption that bonded pairs are screened from the crowder over a wide concentration range. Thus, a simple mean field correction, eq 12, makes it possible to describe the effect of crowding and obtain a quantitative description of the apparent association constant using two adjustable parameters, K0 and v0, that have clear physical meaning. Applying eq 14, we fitted two parameters for each composition line: the effective diameter of the crowding monomers σcr and the association strength (Δnp) at zero content of p-chains. Figure 10 shows that an excellent description of the crowding

(13)

seg

Here a and a are the chain and segment contributions, respectively, to the SAFT residual Helmholtz energy per monomer for the mixture of chains in the solvent.48 For hard chains, aseg is obtained from an equation for the hard sphere mixture.52,53 eqs 9, 10, and 13 then give for the effective constant of the association equilibrium (m)

HS x p K ass(ρn , ρp ) = Δnp(g pp )

(1 − rp)

(m)

HS xn (gnn )

(1 − rn)

HS exp(aseg /kBT )/gnp

(14)

where gHS ij is the contact value of the radial distribution function for the fluid of hard spheres:52,53 3σiiσjj ξ2 1 + 1 − ξ3 σii + σjj (1 − ξ3)2

g ijHS(σij) =

⎡ σiiσjj ⎤2 ξ22 ⎥ + 2⎢ ⎢⎣ σii + σjj ⎥⎦ (1 − ξ3)3

ξk =

πρ 6

(15)

∑ xi(m)riσiik i

and ρ is the total number density of molecules (chains and solvent). Here, the summation is over different ri-mers (including monomeric solvent, rsolv = 1), x(m) is the mole fraction of ri-mer i (not to be confused with xi, the mole fraction of monomer units of type i). The n-decamer chain (rn = 10) is composed of 10 equally sized n-spheres of diameter σn; the p-chain (rp = 10) consists of two terminal p-spheres of diameter σp = σn and eight “crowding” spheres of diameter σcr. Figure 9 shows 1/Kass and 1/K1/2 ass as a function of the density of the inert segments of p-chains. Eqs 11 and 12, respectively, suggest that these functions are linear. The fitted values of K0 and v0 are given in Table SI1. Both equations work fairly well except for solutions poor in solvent and rich in p-chains, where Kass from MD is somewhat lower than eqs 11 or 12 would suggest (1/Kass and 1/K1/2 ass deviate upward from linearity). However, for eq 12, the deviation becomes significant at concentrations much higher than for eq 11. The “effective crowder unit volume” v0 calculated from eq 11 for dilute solutions (where both approximations can be used with reasonable accuracy) is considerably higher than v* (cf. Figure 8), while eq 12 gives v0 that is somewhat lower than v*

Figure 10. “Apparent” association constant in the n-decamer+ p-decamer system vs the fraction of p-chains for constant content of solvent at 365.4 K. Symbols: MD data (for notation, see Figure 9); curves: calculated from eq 14 for a mixture of hard-sphere solvent with tangent hard-sphere chains (n-decamer+p-decamer). The diameters are σp = σn = 0.43 nm for the sticky spheres and σsolv = 0.436759 nm for the solvent. The diameters σcr of the crowding spheres are shown in the graph.

effect can be obtained in this way although the results are rather sensitive to the effective diameter σcr (see Figure SI5). A moderate variation of σcr (close to σp = σn) may lead to opposite effects of the crowder’s concentration on association equilibrium, in line with previous findings.36,37,51 7241

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

The Journal of Physical Chemistry B



4. CONCLUSIONS Our MD results demonstrate the effects of temperature, chain flexibility, and location of the sticker within the chain on the association equilibrium. The terminal stickers are less reactive than the internal stickers; this may be explained by a larger entropy loss upon pinning down a mobile terminal sticker by an associative bond. The difference in reactivity of the terminal and internal stickers diminishes with an increase of the rigidity of chains. The “apparent” association constant has similar dependence on composition for flexible and rigid chains. Increasing the chains’ rigidity hinders association because molecular rearrangement around bonded sites that would produce an energetically favorable local configuration is restrained for less flexible chains. Our data show a significant contribution of the rearrangement of molecular environment to the thermodynamics of association. In dilute solutions, this effect is not likely to favor association, as two sites are partially desolvated when they form an associate. To the contrary, for low content of solvent, the rearrangement of molecular environment involves mostly molecules of the solutes and enhances the association favoring its cooperative character (if two sites are bonded, neighboring sites are more likely to form a bond, too). A minor increase of the “apparent” association constant with dilution in very dilute solutions may be explained by strong fluctuations of concentration: regions of high local concentration of associating sites give a higher total percentage of associated n−p pairs. We show that the presence of nonassociating monomer units in p-chain has a substantial effect on the association equilibrium. This effect disappears in the n-decamer−p-dimer system, in contrast to other manifestations of nonideality, such as cooperativity and preferential solvation, while in n-decamer− p-decamer mixtures, it is much stronger than those nonideality factors that are present also in the noncrowded system. To the best of our knowledge, crowding effects have previously been described only in systems where they are caused by a distinct inert molecular species.36−38,49−51 Such crowding effects can be taken into account by estimating the activity coefficients of the reacting species.38,51 This approach is no longer viable for our system where a multitude of reaction schemes correspond to formation of various polydisperse clusters of chains; an alternative is to account for the activity of stickers in the chains. We show that owing to the symmetry of the reference fluid the standard versions of SAFT may not take into account the crowding observed in our MD simulation. We propose simple empirical modifications of the Semenov−Rubinstein theory that describe these crowding effects. However, a substantial theoretical progress is required to take into account the association cooperativity and composition-dependent conformational effects. Our results may be helpful for better understanding and controlling the behavior of fluids that contain associating chains, e.g., physical gels, because these results show how the association equilibrium depends on the structural details of the molecules. At the same time, it is clear that studies of the crowding effect should be extended and generalized. The chain length plays a role (e.g., the end effects for short molecules). Distribution of associating and nonassociating sites along the chain is likely to be an important factor; particularly, it determines the “effective size” of the crowder that is reflected, e.g., in the value of v0 in eqs 11 or 12. Other issues include local structure and topology of clusters, gelation, and the association kinetics. These issues are among the prospects for future work.

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b12530. “Apparent” association equilibrium constants, v*, energetic characteristics, and θav for n-decamer−p-dimer systems containing flexible (kθ = 85 kJ mol−1) and rigid (kθ = 850 kJ mol−1) decamers at 365.4 K, Kass/G for the five types of n-monomers in the n-decamer−p-dimer system from MD simulation at 365.4 K (r0 = 0.478 nm), K0 and v0 for the n-decamer−p-dimer system, and Kass(xp) according to eq 14 (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work has been supported by St.Petersburg State University (grants ##12.38.199.2014, 12.50.1559.2013) and the Russian Foundation for Basic Research (grant #15-03-04629a). Computational resources of Moscow State University supercomputer center (Lomonosov supercomputer) and of Saint Petersburg State University computational resource center have been used for MD simulations.



REFERENCES

(1) Chiarelli, P.; De Rossi, D. Polyelectrolyte Intelligent Gels: Design and Applications. In Ionic Interactions in Natural and Synthetic Macromolecules; John Wiley & Sons, Inc., 2012; p 581. (2) Intelligent Hydrogels; Sadowski, G., Richtering, W., Eds.; Springer International Publishing, 2013; p 279. (3) Chawla, P.; Srivastava, A. R.; Pandey, P.; Chawla, V. Hydrogels: a journey from diapers to gene delivery. Mini-Rev. Med. Chem. 2014, 14, 154. (4) Manassero, C.; Raos, G.; Allegra, G. Structure of Model Telechelic Polymer Melts by Computer Simulation. J. Macromol. Sci., Part B: Phys. 2005, 44, 855. (5) Tang, S.; Wang, M.; Olsen, B. D. Anomalous Self-Diffusion and Sticky Rouse Dynamics in Associative Protein Hydrogels. J. Am. Chem. Soc. 2015, 137, 3946. (6) Jackson, G.; Chapman, W. G.; Gubbins, K. E. Phase equilibria of associating fluids. Mol. Phys. 1988, 65, 1. (7) Chapman, W. G.; Jackson, G.; Gubbins, K. E. Phase equilibria of associating fluids. Mol. Phys. 1988, 65, 1057. (8) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation-of-state solution model for associating fluids. Fluid Phase Equilib. 1989, 52, 31. (9) Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New reference equation of state for associating liquids. Ind. Eng. Chem. Res. 1990, 29, 1709. (10) Müeller, E. A.; Gubbins, K. E. An Equation of State for Water from a Simplified Intermolecular Potential. Ind. Eng. Chem. Res. 1995, 34, 3662. (11) Economou, I. G.; Donohue, M. D. Chemical, quasi-chemical and perturbation theories for associating fluids. AIChE J. 1991, 37, 1875. (12) Kraska, T. Analytic and Fast Numerical Solutions and Approximations for Cross-Association Models within Statistical Association Fluid Theory. Ind. Eng. Chem. Res. 1998, 37, 4889. (13) Semenov, A. N.; Rubinstein, M. Thermoreversible gelation in solutions of associative polymers. 1. Statics. Macromolecules 1998, 31, 1373.

7242

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243

Article

The Journal of Physical Chemistry B (14) Tanaka, F. Polymer Physics. Application to Molecular Association and Thermoreversible Gelation; Cambridge University Press, 2011. (15) Dormidontova, E.; ten Brinke, G. Phase behavior of hydrogenbonding polymer-oligomer mixtures. Macromolecules 1998, 31, 2649. (16) Tcyrulnikov, S.; Victorov, A. I. Molecular Thermodynamic Modeling of Gelation and Demixing in Solution of Cross-Associating Chains. Macromolecules 2013, 46, 4706. (17) Dormidontova, E. E. Role of competitive PEO-water and waterwater hydrogen bonding in aqueous solution PEO behavior. Macromolecules 2002, 35, 987. (18) Dormidontova, E. E. Influence of End Groups on Phase Behavior and Properties of PEO in Aqueous Solutions. Macromolecules 2004, 37, 7747. (19) Dorn, U.; Enders, S. Mixing enthalpy and liquid−liquid equilibrium of aqueous polyethylene glycol (PEG) solutions. Mol. Phys. 2014, 112, 2310. (20) Fischlschweiger, M.; Enders, S.; Zeiner, T. Solubility calculations of branched and linear amino acids using lattice cluster theory. Mol. Phys. 2014, 112, 2282. (21) Zhang, R.; Shi, T.; An, L.; Sun, Z.; Tong, Z. Conformational Study on Sol−Gel Transition in Telechelic Polyelectrolytes Solutions. J. Phys. Chem. B 2010, 114, 3449. (22) Balabanyan, A. G.; Kramarenko, E. Y.; Ronova, I. A.; Khokhlov, A. R. Monte Carlo study of structure and kinetics of formation of endlinked polymer networks. Polymer 2005, 46, 4248. (23) Baljon, A. R. C.; Flynn, D.; Krawzsenek, D. Numerical study of the gel transition in reversible associating polymers. J. Chem. Phys. 2007, 126, 044907. (24) Cass, M. J.; Heyes, D. M.; Blanchard, R. L.; English, R. J. Simulations and experiments of self-associating telechelic polymer solutions. J. Phys.: Condens. Matter 2008, 20, 335103. (25) Groot, R. D.; Agterof, W. G. M. Monte Carlo study of associative polymer networks. I. Equation of state. J. Chem. Phys. 1994, 100, 1649. (26) Kumar, S. K.; Panagiotopoulos, A. Z. Thermodynamics of Reversibly Associating Polymer Solutions. Phys. Rev. Lett. 1999, 82, 5060. (27) Wang, S.; Chen, C.-C.; Dormidontova, E. E. Reversible association and network formation in 3:1 ligand-metal polymer solutions. Soft Matter 2008, 4, 2039. (28) Hoy, R. S.; Fredrickson, G. H. Thermoreversible associating polymer networks. I. Interplay of thermodynamics, chemical kinetics, and polymer physics. J. Chem. Phys. 2009, 131, 224902. (29) Bianchi, E.; Largo, J.; Tartaglia, P.; Zaccarelli, E.; Sciortino, F. Phase Diagram of Patchy Colloids: Towards Empty Liquids. Phys. Rev. Lett. 2006, 97, 168301. (30) de las Heras, D.; Tavares, J. M.; Telo da Gama, M. M. Phase diagrams of binary mixtures of patchy colloids with distinct numbers of patches: the network fluid regime. Soft Matter 2011, 7, 5615. (31) Tsuchida, E.; Osada, Y.; Ohno, H. Formation of interpolymer complexes. J. Macromol. Sci., Part B: Phys. 1980, 17, 683. (32) Khutoryanskiy, V. V. Hydrogen-bonded interpolymer complexes as materials for pharmaceutical applications. Int. J. Pharm. 2007, 334, 15. (33) Luo, K.; Yin, J.; Song, Z.; Cui, L.; Cao, B.; Chen, X. Biodegradable Interpolyelectrolyte Complexes Based on Methoxy Poly(ethylene glycol)-b-poly(α,l-glutamic acid) and Chitosan. Biomacromolecules 2008, 9, 2653. (34) Xu, D.; Craig, S. L. Scaling Laws in Supramolecular Polymer Networks. Macromolecules 2011, 44, 5465. (35) Minton, A. P. The effect of volume occupancy upon the thermodynamic activity of proteins: some biochemical consequences. Mol. Cell. Biochem. 1983, 55, 119. (36) Hall, D.; Minton, A. P. Macromolecular crowding: qualitative and semiquantitative successes, quantitative challenges. Biochim. Biophys. Acta, Proteins Proteomics 2003, 1649, 127. (37) Zhou, H.-X.; Rivas, G.; Minton, A. P. Macromolecular crowding and confinement: Biochemical, biophysical, and potential physiological consequences. Annu. Rev. Biophys. 2008, 37, 375.

(38) Minton, A. P. Quantitative assessment of the relative contributions of steric repulsion and chemical interactions to macromolecular crowding. Biopolymers 2013, 99, 239. (39) Ellis, R. J. Macromolecular crowding: obvious but underappreciated. Trends Biochem. Sci. 2001, 26, 597. (40) Chebotareva, N. A.; Eronina, T. B.; Roman, S. G.; Poliansky, N. B.; Muranov, K. O.; Kurganov, B. I. Effect of crowding and chaperones on self-association, aggregation and reconstitution of apophosphorylase b. Int. J. Biol. Macromol. 2013, 60, 69. (41) de Gennes, P. G. Scaling Concepts in Polymer Physics; Cornell University Press, 1979. (42) van der Spoel, D.; van Maaren, P. J. The Origin of Layer Structure Artifacts in Simulations of Liquid Water. J. Chem. Theory Comput. 2006, 2, 1. (43) Abraham, M. J.; van der Spoel, D.; Lindahl, E.; Hess, B.and the GROMACS development team. GROMACS User Manual, version 5.0.4; Royal Institution of Technology and Uppsala University: Sweden, 2014; www.gromacs.org. (44) Lee, H.; de Vries, A. H.; Marrink, S.-J.; Pastor, R. W. A CoarseGrained Model for Polyethylene Oxide and Polyethylene Glycol: Conformation and Hydrodynamics. J. Phys. Chem. B 2009, 113, 13186. (45) 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. (46) Voevodin, V. V.; Zhumatiy, S. A.; Sobolev, S. I.; Antonov, A. S.; Bryzgalov, P. A.; Nikitenko, D. A.; Stefanov, K. S.; Voevodin, V. V. Practice of ″Lomonosov″ Supercomputer. Open Systems 2012, 36. (47) Kontogeorgis, G. M.; Voutsas, E. C.; Yakoumis, I. V.; Tassios, D. P. An Equation of State for Associating Fluids. Ind. Eng. Chem. Res. 1996, 35, 4310. (48) Dufal, S.; Lafitte, T.; Haslam, A. J.; Galindo, A.; Clark, G. N. I.; Vega, C.; Jackson, G. The A in SAFT: developing the contribution of association to the Helmholtz free energy within a Wertheim TPT1 treatment of generic Mie fluids. Mol. Phys. 2015, 113, 948. (49) Rivas, G.; Fernández, J. A.; Minton, A. P. Direct observation of the enhancement of noncooperative protein self-assembly by macromolecular crowding: Indefinite linear self-association of bacterial cell division protein FtsZ. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 3150. (50) Zhou, H.-X. Protein folding and binding in confined spaces and in crowded solutions. J. Mol. Recognit. 2004, 17, 368. (51) Hu, Z.; Jiang, J.; Rajagopalan, R. Effects of Macromolecular Crowding on Biochemical Reaction Equilibria: A Molecular Thermodynamic Perspective. Biophys. J. 2007, 93, 1464. (52) Boublík, T. Hard-Sphere Equation of State. J. Chem. Phys. 1970, 53, 471. (53) Henderson, D.; Trokhymchuk, A.; Woodcock, L. V.; Chan, K.-Y. Simulation and approximate formulae for the radial distribution functions of highly asymmetric hard sphere mixtures. Mol. Phys. 2005, 103, 667.

7243

DOI: 10.1021/acs.jpcb.5b12530 J. Phys. Chem. B 2016, 120, 7234−7243