Hydrogels are highly flexible network polymers being developed as scaffolds for tissue engineering and joint replacement. Their mechanical properties depend largely on their water content. To determine the associated mechanical and thermodynamic properties, we apply the new two-phase thermodynamics method (2PT) to short, molecular dynamics (MD) trajectories of solvated carboxybetaine methacrylate (CBMA) hydrogels. The calculated optimum water content agrees well with recent experiments. We find that the thermodynamics is dominated by a competition between the enthalpy of tightly bound water molecules (which enhance the population of low-energy states of the hydrogel) and the entropy-driven formation of a quasi-liquid water phase in the void volume. These new insights into the role of water in stabilizing hydrophilic motifs is expected to guide design strategies aimed at creating hydrogels with improved performance

Goddard, William A., III

He, Yi

Jiang, Shaoyi

Pascal, Tod A.

Caltech Authors - Main

The thermodynamics of water stabilization of Carboxybetaine 
hydrogels from molecular dynamics simulations 
Tod A. Pascal†‡, Yi Heζ, Shaoyi Jiangζ and William A. Goddard III*†‡ 
†Materials and Process Simulation Center, California Institute of Technology, Pasadena, California, USA 
‡Graduate School of EEWS(WCU), Korea Advanced Institute of Science and Technology, Daejeon, Korea 
ζDepartment of Chemical Engineering, University of Washington, Seattle, Washington, USA 
*Corresponding author: Professor William A Goddard III, wag@wag.caltech.edu 
Supporting Information 
I. THE 2PT METHOD 
From molecular dynamics trajectories in the canonical ensemble, we calculate the 
velocity autocorrelation function (VACF) as:  
( ) ( ) ( )
3
1 1
1lim ' '
2
N
k k
j j jtj k
C t m v t t v t dt
τ
ττ→∞= = −
 
= + 
 
∑∑ ∫    1) 
where N is the number of atoms and ( )kjv t is the k-th component of the velocity of atom j at 
time t.  From a Fourier Transform of the VACF, we obtain the total Density of States (DoS):  
( ) 2( ) lim i tDoS v C t e dt
τ
π ν
τ
τ
−
→∞
−
= ∫     2) 
Here DoS(v) is the number of modes of the system at frequency v, including both vibrational 
and diffusional components. Indeed, DoS(0) measures the diffusion coefficient D0: 
( ) 0120 mDDoS
kT
=       3) 
where m is the mass, k is Boltzmann’s constant and T is the temperature.  
A finite DoS(0) would lead to infinite entropy for standard quantum statistical formula1. 
2PT overcomes this limitation by partitioning DoS(v) into two components:  
• DoSdiff(v): the diffusional component is described as a hard sphere diffusing gas, for 
which the VACF decays exponentially with time, leading to: 
( ) ( )
( )
2
0
0
1
2
diff
diff
DoS
DoS v
DoS v
N
π
=
 
+  
  
   4) 
where Ndiff = 3Nf is the total diffusional degrees of freedom (f is the fraction of modes 
that are diffusive). The diffusional contributions to the thermodynamics, Sdiff, Udiff, Adiff, 
are obtained from the Chapman-Enskog hard sphere theory2, as explained in Ref1. For 
the systems considered here, we find that f ranges from 5% (41p) to 22% (91p). For the 
sake of comparison, a box of SPC-E waters has f of 25%. 
• DoSsolid(v): goes smoothly to zero as v→ 0 (no diffusion), representing a vibrating Debye 
crystal.  Here, we evaluate the partition function Q using the standard harmonic 
oscillator expressions from statistical mechanics 2: 
( ) ( )( )0
exp / 2
ln
1 exp / 2solid
hv
Q DoS v dv
hv
β
β
∞ −
=
− −∫   5) 
Since DoSsolid(v) → 0 as v→ 0, there are no singularities in eqn. (5) at v = 0.  
The system thermodynamics of the solid component are then obtained by integrating 
over the DoSsolid(v). Thus the standard molar entropy S0 is: 
( ) ( ) ( )( )
0 1
,
0
lnln
ln 1 exp
exp 1
N V
solid
QS k Q
T
hvk DoS v hv dv
hv
β
β β
β
−
∞
∂ = +  ∂ 
   = − − −  
−    
∫
  6) 
where β = 1/kT. The internal energy E0 is: 
( ) ( )
0
0
1 1
2 exp 1
MD
solidE E DoS v hv dvhvβ
∞  
= − + 
− 
∫    7) 
where EMD is the system energy from MD. E0 is the reference energy of the system with 
all vibrations are in their lowest (zero) vibrational level, explicitly including zero-point 
energies.  Finally, the Helmholtz Free energy A0 is: 
( ) ( )( )
0 0 1
0 1
0
ln
1 exp
ln
exp / 2solid
A E Q
hv
E DoS v dv
hv
β
β
β
β
−
∞
−
= −
−
= +
−∫
   8) 
