- Email: [email protected]
Predict the glass transition temperature of glycerol–water binary cryoprotectant by molecular dynamic simulation
Available online at www.sciencedirect.com
Cryobiology 56 (2008) 114–119 www.elsevier.com/locate/ycryo
Predict the glass transition temperature of glycerol–water binary cryoprotectant by molecular dynamic simulation q Dai-Xi Li a, Bao-Lin Liu
a,* ,
Yi-shu Liu b, Cheng-lung Chen
c
a
University of Shanghai for Science and Technology, Shanghai 200093, China Guiyang College of Traditional Chinese Medicine, Guiyang 550001, China Department of Chemistry, National Sun Yat-Sen University, Kaohsiung 80424, China b
c
Received 22 September 2007; accepted 28 November 2007 Available online 5 December 2007
Abstract Vitrification is proposed to be the best way for the cryopreservation of organs. The glass transition temperature (Tg) of vitrification solutions is a critical parameter of fundamental importance for cryopreservation by vitrification. The instruments that can detect the thermodynamic, mechanical and dielectric changes of a substance may be used to determine the glass transition temperature. Tg is usually measured by using differential scanning calorimetry (DSC). In this study, the Tg of the glycerol-aqueous solution (60%, wt/%) was determined by isothermal-isobaric molecular dynamic simulation (NPT-MD). The software package Discover in Material Studio with the Polymer Consortium Force Field (PCFF) was used for the simulation. The state parameters of heat capacity at constant pressure (Cp), density (q), amorphous cell volume (Vcell) and specific volume (Vspecific) and radial distribution function (rdf) were obtained by NPT-MD in the temperature range of 90–270 K. These parameters showed a discontinuity at a specific temperature in the plot of state parameter versus temperature. The temperature at the discontinuity is taken as the simulated Tg value for glycerol–water binary solution. The Tg values determined by simulation method were compared with the values in the literatures. The simulation values of Tg (160.06– 167.51 K) agree well with the DSC results (163.60–167.10 K) and the DMA results (159.00 K). We drew the conclusion that molecular dynamic simulation (MDS) is a potential method for investigating the glass transition temperature (Tg) of glycerol–water binary cryoprotectants and may be used for other vitrification solutions. 2007 Elsevier Inc. All rights reserved. Keywords: Molecular dynamics simulation; Glass transition temperature; Vitrification; Cryoprotectant
Various types of cells and tissues have been successfully cryopreserved by vitrification since its appearance in cryobiology [23,29,20]. Vitrification is also proposed to be the best way for the cryopreservation of organs [9,10]. The glass transition temperature (Tg) is one of the critical parameters of fundamental importance for cryopreservation by vitrification. Therefore, it is important to know the glass transition temperature for vitrification solutions. q Statement of funding: The research was supported by National Science Foundation of China (50576059), The Shanghai Leading Academic Discipline Project (PO502), The Shanghai Development Funds (05ZZ24) as well as The Shanghai University for Science and Technology (X544). * Corresponding author. Fax: +86 21 55270596. E-mail address: [email protected] (B.-L. Liu).
0011-2240/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.cryobiol.2007.11.003
The glass transition has been investigated experimentally and theoretically in the fields of material and physics for a long time [3,28]. Though the vast amounts of literatures have been already contributed to the investigation of the glass transition and glass transition temperature [3,8,4,5,24,26], the nature of the glass transition still remains an unsolved and challenging problem [7,25]. In fact, up to now, there is no single theory that can account for all the features during the glass transition process [27]. Experimentally, the glass transition can be tracked by differential scanning calorimeter (DSC) [12]. Accordingly, Tg can be determined by DSC characterized by an abrupt change of heat flow in the DSC thermogram [15]. When solution is cooled and changed to glass, the heat capacity
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Cp shows an abrupt change at the glass transition temperature Tg. The molecular simulation has been applied to investigate the glass transition and predict the glass transition temperature in polymer material sciences [1,13,14,18,22,16,2]. Within an equilibrium phase (gas, liquid, or crystal), the thermodynamic parameters (heat capacity, special volume, and so on) will change regularly with temperature. These parameters exhibit a discontinuity when the substance transfers from one equilibrium phase to another. So the plot of these parameters versus temperature will show an abrupt change when the glass transition takes place. To the best of our knowledge, there is still no report on the prediction of glass transition temperature for cryoprotectants by MDS method. In this paper, the glass transition temperature of 60 wt/% glycerol-aqueous solution was determined by molecular dynamic simulation. Theoretical background and methods When a solution is cooled and turned into its glassy state, the heat capacity Cp, density, specific volume, radial distribution function, as well as the formation probability of the H-bonds will show an abrupt change at the glass transition temperature. In theory, the above parameters can be evaluated by using the molecular dynamics simulation. For example, the heat capacity Cp can be expressed as a function of the kinetic and potential energy (j, m), temperature (T), pressure (P) and volume (V) of the system during molecular simulation [19], that is: E 1 D Cp ¼ dðj þ m þ PV Þ2 2 kBT where kB is Boltzmann’s constant (1.38066 · 1023J K1). The notation d(j + m + PV) means (j + m+PV) Æj + m + PVæ, namely energy fluctuation, where Æj + m + PVæ denotes the equilibrium ensemble average value of quantity (j + m + PV). The density, specific volume, radial distribution function, as well as the formation probability of the H-bonds were also investigated in this paper. The software package Discover in Material Studio (Accelrys Software, Inc.) with the Polymer Consortium Force Field (PCFF) is used for predicting the glass transition temperature. Discover is a molecular simulation program in Material Studio. It provides a broad range of simulation methods to study molecular systems. Discover can also be used to perform structural characterization and property prediction for molecules, biological compounds, and materials. The force fields in Discover provide the function of predicting the structure, energetic, and properties of organic, inorganic, and biological systems as well as gas, liquid and solid. Molecular dynamics simulations generate information at the microscopic level, including atomic positions and velocities. The isothermal-isobaric molecular dynamic simulations (NPT-MD) method is employed. In NPT ensemble, tem-
115
perature change is controlled by using the velocity scaling method, pressure (0.0001 GPa) is controlled by the Andersen method and stress is controlled by the Parrinello method. The modules Minimizer, Discover and Amorphous Cell Tools in the Materials Studio (Accelrys, Inc.) are also used here. For the molecular simulation of the glycerol-aqueous solution, the periodic boundary condition with an amorphous cell is applied. As shown in Fig. 1, the cell is composed of 100 glycerol molecules and 341 water molecules and contains 2423 atoms totally. The potential energy of each atom is derived from the PCFF force field. The cutoff distance for van der Waals and electrostatic ˚. forces is 9.50 A The potential energy of the system is minimized by use of the Smart Minimizer method. The temperature range for the molecular simulation is set from 90 K to 270 K with an interval of 20 K. The pre-equilibration of aqueous solution at each temperature is attained from molecular simulations of every 1 fs during 5000 ps. Then molecular simulation is performed thrice or four times by using NPT ensemble for 200 ps in order to get a reliable average value. Each trajectory is sampled every 20 fs, and 10,000 structures are subject to analysis. Based on the data of molecular simulation, the values of the heat capacity at constant pressure (Cp), the density of the solution (q) and the radial distribution function and so on are further extracted and discussed. Results and discussion Prediction of the Tg by the heat capacity Cp Based on the equilibrium thermodynamic theory of fluctuations, a linear law of the heat capacity with temperature in an equilibrium phase is derived: the heat capacity of the system increases linearly with temperature. An abrupt change of the heat capacity at a specific temperature indicates that some solid–liquid transitions of a substance may occur. From 90 K to 150 K, the slope of the curve of the heat capacity vs. temperature (Fig. 2) is 0.01283 kJ mol1 K2 (R = 0.99917). From 170 K to 270 K, the slope of the curve of the heat capacity vs. temperature changes to 0.02056 kJ mol1 K2 (R = 0.99989). When both the curves are extrapolated, they meet with each other at a temperature of 165.50 K. This temperature is the so-called glass transition temperature of the glycerol-aqueous solution. Prediction of the Tg by the density and the specific volume (Vspecific) The curves of the density of the glycerol-aqueous solution vs. temperature are shown in the Fig. 3. Below 150 K, the density increases linearly with the decrement of temperature, the slope of the curve is 4.144 · 104 g cm3 K1 (R = 0.9956). Above 170 K,
116
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Fig. 1. Amorphous cell for the glycerol-aqueous solution (60% wt/%). The amorphous cell is composed of 100 glycerol molecules and 341 water molecules. The configurations of a glycerol molecule and a water molecule are shown on the right of the figure using the ball and stick style. The black balls, gray balls and white balls denote oxygen atoms, carbon atoms and hydrogen atoms, respectively.
Fig. 2. MDS curves of the heat capacity vs. temperature for glycerolaqueous solution (60%, wt/%).
Fig. 3. MDS curves of the density vs. temperature for glycerol-aqueous solution (60%, wt/%).
the slope of the curve is 8.194 · 104 g cm3 K1 (R = 0.9993). The temperature of 164.67 K at the intersection of the two linear curves is taken as the glass transition temperature. Specific volume is mostly used to determine the value of Tg [1,18,16,6]. The curves of the specific volume vs. temperature are shown in Fig. 4. The intersections of the two lines yield the value of Tg (167.51 K), which is 2.84 K higher than 164.67 K.
the rdf peaks implies that the arrangement of atoms and molecules changes gradually from disorder to order with the decrease of the temperature. The first rdf peak ˚ ) denotes the formation probability of the H(r1 = 1.90 A bonds in the glycerol-aqueous solution, which indicates the H-bond strengthens with the decrease of temperature. A plot of P(r) vs. temperature is drawn (Fig. 6). It can be seen from Fig. 6 that the p(r) value for 160 K shows a clearly different change pattern, implying a phase change around 160 K.
Prediction of the Tg by the formation of the H-bonds Verification of the predicted Tg The radial distribution functions (rdfs) between all O atoms and all H on the OH at a temperature range from 90 K to 270 K are shown in the Fig. 5. The pattern of
Gao et al. [11] studied extensively the Tg of glycerol solutions by DSC and obtained the Tg of 60% (wt/%)
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Fig. 4. MDS curves of the specific volume vs. temperature for glycerolaqueous solution (60%, wt/%).
Fig. 5. The radial distribution functions between all O atoms and all H on the OH at various temperatures. The radial distribution functions express the distance distribution between all O atoms and all H on the OH at various temperatures. The rdfs (radial distribution functions) changes with ˚) temperatures (as marked by the arrows). The first peak (r1 1.90A corresponds to the formation probability of the H-bond between all O atoms and all H on the OH. The other peak correspond to the farther distance distribution between all O atoms and all H on the OH. It can be seen that the formation probability of the H-bond between all O atoms and all H on the OH in the solution increases with the decrease of temperature, and simultaneously their bond distances shorten step by step.
glycerol-aqueous solution (Table 1). The Tg of 60% (wt/%) glycerol water solution measured by Luyet and Rasmussen [21] using a DTA is 159 K. As shown in Table 1, the relative errors between the Tg values predicted by MDS and the Tg values in literature [11,21] are all less than 5.5%, which indicates that the NPT-MD method can predict the Tg of glycerol-aqueous solution with relatively high accuracy. Discussion and conclusions Glass transition occurs when a sample goes from a hard, glass like state to a rubber like state. The state parameters
117
Fig. 6. Probability of the first rdf peak vs. temperature for glycerolaqueous solution (60%, wt/%). It is shown that the H-bond becomes stronger with the decrease of temperature. The intersections of the two lines yield the value of Tg (160.06 K), which implies that the microstructures at the temperature above and below Tg are quite different in nature.
of the sample will exhibit an abrupt change within the range of the glass transition temperature. So any instrument that can detect the thermal, mechanical or dielectric changes can theoretically be used to measure the Tg of a material. For example, DSC defines the glass transition as a change in the heat capacity, TMA defines the glass transition in terms of the change in the coefficient of thermal expansion. While Tg can be measured by a variety of methods, it is a time-consuming procedure, especially if the sample is to be kept at subzero temperatures, in anhydrous conditions, or if sampling a portion of the specimen for analysis is cumbersome. Hence, predicting rather than directly measuring Tg as a function of the content of the constituents of a solution can be a powerful tool [17]. Molecular dynamics (MD) is the term used to describe the solution of the classical equation of motion (Newton’s equation) for a set of molecules. In this case, the particles move at constant velocity between perfectly elastic collisions. If the force exerted on each atom is known, it is possible to determine the acceleration of each atom in the system. Integration of the motion equations yields a trajectory that describes the positions, velocities and accelerations of the particles. The average values of properties can be determined from this trajectory. This method is deterministic, that is, once the positions and velocities of each atom are known, the state of the system can be predicted at any time in the future or the past. MD method is different with the Monte Carlo method, which is a statistical technique that involves using random numbers and probability to solve problem. In this study, the heat capacity, the density and specific volume and rdf are used in the software package Discover in Material Studio to predict the Tg of 60% (wt/%) glycerol-aqueous solution. The heat capacity reveals the thermal change when the solution transfers from rubbery state to glass state. While the density and specific volume indicate the mechanical change.
118
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Table 1 The Tg of glycerol (60% wt/%) aqueous solution determined by DSC From literature Gao et al. Tg (K) Relative errors (%)
Predicted in this paper Luyet et al.
163.60 159.00 Compared to Gao et al. Compared to Luyet et al.
The predicted values were verified with the values from literatures. The relative errors are all less than 5.5%. So the calculation of Tg using this software package is relatively accurate. The deviation may be due to differences in glycerol source, method of measurement or the heating rate in DSC or DTA studies. The differences among the three predicted Tg values result from the difference between the dynamics and thermodynamics properties. The vitrification solutions for cryopreservation usually contain several solutes, which are more complex than the 60 wt% glycerol solution. Theoretically, the molecular dynamics simulation method can be applied to calculate the Tg of realistic complex vitrification solutions if the types and the concentration of the components in the solution are known. But it is estimated that half of one year will be cost to finish the calculation because the cell (similar with Fig. 1) will be much bigger for the realistic vitrification solutions. If the direct calculation method is not practical, another way may be used to predict the Tg of realistic vitrification solutions. Firstly, the Tg of each separate component in the vitrification solutions can be calculated using the MD method. Then the Tg of this vitrification solution can be determined by use of the models as described in Ref. [17]. This will be the future work of our study. In conclusion, molecular dynamic simulation (MDS) is a potential way for investigating the glass transition temperature (Tg) of glycerol–water binary cryoprotectants and may be used for other types of cryoprotectants. Acknowledgments The research was supported by National Science Foundation of China (50576059), The Shanghai Leading Academic Discipline Project (PO502), The Shanghai Development Funds (05ZZ24) as well as the Shanghai University for Science and Technology(X544). References [1] B.F. Abu-Sharkh, Glass transition temperature of poly(vinylchloride) from molecular dynamics simulation: explicit atom model versus rigid CH2 and CHCl groups model, Comput. Theor. Polym. Sci. 11 (2001) 29–34. [2] V.L. Alexey, M.A.J. Michels, Simulation of polymer glasses: from segmental dynamics to bulk mechanics, J. Non-Cryst. Solids 352 (2006) 5008–5012. [3] C.A. Angell, Formation of glasses from liquids and biopolymers, Science 267 (5206) (1995) 1924–1935.
q
Cp
Vspecific
rdf
164.67 0.65 3.57
165.50 1.16 4.09
167.51 2.39 5.35
160.06 3.54 1.06
[4] C.A. Angell, K.L. Ngai, G.B. Mackenna, P.F. McMillan, S.W. Martin, Relaxation in glass forming liquids and amorphous solids, J. Appl. Phys. 88 (2000) 3113–3157. [5] P.G. Debenedetti, F.H. Stillinger, Supercooled liquids and the glass transition, Nature 410 (2001) 259–267. [6] V.I. Dimitrov, Theory of fluidity of liquids, glass transition, and melting, J. Non-Cryst. Solids 352 (2006) 216–231. [7] M.D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, Annu. Rev. Phys. Chem. 51 (2000) 99–128. [8] M.D. Ediger, C.A. Angell, S.R. Nagel, Supercooled liquids and glasses, J. Phys. Chem. 100 (1996) 13200–13212. [9] G.M. Fahy, Vitrification, in: J.J. McGrath, K.R. Diller (Eds.), Low Temperature Biotechnology: Engineering Applications and Engineering Contributions, The American Society of Mechanical Engineers, New York, 1988, pp. 113–146. [10] G.M. Fahy, D.R. Macfarlane, C.A. Angell, H.T. Meryman, Vitrification as an approach to cryopreservation, Cryobiology 21 (1984) 407–426. [11] C. Gao, G.-Y. Zhou, Y. Xu, T.C. Hua, Glass transition and enthalpy relaxation of glycerol and its aqueous solution, Thermochim. Acta 435 (2005) 38–43. [12] N.A. Grunina, T.V. Belopolskaya, G.I. Tsereteli, DSC study the glass transition process in humid biopolymers, J. Phys. Conference Series 40 (2006) 105–110. [13] V.R. Gumen, F.R. Jones, D. Attwood, Prediction of the glass transition temperatures for epoxy resins and blends using group interaction modeling, Polymer 42 (2001) 5717–5725. [14] P. Javier, B. Juan, Glass transition temperature of low molecular weight poly (3-aminopropyl methyl siloxane). A molecular dynamics study, Polymer 43 (2002) 6049–6055. [15] M. Jochem, C.H. Korber, Extended phase diagrams for the ternary solutions H2O–NaCl–glycerol and H2O–NaCl–hydroxyethylstarch determined by DSC, Cryobiology 24 (1987) 513– 536. [16] G.W. Karl, M. Martin, K. Andreas, Z. Gerhard, Glass transition temperature of a cationic polymethacrylate dependent on the plasticizer content—simulation vs. experiment, Chem. Phys. Lett. 406 (2005) 90–94. [17] I. Katkov, F. Levine, Prediction of the glass transition temperature of water solutions: comparison of different models, Cryobiology 49 (2004) 62–82. [18] B. Kurt, B. Jorg, W. Paul, Glass transition of polymer melts: test of theoretical concepts by computer simulation, Prog. Polym. Sci. 28 (2003) 115–172. [19] A.R. Leach, Molecular Modeling—Principles and Applications, second ed., Person Education Limited, Harlow (London), England, 2001, pp. 348–349. [20] B.L. Liu, J. McGrath, L. McCabe, M. Crimp, Response of murine osteoblasts and porous hydroxyapatite (HA) scaffolds to two-step, slow freezing and vitrification process, Cell Preservation Technol. 1 (2002) 33–44. [21] B. Luyet, D. Rasmussen, Study by differential thermal analysis of the temperatures of instability of rapidly cooled solutions of glycerol, ethylene glycol, sucrose and glucose, Biodynamica 10 (1968) 167–191. [22] W. Paul, Molecular dynamics simulations of the glass transition in polymer melts, Polymer 45 (2004) 3901–3905.
D.-X. Li et al. / Cryobiology 56 (2008) 114–119 [23] W.F. Rall, G.M. Fahy, Ice-free cryopreservation of mouse embryos at 196 C by vitrification, Nature 313 (1975) 573– 575. [24] D.R. Reichman, P. Charbonneau, Mode-coupling theory, J. Stat. Mech. P05013 (2005) 1–23. [25] R. Richert, C.A. Angell, Dynamics of glass-forming liquids. V. On the link between molecular dynamics and configurational entropy, J. Chem. Phys. 108 (1998) 9016–9026. [26] F. Sciortino, Potential energy landscape description of supercooled liquids and glasses, J. Stat. Mech. P05015 (2005) 1–35.
119
[27] M. Scott Shell, Pablo G. Debenedetti, Athanassios Z. Panagiotopoulos, A conformal solution theory for the energy landscape and glass transition of mixtures, Fluid Phase Equilibr. 241 (2006) 147–154. [28] F.H. Strillinger, A topographic view of supercooled liquids and glass formation, Science 267 (5206) (1995) 1935–1939. [29] M.J. Taylor, Y.C. Song, B.S. Kheirabadi, F.G. Lightfoot, K.G.M. Brockbank (Eds.), Vitrification fulfills its promise as an approach to reducing freeze-induced injury in a multi-cellular tissue, in: Advances in Heat and Mass Transfer in Biotechnology, vol. 44, Plenum Press, NY, 1999, pp. 93–102.
Cryobiology 56 (2008) 114–119 www.elsevier.com/locate/ycryo
Predict the glass transition temperature of glycerol–water binary cryoprotectant by molecular dynamic simulation q Dai-Xi Li a, Bao-Lin Liu
a,* ,
Yi-shu Liu b, Cheng-lung Chen
c
a
University of Shanghai for Science and Technology, Shanghai 200093, China Guiyang College of Traditional Chinese Medicine, Guiyang 550001, China Department of Chemistry, National Sun Yat-Sen University, Kaohsiung 80424, China b
c
Received 22 September 2007; accepted 28 November 2007 Available online 5 December 2007
Abstract Vitrification is proposed to be the best way for the cryopreservation of organs. The glass transition temperature (Tg) of vitrification solutions is a critical parameter of fundamental importance for cryopreservation by vitrification. The instruments that can detect the thermodynamic, mechanical and dielectric changes of a substance may be used to determine the glass transition temperature. Tg is usually measured by using differential scanning calorimetry (DSC). In this study, the Tg of the glycerol-aqueous solution (60%, wt/%) was determined by isothermal-isobaric molecular dynamic simulation (NPT-MD). The software package Discover in Material Studio with the Polymer Consortium Force Field (PCFF) was used for the simulation. The state parameters of heat capacity at constant pressure (Cp), density (q), amorphous cell volume (Vcell) and specific volume (Vspecific) and radial distribution function (rdf) were obtained by NPT-MD in the temperature range of 90–270 K. These parameters showed a discontinuity at a specific temperature in the plot of state parameter versus temperature. The temperature at the discontinuity is taken as the simulated Tg value for glycerol–water binary solution. The Tg values determined by simulation method were compared with the values in the literatures. The simulation values of Tg (160.06– 167.51 K) agree well with the DSC results (163.60–167.10 K) and the DMA results (159.00 K). We drew the conclusion that molecular dynamic simulation (MDS) is a potential method for investigating the glass transition temperature (Tg) of glycerol–water binary cryoprotectants and may be used for other vitrification solutions. 2007 Elsevier Inc. All rights reserved. Keywords: Molecular dynamics simulation; Glass transition temperature; Vitrification; Cryoprotectant
Various types of cells and tissues have been successfully cryopreserved by vitrification since its appearance in cryobiology [23,29,20]. Vitrification is also proposed to be the best way for the cryopreservation of organs [9,10]. The glass transition temperature (Tg) is one of the critical parameters of fundamental importance for cryopreservation by vitrification. Therefore, it is important to know the glass transition temperature for vitrification solutions. q Statement of funding: The research was supported by National Science Foundation of China (50576059), The Shanghai Leading Academic Discipline Project (PO502), The Shanghai Development Funds (05ZZ24) as well as The Shanghai University for Science and Technology (X544). * Corresponding author. Fax: +86 21 55270596. E-mail address: [email protected] (B.-L. Liu).
0011-2240/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.cryobiol.2007.11.003
The glass transition has been investigated experimentally and theoretically in the fields of material and physics for a long time [3,28]. Though the vast amounts of literatures have been already contributed to the investigation of the glass transition and glass transition temperature [3,8,4,5,24,26], the nature of the glass transition still remains an unsolved and challenging problem [7,25]. In fact, up to now, there is no single theory that can account for all the features during the glass transition process [27]. Experimentally, the glass transition can be tracked by differential scanning calorimeter (DSC) [12]. Accordingly, Tg can be determined by DSC characterized by an abrupt change of heat flow in the DSC thermogram [15]. When solution is cooled and changed to glass, the heat capacity
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Cp shows an abrupt change at the glass transition temperature Tg. The molecular simulation has been applied to investigate the glass transition and predict the glass transition temperature in polymer material sciences [1,13,14,18,22,16,2]. Within an equilibrium phase (gas, liquid, or crystal), the thermodynamic parameters (heat capacity, special volume, and so on) will change regularly with temperature. These parameters exhibit a discontinuity when the substance transfers from one equilibrium phase to another. So the plot of these parameters versus temperature will show an abrupt change when the glass transition takes place. To the best of our knowledge, there is still no report on the prediction of glass transition temperature for cryoprotectants by MDS method. In this paper, the glass transition temperature of 60 wt/% glycerol-aqueous solution was determined by molecular dynamic simulation. Theoretical background and methods When a solution is cooled and turned into its glassy state, the heat capacity Cp, density, specific volume, radial distribution function, as well as the formation probability of the H-bonds will show an abrupt change at the glass transition temperature. In theory, the above parameters can be evaluated by using the molecular dynamics simulation. For example, the heat capacity Cp can be expressed as a function of the kinetic and potential energy (j, m), temperature (T), pressure (P) and volume (V) of the system during molecular simulation [19], that is: E 1 D Cp ¼ dðj þ m þ PV Þ2 2 kBT where kB is Boltzmann’s constant (1.38066 · 1023J K1). The notation d(j + m + PV) means (j + m+PV) Æj + m + PVæ, namely energy fluctuation, where Æj + m + PVæ denotes the equilibrium ensemble average value of quantity (j + m + PV). The density, specific volume, radial distribution function, as well as the formation probability of the H-bonds were also investigated in this paper. The software package Discover in Material Studio (Accelrys Software, Inc.) with the Polymer Consortium Force Field (PCFF) is used for predicting the glass transition temperature. Discover is a molecular simulation program in Material Studio. It provides a broad range of simulation methods to study molecular systems. Discover can also be used to perform structural characterization and property prediction for molecules, biological compounds, and materials. The force fields in Discover provide the function of predicting the structure, energetic, and properties of organic, inorganic, and biological systems as well as gas, liquid and solid. Molecular dynamics simulations generate information at the microscopic level, including atomic positions and velocities. The isothermal-isobaric molecular dynamic simulations (NPT-MD) method is employed. In NPT ensemble, tem-
115
perature change is controlled by using the velocity scaling method, pressure (0.0001 GPa) is controlled by the Andersen method and stress is controlled by the Parrinello method. The modules Minimizer, Discover and Amorphous Cell Tools in the Materials Studio (Accelrys, Inc.) are also used here. For the molecular simulation of the glycerol-aqueous solution, the periodic boundary condition with an amorphous cell is applied. As shown in Fig. 1, the cell is composed of 100 glycerol molecules and 341 water molecules and contains 2423 atoms totally. The potential energy of each atom is derived from the PCFF force field. The cutoff distance for van der Waals and electrostatic ˚. forces is 9.50 A The potential energy of the system is minimized by use of the Smart Minimizer method. The temperature range for the molecular simulation is set from 90 K to 270 K with an interval of 20 K. The pre-equilibration of aqueous solution at each temperature is attained from molecular simulations of every 1 fs during 5000 ps. Then molecular simulation is performed thrice or four times by using NPT ensemble for 200 ps in order to get a reliable average value. Each trajectory is sampled every 20 fs, and 10,000 structures are subject to analysis. Based on the data of molecular simulation, the values of the heat capacity at constant pressure (Cp), the density of the solution (q) and the radial distribution function and so on are further extracted and discussed. Results and discussion Prediction of the Tg by the heat capacity Cp Based on the equilibrium thermodynamic theory of fluctuations, a linear law of the heat capacity with temperature in an equilibrium phase is derived: the heat capacity of the system increases linearly with temperature. An abrupt change of the heat capacity at a specific temperature indicates that some solid–liquid transitions of a substance may occur. From 90 K to 150 K, the slope of the curve of the heat capacity vs. temperature (Fig. 2) is 0.01283 kJ mol1 K2 (R = 0.99917). From 170 K to 270 K, the slope of the curve of the heat capacity vs. temperature changes to 0.02056 kJ mol1 K2 (R = 0.99989). When both the curves are extrapolated, they meet with each other at a temperature of 165.50 K. This temperature is the so-called glass transition temperature of the glycerol-aqueous solution. Prediction of the Tg by the density and the specific volume (Vspecific) The curves of the density of the glycerol-aqueous solution vs. temperature are shown in the Fig. 3. Below 150 K, the density increases linearly with the decrement of temperature, the slope of the curve is 4.144 · 104 g cm3 K1 (R = 0.9956). Above 170 K,
116
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Fig. 1. Amorphous cell for the glycerol-aqueous solution (60% wt/%). The amorphous cell is composed of 100 glycerol molecules and 341 water molecules. The configurations of a glycerol molecule and a water molecule are shown on the right of the figure using the ball and stick style. The black balls, gray balls and white balls denote oxygen atoms, carbon atoms and hydrogen atoms, respectively.
Fig. 2. MDS curves of the heat capacity vs. temperature for glycerolaqueous solution (60%, wt/%).
Fig. 3. MDS curves of the density vs. temperature for glycerol-aqueous solution (60%, wt/%).
the slope of the curve is 8.194 · 104 g cm3 K1 (R = 0.9993). The temperature of 164.67 K at the intersection of the two linear curves is taken as the glass transition temperature. Specific volume is mostly used to determine the value of Tg [1,18,16,6]. The curves of the specific volume vs. temperature are shown in Fig. 4. The intersections of the two lines yield the value of Tg (167.51 K), which is 2.84 K higher than 164.67 K.
the rdf peaks implies that the arrangement of atoms and molecules changes gradually from disorder to order with the decrease of the temperature. The first rdf peak ˚ ) denotes the formation probability of the H(r1 = 1.90 A bonds in the glycerol-aqueous solution, which indicates the H-bond strengthens with the decrease of temperature. A plot of P(r) vs. temperature is drawn (Fig. 6). It can be seen from Fig. 6 that the p(r) value for 160 K shows a clearly different change pattern, implying a phase change around 160 K.
Prediction of the Tg by the formation of the H-bonds Verification of the predicted Tg The radial distribution functions (rdfs) between all O atoms and all H on the OH at a temperature range from 90 K to 270 K are shown in the Fig. 5. The pattern of
Gao et al. [11] studied extensively the Tg of glycerol solutions by DSC and obtained the Tg of 60% (wt/%)
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Fig. 4. MDS curves of the specific volume vs. temperature for glycerolaqueous solution (60%, wt/%).
Fig. 5. The radial distribution functions between all O atoms and all H on the OH at various temperatures. The radial distribution functions express the distance distribution between all O atoms and all H on the OH at various temperatures. The rdfs (radial distribution functions) changes with ˚) temperatures (as marked by the arrows). The first peak (r1 1.90A corresponds to the formation probability of the H-bond between all O atoms and all H on the OH. The other peak correspond to the farther distance distribution between all O atoms and all H on the OH. It can be seen that the formation probability of the H-bond between all O atoms and all H on the OH in the solution increases with the decrease of temperature, and simultaneously their bond distances shorten step by step.
glycerol-aqueous solution (Table 1). The Tg of 60% (wt/%) glycerol water solution measured by Luyet and Rasmussen [21] using a DTA is 159 K. As shown in Table 1, the relative errors between the Tg values predicted by MDS and the Tg values in literature [11,21] are all less than 5.5%, which indicates that the NPT-MD method can predict the Tg of glycerol-aqueous solution with relatively high accuracy. Discussion and conclusions Glass transition occurs when a sample goes from a hard, glass like state to a rubber like state. The state parameters
117
Fig. 6. Probability of the first rdf peak vs. temperature for glycerolaqueous solution (60%, wt/%). It is shown that the H-bond becomes stronger with the decrease of temperature. The intersections of the two lines yield the value of Tg (160.06 K), which implies that the microstructures at the temperature above and below Tg are quite different in nature.
of the sample will exhibit an abrupt change within the range of the glass transition temperature. So any instrument that can detect the thermal, mechanical or dielectric changes can theoretically be used to measure the Tg of a material. For example, DSC defines the glass transition as a change in the heat capacity, TMA defines the glass transition in terms of the change in the coefficient of thermal expansion. While Tg can be measured by a variety of methods, it is a time-consuming procedure, especially if the sample is to be kept at subzero temperatures, in anhydrous conditions, or if sampling a portion of the specimen for analysis is cumbersome. Hence, predicting rather than directly measuring Tg as a function of the content of the constituents of a solution can be a powerful tool [17]. Molecular dynamics (MD) is the term used to describe the solution of the classical equation of motion (Newton’s equation) for a set of molecules. In this case, the particles move at constant velocity between perfectly elastic collisions. If the force exerted on each atom is known, it is possible to determine the acceleration of each atom in the system. Integration of the motion equations yields a trajectory that describes the positions, velocities and accelerations of the particles. The average values of properties can be determined from this trajectory. This method is deterministic, that is, once the positions and velocities of each atom are known, the state of the system can be predicted at any time in the future or the past. MD method is different with the Monte Carlo method, which is a statistical technique that involves using random numbers and probability to solve problem. In this study, the heat capacity, the density and specific volume and rdf are used in the software package Discover in Material Studio to predict the Tg of 60% (wt/%) glycerol-aqueous solution. The heat capacity reveals the thermal change when the solution transfers from rubbery state to glass state. While the density and specific volume indicate the mechanical change.
118
D.-X. Li et al. / Cryobiology 56 (2008) 114–119
Table 1 The Tg of glycerol (60% wt/%) aqueous solution determined by DSC From literature Gao et al. Tg (K) Relative errors (%)
Predicted in this paper Luyet et al.
163.60 159.00 Compared to Gao et al. Compared to Luyet et al.
The predicted values were verified with the values from literatures. The relative errors are all less than 5.5%. So the calculation of Tg using this software package is relatively accurate. The deviation may be due to differences in glycerol source, method of measurement or the heating rate in DSC or DTA studies. The differences among the three predicted Tg values result from the difference between the dynamics and thermodynamics properties. The vitrification solutions for cryopreservation usually contain several solutes, which are more complex than the 60 wt% glycerol solution. Theoretically, the molecular dynamics simulation method can be applied to calculate the Tg of realistic complex vitrification solutions if the types and the concentration of the components in the solution are known. But it is estimated that half of one year will be cost to finish the calculation because the cell (similar with Fig. 1) will be much bigger for the realistic vitrification solutions. If the direct calculation method is not practical, another way may be used to predict the Tg of realistic vitrification solutions. Firstly, the Tg of each separate component in the vitrification solutions can be calculated using the MD method. Then the Tg of this vitrification solution can be determined by use of the models as described in Ref. [17]. This will be the future work of our study. In conclusion, molecular dynamic simulation (MDS) is a potential way for investigating the glass transition temperature (Tg) of glycerol–water binary cryoprotectants and may be used for other types of cryoprotectants. Acknowledgments The research was supported by National Science Foundation of China (50576059), The Shanghai Leading Academic Discipline Project (PO502), The Shanghai Development Funds (05ZZ24) as well as the Shanghai University for Science and Technology(X544). References [1] B.F. Abu-Sharkh, Glass transition temperature of poly(vinylchloride) from molecular dynamics simulation: explicit atom model versus rigid CH2 and CHCl groups model, Comput. Theor. Polym. Sci. 11 (2001) 29–34. [2] V.L. Alexey, M.A.J. Michels, Simulation of polymer glasses: from segmental dynamics to bulk mechanics, J. Non-Cryst. Solids 352 (2006) 5008–5012. [3] C.A. Angell, Formation of glasses from liquids and biopolymers, Science 267 (5206) (1995) 1924–1935.
q
Cp
Vspecific
rdf
164.67 0.65 3.57
165.50 1.16 4.09
167.51 2.39 5.35
160.06 3.54 1.06
[4] C.A. Angell, K.L. Ngai, G.B. Mackenna, P.F. McMillan, S.W. Martin, Relaxation in glass forming liquids and amorphous solids, J. Appl. Phys. 88 (2000) 3113–3157. [5] P.G. Debenedetti, F.H. Stillinger, Supercooled liquids and the glass transition, Nature 410 (2001) 259–267. [6] V.I. Dimitrov, Theory of fluidity of liquids, glass transition, and melting, J. Non-Cryst. Solids 352 (2006) 216–231. [7] M.D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, Annu. Rev. Phys. Chem. 51 (2000) 99–128. [8] M.D. Ediger, C.A. Angell, S.R. Nagel, Supercooled liquids and glasses, J. Phys. Chem. 100 (1996) 13200–13212. [9] G.M. Fahy, Vitrification, in: J.J. McGrath, K.R. Diller (Eds.), Low Temperature Biotechnology: Engineering Applications and Engineering Contributions, The American Society of Mechanical Engineers, New York, 1988, pp. 113–146. [10] G.M. Fahy, D.R. Macfarlane, C.A. Angell, H.T. Meryman, Vitrification as an approach to cryopreservation, Cryobiology 21 (1984) 407–426. [11] C. Gao, G.-Y. Zhou, Y. Xu, T.C. Hua, Glass transition and enthalpy relaxation of glycerol and its aqueous solution, Thermochim. Acta 435 (2005) 38–43. [12] N.A. Grunina, T.V. Belopolskaya, G.I. Tsereteli, DSC study the glass transition process in humid biopolymers, J. Phys. Conference Series 40 (2006) 105–110. [13] V.R. Gumen, F.R. Jones, D. Attwood, Prediction of the glass transition temperatures for epoxy resins and blends using group interaction modeling, Polymer 42 (2001) 5717–5725. [14] P. Javier, B. Juan, Glass transition temperature of low molecular weight poly (3-aminopropyl methyl siloxane). A molecular dynamics study, Polymer 43 (2002) 6049–6055. [15] M. Jochem, C.H. Korber, Extended phase diagrams for the ternary solutions H2O–NaCl–glycerol and H2O–NaCl–hydroxyethylstarch determined by DSC, Cryobiology 24 (1987) 513– 536. [16] G.W. Karl, M. Martin, K. Andreas, Z. Gerhard, Glass transition temperature of a cationic polymethacrylate dependent on the plasticizer content—simulation vs. experiment, Chem. Phys. Lett. 406 (2005) 90–94. [17] I. Katkov, F. Levine, Prediction of the glass transition temperature of water solutions: comparison of different models, Cryobiology 49 (2004) 62–82. [18] B. Kurt, B. Jorg, W. Paul, Glass transition of polymer melts: test of theoretical concepts by computer simulation, Prog. Polym. Sci. 28 (2003) 115–172. [19] A.R. Leach, Molecular Modeling—Principles and Applications, second ed., Person Education Limited, Harlow (London), England, 2001, pp. 348–349. [20] B.L. Liu, J. McGrath, L. McCabe, M. Crimp, Response of murine osteoblasts and porous hydroxyapatite (HA) scaffolds to two-step, slow freezing and vitrification process, Cell Preservation Technol. 1 (2002) 33–44. [21] B. Luyet, D. Rasmussen, Study by differential thermal analysis of the temperatures of instability of rapidly cooled solutions of glycerol, ethylene glycol, sucrose and glucose, Biodynamica 10 (1968) 167–191. [22] W. Paul, Molecular dynamics simulations of the glass transition in polymer melts, Polymer 45 (2004) 3901–3905.
D.-X. Li et al. / Cryobiology 56 (2008) 114–119 [23] W.F. Rall, G.M. Fahy, Ice-free cryopreservation of mouse embryos at 196 C by vitrification, Nature 313 (1975) 573– 575. [24] D.R. Reichman, P. Charbonneau, Mode-coupling theory, J. Stat. Mech. P05013 (2005) 1–23. [25] R. Richert, C.A. Angell, Dynamics of glass-forming liquids. V. On the link between molecular dynamics and configurational entropy, J. Chem. Phys. 108 (1998) 9016–9026. [26] F. Sciortino, Potential energy landscape description of supercooled liquids and glasses, J. Stat. Mech. P05015 (2005) 1–35.
119
[27] M. Scott Shell, Pablo G. Debenedetti, Athanassios Z. Panagiotopoulos, A conformal solution theory for the energy landscape and glass transition of mixtures, Fluid Phase Equilibr. 241 (2006) 147–154. [28] F.H. Strillinger, A topographic view of supercooled liquids and glass formation, Science 267 (5206) (1995) 1935–1939. [29] M.J. Taylor, Y.C. Song, B.S. Kheirabadi, F.G. Lightfoot, K.G.M. Brockbank (Eds.), Vitrification fulfills its promise as an approach to reducing freeze-induced injury in a multi-cellular tissue, in: Advances in Heat and Mass Transfer in Biotechnology, vol. 44, Plenum Press, NY, 1999, pp. 93–102.