Article Information

• PDF (306 kb)
• Full Text with Large Figures

PubMed

• Articles by Hwankyu Lee
• Articles by Richard W. Pastor

Related Articles

• Hydrodynamic Recruitment of Rolling Leukocytes In Vitro

  Hydrodynamic Recruitment of Rolling Leukocytes In Vitro Biophysical Journal, Volume 84, Issue 6, 1 June 2003, Pages 4182 Michael R. King and Daniel A. Hammer Full Text | PDF (48 kb)

• Electrophoresis of Positioned Nucleosomes

  Electrophoresis of Positioned Nucleosomes Biophysical Journal, Volume 92, Issue 9, 1 May 2007, Pages 3022-3031 Martin Castelnovo and Sébastian Grauwin Abstract We present in this article an original approach to compute the electrophoretic mobility of rigid nucleo-protein complexes like nucleosomes. This model allows us to address theoretically the influence of complex position along DNA, as well as wrapped length of DNA on the electrophoretic mobility of the complex. The predictions of the model are in qualitative agreement with experimental results on mononucleosomes assembled on short DNA fragments (<400 bp). Influences of additional experimental parameters like gel concentration, ionic strength, and effective charges are also discussed in the framework of the model, and are found to be qualitatively consistent with experiments when available. Based on the present model, we propose a simple semi-empirical formula describing positioning of nucleosomes as seen through electrophoresis. Abstract | Full Text | PDF (360 kb)

• Response to the Letter by King and Hammer

  Response to the Letter by King and Hammer Biophysical Journal, Volume 84, Issue 6, 1 June 2003, Pages 4183 Yi Zhang and Sriram Neelamegham Full Text | PDF (46 kb)

• …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 95, Issue 4, 1590-1599, 15 August 2008

doi:10.1529/biophysj.108.133025

Biophysical Theory and Modeling

Molecular Dynamics Studies of Polyethylene Oxide and Polyethylene Glycol: Hydrodynamic Radius and Shape Anisotropy

Hwankyu Lee^*, Richard M. Venable^*, Alexander D. MacKerell^† and Richard W. Pastor^*, ,  

^* Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Bethesda, Maryland 20892
^† Department of Pharmaceutical Sciences, School of Pharmacy, University of Maryland, Baltimore, Maryland 21201

Address reprint requests to Richard W. Pastor.

Abstract

A revision (C35r) to the CHARMM ether force field is shown to reproduce experimentally observed conformational populations of dimethoxyethane. Molecular dynamics simulations of 9, 18, 27, and 36-mers of polyethylene oxide (PEO) and 27-mers of polyethylene glycol (PEG) in water based on C35r yield a persistence length λ=3.7Å, in quantitative agreement with experimentally obtained values of 3.7Å for PEO and 3.8Å for PEG; agreement with experimental values for hydrodynamic radii of comparably sized PEG is also excellent. The exponent υ relating the radius of gyration and molecular weight of PEO from the simulations equals 0.515±0.023, consistent with experimental observations that low molecular weight PEG behaves as an ideal chain. The shape anisotropy of hydrated PEO is 2.59:1.44:1.00. The dimension of the middle length for each of the polymers nearly equals the hydrodynamic radius R_h obtained from diffusion measurements in solution. This explains the correspondence of R_h and R_p, the pore radius of membrane channels: a polymer such as PEG diffuses with its long axis parallel to the membrane channel, and passes through the channel without substantial distortion.

Polyethylene oxide (PEO) and polyethylene glycol (PEG) are polymers with the subunit C-O-C (Fig. 1). They are well-known encapsulating agents for drug delivery 1, solvents for low temperature crystallography 2,3,4, and modulators of osmotic pressure 5,6,7,8. Low molecular weight PEG readily passes through the pores of membrane proteins 9, and, in fact, can be sized by single channels 10. Comparisons with crystal structures and electron micrographs indicate that the pore radius, R_p, is close to the effective hydrodynamic radius in solution, R_h, of the largest PEG able to diffuse through the pore or to block ion conductance 11,12,13. This simple correspondence has allowed estimates of R_p for a variety of toxins 11,12,14, porins 13,15, and other ion channels 16,17.

Display large version of this figure

Figure 1

Structures of DME, PEO, and PEG.

