Hybrid Atomistic and Coarse-Grained Molecular Dynamics

Article pubs.acs.org/JPCB

Hybrid Atomistic and Coarse-Grained Molecular Dynamics Simulations of Polyethylene Glycol (PEG) in Explicit Water Francesca Stanzione† and Arthi Jayaraman*,†,‡ †

Department of Chemical and Biomolecular Engineering, University of Delaware, Newark, Delaware 19716, United States Department of Materials Science and Engineering, University of Delaware, Newark, Delaware 19716, United States

‡

S Supporting Information *

ABSTRACT: In-silico design of polymeric biomaterials requires molecular dynamics (MD) simulations that retain essential atomistic/molecular details (e.g., explicit water around the biofunctional macromolecule) while simultaneously achieving large length and time scales pertinent to macroscale function. Such large-scale atomistically detailed macromolecular MD simulations with explicit solvent representation are computationally expensive. One way to overcome this limitation is to use an adaptive resolution scheme (AdResS) in which the explicit solvent molecules dynamically adopt either atomistic or coarse-grained resolution depending on their location (e.g., near or far from the macromolecule) in the system. In this study we present the feasibility and the limitations of AdResS methodology for studying polyethylene glycol (PEG) in adaptive resolution water, for varying PEG length and architecture. We first validate the AdResS methodology for such systems, by comparing PEG and solvent structure with that from all-atom simulations. We elucidate the role of the atomistic zone size and the need for calculating thermodynamic force correction within this AdResS approach to correctly reproduce the structure of PEG and water. Lastly, by varying the PEG length and architecture, we study the hydration of PEG, and the effect of PEG architectures on the structural properties of water. Changing the architecture of PEG from linear to multiarm star, we observe reduction in the solvent accessible surface area of the PEG, and an increase in the order of water molecules in the hydration shells.

I. INTRODUCTION Molecular dynamics (MD) simulation is a well-established computational technique used to study molecular interactions and dynamics in organic and hybrid (inorganic and organic), synthetic, and biological materials.1−4 In terms of computational complexity, atomistically detailed MD simulations of macromolecules (e.g., proteins, polymers) are particularly challenging because of the chemical heterogeneity (e.g., small molecule solvent and macromolecular solute), and the broad range of length and time scales of the molecular driving forces, and the resulting macroscale structure and dynamics seen in these systems. To overcome these limitations of atomistic MD simulations, coarse-grained (CG) approaches have been developed.5 In the CG MD simulations the number of particles/beads and the corresponding degrees of freedom are reduced compared to atomistic MD simulations, but the essential interactions that drive the property of interest are still retained via effective interactions between the CG beads. Such CG approaches enable larger scale simulations, e.g. membrane self-assembly, biopolymer assembly, and protein−protein recognition, than would be impossible with atomistic MD simulation.5−12 However, coarse-graining through implicit representation of the solvent has a major limitation in biologically relevant systems where the solvent (i.e., water) plays a critical role in the behavior of the biomolecules. Water has been shown to affect the structure, dynamics, and functional properties of the biomacromolecules through both © XXXX American Chemical Society

direct hydrogen-bonding and indirect hydrodynamic interactions.13,14 Furthermore, biomacromolecules also influence the local structure and dynamics of the water molecules.15−17 For example, in polyethylene glycol (PEG)-based biomaterials, the high affinity of PEG for water molecules not only increases the solubility of the biomaterial in water, but also creates a hydrated layer around PEG that impacts the interaction of PEG with other biologically relevant proteins.18−20 Therefore, when conducting simulations of biologically relevant macromolecules, water needs to be represented explicitly, and as a result the computational intensity increases because the majority of the simulation time is spent on calculations involving the many water molecules instead of the macromolecule(s) of interest. One way to reduce the computational intensity while representing the solvent explicitly is to conduct simulations at multiple resolutions (e.g., atomistic and coarse-grained) where the molecules either stay fixed in separate domains with different resolutions21−24 or move freely from one resolution to another resolution.12,22,24−26 One way to conduct multiple resolution simulations is to simulate the atoms of the biomacromolecule and the solvent in the region of interest with atomistic resolution and the atoms far away from the region of interest in a CG manner. This is the philosophy Received: March 4, 2016 Revised: April 8, 2016

A

DOI: 10.1021/acs.jpcb.6b02327 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. (a) Schematic representation of AdResS implementation of a linear PEG chain in water. The three regions in this implementation are atomistic (AA, innermost sphere), hybrid (HY, intermediate spherical shell), and coarse-grained (CG, outermost shell) regions. We note that the sizes of oxygen atoms in water and in the PEG are not the same in this schematic for aesthetic reasons, and to visually emphasize the hybrid resolution representation of water. (b) The representations of the molecules in the three regions are assigned based on their positions. In the middle is the transition (hybrid) region HY with the switching function w (curve in black), as described in eqs 1−3 (This schematic is inspired by a similar figure presented in the AdRes work of Fritsch et al. in ref 28.).

underlying the adaptive resolution scheme (AdResS) methodology.26−30 In AdResS, the particle’s resolution depends on its instantaneous location in the predefined atomistic (AA) or coarse-grained (CG) regions. A hybrid (HY) or transition region between AA and CG regions allows for smooth coupling between resolutions and free diffusion across their interface. By maintaining the atomistic details only in the “interesting” regions of the system and reducing the degrees of freedom in the remaining regions of the simulation box, the relevant macromolecule−water and intramacromolecular interactions are still retained with reduction in the computational cost. AdResS methodology has already been used successfully to simulate bulk water27,31,32 or fluids in general26,33,34 and protein and DNA in water.35−37 These previous studies have shown that the AdResS methodology captures the structural properties of the biomolecular system (e.g., folding of proteins) correctly in adaptive resolution water and provides the additional (otherwise inaccessible) insight into the hydration/solvation phenomena.33,38 While the goal of introducing the AdResS method is to reduce the heavy computational costs associated with all-atom explicit solvent simulations, these studies have shown that for small systems the computational costs of AdResS simulations are higher or comparable with that of allatom simulations.33 For larger systems (e.g., proteins or DNA), where the advantage of the adaptive resolution water is fully exploited, a considerable speed up is gained.34,36,37 In this paper we adapt the existing AdResS methodology to investigate polyethylene glycol (PEG) chains of varying polymer lengths and architectures in water to demonstrate the applicability of this AdResS methodology to study polymer based biomaterials in explicit water. In our protocol, an atomistically represented PEG chain with different lengths and architectures is solvated in water, with the water molecules adopting atomistic/coarse-grained resolution depending on their coordinates with respect to the PEG chain. The water molecules in a fixed shell around the PEG chain (termed atomistic zone) is modeled atomistically while the water in the regions outside that fixed shell from the PEG chain is modeled in a CG manner. We compare the structure of the polymer and

water from AdResS approach to that from fully atomistic (allatom) reference simulations to validate the AdResS approach. We then use this AdResS approach to investigate the effect of the PEG length and architecture on PEG hydration. We conclude with a direct comparison of the computational speeds of the two approaches- AdResS and all-atom- for studying the same system. The outline of the paper is as follows. In Section II we describe the AdResS methodology and specific details of the simulations used in this paper. In Section III we present the AdResS simulation results of pure water, linear PEG chains with different lengths (5mer, 10mer, and 20mer), and two different PEG architectures (PEG star and linear PEG), and we compare those to all-atom simulations results. Lastly we compare the computational performance of the AdResS to that of all-atom simulations, followed by the general conclusions in Section IV.

II. METHODS II.A. Adaptive Resolution Simulations. In this paper the AdResS methodology couples atomistic and coarse-grained resolutions for the solvent in a system consisting of a PEG polymer chain in explicit water. The system is divided into predefined regions: an atomistic region (AA), a coarse-grained region (CG), and a hybrid region or transition region (HY) between the AA and CG regions serving to couple these two levels of resolution (Figure 1). To model a fully atomistic PEG polymer in AdResS water, the atomistic region is defined as a sphere of radius dAA centered on a point within the PEG polymer, with dAA greater than the size (e.g., end-end distance) of the PEG polymer (Figure 1). During the simulation, the atomistic region moves with the freely diffusing atomistic PEG. The AA region is surrounded by a spherical shell of width dHY, termed the hybrid (HY) region. The remaining region in the simulation box that is outside the HY shell is the CG region which serves as a particle reservoir for the atomistic region (Figure 1). Water molecules’ representation (AA or CG) depends on its spatial location within these regions, and the B