We have developed a post trajectory code that performs the 2PT calculations. The code is 
available upon request from the authors (wag@wag.caltech.edu, tpascal@wag.caltech.edu). 
II. COMPUTATIONAL METHODS 
II.a. Molecular Dynamics simulations 
The zwitterionic CBMA hydrogel model was built using an approach similar to that 
reported by Chiessi, et al. 3. The polymer chains of the network were built with Accelrys 
Material Studio (version 3.2). Each of these chains has 15 repeating units in a fully extended 
conformation and a parallel orientation with respect to the Cartesian axes. The polymer 
network was then solvated with various amounts of water molecules to obtain hydrogels with 
various water contents. The all-atom Consistent Valence Force Field (CVFF)4 was used to 
describe the interactions in the polymer networks. The partial charges of each atom in the 
CBMA polymer were calculated using the Jaguar program 5. These calculations were carried out 
at the Hartree-Fock (HF) level using the 6-31G** basis set. In this procedure, the electrostatic 
field at a grid of points was calculated from the HF wave function. Using the grid points outside 
the VDW radii, atom-centered charges were derived so as to match the HF potential while 
reproducing the dipole moment from HF. The water molecules are described with the Extended 
Simple Point Charge (SPC/E) model.6  
Our MD simulations were performed using the LAMMPS7,8 simulation engine, which 
affords the flexibility of using various forcefields in a common framework.  Long-range 
coulombic interactions were calculated using the particle-particle particle-mesh Ewald method9 
(with a precision of 10-5 kcal/mol), while the van der Waals interaction were computed with a 
cubic spline (inner cutoff of 11Å and an outer cutoff of 12Å). We used the spline to guarantee 
that the energies and forces go smoothly to zero at the outer cutoff, preventing energy drifts 
that might arise from to inconsistent forces.  We also tested the effect of the cutoff by 
computing the energy of benzene with cutoffs ranging from 8 to 20Å and found converged 
results at 12Å. 
For each system, we used the Continuous Configurational Boltzmann Biased (CCBB) 
Monte Carlo (MC) method10,11 to generate a random starting structure of 512 solvent molecules 
packed to minimize the system interaction energy. To rapidly equilibrate the systems, we used 
our standard procedure12-14: after an initial conjugant gradient minimization to an RMS force of 
10-4kcal/mol/Å, the system was slowly heated from 0K to 298K over a period of 100 ps using a 
Langevin thermostat in the constant temperature, constant volume canonical (NVT) ensemble. 
The temperature coupling constant was 0.1 ps and the simulation timestep was 1.0 fs.  
This equilibration was followed by 20ns of constant-pressure(iso-baric), constant-
temperature (NPT) dynamics at 298K and 1 atm. The temperature coupling constant was 0.1 ps 
while the pressure piston constant was 2.0 ps. The equations of motion used are those of 
Shinoda et al.15, which combine the hydrostatic equations of Martyna et al.16with the strain 
energy proposed by Parrinello and Rahman17. The time integration schemes closely follow the 
time-reversible measure-preserving Verlet integrators derived by Tuckerman et al.18.  
Production dynamics was then run for a further 5.0 ns in the NPT ensemble. All data 
reported in this paper are the statistical averages of two individual MD simulations from two 
different initial configurations.  
II.b. Free energy calculations 
During the final 5ns of the production dynamics simulations outlined above, we selected 
snapshots of the system (coordinates and velocities) every 0.5 ns. Each of the 10 snapshots was 
then simulated for 20 ps of MD using the Gibbs (constant particle, constant volume, constant 
temperature – NVT) ensemble.  The velocities and coordinates were saved every 4 fs (must be 
shorter than the fastest vibrational levels, which have periods of ~10 fs for the 3000 cm-1 C-H 
vibrations).  
 III. TABLES 
Table S1: Description of CBMA – water systems considered in this study. The quoted average 
volumes and solvent accessible surface areas are obtained by averaging 1000 structure during 
the last 4ns of MD. 
% water # num a<SASA> Å2 <Volume> Å3  <SASA/V> 
41 592 17877.5 37481.0 0.477 
49 819 20764.6 43960.8 0.472 
62 1396 16668.4 61287.4 0.272 
71 2120 14987.1 83048.0 0.180 
81 3630 15739.8 126859.0 0.124 
91 8200 16867.4 248168.0 0.068 
aSolvent accessible surface area, obtained according to the method of Lee-Richards19 with a 
1.4Å probe radius 
 
Table S2: Thermodynamics of the CBMA hydrogel from 2PT analysis. 
% water A0 (kJ/mol) E0 (kJ/mol) ZPE (kJ/mol)a Emd (kJ/mol)b S0 (J/mol/K) 
 avg ± avg ± avg ± avg ± avg ± 