The similarity of R_p and R_h is clearly useful (crystal structures are available for only a small number of membrane proteins), though somewhat surprising. Polymer theory predicts for a random coil polymer in a θ solvent (i.e., an ideal random flight chain) that

(1)

This study, based on a combination of molecular dynamics (MD) simulations of PEO and PEG in solution and hydrodynamic bead model calculations, examines anisotropy of the shape and translational diffusion tensors of these polymers as a step toward understanding their interaction with membrane pores. The study consists of two parts. The first validates a revision (C35r) of the recently developed CHARMM ether force field (C35) 24. This involved MD simulations of solutions of dimethoxyethane (DME; Fig. 1) and water at different mole fractions, and of 9, 18, 27, and 36-mers of PEO, and 27-mers of PEG in water. The populations of the dominant DME conformations for C35 and C35r are compared with Raman data 25. R_g for the polymers for the two parameter sets are then shown to have the molecular weight dependence of an ideal chain, as expected from experimental measurements of short length PEG 26,27. This finding enables a determination of the persistence length, λ, from each parameter set using the worm-like chain model 18, and a comparison with λ obtained from experiments on high molecular weight PEO 28 and PEG 29.

The second part focuses on the diffusion tensors, hydrodynamic radii, and other measures of polymer shape. It is difficult to obtain unambiguous values of the diffusion constant (and hence R_h) from MD simulations of the present systems for two reasons: the statistical error is substantial, and the viscosity of the TIP3P water model 30 must be scaled to compare with experiment 31. Consequently, an alternative estimate of the diffusion constant was obtained using bead model hydrodynamic calculations 32 based on the conformations generated by the MD. The consistency of both approaches lends confidence to the calculated values of R_h, examination of the assumptions inherent in Eq. (1), and a detailed analysis of the anisotropy.

All the simulations and analyses were performed using CHARMM c33b2 33. Ether parameters were from the CHARMM force field (FF) developed by Vorobyov et al. (C35) 24 and the present revision (C35r), and TIP3P model [30,34] was used for water. Using the new velocity Verlet integrator 35 implemented in CHARMM, the temperature was maintained by applying a Nose-Hoover thermostat 36, and the pressure was maintained at 1 atm by applying an Andersen-Hoover barostat 37,38. Electrostatic forces were evaluated using a particle mesh Ewald 39 summation with κ=0.34Å⁻¹, and a real space cutoff of 12Å; the Lennard-Jones forces were switched to zero between 8Å and 12Å, and an isotropic long-range correction 40 was applied. All hydrogen-carbon bond lengths were constrained by SHAKE 41. Simulations were carried out with a time step of 1 fs, and coordinates were saved every picosecond. The parameter and topology files in CHARMM readable format may be obtained from http://mackerell.umaryland.edu/CHARMM_ff_params.html.

Three hundred DME molecules were positioned in a cubic periodic cell of size 38Å/side and energy minimized by 100 steps of steepest descent (SD) and 500 steps of adopted basis Newton Ralphson (ABNR). Mole fractions of 0.6 and 0.3 were constructed by adding 200 and 700 water molecules, respectively, and each system was energy minimized with 100 steps of SD and 500 steps of ABNR. Six 2 ns trajectories (three for each parameter set) were generated at 318K, the temperature of the Raman experiment 25, with conformational averages obtained over the last 1 ns; 2 ns simulations of pure DME were also carried out at 298K to compare with experimental densities, dielectric constant, and heat of vaporization. Gas phase energies at 298K (necessary for calculating the heat of vaporization) were obtained from 2 ns Langevin dynamics simulations at 298K.

A fully extended PEO molecule was used as an initial configuration for 9-mers, and the initial configurations for 18, 27, and 36-mers were randomly obtained from the Langevin dynamics simulations at T=296K. Two simulations with different initial configurations were performed for each molecular size and force field. The simulated system with 9, 18, or 27-mers consists of a PEO or PEG molecule and ∼2900 water molecules in a periodic box of size 44Å/side, and the system with 36-mers molecules consists of a PEO molecule and ∼6500 water molecules in a box of dimension 58Å/side. For each system, the energy was minimized with 30 steps of SD and 100 steps of ABNR. Simulations were performed for 20 ns at 296K, with averages calculated over the last 18 ns.