Article


straightforward in the AdResS setup. In addition, such thermostats do not sample the canonical ensemble. Since the interaction potentials in the CG region are meant to match the radial distribution functions between the center of masses of the molecules in the atomistic simulations (Supporting Information section S.1), the compressibilities of AA and CG remain same40 by construction. However, the pressures in the AA and CG regions can be different due to different AA and CG interaction potentials.41 The inequality in AA and CG pressures can lead to density variations throughout the whole simulation box. Different pressures at the boundaries of the HY region will produce a drift force on the molecules, leading to an unphysical inhomogeneous density profile. An alternative option would be to apply a pressure correction terms to the CG potential during the IBI.42,43 If pressure corrections were implemented in the derivation of the CG potentials, then the pressures of AA and CG regions would match, however CG and AA compressibilities would be different. Different compressibilities would then affect the fluctuations in the number of molecules in the AA and CG regions leading to a barrier to molecule exchange between regions.28 Therefore, to compensate for the pressure generated drift force, while keeping the same compressibilities, we used an external force known as the thermodynamic force (TF), obtained iteratively via an expression based on the gradient of the density difference between AA and CG regions. This force acts on the center of mass of molecules in the HY region and ensures thermodynamic equilibrium between AA and CG regions. Since the water molecules go through the change in resolution in our system, to calculate the thermodynamic force (TF) we conduct a series of short simulations of pure water or systems of interest, where at the i+1 iterative step the thermodynamic force is

water molecules move freely between these AA, CG and HY regions (Figure 1b). In the AA region, the system is described in atomistic detail with standard atomistic force fields for the polymer and water. The CG reservoir of water molecules is described with a suitable CG interaction potential at the appropriate density or pressure obtained from a reference atomistic simulation, through iterative Boltzmann inversion method (IBI) (see Supporting Information section S.1). The two resolutionsAA and CGare coupled via a force-interpolation scheme, in which the intermolecular force between the centers of mass of molecules α and β is given by AA Fαβ = w(rα)w(rβ)Fαβ + [1 − w(rα)w(rβ)]FCG αβ

(1)

where AA Fαβ =

∑ ∑ FijAA i∈α j∈β

(2)

FAA ij is the force between atoms i and j of molecules α and β, respectively, using the atomistic force field, and FCG αβ is the force between the centers of CG α and β molecules. The function w(r) in eq 1 is a transition function varying smoothly and monotonically between 1 in the AA region and 0 in the CG region, and in the HY region it takes the form ⎞ ⎛ π (|rcentre − rα| − dAA)⎟ w(rα) = cos2⎜ ⎠ ⎝ 2dHY

(3)

This way in eq 1, the force between two molecules that are both in the atomistic region or both in the CG region simplifies CG to FAA αβ or Fαβ , respectively. By this choice of forces in the AA− CG interface region, the hybrid resolution molecule experiences forces from molecules in the CG region on a coarse-grained level, and from molecules in the AA region via a combination of atomistic−atomistic and coarse-grained−coarse-grained forces.30 The transition HY region also facilitates the exchange of molecules and ensures equilibrium between the AA and CG regions. An important characteristic of the AdResS approach is that while the overall number of molecules is conserved, the number of the molecules in the individual AA, CG, and HY regions is not. However, a flat density profile throughout the whole simulation box, including the HY region, ensures that the molecules move between resolutions without any (kinetic) barriers. The temperature in the simulation box has to be welldefined and equal throughout the system. While this is simple to achieve in the bulk regions, the transition region (HY) requires some care, since in this region a change in the number of degrees of freedom (DOF) occurs. By switching from CG to atomistic or vice versa, DOFs are reintroduced or removed. Thus, equilibration is accomplished by a local thermostatting procedure30 in the HY region that must be employed together with a system thermostat.39 However, using the Langevin thermostat in the HY and CG regions, we found it necessary to maintain the HY and CG region thermostat at 299 K in order to achieve a temperature in the AA of 300 K. This is due to the fact that particles gain excess heat in the HY region and may not reach the temperature maintained in the HY region, before they cross into the AA region. This artifact does not occur with the use of, e.g., a velocity-rescaling thermostat in the HY region.36 However, thermostats which employ velocity rescaling require the calculation of a global temperature, or at least a temperature over some reasonably large area, which is not

+1 FiTF (r) = FiTF(r) −

M water ρ02 κTAA

∇r ρi (|r|)

(4)

where Mwater is the molecular weight of water, ρ0 is the reference uniform density defined by the NPT ensemble, ρi is the normalized radial density measured from the center of the atomistic zone in the ith iteration, and κAA T is the all-atom isothermal compressibility.30 The thermodynamic force is considered to be converged when the system’s density profile evolves to the target uniform density, i.e. when the density is homogeneous between AA, HY, and CG regions. II.B. Simulation Details. All of the simulations in this paper were carried out using the GROMACS 4.6.744 package, which allows for easy implementation of the AdResS algorithm. To validate the AdResS results, we compare them to the results from purely atomistic (all-atom) simulations. The systems we studied include (a) pure water, (b) a linear PEG chain in water, with varying PEG lengths: 5mer, 10mer, and 20mer (or PEG5, PEG10 and PEG20), and (c) a 4-arm PEG5 star and 8-arm PEG5 star molecule in water. II.B.1. All-atom simulations of pure water. A system composed of 2180 SPC water molecules45 was simulated in a cubic box of dimensions 4 × 4 × 4 nm. The nonbonded interactions between an atom and its neighbors stored in the neighbor list were modeled using Lennard-Jones interactions with a cutoff of 1 nm. The electrostatic interactions between the partially charged atoms in the water molecules were calculated using the reaction fields method46 with a dielectric constant of 80 for water−water interactions.47 Standard C


Article


(see Section II.A for details on CG-HY thermostat effects).36 As done in the all-atom simulations (Section II.B.1), electrostatics were treated using the reaction field method with a dielectric constant of 80 for water−water interactions.46 For all nonbonded Lennard-Jones interactions the cutoff was 1 nm, and simulations were performed in periodic boundary conditions with minimum image convention. The bonds were constrained via the LINCS algorithm. Equilibration runs were 2 ns long, while the production runs had a total length of 40 ns and only the configurations from the last 20 ns were used for analysis. In the AdResS simulations of PEG-water systems, the atomistic (AA) region was centered on the carbon closest to the center of mass of linear PEG 5 mer (PEG5), of linear PEG 10 mer (PEG10) and of linear PEG 20 mer (PEG20), and on the central atom forming the core of the 4-arm and 8-arm PEG5 star. The AA region moves with the freely diffusing PEG while the water molecules change in resolution from AA to CG and vice versa depending on their spatial location in the box. The AA region was sized as dAA= 1.5 nm and the HY region a width of 1.5 nm. For linear PEG5 system, we performed additional simulations to study the behavior of the system as a function of dAA in the range of 1.1−1.5 nm. The thickness of the hybrid zone was kept constant at 1.5 nm in most cases except for the last study where we evaluated the influence of the size of the hybrid zone on the computational efficiency of the simulations. II.C. Analyses. We quantified various properties of both PEG chains and water in all-atom and in AdResS simulations. The structural properties of PEG were evaluated by calculating: (i) Average root mean squared end-end distance (⟨Ree2⟩1/2). This represents the distance between the first and the last heavy atom of the PEG chain. The error bars in the < Ree2>1/2 are the standard deviation calculated from 4000 configurations collected every 5 ps during the last 20 ns of the production run in the simulations. (ii) Radius of gyration of PEG. Rg2, is the second moment around the center of mass of the chain, and is defined as the square of the distance between each atom in the molecule and the center of mass of the molecule. We calculated the average root-mean-square distance denoted as < Rg2>1/5 with the error bars being the standard deviation calculated from 4000 configurations collected every 5 ps during the last 20 ns of the simulations. (iii) Persistence length of PEG chain. The persistence length, P, is calculated by fitting to the equation

