|Home | About | Journals | Submit | Contact Us | Français|
A coarse-grained (CG) model for polyethylene oxide (PEO) and polyethylene glycol (PEG) developed within the framework of the MARTINI CG force field (FF) using the distributions of bonds, angles, and dihedrals from the CHARMM all-atom FF is presented. Densities of neat low molecular weight PEO agree with experiment, and the radius of gyration Rg = 19.1 Å±0.7 for 76-mers of PEO (Mw ≈ 3400), in excellent agreement with neutron scattering results for an equal sized PEG. Simulations of 9, 18, 27, 36, 44, 67, 76, 90, 112, 135, and 158-mers of the CG PEO (442 < Mw < 6998) at low concentration in water show the experimentally observed transition from ideal chain to real chain behavior at 1600 < Mw < 2000, in excellent agreement with the dependence of experimentally observed hydrodynamic radii of PEG. Hydrodynamic radii of PEO calculated from diffusion coefficients of the higher Mw PEO also agree well with experiment. Rg calculated from both all-atom and CG simulations of PEO76 at 21 and 148 mg/cm3 are found to be nearly equal. This lack of concentration dependence implies that apparent Rg from scattering experiments at high concentration should not be taken to be the chain dimension. Simulations of PEO grafted to a nonadsorbing surface yield a mushroom to brush transition that is well described by the Alexander-de Gennes formalism.
Polyethylene oxide (PEO) and polyethylene glycol (PEG) are polymers with the formulas H3C-O-(CH2-CH2-O)n-CH3 and HO-(CH2-CH2-O)n-H, respectively. Due to their low toxicity and high solubility in water, they have been conjugated to an array of pharmaceuticals to overcome limitations of low solubility, short circulating life time, and immunogenicity.1,2 Initially, these pharmaceuticals mainly included peptide drugs3,4 and lipid liposomes,5–9 but recently PEGylation has been also applied to oligonucleotides,10–13 saccharides,14 biodegradable hydrogels,15,16 and dendrimers.17–20 A critical issue regarding these assemblies is the extent to which conformations of conjugated and unconjugated PEG or PEO are similar. Good correspondence implies that measurements of free polymers in solution provide insight to the complex (e.g., the height of the PEG layer above the surface of a PEGylated lipid bilayer). In contrast, strong lipid-PEG interactions will lead to deviations with predictions of theories based on nonadsorbing surfaces.21,22 Molecular dynamics (MD) simulation is a natural method to apply to this question. For example, simulations of PEO suggested how the solution hydrodynamic radius Rh is related to the pore radius of membrane proteins.23 However, many important PEG-containing assemblies are too large for all-atom MD simulation. To overcome this limitation, coarse-grained (CG) treatments are required.
CG models for PEO have been developed for implicit solvent, and yield good agreement with experiment for chain dimensions,24–26 and aggregation number and critical micelle concentration.26 Although not computationally demanding, the transferability of implicit solvent models to multicomponent mixtures is limited. Models with explicit solvent are more directly applicable to the interaction of PEG with macromolecular assemblies of lipids and proteins. For example, the CG PEO model with explicit solvent developed by Klein and coworkers showed self-assembly of diblock copolymers in explicit water, and strong interaction with a lipid bilayer.27,28 This paper presents a PEO/PEG model suitable for simulations in the MARTINI CG force field (FF) due to Marrink et al.29,30 The MARTINI CG FF, originally designed for lipids,30 has been extended to proteins31 and successfully applied to large assemblies including proteins, nanoparticles, and membranes.32–37
Particles in the MARTINI FF typically consist of 3 or 4 heavy atoms. For PEG/PEO chains, an intuitive building block defines the sequence C-O-C as a CG bead. Then, as illustrated in Figure 1 for n=2, an n-mer of PEO (denoted PEOn) contains n+1 beads, while an n-mer of PEG has n beads. In this study, CG molecules are consistently named PEOn instead of PEGn+1, though can model for either. The parameterization of the bonded interactions between CG beads is primarily based on all-atom simulations of 9, 18, 27, and 36-mers of PEO (442 < molecular weight Mw < 1,630) in water.23 These all-atom simulations, carried out with the CHARMM ether force field (C35r), show excellent agreement with experiment, as demonstrated for the persistence length λ = 3.7 Å (experimental values are 3.7 Å for PEO38 and 3.8 Å for PEG39), the radius of gyration (ideal chain behavior, as observed40), and the shape anisotropy of 2.59:1.44:1.00 (in reasonable agreement with the axial ratio of 2.2 obtained experimentally for DNA41). Following the strategy used in the development of the MARTINI model for proteins,35 the configurations from the all-atom simulations were mapped to CG pseudo-atoms, and the CG potentials were adjusted to match the bond, angle, and dihedral angle distributions of the pseudo-atoms. Nonbonded interactions between CG beads were parameterized using experimental densities of low molecular weight PEO as targets. End-to-end distributions of the all-atom and CG simulations were compared, but not used explicitly for the parameterization.
The CG model was tested by comparison to the experimentally determined Rg of PEG77 (Mw ≈ 3400), hydrodynamic radii for a range of Mw, and by the transition point of ideal to real chain behavior as measured by Rg. Specifically, theoretical considerations indicate that the coefficient ν in the relation equals 0.5 (ideal chain, no excluded volume) at low molecular weight, and 0.588 (real chain) at high Mw for polymers in good solvents such as water.42 Light scattering of high molecular weight PEO yields ν = 0.583±0.031 in .43,44 An experimental value of ν at low-to-medium Mw can be estimated from hydrodynamic radii Rh obtained from measurements of viscosities and diffusion constants of PEG by Kuga.40 Fitting this data to the relation yields ν = 0.505±0.004 for Mw < 2000 (45 monomers) and 0.571±0.001 for 2000 < Mw < 7500. Assuming that Rh Rg,45 a cross-over from ideal to real chain behavior at Mw ≈ 2000 is a valuable test for the model.
This paper is organized as follows. Sections 2.1 and 2.2 describe the simulation conditions for the coarse-grained and all-atom simulations, respectively, and Section 2.3 outlines the steps required to evaluate Rh. The Results and Discussion contains five subsections. Section 3.1 details the parameterization of CG PEO/PEG, and the comparison with all-atom simulations of low molecular weight PEO. Sections 3.2 and 3.3 compare the molecular weight dependence of Rg and Rh respectively from CG simulations of PEO9 to PEO158 at low concentration. Section 3.4 considers the concentration dependence of Rg for PEO76 based on both all-atom and CG simulations. Experimental measurements46 indicate a large drop in the apparent Rg as the concentration of PEG77 increases from 30 mg/cm3 to 160 mg/cm3. The contributions from interchain scattering to this drop are not known; i.e., the actual Rg of PEG77 may be equal at the preceding concentrations. Hence, the simulations provide another comparison of the CG and all-atom FF, and provide a direct estimate of Rg for the individual chains at high concentration. As an application of the model, Section 3.5 presents the results of simulations of nonadsorbing membrane-like surfaces with different densities of grafted PEO44, and analysis using the Alexander-de Gennes model21,22 of the mushroom to brush transition. This is a preliminary exploration of PEGylated lipid bilayers, with a focus on determining a consistent measure of the height above the surface (the thickness of the brush) for the two regimes. Section 4 presents the conclusions.
Simulations and analyses were performed using the GROMACS simulation package47 with the MARTINI CG force field developed by Marrink et al.29,30 (downloaded from http://md.chem.rug.nl/~marrink/coarsegrain.html). A cutoff of 12 Å was set for Lennard-Jones (LJ) and electrostatic interactions. The LJ potential was smoothly shifted to zero between 9 and 12 Å, and the Coulomb potential was smoothly shifted to zero between 0 and 12 Å. The pressure was maintained at 1 bar and temperature at 296 K (solutions, and PEO grafted to nonadsorbing surface) or 293 K (neat liquids) by the weak-coupling algorithm.48 A time step of 10 fs was used for neat PEO and for the low concentration solutions; an 8 fs time step was employed for the high concentration solutions. Coordinates were saved every 20 ps for analysis. As a final general comment, freezing of pure MARTINI water has been reported at temperatures close to that of the present simulation.29 Freezing has been avoided by applying a 21 % radius increase to 10 % of the water particles.29 No freezing was observed in any of the PEO/water or PEO/lipid/water simulations carried out for this study, though all waters were identical; i.e., the presence of PEO appears to prevent nucleation. The following subsections provide further details for each system.
Neat PEO1, PEO2, PEO3, PEO4, and PEO5 were simulated to calculate densities. 300 PEO molecules were randomly positioned in an initial periodic box of size 51 Å/side, and simulations were performed for 100 ns, with averages calculated over the last 50 ns.
Single chains of PEO were simulated in water for the following lengths and times: PEO9, PEO18, PEO27, PEO36 (400 ns); PEO44, PEO67 (600 ns); PEO76, PEO90, PEO112, PEO135, PEO158 (800 ns). Four trajectories for each molecular size were generated. To avoid interactions between PEO molecules through the periodic boundary conditions, simulation systems with PEO9-PEO36 included ~9,200 water beads (equivalent to ~36,800 real waters) in box sizes of 100 Å/side; and systems with PEO44-PEO158 included ~16,400 water beads (~65,600 real waters) in 125 Å/side boxes. Analyses were performed with the first 20 ns deleted for PEO9 to PEO36, and the first 100 ns deleted for PEO44 to PEO158.
8 copies of PEO76 were randomly positioned in a box of 130 Å/side (18,544 water beads), leading to a concentration of 21 mg/cm3. Similarly, 72 copies of PEO76 in a 140 Å/side box (19,928 water beads) yielded 148 mg/cm3. These are close to ~30 mg/cm3 and ~160 mg/cm3, the concentrations for experimentally measured Rg.46 The simulations were performed for 800 ns, and the last 700 ns were averaged for analysis.
A model nonadsorbing bilayer was constructed as a mixture of modified DSPC (distearoylphosphatidylcholine) and DSPE-PEO44 (PEO44 attached to the ethanolamine bead in the head group of distearoylphosphatidylethanolamine). Glycerol and head group beads of DSPC and lipid parts of DSPE-PEO44 were changed into hydrophobic C1-type beads which are usually designed for the lipid tail group in the MARTINI FF. The lipid-PEO interaction was also modified to prevent adsorption (LJ: ε = 2.0, σ = 6.2) without any change of the interactions of PEO-water and PEO-PEO.
Seven systems at grafting densities ranging from 1 to 100% were equilibrated at constant pressure and temperature (NPT), and then simulated at fixed surface area and constant normal pressure (NPzAT). Numbers of DSPE-PEO44 and total lipids are as follows: 2, 200; 72, 1820; 98, 1568; 200, 1800; 450, 1800; 900, 1800; 1800, 1800. Each system contained approximately 40 CG waters/lipid, and the surface area ranged from 66 (lowest grafting density) to 80 (highest) Å2/lipid. Simulations were performed for 800 ns for the highest grafting density, and 300 for the others; the last 100 ns were used for analysis.
All-atom simulations and analyses were performed using CHARMM c33b2.49 The CHARMM C35r ether force field was used for PEO parameters,23,50 and TIP3P model was used for water.51,52 Using the velocity Verlet integrator53 with a time step of 2 fs, the temperature of 296 K was maintained by the Nose-Hoover thermostat,54,55 and the pressure was maintained at 1 atm by the Andersen-Hoover barostat.56 Electrostatic interactions were calculated using particle mesh Ewald57 and a real space cutoff of 12 Å; van der Waals interactions were switched to zero between 8 Å and 12 Å, and an isotropic long-range correction was applied.58 Coordinates were saved every picosecond for analysis.
A single PEO76 was solvated in a box of 74 Å/side (~13,600 water molecules). Two solvated systems were generated with different configurations of PEO76 having initial Rg of 9 and 24 Å, respectively. Simulations were performed for 70 ns, and the last 50 ns were used for analysis.
5 copies of PEO76 were positioned in a box of 100 Å/side (~35,200 water molecules) for a concentration of 26 mg/cm3; 27 copies of PEO76 were positioned in 100 Å/side (~30,700 water molecules, 146 mg/cm3). These PEO76 concentrations are close to those used in the CG simulation and experiment.46 Initial conditions were generated with the Rg of each PEO ≈ 9, 14 and 27 Å. The experimental Rg for PEG77 at low concentration is 19.7 Å, so the preceding initial conditions span the range from compact to relatively extended. The six simulations (two concentrations, three initial conditions) are denoted “L9”, “L14”, “L27”, “H9”, “H14”, and “H27”, where the number is the initial Rg, and the letter specifies low or high concentration. Trajectories were generated for 24 and 21 ns for L9-L27 and H9-H27, with the first 10 and 15 deleted for analysis, respectively.
The effective hydrodynamic radius Rh is a construct obtained from the Stokes-Einstein relationship42
where D is the translational diffusion constant, T is the absolute temperature, kB is Boltzmann’s constant, and f is the translational friction constant. Stokes’ law for a sphere of radius a with stick boundary conditions in a medium of viscosity η is f = 6πηa.59 Equating the hydrodynamic radius Rh to the Stokes’ radius a, Rh is then obtained from the diffusion constant as
Rh is a useful quantity for comparing models because the simulated diffusion constant Dsim is inversely proportional to the viscosity of the medium. For example, the all-atom low molecular weight PEO discussed here were simulated in TIP3P water,51 where η ≈ 0.32 cP at room temperature.60 This viscosity differs substantially from experiment, and the viscosity of MARTINI CG water. Consequently, D from different simulations and from experiment will likely be very different, whereas Rh, which is both a measure of conformation and interaction with solvent, may be similar.
Calculation of Dsim for each polymer was a multistep process, beginning with DPBC, the diffusion constant obtained in a simulation with periodic boundary conditions. The four individual trajectories for each polymer size were divided in half yielding eight independent segments. The mean square displacement (MSD) of the centers of mass versus time was evaluated for each segment, and a standard deviation was determined for each point of the average MSD. Linear regression61 was then carried out on the averaged MSD with (inverse) weighting equal to the standard deviations over the following ranges: 10 to 100 ns for PEO9,18,27,36; 10 to 150 ns for PEO44,67; 10 to 250 ns for PEO76,158. The slopes were divided by 6, leading to DPBC. Slopes obtained from the 8 individual MSD yielded an estimate of the standard error for DPBC.
DPBC was then corrected for finite size effects using the formula derived by Yeh and Hummer:62
where L is the cubic box length, ξ = 2.837297, and η is the viscosity of the medium. The viscosity increases with higher polymer concentration, and hence the solution viscosity was corrected using the Einstein formula45
where is the volume fraction of the particles, and ηw is the viscosity of pure water (taken to be 0.75 cP from the value obtained from simulations of the present model at 298 K).63 This correction is small, given that ranged from 0.00186 to 0.00346. The final diffusion constant is given by
This formula differs from Eq. (6) of ref. 23, where Dsim was further scaled to the experimental viscosity.
The energy function was designed to be consistent with the MARTINI force field (FF).29 In general, the MARTINI potential consists of bond, angle, Lennard-Jones (LJ), electrostatic, and torsional terms. The following functions were parameterized:
where Kb is the bond force constant, b is the instantaneous bond length, and b0 is the bond length at the minimum energy;
(parameters are defined as for bonds);
where ni and i are the multiplicities and offsets, respectively, of the m individual dihedral terms (m=4 for the present parameterization);
where σij is the zero point of the potential and εij is the well depth at the minimum; εij = 0 for particles bonded to each other.
Parameters for the CG potential energy function (Table 1) were obtained by matching bond, angle, and dihedral distributions, from all-atom simulations23 of PEO9, PEO18, PEO27, and PEO36, and experimental densities of very low molecular weight neat PEO. End-to-end distributions and Rg were used to evaluate the parameters. The following paragraphs describe the critical considerations for refining each of the terms.
In the MARTINI FF, values of σij (in Å) are set to 4.7 for particle types representing approximately 4 heavy atoms, and 4.3 for smaller particles representing 2–3 atoms including those in rings. For particles with σij = 4.7 Å, εij (in kJ/mol) takes values that represent different strengths of interaction. From most attractive to least attractive, the values are: 5.6 (bead types: P5, Qda), 5.0 (P4, P3, Qd, Qa), 4.5 (P2, P1, Nda), 4.0 (Nd, Na), 3.5 (Q0, N0, C5, C4, C3, C2, C1), 3.1, 2.7, 2.3, and 2.0. For interaction between particles with σij = 4.3, the εij are scaled to 75% of those for σij = 4.7, and the bead types are renamed by adding the prefix “S” (e.g., SNda from Nda); εij are not scaled when σij = 4.3 and σij = 4.7 particles interact, and σij is reset to 4.7. Since a monomer of PEO includes 3 heavy atoms (COC), σij = 4.3 and εij = 3 were initially chosen; these are typical parameters for a neutral bead with H-bonding acceptor properties in the MARTINI FF. However, after matching bond, angle, and dihedral distributions to the all-atom results, the Rg and the root mean square end-to-end distances, h21/2, from the CG simulations were too high. This implied that either the PEO-PEO interaction was too weak, or the PEO-water interaction was too strong. Other bead types for PEO were tested, and σij = 4.3 and εij = 3.375 yielded the best results. Thus, the PEO bead is more attractive to itself than expected on the basis of the small molecule parameterization of the standard MARTINI FF. In effect, as far as the PEO-PEO interaction and the interaction with water are concerned, the PEO bead is treated as a SNda type. We anticipate that for interactions with water and other beads the SNa type will be more appropriate.
Parameters for the bond, angle and dihedral potentials were obtained by comparing distributions from our previous all-atom simulations of PEO36,23 after mapping them to the CG representation. Following the conventions for development of the MARTINI FF, distributions were evaluated from the centers of mass for each monomer (C-O-C) in the all-atom model, leading to unimodal distributions for bonds and angles, and a bimodal one for dihedrals (Figure 2, dashed lines). The solid lines in Figure 2, obtained from simulations of PEO36 based on the CG parameters listed in Table 1, agree very well with the all-atom set, demonstrating the success of the parameterization of the bond and angle terms. The standard MARTINI approach does not include proper torsional terms. CG dihedral distributions for alkanes and alkenes can be matched to atomistic distributions sufficiently well using bond and angle potentials in combination with the nonbonded terms that are excluded only between nearest neighbors. For peptides and proteins, elastic bonds have been used to match the specific dihedral distributions for secondary structure elements (extended versus helical arrangements). The preceding strategies failed here, and explicit proper torsional terms were required. Figure 2 shows the very good correspondence of CG and the mapped atomistic dihedral distributions. A similar result for a different CG model of PEO was reported by Fischer et al.25
Nonbonded parameters were optimized on the basis of comparison of densities of low molecular weight PEO to experimental values as well as comparison of the radius of gyration and end-to-end distance distributions of single PEO chains in water to mapped atomistic simulation results. Table 2 compares simulated and experimental densities of low molecular weight PEO. Differences of <5% in PEO3 to PEO5 are comparable to those obtained in the parameterization of CG alkanes.30
Table 3 lists Rg and h21/2 for the single molecule simulations of CG PEO9 to PEO158, all-atom PEO76, and all-atom of PEO9 to PEO36 from reference 23. Figure 3 compares the end-to-end distance distributions Prob (h) at lower molecular weights. Overall agreement of CG and all-atom results is excellent, including PEO76 which was not part of the parameter optimization. Rg for PEO76 equals 19.1 ± 0.7 and 20.4 ± 0.8 Å for the CG and all-atom systems respectively, in excellent agreement with the value of 19.7 Å obtained from neutron scattering of PEG77 at low concentration46 (recall from Figure 1 that PEOn contains the same number of CG beads as PEGn+1). The close agreement of distributions from all-atom PEO27, all-atom PEG28, and CG PEO27/PEG28 (third panel of Figure 3) demonstrates that the CG model can be applied to both PEG and PEO. Figure 3 also includes the analytic result for the worm-like chain model,64
where L is the length of the fully extended (all-atom) polymer. , λ is the persistence length, and α = 3L/4λ. Agreement of the simulated and analytic distributions is very good for the 4 lengths shown for λ = 3.7 Å (the value of the all-atom simulations and experiment).
A minor peak at approximately 5 Å (close to σij) is evident for PEO9 (Fig. 3). This implies that a small subset of the lower molecular weight PEO forms ring-like conformations not observed in the MD simulations or in the worm-like chain model. The peak gradually diminishes by PEO36. Further refinement of the CG model to include a separate bead type for the terminal groups could likely eliminate this conformation. However, given that this artifact is restricted to the low Mw, no changes were made in order to avoid complicating the model.
Figure 4 plots the instantaneous Rg and h of PEO158 over the 800 ns trajectory, and demonstrates the good convergence and stability of the system. When comparisons were possible, the autocorrelation functions of Rg for CG models relaxed more than twice as slowly as those from the all atom models. A close correspondence would not necessarily be expected because the potentials are not the same; however, as established in Section 2.3, the viscosity of CG water is almost 3 times that of TIP3P. Hence, the difference in solvents explains much of the difference in relaxation times.
The present range of PEO allows an examination of the transition from ideal to real chain behavior; i.e., a determination of where the coefficient ν in shifts from ν = 0.5 to the experimental value of 0.583. Figure 5 plots log Rg vs. log Mw. Least squares fits over 442 < Mw < 1630 and 1630 < Mw < 6998 yield ν =0.51 ± 0.02 and 0.57 ± 0.02, respectively. Similarly, those over 442 < Mw < 1982 and 1982 < Mw < 6998 yield ν =0.52 ± 0.01 and 0.57 ± 0.02. Higher break points result in higher ν for the ideal chain regime. Hence, the break in slope (Mw ≈ 1600–2000) occurs near the experimental value Mw = 2000. Neither the value of Rg for PEO76 nor the transition point of ν was explicitly parameterized, and the excellent agreement with experiment thereby provides strong support for the model.
Columns 2 and 3 of Table 4 list DPBC and Dsim for the single particle CG simulations, and those of 8 PEO76 at 21 mg/cm3 (see Section 2.1.c and 3.4). The corrections for periodic boundary conditions are substantial, up to approximately half Dsim. The remaining columns compare Rh for the CG and all-atom models, and experiment. Results from Rh from larger PEO (PEO67, PEO76, and PEO158) agree very well with experiment (average error of 6%). While the trend is captured for PEO9 to PEO44, the simulated values systematically underestimate the experimental values by approximately 20%. Such a limitation might well be expected. Details of the surface become less important in hydrodynamics as size increases. However, the preceding results indicate that caution is required when attempting to model hydrodynamic results for N < 67 for the present model.
At low polymer concentration, excluded-volume effects are predominantly intramolecular, and cause the coil to swell. As concentration increases, intermolecular excluded-volume interactions reduce the size of each coil.46,65 The cross-over point for PEO, however, is not well established. Small angle neutron scattering (SANS) measurments46 on PEG77 indicate a decrease in the apparent Rg from 19.7 Å at 30 mg/cm3 to 8 Å at 160 mg/cm3. The scattering intensity from SANS spectra is only directly related Rg at very low concentration,46,66 so the extent of decrease of Rg, if any, is not clear. Recent all-atom simulations of PEO3 to 30 showed no concentration dependence of Rg,25 although larger PEO were not simulated. To directly compare with experimental values, both all-atom and CG simulations of multiple copies of PEO76 were carried out at low (21–26 mg/cm3) and high concentrations (146–148 mg/cm3).
Figure 6 shows the average Rg of the sets of PEO76 for the all-atom simulations at low and high concentrations (see Section 2.2.b for nomenclature). The average Rg of L9, L14, and L27 are reasonably equilibrated by 10 ns. Similarly, H9, H14, and H27 are equilibrated at 15 ns. The far less costly CG simulations were carried out for 800 ns at both concentrations and the first 100 ns were not included in the final averages. Table 5 lists the results. There is no compelling evidence for a reduction of Rg at the higher concentration; i.e., the apparent Rg obtained from scattering experiments at 160 mg/cm3 likely arise from interchain scattering.
To increase solubility and circulating time of drug molecules, PEG has been attached to the drug transporters such as vesicles, micelles, and nanoparticles,2 which has motivated theoretical and experimental studies of the behavior of PEG grafted on various surfaces.67–71 The theoretical treatment by Alexander21 and de Gennes22 defines two regimes for a polymer grafted on a nonadsorbing surface. At very low grafting density, the chain behaves much like an isolated chain in solution (with the obvious proviso that half of the conformational space is excluded by the surface). Consequently it traces out a hemisphere (or “mushroom”) with a size given by the Flory radius
where N is the degree of polymerization and a is the monomer size; for this study, a is set equal to the bond length b0 = 3.3 Å. RF is proportional to, though not necessarily exactly equal to, h21/2 of the polymer in solution. At high grafting densities, the polymer enters the “brush” regime, with the thickness of the brush L given by
where D is the distance between the grafting points of polymers in the lateral plane. The mushroom and brush states are sketched on the top row of Figure 7. This subsection applies Alexander-de Gennes theory to the results of simulations of PEO44 (N=45) grafted to a bilayer modified to model a hydrophobic nonadsorbing surface. Of particular interest is a consistent specification of L that can be applied to simulations of PEG on more complicated surfaces.
As described in Section 2.1.d, lipid bilayers consisting a mixtures of modified DSPC and DSPE-PEO44 were simulated; Table 6 lists the mole fractions of DSPE-PEO44 and D-values for the 7 systems. The middle and bottom rows of Figure 7 show the side and top-down views, respectively, of the final conformations of the grafted PEO44. It is clear that the polymer is in mushroom and brush-like states at D=44 and 12 Å, respectively. Figure 8 plots the densities of PEO44 normal to the surface for D=9 to 81 Å. The distributions are reasonably similar at high D, and become more extended as D decreases and the chains overlap.
Table 6 lists <h2>1/2 for each value of D. For D =35–81 Å, <h2>1/2 = 28–33 Å; i.e., nearly equal to RF = 32 Å. For D=9–27 Å, <h2>1/2 = 34–43 Å, indicating a departure from the mushroom regime. Alexander-de Gennes theory only defines the brush thickness precisely in the limit a=D from Eq (12); the intermediate values are left unspecified. Here the density profiles were used to define brush thickness. For the mushroom regime, it is reasonable to equate <h2>1/2 to an effective L for a nonadsorbing surface. From the density profiles in this range of D (Figure 8), distances of 28–33 Å correspond to 97–99% of densities. If this same 97–99% density criterion is applied to the polymer at D = 17 Å, L ranges from 45 to 53 Å. These values bracket L = 50 Å obtained for a brush (Eq. 12), and differ substantially from RF = 32 Å. Figure 9 shows the results for all D analyzed in this manner, and compares them with results from Eqns. (11) and (12). There is good agreement between simulation and theory except at the very highest grafting density; this may arise from chain packing interactions in the simulation. The transition from mushroom to brush occurs between D = 27 and 35 Å; i.e., about the value of RF for N=45. An examination of the relation of RF and <h2>1/2 for other values of N, and for adsorbing surfaces will be considered elsewhere.
A coarse-grained (CG) model for polyethylene oxide (PEO) and polyethylene glycol (PEG) was developed within the framework of the MARTINI CG force field. Densities of neat PEO1 to PEO5, and distributions of bond, angle, dihedral, and end-to-end distances of PEO9 to PEO36 in water are in excellent agreement with those from experiments and all-atom simulations. Simulations of PEO9 to PEO158 (442 < Mw < 6998) in explicit water then yielded 1600 < Mw < 2000 for the turnover from ideal to real chains in the relation , and Rg = 19.1 Å for PEO76. Both are in excellent agreement with the experimental targets Mw = 2000 and Rg = 19.7 Å, respectively.
The hydrodynamic radii Rh for the longer chains (>PEO67) are in very good agreement with those from experiment; the lower Mw underestimate experiment by 20% (on average). This observation, along with the small peak at 5 Å in the end-to-end distributions for PEO9 to PEO27, indicates that the model becomes more accurate as molecular weight increases.
Average Rg calculated from CG simulations of multiple copies of PEO76 are 18.9 ± 1.1 Å at low concentration (21 mg/cm3) and 18.6 ± 0.1 Å at high concentration (148 mg/cm3), respectively, close to those calculated from all-atom simulations. These values agree well with the experimental value of 19.7 Å at 30 mg/cm3, but not 8 Å at 160 mg/cm3. Thus, simulations show no concentration dependence of PEO76.
The effective lengths (set to 97–99% of the density above the surface) of PEO44 grafted on hydrophobic (nonadsorbing) surface are described well by Alexander-de Gennes theory. The transition between mushroom and brush occurs at D = 27–35 Å, where D becomes close to the Flory radius RF=32 Å for this length polymer. The agreement of theory and simulation, after an effective brush thickness is defined, provides the basis for analysis of more complex systems including PEGylated bilayers. However, the PEO-lipid interaction parameters must be developed before the present model can be applied to such systems in the MARTINI FF.
We thank Richard Venable, Edward O’Brien, Monica Pate, and P. Thiyagarajan for helpful discussions. This research was supported by the Intramural Research Program of the NIH, National Heart, Lung and Blood Institute. This study utilized the high-performance computational capabilities of the CIT Biowulf/LoBoS3 and NHLBI LOBOS clusters at the National Institutes of Health, Bethesda, MD.