The persistence length λ was calculated for PEO by a nonlinear least-squares fit 42 of mean-squared end-to-end distance h for the polymer series based on the worm-like chain model 43,44

(2)

(3)

Diffusion constants from the trajectories D_PBC (the subscript denotes periodic boundary conditions) were evaluated from the slopes of the long-time average mean-squared of displacements, MSD(t), of the centers of mass (CM) versus time. As expected for a system consisting of a single solute, MSD(t) was quite irregular at longer times. To improve sampling, the last 18 ns for each simulation was divided into two trajectories of 9 ns, and 4 MSD(t) for each FF, and length of PEO were calculated for 10 ps intervals; these were averaged and a standard deviation σ was calculated for each time point. Displacement of the CM at short times arises from isomerization of dihedral angles, and must be removed when calculating the long-time slope. In this case, the “short-time” cutoff, t₀, was defined as the statistically independent block size determined for the averaging of R_g; t₀=250, 500, 1000, and 2000 ps for the 9, 18, 27, and 36-mers, respectively. Linear regression 42 on the averaged MSD(t) with weighting based on σ[MSD(t)] was then carried out over the range t₀ to 5 ns, and the slopes were divided by 6 to obtain D_PBC for each FF and polymer length. The standard deviations of the slopes evaluated for the individual MSD(t) (over the same time intervals and with the same σ[MSD(t)]) were used to estimate standard errors for each diffusion constant.

The preceding diffusion constants require two adjustments before they can be usefully compared with experiment. The viscosity, η, of pure TIP3P water at 293K equals 0.0035 P 31, which is substantially less than the experimental value of 0.01 P 46. Assuming that the underestimation is the same for the polymer solutions at 296K, where the viscosity of pure water is 0.00933 P, the diffusion constants were scaled by a factor of 0.0035/0.00933=0.375. D_PBC also requires the following finite size correction developed by Yeh and Hummer 47,

(4)

(5)

(6)

(7)

Translational diffusion constants of PEO and PEG were also obtained using the hydrodynamic bead model developed by Garcia de la Torre and Bloomfield 32. In this model, atom centers are considered point sources of friction with hydrodynamic interactions described by a specific shielding tensor, in this case the Oseen tensor 32,43. Inversion of the matrix of interactions and subsequent summing of forces yield a rotational and translational diffusion tensor for the array.

For the present application, friction constants ς for carbons and oxygens of the polymers were calculated from ς=6πηa, where a=0.732Å is the effective bead radius, and η equaled the viscosity of water at 296 K corrected for polymer concentration (Eq. (5)). The preceding value for a is the geometric mean of the C-O and C-C bond distances from the parameter set. It is analogous to the value of 0.77Å (1/2 of the C-C bond length) found appropriate for alkanes 49,50 when calculating diffusion constants using the Kirkwood-Riseman equation 51. It is necessary to explicitly include bound water to the array of beads when modeling diffusion of proteins 52,53, micelles 54,55, and carbohydrates 56,57 in water, so it is expected that some level of bound water is required for solutions of PEO and PEG. Diffusion tensors were evaluated for configurations with no added water, and with waters bound at k_BT, 2k_BT, and 3k_BT generated at 500 ps intervals over the 2–20 ns range for each trajectory. The components of the diffusion tensor were averaged to yield 〈D_xx〉, 〈D_yy〉, and 〈D_zz〉 and average diffusion constant D=(1/3)(〈D_xx〉+〈D_yy〉+〈D_zz〉). Because the radius of gyration for the largest polymers indicated relaxation within 2 ns as an upper limit, sets of four consecutive data points were averaged to obtain independent estimates of the translational diffusion constants for computing the standard errors. Translational diffusion constants evaluated using the bead model without and with bound water are denoted D_bead and .

Raman spectroscopy of DME/water solutions at mole fractions 0.3, 0.6, and 1.0 25 indicates that the principal conformers are TGT, TGG′, TTT, TGG, and TTG, where T and G respectively denote trans, and gauche in the order of C-O-C-C, O-C-C-O and C-C-O-C; GG′ denotes a pair of gauche of opposite signs. As is evident in Table 1 and Fig. 2, trajectories generated by C35 underestimate the populations of TGT and TGG′ and overestimate TTT and TTG. This suggested that stabilization of the gauche state of the central torsion could improve agreement with experiment. The O-C-C-O dihedral potential energy term was refit against the previously calculated gas-phase ab initio surface 24, yielding the revision