periodic boundary conditions and the minimum image convention were implemented. The bonds were constrained via the LINCS algorithm.48 We used a Langevin thermostat39 with a time constant of τ = 2 ps−1 to maintain the temperature at 300 K. Before the production run in the NVT ensemble, the equilibrium density was obtained via 1000 ps of equilibration with a Parrinello−Rahman barostat.49 The production run in the NVT ensemble consisted of 40 ns, and only the last 20 ns were used for the data analyses. During these last 20 ns, the configurations were stored every 2500 steps (5 ps), ensuring negligible statistical correlation between successive stored configurations. We used these atomistic simulations to obtain both the structure and dynamics of atomistic water as well as to obtain the CG potential using IBI techniques, as described in the Supporting Information (Supporting Information section S.1). The water−water radial distribution function from these all-atom simulations and purely coarse-grained simulations using the effective CG potentials are shown in Supporting Information Figure S1. II.B.2. All-atom simulations of linear PEG and star PEG in water. The details of the simulations composed of a single linear PEG chain of varying lengths (5mer, 10mer, and 20mer) in water, and systems composed of a single PEG star (4-arm PEG5 star and 8-arm PEG5 star) are listed in the Table 1. The Table 1. System Details for the All-Atom Simulations of Linear and Star PEG Molecules in Water Model

No. of polymers

No. of water molecules

PEG 5mer (PEG5) PEG 10mer (PEG10) PEG 20mer (PEG20) 4-arm PEG5 star 8-arm PEG5 star

1 1 1 1 1

17116 17102 17156 17128 17021

Box size (in nm) 8 8 8 8 8

× × × × ×

8 8 8 8 8

× × × × ×

8 8 8 8 8

parameters for water simulation were the same as discussed in Section II.B.1. The interaction parameters for PEG atoms were derived from a modified version of the Gromos53a6 force field (Gromos53a6_OXYD),50−53 which has been specifically parametrized to reproduce both the structural and thermodynamic properties of PEG oligomers and polymers in polar and weakly polar environments.51 The simulation protocols used for these systems were the same as that used for the pure water system described in section II.B.1. II.B.3. AdResS simulations. In the adaptive resolution MD simulations, we used a cubic box of length 8 nm in an NVT ensemble at 300 K. The cubic box of size 8 nm was chosen to produce a density determined in an NPT ensemble at 300 K and 1 atm. For the water system, the atomistic (AA) spherical zone was centered on the center of the simulation box with a radius (dAA) of 1.5 nm. The hybrid zone was a spherical shell of thickness (dHY) 1.5 nm. The atomistic water simulation details was described in section II.B.1; the interaction potential of CG water was calculated as described in Supporting Information section S.1. The thermodynamic force was calculated as described in Section II.A. A time step of 0.5 fs was used to capture the fastest time scale in the AA region, namely the vibrations of hydrogen bond. In the AdResS implementation, the particles in the HY and CG regions were coupled to a Langevin thermostat maintained at 299 K using τ = 30 ps−1. This lead to the atomistic region temperature to be 300 ± 1 K

cos θ = e−(L / P)

(5)

where θ is the angle between a vector that is tangent to the polymer chain at position 0 and a tangent vector at a distance L away from position 0, along the contour of the polymer chain. L was defined in integral multiples of the geometric mean bond length of the repeat unit in PEG (=1.464 Å).54 (iv) PEG−water radial distribution f unction. The radial distribution function (or pair correlation function) between the center of mass of the PEG chain (in the case of linear PEG) or PEG arm (in the case of star PEG) and the center of mass of water molecules was calculated to understand their average spatial distribution. For the all-atom simulations we averaged over all the water molecules in the box, whereas for the AdResS simulations only the water molecules in the AA region were considered. D


Article


Figure 2. Water−water radial distribution function from an all-atom simulation (blue line), the atomistic zone of the AdResS simulation (AdResSAA) (black line), and the coarse-grained zone of the AdResS simulation (AdResS-CG) (red line). All-atom and AdResS simulations are performed in NVT at 300 K. The dAA and dHY regions of the AdResS simulation are fixed at 1.5 nm.

sphere in all-atom systems, i.e. the sphere with the same center and radius as the AA region in AdResS. The calculations were limited to subpopulations of water, namely, first hydration shell, second hydration shell and the bulk atomistic water. In those calculations limited to water molecules in the first and second hydration shells, the first hydration shell was defined as those water molecules whose oxygen atom was within a cutoff distance from the nearest PEG heavy atom, the cutoff being taken as the first minimum in the PEG heavy atom-water oxygen RDF. The second hydration shell was defined as the water molecules at 0.35 nm from the first hydration shell. For bulk atomistic water, the water molecules included in the calculation were those within 0.35 nm from the second hydration shell.

For the water molecules we also evaluated in all-atom simulations and in the AA and CG regions of the AdResS simulations, the following quantities: (i) Density of water. To validate the efficacy of the thermodynamic force in compensating the pressure imbalance between AA and CG models we verified that the density is homogeneous across the different regions of the simulation boxes. The water density was evaluated as a function of distance from the center of the AA region. The error bars in the plots are the standard deviation calculated by dividing the last 20 ns of the simulation into 5 blocks of 5 ns time steps and taking the average over these blocks. (ii) Water−water radial distribution f unction. We computed the radial distribution function of the centers of masses of the water molecules. In all-atom simulations we averaged over all the water molecules in the box. To check the local structure, in the AdResS simulations we averaged over the molecules either in the AA or in the CG regions. (iii) Tetrahedral order parameter of water (Sg). To specifically probe the order in the bulk and in the hydration shells, we calculated the tetrahedral order parameter. This order parameter quantifies how well the four nearest neighbors of a water molecule form a tetrahedron. It has a value of 0 in a perfectly tetrahedral arrangement, and a value of 1 in an ideal gas and it is defined as55 Sg =

3 32

3

2 ⎛ 1⎞ ⎜cos ψj , k + ⎟ 3⎠ k=j+1 ⎝

III. RESULTS AND DISCUSSION III.A. Validation of AdResS Results for a Pure Water System. First we validated the AdResS methodology by simulating an adaptive resolution system composed only of water where the water molecules freely diffuse between the AA and CG regions of the box. The AA region was centered on the center of the box and the dAA and dHY were set both to 1.5 nm. In Figure 2 we compare the results from these AdResS simulations to that from the all-atom simulation of pure water, both conducted in an NVT ensemble at 300 K. For the AdResS simulation, two RDFs are shown, one an average over molecules in the AA region and another as an average over molecules in the CG region. The reference all-atom and the AdResS simulation radial distribution functions are in fairly good agreement with only small deviations at the first peak. Furthermore, the agreement in AA and CG regions also shows that the effective potentials in the CG region (Supporting Information Figure S1) calculated through iterative Boltzmann inversion (IBI) are reasonable.

4

∑ ∑ j=1

(6)

where ψj,k is the angle formed by a water oxygen i and its nearest neighbors j and k. Sg (or Sg,w) is calculated considering only the oxygen atoms in the water as potential nearest neighbors. For a fair comparison between AdResS and all-atom systems, analyses were performed on the water molecules in the AA part of the AdResS systems and in the corresponding E


Article


ensure a homogeneous density profile in the different regions of the simulation box (Figure 3b). Without the thermodynamic force we can expect a disparity in densities in the AA, HY and CG regions, as described in later sections. With this agreement in water−water pair correlations between the AdResS (AA and CG zones) and all-atom simulations of the water in NVT ensemble at 300 K, we were justified to use this CG water potential to study infinitely dilute solution of PEG chains in water in NVT ensemble at 300 K using the AdResS methodology. As noted in later sections we recalculated a new thermodynamic forces as needed if we observed inhomogeneity in water density profile in the AA, HY and CG regions using the thermodynamic force of pure water. III.B. AdResS Simulations of Linear and Star PEG Molecules in Adaptive Resolution Water. We calculate the average root mean squared radius of gyration (⟨Rg2⟩1/2), average root mean squared end-end distance (⟨Ree2⟩1/2) and persistence length to show the structural properties of the PEGs in all-atom and AdResS simulations. The structure of PEG in water plays an important role in its application in biomaterials,20,56,57 and being able to reproduce the all-atom structure of PEG in this adaptive resolution water model is important to prove the applicability of the AdResS methodology to study PEG-based biomaterials. First, we present the results for linear PEG 5mer in adaptive resolution water. In Figures 4a−c we present the evolution of these three structural properties of linear 5mer PEG chain (PEG5) during the all-atom and AdResS simulations. The structural properties of linear PEG5 seen in all-atom simulations are retained in the AdResS with adaptive resolution representation of water, and are not affected by a finite size (spherical region of dAA= 1.5 nm) of the atomistically represented water. Furthermore, the PEG-water pair correlation

