Molecular dynamics simulations on aqueous two-phase systems - Single PEG-molecules in solution
© Oelmeier et al.; licensee BioMed Central Ltd. 2012
Received: 31 January 2012
Accepted: 25 July 2012
Published: 8 August 2012
Skip to main content
© Oelmeier et al.; licensee BioMed Central Ltd. 2012
Received: 31 January 2012
Accepted: 25 July 2012
Published: 8 August 2012
Molecular Dynamics (MD) simulations are a promising tool to generate molecular understanding of processes related to the purification of proteins. Polyethylene glycols (PEG) of various length are commonly used in the production and purification of proteins. The molecular mechanisms behind PEG driven precipitation, aqueous two-phase formation or the effects of PEGylation are however still poorly understood.
In this paper, we ran MD simulations of single PEG molecules of variable length in explicitly simulated water. The resulting structures are in good agreement with experimentally determined 3D structures of PEG. The increase in surface hydrophobicity of PEG of longer chain length could be explained on an atomic scale. PEG-water interactions as well as aqueous two-phase formation in the presence of PO4 were found to be correlated to PEG surface hydrophobicity.
We were able to show that the taken MD simulation approach is capable of generating both structural data as well as molecule descriptors in agreement with experimental data. Thus, we are confident of having a good in silico representation of PEG.
Polyethylene glycol (PEG) is among the most commonly used chemicals in protein purification processes. Its use ranges from precipitating agent [1, 2], phase forming component in liquid-liquid extraction [3–6], displacer in hydrophobic interaction chromatography  to potential additive in protein formulation [8, 9] and drug modifying agent [10, 11]. Despite its frequent application, the exact nature of the polymer solvent interactions as well as the structural dynamics in solution governing its role in protein purification processes remain unclear.
One approach to understand the nature of molecules on an atomic scale are molecular dynamics (MD) simulations. MD simulations have proven a useful tool to understand the binding and elution behavior of proteins in ion exchange chromatography . Its results could be correlated to experimental data and thus used predictively. Besides HTS process development methodologies, MD simulations represent another promising approach to speed up process development in protein purification.
Several molecular dynamics studies have been performed focusing on the tertiary structure of polyethylene glycol or polyethylene oxide in aqueous solutions. Lee et al. in 2008  used a revision of the CHARMM ether force field (C35)  to establish a molecular understanding of the relation between hydrodynamic radius and radius of gyration of PEG, which differs in behavior from polymer theory under certain conditions. They extended their study in 2009  by using a coarse grained model and were thus able to extend the simulated timespan from 20 to 800 ns. These studies were again focused on investigating the hydrodynamic properties of the polymer molecules, their hydrodynamic radius in relation to their radius of gyration and the coherence of the simulation results with polymer theory. One major goal was to understand - on a molecular scale - the influence of the three dimensional structure of polymer molecules in solution on macroscopic, hydrodynamic properties of these solutions e.g. viscosity. Other studies, such as the one conducted by Borodin et al. in 2001  are geared towards polymer dynamics of pure polymers and the effect of added water. One main goal of this study was again to understand the physicochemical properties of polymer solutions such as viscosity.
Tasaki et al. performed MD simulation of a single PEG molecule of 722 Da in 1996  focused more on the secondary structure and polymer-water interaction. While their study yielded interesting insights into the 3D structure of single PEG molecules in solution, the advance in computational power since their publication now allows for longer simulations of larger PEG molecules. Another, more recent, MD simulation study using a PEG-derivate focused the importance of choosing a proper force field  to obtain secondary structure results in agreement with experimental data. One force field, not tested in their study, but commonly used for protein simulations, is the amber03 force field . While this force field is designated to proteins, a self parameterizing algorithm can be applied to adjust for non proteinaceous molecules . The advantage of using a force field designated to proteins is that combined MD simulations of proteins and PEG could be readily run, once it is established that the self parameterizing algorithm yield reasonable results for PEG. In order to verify the structural data resulting from such MD simulations, experimentally determined structures are needed. There have been several reports of 3D structural data of PEG  with the recent addition of a high quality crystallographic record .
As stated above, MD simulations focused on the secondary structure of single PEG molecules of fixed molecular weight as well as studies concentrating on the three dimensional structure and its implications for macroscopic solution properties have been conducted. Simulations of larger PEG molecules and an investigation of the effects of polymer molecular weight on properties of PEG such as secondary structure, surface hydrophobicity and H-bond network with the solvent are lacking. These properties most likely govern the characteristics of PEG molecules in protein purification, thus their deeper understanding on a molecular scale would be beneficial. PEG polymer molecular weight influences the effect PEG has on protein solubility or aqueous two-phase formation. To understand the influence of molecular weight of PEG molecules on their structural and surface characteristics on an atomic scale would thus be valuable.
This study’s first aim is to confirm the validity of all-atom molecular dynamics simulations of PEG in water based on the amber03 force field by comparing the resulting structure to the recently published crystallographic data and previous MD results. Results from 10-30 ns long MD simulations of PEG molecules ranging in molecular weight from 300 Da to 3500 Da are shown. This study puts a focus on the effect PEG chain length has on geometric parameters, secondary structure, polymer surface characteristics and polymer solvent interactions. We aim to understand the influence on PEG molecular weight on polymer properties important for its use in protein purification. Correlations to experimental data on the phase forming behavior of 4 different PEGs in the presence of phosphate are given. Correlation between solvent polarity determined using a solvatochromic dye to hydrophobicity determined via MD simulations are shown.
Molecular dynamics simulations were performed using the Yasara Structure software package , version 10.2.1. The software was installed in a cluster computer environment running Suse Linux Enterprise 10 as operating system. Each simulation was run using a single node of the cluster computer. Each node was equipped with 2 Intel Xeon Quad Core (X5355) processors at 2.66 GHz and 16GB local memory (2GB per processor core).
The software package used for the simulations employs an automatic parametrization algorithm (termed “AutoSMILES”) to generate force field parameters for unknown structures. The method is described in detail in [19, 20]. This method was used to generate force field parameters for the polymer molecules. Van der Vaals forces were truncated at a cutoff of 10 Å. Long range Coulomb interactions were calculated using the Particle Mesh Ewald algorithm detailed by Essmann et. al. Grid point for the PME evaluation were evenly spaced in each dimension. The amber03 force field  was used for all molecular dynamics simulations. As water model, TIP3P was used. The amber03 force field was chosen as following studies were to include protein molecules in the simulation and amber03 is a well established force field for protein simulations.
The PEG molecule in each simulation snapshot was analyzed for a series of parameters. The following geometric parameters were recorded: dihedral angle between adjacent oxygen atoms (“CC-dihedral”), dihedral angle between adjacent oxygen and carbon atoms (“OC-dihedral”), bond angles (“COC”-, “OCC”-, “CCO”-, and “HCH”-angle), bond length (“OC”-, “CC”-, and “CH”-bond length). Additionally, H-bonds (total number, number donated, number accepted, and total energy), surface characteristics (solvent accessible surface (“SAS”) of CH-groups, O-atoms, and OH-groups), and solvent accessible volume (“SAV”, 1.4Å probe radius) were measured.
Recently published crystallographic data of a PEG molecule  was used to build reference structures. The geometric parameters (dihedral angles, bond angles, bond lengths) of the crystallographic record were measured and PEG-molecules with 6 to 81 subunits constructed accordingly. Surface characteristics and solvent accessible volume of these molecules were measured as described for the MD snapshots.
Data evaluation was performed using Matlab®; (The Mathworks™, Natick, ME, USA), version 2010a. All geometric parameters were first analyzed for a trend over simulation time. If no dependency on simulation time was found, the mean value over all snapshots was used. If no dependency on the position within the molecule was found, the mean value across the molecule was used. Finally, if no dependency on PEG-length was found, the mean value over all PEG-length was calculated. To assess the degree of secondary structure formation of the simulated PEG molecules, regions of four and more consecutive CC-dihedral angles within the “gauche-conformation range” of 50° to 100° or -50° to -100° were considered to form a helical structure. The ratio of CC dihedrals forming a helical structure over the total number of CC dihedrals in the molecule was taken as a measure for the helicality of the molecule. Mean helicality over all snapshots were calculated for each PEG-length simulated. Curve fitting was done using the non-linear-least-square method for fitting and the least-average-residual (“LAR”) method for robustness.
Experimentally determined binodal curves were obtained using the approach previously described . In short, a Tecan (Crailsheim, Germany) Freedom EVO®; 200 robotic platform was used to determine phase transition points using a “cloud point method” modified for liquid handling platforms. PEG used were obtained from Sigma-Adrich (Product numbers were: PEG300: 202371 - PEG600: 202401 - PEG1000: 202428 - PEG1500: 81214 ). 40%[w/w] PO4 stock solution at pH 7.0 was composed of 11.9 g sodiumdihydrogenphosphate and 28.1 g di-potassiumhydrogenphosphate per 100 g of solution.
Determination of empirical E T (30)values were done using the solvatochromic dye nile red (Sigma Aldrich, Steinheim) by measuring the shift of the last absorbance maximum between 500 and 650 nm with a Perkin Elmer lambda 35 spectrophotometer . The following PEGs were investigated: PEG 200, PEG 400, PEG 600, PEG 1000 and PEG 1450 (all purchased from Sigma Aldrich, Steinheim as stated in the previous section). Samples were prepared for different molar fractions of PEG in water with the investigated range depending on the PEG molecular weight (0-44.2% for PEG 200, 0-21% for PEG 400, 0-15.8% for PEG 600, 0-7.8% for PEG 1000 and 0-3.3% for PEG 1450). 2 μl of a saturated solution of nile red in 75% v/v acetonitrile were added to 1ml of sample solution. Calibration of the E T (30) scale was done by plotting E T (30)values of organic solvents (ethanol and different dilutions of ethanol in water, isopropanol, acetone, toluol and octanol) over the wavelength of the absorbance maximum of nile red. E T (30) values were taken from literature . Data was fitted with a 2nd degree polynomial function.
Comparison of geometric parameter between crystal structure of PEG[ and MD simulation results
HCH angle [°]
CCO angle [°]
COC angle [°]
OCC angle [°]
OC distance [Å]
CC distance [Å]
CC dihedral gauche [°]
OC dihedral trans [°]
As mentioned in the previous section, the PEG molecules were found to settle into a random-coil conformation in agreement with previously published studies. While the fluctuation of total system energy suggested to have reached an equilibrium after 1 ns of simulation time, values of R g reached equilibrium only after 25 ns for the longest PEG molecules. In Figure 2 the log R G is plotted over log PE G MW . A linear fit using the least average residual algorithm was calculated and is shown in Figure 2.
The radius of gyration R g was found to correlate with the molecular weight with an exponent v for not significantly different from 0.5. This is in agreement with polymer theory for molecules within the tested range of molecular weights [29, 30] when ideal chain behavior, or a theta solvent respectively, are assumed. While the exponent v was found not to be significantly different from 0.5 the 95% confidence interval accommodates values as high as 0.535 thus also allowing for a divergence towards real-chain behavior. As helical coils were observed as a structural element in the secondary structure, a divergence from the ideal chain model would be reasonable. Longer simulation times might be needed to average out structural fluctuations of R g thus narrowing the confidence interval around v. However, as mentioned above, as the aim of this study was not to infer from structure to macroscopic solution properties, the fluctuations of R G are acceptable within the scope of this work. Solvent accessible volume (“SAV”) of the artificial linear PEG helices and the equivalent average volume of the PEG molecules from MD simulations did not differ significantly. SAVs of simulated PEG molecules were normally distributed around their mean with an average coefficient of variation (“CV”) of 2.4 %. It was concluded, that, while the volumes of PEG molecules in our simulations change dynamically with time, no stable tertiary structures excluding the solvent were formed. While the observed helical regions influenced surface hydrophobicity, they did not significantly change the accessible volume of the molecules.
Figure 1 highlights three aspects of the structural dynamics observed over the course of the simulations. Figure 1-a shows the value of two CC-dihedrals arbitrarily selected from a PEG1162 (n=25) over the entire course of a 5 ns simulation. The two angles are in gauche(-) or gauche(+) conformation in the vast majority of the snapshots. The angles flip irregularly between the (-) and (+) conformation. While the angles are either gauche(-) or gauche(+) most of the time, regions where both angles have the same conformation are more limited. Helical structures form in regions where multiple consecutive CC-dihedrals are of the same conformation. Figure 1-b shows the normalized PEG dihedral energy of three PEGs of different molecular weight over the course of 0.5 ns of simulation. The higher the molecular weight of the PEG, the longer it took to reach equilibrium. Figure 1-c shows the normalized system energy of three PEGs of different molecular weight over the course of 1 ns of simulation. All three simulations show the same trend and reach equilibrium within this time frame. Scattering of the data is lower, the higher the molecular weight of the PEG. Figure 1-d shows the radius of gyration plotted over the simulation time of four exemplary PEG molecules. It can be seen, that the tertiary structure of the PEG reached equilibrium within the simulated timespan. Longer PEG molecules needed longer simulation times. The range over which average properties were calculated is detailed in Figure 2-d. It can be seen, that the tertiary structure is in equilibrium in the time span which was used for the calculation of average PEG properties. Both Figure 1-b and 1-c suggest, that the flexibility of the PEG molecule is dependent on its molecular weight.
It should be noted that due to the compaction of the overall structure, both the hydrophobic and the hydrophilic surface area decrease in the process of helix formation. However, hydrophilic surface area decreases more than the hydrophobic surface area thus making the overall surface more hydrophobic. Second, all three curves show a dependency of the surface hydrophobic fraction on the number of PEG subunits. This effect is mostly independent of the secondary structure. By comparing the hydrophobicity of the end- and mid-groups of the various PEGs, the effect could be ascribed to the “dilution” of the end group with increasing PEG chain length. The end groups are less hydrophobic than the repetitive units in the middle of the molecule. Thus, surface hydrophobicity increases the more of these repetitive units exist. Third, surface hydrophobic fraction of the MD simulations are generally below those of the entirely helical structures. At higher numbers of subunits, the values start to converge. To explain this trend, several underlying effects need to be looked at. As helicality of a PEG molecule has a major effect on its surface hydrophobicity, and helicality increases with number of PEG subunits (see Figure 3-c), it is reasonable that MD simulation results close in on the surface characteristics of a perfectly helical molecule with increasing molecular weight. While helicality increases, helix angle also increases (Figure 3-b) and fraction of CC dihedrals in gauche conformation decreases (Figure 3-d), with both effects decreasing surface hydrophobicity. Additionally, as shown in Figure 4-b, during the MD simulations, temporary structures can be formed that hardly change the solvent accessible volume but strongly reduce the solvent accessible surface of certain parts of the molecules. These regions were found randomly distributed over time and place in the MD simulations and contribute to a decrease in surface hydrophobicity as they only affect mid-groups. The relation of surface hydrophobicity to molecular weight were the result of the combination of all the underlying effects mentioned above.
Starting from a completely linear structure with all geometric parameters set to standard values by the importing plug-in, our molecular simulations resulted in average PEG structures with geometric parameters in good agreement with both the crystal structure  and previous MD simulations .
While Tasaki et al. found a ratio of gauche(+) to gauche(-) of 0.33 to 0.67 for the CC dihedral in their 0.5-2ns long simulations, we found equal distribution of the two conformations. As there is no molecular constraint that would lead to an uneven distribution of the conformation, we conclude that the chosen simulation time of 5 ns is sufficient to yield a representative, average structure. System total energy reaching equilibrium within 1 ns supports this conclusion.
A clear dependency of CC dihedral and the OC distance on PEG chain length was found. The increase in CC dihedral (1.02%) was significantly larger than the increase in OC bond length (0.12%). While the dependency was obvious, we consider the degree of change of OC bond length too low to have a significant effect on the resulting surface characteristics. The reason for this dependency remains unclear.
In a recently published paper, Winger et al. investigated the influence of force-fields on the resulting structure of a modified PEG molecule. They concluded that the development of a helical structure of the simulated PEG molecule in aqueous solution is dependent on the force-field used. Improper force fields or simulations in vacuo resulted in the PEG molecule collapsing into a random coil. In our study, PEG molecules formed helical structures and stayed elongated. The resulting geometric parameters are in good agreement with the results from previous MD simulations  as well as crystallographic data  and structural data from other sources (as summarized in ). We thus conclude, that the amber03 force-field in concert with the employed parametrization algorithm is well suited to run MD simulations of PEG molecules.
Structure dynamics were looked at in terms of time to reach equilibrium of total system energy, scattering of total system energy and time to reach dihedral energy equilibrium. The general conclusion drawn from these observations was, that longer PEG chains are less flexible. It was further concluded, that the increase in helicality with increasing PEG chain length is a direct consequence of the decreased flexibility. Helical regions are less prone to be broken up if the structure is changing less dynamically. Tasaki et al. concluded that their simulation times (0.5-2 ns) were not sufficient to reach an equilibrium. In contrast, system total energy in our simulations reached equilibrium within 1 ns of simulations. We conclude, that simulation times of 1 ns are sufficient when using the simulation protocol employed herein and when the surface of the molecule is the target of the investigation. Reaching an equilibrium tertiary structure, measured as reaching an equilibrium in the radius of gyration, needed significantly longer simulations times as high as 25ns for the largest molecules investigated. Lee et. al in 2009  used a different modeling approach to realize simulation times of up to 800 ns in order to investigate PEG tertiary structure and its hydrodynamic properties in water. Such long simulation times are currently outside the reach of all-atom MD simulations. Again, it should be pointed out that the focus of this study was on the secondary structure and PEG solvent interactions.
In 1970 Koenig et. al concluded from Raman spectroscopy studies that PEG molecules display a more highly ordered structure in water than in methanol. They hypothesized that the helical nature of the solid state is partly retained in the solubilized state . This was also discussed by Devanand et. al and is supported by the results presented herein. While the radius of gyration observed did not support a large deviation from ideal chain behavior, helical regions were found to form within the molecule.
Surface hydrophobicity of the PEG molecules was looked at in terms of surface contribution of CH groups. Two main effects were identified. First, formation of helical regions within the molecule increased the overall hydrophobicity. Second, the effect of the hydrophilic end-group got diluted with increasing PEG chain length. These two effects could explain the overall trend in surface hydrophobicity of the simulated PEG molecules. Additionally, several underlying effects were identified, that influenced surface hydrophobicity and, in part, counteracted one another. For example, both helix angle and overall helicality increased with PEG chain length, with the former decreasing and the latter increasing surface hydrophobicity. While the interplay of the various effects is complex, the overall trend, increasing surface hydrophobicity with increasing PEG chain length showed clearly in our simulations. Zaslavsky et al. measured the hydrophobicity of tritium labeled PEG molecules (1.5 to 40 kDa) in terms of distribution in Ficoll-400-dextran-70 ATPS. The distribution was found to be equal within experimental certainty. The authors concluded that the investigated PEG molecules had equal hydrophobicity. While results from the MD simulations described herein suggest an increase in surface hydrophobicity, the smallest PEG molecule investigated by Zaslavsky et al. had a MW of 1.5 kDa. At this molecular weight (no. of subunits = 34) we found surface hydrophobicity to become less dependent on PEG chain length. Surface hydrophobicity in our simulations changed most in the range 300Da to 1.1kDa. Thus, it would be most interesting to obtain distribution data for PEG molecules within this range.
In this study, the direct interactions of PEG molecules and the surrounding water molecules via H-bonds were quantified. It was found that the number of H-bonds per PEG subunit is a function of PEG chain length and decreases with increasing number of PEG subunits. This is in agreement with the increased hydrophobic surface fraction and the helical structure formed by the PEG molecule, in which the oxygen atoms face inwards and are thus excluded from interactions with the solvent. Tasaki et al. discussed the average number of water molecules associated with the PEG molecule. The association of water molecules with PEG was determined by a water density function around the PEG molecule and found to be 2.9. There are several studies in which water associated with PEG molecules was quantified, the results ranging from 1 to 5 water molecules per PEG subunit. Our results are at the lower end of this range. However, it should be noted, that the interaction via H-bonds quantified here does not necessarily translate directly into the association measured in the before mentioned studies, as the nature of this association remains unclear and is governed by proximity rather than direct interaction.
In general, H-bond formation is an exothermic process. Less formation of H-bond with increasing PEG chain length might create the need for a higher ratio of order water structure around the PEG molecule, which is an entropically unfavorable process. This might have implications on the effect PEG molecules have on other dissolved entities or its behavior on aqueous two-phase systems. A direct correlation of the surface hydrophobic fraction to the experimentally determined PEG concentration needed for two-phase formation in the presence of fixed concentrations of PO4 was found. We concluded that the simulation results are a good representation of the actual behavior of PEG molecules in solution and that the hydrophobic character of the molecule might be a driver of phase separation in PEG ATPSs. Further investigations will be based on this approach focusing on phase formation in PEG-PO4 ATPSs.
While this article focuses on ATPSs, it should be noted that the results might also be applicable to protein precipitation by PEG. Both its structure and thus the space occupied by the PEG molecule in the solvent [1, 34] as well as PEG solvent interactions and thus its influence on solvation energy of the protein might be associated with protein solubility. This will be the scope of following investigations.
To approve the calculation of relative hydrophobicity of PEG based on structural information, we compared E T (30) as an empirical measure for solvent polarity with calculation results based on PEG structure information. In accordance with theory, PEG solutions became more hydrophobic with increasing molar fraction of PEG, indicating that PEG itself is less polar than water (see Figure 7-a). Over a wide range of PEG concentrations this relationships was found to be linear, at higher PEG concentrations further addition of PEG led to an disproportional decrease of E T (30). For that reason the slope for the linear part was used as a measure for PEG polarity and correlated with relative PEG hydrophobicity based surface contributions of unpolar −CH and polar −O− and −OH groups of 1) linear PEG molecules and 2) partly helical PEG molecules from MD simulations. In both cases a linear correlation was obtained with similar R2values (see Figure 7-b), indicating that the structure-based hydrophobicity measure was suitable to describe changes in PEG hydrophobicity. Nevertheless, the experimental data could not be used to differentiate between linear and partly helical PEG structures.
In this study, the validity of using the amber03 force field in combination with the AutoSMILES self parameterizing algorithm for MD simulations of PEG was confirmed. MD simulations were run on a series of PEG molecules ranging in molecular weight from 300 Da to 3500 Da. 3D data from these simulations were found in good agreement with recently published crystallographic data and published MD simulation results. It was found that PEG chain length has a major influence on the surface characteristics and solvent interaction of PEG. Surface hydrophobicity values derived from these simulations could be correlated to experimentally determined minimum PEG concentrations needed to establish two-phase systems in the presence of PO4. Surface hydrophobicity values from MD simulations correlated linearly with solution polarity experimentally determined via a solvatochromic dye. This work forms the basis for conducting further MD studies on phase behavior of PEG as well as simulations including both PEG and proteins such as proteins in ATPSs or simulations of PEGylated proteins.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.