(8)

Display large version of this figure

Figure 2

Conformer populations for DME for three mole fractions from experiment (squares), and simulations with C35 (triangles) and C35r (solid circles).

All other terms in C35r are the same as those in C35. Fig. 3 plots selected torsional scans to illustrate the difference between the two parameter sets. Table 1 and Fig. 2 show improved agreement with experiment at all three mole fractions for DME conformers generated with C35r. The molecular volume, dielectric constant, and heat of vaporization also agree well with experiment (Table 1).

Display large version of this figure

Figure 3

Potential energies of O-C-C-O with adjacent C-C-O-C in trans (top); O-C-C-O with C-C-O-C gauche and C-C-O-C trans (middle); C-O-C-C with O-C-C-O and C-O-C-C trans (bottom).

To examine the effect of the FF revision on properties of longer chains, simulations were carried out on 9, 18, 27, and 36-mers of PEO in water and 27-mers of PEG using both C35 and C35r. Fig. 4 (top) plots the instantaneous values of the end-to-end distance and the radius of gyration for PEO36. Although fluctuations are clearly higher for h (Table 2), their autocorrelation functions are similar as expected for these closely related quantities (Fig. 4, bottom). The ∼1 ns decay time justifies the 2 ns block size used to calculate statistical errors, and not averaging over the first 2 ns of the 20 ns trajectories to account for equilibration. Values for 〈h²〉 and R_g from C35r are slightly lower than those with C35, as expected from the increased population of gauche conformers in the O-C-C-O dihedral angle. Entries in Table 2 for PEO27 and PEG27 are statistically equal for each parameter set, indicating that the differences of the terminal groups do not significantly change their properties.

Display large version of this figure

Figure 4

Time series of the end-to-end distance (h) and radius of gyration (R_g) for 36-mers of PEO with C35r (top), and their autocorrelation functions (bottom).

Table 2 Mean-squared end-to-end distance 〈h2〉 and radius of gyration , and their root mean-square fluctuations, σ[h2] and , (in Å2) of PEO and PEG, averaged over the last 18 ns of each simulation

	Force field	Name (Mw)	〈h2〉	σ[h2]
	C35	PEO9	258±16	151	44±1	11
		(442)	289±15	159	46±1	12
		PEO18	542±45	324	91±5	29
		(838)	429±47	316	76±5	29
		PEO27	973±70	515	149±10	46
		(1234)	700±114	494	126±10	46
		PEO36	1316±145	785	204±17	67
		(1630)	1230±259	898	183±26	84
		PEG27	798±100	491	133±9	46
		(1250)	910±195	658	144±21	60
	C35r	PEO9	229±16	147	40±1	11
		(442)	254±16	148	43±1	11
		PEO18	386±40	280	74±4	25
		(838)	418±55	295	79±5	27
		PEO27	765±105	523	125±9	42
		(1234)	754±147	480	129±12	41
		PEO36	1098±295	760	181±65	68
		(1630)	967±198	696	169±21	70
		PEG27	878±156	554	141±24	49
		(1250)	614±107	463	115±11	46

Molecular weights (in daltons) of each polymer are listed in parenthesis below the name.

Fig. 5 plots the molecular weight (M_w) dependence of R_g for PEO, and linear fits yield the coefficient υ in equal to 0.5 within statistical error. Hence, the PEO behaves as an ideal chain 18. This result is consistent with polymer theory. For high molecular weight polymers in “good solvents” (such as water for PEO), mean field and renormalization group treatments of excluded volume interactions yield υ=0.6 and 0.588, respectively 58; υ=0.583 has been experimentally determined for PEO in water for 80,000 < M_w < 10⁶59. However, effects of excluded volume interactions diminish at low M_w, and υ→0.5 58. Solvent systems can also be developed for real polymers for which υ=0.5; these are so-called “θ-conditions”.

Display large version of this figure

Figure 5

log R_g versus log M_w from simulations of PEO.