The details of the IBI method for calculating the potentials used in the coarse-grained model for water are described in Supporting Information section S.1. As stated earlier, the best agreement between the AdResS simulation and all-atom pair correlations would be achieved through the use of the correct thermodynamic force. In Figure 3a we present the thermody-

Figure 3. (a) Thermodynamic force (TF) used and the resulting (b) density profile of water as a function of distance from the center of the atomistic zone in an AdResS simulation composed of 2180 water molecules, performed in an NVT ensemble at 300 K, with the dAA and dHY fixed at 1.5 nm.

namic force that was calculated for this system of water in a cubic box of 8 nm length at 300 K as described in Section II.A. We applied this thermodynamic force in the hybrid zone to

Figure 4. Structural properties of PEG calculated from all-atom (blue lines) and AdResS (red lines) simulations. (a) Average root mean squared end−end distance (⟨Ree2⟩1/2), (b) average root mean squared radius of gyration (⟨Rg2⟩1/2), and (c) the persistence length of PEG. (d) PEG−water radial distribution function from the all-atom simulation (blue line) and atomistic zone of the AdResS simulation (red line). The system is composed of a single linear PEG 5mer immersed in a cubic box of 8 nm in 17116 water molecules. All-atom and AdResS simulations are performed in NVT at 300 K, with the dAA and dHY regions of the AdResS simulation fixed at 1.5 nm. F


Article

The Journal of Physical Chemistry B function from AdResS simulations shows a perfect agreement with that from the all-atom simulation (Figure 4d). To compare the water structure in the bulk, first and second hydration shells for the AdResS and all atom simulations, in Figure 5 we show the probability distribution of the tetrahedral

Figure 6. Water density profile across the simulation box from the allatom simulation (blue line), from the AdResS simulation with the thermodynamic force for pure water (solid red line with solid circles) and from the AdResS simulation using a recalculated thermodynamic force for a linear 5mer PEG chain in water (dashed red line with dashed red circle). The system is composed of a single linear PEG 5mer immersed in a cubic box of 8 nm in 17116 water molecules. Allatom and AdResS simulations are performed in NVT at 300 K, with the dAA and dHY regions of the AdResS simulation fixed at 1.5 nm.

Figure 5. Probability distribution of the tetrahedral order parameter (Sg,w) for water molecules in the first (black line) and second (red line) hydration shells, and in the bulk region (blue line) for a system that is composed of a single linear PEG 5mer (PEG5) immersed in 17116 water molecules placed in a cubic box of 8 nm. All-atom and AdResS simulations are performed in the NVT ensemble at 300 K with the dAA and dHY regions of the AdResS simulation both fixed at 1.5 nm.

thermodynamic force was primarily on the water density profile and was not evident in PEG-water pair correlation function (Supporting Information Figure S2). We expect that the effect of an incorrect thermodynamic force can be more drastic upon decreasing the atomistic zone size (dAA) or hybrid zone thickness (dHY) or upon varying the linear PEG length, PEG architecture, and PEG concentration in water. To test this, we evaluated the validity of the thermodynamic force calculated for a linear 5mer PEG in water with dAA and dHY of 1.5 nm in the NVT ensemble at 300 K, for other closely related PEG systems, i.e. linear PEG chains with varying PEG lengths: 10mer and 20mer (noted as PEG10 and PEG20) at similar dilute concentrations (single chain in ∼17100 water molecules in a cubic box of 8 nm and water density of 1000 kg m−3). We first show the results for PEG10. In Figure 7, we present the density of water across the three regions of the AdResS simulation, and the PEG−water and water−water pair correlation functions. As shown, the thermodynamic force calculated for linear PEG5 works well for linear PEG10 in reproducing the all-atom water density profile (Figure 7a). Furthermore, the AdResS and all-atom PEG−water pair correlation functions show agreement (Figure 7b). Small discrepancies in the water−water pair correlation function (Figure 7c and Supporting Information Figure S3) between allatom and AdResS simulations are observed because the CG water potential used in the AdResS simulations is the effective CG water potential derived from the pure water system, and thus does not capture the effect of PEG on the water−water pair correlation function. As shown in Table 2, the structural properties of PEG10 in adaptive resolution water match the structural properties of PEG10 in water in the all-atom simulation. While the thermodynamic force calculated for PEG5 in water worked reasonably well for PEG10, when we simulated a single linear PEG20 in water with the thermodynamic force calculated for PEG5 in water, we observed that the results from the AdResS did not match all-atom simulations (Supporting Information Figure S4 and Table 3). While the incorrect thermodynamic force is still able to reproduce the same distribution of water across the different regions of the box (Supporting Information Figure S4), Table 3 shows that the

order parameter (Sg,w) for water molecules in the first, second hydration shells and in the bulk region. The all-atom (solid lines) probability distribution of tetrahedral order parameter (Sg,w) in the second hydration shells and the bulk water are reproduced in the AdResS (dashed lines) approach, with only minor deviations in the first hydration shell. Overall, based on these curves we conclude that the tetrahedral packing of water around the linear PEG5, both in the hydration shells as well as in bulk, is captured reasonably well in the AdResS (dashed lines) approach. Thus, the AdResS methodology not only reproduces the structural properties of PEG, but it also captures the solvation of PEG molecule and the hydration shells around the PEG. This means that AdResS is a viable approach for studying systems where the water/solvent arrangement around the macromolecule needs to be treated with chemical detail. We will discuss the solvation of PEG as a function of architecture in Section III.C. We note that all of the above quantities from the AdResS simulations in Figures 4 and 5 were calculated using the thermodynamic force recalculated for linear PEG5 in water. The thermodynamic force calculated for the pure water system and the thermodynamic force calculated for the linear PEG5−water system are shown in Supporting Information Figure S2. The thermodynamic force calculated for the pure water system was not sufficient to reproduce the all-atom results for the system composed of linear PEG5 in water, likely due to the high affinity of PEG to water. The (incorrect) thermodynamic force resulted in slightly inhomogeneous density of the water in the box, as shown in Figure 6. The density profile of water across the simulation box obtained from all-atom (blue line), from AdResS simulations with the thermodynamic force for pure water (solid red circles) and from AdResS simulations using a recalculated thermodynamic force for linear PEG5 in water (dashed red circles). The water density profile from AdResS simulations with a recalculated thermodynamic force better mimics the all-atom water density profile in AA, CG, and HY regions. We note however that this small effect of the incorrect G


Article


Table 3. Structural Properties of PEG20: ⟨Ree2⟩1/2, Persistence Length, and ⟨Rg2⟩1/2 for a System Composed of a Single Linear PEG20 in a Cubic Box of 8 nm in an NVT Ensemble at 300 K (0.17 M)a Simulation approach

⟨Ree2⟩1/2 (nm)

Persistence length (Å)

⟨Rg2⟩1/2 (nm)

All-atom simulation AdResS simulation

2.3 ± 0.94 1.6 ± 0.52

3.38 ± 0.44 3.54 ± 0.43

0.92 ± 0.19 0.78 ± 0.09

The AdResS simulation dAA and dHY are fixed at 1.5 nm, and the thermodynamic force used is the one calculated for PEG5 in water. a

