Dynamics of Dilute Solutions of Poly(aspartic acid) - American

Oct 7, 2013 - Computational Simulations Group, Samsung R&D Institute India-Bangalore, Block-B, Bagmane Tridib, Bagmane Tech Park, C. V.. Raman Nagar, ...
0 downloads 12 Views 937KB Size
Article pubs.acs.org/JPCB

Dynamics of Dilute Solutions of Poly(aspartic acid) and Its Sodium Salt Elucidated from Atomistic Molecular Dynamics Simulations with Explicit Water Sanoop Ramachandran,*,† Anki Reddy Katha,† Subramanya Mayya Kolake,† Bokyung Jung,‡ and Sungsoo Han‡ †

Computational Simulations Group, Samsung R&D Institute India-Bangalore, Block-B, Bagmane Tridib, Bagmane Tech Park, C. V. Raman Nagar, Bangalore 560 093, India ‡ Organic Electronics Materials Lab, Samsung Advanced Institute of Technology (SAIT), Gyeonggi-do 446 712, Republic of Korea S Supporting Information *

ABSTRACT: The use of forward osmosis (FO) process for seawater desalination has attracted tremendous interest in recent years. Besides the manufacture of suitable membranes, the major technical challenge in the efficient deployment of the FO technology lies in the development of a suitable “draw solute”. Owing to its inherent advantages, poly(aspartic acid) has arisen to be an attractive candidate for this purpose. However, an investigation of its molecular level properties has not been studied in detail. In this paper, the dynamics of poly(aspartic acid) and its sodium salt in the dilute concentration regime have been reported. The quantification of the polymer conformational properties, its solvation behavior, and the counterion dynamics are studied. The neutral polymer shows a preferentially coiled structure whereas the fully ionized polymer has an extended structure. Upon comparing with poly(acrylic acid) polymer, another polymer which has been used as a draw solute, poly(aspartic acid) forms more number of hydrogen bonds as well as fewer ion pairs.



INTRODUCTION Polyelectrolytes are polymers with ionizable groups which dissociate in aqueous solution.1 The ions which are released into the solution are termed counterions. The resulting polymer backbone is charged and is responsible for the deviation of both physical and chemical properties from that of neutral polymers. The competition between the relatively shortranged solvation interactions and long-ranged electrostatic interactions can lead to phenomena such as counterion condensation. Such charged polymeric systems also play important roles in a variety of fields such as the chemical and pharmaceutical industry, material science, and biology.2,3 Of particular interest is in the field of seawater desalination via the exploitation of the phenomenon of osmosis. Osmosis is the flow of water across a semipermeable membrane from a region of higher chemical potential to a lower one. In the context of seawater desalination and water treatment, the process of reverse osmosis (RO) has been the most favored method till date.4 In this method, an external pressure is applied against the gradient in chemical potential of water so as to obtain purified water. In recent years, the use of direct osmosis or forward osmosis for the same purpose has gained acceptance.5,6 In this method, the normal osmosis process is exploited by using a “draw solution” with a higher osmotic © 2013 American Chemical Society

pressure than the feed (solution to be purified). During the process, the draw solution becomes diluted and pure water is subsequently separated out via a secondary filtration method. The FO process operates at substantially lower hydraulic pressures which is a major advantage. This leads to a considerable decrease in energy costs as well as a reduction in the membrane fouling (degradation owing to bacterial growth, accumulation of salts etc) observed during the standard RO process. Apart from the research into the design and manufacture of efficient semipermeable membranes, the key focus in this area is in the development of suitable draw solutes.5−7 There are several criteria which must be simultaneously satisfied by a candidate compound in order for it to be a suitable draw solute. Some of the most important considerations are the following: (a) generation of high osmotic pressure (at lower concentration of draw solute), (b) high solubility in water, (c) low viscosity and high diffusivity, (d) relatively low molecular weight, and (e) ambient operating temperature for both FO and draw solute recovery. Several Received: July 9, 2013 Revised: September 30, 2013 Published: October 7, 2013 13906

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B

Article

studies on draw solutes based on inorganic salts,8 organic ionic salts,9 and polyelectrolyte solutions from poly(acrylic acid) polymer10 have appeared recently. A more comprehensive survey of draw solutes can be found in the recent review by Chekli et al.7 Although individual experiments can be used to quantify all aspects described above, a microscopic understanding of the mechanisms will provide key insights for the selection criteria. Despite the need for several aspects to be simultaneously satisfied, in this article, we attempt to investigate the factors that may be responsible for the generation of high osmotic pressure. Theoretical approaches based on electrostatic theory are usually very successful in explaining the global properties of charged polymers. However, they are found wanting in the atomistic length and time scale regime.2,3,11,12 Theories generally address universal behavior and focus on low ionic strength solutions.13−15 For more complicated systems with high ionic strength and the presence of multivalent ions, the properties are not well understood. Most of the computational works on polyelectrolyte systems also focus on generic behavior via coarse-grained models.16−21 One of the aspects that are neglected in theoretical works is the phenomenon of hydrogen bonding. This can be ascertained easily via atomic level simulations. Atomistic simulations (although with their own set of challenges) can therefore provide much needed insights. Recently, there have been a number of computational studies using atomistic models on polyelectrolytes such as poly(acrylic acid),22,23 and poly(styrenesulfonate).24−26 However, to the best of our knowledge, there have not been any notable study on poly(aspartic acid) (PAsp) systems till date. PAsp-Na is the sodium salt of poly(aspartic acid) and is a synthetic polypeptide which can be manufactured on an industrial scale. It is biodegradable27 and poses no notable environmental risks, an important consideration if the FO process is used to generate potable water. The present study explores the characteristics and dynamics of a single-chain poly(aspartic acid) and its sodium salt to gain a microscopic understanding using molecular dynamics (MD) simulations with explicit water. Further, a comparison is also made between the sodium salts of poly(acrylic acid) and poly(aspartic acid) polymers. This will then set the foundation for future studies on concentrated systems and help in gaining further insights into the development of more optimum draw solutes for the FO process.