It is difficult to predict analytically the molecular weight for which nonideality becomes important for a particular polymer, and the following subsection compares hydrodynamic radii obtained from simulation and experiment to demonstrate that υ=0.5 is correct for the present set. For this subsection, the preceding result justifies the use of the worm-like chain model (Eq. (2)) to estimate persistence lengths (λ) from simulation, and to compare them to those obtained experimentally for much higher molecular weight PEO. Nonlinear fitting of the simulated 〈h²〉 to Eq. (2) (see Methods) yields λ=4.3±0.3Å for C35 and 3.7±0.4Å for C35r. The slightly larger value for C35 is consistent with the larger values of 〈h²〉 obtained for each molecular weight.

Mark and Flory 28 obtained a characteristic ratio C_∞=4.1±0.4 for high M_w PEO in aqueous K₂SO₄ at 308K (θ-conditions). Correcting for temperature and applying Eq. (3) yields λ=3.73±0.29Å. More recently, Kienberger et al. 29 measured the persistence length of 3.80±0.02 for PEG in saline solution using atomic force microscopy. Both experimental measurements are in quantitative agreement with the value from C35r.

Chain properties obtained here for PEO are similar, though not identical, to those obtained with the parameter set of Bedrov, Borodin, and Smith (BBS) 60,61. Specifically, extrapolating results for the 11-mer reported by them 60 to the temperature of the present simulations yields ; interpolating between results the PEO-9 and PEO-18 yields for C35 and 47Ų for C35r. To understand this difference, P_OCCO(T) and P_COCC(T), the probabilities of trans in the central and terminal dihedrals of DME, respectively, were evaluated for the three FF. Interpolating from the mole fractions presented by BBS 61, P_OCCO(T) for BBS is between those of C35 and C35r. In contrast, for mole fractions 1, 0.6, and 0.3, P_COCC(T)=0.825, 0.835, and 0.862, respectively, for BBS, whereas P_COCC(T)=0.767–0.768, 0.792–0.794, and 0.815–0.819 for C35 and C35r. This increased P_COCC(T) of DME is the likely reason for the higher radius of gyration for PEO in simulations carried out with the BBS force field.

Columns 3–5 of Table 3 list D_PBC (the translational diffusion constants directly from the trajectories), D_sim (the correction from Eq. (6)), and R_h (the hydrodynamic radius calculated from Eq. (7)). As expected, the diffusion constants decrease and R_h increases with molecular weight, and the statistical error is substantial (as high as 45% for PEO-36). As a check of the validity of these values and for characterizing the anisotropy, diffusion constants were calculated using the hydrodynamic bead model. The last two columns of Table 3 list the diffusion constants based on 40 conformations for each length of PEO and PEG, with and without (D_bead) bound water. Diffusion constants are very similar for C35 and C35r for each case, indicating that the differences in D_sim between the two FF mostly arise from statistical error. This is not surprising. The diffusion constant may be calculated analytically for a sphere in a viscous solution, whereas a diffusion constant calculated from a short simulation of the same system would contain a huge statistical error. In this case, the conformational populations of the polymers are reasonably averaged, so bead model calculations using these conformations (37 for each trajectory) yield standard errors of 7% or less for all systems.

Table 3 Diffusion coefficients (10−6 cm2s−1) from simulations, and the hydrodynamic bead model

			MD simulation			Hydrodynamic bead model
		Name	DPBC	Dsim	Rh	Dbead
	C35	PEO9	5.29±0.97	3.44	6.6±1.2	5.11±0.03	3.91±0.04
		PEO18	4.45±0.97	3.09	7.2±1.6	3.51±0.07	2.85±0.09
		PEO27	2.33±0.15	2.29	9.7±0.6	2.70±0.10	2.24±0.13
		PEO36	1.42±0.37	1.64	13.8±3.6	2.36±0.11	1.98±0.14
		PEG27	3.87±0.66	2.86	7.7±1.3	2.71±0.08	2.22±0.11
	C35r	PEO9	8.88±1.16	4.79	4.8±0.6	5.16±0.03	3.94±0.03
		PEO18	4.12±0.30	2.97	7.5±0.5	3.49±0.05	2.80±0.06
		PEO27	2.21±0.32	2.24	9.9±1.4	2.73±0.08	2.26±0.10
		PEO36	1.64±0.74	1.72	13.1±5.9	2.39±0.12	2.01±0.14
		PEG27	1.57±0.53	2.00	11.0±3.7	2.71±0.08	2.23±0.12