While the incorrect thermodynamic force is likely the cause of the disagreement in < Ree2>1/2 of PEG20 between all-atom and AdResS, we also anticipate the atomistic zone size could be tight for this large linear PEG chain. The < Ree2>1/2 of PEG20 in all-atom simulation is ∼2.3 nm and in the AdResS simulations the PEG20 is confined to atomistic zone diameter of 2dAA = 3 nm where the center of this atomistic zone is centered on the central carbon atom of PEG20. Thus, the atomistic zone size used in the AdResS could constrain the PEG20, leading to smaller sizes of PEG (Table 3) in AdResS than in all-atom simulations. A simulation conducted with a larger dAA (dAA= 2 nm) but with the same dHY (dHY = 1.5 nm) using a correct thermodynamic force derived for linear PEG20 in water (Supporting Information Figure S5) results in a better agreement between the all-atom and the AdResS simulation. Using the structural results of all three PEG oligomer lengths from all-atom and AdResS simulations we also calculated the scaling exponent ν relating the radius of gyration to chain length, Rg ∼ Nν. For all-atom simulations we obtained a scaling exponent of ν = 0.52 ± 0.09, and for the AdResS simulations we obtained a scaling exponent of ν = 0.48 ± 0.07, both of which are consistent with each other and the values of ν obtained previously in other studies conducted on PEG chains in water using a different force field.51,58,59 These data are also consistent with the experimental observation that low molecular weight PEG behaves as an ideal chain in water.60,61 Lastly, we evaluate the validity of the thermodynamic force calculated for linear PEG5 in water in a cubic box of 8 nm in NVT ensemble at 300 K and with dAA and dHY of 1.5 nm, for a 4-arm PEG5 star with 5mer arms (PEG5-star) at similar dilute concentrations (i.e., a single star in a cubic box of 8 nm with a density of water of 1000 kg m−3). The thermodynamic force calculated for a system of a linear PEG5 in water, is incorrect for a system of single PEG5 star with 4 arms of PEG5 in water. As seen in Figure 8a and 8b, the density profile of water seen in the AdResS simulation with the incorrect thermodynamic force does not match that of all-atoms. Furthermore, it also influences the PEG−water pair correlation function. The 4arm star architecture of PEG is effectively at a higher concentration (0.26 M) than a linear PEG5 (0.066 M) and consequently impacts the water, and the PEG−PEG and PEG− water pair correlations. So, we can deduce that the thermodynamic force calculated for linear PEG5 is not sufficient to provide the correct distribution of the water in the atomistic zone in order to properly hydrate the 4-arm PEG5 star. A recalculated thermodynamic force for PEG5-star system is able to correctly reproduce the distribution of water across the box and the PEG-water pair correlation function (Figure 8c and 8d) as seen in all-atom simulations. Additionally the structural properties of the 4-arm PEG5 star agree between all-atom and

Figure 7. (a) Water density profile of the all-atom (blue line) and the AdResS simulation (red line) performed with the thermodynamic force calculated in a PEG5−water system. The system for all parts of this figure is composed of a single linear PEG 10mer (PEG10) immersed in 17102 water molecules in a cubic box of 8 nm. (b) PEG− water radial distribution function (RDF) from the all-atom simulation (blue line), and the atomistic zone of the AdResS simulation (red line). (c) Water−water RDF from the all-atom simulation (blue line), the atomistic zone of the AdResS simulation (AdResS-AA) (black line), and the coarse-grained zone of the AdResS simulation (AdResSCG) (red line). Supporting Information Figure S3 shows the zoomedin versions of the water−water RDF. All-atom and AdResS simulations are performed in the NVT ensemble at 300 K with the dAA and dHY regions of the AdResS simulation both fixed at 1.5 nm.

Table 2. Structural Properties of PEG10: ⟨Ree2⟩1/2, Persistence Length, and ⟨Rg2⟩1/2 for a System Composed of a Single Linear PEG10 in a Cubic Box of Size 8 nm in an NVT Ensemble at 300 K (resulting in a PEG concentration of 0.1 M)a Simulation approach

⟨Ree2⟩1/2 (nm)


All-atom simulation AdResS simulation

1.43 ± 0.51

3.38 ± 0.47

0.615 ± 0.094

1.46 ± 0.52

3.4 ± 0.64

0.62 ± 0.093

a

⟨Rg2⟩1/2

(nm)

For the AdResS simulations, the dAA and dHY are fixed at 1.5 nm.

size of linear PEG20 is underestimated by AdResS simulation when using the thermodynamic force calculated for PEG5 in water. H


Article


Figure 8. (a) Water density profile and (b) PEG−water radial distribution function from all-atom simulation (blue line) and AdResS simulation (red line) for a system composed of a single 4-arm PEG5 star immersed in a cubic box of 8 nm with ∼17100 water molecules, with the AdResS simulation performed using the thermodynamic force calculated for a linear PEG5−water system. (c) Water density profile and (d) PEG−water radial distribution function from the all-atom simulation (blue line) and the atomistic zone of the AdResS simulation (red line) using a recalculated TF for a system composed of a single 4-arm PEG5 star immersed in a cubic box of 8 nm with ∼17100 water molecules. In all cases, the all-atom and AdResS simulations are performed in the NVT ensemble at 300 K with dAA and dHY regions of AdResS simulation fixed at 1.5 nm.

chains in the hydration shell are structurally perturbed by the PEG polymer as a function of PEG chain architecture. Past studies have shown that PEG chains in water62−64 exhibit a “water sponge” effect of PEG, where the PEG monomers form H-bonds with water and effectively retain water molecules. This is in contrast to the hydrophobic effects seen in hydrophobic polymers or proteins with hydrophobic residues in water, where one sees depletion of water from the hydrophobic regions as well as increased density fluctuations of water near the hydrophobic regions.65 To first understand the distance from PEG up to which the explicit water needs to be represented atomistically within the AdResS simulation setup to capture correct (all-atom) PEG hydration, we vary the size of the atomistic zone in the AdResS simulation. In addition to the simulation of linear PEG5 in ∼17100 water molecules with dAA= 1.5 nm, presented in previous sections, we performed AdResS MD simulations of linear PEG5 in the same number of water molecules but with dAA set to 1.25 and 1.1 nm, and same dHY, as depicted with visuals in Supporting Information Figure S6. We reuse the thermodynamic force calculated from PEG5 − water simulation with dAA and dHY set to 1.5 nm. This simulation is performed in NVT ensemble at 300 K and cubic box of 8 nm. As seen in Figure 9a the water density across the AA, HY and CG regions is comparable for all the dAA values. In Figure 9b we show the probability distribution of the tetrahedral order parameter (Sg,w) for water molecules in bulk, and in the first and second hydration shells of PEG5. The distribution of the Sg,w values in the hydration shells and bulk water are similar to the all-atom values for dAA values, with only small shifts observed in the first and second hydration shells for the smallest dAA (1.1 nm).

AdResS simulations with the recalculated thermodynamic force (Table 4). Table 4. Structural Properties of 4-Arm PEG5 Star (PEG5star)a

All-atom simulation AdResS simulation with TF from linear PEG5 in water AdResS simulation with recalculated TF for 4-arm PEG5 star in water

⟨Ree2⟩1/2 (nm)


⟨Rg2⟩1/2 (nm)

1.18 ± 0.38 1.21 ± 0.35

2.9 ± 0.32 2.97 ± 0.32

0.85 ± 0.32 0.8 ± 0.06

1.18 ± 0.33

3.0 ± 0.32

0.82 ± 0.06

a The ⟨Ree2⟩1/2, persistence length, and ⟨Rg2⟩1/2 for a system composed of a single 4-arm PEG5 star (each arm being a linear PEG5) in a cubic box of size 8 nm in an NVT ensemble at 300 K; the AdResS system has dAA and dHY fixed at 1.5 nm.

Based on these data we can deduce that the thermodynamic force can be applied to another system with similar geometry, size, and concentration. However, changing the geometry of the polymer (from linear to star architecture) drastically changes the density of the water across the different (AA, HY, CG) regions of the box, thus warranting a recalculation of the thermodynamic force. Changing the concentration and/or the geometry also impacts the hydration of the polymer in the atomistic zone, and the thermodynamic force has to be recalculated to overcome the imbalance of the water density in the different resolution regions of the box. In the next section we discuss this hydration phenomena in more detail. III.C. Hydration of Linear and Star PEG Polymers. In this section we quantify how the water molecules near the PEG I


Article


Table 6. Structural Properties of Linear PEG5 in All-Atom and AdResS Simulations with Different Atomistic Sizes, dAA = 1.5, 1.25, and 1.1 nm Linear PEG5 All-atom simulation AdResS simulation dAA= 1.5 nm AdResS simulation dAA= 1.25 nm AdResS simulation dAA= 1.1 nm

Furthermore, we show in Table 5 that the average number of water populating the solvation shells is also comparable between all-atom and the AdResS simulations for all dAA. Table 5. Number of Waters in the First Hydration Shell (1SS), in the Second Hydration Shell (2SS), and in the Bulk Region of All-Atoms and AdResS Simulations of Linear PEG5 with Different Sizes of the Atomistic Zone (dAA = 1.5, 1.25, and 1.1 nm) ⟨1SS⟩

⟨2SS⟩

⟨Bulk⟩

All-atom simulation AdResS simulation with dAA = 1.5 nm AdResS simulation with dAA = 1.25 nm AdResS simulation with dAA = 1.1 nm

20 ± 5 19 ± 3 19 ± 4