41 56630.1 210.4 62712.5 188.4 53815.9 158.0 23716.7 86.9 20274.7 130.6 
49 56501.6 161.8 62563.9 151.1 53717.5 141.9 23439.3 78.7 20207.7 120.5 
62 56218.3 157.1 62194.7 176.8 53380.9 129.3 23232.1 100.8 19921.3 55.7 
71 55948.8 114.9 61900.9 126.9 53341.6 104.8 22912.3 90.9 19840.3 84.8 
81 56589.8 182.1 62475.6 171.1 53355.1 169.0 22420.1 57.8 19619.3 101.7 
91 56780.1 186.8 62756.5 154.8 53618.2 104.5 23408.7 173.2 19921.2 152.1 
dry: 0%           
aZero point energy determined using equation 7. 
bTotal Energy of hydrogel from MD. The per-atom long range electrostatic contributions were 
obtained by 1000 single point energy calculations using Ewald summations along the 20ps 
trajectory, as detailed in ref20.  
  
Table S3: Thermodynamics per molecule of waters in CBMA hydrogel from 2PT analysis 
% water A0 (kJ/mol) E0 (kJ/mol) Cv (kJ/mol/K)a S0 (J/mol/K) 
 avg ± avg ± avg ± avg ± 
41 -51.74 0.08 -35.99 0.08 35.93 0.04 50.67 0.23 
49 -52.03 0.11 -35.73 0.08 35.88 0.06 52.64 0.21 
62 -52.57 0.10 -35.17 0.09 35.92 0.07 56.42 0.38 
71 -52.81 0.06 -34.94 0.06 35.64 0.04 58.30 0.23 
81 -52.94 0.03 -34.80 0.03 35.35 0.05 59.53 0.13 
91 -52.86 0.04 -34.77 0.01 34.96 0.05 59.42 0.08 
bulkb -53.2 0.03 -34.61 0.02 34.79 0.04 62.76 0.11 
aConstant volume heat capacity 
bBulk water values obtained by 100 separate calculations of a 1200 SPC-E water molecules 
simulation over 5ns 
 
Table S4: Components of entropy (J/mol/K/molecule) of waters in CBMA hydrogel from 2PT 
analysis 
% water Diffusional Librational Rotational 
 avg ± avg ± avg ± 
41 0.27 0.17 41.94 0.47 8.46 0.05 
49 1.77 0.38 41.98 0.06 8.89 0.06 
62 3.10 0.22 43.51 0.24 9.81 0.11 
71 7.00 0.41 41.15 0.31 10.15 0.05 
81 8.43 0.11 40.72 0.20 10.38 0.03 
91 9.07 0.24 39.84 0.35 10.52 0.04 
bulk 11.10 0.04 38.99 0.10 10.67 0.03 
 
IV. REFERENCES 
  (1) Lin, S. T.; Blanco, M.; Goddard, W. A. J. Chem. Phys. 2003, 119, 11792.  (2) McQuarrie, D. A. Statistical mechanics; University Science Books: Sausalito (Calif.), 2000.  (3) Chiessi, E.; Cavalieri, F.; Paradossi, G. J Phys Chem B 2007, 111, 2820.  (4) Dauberosguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins 1988, 4, 31.  (5) Schrödinger, I.  Portland, OR, 1991-2000.  (6) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 
1987, 91, 6269.  (7) Plimpton, S. J.; Pollock, R.; Stevens, M. In Proc of the Eighth SIAM 
Conference on Parallel Processing for Scientific Computing Minneapolis, MN 1997.  (8) Plimpton, S. J. Comput. Phys. 1995, 117, 1.  (9) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using 
Particles; Taylor & Francis: New York, 1989.  (10) Maiti, P. K.; Çaǧın, T.; Wang, G.; Goddard, W. A. Macromolecules 
2004, 37, 6236.  (11) Cagin, T.; Wang, G. F.; Martin, R.; Breen, N.; Goddard, W. A. 
Nanotechnology 2000, 11, 77. 
 (12) Maiti, P. K. P., Tod A.; Vaidehi, Nagarajan; Goddard, William A. J. 
Nanosci. Nanotechnol. 2007, 7, 1712.  (13) Maiti, P. K.; Pascal, T. A.; Vaidehi, N.; Heo, J.; Goddard, W. A. 
Biophys. J. 2006, 90, 1463.  (14) Maiti, P. K.; Pascal, T. A.; Vaidehi, N.; Goddard, W. A. NUCLEIC 
ACIDS RES 2004, 32, 6047.  (15) Shinoda, W.; Shiga, M.; Mikami, M. Physical Review B 2004, 69.  (16) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177.  (17) Parrinello, M.; Rahman, A. Journal of Applied Physics 1981, 52, 7182.  (18) Tuckerman, M. E.; Alejandre, J.; Lopez-Rendon, R.; Jochim, A. L.; Martyna, G. J. J Phys a-Math Gen 2006, 39, 5629.  (19) Lee, B.; Richards, F. M. J. Mol. Biol. 1971, 55, 379.  (20) Kumar, H.; Mukherjee, B.; Lin, S. T.; Dasgupta, C.; Sood, A. K.; Maiti, P. K. J. Chem. Phys. 2011, 134.  
 


Thermodynamics of Water Stabilization of Carboxybetaine Hydrogels from Molecular Dynamics Simulations

https://authors.library.caltech.edu/24832/2/jz200760n_si_001.pdf

Thermodynamics of Water Stabilization of Carboxybetaine Hydrogels from Molecular Dynamics Simulations

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main