DPBC is obtained directly from the trajectories, and Dsim is corrected for finite size effects and the low viscosity of the TIP3P water model using Eq. (6). Hydrodynamic radii, Rh (in Å), are calculated from Dsim and Eq. (7). Diffusion constants for the bead model are evaluated without bound water, Dbead, and with hydration at 2 kBT level (∼1.5 waters per oxygen), .

D_bead are, on average, 24% higher than D_sim. Although this could be ascribed to an error in scaling the viscosities, as explained in the Methods a hydration shell must be included when calculating diffusion constants of water-soluble polymers. Coordinate sets were therefore generated with increasing levels of bound water based on the polymer-water interaction energies. Cutoffs 3 k_BT (1.1 waters per monomer), 2 k_BT (1.5 waters), and k_BT (3.2 waters) yielded average differences with D_sim of 6.5%, 0.8%, and −15%, respectively. Hence, the 2 k_BT level provides the best match. This level of hydration is consistent with previous ones for carbohydrates 56,57, where ∼1 bound water per hydroxyl group provided optimal agreement with experimental rotational and translational diffusion measurements. The radial distribution functions plotted in Fig. 6 demonstrate that the bound waters are closely associated with the oxygens of the PEO chain. Consequently, the diffusion constants and hydrodynamic radii obtained from the MD simulations are reasonable, and may be used for further analysis.

Display large version of this figure

Figure 6

Radial PEO-water distribution functions from simulations of the 9-mer (solid) and 36-mer (dashed) based on C35r for all neighboring waters (unrestricted) and only those whose binding energy equaled 2 k_BT or greater.

Fig. 7 shows the good agreement of R_h obtained from simulation and experiment 26. A linear fit of the experimental logR_h versus logM_w yields υ=0.506±0.004 for 200<M_w<2000, and υ=0.571±0.003 for 2000<M_w<7500. This further supports the simulation result υ=0.528±0.019 (C35) and 0.515±0.023 (C35r) for 442<M_w<1630 (Fig. 5); i.e., PEO at a low molecular weight is appropriately modeled as a worm-like chain in a θ-solvent.

Display large version of this figure

Figure 7

Hydrodynamic radii (R_h) of PEO versus molecular weight for C35 (diamonds), C35r (squares); experimental values 26 (triangles) are for PEG.

The excellent agreement of simulation and experiment for λ, R_h, and υ indicates that the parameter set C35r yields structures that accurately reflect the conformation of PEO and PEG in solution. The top and middle panels of Fig. 8 compare two such structures of PEO-36 (with waters bound at the 2 k_BT level) with the hydrodynamic radius evaluated from the average of all structures. The bottom panel shows 20 (from 2 ns intervals) in a common alignment. As anticipated from Eq. (1), the length along the long axis (Z) is substantially larger than the hydrodynamic diameter. However, dimensions along the other two axes are clearly comparable to the hydrodynamic diameter, 2 R_h. Table 4 lists the averages for all four PEO. The length along X (the middle axis) is ∼1Å smaller than 2 R_h for each of the four polymers; the length along Z is comparable to 〈h²〉^1/2. Values of the root mean-squared fluctuations are 21%, 17%, and 16% of the lengths along the Z, X, and Y axes, respectively. This result explains the curious relationship between the pore radius R_p of membrane proteins and R_h. Polymers such as PEG in all likelihood diffuse through the pore aligned along their long axes, so the relevant size of the polymer equals 2 R_h to within an angstrom.

Display large version of this figure

Figure 8

Side (left) and end-on (right) views of extended (top) and compact (middle) configurations from a trajectory of PEO36 simulated with C35r. PEO carbons and oxygens are blue and red, respectively; waters interacting at the 2 k_BT level are green spheres; and the hydrodynamic radius of PEO36 evaluated from hydrodynamics of the chain conformations is represented as a transparent sphere positioned in the center of mass. The bottom panel shows the superposition of 20 snapshots at 2 ns intervals, each aligned similarly. The images were created with Visual Molecular Dynamics 67.