67 ± 9 64 ± 12 68 ± 10

165 ± 21 160 ± 27 171 ± 24

22 ± 2

72 ± 9

175 ± 20


⟨Rg2⟩1/2 (nm)

1.14 ± 0.37 1.14 ± 0.36

3.35 ± 1.14 3.32 ± 1.07

0.47 ± 0.05 0.46 ± 0.05

1.12 ± 0.38

3.37 ± 1.15

0.46 ± 0.05

1.1 ± 0.32

3.22 ± 0.83

0.46 ± 0.05

simulations are in agreement with that of the all-atom simulations, with largest deviation seen for the smallest dAA. Next, we use AdResS methodology to study the effect of PEG architecture on the hydration of PEG. In Figure 10 we compare the AdResS simulations of PEG with different length and architecture (linear PEG5, 4-arm PEG5 star and 8-arm PEG5 star). The comparison of all-atom and AdResS simulation of 8-arm PEG5 star are shown in Supporting Information Figures S7 and Table S.1. These three systems in Figure 10 effectively have increasing number of PEG5 arms; 1 arm in linear PEG5 case to 8 arms in 8-arm PEG5 star case, and are representative of PEG chains with increasing PEG crowding going from linear PEG5 to 8-arm PEG5 star. Figure 10a shows that changing PEG architecture does not impact water−water pair correlation function as the majority of the water in this calculation being beyond the hydration shells remain unperturbed, and dominate the correlation function to be same as bulk atomistic water. The PEG-water pair correlation decreases as PEG crowding increases going from linear chain (or 1 arm) to 8 arm star (Figure 10b) at PEG-water distances r < 1.75 nm. The solvent accessible surface area per PEG5 arm also shows the decreasing hydration of PEG chains with increasing crowding going from linear chain to 8 arm star (Figure 10c). Figure 10d-f shows that with increasing crowding of PEG monomers going from linear to star architecture there is increasing tetrahedral order of the water molecules, especially in the first and second hydration shells. We expect that increasing PEG monomer crowding is able to confine the water molecules to a more ordered state. III.D. Computational Performance of AdResS versus All-Atom Simulations. One of the major motivations behind conducting a simulation with an adaptive resolution for the explicit solvent molecules is to reduce the heavy computational costs associated with all atom representation of explicit solvent in macromolecular systems. To test if and how the AdResS implementation saves computational expense, we compare the computational performance of all-atom simulations and AdResS simulations conducted in same simulation box sizes on same cores, as described in Table 7. Considering the same simulation box size, the speedup for the AdResS set up over allatom simulations would arise from the coarse-grained (CG) region which have fewer calculations due to fewer degrees of freedom. Thus, one can also expect that the actual speedup of AdResS simulations depends on the ratio of the atomistic and coarse-grained regions For small simulation boxes with a large atomistic zone (dAA = 1.5 nm) the AdResS is slightly slower than all-atom simulation because the additional costs from calculating the forces in the HY region (interpolation of the forces) outweighs the gain in computational speed from the CG zone. For larger CG region achieved via a choice of a smaller

Figure 9. (a) Water density profile of the AdResS simulations with dAA of 1.5 nm (red), dAA of 1.25 nm (blue), and dAA of 1.1 nm (green). (b) Probability distribution of the tetrahedral order parameter (Sg,w) for water molecules in the first (black line) and second (red line) hydration shells and in the bulk region (blue line). The system is composed of a single linear PEG5 immersed in a cubic box of 8 nm with ∼17100 water molecules. In all cases, the all-atom and AdResS simulations are performed in an NVT ensemble at 300 K.

Linear PEG 5mer

⟨Ree2⟩1/2 (nm)

We expect that the lack of differences in hydration patterns with varying dAA will manifest itself in consistent structural properties of the linear PEG5 for varying dAA. Table 6 shows that the ⟨Ree2⟩1/2 for dAA = 1.5 nm matches that of all-atom simulation, and then decreases as dAA decreases further to 1.25 and 1.1 nm. Even though these differences are not statistically significant (within standard deviations), it is important to note the decreasing mean value which could arise likely due to confinement arising from a small atomistic zone. The persistence length and radius of gyration of all AdResS J


Article


Figure 10. (a) Water−water radial distribution function from AdResS simulation of linear PEG5 (black line), 4-arm PEG5 star (red line), and 8-arm PEG5 star (green line). (b) PEG−water radial distribution function from AdResS simulation of linear PEG5 (black line), 4-arm PEG5 star (red line), and 8-arm PEG5 star (green line). (c) Solvent accessible surface area (SASA) per PEG5 arm, and probability distribution of the tetrahedral order parameter (Sg,w) for water molecules in the (d) first and (e) second hydration shells and in the (f) bulk from the AdResS simulation of linear PEG5 (black line), 4-arm PEG5 star (red line), and 8-arm PEG5 star (green line).

IV. CONCLUSION In this paper we adapted the AdResS methodology, which couples coarse-grained and atomistic representations of molecules in a molecular dynamics simulation, to study poly(ethylene glycol) or PEG polymers in explicit water. In our setup, the explicit water molecules underwent a change in resolution from atomistic to coarse-grained (CG) depending on their location with respect to the polymer that is represented atomistically throughout the simulation. In other words, we simulated an atomistic (linear or star) PEG chain with an atomistic hydration shell coupled to a coarse-grained water reservoir. The CG water representation was such that one CG bead represents one water molecule, and the CG interactions were obtained from the iterative Boltzmann inversion (IBI) method to reproduce the all-atom water−water radial distribution function. In the AdResS methodology the hybrid region in between the atomistic and coarse-grained regions ensures a smooth transition of the water molecules from one resolution to the other. We first demonstrated the importance of obtaining a corrective thermodynamic force that ensures homogeneous water density between the atomistic, hybrid, and coarse-grained regions. With the correct thermodynamic force, the AdResS set up produced the same atomistic structure of water and PEG as seen in all-atom simulations. We showed that the thermodynamic force calculated for a system at a state point can be used for systems with similar concentration, and polymer architecture. However, changing the concentration or the geometry of the PEG drastically impacts the hydration of the polymer in the atomistic region and consequently the density of the water across the box, thus needing a recalculation of the thermodynamic force.

atomistic zone (system with dAA = 1.25 or 1.1 nm), we observe a speedup with AdResS simulations over that of all-atom simulations. This speedup is further improved by decreasing the size of the HY zone (dHY) as seen in Table 7. In addition to the time needed for simulation we have to consider the additional time required for the calculation of the thermodynamic force. The calculation of the thermodynamic force requires additional time where the density of the system has to be updated for several steps until the density profile converges with the target (all-atom) density profile. This additional computational cost for the thermodynamic force calculation is a one-time cost for a specific system, and can be reused for multiple trials of the same system. Additionally, as discussed in the previous sections there is also a possibility of using the same thermodynamic force for closely related systems with similar size/architecture and concentration. As discussed in Section III.C the different solvation of the PEG-star architecture warrants recalculation of the thermodynamic force to provide the correct density profile of the water across the box, and a homogeneous distribution of the water molecules. However, as shown in Table 7, the star architecture requires less time for the thermodynamic force calculations due to fewer number of steps needed for convergence of the thermodynamic force calculations than linear PEG. We note that the modest speed-up we have seen with AdResS simulations over all-atom simulation could be further improved by increasing the time step in the CG region and combining it with multiple time step techniques (not implemented yet, but suggested in ref.35) or by using a heavier coarse-graining for the solvent/using an implicit solvent representation in the coarse-grained region.66 K



■

Table 7. Computational Efficiency of All-Atom and AdResS Simulations Discussed in This Studya Simulation

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b02327. Details about iterative the Boltzmann inversion (IBI), method and additional figures, e.g. zoomed in versions of some of the plots in the main paper (PDF)

Total Simulation time as hours per simulation (in brackets are additional hours needed to calculate the TF)

Single 5mer PEG (PEG5) All-atom AdResS dAA = 1.5 nm and dHY = 1.5 nm dAA = 1.25 nm and dHY = 1.5 nm dAA = 1.1 nm and dHY = 1.5 nm dAA = 1.5 nm and dHY = 1 nm dAA = 1.25 nm and dHY = 1 nm Single 10mer PEG (PEG10) All-atom AdResS Single 20mer PEG (PEG20) All-atom AdResS 4-arm PEG5 star (PEG5-star) All-atom AdResS 8-arm PEG5 star All-atom AdResS