Figure 1. (a) PAsp and (b) PAsp-Na repeat units.

coarse-grained into one effective carbon atom) and the polar hydrogens are treated explicitly. The water molecules are also treated explicitly using the SPC/E (extended simple point charge) water model29 whose bond lengths are held constant using the SETTLE algorithm.30 All the simulations were performed using the GROMACS simulation engine (version 4.5.5).31 The simulation box is cubic with periodic boundary conditions applied in all the three independent directions. The van der Waals interactions represented by the Lennard-Jones potential and were truncated and shifted at rc = 1.2 nm. The electrostatic interactions were calculated using the particle mesh Ewald technique (PME).32 The trajectories were calculated using the leapfrog integrator with an integration time step of 2 fs. The covalent bonds in the polymer chain were constrained using the LINCS algorithm.33 Poly(acrylic acid) has been previously simulated using the same force field, and we adopt those bonded and nonbonded parameters in our simulations.23 In the GROMOS 53a6 force field, the intramolecular bonded interactions are given by the following set of potentials 1 Vb(rij) = kijb(rij2 − b0 2)2 (1) 4



Va(θijk) =

ATOMISTIC MODEL AND SIMULATION METHODOLOGY The chemical structures of the polyelectrolytes, namely PAsp and PAsp-Na, are shown in Figure 1. In the simulations, they are represented by a single polymer strand in the syndiotactic configuration with degrees of polymerization N = 10 and N = 30 (in individual monomer units). The syndiotactic configuration was found to have the lowest energy verified via ab initio calculations using short polymer strands. As a measure of ionization of the polyelectrolytes we define, f = Nc/N, where Nc is the count of the number of charged monomers. In this study, we have considered f = 0 (neutral) and f = 1 (fully ionized). The ends of the polyelectrolyte chains were capped with NH2 and hydroxyl groups. For our simulations, we make use of the GROMOS 53a6 force field.28 In this force field the carbon atoms in the aliphatic chain are treated as united atoms (the aliphatic CHn atoms are

1 θ kijk[(cos(θijk) − cos(θijk0 )]2 2

Vd(ϕijkl) = kϕ(1 + cos(nϕ − ϕs))

(2) (3)

1 kξ(ξijkl − ξ0)2 (4) 2 Here, Vb is the bond interaction between atoms i and j separated by a distance rij with equilibrium bond distance b0 and spring constant kbij for the particular pair of bonded atoms. The angular interactions between two bonds connecting atoms, labeled by the indices i, j, and k with angle θijk are given by Va. θ0ijk is the equilibrium angle for the same set of bonds and kθijk is the force constant. The dihedral interaction energy due to the dihedral angle ϕijkl between the planes containing the atoms with indices i,j,k and j,k,l is given by Vd. ϕs is the equilibrium angle, n the dihedral multiplicity, and kϕ the force constant. Finally, Vid forms the improper dihedral interaction needed to prevent out-of-plane bending of atoms. The angle between the Vid(ξijkl) =

13907

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B



RESULTS AND DISCUSSION In this section, we discuss the polymer chain configurations, hydrogen bond formation, hydration behavior, and the counterion dynamics. For the N = 10 systems, the polymer weight percentage is ≈0.8%, and for N = 30 it is ≈0.2%. This puts the systems under study in the very dilute solution regime. In order to get an estimate of the relevant length scales involved, we appeal to mean-field electrostatic theory based on the linearized Poisson−Boltzmann equation, although, considering the atomistic length scale and local nonuniformity of the system, the validity of mean-field values can be questioned. However, at the very least it can be used to obtain an initial estimate. The distance over which the electrostatic energy becomes comparable to the thermal energy is defined as the Bjerrum length, lB = e2/(4πε0εrkBT), where e is the elementary charge (1.6 × 10−19 C), ε0 is the vacuum permittivity (8.85 × 10−12 C2 N−1 m2), εr is the dielectric constant of the medium (for SPC/E water a range of values between 68 and 81 have been reported36,37), kB is the Boltzmann constant (1.38 × 10−23 J K−1), and T is the temperature (300 K). Thus, we can approximate 0.7 < lB < 0.8 nm. Another relevant length scale is the Debye−Hü c kel screening length given by, l DH = (8πlBcs)−1/2, where cs is the concentration of additional salt added to the system. Since we have a salt-free system, lDH ≫ lB and is of the order of micrometers for pure water. This indicates that the electrostatic interactions are important at the length scales of the system studied and should be properly taken into account. Chain Conformations. As mentioned in the Introduction, the presence of charged groups can alter the equilibrium conformation of the polymers. In Figures 2 and 3, the

planes containing atoms with indices i,j,k and j,k,l is given by ξijkl, the equilibrium angle is given by ξ0 and force constant kξ. There are two types of nonbonded interactions used, namely van der Waals interactions represented by the Lennard-Jones (LJ) potential and Coulombic interactions. The LJ interaction potential is given by VLJ(rij) =

Cij(12) rij12



Cij(6) rij6

(5)

The parameters C(12) and C(6) ij ij depend on pairs of atom types with the following combination rule employed for interaction between different atom types 1/2 Cij(6) = (Cii(6)C(6) jj )

(6)

1/2 Cij(12) = (Cii(12)C(12) jj )

(7)

The Coulombic interaction is given by 1 qiqj VC(rij) = 4πε0 εrrij

Article

(8)

Here, qi and qj are the charges on particles labeled by i and j separated by a distance rij, and ε0 and εr are the vacuum and relative dielectric permittivities, respectively. For completeness, the force-field parameters used in the simulations are provided as Supporting Information. For the initial configuration of the polymers, we used an extended syndiotactic structure centered in the simulation box with the minimum distance from the boundary set at rc = 1.2 nm. The water molecules and the Na+ counterions (for the f = 1 case) were distributed randomly in the simulation box. The un-ionized, f = 0, case was solvated with 7557 water molecules in a simulation box of size 6.1 × 6.1 × 6.1 nm3 for the N = 10. For N = 30, 79 756 water molecules in a 13.4 × 13.4 × 13.4 nm3 box were used. For f = 1, to maintain charge neutrality, an equal number of water molecules were replaced by Na+ counterions. The energy of the solvated system was first minimized using the steepest descent method without any constraints. This was followed by a short 100 ps simulation under isothermic− isochoric (NVT) condition at temperature T = 300 K. Subsequently, a 1 ns simulation under isothermal−isobaric (NPT) conditions at pressure, p = 1 atm and at temperature T = 300 K was carried out. Both these simulations involved the application of the position restraints on the polymer chains so as to evenly thermalize the system and also allow the water molecules to equilibrate properly. The final equilibration step involved a 10 ns NVT simulation without any position restraints on the polypeptide chains. The production runs involved 50 and 30 ns NVT simulations for the chain lengths 10 and 30 monomers, respectively. The trajectories and energies are stored every 2 ps. For chain length 10, the last 20 ns were used for the calculations of the structural quantities, and for chain length 30, the last 10 ns were used. The physical quantities were then averaged over four simulation trials with independent initial configurations. For comparison, we also simulated the poly(aspartic acid) systems using the all-atom CHARMM-22 force field.34,35 The results obtained were qualitatively similar to those from the united atom GROMOS force field employed in the present study and hence the choice.

Figure 2. Polymer conformation for N = 10 with (a) PAsp and (b) PAsp-Na. The aliphatic carbon chain backbone is shaded cyan, the nitrogen atoms are blue, the oxygen atoms are red, and the hydrogen atoms are shown in white. Sodium ions are shaded yellow in color. Please note that only those ions within close proximity of the polymer are shown. The whole system is strictly electroneutral.

representative conformations for chain length N = 10 and N = 30 respectively are shown for both PAsp and PAsp-Na. The structures were visualized using the Visual Molecular Dynamics software.38 For the shorter chains (N = 10), the polymer conformations do not show any signs of collapse (even for the neutral chain) and the polymer configurations appear similar. The N = 30 snapshots depicted in Figure 3 highlight the salient differences. The neutral chain, shown in Figure 3a, has a collapsed or globular structure. This conformation is adopted in order to minimize the hydrophobic interactions between the aliphatic backbone and the water molecules. In Figure 3b, the conformation of PAsp-Na is shown. As expected, the competition between the electrostatic interactions between the charged groups, hydrophobic interaction between the 13908

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B

Article

Figure 4. (a) Normalized probability distribution of the radius of gyration, P(Rg), as a function of Rg for chain length N = 10 and (b) for N = 30. (c) Normalized probability distribution of the end-to-end distance, P(R), as a function of R for chain length N = 10 and (d) for N = 30. The solid and dashed curves indicate PAsp and PAsp-Na, respectively.

Figure 3. Polymer conformation for N = 30 with (a) PAsp and (b) PAsp-Na. The aliphatic carbon chain backbone is shaded cyan, the nitrogen atoms are blue, the oxygen atoms are red, and the hydrogen atoms are shown in white. Sodium ions are shaded yellow in color. Please note that only those ions within close proximity of the polymer are shown. The whole system is strictly electroneutral.

The probability distribution of the end-to-end distance is plotted in Figure 4, c and d, for N = 10 and N = 30, respectively. The neutral polymer, the N = 10 case, shows a major peak at ≈1 nm and a second one at 2.6 nm. This is consistent with the twin peaks seen in the plot of the distribution of the radius of gyration. For PAsp-Na, the major peak is observed at ≈2 nm. However, it is now broadened out, implying that the polymer explores a much broader range of R values between 1 and 3 nm. Proceeding to N = 30, the neutral chain also has a peak value of ≈1.1 nm and a subsequent second peak at 2 nm. It is indicative of a coiled structure. Similar to the shorter chain case, PAsp-Na shows a broader distribution with a peak at R ≈ 7 nm. From the data so far, it is clear that the neutral chain has a preferentially collapsed structure while the sodium salt case has an extended structure. There are limited experimental studies on poly(aspartic acid) systems. However, it has been reported that the PAsp adopts a coiled structure under normal physiological conditions.39 The observations from the simulations described so far are in qualitative agreement with the experiments. The averaged values are summarized in Table 1. Also provided are the averaged eigenvalues of the diagonalized radius of gyration tensor. The relative magnitudes of the gyration tensor gives a measure of the shape of the polymer. For a spherical shape all the three eigenvalues will be similar, as seen for the neutral polymers.

polymer backbone, and thermal energy leads to an extended conformation. In order to better quantify the equilibrium structural conformations of the polymers, we have calculated the radius of gyration, Rg, as well as the end-to-end distance, R. Experimentally, the radius of gyration can be determined via dynamic light scattering techniques. In the context of the simulations it is defined as N

Rg2 =

tot 1 ⟨∑ |ri − rCM|2 ⟩ Ntot i = 1

(9)

where rCM is the position vector of the center of mass of the polymer chain and i stands for the atom index with position vector ri. Please note Ntot is the total number of atoms in the polymer and should not be confused with N, which represents the monomer or polymerization index. In Figure 4, a and b, we plot probability distribution of the radius of gyration, P(Rg) for N = 10 and N = 30, respectively. The solid and dashed curves represent PAsp and PAsp-Na, respectively. Although the simulation snapshots of the polymers for N = 10 appeared to look similar, the probability distributions highlight the differences. The smallest Rg is for the neutral case whose distribution is bimodal with the largest peak value at ≈0.6 nm and the smaller peak at 0.9 nm. PAsp-Na has a higher value at ≈0.81 nm. For the longer chains, N = 30, the distinction between the neutral and sodium salt of the polymers is clearer. The solid line is the neutral case and has a sharp peak at Rg at 0.9 nm. The dashed line is for PAsp-Na showing a broad distribution with a peak value ≈2.3 nm. Another related quantity which can be readily determined from simulations is the end-to-end distance. The end-to-end distance, R, is calculated as the distance between the center of masses of first and last monomer unit of the polymer chain and is given by

R = ⟨|r1 − rN |⟩

Table 1. Average Values of the End-to-End Distance R (nm), the Radius of Gyration Rg (nm) and the Three Eigenvalues of the Radius of Gyration Tensora N

polymer

10

a b a b

30

(10)

a

13909

⟨R⟩ 1.18 2.0 1.28 6.63

± ± ± ±

0.68 0.39 0.47 0.89

⟨Rg⟩

⟨Rxx g ⟩

⟨Ryy g⟩

⟨Rzz g⟩

± ± ± ±

0.57 0.76 0.68 2.16

0.30 0.29 0.48 0.59

0.20 0.19 0.34 0.31

0.69 0.84 0.91 2.27

0.13 0.07 0.07 0.09

The polymers PAsp and PAsp-Na are labeled as a and b, respectively. dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B

Article

The structural dynamics of the polymers can be probed by measuring the appropriate time-correlation functions as well as time-series measurements, although owing to the relatively short as well as linear chain structures of the polymers we expect simple exponential decay. This has been verified by the calculations of time correlation of both the radius of gyration and the end-to-end distance. By fitting the curves to an exponential function and calculating the integral of the exponential fit, an approximate correlation time can be calculated. For example, for PASP-Na with N = 10, the correlation time of the radius of gryation is found to be close to 90 ps. Hydration Behavior. The PAsp molecule is a linear molecule with short side chains. Therefore, a simple configuration of the solvation shell around molecule is expected, with chain ends attracting more water molecules than the central region. The pair correlation function, g(r), between the two carboxylate oxygens of the polymer and the oxygen of water as a function of distance r, is plotted in Figure 5a for N = 30.

Figure 6. Spatial density isosurface of water molecules within 0.35 nm of the carboxylate oxygens of the first residue of PAsp-Na with N = 10.

figure illustrates that the water molecules are almost spherically distributed around the carboxylate groups. Hydrogen Bond Formation. The formation of hydrogen bonds (H-bonds) play a very significant role in manipulating the dynamics of biopolymers and charged polymers in aqueous solutions. On the atomic scales, both inter- as well as intrapolymer H-bond formation can take place via the carboxylate groups. A priori, the definition of H-bonds in the context of this work should be made clear. There are several conventions followed in prescribing a geometric criteria for determining the formation of H-bonds.42−44 A H-bond is defined to be formed if the distance between the donor and acceptor is less than 0.35 nm and the angle between hydrogen donor−acceptor is less than 30°.45 This suitably reproduces the results already published in the literature.46,47 A geometric criteria was followed over an energetic one as the framework of classical dynamics does not fully capture the hydrogen bond dynamics. The important criterion with respect to this paper is the relative differences in the H-bond formation between the different polymers. In the system under consideration, polymer−polymer, polymer−water, and water−water Hbonds can be formed. However, we will focus on the H-bond formation only involving the polymer and water. In Figure 7, we plot the distribution of the number of polymer−solvent H-bonds, P(Nhb). For N = 10, there is an

Figure 5. Pair correlation function between the two carboxylate oxygens of the polymer and (a) the oxygen of water and (b) both the hydrogens of water for N = 30. The solid and dashed curves indicate PAsp and PAsp-Na, respectively.

For PAsp, seen from the solid curve, the g(r) shows a value less than unity for r < 2 nm. This indicates a preferential expulsion of water molecules from regions close to the polypeptide residue. For PAsp-Na the carboxylate group is fully ionized. The first solvation shell, whose radius is defined by first minimum of the g(r), is found to be at r = 0.35 nm. The sharp peaks below r = 0.35 nm (dashed curves) indicate the possible formation of hydrogen bonds between the fully ionized carboxylate groups and water molecules. This will be further clarified in the next section. It can be seen that, for every COO− residue in the polymer, approximately 2.25 water molecules are present. The globular structure of PAsp, as seen by the low Rg and R values, can be thus be rationalized by the hydrophobic nature the polypeptide backbone. The pair correlation function, g(r), between the two carboxylate oxygens of the polymer and the two hydrogens of water is shown in Figure 5b for the N = 30 case. The trends are similar to Figure 5a with a linear shift, owing to the closer distance of approach of the hydrogen atoms to the carboxylate oxygens. For more insight into the local solvation structures of water around the PAsp-Na polymer, the spatial distribution functions (SDF) are computed (for example, see ref 40). In contrast to the average numbers provided by the radial distribution functions, the SDF provides the probability of finding an atom in the three-dimensional space around the central molecule.41 In Figure 6, the SDF is shown for water molecules at a distance of less than 0.35 nm from the carboxylate oxygens for PAsp-Na with chain length N = 10. The

Figure 7. Distribution of the number of hydrogen bonds between the polymer and water for chain length (a) N = 10 and (b) N = 30 (the solid and patterned shades represent PAsp and PAsp-Na, respectively).

average of 31 H-bonds for PAsp and it is 56 for PAsp-Na. For N = 30 the average total number of H-bonds formed by PAsp is 76 while PAsp-Na forms approximately 165. It is thus observed that the ionized polymer forms about twice the number of the H-bonds as that of the neutral polymer. 13910

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B

Article

of the measured diffusion coefficient appears to match that of polymers with similar molecular weight. However, due to the lack of experimental data, a more quantitative comparison could not be made. It can also be seen that the diffusion coefficient of the neutral chain is higher than that of the fully charged case. This can be rationalized by considering the more compact nature of the neutral chains (corroborated by the Rg data) as compared to the fully ionized sodium salt version. Distribution of Counterions. In Figure 8a, the radial distribution of the number of sodium ions around the

The hydrogen bond lifetime can be calculated using the theory of Luzar and Chandler.48−50 The H-bond correlation function is defined as c(t ) =

⟨h(0)h(t )⟩ ⟨h⟩

(11)

where the variable h(t) is unity if a tracked hydrogen bond exists and zero otherwise. The derivative −

dc(t ) ≡ k(t ) dt

(12)

describes the decay rate of the H-bonds from which the average bond lifetime, τHB = 1/k, can be calculated. By storing the configurations every 2 ps for a total time of 5 ns, we obtained the average lifetime of polymer−water H-bond to be as shown in Table 2. Table 2. Hydrogen Bond Lifetime for PAsp-Na for Two Chain Lengths N

lifetime (ps)

10 30

1.2 1.1

Figure 8. (a) Radial distribution function between the carboxylate oxygens and the sodium ion for N = 10. (b) Coordination number of the two carboxylate oxygens and sodium ions. The solid and dashed curves indicate PAA-Na and PAsp-Na, respectively.

carboxylate oxygens is shown for N = 10. There is a sharp first peak at r = 0.25 nm followed by a second one at r = 0.5 nm. Although the absolute values of the peaks are very high, the number of sodium ions in the near proximity are small. This is more clearly represented by the coordination number, cN(r) (equal to ∫ r04πx2g(x)ρ(x) dx with ρ(x) being the density in the shell of thickness dx at a radial distance x), of the two carboxylate oxygens of the polymer and the Na+ ions in Figure 8b for N = 10 case. An ion pair is defined to be formed if the atoms are closer than 0.35 nm (taken from the first minima found in the radial distribution function). Every possible ion pair that can be formed between carboxylate oxygen and sodium ions can be uniquely labeled and, subsequently, tracked as a function of time over the course of the simulations. The lifetime of the ion pairs can then be obtained by averaging over all the ion pairs formed. In Table 4, the average number of ion pairs ⟨Nip⟩, their

Due to the more favorable interactions with water, the PAspNa polymer forms a larger number of H-bonds with water molecules. This also contributes to the relatively open conformations of the polymer. Diffusion Coefficient. By tracking the center of mass position of the polymer as a function of time, an estimate of the polymer diffusion coefficient can be made via the calculation of the mean-squared displacement 1 D = lim ⟨ΔrCM 2⟩ (13) t →∞ 6t in three dimensions. Here ΔrCM2 ≡ |rCM(t) − rCM(0)|2 is the square of the displacement of the center of mass, D is the diffusion coefficient in three dimensions and t is the time interval. In order to calculate the diffusion coeffient in a more efficient way, a modified version of the order-n algorithm was used.51,52 This algorithm allows the simultaneous sampling of short and long time correlations at minimal computational cost (the conventional algorithm scales as t2 while the modified one scales as t). Briefly, the algorithm makes use of multiple time interval blocks, each with a different sampling frequency. If the sampling frequency between n blocks differ by a multiple of 10, then the block n samples every 10n time steps. The whole mean-squared displacement plot can then be reconstructed by stacking the accumulated blocks side by side. Further details as well as the computer code can be found in ref 52. The results of the mean-squared displacement of the polymer center of mass are tabulated in Table 3. The order of magnitude

Table 4. Comparison of Average Number (⟨Nip⟩), Distance (⟨dip⟩), and Lifetime (⟨tip⟩) Ion Pairs Formed for PAsp-Na and PAA-Na

N

PAsp

PAsp-Na

258.2 ± 13.3 237.7 ± 12.3

202.6 ± 14.2 97.5 ± 12.0

⟨Nip⟩

⟨dip⟩ (nm)

⟨tip⟩ (ps)

0.9 1.8

0.7 1.6

15.7 29.0

average distance ⟨dip⟩, and lifetime ⟨tip⟩ calculated from the simulations are listed. PAsp-Na forms less than one ion pair as well as the ions are relatively far away from the polymer backbone (0.7 nm). This indicates the negligible formation of ion pairs and also that the carboxylate groups are fully ionized. The ion-pair lifetime can be better appreciated when compared with poly(acrylic acid) sodium salt, which will be presented in the next section. Owing to the relatively fast dissociation of the ionic groups, the sodium ions are free to explore the full volume of the system leading to a large favorable entropic contribution to the system energy. This can lead to a positive contribution to the solution osmotic pressure.

Table 3. Diffusion Coefficient in Units of 10−12 m2 s−1, Measured from the Mean-Squared Displacement of the Polymer Center of Mass for PAsp and PAsp-Na with Chain Lengths N = 10 and 30

10 30

polymer PAsp-Na PAA-Na

13911

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B

Article

Comparison with Poly(acrylic acid) Sodium Salt. A recent study involved the use of the poly(acrylic acid) sodium salt (PAA-Na) as a draw solute for the FO process.10 It forms the obvious benchmark to which comparisons must be made. We now make comparisons between PAsp-Na and PAA-Na polymers. In Figure 9a, the radiuses of gyration of PAsp-Na and PAANa for chain length N = 10 are compared. The plot highlights a

contributions to the osmotic pressure, namely the entropic contribution due to free counterions in solution and a polymeric contribution.2,53 When compared to PAsp-Na, due to the formation of ion pairs (which implies attractive polymer−counterion interactions), the average number of free sodium ions in solution is reduced and thus could result in a relatively lower osmotic pressure of PAA-Na solutions. However, at this stage, the argument remains to be verified through more quantitative studies. Investigations on the concentrated regime of poly(aspartic acid) are ongoing and the results will be reported elsewhere. However, based on the results of the present study, the PAspNa shows promise to be a good candidate for an efficient draw solute in the forward osmosis process.



ASSOCIATED CONTENT

S Supporting Information *

The full parameter set for the force field used in the simulations. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 9. Comparison between PAsp-Na and PAA-Na for chain length N = 10: (a) the distribution of the radius of gyration; (b) the distribution of the number of hydrogen bonds between polymer− water. The dashed curve indicates PAsp-Na and the solid curve denotes PAA-Na.



AUTHOR INFORMATION

Corresponding Author

*Phone: +91-80-4181-9999. E-mail: [email protected].

clear distinction between the two polymer species. The PAANa shows a narrow peak about 0.7 nm whereas PAsp-Na has a very broad distribution. This indicates that PAsp-Na has a very dynamic behavior with constant fluctuations between collapsed and extended structures. In Figure 9b, we compare the total number of hydrogen bonds formed between the polymer and water for the N = 10 case, based on the previously mentioned distance and angle criteria. From the plot, it can be observed that PAsp-Na forms a slightly larger number of of H-bonds as compared to the PAANa. Quantitatively, the average number of H-bonds for PAA-Na is 48 while it is 56 for PAsp-Na. Also, a H-bond correlation analysis resulted in the average H-bond lifetime to be 1.42 ps which is approximately 15% higher than the same for PAsp-Na. In Table 4 we also compare the lifetimes of the ion-pair formation between the carboxylate oxygens and sodium ions. We observe that the average lifetime of the ion pairs for PAANa is double that of PAsp-Na.

Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS Thanks to Shanthi Pandian and Srinivasa Raju G. for helpful discussions. REFERENCES

(1) Oosawa, F. Polyelectrolytes; Marcel Dekker: New York, 1971. (2) Dobrynin, A. V. Polymer Science: A Comprehensive Reference; Elsevier BV: Amsterdam, 2012; Vol. 1, pp 81−132. (3) Dobrynin, A. V. Theory and Simulations of Charged Polymers: From Solution Properties to Polymeric Nanomaterials. Curr. Opin. Colloid Interface Sci. 2008, 13, 376−388. (4) Sourirajan, S. Reverse Osmosis; Academic Press Inc.: New York, 1970. (5) Cath, T. Y.; Childress, A. E.; Elimelech, M. Forward Osmosis: Principles, Applications and Recent Developments. J. Membr. Sci. 2006, 281, 70−87. (6) Zhao, S.; Zou, L.; Tang, C. Y.; Mulcahy, D. Recent Developments in Forward Osmosis: Opportunities and Challenges. J. Membr. Sci. 2012, 396, 1−21. (7) Chekli, L.; Phuntsho, S.; Shon, H. K.; Vigneswaran, S.; Kandasamy, J.; Chanan, A. A Review of Draw Solutes in Forward Osmosis Process and Their Use in Modern Applications. Desalin. Water Treat. 2012, 43, 167−184. (8) Achilli, A.; Cath, T. Y.; Childress, A. E. Selection of Inorganicbased Draw Solutions for Forward Osmosis Applications. J. Membr. Sci. 2010, 364, 233−241. (9) Bowden, K. S.; Achilli, A.; Childress, A. E. Organic Ionic Salt Draw Solutions for Osmotic Membrane Bioreactors. Bioresour. Technol. 2012, 122, 207−216. (10) Ge, Q.; Su, J.; Amy, G. L.; Chung, T. Exploration of Polyelectrolytes as Draw Solutes in Forward Osmosis Processes. Water Res. 2012, 46, 1318−1326. (11) Grosberg, A. Y.; N., T. T.; Shklovskii, B. I. The Physics of Charge Inversion in Chemical and Biological systems. Rev. Mod. Phy. 2002, 74, 329−345. (12) Barrat, J. L.; Joanny, J. F. In Advances in Chemical Physics: Polymeric Systems; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons, Inc.: Hoboken, NJ, 2007; pp 1−66.



CONCLUSIONS In this paper, we considered conformational and solvation behavior of poly(aspartic) acid and its sodium salt using molecular dynamics simulations with explicit solvent. We focused on the neutral poly(aspartic acid) polymer ( f = 0) and the fully ionized version of its sodium salt ( f = 1). The neutral polymers are found to favor a coiled or globular conformation whereas the fully ionized polymers prefer extended conformations. It can be concluded that water is a good solvent for the f = 1 case and a bad solvent for the neutral case. Upon comparison with PAA-Na, the significant difference appears to arise in the conformational state of the polymers with PAsp-Na showing a broad distribution of Rg values indicating a dynamic behavior. Further, PAA-Na forms fewer H-bonds but with a marginally higher lifetime as compared to PAsp-Na. The probability to form larger number of ion pairs (which are also longer lived) is higher for PAA-Na. In order to qualitatively estimate the osmotic pressure, we invoke the understanding gained from the theory of dilute solutions of saltfree polyelectrolytes.2,53 It is accepted that there are two major 13912

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913

The Journal of Physical Chemistry B

Article

(13) de Gennes, P. G.; Pincus, P.; Velasco, R. M.; Brochard, F. Remarks on Polyelectrolyte Conformation. J. Phys. (Paris) 1976, 37, 1461−1473. (14) Muthukumar, M. Double Screening in Polyelectrolyte Solutions: Limiting Laws and Crossover Formulas. J. Chem. Phys. 1996, 105, 5183−5199. (15) Yethiraj, A. Conformational Properties and Static Structure Factor of Polyelectrolyte Solutions. Phys. Rev. Lett. 1997, 78, 3789− 3792. (16) Stevens, M. J.; Kremer, K. The Nature of Flexible Linear Polyelectrolytes in Salt Free Solutions: A Molecular Dynamics Study. J. Chem. Phys. 1995, 103, 1669−1690. (17) Winkler, R. G.; Gold, M.; Reineker, P. Collapse of Polyelectrolyte Macromolecules by Counterion Condensation and Ion Pair Formation: A Molecular Dynamics Simulation Study. Phys. Rev. Lett. 1998, 80, 3731−3734. (18) Chang, R.; Yethiraj, A. Brownian Dynamics Simulations of Saltfree Polyelectrolyte Solutions. J. Chem. Phys. 2002, 116, 5284−5298. (19) Chang, R.; Yethiraj, A. Brownian Dynamics Simulations of Polyelectrolyte Solutions with Divalent Counterions. J. Chem. Phys. 2003, 118, 11315−11325. (20) Liu, S.; Muthukumar, M. Langevin Dynamics Simulation of Counterion Distribution Around Isolated Flexible Polyelectrolyte Chains. J. Chem. Phys. 2002, 116, 9975−9982. (21) Liu, S.; Ghosh, K.; Muthukumar, M. Polyelectrolyte Solutions with Added Salt: A Simulation Study. J. Chem. Phys. 2003, 119, 1813− 1823. (22) Molnar, F.; Rieger, J. Like-charge Attraction Between Anionic Polyelectrolytes: Molecular Dynamics Simulations. Langmuir 2005, 21, 786−789. (23) Sulatha, M. S.; Natarajan, U. Origin of the Difference in Structural Behavior of Poly(acrylic acid) and Poly(methacrylic acid) in Aqueous Solution Discerned by Explicit-solvent Explicit-ion MD Simulations. Ind. Eng. Chem. Res. 2011, 50, 11785−11796. (24) Chialvo, A. A.; Simonson, J. M. Solvation Behavior of Shortchain Polystyrene Sulfonate in Aqueous Electrolyte Solutions: A Molecular Dynamics Study. J. Phys. Chem. B 2005, 109, 23031−23042. (25) Carrillo, J. M.; Dobrynin, A. V. Detailed Molecular Dynamics Simulations of a Model NaPSS in Water. J. Phys. Chem. B 2010, 114, 9391−9399. (26) Park, S.; Zhu, X.; Yethiraj, A. Atomistic Simulations of Dilute Polyelectrolyte Solutions. J. Phys. Chem. B 2012, 116, 4319−4327. (27) Obst, M.; Steinbüchel, A. Microbial Degradation of Poly(Amino Acid)s. Biomacromolecules 2004, 5, 1166−1176. (28) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force Field Parameter Sets 53a5 and 53a6. J. Comput. Chem. 2004, 25, 1656−1676. (29) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (30) Miyamoto, S.; Kollman, P. A. SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (31) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-balanced and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (32) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8594. (33) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472. (34) Foloppe, N.; MacKerell, A. D., Jr. All-atom Empirical Force Field for Nucleic Acids: I. Parameter Optimization Based on Small Molecule and Condensed Phase Macromolecular Target Data. J. Comput. Chem. 2000, 21, 86−104. (35) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.;

et al. All-atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (36) Hö chtl, P.; Boresch, S.; Bitomsky, W.; Steinhauser, O. Rationalization of the Dielectric Properties of Common Three-site Water Models in Terms of Their Force Field Parameters. J. Chem. Phys. 1998, 109, 4927−4937. (37) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. A Systematic Study of Water Models for Molecular Simulation: Derivation of Water Models Optimized for Use with a Reaction Field. J. Chem. Phys. 1998, 108, 10220−10230. (38) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (39) Horkay, F.; Hecht, A. M.; Geissler, E. Similarities Between Polyelectrolyte Gels and Biopolymer Solutions. J. Polym. Sci., Part B: Polym. Phys. 2006, 44, 3679−3686. (40) Meng, S.; Ma, J. Solvatochromic Shift of Donor-Acceptor Substituted Bithiophene in Solvents of Different Polarity: Quantum Chemical and Molecular Dynamics Simulations. J. Phys. Chem. B 2008, 112, 4313−4322. (41) Liu, Z.; Huang, S.; Wang, W. A Refined Force Field for Molecular Simulation of Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2004, 108, 12978−12989. (42) Haughney, M.; Ferrario, M.; McDonald, I. R. MolecularDynamics Simulation of Liquid Methanol. J. Phys. Chem. 1987, 91, 4934−4940. (43) Wernet, Ph.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Näslund, L. Å.; Hirsch, T. K.; Ojamäe, L.; Glatzel, P.; et al. The structure of the first coordination shell in liquid water. Science 2004, 304, 995−999. (44) Chen, B.; Siepmann, J. I. Microscopic Structure and Solvation in Dry and Wet Octanol. J. Phys. Chem. B 2006, 110, 3555−3563. (45) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Timneanu, N. Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media. J. Phys. Chem. B 2006, 110, 4393−4398. (46) Starr, F. W.; Nielsen, J. K.; Stanley, H. E. Hydrogen-bond Dynamics for the Extended Simple Point Charge Model of Water. Phys. Rev. E 2000, 62, 579−587. (47) Modig, K.; Pfrommer, B. G.; Halle, B. Temperature-dependent Hydrogen Bond Geometry in Liquid Water. Phys. Rev. Lett. 2003, 90, 075502. (48) Luzar, A.; Chandler, D. Hydrogen-bond Kinetics in Liquid Water. Nature 1996, 379, 55−57. (49) Luzar, A.; Chandler, D. Effect of Environment on Hydrogenbond Dynamics in Liquid Water. Phys. Rev. Lett. 1996, 76, 928−931. (50) Luzar, A. Resolving the Hydrogen-bond Dynamics Conundrum. J. Chem. Phys. 2000, 113, 10663−10675. (51) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 2002. (52) Dubbeldam, D.; Ford, D. C.; Ellis, D. E.; Snurr, R. Q. A New Perspective on the Order-n Algorithm for Computing Correlation Functions. Mol. Simul. 2009, 35, 1084−1097. (53) Netz, R. R.; Andelman, D. Encyclopedia of Electrochemistry; Wiley-VCH: Weinheim, Germany, 2002; Vol. 1, pp 282−322.

13913

dx.doi.org/10.1021/jp406760v | J. Phys. Chem. B 2013, 117, 13906−13913