Table 4 Average lengths from independent alignments of C35r PEO coordinate sets (including bound water at 2 kBT hydration level) along the three Cartesian axis, assigned in the order Z > X > Y

	Name	Z	X	Y	〈h2〉1/2	2 Rh
	PEO9	18.2±0.1	9.9±0.04	7.1±0.04	15.5±0.7	11.6±0.1
	PEO18	25.3±0.4	14.4±0.1	10.2±0.1	20.0±1.8	15.9±0.3
	PEO27	33.4±0.6	18.9±0.2	12.6±0.1	27.6±3.2	19.9±0.9
	PEO36	39.5±1.1	21.8±0.3	14.8±0.1	32.1±5.4	22.4±1.6

Mean-squared end-to-end distance from Table 2; and hydrodynamic diameters from (Table 3) and Eq. (7). All lengths in angstroms.

The ratios Z/X/Y averaged over the PEO listed in Table 4 are 2.59:1.44:1.00. Further averaging X and Y leads to an aspect ratio of 2.13, in agreement with the value of 2.2 obtained experimentally for DNA 23. The results presented here do not speak to the diffusion constant of PEG or PEO while diffusing through a pore 62,63. Interaction with the pore atoms would be expected to modulate the diffusion. Nevertheless, it is useful to note that the bead calculations yield individual elements of the diffusion tensors, and hence the anisotropy. These are D_zz/D_xx/D_yy=1.18:1.05:1.00. This relatively small diffusional anisotropy is consistent with the shape.

The final topic of this article regards Eq. (1). This equation is derived from the Kirkwood-Riseman (KR) model for polymer diffusion 51, which assumes preaveraging of the Oseen tensor; i.e., the hydrodynamic interactions between beads are evaluated from the average, not instantaneous, distances. This highly useful assumption has been examined by a number of theoretical methods 64,65,66, leading to the conclusion that the exact diffusion constant is ∼10% lower than that predicted by the KR equation. This implies that the coefficient b in R_h=b×R_g should be >0.665. The present bead model calculations were carried out on each conformations (i.e., there is no preaveraging) and so provide another test of the KR equation. A comparison of results with and without hydration also speaks to the applicably of Eq. (1). Fig. 9 plots R_h (with no bound waters and with hydration at the 2 k_BT level) versus R_g (for the polymer only); R_h was taken from the hydrodynamic calculations to reduce scatter. A linear fit to the function R_h=b×R_g yields a slope b=0.70 and 0.85 for unhydrated and hydrated PEO, respectively. As expected, the former slope is slightly higher than 0.665 in Eq. (1). However, that hydration is the more important effect to consider, and may lead to nonlinearity at higher molecular weights.

Display large version of this figure

Figure 9

Hydrodynamic radius versus radius of gyration from simulations based on C35r. The dashed and solid lines represent linear fits of R_h obtained from bead model calculations on hydrated and unhydrated PEO, respectively.

The CHARMM ether force field C35 developed by Vorobyov et al. 24 was revised by adjusting the dihedral potentials to match experimentally measured conformational populations for dimethoxyethane at different mole fractions and the previously published ab initio potential energy surface. MD simulations of 9, 18, 27, and 36-mers of PEO, and 27-mers of PEG based on the revised FF C35r, yield excellent agreement with experiment for persistence lengths and hydrodynamic radii at high and low molecular weights, respectively. In the M_w regime studied, PEO behaves much like an ideal chain; υ in equals 0.5 within statistical error. The properties of the 27-mers for PEO and PEG are statistically equivalent, and conclusions presented are applicable to both polymers.

The shape anisotropy of hydrated PEO simulated here is 2.59:1.44:1.00. Somewhat remarkably, the middle dimension is only ∼1Å larger than twice the hydrodynamic radius R_h. This explains the correspondence of R_h and R_p, the radius of membrane pores: a polymer such as PEG diffuses with its long axis parallel to the membrane channel, so that the effective radius is approximately R_h. It also implies that the polymer does not substantially distort from its solution shape while diffusing through relatively cylindrical membrane channels. Conversely, a breakdown in the relationship R_h≈R_p is expected to occur for irregularly shaped pores.

Lastly, the present simulations and analyses indicate that R_h≈0.85R_g for low molecular PEG and PEO in solution. This increase in the constant of proportionally from 0.665 18 deduced from the Kirkwood-Riseman equation 51 primarily arises from waters of hydration.