Article

13.5

■

13.7 (150) 12.6

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: (302) 831 8682.

10.4

Notes

10.8

The authors declare no competing financial interest.

■

9.7

ACKNOWLEDGMENTS This project was financially supported by the Delaware COBRE pilot program, with a grant from the National Institute of General Medical Sciences−NIGMS (1 P30GM110758-01) from the National Institutes of Health. The authors thank Teragrid/XSEDE (allocation MCB100140) for supercomputing time and the Farber supercomputer supported by the University of Delaware.

15.0 13.7 15.2 13.7 13.9 12.9 (51)

■

15.0 14.8 (24)

REFERENCES

(1) Karplus, M.; McCammon, J. A. Molecular Dynamics Simulations of Biomolecules. Nat. Struct. Biol. 2002, 9 (9), 646−652. (2) Steinhauser, M. O.; Hiermaier, S. A Review of Computational Methods in Materials Science: Examples from Shock-Wave and Polymer Physics. Int. J. Mol. Sci. 2009, 10, 5135−5216. (3) Redondo, A.; LeSar, R. MODELING AND SIMULATION OF BIOMATERIALS. Annu. Rev. Mater. Res. 2004, 34, 279−314. (4) Costa, D.; Garrain, P. A.; Baaden, M. Understanding Small Biomolecule-Biomaterial Interactions: A Review of Fundamental Theoretical and Experimental Approaches for Biomolecule Interactions with Inorganic Surfaces. J. Biomed. Mater. Res., Part A 2013, 101A, 1210−1222. (5) Ingólfsson, H. I.; Lopez, C. a.; Uusitalo, J. J.; de Jong, D. H.; Gopal, S. M.; Periole, X.; Marrink, S. J. The Power of Coarse Graining in Biomolecular Simulations. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (3), 225−248. (6) Clementi, C. Coarse-Grained Models of Protein Folding: Toy Models or Predictive Tools? Curr. Opin. Struct. Biol. 2008, 18, 10−15. (7) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111 (27), 7812− 7824. (8) Riniker, S.; Allison, J. R.; van Gunsteren, W. F. On Developing Coarse-Grained Models for Biomolecular Simulation: A Review. Phys. Chem. Chem. Phys. 2012, 14, 12423. (9) Wang, Q.; Keffer, D. J.; Nicholson, D. M.; Wang, Q.; Keffer, D. J.; Nicholson, D. M. A Coarse-Grained Model for Polyethylene Glycol Polymer A Coarse-Grained Model for Polyethylene Glycol Polymer. J. Chem. Phys. 2011, 135, 214903. (10) Wu, C.; Shea, J.-E. Coarse-Grained Models for Protein Aggregation. Curr. Opin. Struct. Biol. 2011, 21 (2), 209−220. (11) Noid, W. G. Perspective: Coarse-Grained Models for Biomolecular Systems. J. Chem. Phys. 2013, 139, 090901. (12) Kamerlin, S. C. L.; Vicatos, S.; Dryga, A.; Warshel, A. CoarseGrained (multiscale) Simulations in Studies of Biophysical and Chemical Systems. Annu. Rev. Phys. Chem. 2011, 62, 41−64. (13) Chopra, G.; Summa, C. M.; Levitt, M. Solvent Dramatically Affects Protein Structure Refinement. Proc. Natl. Acad. Sci. U. S. A. 2008, 105 (51), 20239−20244.

a The computational efficiency is defined as hours required to perform a single simulation of 40 ns. The additional hours required to calculate the thermodynamic force (TF) are shown in parentheses. The simulations were performed using 40 cores while the TF calculation was performed using 4 cores, all on University of Delaware’s Farber supercomputer, which has 100 compute nodes, each node with 20 Intel “Ivy Bridge” cores, with 6.4TB RAM, 256TB Lustre.

With the validated implementation of the AdResS methodology for PEG in water systems, we then presented the hydration phenomena of PEG−water systems as a function of PEG architecture. By decreasing the size of the atomistic zone from dAA = 1.5 to 1.1 nm, we determined the minimum possible size of the atomistic region in such a simulation setup where the properties of the PEG and water are still preserved and comparable with the all-atom simulations. We showed that changing the PEG architecture from linear to star increases the local PEG monomer crowding, and as a result decreases the hydration of the PEG monomers. In addition we showed that the tetrahedral order of the water around the PEG increases with increasing PEG monomer crowding. We concluded the paper with a comparison of the computational expense of the AdResS and all-atom simulations. For the systems studied here, we did not see a speed-up of the AdResS simulations over the all-atom simulations when the atomistic and hybrid zones were comparable in volume to the coarse-grained zone (dAA and dHY of 1.5 and 1.5 nm, respectively, in a 8 nm cubic box). However, for small atomistic and hybrid zones (dAA and dHY of 1.1 and 1 nm, respectively), the AdResS simulations provide a consistent speed up compared to the all-atom reference simulation. The thermodynamic force calculation was an additional computational cost that was incurred when doing the AdResS simulations that is not needed in all-atom simulations. L


Article

The Journal of Physical Chemistry B (14) Vogler, E. A. Role of Water in Biomaterials. Biomaterials Science: An Introduction to Materials in Medicine; The Pennsylvania State University: 2012. (15) Agarwal, M.; Kushwaha, H. R.; Chakravarty, C. Local Order, Energy, and Mobility of Water Molecules in the Hydration Shell of Small Peptides. J. Phys. Chem. B 2010, 114 (1), 651−659. (16) Bhattacharjee, N.; Biswas, P. Structure of Hydration Water in Proteins: A Comparison of Molecular Dynamics Simulations and Database Analysis. Biophys. Chem. 2011, 158 (1), 73−80. (17) Henchman, R. H.; McCammon, J. A. Extracting Hydration Sites around Proteins from Explicit Water Simulations. J. Comput. Chem. 2002, 23 (9), 861−869. (18) Zhang, Z.; Berns, A. E.; Willbold, S.; Buitenhuis, J. Synthesis of Poly(ethylene Glycol) (PEG)-Grafted Colloidal Silica Particles with Improved Stability in Aqueous Solvents. J. Colloid Interface Sci. 2007, 310 (2), 446−455. (19) Mu, Q.; Hu, T.; Yu, J. PLoS One 2013, 8 (7), 1−10. (20) Park, K. D.; Kim, Y. S.; Han, D. K.; Kim, Y. H.; Lee, E. H. B.; Suh, H.; Choi, K. S. Bacterial Adhesion on PEG Modified Polyurethane Surfaces. Biomaterials 1998, 19 (7−9), 851−859. (21) Sokkar, P.; Choi, S. M.; Rhee, Y. M. Simple Method for Simulating the Mixture of Atomistic and Coarse-Grained Molecular Systems. J. Chem. Theory Comput. 2013, 9 (8), 3728−3739. (22) Rzepiela, A. J.; Louhivuori, M.; Peter, C.; Marrink, S. J. Hybrid Simulations: Combining Atomistic and Coarse-Grained Force Fields Using Virtual Sites. Phys. Chem. Chem. Phys. 2011, 13 (22), 10437− 10448. (23) Wassenaar, T. A.; Ingólfsson, H. I.; Prieß, M.; Marrink, S. J.; Schäfer, L. V. Mixing MARTINI: Electrostatic Coupling in Hybrid Atomistic-Coarse-Grained Biomolecular Simulations. J. Phys. Chem. B 2013, 117 (13), 3516−3530. (24) Nielsen, S. O.; Bulo, R. E.; Moore, P. B.; Ensing, B. Recent Progress in Adaptive Multiscale Molecular Dynamics Simulations of Soft Matter. Phys. Chem. Chem. Phys. 2010, 12 (39), 12401−12414. (25) Praprotnik, M.; Delle Site, L.; Kremer, K. A Macromolecule in a Solvent: Adaptive Resolution Molecular Dynamics Simulation. J. Chem. Phys. 2007, 126 (13), 134902. (26) Potestio, R.; Fritsch, S.; Español, P.; Delgado-Buscalioni, R.; Kremer, K.; Everaers, R.; Donadio, D. Hamiltonian Adaptive Resolution Simulation for Molecular Liquids. Phys. Rev. Lett. 2013, DOI: 10.1103/PhysRevLett.110.108301. (27) Matysiak, S.; Clementi, C.; Praprotnik, M.; Kremer, K.; Delle Site, L. Modeling Diffusive Dynamics in Adaptive Resolution Simulation of Liquid Water. J. Chem. Phys. 2008, 128 (2), 024503. (28) Fritsch, S.; Poblete, S.; Junghans, C.; Ciccotti, G.; Site, L. D.; Kremer, K. Adaptive Resolution Molecular Dynamics Simulation through Coupling to an Internal Particle Reservoir. Phys. Rev. Lett. 2012, 108 (April), 170602. (29) Praprotnik, M.; Delle Site, L.; Kremer, K. Adaptive Resolution Molecular-Dynamics Simulation: Changing the Degrees of Freedom on the Fly. J. Chem. Phys. 2005, 123 (22), 224106. (30) Wang, H.; Schütte, C.; Delle Site, L. Adaptive Resolution Simulation (AdResS): A Smooth Thermodynamic and Structural Transition from Atomistic to Coarse Grained Resolution and Vice Versa in a Grand Canonical Fashion. J. Chem. Theory Comput. 2012, 8 (8), 2878−2887. (31) Praprotnik, M. Adaptive Resolution Simulation of Liquid Water. J. Phys.: Condensed Matter 2007, 19 (29)292201. (32) Zavadlav, J.; Melo, M. N.; Cunha, A. V.; De Vries, A. H.; Marrink, S. J.; Praprotnik, M. Adaptive Resolution Simulation of MARTINI Solvents. J. Chem. Theory Comput. 2014, 10 (6), 2591− 2598. (33) Fritsch, S.; Junghans, C.; Kremer, K. Structure Formation of Toluene around C60: Implementation of the Adaptive Resolution Scheme (AdResS) into GROMACS. J. Chem. Theory Comput. 2012, 8 (2), 398. (34) Bevc, S.; Junghans, C.; Kremer, K.; Praprotnik, M. Adaptive Resolution Simulation of Salt Solutions. New J. Phys. 2013, 15 (10), 105007.

