- Email: [email protected]
Phase behavior of aqueous polyacrylic acid solutions using atomistic molecular dynamics simulations of model oligomers
Accepted Manuscript Phase behavior of aqueous polyacrylic acid solutions using atomistic molecular dynamics simulations of model oligomers Ratna S. Katiyar, Prateek K. Jha PII:
S0032-3861(17)30254-9
DOI:
10.1016/j.polymer.2017.03.007
Reference:
JPOL 19503
To appear in:
Polymer
Received Date: 15 December 2016 Revised Date:
3 March 2017
Accepted Date: 4 March 2017
Please cite this article as: Katiyar RS, Jha PK, Phase behavior of aqueous polyacrylic acid solutions using atomistic molecular dynamics simulations of model oligomers, Polymer (2017), doi: 10.1016/ j.polymer.2017.03.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
ACCEPTED MANUSCRIPT
Phase behavior of aqueous polyacrylic acid solutions using atomistic molecular dynamics simulations of model oligomers
RI PT
Ratna S. Katiyar and Prateek K. Jha* Department of Chemical Engineering, IIT Roorkee, Uttarakhand, India, 247667
*Corresponding author, Email: [email protected], Telephone: +91 1332 284810, Fax: +91 1332 276535, 273560
M AN U
§ Electronic supporting information available
SC
5
Abstract
We have performed fully atomistic molecular dynamics simulations of aqueous solutions of a weak, pH-responsive polyelectrolyte, polyacrylic acid (PAA). Model oligomers of PAA of
TE D
10
different tacticities, molecular weights, degrees of deprotonation, and deprotonation patterns are simulated with water molecules. Deprotonation of PAA chains that occurs with an
EP
increase in pH results in an increase in Coulomb repulsion between chain segments on one
15
AC C
hand, and a non-monotonic change in the hydrogen bonding between chain segments on the other hand. Consequently, at the single chain level, PAA chains are stretched at higher pH values, where the amount of stretching varies with chain tacticity. For the multiple chains case, PAA forms aggregates at higher concentrations, which are relatively denser and contain lesser water (solid-like) at lower pH than compared to higher pH (liquid-like). Such phase transitions of PAA aggregates with pH has possible implications in the design of pH20
responsive polyelectrolytes for applications in drug delivery. Keywords: atomistic simulations; polyacrylic acid; phase behavior 1
ACCEPTED MANUSCRIPT
1. Introduction According to the IUPAC definition, polyelectrolyte is a “polymer composed of macromolecules in which a substantial portion of the constitutional units contains ionic or
5
RI PT
ionisable groups, or both”[1]. Depending on the strength of the ionisable group, they can be classified as strong or weak. Strong polyelectrolytes contain strong acidic/basic groups (e.g. sulfonate, hydroxyl, etc.) that remain ionised at all practical pH conditions. On the other
SC
hand, weak polyelectrolytes contain weak acidic/basic groups (e.g., carboxyl, amide, etc.) that undergo a transition from a neutral state to a charged state with change in pH. This
10
M AN U
transition from neutral to charged state occurs due to deprotonation of acidic groups with an increase in pH or protonation of basic groups with a decrease in pH. In general, polyelectrolytes have better aqueous solubility than neutral polymers. Moreover, they exhibit interesting self-assembly behavior driven by electrostatic forces.[2,3] For example,
TE D
oppositely charged polyelectrolytes self-assemble to form polyelectrolyte complexes, which can either be solid-like (precipitate) or a liquid/gel-like (coacervate) depending on the pH and 15
ionic strength (salt concentration) of the medium[4–8]. Similar self-assembly behavior can
EP
also be observed in systems containing identical polyelectrolytes, by either varying the degree of ionization by changing pH[9,10] or screening the Coulomb repulsion between the by increasing
salt
AC C
chains
concentration[11,12].
Such
self-assembly behavior
of
polyelectrolytes finds use in a range of applications; some recent areas of applications include 20
fuel cell technology[13,14], drug delivery[15–17], tissue engineering[18,19], and optoelectronic devices[20]. Physical properties of polyelectrolytes are fundamentally interesting and have intrigued scientists working in the area of biology and soft matter from last several decades.[2,21,22] At the single chain level (dilute concentrations), polyelectrolyte chains are stretched and have 2
ACCEPTED MANUSCRIPT high persistence length due to Coulomb repulsion between chain segments. For charge densities beyond a threshold value or for the case of multivalent salts, counterions condense on the polyelectrolyte backbone (Manning condensation[23]), resulting in chain collapse. Concentrated solutions of polyelectrolytes possess strong ion-ion correlations resulting in the formation of strongly segregated structures. These behaviors have been captured in field-
RI PT
5
theoretical studies[24,25] and coarse-grained simulations[26,27], which however have mainly focused on strong polyelectrolytes. Such studies of pH-responsive, weak polyelectrolytes are
SC
relatively few[28,29]. Although the field-theoretical and coarse-grained simulation approaches succeed in achieving a qualitative understanding of polyelectrolyte behavior, they do not adequately account for the chemical details of polyelectrolytes and thus do not provide
M AN U
10
a molecular-level insight of the underlying physics. In particular, physical interactions that vary strongly with the polymer chemistry (e.g., hydrogen-bonding) are only captured at a rudimentary level and often requires a priori knowledge about the magnitudes of such
15
TE D
interactions. Atomistic simulations, on the other hand, include the necessary chemical details and thus capture such interactions in a predictive manner, that is, without a prior knowledge about these interactions. Despite this obvious advantage, studies of polyelectrolyte systems
EP
using atomistic simulations are relatively scarce, mainly because of the computational
AC C
expense of these simulations. While atomistic simulations of realistic large molecular weight polymers remains impossible, simulations of reasonably long oligomers that mimic the 20
behavior of real polymers are only now becoming feasible[30–33]. Atomistic simulations of long oligomers coupled with systematic coarse-graining approaches[34–36] now provide a pathway to study polyelectrolyte systems at atomic resolution, and is being actively pursued. The degree of deprotonation of polyelectrolytes at a given pH can be obtained using the Henderson-Hasselbalch equation[37], that is, in the case of polyacids,
3
ACCEPTED MANUSCRIPT =
1 . 1 + 10
(1)
Here, is the fraction of deprotonated groups (degree of deprotonation) and is the
negative logarithm (base 10) of the acid dissociation constant, . There are two major issues
RI PT
in the use of this equation for polyelectrolytes. First, the values of polyelectrolyte
solutions are a function of the polyelectrolyte molecular weight and salt concentration[5] and are generally not available. Second, the Henderson-Hasselbalch equation can only provide
SC
5
the average fraction of charged groups and provide no information on the actual positions of
M AN U
these charged groups on the polyelectrolyte chain, which further keep changing due to dynamic nature of protonation-deprotonation equilibrium. This is demonstrated in fig. 1, where the two PAA chains have the same deprotonation fraction but different deprotonation 10
patterns. groups are closer to each other in fig. 1(b) compared to fig. 1(a), which results in relatively larger intra-chain electrostatic repulsion. Therefore, the structure and
TE D
properties of these two PAA chain conformations, both of which can be obtained at the same pH, may be strikingly different. Further, the Henderson-Hasselbalch equation provide no clue about which of these conformations would be more likely. It ignores the fact that the deprotonation free energy of a group is generally expected to be larger if a
EP
15
AC C
group is present on the neighbouring monomer(s), since the deprotonation would lead to high electrostatic repulsion. Thus, we need to resort to molecular simulations that provide detailed insights into the effect of pH on the polyelectrolyte behavior. Yet another advantage of molecular simulations is their predictive nature; simulations can be performed to study 20
behavior of novel molecules that have not been synthesized to access their potential for a given application.[33] This can save experimental effort and expense, while screening molecules for different applications.
4
RI PT
ACCEPTED MANUSCRIPT
Two distinct approaches can be used in the atomistic simulations of weak polyelectrolytes.
M AN U
5
SC
Figure 1: PAA chains with same fraction of deprotonation, but different deprotonation patterns.
First is the conventional molecular dynamics (MD) simulations performed in the canonical (NVT) ensemble or isothermal-isobaric (NPT) ensemble simulations of polyelectrolytes, performed for different degrees of deprotonation and deprotonation patterns[38,39], which
10
TE D
include the pH-induced changes in the average deprotonation but do not account for the dynamic nature of protonation-deprotonation equilibrium. Second is the constant-pH molecular dynamics approaches[40,41] that have originally been developed for proteins and
EP
can be extended for simulations of weak polyelectrolytes as done in a recent study[9]. Constant-pH simulations also accounts for the dynamic protonation-deprotonation processes
15
AC C
and therefore can predict the values and the titration curve of polyelectrolytes. However, there are several issues in the practical implementation of these methods. For instance, the use of an artificial barrier potential, a reference “chemically similar” compound, incorporation of explicit titratable water, and presence of multiple titration sites on a molecule remains problematic from a fundamental standpoint and computationally difficult.[40] More importantly, such methods should be adequately able to sample 2 20
possible deprotonation patterns for a polyelectrolyte chain containing ionisable groups, 5
ACCEPTED MANUSCRIPT which results in a very large number of possible conformations in a solution of polyelectrolytes. We therefore adopt the first approach and vary both the degree of deprotonation and deprotonation patterns to understand the effect of these on the phase behavior. We expect that the main conclusions of our study will be the same for constant-pH simulations, provided they are able to adequately sample all possible conformations.
RI PT
5
However, the effects of dynamic changes in the deprotonation patterns are not captured in our
SC
study.
In this paper, we discuss atomistic MD simulations of aqueous solutions of PAA oligomers
10
M AN U
(chains) of different tacticities, concentrations, degrees of deprotonation, and deprotonation patterns. Simulation protocol and methods employed in the analysis of simulation data are discussed in the next section. This is followed by the Results and Discussion section, which include a detailed discussion of the effects of various parameters (tacticity, PAA concentration, degree of deprotonation, deprotonation patterns) on the single-chain
study and an outlook for future studies.
2. Methods
EP
15
TE D
characteristics and the solution phase behavior. We conclude with the main outcomes of this
AC C
Simulations are performed in the GROMACS 4.6.7 simulation package (Groningen Machine for Chemical Simulation)[42] using the CHARMM27 force field[43]. The potential energy of the system in the CHARMM27 force field includes both bonded and non-bonded 20
contributions. Bonded contributions are of the form:
6
ACCEPTED MANUSCRIPT !"#
!&'(#
− , + $ % − % , +
8+ (+#
) ,1 + cos01 2 − 3 45
6 7 − 7 , ,
(2)
RI PT
+
"*("+'#
i.e., contains bond stretching, angle bending, proper dihedrals, and improper dihedral energies, summed over all possible bonds, angles, proper dihedrals, and improper dihedrals of
SC
PAA chains and water molecules. Values of force constants 0, , $, , ), , 6, 4,
5
M AN U
multiplicities (1 ) and phase shifts (3 ), and equilibrium values of bond distances ( , ), angles (% , ), and out-of-plane angles (7 , ) are provided in the topology file. All the above quantities remain constant throughout the simulations. However, magnitudes of the bond distances , angles % , dihedral angle 2 , and out-of-plane angles 7 are computed on-the-fly
of the form: B 8# B 8#
:C
@
B 8# B 8#
=8!,: 1 − 2< > A+ =: 4EF
:C
G G: , =:
(3)
and apply between all pairs of atoms in the simulation box, excluding pairs of atoms on the
AC C
10
?
=8!,: 9: ;< > =:
EP
TE D
during simulations and vary with time. Nonbonded contributions to the potential energy are
same chain separated by less than three bonds. The first term is the van der Waals contribution represented using a 6-12 Lennard Jones (LJ) potential, where energy well depth 9: and “minimum distance” =8!,: (defined as the finite distance at which LJ potential is zero) values of all pairs are provided in the input topology and remain constant throughout 15
the simulation. The distance between pairs, =: , vary in the course of the simulation. The
7
ACCEPTED MANUSCRIPT
second term is the electrostatic energy of pairs of atoms of partial charges G and G: at distance =: , where F is vacuum permittivity.
PAA oligomers (chains) of different tacticities (isotactic, syndiotactic, and atactic) containing
5
RI PT
20 monomers (repeating units) and various degrees of deprotonation and deprotonation patterns are constructed using Gaussian’09 software. Two different kinds of deprotonation patterns are studied for partially ionized PAA. In random deprotonation, groups are
SC
randomly distributed on the PAA chain. In end deprotonation, groups are
deprotonated starting from one end of PAA chain. GROMACS compatible topologies of
10
M AN U
model oligomers for force calculations are generated using the automatic topology building tool Swiss Param[44] (Molecular modelling group, Swiss Institute of Bioinformatics, Lausanne, Switzerland)[42]. Single chains and multiple chains (for various concentrations) of PAA are simulated with simple point charge (SPC) water in a cubic simulation box with
TE D
periodic boundary conditions. Appropriate number of HI counterions are added, as needed, to maintain the overall electroneutrality of the simulation box. The simulation box size (6 1L 15
for all simulations) is significantly larger than the chain contour length (≈ 5 1L) in order to
simulations.
EP
avoid periodicity artefacts. The density of the simulation box is approximately 1 O/QQ for all
AC C
Simulations have been performed in three steps. In the first step, energy minimization is performed using the steepest descent technique, with an initial step size of 0.001 1L and 20
adaptive step size control, until the maximum force on each atom become smaller than 10 R/01L ⋅ LTU4. In the second step, MD equilibration is performed for ∼ 60 1W, beyond which the properties of interest converge around an average value. Finally, in the third step, MD production is performed for 5 1W with a typical sampling frequency of 0.5 W. Both equilibration and production simulations are performed in the NVT ensemble with a
8
ACCEPTED MANUSCRIPT reference temperature of 298 K. Nosé-Hoover thermostat is used for temperature control, with a coupling constant of 0.5 W in both equilibration and production simulations. Electrostatic calculations are performed using the particle mesh Ewald (PME) scheme with a Fourier spacing of 0.1 1L, van der Waals cut-off of 1 1L, and cubic interpolation. Bonds containing hydrogen atoms are constrained using a fourth-order LINCS (Linear Constraint
RI PT
5
Solver) algorithm.
SC
Number of oligomer-oligomer contacts with time, H 0X4 (normalized by the total number of atoms of oligomers) is monitored to track the equilibrium process. Here, a “contact” is
10
M AN U
counted if any atom of any oligomer is within a chosen threshold distance (0.6 1L) from any atom of another oligomer. Similar concept is used to obtain the fraction of “condensed” counterions (Y ), which is defined as the total number of “contacts” between the oxygen
atoms of groups and HI counterions, divided by the total number of oxygen atoms of
groups. In this case, a “contact” is counted if any oxygen atom of any group is
15
TE D
within a threshold distance of 0.22 1L from any H I counterion. Please note that the choice of the above threshold distances are somewhat arbitrary due to inherent ambiguities in the
EP
definition of a “contact” or a “condensed counterion”. The number of hydrogen bonds between the two species is defined as the number of donor-acceptor pairs that are within a
AC C
threshold distance 0.35 1L and the hydrogen-donor-acceptor angle less than 30 degrees.
Inter-molecular PAA-PAA radial distribution functions (RDFs), O0=4, are defined as the
20
normalized probability of finding any atom of a PAA chain at a distance between = and = + 3= from any atom of another PAA chain, where bin width, 3= = 0.002 1L, and
normalization is performed with such probability of an ideal gas. Radius of gyration ([& ) of
PAA is computed in two ways – average [& of individual chains and the [& of all chains. While the former represents the size of individual chains, the latter represents the size of PAA
9
ACCEPTED MANUSCRIPT
aggregates. Water content of PAA aggregates are quantified by computing the [& of all chains (entire PAA aggregate) and the solvent-accessible surface area[45] (SASA) computed using a 0.14 1L probe radius. A decrease in [& of all chains signifies compaction of
aggregates, i.e., an increase in the local density of aggregates (∝ H/[&] , H being the total number of atoms in PAA chains). Equivalently, a decrease in SASA/monomer (SASA
RI PT
5
divided by the total number of monomers in PAA chains) indicates a decrease in water content or an increase in the compaction (local density) of aggregates. Except H 0X4, all the
10
3.1. Single-Chain Simulations
M AN U
3. Results and Discussion
SC
above quantities are time-averaged over trajectories sampled during production simulations.
In order to validate the force-field employed in our simulations, we first compare the results
TE D
of single-chain simulations of atactic PAA solutions for different degrees of deprotonation () with that reported by Sulatha and Natarajan[46]. Unlike the other simulations reported in this work that uses the simulation protocol given in section 2, we use the simulation protocol used in their study[46] for these single-chain simulations. Since the actual tacticity patterns
EP
15
used in the simulations are not provided in their paper, we have performed the simulations for
AC C
three different tacticity patterns (Tables S1, S2, S3 in the Supporting Information). We observe that the [& averaged over the three chains increases with the increase in the degree of deprotonation (fig. S1 in the Supporting Information), though the magnitude of [&
20
significantly vary for the three chains. Further, the typical chain conformations at different values of show a transition from a collapsed to a stretched conformation with an increase in (fig. S2 in the Supporting Information). All these results compare well with the results
reported in their paper, but we do not attempt a direct comparison due to differences between
10
ACCEPTED MANUSCRIPT tacticity patterns in the two studies. It is worth noting that the force-field employed in our study has also been used in a previous molecular dynamics study of PAA-PDMAEMA complexation.[5]
5
RI PT
In order to investigate the effects of tacticity on the single-chain statistics in more detail, we then perform single chain simulations for PAA of different tacticities (isotactic, syndiotactic, and atactic) for the neutral ( = 04 and fully deprotonated ( = 1) cases. As shown in fig. 2 and Table 1, while the PAA chains are always more stretched for the fully deprotonated case
SC
than compared to the neutral case, [& for a given varies with tacticity in the order syndiotactic > atactic > isotactic. This difference between [& values for different tacticity is
explained using the fact that all groups are on the same side of PAA chain in the
M AN U
10
isotactic case, and adjacent groups are on the opposite sides of PAA chain in the
syndiotactic case. Therefore, the magnitude of short-ranged hydrophobic attraction for the
TE D
neutral case is larger for isotactic PAA than compared to syndiotactic PAA, resulting in more collapsed conformations. If a similar reasoning is applied to the case of fully deprotonated 15
chains, electrostatic repulsion should be higher in the isotactic case than compared to
EP
syndiotactic case, which would result in larger stretching in the isotactic case than compared to syndiotactic case. This is in contradiction with the observed result that syndiotactic PAA is
AC C
more stretched compared to isotactic PAA at = 1. It can be justified using the fact that the
presence of groups on the same side of PAA chain in the isotactic case results in the
20
bending of chain in order to keep the groups as far away as possible. Such bending is
not favored in the syndiotactic case as the groups are on both sides of the chain, resulting in a straight chain. Yet another explanation for larger stretching in the syndiotactic case compared to the isotactic case is the difference in the magnitudes of counterion condensation in the two cases. Since the fraction of condensed counterions in the isotactic 25
case is significantly higher than that in the syndiotactic case (Table 1), it would have 11
ACCEPTED MANUSCRIPT relatively lower electrostatic repulsion between chain segments due to the charge screening effect of counterions. The former reason appears more convincing, given the relatively large differences in the values of [& in the isotactic and syndiotactic cases. It is useful to compare
these [& values with their upper bound, that is, the theoretical estimate in the rodlike limit.
The − equilibrium bond distances along the backbone is ≈ 0.15 1L, which results in a theoretical estimate of [&,+
"
RI PT
5
= 39 /√12 ≈ 1.69 1L in the rodlike limit[47]. As
expected, all values of [& shown in Table 2 are lower than [&,+
"
and maximum [& is
SC
obtained for = 1 case of syndiotactic PAA that has most stretched configuration (fig. 2d) in
AC C
EP
TE D
M AN U
fig. 2.
10
Figure 2: Typical single chain conformations of PAA for neutral and fully deprotonated cases for three types of tacticity: (a)-(b) isotactic, (c)-(d) syndiotactic, (e)-(f) atactic PAA chain. Carbon, oxygen, and hydrogen atoms in PAA are shown in grey, red, and white colors, respectively. Water molecules and counterions are not shown. 15
12
ACCEPTED MANUSCRIPT
Table 1: Comparison of [& values of neutral ( = 0) and fully deprotonated states ( = 1) of single PAA chain of different tacticity. Fraction of “condensed” counterions (Y ) is also shown for the = 1 case. The values in brackets are the standard deviation of the last decimal. =1
[& (nm)
[& (nm)
Isotactic
0.83(4)
0.97(1)
0.28(6)
Syndiotactic
1.26(7)
1.42(2)
0.20(6)
Atactic
1.06(1)
1.16(1)
0.27(6)
SC
3.2. Multiple-Chain Simulations
Y
RI PT
Tacticity
5
=0
Four different concentrations of PAA are simulated by varying the number of PAA chains in
M AN U
a simulation box of 6 1L size containing SPC water. PAA concentration in the four systems, defined in moles per liter (M) of monomers are: (a) 0.615 M (4 PAA chains), (b) 1.223 M (8 PAA chains), (c) 2.456 M (16 PAA chains), and (d) 3.689 M (24 PAA chains). 10
To decide the amount of equilibration time required for the simulations, we monitor the
TE D
number of oligomer-oligomer contacts with time, which appear to converge after 60 1W, as established by simulations performed for a longer time of 90 1W in fig. 3 for the neutral
EP
syndiotactic case. Equilibration times of this order had also been reported in earlier atomistic simulation studies on polymer aggregation for systems of similar sizes.[33,36,48] Note that our estimate of equilibrium time is based on simulations performed over significantly long
AC C
15
time permitted within the current computational capabilities of atomistic simulations. The possibility of further aggregation or crystallization occurring over time scales that are much longer than that permitted by existing computational capabilities (e.g., time scales higher than 1 `W) can neither be established nor be ruled out. Interestingly, PAA aggregates span the 20
simulation box after equilibration in fig. 3, which at first glance, appear to be an artefact of periodic boundary condition used in our simulations. However, since the chain contour length is significantly smaller than the box size, PAA spanning the simulation box is less likely to be 13
ACCEPTED MANUSCRIPT a periodicity artifact, and is more likely to be reminiscent of a “gel”-like structure reported in
M AN U
SC
RI PT
a previous study.[33]
TE D
5
Figure 3: Number of oligomer-oligomer contacts with time for neutral syndiotactic PAA of concentration 3.689 M. Red dotted line indicate the convergence of oligomer-oligomer contacts around an average value. Simulation snapshots at the start of simulations (left) and after equilibration (right) are also shown. PAA chains are shown in red color; water molecules are not shown.
fully deprotonated cases. As expected, neutral cases for all tacticities formed more compact
AC C
10
EP
Next, we compare the aggregation behavior of PAA of different tacticity for the neutral and
aggregates than compared to the fully deprotonated cases, as evidenced by larger intermolecular RDF peaks for the neutral cases in fig. 4. However, the difference between the RDF for different tacticities at a given are not as substantial as the difference in [& of
single-chains for different tacticities (Table 1). This is because the groups of a PAA 15
chain in an aggregate interact not only with the groups on the same chain, but also
with groups on other PAA chains, which in turn may be oriented in any direction with respect to the PAA chain. Therefore, changes in the single-chain conformations due to
14
ACCEPTED MANUSCRIPT tacticity play little role in determining the aggregation behavior. It is worth pointing out that the first peak around = ∼ 0.2 1L for the neutral case in fig. 4 is due to hydrogen bonding between groups. Roughly, ∼ 0.25 − hydrogen bonds/monomer were observed for all tacticities. This peak is absent for the fully deprotonated cases, as the
groups do not hydrogen bond. The second peak occurring around = ∼ 1 1L represent an
RI PT
5
aggregation length scale of the order of [& (Table 1). This peak is much broader for the neutral case than compared to the fully deprotonated case, implying more PAA aggregation
Figure 4: Intermolecular PAA-PAA RDF for neutral ( = 0) and fully deprotonated ( = 1) PAA of different tacticities at a concentration of 3.689 M.
AC C
10
EP
TE D
M AN U
SC
in the neutral case.
Figure 5 shows the equilibrium configurations of syndiotactic PAA chains at various concentrations for = 0.4 with random deprotonation (Table S4 in the Supporting 15
Information). As expected, aggregation increases with increase in PAA concentration. Although the counterions have higher mobility than compared to PAA, they mostly reside within the PAA aggregate, thus neutralizing its overall charge. Further, when compared to = 0 case in fig. 3 (fig. 5d has the same concentration as fig. 3), the aggregate appears less 15
ACCEPTED MANUSCRIPT compact. The effect of PAA concentration on the aggregation behavior is quantified in terms of several measurable parameters in Table 2. The average [& of individual chains increases with increase in PAA concentration, implying an increase in the stretching of individual chains. On the other hand, the overall [& of all chains first increases (until 2.456 M in Table
2) and then decrease (for 3.689 M in Table 2). The initial increase in [& can be attributed to
RI PT
5
an increase in the size of PAA aggregates with increase in the number of PAA chains, and the final decrease (after [& ≈ 3 1L) can be understood to result from the compaction of
SC
aggregates. Although the [& values show a slightly non-monotonic trend with concentration,
increases uniformly with increase in concentration. This is further confirmed by the fact that the water content of aggregate, characterized by the SASA/monomer roughly decreases with increase in concentration, with small anomaly at the two intermediate concentrations where
EP
TE D
the SASA/monomer is nearly same.
AC C
10
M AN U
the local density of aggregates (∝ H/[&] , H being the total number of atoms in PAA chains)
16
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
TE D
5
Figure 5: Simulation snapshots (after equilibration) of syndiotactic PAA with = 0.4 for different concentration with random deprotonation (Table S4 in the Supporting Information). PAA chains and counterions are shown in red and green colors, respectively; water molecules are not shown.
AC C
EP
Table 2: Radius of gyration (for single chain and all chains), SASA/monomer, and fraction of condensed ions, of syndiotactic PAA with = 0.4 for different concentration with random deprotonation (Table S4 in the Supporting Information). The values in brackets are the standard deviation of the last decimal. PAA
[& (nm)
Concentration Individual
SASA/monomer All
2
(nm )
Y
chains
chains
0.615 M
1.05 (2)
2.0(4)
0.65(1)
0.26(5)
1.223 M
1.185(9)
2.8(2)
0.567(6)
0.30(4)
2.456 M
1.296(5)
3.06(6)
0.573(6)
0.32(3)
3.689 M
1.309(3)
2.68(3)
0.542(4)
0.32(2)
10
17
ACCEPTED MANUSCRIPT
The effect of on PAA aggregation behavior is more clearly demonstrated in fig. 6, which
shows that the compaction of PAA aggregate decreases with increase in . Electrostatic
repulsion between PAA chains increases with increase in resulting in the formation of less
compact structures. This is further confirmed by a monotonic increase in SASA/monomer with an increase in at the given concentration (Table 3). Such monotonic trend is not
RI PT
5
observed in the [& of individual chains and [& of all chains in Table 3. However, roughly
speaking, both the [& of individual chains [& of all chains increase with increase in ,
SC
signifying an increase in the stretching of individual chains and a decrease in overall
M AN U
compaction of aggregates, respectively.
AC C
EP
TE D
10
18
5
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 6: Simulation snapshots (after equilibration) for syndiotactic PAA of concentration 3.689 M for different values of with random deprotonation (Table S4 in the Supporting Information). PAA chains and counterions are shown in red and green colors, respectively; water molecules are not shown.
19
ACCEPTED MANUSCRIPT
[& (nm)
Y
SASA/monomer 2
(nm )
All
chains
chains
0
1.272(5)
2.66(4)
0.375(4)
0.2
1.268(4)
2.76(5)
0.511(4)
0.4
1.309(3)
2.68(3)
0.542(4)
0.6
1.342(2)
2.70(1)
0.569(4)
0.39(2)
0.8
1.320(3)
2.91(4)
0.603(3)
0.43(2)
1.0
1.357(2)
2.99(1)
RI PT
Individual
-
0.19(2) 0.32(2)
SC
M AN U
5
Table 3: Radius of gyration (for single chain and all chains), SASA/monomer, and fraction of condensed ions, of syndiotactic PAA at different degrees of deprotonation with random deprotonation (Table S4 in the Supporting Information). The values in brackets are the standard deviation of the last decimal.
0.651(3)
0.45(2)
More information about the structure of PAA aggregates is revealed by the intermolecular PAA-PAA RDF (fig. 7). Second peak of this RDF occurring around ∼ 1 1L is higher and
10
TE D
broader for smaller values of indicating an increase in PAA aggregation, as already observed in fig. 6. However, the first peak occurring around ∼ 0.15 1L (inset of fig. 7) do
not exhibit the same trend due to rather non-monotonic changes in intermolecular hydrogen-
EP
bonding with . While cannot hydrogen bond with , it can hydrogen bond with , which contribute to the first peak along with − hydrogen bonding. In
15
AC C
the extreme case of fully deprotonated chain ( = 1 case in fig. 7), the first peak is missing as there is no possibility of hydrogen bonding between groups. For partially ionized cases
( = 0.2,0.4,0.6,0.8 case in fig. 7), we observe a peak at = ≈ 0.16 1L due to the combined
effect of hydrogen bonding between and groups, along with −
hydrogen bonding. In the other extreme of a neutral chain ( = 0 case in fig. 7), the first peak
occurs at a slightly larger value of = ≈ 0.17 1L as the groups are not present. The
20
hydrogen bonding data in Table 4 shows that the − hydrogen bonding 20
ACCEPTED MANUSCRIPT
decreases with due to a reduction in the number of groups. On the other hand,
− hydrogen bonding first increases and reaches a maximum at an intermediate value of , followed by a decrease at even higher value of . This is because the number of
RI PT
possible − pairs is maximum at an intermediate value.
M AN U
SC
5
Figure 7: Intermolecular PAA-PAA RDF for syndiotactic PAA of 3.689 M concentration for different values of with random deprotonation. O0=4 for small values of = are shown in the inset.
and
and
and
0.2 0.13(1)
0.01(1)
0.4 0.18(2)
and
2.456 M
and
0.18(1)
0.05(1)
0.03(1)
0.21(1)
0.6 0.17(1)
< 0.01
0.8 0.12(1)
< 0.01
0
1.223 M
EP
0.615 M
TE D
Table 4: Hydrogen bonding per monomer for syndiotactic PAA of different concentration for different f with random deprotonation. The values in brackets are the standard deviation of the last decimal. Values c 0.01 are not reported, since the standard deviation becomes comparable or higher than the mean.
AC C
10
−
0.05(3)
3.689 M
and
and
and
0.23(1)
0.04(1)
0.24(1)
0.05(1)
0.02(0)
0.22(1)
0.02(1)
0.25(1)
< 0.01
0.20(1)
< 0.01
0.19(0)
< 0.01
0.21(1)
< 0.01
0.07(1)
< 0.01
0.08(0)
< 0.01
0.11(0)
< 0.01
−
0.19(1)
−
0.25(2)
−
0.22(1)
21
ACCEPTED MANUSCRIPT Further insight into the intermolecular interactions can be obtained by a closer inspection of two nearby PAA chains in the concentrated system. Figures 8(a), 8(b), 8(c) show the zoomed view of simulation snapshots of the figs. 6(a), 6(c), and 6(f), for the neutral, partially deprotonated ( = 0.4), and fully deprotonated cases, respectively. As evident from fig. 8,
the distance between PAA chains at a given concentration increases with an increase in due
RI PT
5
to an increase in electrostatic repulsion (or a decrease in hydrophobic attraction) between PAA chains. However, as noted above, hydrogen bonding between and groups
SC
in the partially ionized case (fig. 8b) works together with − hydrogen bonding,
resulting in substantial attraction between PAA chains. Finally, while such hydrogen bonding is absent in the fully deprotonated case (fig. 8c), segments of PAA chain can still show an
M AN U
10
effective attraction due to “counterion bridging” as the HI counterions sandwiched between the two PAA chains experience electrostatic attraction from both the PAA chains. The likelihood of counterion bridging increases with increasing either the concentration or the
15
TE D
degree of deprotonation of PAA, as the fraction of condensed ions (Y ) increases both with increasing concentration (Table 2) and with increase in (Table 3). However, this
“counterion bridging” (or “salt-bridging”) effect is much weaker compared to hydrogen
EP
bonding interactions as the counterions have high entropy and prefers to disperse in the
AC C
solution. Therefore, PAA chains are mostly separated from each other in the fully ionized case (fig. 6f). 20
22
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
TE D
EP
5
Figure 8: Zoomed view of simulation snapshot for syndiotactic PAA of 3.689 M concentration for different values of with random deprotonation, showing the typical interactions between two nearby chains. Average distance between the center of mass of chains are 2.31 d 0.02 1L, 2.60 d 0.02 1L and 3.20 d 0.01 1L for = 0, 0.4, and 1, respectively. Carbon, oxygen, and hydrogen atoms in PAA are shown in grey, red, and white colors, respectively. Counterions are shown in green color. Water molecules are not shown.
10
AC C
Next, we compare the aggregation behavior of the partially deprotonated case for random deprotonation and end deprotonation at the same concentration. As shown in figs. 9(a) and 9(b), the overall aggregation behavior appear similar. This is further confirmed by the comparison of their intermolecular PAA-PAA RDF in fig. 9c. However, the dense regions formed by the hydrogen bonding of groups is more populated and prominent in the end deprotonation case (fig. 9b) compared to the random deprotonation case (fig. 9a). In fact, 15
for realistic chains (much longer than the oligomers studied here), a micelle-like formation (appearance of a dense core of hydrophobic segments) is expected in the end deprotonation 23
ACCEPTED MANUSCRIPT case due to the hydrogen bonding between the protonated ends of chains. On the other hand, a network-like formation with no such dense core is expected in the random deprotonation case of realistic chains, as the protonated groups are randomly distributed along the chains. Such differences will however be only apparent at intermediate values of (e.g., = 0.2 or
0.4). At lower (higher) values of , most of the chain is hydrophobic (hydrophilic) resulting
RI PT
5
in the formation of a dense precipitate (dilute solution). This can be substantiated by comparing the magnitudes of SASA/monomer for the random deprotonation case (Table 3)
SC
and end deprotonation case (Table 5). SASA/monomer magnitudes are significantly lower for the end deprotonation case than compared to the random deprotonation case for = 0.2 and
M AN U
0.4, which implies the formation of denser structures in the end deprotonation case. However, the system size and the length of chain is not sufficient enough to obtain a clear distinction of
EP
TE D
micelle-like and network-like structures.
AC C
10
24
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
10
EP
AC C
5
Figure 9: Simulation snapshot of syndiotactic PAA of 3.689 M concentration for (a) random deprotonation (Table S4 in the supporting information) and (b) end deprotonation (Table S5 in the supporting information), at fixed = 0.4, Water molecules and counterions are not shown here; aliphatic backbone carbon, group, and group are shown in red, blue, and green color, respectively. (c) Intermolecular PAA-PAA RDF for random (dashed line) and end deprotonation (bold line).
25
ACCEPTED MANUSCRIPT Table 5: Radius of gyration (for single chain and all chains), SASA/monomer (hydrophobic, hydrophilic, and total), and fraction of condensed ions, of syndiotactic PAA at different degrees of deprotonation with end deprotonation (Table S5 in the Supporting Information). The values in brackets are the standard deviation of the last decimal.
[& (nm)
Y
SASA/monomer 2
(nm )
All
chains
chains
0
1.272(5)
2.66(4)
0.375(4)
0.2
1.258(5)
2.83(7)
0.467(4)
0.4
1.271(4)
2.96(6)
0.522(4)
0.6
1.331(3)
2.79(4)
0.585(3)
0.8
1.339(3)
2.60(3)
0.624(4)
0.42(2)
1.0
1.357(2)
2.99(1)
0.651(3)
0.45(2)
RI PT
Individual
-
0.28(3) 0.38(2)
M AN U
SC
0.41(2)
5
The differences in phase behavior observed with changes in (fig. 6) and deprotonation pattern (fig. 8) has important implication in the use of polyelectrolytes as stimuli-responsive
10
TE D
materials for various applications. For example, for pH-responsive polyelectrolytes used as drug carriers, it is critical to attain substantial change in the drug release rate with change in pH. This can be realized due to an increase in the drug diffusivity within the carrier with a
EP
decrease in polyelectrolyte volume fraction, or equivalently, an increase in the water volume fraction. The water content of the PAA aggregates can be estimated by the solvent-accessible
15
AC C
surface area (SASA). In figure 10, we plot the SASA normalized by its value at = 0 (measure of change in water content with ) against . As shown in the plot, SASA changes
by ∼ 3 times in the entire range of , which should give rise to substantial change in the drug diffusivity within the aggregate with pH.
26
RI PT
ACCEPTED MANUSCRIPT
SC
M AN U
5
Figure 10: Solvent accessible surface area (SASA) divided by its value at = 0 against for random deprotonation and end deprotonation cases for syndiotactic PAA of 3.689 M concentration.
4. Conclusion
In summary, the results of atomistic MD simulations provide a detailed molecular picture of
10
TE D
the PAA aggregation behavior as a function of the degree of deprotonation or pH. PAA aggregation increases with PAA concentration and an increase in its degree of deprotonation (pH). Further, the nature of aggregation depends on the deprotonation pattern; end
EP
deprotonation and random deprotonation showed the formation of micelle-like and networklike aggregates, respectively. Interactions that were found relevant to the aggregation
15
AC C
behavior included the hydrogen bonding between − and −
groups, electrostatic repulsion between groups, and counterion (HI ) bridging
between the groups. Finally, the water content of aggregates expressed in terms of solvent-accessible surface area showed a 3-fold increase between the neutral and fully deprotonated limits.
20
Several important aspects studied in this paper require further investigation. First and foremost, our simulations are limited to model oligomers of PAA containing 20 monomers. 27
ACCEPTED MANUSCRIPT Studying the effect of chain flexibility on the aggregation behavior will require the simulations of much longer chains, well beyond current computational limitations of atomistic simulations. A useful compromise can be made using a systematic coarse-grained model parametrized using the atomistic simulation results described here, which may then be applied to study longer realistic chains. Second, although an increase in water content with
RI PT
5
pH signify a decrease in the drug diffusivity when PAA is used as a drug carrier, simulations using actual drug molecules are required to also include the PAA-drug interactions on the
SC
drug release. Third, we have confined our model to aqueous PAA solutions with no additional salt. Effects of the inclusion of salt species, particularly those present in the body fluids, are thus not captured. We are currently working along these lines of investigation.
M AN U
10
Acknowledgments
TE D
This work is supported by the INSPIRE award from the Department of Science and Technology (DST), India (no. IFA14/ENG-72) and a faculty initiation grant from IIT Roorkee obtained by P.K.J. References
M. Nič, J. Jirát, B. Košata, A. Jenkins, A. McNaught, eds., IUPAC Compendium of
AC C
[1]
EP
15
Chemical Terminology, IUPAC, Research Triagle Park, NC, 2009. doi:10.1351/goldbook.
20
[2]
A. V. Dobrynin, Theory and simulations of charged polymers: From solution properties to polymeric nanomaterials, Curr. Opin. Colloid Interface Sci. 13 (2008) 376–388. doi:10.1016/j.cocis.2008.03.006.
28
ACCEPTED MANUSCRIPT [3]
M.A. Cohen Stuart, B. Hofs, I.K. Voets, A. de Keizer, Assembly of polyelectrolytecontaining block copolymers in aqueous media, Curr. Opin. Colloid Interface Sci. 10 (2005) 30–36. doi:10.1016/j.cocis.2005.04.004. E. Kizilay, A.B. Kayitmazer, P.L. Dubin, Complexation and coacervation of
RI PT
[4]
polyelectrolytes with oppositely charged colloids, Adv. Colloid Interface Sci. 167
5
(2011) 24–37. doi:10.1016/j.cis.2011.06.006.
P.K. Jha, P.S. Desai, J. Li, R.G. Larson, pH and salt effects on the associative phase
SC
[5]
separation of oppositely charged polyelectrolytes, Polymers (Basel). 6 (2014) 1414–
10
[6]
M AN U
1436.
D. Priftis, M. Tirrell, Phase behaviour and complex coacervation of aqueous polypeptide solutions, Soft Matter. 8 (2012) 9396–9405. doi:10.1039/c2sm25604e. J. van der Gucht, E. Spruijt, M. Lemmers, M.A. Cohen Stuart, Polyelectrolyte
TE D
[7]
complexes: Bulk phases and colloidal systems, J. Colloid Interface Sci. 361 (2011) 407–422. doi:10.1016/j.jcis.2011.05.080. [8]
A. Kudlay, A. V. Ermoshkin, M.O. De La Cruz, Complexation of oppositely charged
EP
15
AC C
polyelectrolytes: Effect of ion pair formation, Macromolecules. 37 (2004) 9231–9241. doi:10.1021/ma048519t.
[9]
B.H. Morrow, G.F. Payne, J. Shen, pH-Responsive Self-Assembly of Polysaccharide through a Rugged Energy Landscape, J. Am. Chem. Soc. 137 (2015) 13024–13030.
20
doi:10.1021/jacs.5b07761.
29
ACCEPTED MANUSCRIPT [10] M.O. de la Cruz, L. Belloni, M. Delsanti, J.P. Dalbiez, O. Spalla, M. Drifford, Precipitation of highly-charged polyelectrolyte solutions in the presence of multivalent salts, J. Chem. Phys. 103 (1995) 5781–5791. doi:10.1063/1.470459.
5
RI PT
[11] A. Kudlay, M.O. De la Cruz, Precipitation of oppositely charged polyelectrolytes in salt solutions, J. Chem. Phys. 120 (2004) 404–412. doi:10.1063/1.1629271.
[12] A. Kundagrami, M. Muthukumar, Theory of competitive counterion adsorption on
SC
flexible polyelectrolytes: Divalent salts, J. Chem. Phys. 128 (2008) 244901.
M AN U
doi:10.1063/1.2940199.
[13] X. Guo, J. Fang, T. Watari, K. Tanaka, H. Kita, K.I. Okamoto, Novel sulfonated 10
polyimides as polyelectrolytes for fuel cell application. 2. Synthesis and proton conductivity of polyimides from 9,9-bis(4-aminophenyl)fluorene-2,7-disulfonic acid,
TE D
Macromolecules. 35 (2002) 6707–6713. doi:10.1021/ma020260w. [14] V.K. Sachan, A. Devi, R.S. Katiyar, R.K. Nagarale, P.K. Bhattacharya, Proton transport properties of sulphanilic acid tethered poly(methyl vinyl ether-alt-maleic anhydride)-PVA blend membranes, Eur. Polym. J. 56 (2014) 45–58.
EP
15
AC C
doi:10.1016/j.eurpolymj.2014.04.009. [15] C.S. Peyratout, L. Dahne, Tailor-made polyelectrolyte microcapsules: From multilayers to smart containers, Angew. Chemie - Int. Ed. 43 (2004) 3762–3783. doi:10.1002/anie.200300568.
20
[16] S.Y. Nie, W.J. Lin, N. Yao, X.D. Guo, L.J. Zhang, Drug Release from pH-Sensitive Polymeric Micelles with Different Drug Distributions: Insight from Coarse-Grained Simulations, ACS Appl. Mater. Interfaces. 6 (2014) 17668–17678. doi:10.1021/am503920m. 30
ACCEPTED MANUSCRIPT [17] D. Schmaljohann, Thermo- and pH-responsive polymers in drug delivery, Adv. Drug Deliv. Rev. 58 (2006) 1655–1670. doi:10.1016/j.addr.2006.09.020. [18] P. Coimbra, P. Ferreira, H.C. de Sousa, P. Batista, M.A. Rodrigues, I.J. Correia, M.H.
5
RI PT
Gil, Preparation and chemical and biological characterization of a pectin/chitosan polyelectrolyte complex scaffold for possible bone tissue engineering applications, Int. J. Biol. Macromol. 48 (2011) 112–118. doi:10.1016/j.ijbiomac.2010.10.006.
SC
[19] H. Almeida, M.H. Amaral, P. Lobão, Temperature and pH stimuli-responsive
polymers and their applications in controlled and selfregulated drug delivery, J. Appl.
10
M AN U
Pharm. Sci. 2 (2012) 01–10. doi:10.7324/JAPS.2012.2609.
[20] C. V. Hoven, A. Garcia, G.C. Bazan, T.Q. Nguyen, Recent applications of conjugated polyelectrolytes in optoelectronic devices, Adv. Mater. 20 (2008) 3793–3810.
TE D
doi:10.1002/adma.200800533.
[21] A. V. Dobrynin, M. Rubinstein, Theory of polyelectrolytes in solutions and at surfaces, Prog. Polym. Sci. 30 (2005) 1049–1118. doi:10.1016/j.progpolymsci.2005.07.006. [22] C. Holm, J.F. Joanny, K. Kremer, R.R. Netz, P. Reineker, C. Seidel, T.A. Vilgis, R.G.
EP
15
AC C
Winkler, Polyelectrolyte Theory, Adv. Polym. Sci. 166 (2004) 3–17. doi:10.1007/b11349.
[23] G.S. Manning, Limiting Laws and Counterion Condensation in Polyelectrolyte Solutions I. Colligative Properties, J. Chem. Phys. 51 (1969) 924.
20
doi:10.1063/1.1672157. [24] M. Muthukumar, Theory of counter-ion condensation on flexible polyelectrolytes: Adsorption mechanism, J. Chem. Phys. 120 (2004) 9343. doi:10.1063/1.1701839.
31
ACCEPTED MANUSCRIPT [25] Q. Wang, T. Taniguchi, G.H. Fredrickson, Self-consistent field theory of polyelectrolyte systems, J. Phys. Chem. B. 108 (2004) 6733–6744. doi:10.1021/jp037053y.
5
RI PT
[26] F. Alarcón, G. Pérez-Hernández, E. Pérez, A. Gama Goicochea, Coarse-grained simulations of the salt dependence of the radius of gyration of polyelectrolytes as
models for biomolecules in aqueous solution, Eur. Biophys. J. 42 (2013) 661–672.
SC
doi:10.1007/s00249-013-0915-z.
[27] S. Liu, M. Muthukumar, Langevin dynamics simulation of counterion distribution
10
doi:10.1063/1.1476930.
M AN U
around isolated flexible polyelectrolyte chains, J. Chem. Phys. 116 (2002) 9975.
[28] A. Laguecir, S. Ulrich, J. Labille, N. Fatin-Rouge, S. Stoll, J. Buffle, Size and pH effect on electrical and conformational behavior of poly(acrylic acid): Simulation and
TE D
experiment, Eur. Polym. J. 42 (2006) 1135–1144. doi:10.1016/j.eurpolymj.2005.11.023.
[29] G. Berghold, P. van der Schoot, C. Seidel, Equilibrium charge distribution on weak
EP
15
AC C
polyelectrolytes, J. Chem. Phys. 107 (1997) 8083–8088. doi:10.1063/1.475071. [30] R. Chockalingam, U. Natarajan, Self-association behaviour of atactic polymethacrylic acid in aqueous solution investigated by atomistic molecular dynamics simulations, Mol. Simul. 41 (2015) 1110–1121. doi:10.1080/08927022.2014.947481.
20
[31] S.H. Min, S.K. Kwak, B.-S. Kim, Atomistic simulation for coil-to-globule transition of poly(2-dimethylaminoethyl methacrylate), Soft Matter. 11 (2015) 2423–2433. doi:10.1039/C4SM02242D.
32
ACCEPTED MANUSCRIPT [32] Y. Liu, V.S. Bryantsev, M.S. Diallo, W.A. Goddard III, PAMAM Dendrimers Undergo pH Responsive Conformational Changes without Swelling, J. Am. Chem. Soc. 131 (2009) 2798–2799. doi:10.1021/ja8100227.
5
RI PT
[33] P.K. Jha, R.G. Larson, Assessing the efficiency of polymeric excipients by atomistic molecular dynamics simulations, Mol. Pharm. 11 (2014) 1676–1686. doi:10.1021/mp500068w.
SC
[34] D. Reith, B. Müller, F. Müller-Plathe, S. Wiegand, How does the chain extension of poly (acrylic acid) scale in aqueous solution? A combined study with light scattering
10
M AN U
and computer simulation, J. Chem. Phys. 116 (2002) 9100. doi:10.1063/1.1471901. [35] D.D. Hsu, W. Xia, S.G. Arturo, S. Keten, Systematic method for thermomechanically consistent coarse-graining: A universal model for methacrylate-based polymers, J.
TE D
Chem. Theory Comput. 10 (2014) 2514–2527. doi:10.1021/ct500080h. [36] W. Huang, R. Ramesh, P.K. Jha, R.G. Larson, A Systematic Coarse-Grained Model for Methylcellulose Polymers: Spontaneous Ring Formation at Elevated Temperature, Macromolecules. 49 (2016) 1490–1503. doi:10.1021/acs.macromol.5b02373.
EP
15
AC C
[37] P.W. Atkins, J. De Paula, Physical chemistry for the life sciences, W.H. Freeman and Co., New York, 2011.
[38] N. Hoda, R.G. Larson, Explicit- and Implicit-Solvent Molecular Dynamics Simulations of Complex Formation between Polycations and Polyanions,
20
Macromolecules. 42 (2009) 8851–8863. doi:10.1021/ma901632c.
33
ACCEPTED MANUSCRIPT [39] R. Chockalingam, U. Natarajan, Self-association behaviour of atactic polymethacrylic acid in aqueous solution investigated by atomistic molecular dynamics simulations, Mol. Simul. 41 (2015) 1110–1121. doi:10.1080/08927022.2014.947481.
5
RI PT
[40] S. Donnini, F. Tegeler, G. Groenhof, H. Grubmüller, Constant pH molecular dynamics in explicit solvent with λ-dynamics, J. Chem. Theory Comput. 7 (2011) 1962–1978. doi:10.1021/ct200061r.
SC
[41] U. Börjesson, P.H. Hünenberger, Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines, J. Chem. Phys. 114 (2001)
10
M AN U
9706–9719. doi:10.1063/1.1370959.
[42] S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M.R. Shirts, J.C. Smith, P.M. Kasson, D. Van Der Spoel, B. Hess, E. Lindahl, GROMACS 4.5: A highthroughput and highly parallel open source molecular simulation toolkit,
TE D
Bioinformatics. 29 (2013) 845–854. doi:10.1093/bioinformatics/btt055. [43] E. Lindahl, P. Bjelkmar, P. Larsson, M.A. Cuendet, B. Hess, Implementation of the charmm force field in GROMACS: Analysis of protein stability effects from
EP
15
correction maps, virtual interaction sites, and water models, J. Chem. Theory Comput.
AC C
6 (2010) 459–466. doi:10.1021/ct900549r.
[44] V. Zoete, M.A. Cuendet, A. Grosdidier, O. Michielin, SwissParam: A fast force field generation tool for small organic molecules, J. Comput. Chem. 32 (2011) 2359–2368.
20
doi:10.1002/jcc.21816.
34
ACCEPTED MANUSCRIPT [45] F. Eisenhaber, P. Lijnzaad, P. Argos, C. Sander, M. Scharf, The double cubic lattice method: Efficient approaches to numerical integration of surface area and volume and to dot surface contouring of molecular assemblies, J. Comput. Chem. 16 (1995) 273–
5
RI PT
284. doi:10.1002/jcc.540160303. [46] M.S. Sulatha, U. Natarajan, Origin of the difference in structural behavior of
poly(acrylic acid) and poly(methacrylic acid) in aqueous solution discerned by
SC
explicit-solvent explicit-ion MD simulations, Ind. Eng. Chem. Res. 50 (2011) 11785– 11796. doi:10.1021/ie2014845.
Sons, Inc., New York, USA, 2002. doi:10.1002/0471224510.fmatter_indsub. [48] W. Huang, I.S. Dalal, R.G. Larson, Analysis of solvation and gelation behavior of methylcellulose using atomistic molecular dynamics simulations, J. Phys. Chem. B.
EP
TE D
118 (2014) 13992–14008. doi:10.1021/jp509760x.
AC C
10
M AN U
[47] I. Teraoka, Polymer Solutions: An Introduction to Physical Properties, John Wiley &
35
ACCEPTED MANUSCRIPT
Highlights
EP
TE D
M AN U
SC
RI PT
Atomistic MD simulations of polyacrylic acid (PAA) have been performed. We study the effects of tacticity, concentration, and pH on the phase behavior. Aggregates form at high concentrations and their water content increase with pH. Results of this study may be useful for the design of pH-responsive drug carriers.
AC C
• • • •
S0032-3861(17)30254-9
DOI:
10.1016/j.polymer.2017.03.007
Reference:
JPOL 19503
To appear in:
Polymer
Received Date: 15 December 2016 Revised Date:
3 March 2017
Accepted Date: 4 March 2017
Please cite this article as: Katiyar RS, Jha PK, Phase behavior of aqueous polyacrylic acid solutions using atomistic molecular dynamics simulations of model oligomers, Polymer (2017), doi: 10.1016/ j.polymer.2017.03.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
ACCEPTED MANUSCRIPT
Phase behavior of aqueous polyacrylic acid solutions using atomistic molecular dynamics simulations of model oligomers
RI PT
Ratna S. Katiyar and Prateek K. Jha* Department of Chemical Engineering, IIT Roorkee, Uttarakhand, India, 247667
*Corresponding author, Email: [email protected], Telephone: +91 1332 284810, Fax: +91 1332 276535, 273560
M AN U
§ Electronic supporting information available
SC
5
Abstract
We have performed fully atomistic molecular dynamics simulations of aqueous solutions of a weak, pH-responsive polyelectrolyte, polyacrylic acid (PAA). Model oligomers of PAA of
TE D
10
different tacticities, molecular weights, degrees of deprotonation, and deprotonation patterns are simulated with water molecules. Deprotonation of PAA chains that occurs with an
EP
increase in pH results in an increase in Coulomb repulsion between chain segments on one
15
AC C
hand, and a non-monotonic change in the hydrogen bonding between chain segments on the other hand. Consequently, at the single chain level, PAA chains are stretched at higher pH values, where the amount of stretching varies with chain tacticity. For the multiple chains case, PAA forms aggregates at higher concentrations, which are relatively denser and contain lesser water (solid-like) at lower pH than compared to higher pH (liquid-like). Such phase transitions of PAA aggregates with pH has possible implications in the design of pH20
responsive polyelectrolytes for applications in drug delivery. Keywords: atomistic simulations; polyacrylic acid; phase behavior 1
ACCEPTED MANUSCRIPT
1. Introduction According to the IUPAC definition, polyelectrolyte is a “polymer composed of macromolecules in which a substantial portion of the constitutional units contains ionic or
5
RI PT
ionisable groups, or both”[1]. Depending on the strength of the ionisable group, they can be classified as strong or weak. Strong polyelectrolytes contain strong acidic/basic groups (e.g. sulfonate, hydroxyl, etc.) that remain ionised at all practical pH conditions. On the other
SC
hand, weak polyelectrolytes contain weak acidic/basic groups (e.g., carboxyl, amide, etc.) that undergo a transition from a neutral state to a charged state with change in pH. This
10
M AN U
transition from neutral to charged state occurs due to deprotonation of acidic groups with an increase in pH or protonation of basic groups with a decrease in pH. In general, polyelectrolytes have better aqueous solubility than neutral polymers. Moreover, they exhibit interesting self-assembly behavior driven by electrostatic forces.[2,3] For example,
TE D
oppositely charged polyelectrolytes self-assemble to form polyelectrolyte complexes, which can either be solid-like (precipitate) or a liquid/gel-like (coacervate) depending on the pH and 15
ionic strength (salt concentration) of the medium[4–8]. Similar self-assembly behavior can
EP
also be observed in systems containing identical polyelectrolytes, by either varying the degree of ionization by changing pH[9,10] or screening the Coulomb repulsion between the by increasing
salt
AC C
chains
concentration[11,12].
Such
self-assembly behavior
of
polyelectrolytes finds use in a range of applications; some recent areas of applications include 20
fuel cell technology[13,14], drug delivery[15–17], tissue engineering[18,19], and optoelectronic devices[20]. Physical properties of polyelectrolytes are fundamentally interesting and have intrigued scientists working in the area of biology and soft matter from last several decades.[2,21,22] At the single chain level (dilute concentrations), polyelectrolyte chains are stretched and have 2
ACCEPTED MANUSCRIPT high persistence length due to Coulomb repulsion between chain segments. For charge densities beyond a threshold value or for the case of multivalent salts, counterions condense on the polyelectrolyte backbone (Manning condensation[23]), resulting in chain collapse. Concentrated solutions of polyelectrolytes possess strong ion-ion correlations resulting in the formation of strongly segregated structures. These behaviors have been captured in field-
RI PT
5
theoretical studies[24,25] and coarse-grained simulations[26,27], which however have mainly focused on strong polyelectrolytes. Such studies of pH-responsive, weak polyelectrolytes are
SC
relatively few[28,29]. Although the field-theoretical and coarse-grained simulation approaches succeed in achieving a qualitative understanding of polyelectrolyte behavior, they do not adequately account for the chemical details of polyelectrolytes and thus do not provide
M AN U
10
a molecular-level insight of the underlying physics. In particular, physical interactions that vary strongly with the polymer chemistry (e.g., hydrogen-bonding) are only captured at a rudimentary level and often requires a priori knowledge about the magnitudes of such
15
TE D
interactions. Atomistic simulations, on the other hand, include the necessary chemical details and thus capture such interactions in a predictive manner, that is, without a prior knowledge about these interactions. Despite this obvious advantage, studies of polyelectrolyte systems
EP
using atomistic simulations are relatively scarce, mainly because of the computational
AC C
expense of these simulations. While atomistic simulations of realistic large molecular weight polymers remains impossible, simulations of reasonably long oligomers that mimic the 20
behavior of real polymers are only now becoming feasible[30–33]. Atomistic simulations of long oligomers coupled with systematic coarse-graining approaches[34–36] now provide a pathway to study polyelectrolyte systems at atomic resolution, and is being actively pursued. The degree of deprotonation of polyelectrolytes at a given pH can be obtained using the Henderson-Hasselbalch equation[37], that is, in the case of polyacids,
3
ACCEPTED MANUSCRIPT =
1 . 1 + 10
(1)
Here, is the fraction of deprotonated groups (degree of deprotonation) and is the
negative logarithm (base 10) of the acid dissociation constant, . There are two major issues
RI PT
in the use of this equation for polyelectrolytes. First, the values of polyelectrolyte
solutions are a function of the polyelectrolyte molecular weight and salt concentration[5] and are generally not available. Second, the Henderson-Hasselbalch equation can only provide
SC
5
the average fraction of charged groups and provide no information on the actual positions of
M AN U
these charged groups on the polyelectrolyte chain, which further keep changing due to dynamic nature of protonation-deprotonation equilibrium. This is demonstrated in fig. 1, where the two PAA chains have the same deprotonation fraction but different deprotonation 10
patterns. groups are closer to each other in fig. 1(b) compared to fig. 1(a), which results in relatively larger intra-chain electrostatic repulsion. Therefore, the structure and
TE D
properties of these two PAA chain conformations, both of which can be obtained at the same pH, may be strikingly different. Further, the Henderson-Hasselbalch equation provide no clue about which of these conformations would be more likely. It ignores the fact that the deprotonation free energy of a group is generally expected to be larger if a
EP
15
AC C
group is present on the neighbouring monomer(s), since the deprotonation would lead to high electrostatic repulsion. Thus, we need to resort to molecular simulations that provide detailed insights into the effect of pH on the polyelectrolyte behavior. Yet another advantage of molecular simulations is their predictive nature; simulations can be performed to study 20
behavior of novel molecules that have not been synthesized to access their potential for a given application.[33] This can save experimental effort and expense, while screening molecules for different applications.
4
RI PT
ACCEPTED MANUSCRIPT
Two distinct approaches can be used in the atomistic simulations of weak polyelectrolytes.
M AN U
5
SC
Figure 1: PAA chains with same fraction of deprotonation, but different deprotonation patterns.
First is the conventional molecular dynamics (MD) simulations performed in the canonical (NVT) ensemble or isothermal-isobaric (NPT) ensemble simulations of polyelectrolytes, performed for different degrees of deprotonation and deprotonation patterns[38,39], which
10
TE D
include the pH-induced changes in the average deprotonation but do not account for the dynamic nature of protonation-deprotonation equilibrium. Second is the constant-pH molecular dynamics approaches[40,41] that have originally been developed for proteins and
EP
can be extended for simulations of weak polyelectrolytes as done in a recent study[9]. Constant-pH simulations also accounts for the dynamic protonation-deprotonation processes
15
AC C
and therefore can predict the values and the titration curve of polyelectrolytes. However, there are several issues in the practical implementation of these methods. For instance, the use of an artificial barrier potential, a reference “chemically similar” compound, incorporation of explicit titratable water, and presence of multiple titration sites on a molecule remains problematic from a fundamental standpoint and computationally difficult.[40] More importantly, such methods should be adequately able to sample 2 20
possible deprotonation patterns for a polyelectrolyte chain containing ionisable groups, 5
ACCEPTED MANUSCRIPT which results in a very large number of possible conformations in a solution of polyelectrolytes. We therefore adopt the first approach and vary both the degree of deprotonation and deprotonation patterns to understand the effect of these on the phase behavior. We expect that the main conclusions of our study will be the same for constant-pH simulations, provided they are able to adequately sample all possible conformations.
RI PT
5
However, the effects of dynamic changes in the deprotonation patterns are not captured in our
SC
study.
In this paper, we discuss atomistic MD simulations of aqueous solutions of PAA oligomers
10
M AN U
(chains) of different tacticities, concentrations, degrees of deprotonation, and deprotonation patterns. Simulation protocol and methods employed in the analysis of simulation data are discussed in the next section. This is followed by the Results and Discussion section, which include a detailed discussion of the effects of various parameters (tacticity, PAA concentration, degree of deprotonation, deprotonation patterns) on the single-chain
study and an outlook for future studies.
2. Methods
EP
15
TE D
characteristics and the solution phase behavior. We conclude with the main outcomes of this
AC C
Simulations are performed in the GROMACS 4.6.7 simulation package (Groningen Machine for Chemical Simulation)[42] using the CHARMM27 force field[43]. The potential energy of the system in the CHARMM27 force field includes both bonded and non-bonded 20
contributions. Bonded contributions are of the form:
6
ACCEPTED MANUSCRIPT !"#
!&'(#
− , + $ % − % , +
8+ (+#
) ,1 + cos01 2 − 3 45
6 7 − 7 , ,
(2)
RI PT
+
"*("+'#
i.e., contains bond stretching, angle bending, proper dihedrals, and improper dihedral energies, summed over all possible bonds, angles, proper dihedrals, and improper dihedrals of
SC
PAA chains and water molecules. Values of force constants 0, , $, , ), , 6, 4,
5
M AN U
multiplicities (1 ) and phase shifts (3 ), and equilibrium values of bond distances ( , ), angles (% , ), and out-of-plane angles (7 , ) are provided in the topology file. All the above quantities remain constant throughout the simulations. However, magnitudes of the bond distances , angles % , dihedral angle 2 , and out-of-plane angles 7 are computed on-the-fly
of the form: B 8# B 8#
:C
@
B 8# B 8#
=8!,: 1 − 2< > A+ =: 4EF
:C
G G: , =:
(3)
and apply between all pairs of atoms in the simulation box, excluding pairs of atoms on the
AC C
10
?
=8!,: 9: ;< > =:
EP
TE D
during simulations and vary with time. Nonbonded contributions to the potential energy are
same chain separated by less than three bonds. The first term is the van der Waals contribution represented using a 6-12 Lennard Jones (LJ) potential, where energy well depth 9: and “minimum distance” =8!,: (defined as the finite distance at which LJ potential is zero) values of all pairs are provided in the input topology and remain constant throughout 15
the simulation. The distance between pairs, =: , vary in the course of the simulation. The
7
ACCEPTED MANUSCRIPT
second term is the electrostatic energy of pairs of atoms of partial charges G and G: at distance =: , where F is vacuum permittivity.
PAA oligomers (chains) of different tacticities (isotactic, syndiotactic, and atactic) containing
5
RI PT
20 monomers (repeating units) and various degrees of deprotonation and deprotonation patterns are constructed using Gaussian’09 software. Two different kinds of deprotonation patterns are studied for partially ionized PAA. In random deprotonation, groups are
SC
randomly distributed on the PAA chain. In end deprotonation, groups are
deprotonated starting from one end of PAA chain. GROMACS compatible topologies of
10
M AN U
model oligomers for force calculations are generated using the automatic topology building tool Swiss Param[44] (Molecular modelling group, Swiss Institute of Bioinformatics, Lausanne, Switzerland)[42]. Single chains and multiple chains (for various concentrations) of PAA are simulated with simple point charge (SPC) water in a cubic simulation box with
TE D
periodic boundary conditions. Appropriate number of HI counterions are added, as needed, to maintain the overall electroneutrality of the simulation box. The simulation box size (6 1L 15
for all simulations) is significantly larger than the chain contour length (≈ 5 1L) in order to
simulations.
EP
avoid periodicity artefacts. The density of the simulation box is approximately 1 O/QQ for all
AC C
Simulations have been performed in three steps. In the first step, energy minimization is performed using the steepest descent technique, with an initial step size of 0.001 1L and 20
adaptive step size control, until the maximum force on each atom become smaller than 10 R/01L ⋅ LTU4. In the second step, MD equilibration is performed for ∼ 60 1W, beyond which the properties of interest converge around an average value. Finally, in the third step, MD production is performed for 5 1W with a typical sampling frequency of 0.5 W. Both equilibration and production simulations are performed in the NVT ensemble with a
8
ACCEPTED MANUSCRIPT reference temperature of 298 K. Nosé-Hoover thermostat is used for temperature control, with a coupling constant of 0.5 W in both equilibration and production simulations. Electrostatic calculations are performed using the particle mesh Ewald (PME) scheme with a Fourier spacing of 0.1 1L, van der Waals cut-off of 1 1L, and cubic interpolation. Bonds containing hydrogen atoms are constrained using a fourth-order LINCS (Linear Constraint
RI PT
5
Solver) algorithm.
SC
Number of oligomer-oligomer contacts with time, H 0X4 (normalized by the total number of atoms of oligomers) is monitored to track the equilibrium process. Here, a “contact” is
10
M AN U
counted if any atom of any oligomer is within a chosen threshold distance (0.6 1L) from any atom of another oligomer. Similar concept is used to obtain the fraction of “condensed” counterions (Y ), which is defined as the total number of “contacts” between the oxygen
atoms of groups and HI counterions, divided by the total number of oxygen atoms of
groups. In this case, a “contact” is counted if any oxygen atom of any group is
15
TE D
within a threshold distance of 0.22 1L from any H I counterion. Please note that the choice of the above threshold distances are somewhat arbitrary due to inherent ambiguities in the
EP
definition of a “contact” or a “condensed counterion”. The number of hydrogen bonds between the two species is defined as the number of donor-acceptor pairs that are within a
AC C
threshold distance 0.35 1L and the hydrogen-donor-acceptor angle less than 30 degrees.
Inter-molecular PAA-PAA radial distribution functions (RDFs), O0=4, are defined as the
20
normalized probability of finding any atom of a PAA chain at a distance between = and = + 3= from any atom of another PAA chain, where bin width, 3= = 0.002 1L, and
normalization is performed with such probability of an ideal gas. Radius of gyration ([& ) of
PAA is computed in two ways – average [& of individual chains and the [& of all chains. While the former represents the size of individual chains, the latter represents the size of PAA
9
ACCEPTED MANUSCRIPT
aggregates. Water content of PAA aggregates are quantified by computing the [& of all chains (entire PAA aggregate) and the solvent-accessible surface area[45] (SASA) computed using a 0.14 1L probe radius. A decrease in [& of all chains signifies compaction of
aggregates, i.e., an increase in the local density of aggregates (∝ H/[&] , H being the total number of atoms in PAA chains). Equivalently, a decrease in SASA/monomer (SASA
RI PT
5
divided by the total number of monomers in PAA chains) indicates a decrease in water content or an increase in the compaction (local density) of aggregates. Except H 0X4, all the
10
3.1. Single-Chain Simulations
M AN U
3. Results and Discussion
SC
above quantities are time-averaged over trajectories sampled during production simulations.
In order to validate the force-field employed in our simulations, we first compare the results
TE D
of single-chain simulations of atactic PAA solutions for different degrees of deprotonation () with that reported by Sulatha and Natarajan[46]. Unlike the other simulations reported in this work that uses the simulation protocol given in section 2, we use the simulation protocol used in their study[46] for these single-chain simulations. Since the actual tacticity patterns
EP
15
used in the simulations are not provided in their paper, we have performed the simulations for
AC C
three different tacticity patterns (Tables S1, S2, S3 in the Supporting Information). We observe that the [& averaged over the three chains increases with the increase in the degree of deprotonation (fig. S1 in the Supporting Information), though the magnitude of [&
20
significantly vary for the three chains. Further, the typical chain conformations at different values of show a transition from a collapsed to a stretched conformation with an increase in (fig. S2 in the Supporting Information). All these results compare well with the results
reported in their paper, but we do not attempt a direct comparison due to differences between
10
ACCEPTED MANUSCRIPT tacticity patterns in the two studies. It is worth noting that the force-field employed in our study has also been used in a previous molecular dynamics study of PAA-PDMAEMA complexation.[5]
5
RI PT
In order to investigate the effects of tacticity on the single-chain statistics in more detail, we then perform single chain simulations for PAA of different tacticities (isotactic, syndiotactic, and atactic) for the neutral ( = 04 and fully deprotonated ( = 1) cases. As shown in fig. 2 and Table 1, while the PAA chains are always more stretched for the fully deprotonated case
SC
than compared to the neutral case, [& for a given varies with tacticity in the order syndiotactic > atactic > isotactic. This difference between [& values for different tacticity is
explained using the fact that all groups are on the same side of PAA chain in the
M AN U
10
isotactic case, and adjacent groups are on the opposite sides of PAA chain in the
syndiotactic case. Therefore, the magnitude of short-ranged hydrophobic attraction for the
TE D
neutral case is larger for isotactic PAA than compared to syndiotactic PAA, resulting in more collapsed conformations. If a similar reasoning is applied to the case of fully deprotonated 15
chains, electrostatic repulsion should be higher in the isotactic case than compared to
EP
syndiotactic case, which would result in larger stretching in the isotactic case than compared to syndiotactic case. This is in contradiction with the observed result that syndiotactic PAA is
AC C
more stretched compared to isotactic PAA at = 1. It can be justified using the fact that the
presence of groups on the same side of PAA chain in the isotactic case results in the
20
bending of chain in order to keep the groups as far away as possible. Such bending is
not favored in the syndiotactic case as the groups are on both sides of the chain, resulting in a straight chain. Yet another explanation for larger stretching in the syndiotactic case compared to the isotactic case is the difference in the magnitudes of counterion condensation in the two cases. Since the fraction of condensed counterions in the isotactic 25
case is significantly higher than that in the syndiotactic case (Table 1), it would have 11
ACCEPTED MANUSCRIPT relatively lower electrostatic repulsion between chain segments due to the charge screening effect of counterions. The former reason appears more convincing, given the relatively large differences in the values of [& in the isotactic and syndiotactic cases. It is useful to compare
these [& values with their upper bound, that is, the theoretical estimate in the rodlike limit.
The − equilibrium bond distances along the backbone is ≈ 0.15 1L, which results in a theoretical estimate of [&,+
"
RI PT
5
= 39 /√12 ≈ 1.69 1L in the rodlike limit[47]. As
expected, all values of [& shown in Table 2 are lower than [&,+
"
and maximum [& is
SC
obtained for = 1 case of syndiotactic PAA that has most stretched configuration (fig. 2d) in
AC C
EP
TE D
M AN U
fig. 2.
10
Figure 2: Typical single chain conformations of PAA for neutral and fully deprotonated cases for three types of tacticity: (a)-(b) isotactic, (c)-(d) syndiotactic, (e)-(f) atactic PAA chain. Carbon, oxygen, and hydrogen atoms in PAA are shown in grey, red, and white colors, respectively. Water molecules and counterions are not shown. 15
12
ACCEPTED MANUSCRIPT
Table 1: Comparison of [& values of neutral ( = 0) and fully deprotonated states ( = 1) of single PAA chain of different tacticity. Fraction of “condensed” counterions (Y ) is also shown for the = 1 case. The values in brackets are the standard deviation of the last decimal. =1
[& (nm)
[& (nm)
Isotactic
0.83(4)
0.97(1)
0.28(6)
Syndiotactic
1.26(7)
1.42(2)
0.20(6)
Atactic
1.06(1)
1.16(1)
0.27(6)
SC
3.2. Multiple-Chain Simulations
Y
RI PT
Tacticity
5
=0
Four different concentrations of PAA are simulated by varying the number of PAA chains in
M AN U
a simulation box of 6 1L size containing SPC water. PAA concentration in the four systems, defined in moles per liter (M) of monomers are: (a) 0.615 M (4 PAA chains), (b) 1.223 M (8 PAA chains), (c) 2.456 M (16 PAA chains), and (d) 3.689 M (24 PAA chains). 10
To decide the amount of equilibration time required for the simulations, we monitor the
TE D
number of oligomer-oligomer contacts with time, which appear to converge after 60 1W, as established by simulations performed for a longer time of 90 1W in fig. 3 for the neutral
EP
syndiotactic case. Equilibration times of this order had also been reported in earlier atomistic simulation studies on polymer aggregation for systems of similar sizes.[33,36,48] Note that our estimate of equilibrium time is based on simulations performed over significantly long
AC C
15
time permitted within the current computational capabilities of atomistic simulations. The possibility of further aggregation or crystallization occurring over time scales that are much longer than that permitted by existing computational capabilities (e.g., time scales higher than 1 `W) can neither be established nor be ruled out. Interestingly, PAA aggregates span the 20
simulation box after equilibration in fig. 3, which at first glance, appear to be an artefact of periodic boundary condition used in our simulations. However, since the chain contour length is significantly smaller than the box size, PAA spanning the simulation box is less likely to be 13
ACCEPTED MANUSCRIPT a periodicity artifact, and is more likely to be reminiscent of a “gel”-like structure reported in
M AN U
SC
RI PT
a previous study.[33]
TE D
5
Figure 3: Number of oligomer-oligomer contacts with time for neutral syndiotactic PAA of concentration 3.689 M. Red dotted line indicate the convergence of oligomer-oligomer contacts around an average value. Simulation snapshots at the start of simulations (left) and after equilibration (right) are also shown. PAA chains are shown in red color; water molecules are not shown.
fully deprotonated cases. As expected, neutral cases for all tacticities formed more compact
AC C
10
EP
Next, we compare the aggregation behavior of PAA of different tacticity for the neutral and
aggregates than compared to the fully deprotonated cases, as evidenced by larger intermolecular RDF peaks for the neutral cases in fig. 4. However, the difference between the RDF for different tacticities at a given are not as substantial as the difference in [& of
single-chains for different tacticities (Table 1). This is because the groups of a PAA 15
chain in an aggregate interact not only with the groups on the same chain, but also
with groups on other PAA chains, which in turn may be oriented in any direction with respect to the PAA chain. Therefore, changes in the single-chain conformations due to
14
ACCEPTED MANUSCRIPT tacticity play little role in determining the aggregation behavior. It is worth pointing out that the first peak around = ∼ 0.2 1L for the neutral case in fig. 4 is due to hydrogen bonding between groups. Roughly, ∼ 0.25 − hydrogen bonds/monomer were observed for all tacticities. This peak is absent for the fully deprotonated cases, as the
groups do not hydrogen bond. The second peak occurring around = ∼ 1 1L represent an
RI PT
5
aggregation length scale of the order of [& (Table 1). This peak is much broader for the neutral case than compared to the fully deprotonated case, implying more PAA aggregation
Figure 4: Intermolecular PAA-PAA RDF for neutral ( = 0) and fully deprotonated ( = 1) PAA of different tacticities at a concentration of 3.689 M.
AC C
10
EP
TE D
M AN U
SC
in the neutral case.
Figure 5 shows the equilibrium configurations of syndiotactic PAA chains at various concentrations for = 0.4 with random deprotonation (Table S4 in the Supporting 15
Information). As expected, aggregation increases with increase in PAA concentration. Although the counterions have higher mobility than compared to PAA, they mostly reside within the PAA aggregate, thus neutralizing its overall charge. Further, when compared to = 0 case in fig. 3 (fig. 5d has the same concentration as fig. 3), the aggregate appears less 15
ACCEPTED MANUSCRIPT compact. The effect of PAA concentration on the aggregation behavior is quantified in terms of several measurable parameters in Table 2. The average [& of individual chains increases with increase in PAA concentration, implying an increase in the stretching of individual chains. On the other hand, the overall [& of all chains first increases (until 2.456 M in Table
2) and then decrease (for 3.689 M in Table 2). The initial increase in [& can be attributed to
RI PT
5
an increase in the size of PAA aggregates with increase in the number of PAA chains, and the final decrease (after [& ≈ 3 1L) can be understood to result from the compaction of
SC
aggregates. Although the [& values show a slightly non-monotonic trend with concentration,
increases uniformly with increase in concentration. This is further confirmed by the fact that the water content of aggregate, characterized by the SASA/monomer roughly decreases with increase in concentration, with small anomaly at the two intermediate concentrations where
EP
TE D
the SASA/monomer is nearly same.
AC C
10
M AN U
the local density of aggregates (∝ H/[&] , H being the total number of atoms in PAA chains)
16
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
TE D
5
Figure 5: Simulation snapshots (after equilibration) of syndiotactic PAA with = 0.4 for different concentration with random deprotonation (Table S4 in the Supporting Information). PAA chains and counterions are shown in red and green colors, respectively; water molecules are not shown.
AC C
EP
Table 2: Radius of gyration (for single chain and all chains), SASA/monomer, and fraction of condensed ions, of syndiotactic PAA with = 0.4 for different concentration with random deprotonation (Table S4 in the Supporting Information). The values in brackets are the standard deviation of the last decimal. PAA
[& (nm)
Concentration Individual
SASA/monomer All
2
(nm )
Y
chains
chains
0.615 M
1.05 (2)
2.0(4)
0.65(1)
0.26(5)
1.223 M
1.185(9)
2.8(2)
0.567(6)
0.30(4)
2.456 M
1.296(5)
3.06(6)
0.573(6)
0.32(3)
3.689 M
1.309(3)
2.68(3)
0.542(4)
0.32(2)
10
17
ACCEPTED MANUSCRIPT
The effect of on PAA aggregation behavior is more clearly demonstrated in fig. 6, which
shows that the compaction of PAA aggregate decreases with increase in . Electrostatic
repulsion between PAA chains increases with increase in resulting in the formation of less
compact structures. This is further confirmed by a monotonic increase in SASA/monomer with an increase in at the given concentration (Table 3). Such monotonic trend is not
RI PT
5
observed in the [& of individual chains and [& of all chains in Table 3. However, roughly
speaking, both the [& of individual chains [& of all chains increase with increase in ,
SC
signifying an increase in the stretching of individual chains and a decrease in overall
M AN U
compaction of aggregates, respectively.
AC C
EP
TE D
10
18
5
AC C
EP
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 6: Simulation snapshots (after equilibration) for syndiotactic PAA of concentration 3.689 M for different values of with random deprotonation (Table S4 in the Supporting Information). PAA chains and counterions are shown in red and green colors, respectively; water molecules are not shown.
19
ACCEPTED MANUSCRIPT
[& (nm)
Y
SASA/monomer 2
(nm )
All
chains
chains
0
1.272(5)
2.66(4)
0.375(4)
0.2
1.268(4)
2.76(5)
0.511(4)
0.4
1.309(3)
2.68(3)
0.542(4)
0.6
1.342(2)
2.70(1)
0.569(4)
0.39(2)
0.8
1.320(3)
2.91(4)
0.603(3)
0.43(2)
1.0
1.357(2)
2.99(1)
RI PT
Individual
-
0.19(2) 0.32(2)
SC
M AN U
5
Table 3: Radius of gyration (for single chain and all chains), SASA/monomer, and fraction of condensed ions, of syndiotactic PAA at different degrees of deprotonation with random deprotonation (Table S4 in the Supporting Information). The values in brackets are the standard deviation of the last decimal.
0.651(3)
0.45(2)
More information about the structure of PAA aggregates is revealed by the intermolecular PAA-PAA RDF (fig. 7). Second peak of this RDF occurring around ∼ 1 1L is higher and
10
TE D
broader for smaller values of indicating an increase in PAA aggregation, as already observed in fig. 6. However, the first peak occurring around ∼ 0.15 1L (inset of fig. 7) do
not exhibit the same trend due to rather non-monotonic changes in intermolecular hydrogen-
EP
bonding with . While cannot hydrogen bond with , it can hydrogen bond with , which contribute to the first peak along with − hydrogen bonding. In
15
AC C
the extreme case of fully deprotonated chain ( = 1 case in fig. 7), the first peak is missing as there is no possibility of hydrogen bonding between groups. For partially ionized cases
( = 0.2,0.4,0.6,0.8 case in fig. 7), we observe a peak at = ≈ 0.16 1L due to the combined
effect of hydrogen bonding between and groups, along with −
hydrogen bonding. In the other extreme of a neutral chain ( = 0 case in fig. 7), the first peak
occurs at a slightly larger value of = ≈ 0.17 1L as the groups are not present. The
20
hydrogen bonding data in Table 4 shows that the − hydrogen bonding 20
ACCEPTED MANUSCRIPT
decreases with due to a reduction in the number of groups. On the other hand,
− hydrogen bonding first increases and reaches a maximum at an intermediate value of , followed by a decrease at even higher value of . This is because the number of
RI PT
possible − pairs is maximum at an intermediate value.
M AN U
SC
5
Figure 7: Intermolecular PAA-PAA RDF for syndiotactic PAA of 3.689 M concentration for different values of with random deprotonation. O0=4 for small values of = are shown in the inset.
and
and
and
0.2 0.13(1)
0.01(1)
0.4 0.18(2)
and
2.456 M
and
0.18(1)
0.05(1)
0.03(1)
0.21(1)
0.6 0.17(1)
< 0.01
0.8 0.12(1)
< 0.01
0
1.223 M
EP
0.615 M
TE D
Table 4: Hydrogen bonding per monomer for syndiotactic PAA of different concentration for different f with random deprotonation. The values in brackets are the standard deviation of the last decimal. Values c 0.01 are not reported, since the standard deviation becomes comparable or higher than the mean.
AC C
10
−
0.05(3)
3.689 M
and
and
and
0.23(1)
0.04(1)
0.24(1)
0.05(1)
0.02(0)
0.22(1)
0.02(1)
0.25(1)
< 0.01
0.20(1)
< 0.01
0.19(0)
< 0.01
0.21(1)
< 0.01
0.07(1)
< 0.01
0.08(0)
< 0.01
0.11(0)
< 0.01
−
0.19(1)
−
0.25(2)
−
0.22(1)
21
ACCEPTED MANUSCRIPT Further insight into the intermolecular interactions can be obtained by a closer inspection of two nearby PAA chains in the concentrated system. Figures 8(a), 8(b), 8(c) show the zoomed view of simulation snapshots of the figs. 6(a), 6(c), and 6(f), for the neutral, partially deprotonated ( = 0.4), and fully deprotonated cases, respectively. As evident from fig. 8,
the distance between PAA chains at a given concentration increases with an increase in due
RI PT
5
to an increase in electrostatic repulsion (or a decrease in hydrophobic attraction) between PAA chains. However, as noted above, hydrogen bonding between and groups
SC
in the partially ionized case (fig. 8b) works together with − hydrogen bonding,
resulting in substantial attraction between PAA chains. Finally, while such hydrogen bonding is absent in the fully deprotonated case (fig. 8c), segments of PAA chain can still show an
M AN U
10
effective attraction due to “counterion bridging” as the HI counterions sandwiched between the two PAA chains experience electrostatic attraction from both the PAA chains. The likelihood of counterion bridging increases with increasing either the concentration or the
15
TE D
degree of deprotonation of PAA, as the fraction of condensed ions (Y ) increases both with increasing concentration (Table 2) and with increase in (Table 3). However, this
“counterion bridging” (or “salt-bridging”) effect is much weaker compared to hydrogen
EP
bonding interactions as the counterions have high entropy and prefers to disperse in the
AC C
solution. Therefore, PAA chains are mostly separated from each other in the fully ionized case (fig. 6f). 20
22
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
TE D
EP
5
Figure 8: Zoomed view of simulation snapshot for syndiotactic PAA of 3.689 M concentration for different values of with random deprotonation, showing the typical interactions between two nearby chains. Average distance between the center of mass of chains are 2.31 d 0.02 1L, 2.60 d 0.02 1L and 3.20 d 0.01 1L for = 0, 0.4, and 1, respectively. Carbon, oxygen, and hydrogen atoms in PAA are shown in grey, red, and white colors, respectively. Counterions are shown in green color. Water molecules are not shown.
10
AC C
Next, we compare the aggregation behavior of the partially deprotonated case for random deprotonation and end deprotonation at the same concentration. As shown in figs. 9(a) and 9(b), the overall aggregation behavior appear similar. This is further confirmed by the comparison of their intermolecular PAA-PAA RDF in fig. 9c. However, the dense regions formed by the hydrogen bonding of groups is more populated and prominent in the end deprotonation case (fig. 9b) compared to the random deprotonation case (fig. 9a). In fact, 15
for realistic chains (much longer than the oligomers studied here), a micelle-like formation (appearance of a dense core of hydrophobic segments) is expected in the end deprotonation 23
ACCEPTED MANUSCRIPT case due to the hydrogen bonding between the protonated ends of chains. On the other hand, a network-like formation with no such dense core is expected in the random deprotonation case of realistic chains, as the protonated groups are randomly distributed along the chains. Such differences will however be only apparent at intermediate values of (e.g., = 0.2 or
0.4). At lower (higher) values of , most of the chain is hydrophobic (hydrophilic) resulting
RI PT
5
in the formation of a dense precipitate (dilute solution). This can be substantiated by comparing the magnitudes of SASA/monomer for the random deprotonation case (Table 3)
SC
and end deprotonation case (Table 5). SASA/monomer magnitudes are significantly lower for the end deprotonation case than compared to the random deprotonation case for = 0.2 and
M AN U
0.4, which implies the formation of denser structures in the end deprotonation case. However, the system size and the length of chain is not sufficient enough to obtain a clear distinction of
EP
TE D
micelle-like and network-like structures.
AC C
10
24
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
10
EP
AC C
5
Figure 9: Simulation snapshot of syndiotactic PAA of 3.689 M concentration for (a) random deprotonation (Table S4 in the supporting information) and (b) end deprotonation (Table S5 in the supporting information), at fixed = 0.4, Water molecules and counterions are not shown here; aliphatic backbone carbon, group, and group are shown in red, blue, and green color, respectively. (c) Intermolecular PAA-PAA RDF for random (dashed line) and end deprotonation (bold line).
25
ACCEPTED MANUSCRIPT Table 5: Radius of gyration (for single chain and all chains), SASA/monomer (hydrophobic, hydrophilic, and total), and fraction of condensed ions, of syndiotactic PAA at different degrees of deprotonation with end deprotonation (Table S5 in the Supporting Information). The values in brackets are the standard deviation of the last decimal.
[& (nm)
Y
SASA/monomer 2
(nm )
All
chains
chains
0
1.272(5)
2.66(4)
0.375(4)
0.2
1.258(5)
2.83(7)
0.467(4)
0.4
1.271(4)
2.96(6)
0.522(4)
0.6
1.331(3)
2.79(4)
0.585(3)
0.8
1.339(3)
2.60(3)
0.624(4)
0.42(2)
1.0
1.357(2)
2.99(1)
0.651(3)
0.45(2)
RI PT
Individual
-
0.28(3) 0.38(2)
M AN U
SC
0.41(2)
5
The differences in phase behavior observed with changes in (fig. 6) and deprotonation pattern (fig. 8) has important implication in the use of polyelectrolytes as stimuli-responsive
10
TE D
materials for various applications. For example, for pH-responsive polyelectrolytes used as drug carriers, it is critical to attain substantial change in the drug release rate with change in pH. This can be realized due to an increase in the drug diffusivity within the carrier with a
EP
decrease in polyelectrolyte volume fraction, or equivalently, an increase in the water volume fraction. The water content of the PAA aggregates can be estimated by the solvent-accessible
15
AC C
surface area (SASA). In figure 10, we plot the SASA normalized by its value at = 0 (measure of change in water content with ) against . As shown in the plot, SASA changes
by ∼ 3 times in the entire range of , which should give rise to substantial change in the drug diffusivity within the aggregate with pH.
26
RI PT
ACCEPTED MANUSCRIPT
SC
M AN U
5
Figure 10: Solvent accessible surface area (SASA) divided by its value at = 0 against for random deprotonation and end deprotonation cases for syndiotactic PAA of 3.689 M concentration.
4. Conclusion
In summary, the results of atomistic MD simulations provide a detailed molecular picture of
10
TE D
the PAA aggregation behavior as a function of the degree of deprotonation or pH. PAA aggregation increases with PAA concentration and an increase in its degree of deprotonation (pH). Further, the nature of aggregation depends on the deprotonation pattern; end
EP
deprotonation and random deprotonation showed the formation of micelle-like and networklike aggregates, respectively. Interactions that were found relevant to the aggregation
15
AC C
behavior included the hydrogen bonding between − and −
groups, electrostatic repulsion between groups, and counterion (HI ) bridging
between the groups. Finally, the water content of aggregates expressed in terms of solvent-accessible surface area showed a 3-fold increase between the neutral and fully deprotonated limits.
20
Several important aspects studied in this paper require further investigation. First and foremost, our simulations are limited to model oligomers of PAA containing 20 monomers. 27
ACCEPTED MANUSCRIPT Studying the effect of chain flexibility on the aggregation behavior will require the simulations of much longer chains, well beyond current computational limitations of atomistic simulations. A useful compromise can be made using a systematic coarse-grained model parametrized using the atomistic simulation results described here, which may then be applied to study longer realistic chains. Second, although an increase in water content with
RI PT
5
pH signify a decrease in the drug diffusivity when PAA is used as a drug carrier, simulations using actual drug molecules are required to also include the PAA-drug interactions on the
SC
drug release. Third, we have confined our model to aqueous PAA solutions with no additional salt. Effects of the inclusion of salt species, particularly those present in the body fluids, are thus not captured. We are currently working along these lines of investigation.
M AN U
10
Acknowledgments
TE D
This work is supported by the INSPIRE award from the Department of Science and Technology (DST), India (no. IFA14/ENG-72) and a faculty initiation grant from IIT Roorkee obtained by P.K.J. References
M. Nič, J. Jirát, B. Košata, A. Jenkins, A. McNaught, eds., IUPAC Compendium of
AC C
[1]
EP
15
Chemical Terminology, IUPAC, Research Triagle Park, NC, 2009. doi:10.1351/goldbook.
20
[2]
A. V. Dobrynin, Theory and simulations of charged polymers: From solution properties to polymeric nanomaterials, Curr. Opin. Colloid Interface Sci. 13 (2008) 376–388. doi:10.1016/j.cocis.2008.03.006.
28
ACCEPTED MANUSCRIPT [3]
M.A. Cohen Stuart, B. Hofs, I.K. Voets, A. de Keizer, Assembly of polyelectrolytecontaining block copolymers in aqueous media, Curr. Opin. Colloid Interface Sci. 10 (2005) 30–36. doi:10.1016/j.cocis.2005.04.004. E. Kizilay, A.B. Kayitmazer, P.L. Dubin, Complexation and coacervation of
RI PT
[4]
polyelectrolytes with oppositely charged colloids, Adv. Colloid Interface Sci. 167
5
(2011) 24–37. doi:10.1016/j.cis.2011.06.006.
P.K. Jha, P.S. Desai, J. Li, R.G. Larson, pH and salt effects on the associative phase
SC
[5]
separation of oppositely charged polyelectrolytes, Polymers (Basel). 6 (2014) 1414–
10
[6]
M AN U
1436.
D. Priftis, M. Tirrell, Phase behaviour and complex coacervation of aqueous polypeptide solutions, Soft Matter. 8 (2012) 9396–9405. doi:10.1039/c2sm25604e. J. van der Gucht, E. Spruijt, M. Lemmers, M.A. Cohen Stuart, Polyelectrolyte
TE D
[7]
complexes: Bulk phases and colloidal systems, J. Colloid Interface Sci. 361 (2011) 407–422. doi:10.1016/j.jcis.2011.05.080. [8]
A. Kudlay, A. V. Ermoshkin, M.O. De La Cruz, Complexation of oppositely charged
EP
15
AC C
polyelectrolytes: Effect of ion pair formation, Macromolecules. 37 (2004) 9231–9241. doi:10.1021/ma048519t.
[9]
B.H. Morrow, G.F. Payne, J. Shen, pH-Responsive Self-Assembly of Polysaccharide through a Rugged Energy Landscape, J. Am. Chem. Soc. 137 (2015) 13024–13030.
20
doi:10.1021/jacs.5b07761.
29
ACCEPTED MANUSCRIPT [10] M.O. de la Cruz, L. Belloni, M. Delsanti, J.P. Dalbiez, O. Spalla, M. Drifford, Precipitation of highly-charged polyelectrolyte solutions in the presence of multivalent salts, J. Chem. Phys. 103 (1995) 5781–5791. doi:10.1063/1.470459.
5
RI PT
[11] A. Kudlay, M.O. De la Cruz, Precipitation of oppositely charged polyelectrolytes in salt solutions, J. Chem. Phys. 120 (2004) 404–412. doi:10.1063/1.1629271.
[12] A. Kundagrami, M. Muthukumar, Theory of competitive counterion adsorption on
SC
flexible polyelectrolytes: Divalent salts, J. Chem. Phys. 128 (2008) 244901.
M AN U
doi:10.1063/1.2940199.
[13] X. Guo, J. Fang, T. Watari, K. Tanaka, H. Kita, K.I. Okamoto, Novel sulfonated 10
polyimides as polyelectrolytes for fuel cell application. 2. Synthesis and proton conductivity of polyimides from 9,9-bis(4-aminophenyl)fluorene-2,7-disulfonic acid,
TE D
Macromolecules. 35 (2002) 6707–6713. doi:10.1021/ma020260w. [14] V.K. Sachan, A. Devi, R.S. Katiyar, R.K. Nagarale, P.K. Bhattacharya, Proton transport properties of sulphanilic acid tethered poly(methyl vinyl ether-alt-maleic anhydride)-PVA blend membranes, Eur. Polym. J. 56 (2014) 45–58.
EP
15
AC C
doi:10.1016/j.eurpolymj.2014.04.009. [15] C.S. Peyratout, L. Dahne, Tailor-made polyelectrolyte microcapsules: From multilayers to smart containers, Angew. Chemie - Int. Ed. 43 (2004) 3762–3783. doi:10.1002/anie.200300568.
20
[16] S.Y. Nie, W.J. Lin, N. Yao, X.D. Guo, L.J. Zhang, Drug Release from pH-Sensitive Polymeric Micelles with Different Drug Distributions: Insight from Coarse-Grained Simulations, ACS Appl. Mater. Interfaces. 6 (2014) 17668–17678. doi:10.1021/am503920m. 30
ACCEPTED MANUSCRIPT [17] D. Schmaljohann, Thermo- and pH-responsive polymers in drug delivery, Adv. Drug Deliv. Rev. 58 (2006) 1655–1670. doi:10.1016/j.addr.2006.09.020. [18] P. Coimbra, P. Ferreira, H.C. de Sousa, P. Batista, M.A. Rodrigues, I.J. Correia, M.H.
5
RI PT
Gil, Preparation and chemical and biological characterization of a pectin/chitosan polyelectrolyte complex scaffold for possible bone tissue engineering applications, Int. J. Biol. Macromol. 48 (2011) 112–118. doi:10.1016/j.ijbiomac.2010.10.006.
SC
[19] H. Almeida, M.H. Amaral, P. Lobão, Temperature and pH stimuli-responsive
polymers and their applications in controlled and selfregulated drug delivery, J. Appl.
10
M AN U
Pharm. Sci. 2 (2012) 01–10. doi:10.7324/JAPS.2012.2609.
[20] C. V. Hoven, A. Garcia, G.C. Bazan, T.Q. Nguyen, Recent applications of conjugated polyelectrolytes in optoelectronic devices, Adv. Mater. 20 (2008) 3793–3810.
TE D
doi:10.1002/adma.200800533.
[21] A. V. Dobrynin, M. Rubinstein, Theory of polyelectrolytes in solutions and at surfaces, Prog. Polym. Sci. 30 (2005) 1049–1118. doi:10.1016/j.progpolymsci.2005.07.006. [22] C. Holm, J.F. Joanny, K. Kremer, R.R. Netz, P. Reineker, C. Seidel, T.A. Vilgis, R.G.
EP
15
AC C
Winkler, Polyelectrolyte Theory, Adv. Polym. Sci. 166 (2004) 3–17. doi:10.1007/b11349.
[23] G.S. Manning, Limiting Laws and Counterion Condensation in Polyelectrolyte Solutions I. Colligative Properties, J. Chem. Phys. 51 (1969) 924.
20
doi:10.1063/1.1672157. [24] M. Muthukumar, Theory of counter-ion condensation on flexible polyelectrolytes: Adsorption mechanism, J. Chem. Phys. 120 (2004) 9343. doi:10.1063/1.1701839.
31
ACCEPTED MANUSCRIPT [25] Q. Wang, T. Taniguchi, G.H. Fredrickson, Self-consistent field theory of polyelectrolyte systems, J. Phys. Chem. B. 108 (2004) 6733–6744. doi:10.1021/jp037053y.
5
RI PT
[26] F. Alarcón, G. Pérez-Hernández, E. Pérez, A. Gama Goicochea, Coarse-grained simulations of the salt dependence of the radius of gyration of polyelectrolytes as
models for biomolecules in aqueous solution, Eur. Biophys. J. 42 (2013) 661–672.
SC
doi:10.1007/s00249-013-0915-z.
[27] S. Liu, M. Muthukumar, Langevin dynamics simulation of counterion distribution
10
doi:10.1063/1.1476930.
M AN U
around isolated flexible polyelectrolyte chains, J. Chem. Phys. 116 (2002) 9975.
[28] A. Laguecir, S. Ulrich, J. Labille, N. Fatin-Rouge, S. Stoll, J. Buffle, Size and pH effect on electrical and conformational behavior of poly(acrylic acid): Simulation and
TE D
experiment, Eur. Polym. J. 42 (2006) 1135–1144. doi:10.1016/j.eurpolymj.2005.11.023.
[29] G. Berghold, P. van der Schoot, C. Seidel, Equilibrium charge distribution on weak
EP
15
AC C
polyelectrolytes, J. Chem. Phys. 107 (1997) 8083–8088. doi:10.1063/1.475071. [30] R. Chockalingam, U. Natarajan, Self-association behaviour of atactic polymethacrylic acid in aqueous solution investigated by atomistic molecular dynamics simulations, Mol. Simul. 41 (2015) 1110–1121. doi:10.1080/08927022.2014.947481.
20
[31] S.H. Min, S.K. Kwak, B.-S. Kim, Atomistic simulation for coil-to-globule transition of poly(2-dimethylaminoethyl methacrylate), Soft Matter. 11 (2015) 2423–2433. doi:10.1039/C4SM02242D.
32
ACCEPTED MANUSCRIPT [32] Y. Liu, V.S. Bryantsev, M.S. Diallo, W.A. Goddard III, PAMAM Dendrimers Undergo pH Responsive Conformational Changes without Swelling, J. Am. Chem. Soc. 131 (2009) 2798–2799. doi:10.1021/ja8100227.
5
RI PT
[33] P.K. Jha, R.G. Larson, Assessing the efficiency of polymeric excipients by atomistic molecular dynamics simulations, Mol. Pharm. 11 (2014) 1676–1686. doi:10.1021/mp500068w.
SC
[34] D. Reith, B. Müller, F. Müller-Plathe, S. Wiegand, How does the chain extension of poly (acrylic acid) scale in aqueous solution? A combined study with light scattering
10
M AN U
and computer simulation, J. Chem. Phys. 116 (2002) 9100. doi:10.1063/1.1471901. [35] D.D. Hsu, W. Xia, S.G. Arturo, S. Keten, Systematic method for thermomechanically consistent coarse-graining: A universal model for methacrylate-based polymers, J.
TE D
Chem. Theory Comput. 10 (2014) 2514–2527. doi:10.1021/ct500080h. [36] W. Huang, R. Ramesh, P.K. Jha, R.G. Larson, A Systematic Coarse-Grained Model for Methylcellulose Polymers: Spontaneous Ring Formation at Elevated Temperature, Macromolecules. 49 (2016) 1490–1503. doi:10.1021/acs.macromol.5b02373.
EP
15
AC C
[37] P.W. Atkins, J. De Paula, Physical chemistry for the life sciences, W.H. Freeman and Co., New York, 2011.
[38] N. Hoda, R.G. Larson, Explicit- and Implicit-Solvent Molecular Dynamics Simulations of Complex Formation between Polycations and Polyanions,
20
Macromolecules. 42 (2009) 8851–8863. doi:10.1021/ma901632c.
33
ACCEPTED MANUSCRIPT [39] R. Chockalingam, U. Natarajan, Self-association behaviour of atactic polymethacrylic acid in aqueous solution investigated by atomistic molecular dynamics simulations, Mol. Simul. 41 (2015) 1110–1121. doi:10.1080/08927022.2014.947481.
5
RI PT
[40] S. Donnini, F. Tegeler, G. Groenhof, H. Grubmüller, Constant pH molecular dynamics in explicit solvent with λ-dynamics, J. Chem. Theory Comput. 7 (2011) 1962–1978. doi:10.1021/ct200061r.
SC
[41] U. Börjesson, P.H. Hünenberger, Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines, J. Chem. Phys. 114 (2001)
10
M AN U
9706–9719. doi:10.1063/1.1370959.
[42] S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M.R. Shirts, J.C. Smith, P.M. Kasson, D. Van Der Spoel, B. Hess, E. Lindahl, GROMACS 4.5: A highthroughput and highly parallel open source molecular simulation toolkit,
TE D
Bioinformatics. 29 (2013) 845–854. doi:10.1093/bioinformatics/btt055. [43] E. Lindahl, P. Bjelkmar, P. Larsson, M.A. Cuendet, B. Hess, Implementation of the charmm force field in GROMACS: Analysis of protein stability effects from
EP
15
correction maps, virtual interaction sites, and water models, J. Chem. Theory Comput.
AC C
6 (2010) 459–466. doi:10.1021/ct900549r.
[44] V. Zoete, M.A. Cuendet, A. Grosdidier, O. Michielin, SwissParam: A fast force field generation tool for small organic molecules, J. Comput. Chem. 32 (2011) 2359–2368.
20
doi:10.1002/jcc.21816.
34
ACCEPTED MANUSCRIPT [45] F. Eisenhaber, P. Lijnzaad, P. Argos, C. Sander, M. Scharf, The double cubic lattice method: Efficient approaches to numerical integration of surface area and volume and to dot surface contouring of molecular assemblies, J. Comput. Chem. 16 (1995) 273–
5
RI PT
284. doi:10.1002/jcc.540160303. [46] M.S. Sulatha, U. Natarajan, Origin of the difference in structural behavior of
poly(acrylic acid) and poly(methacrylic acid) in aqueous solution discerned by
SC
explicit-solvent explicit-ion MD simulations, Ind. Eng. Chem. Res. 50 (2011) 11785– 11796. doi:10.1021/ie2014845.
Sons, Inc., New York, USA, 2002. doi:10.1002/0471224510.fmatter_indsub. [48] W. Huang, I.S. Dalal, R.G. Larson, Analysis of solvation and gelation behavior of methylcellulose using atomistic molecular dynamics simulations, J. Phys. Chem. B.
EP
TE D
118 (2014) 13992–14008. doi:10.1021/jp509760x.
AC C
10
M AN U
[47] I. Teraoka, Polymer Solutions: An Introduction to Physical Properties, John Wiley &
35
ACCEPTED MANUSCRIPT
Highlights
EP
TE D
M AN U
SC
RI PT
Atomistic MD simulations of polyacrylic acid (PAA) have been performed. We study the effects of tacticity, concentration, and pH on the phase behavior. Aggregates form at high concentrations and their water content increase with pH. Results of this study may be useful for the design of pH-responsive drug carriers.
AC C
• • • •