(35) Zavadlav, J.; Melo, M. N.; Marrink, S. J.; Praprotnik, M. Adaptive Resolution Simulation of an Atomistic Protein in MARTINI Water. J. Chem. Phys. 2014, 140 (5), 054114. (36) Fogarty, A. C.; Potestio, R.; Kremer, K. Adaptive Resolution Simulation of a Biomolecule and Its Hydration Shell: Structural and Dynamical Properties. J. Chem. Phys. 2015, 142 (19), 195101. (37) Zavadlav, J.; Podgornik, R.; Praprotnik, M. Adaptive Resolution Simulation of a DNA Molecule in Salt Solution. J. Chem. Theory Comput. 2015, 11 (10), 5035−5044. (38) Fogarty, A. C.; Potestio, R.; Kremer, K. Adaptive Resolution Simulation of a Biomolecule and Its Hydration Shell: Structural and Dynamical Properties. J. Chem. Phys. 2015, 142 (19), 195101. (39) Van Gunsteren, W. F.; Berendsen, H. J. C. A Leap-Frog Algorithm for Stochastic Dynamics. Mol. Simul. 1988, 1, 173−185. (40) Hansen, J.-P.; McDonald, I. Theory of Simple Liquids 2006, 78. (41) Izvekov, S.; Voth, G. A. Multiscale Coarse Graining of LiquidState Systems. J. Chem. Phys. 2005, 123 (13).13410510.1063/ 1.2038787 (42) Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24 (13), 1624−1636. (43) Wang, H.; Junghans, C.; Kremer, K. Comparative Atomistic and Coarse-Grained Study of Water: What Do We Lose by CoarseGraining? Eur. Phys. J. E: Soft Matter Biol. Phys. 2009, 28 (2), 221−229. (44) Hess, B.; Uppsala, S.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (45) Mark, P.; Nilsson, L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J. Phys. Chem. A 2001, 105 (43), 9954−9960. (46) Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. A Generalized Reaction Field Method for Molecular Dynamics Simulations. J. Chem. Phys. 1995, 102 (13), 5451. (47) Neumann, M. Dipole Moment Fluctuation Formulas in Computer Simulations of Polar Systems. Mol. Phys. 1983, 50, 841− 858. (48) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18 (12), 1463−1472. (49) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (50) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6. J. Comput. Chem. 2004, 25 (13), 1656−1676. (51) Fuchs, P. F. J.; Hansen, H. S.; Hünenberger, P. H.; Horta, B. A. C. A GROMOS Parameter Set for Vicinal Diether Functions: Properties of Polyethyleneoxide and Polyethyleneglycol. J. Chem. Theory Comput. 2012, 8 (10), 3943−3963. (52) Daura, X.; Mark, A. E.; Van Gunsteren, W. F. Parametrization of Aliphatic CHn United Atoms of GROMOS96 Force Field. J. Comput. Chem. 1998, 19 (5), 535−547. (53) Oostenbrink, C.; Soares, T. A.; Van Der Vegt, N. F. A.; Van Gunsteren, W. F. Validation of the 53A6 GROMOS Force Field. Eur. Biophys. J. 2005, 34 (4), 273−284. (54) Fredrickson, G. H. The Theory of Polymer Dynamics. Curr. Opin. Solid State Mater. Sci. 1996, 1 (6), 812−816. (55) Chau, P.-L.; Hardwick, a. J. A New Order Parameter for Tetrahedral Configurations. Mol. Phys. 1998, 93 (3), 511−518. (56) Kim, J.; Kim, P.; Lee, S. H.; Suh, G. Y.; Yoon, J. Effect of the PEG or PMMA Micro-Patterned Surface Roughness on Bacterial Adhesion. In 2nd ASM-IEEE EMBS Conference on Bio-, Micro- and Nanosystems; IEEE: San Francisco, CA, 2006; pp 97−99. (57) Wagner, V. E.; Koberstein, J. T.; Bryers, J. D. Protein and Bacterial Fouling Characteristics of Peptide and Antibody Decorated Surfaces of PEG-Poly(acrylic Acid) Co-Polymers. Biomaterials 2004, 25 (12), 2247−2263. M


Article

The Journal of Physical Chemistry B (58) Mondal, J.; Choi, E.; Yethiraj, A. Atomistic Simulations of Poly(ethylene Oxide) in Water and an Ionic Liquid at Room Temperature. Macromolecules 2014, 47 (1), 438−446. (59) Lee, H.; Venable, R. M.; Mackerell, A. D.; Pastor, R. W. Molecular Dynamics Studies of Polyethylene Oxide and Polyethylene Glycol: Hydrodynamic Radius and Shape Anisotropy. Biophys. J. 2008, 95 (4), 1590−1599. (60) Couper, A.; Stepto, R. F. T. Diffusion of Low-Molecular Weight Poly(ethylene Oxide) in Water. Trans. Faraday Soc. 1969, 65, 2486. (61) Mark, J.; Flory, P.; American, O. F. T. H. E.; Society, C. The Configuration of the Polyoxyethylene Chain. J. Am. Chem. Soc. 1965, 87 (7), 1415−1423. (62) Wahab, S. A.; Matsuura, H. IR Spectroscopic Study of Conformational Properties of a Short-Chain Poly (Oxyethylene) (C1E3C1) in Binary Mixtures with Liquids of Different HydrogenBonding Abilities. Phys. Chem. Chem. Phys. 2001, 3, 4689−4695. (63) Tasaki, K. Poly (Oxyethylene) -Water Interactions: A Molecular Dynamics Study. J. Am. Chem. Soc. 1996, 118, 8459−8469. (64) Oesterhelt, F.; Rief, M.; Gaub, H. E. Single Molecule Force Spectroscopy by AFM Indicates Helical Structure of Poly (EthyleneGlycol) in Water. New J. Phys. 1999, 1, 1−11. (65) Jamadagni, S. N.; Godawat, R.; Garde, S. Hydrophobicity of Proteins and Interfaces: Insights from Density Fluctuations. Annu. Rev. Chem. Biomol. Eng. 2011, 2 (1), 147−171. (66) Kreis, K.; Fogarty, A. C.; Kremer, K.; Potestio, R. Advantages and Challenges in Coupling an Ideal Gas to Atomistic Models in Adaptive Resolution Simulations. Eur. Phys. J.: Spec. Top. 2015, 224 (12), 2289−2304.

N


Hybrid Atomistic and Coarse-Grained Molecular Dynamics

Recommend Documents