A polarizable coarse-grained water model for dissipative particle dynamics Emanuel K. Peter and Igor V. Pivkin Citation: The Journal of Chemical Physics 141, 164506 (2014); doi: 10.1063/1.4899317 View online: http://dx.doi.org/10.1063/1.4899317 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/16?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Equilibration of complexes of DNA and H-NS proteins on charged surfaces: A coarse-grained model point of view J. Chem. Phys. 141, 115102 (2014); 10.1063/1.4895819 Coarse-grained modeling of DNA oligomer hybridization: Length, sequence, and salt effects J. Chem. Phys. 141, 035102 (2014); 10.1063/1.4886336 A nucleotide-level coarse-grained model of RNA J. Chem. Phys. 140, 235102 (2014); 10.1063/1.4881424 Coarse-grained simulations of DNA overstretching J. Chem. Phys. 138, 085101 (2013); 10.1063/1.4792252 A coarse-grain three-site-per-nucleotide model for DNA with explicit ions J. Chem. Phys. 135, 165104 (2011); 10.1063/1.3652956

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

THE JOURNAL OF CHEMICAL PHYSICS 141, 164506 (2014)

A polarizable coarse-grained water model for dissipative particle dynamics Emanuel K. Peter1 and Igor V. Pivkin1,2,a) 1 2

Institute of Computational Science, Faculty of Informatics, University of Lugano, Lugano, Switzerland Swiss Institute of Bioinformatics, Lausanne, Switzerland

(Received 7 August 2014; accepted 14 October 2014; published online 28 October 2014) We present a polarizable water model for the Dissipative Particle Dynamics (DPD) method. Employing long-range electrostatics and Drude oscillators, we calibrate the model using the compressibility and the dielectric constant of water. We validate the model by sampling the dielectric properties of solutions of sodium chloride at various concentrations. Additionally, we apply our model in equilibrium and electroporation simulations of a pure dipalmitoylphosphatidylcholine (DPPC) bilayer, a pure cholesterol domain and a mixed DPPC-cholesterol membrane in polarizable water. Finally, we simulate the transport of a short DNA segment through a DPPC bilayer driven by an external electric field. The new water model is suitable for the DPD simulations of systems where polarization effects play an essential role. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4899317] I. INTRODUCTION

Coarse-graining in biomolecular simulation is a well established and very popular method to reduce the computational cost of simulations through a decrease in the model complexity of proteins, sugars, lipids, and water.1–3 Many biomolecular processes involve electrostatic effects and the influence of polarity. For instance in membrane systems differences in the electrostatic potential play an essential role in their biological activity. At present, the most popular all-atom force fields lack electronic polarizability, which is a significant drawback in simulations of systems with ions and polarization effects.4, 5 Polarizable all-atom force fields have been developed, using inducible point dipoles, Drude oscillators, and fluctuating charge models.6–9 Polarizability has been introduced into coarse-grained molecular dynamics simulations,10, 11 while soft sphere dipole models have been developed for the description of polarization effects in enzymatic processes on a coarser level by Warshel et al.,12, 13 where water is represented by point dipoles which reorient in response to an electric field. Recently, inducible dipoles have also been included into coarse-grained force fields,10, 14, 15 while a water model based on Drude oscillators has been implemented into the popular MARTINI force field by Yesylevskyy et al.,2, 3 which was explicitly validated by experimental data. Dissipative Particle Dynamics (DPD) is an efficient method for modeling the mesoscopic behavior of fluids, describing the dynamics of complex systems on long time and length scales,16, 17 while preserving the correct hydrodynamics.18 In prior works, DPD has been extended by adding an electrostatic description in coarse-grained systems using smeared-out charges across a three-dimensional mesh and resolving long range electrostatics through solving the Poisson equation19, 20 and Ewald sums.21 It is worth mentioning, that DPD also has been extended to the description of systems under isobaric, isoenergetic and isoenthalpic conditions.22

In this paper, we present a polarizable water model for DPD, which is based on Drude oscillators. By contrast to other coarse-grained approaches using this Drude model, our approach aims on a DPD description of polarizable water, which enables to simulate larger timescales through the use of the soft potentials in the DPD formulation. To define the model parameters, we chose the compressibility and the dielectric constant of water. We validate the model by comparing the simulation results with experimental data on pure water and ionic solutions. To illustrate possible applications of the model, we use it for equilibrium simulations of model systems which represent dipalmitoylphosphatidylcholine (DPPC), cholesterol, and mixed DPPC-cholesterol membrane systems. We study electroporation effects under the influence of an electric field, as well as DNA transport through a DPPC bilayer, on time and length scales larger than currently accessible by conventional molecular dynamics simulations. II. METHODS

We use the Dissipative Particle Dynamics method16 in the simulations presented here. The time evolution of the Nparticle system is described by Newton’s equations of motion as (1)

dvi = fi , dt

(2)

and

where ri and vi are position and velocity of particle i, respectively. The force fi acting on particle i consists of three additive parts,   R FCij + FD (3) fi = ij + Fij , j =i

a) Electronic mail: [email protected]

0021-9606/2014/141(16)/164506/10/$30.00

dri = vi , dt

141, 164506-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-2

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

which are non-zero within a cutoff radius Rc = 1, defining the length scale of the DPD system. The conservative force FCij is a soft repulsion force, acting along the vector between particles i and j, with parameter aij defining the maximum repulsion between the 2 particles:  aij (1 − rij /Rc )ˆrij , rij < Rc C , (4) Fij = 0, rij ≥ Rc with rij = ri − rj , rij = |rij |, and rˆ ij = rij /|rij |. The two other forces are the dissipative force, FD ij , and the random R force, Fij . These forces are given by D FD rij vij )ˆrij , ij = −γ w (rij )(ˆ

written as U (rN )Coul

⎛ erfc(αrij ) 1 ⎝  = ρi ρj 4π 0 r rij i j >i ⎞ ∞ N 2π  α  2⎠ + Q(k)S(k)S(−k) − √ ρ , V π i i k=0 (11)

with e−k /4α Q(k) = , k2 2

(5)

and

S(k) = FRij

= σ w (rij )θij rˆ ij , R

FBij = dU B /drij ,

with U B =

R

2

(7)

1 k (θ − θ0 )2 , 2 θ ij k

(8)

and A FA ij = dU /dθij k ,

with U A =

where θ ijk stands for the instantaneous angle between 3 particles, while θ 0 is the equilibrium angle, and d0 is the equilibrium bond length. Finally, for the description of electrostatics a non-bonded Coulombic force is added, = dU Coul /drij . FCoul ij

(9)

For the electrostatic interactions, we use the Particle ParticleParticle Mesh Ewald (P3ME) method in combination with a Slater type smearing out of charges.19–21 In the following, we will briefly describe this method. The electrostatic interactions between N point charges are described by U (rN )Coul =

1   ρi ρj , 4π 0 r i j >i |rij |

k=

R

1 k (|r | − d0 )2 , 2 D ij

(10)

where ρ i and ρ j are the charge density values of a pair of DPD particles, |rij | is the distance between the charges,  0 and  r are the dielectric constants of vacuum and water at room temperature. In the P3ME treatment the long-range electrostatic energy given by Eq. (10) is decomposed into a real space and a reciprocal space contribution.24, 25 This expression is

qi eikri ,

(12)

(13)

i=1

(6)

where vij = vi − vj , w (rij ) = [w (rij )] and w (rij ) = 1 − rij /Rc are weight functions that depend on the interparticle distance rij and vanish at rij > Rc , σ 2 = 2γ kB T, and θ ij (t) is a randomly fluctuating variable with Gaussian statistics.16 Dissipative and random forces form the DPD thermostat and are linked through the fluctuation-dissipation theorem.23 We note that all forces in DPD act along the line which connects the centers of interacting particles, conserving the linear and angular momentum of the system. The Velocity Verlet algorithm is used to advance the system in time.16 Additionally, we consider two harmonic bonded forces: the force between bonded particles, FBij , and the harmonic angular force, FA ij , given by D

N 

2

2π (m , m , m ), L x y z k = |k|,

(14) (15)

where α is the parameter that controls the contribution in real space, k is the magnitude of the reciprocal vector k, while mx , my , mz are integer numbers. In DPD, the soft conservative force would allow a full overlap of particles at |rij | = 0. This effect leads to infinite ion pair formation at zero distance. We use the approach of Gonzalez-Melchor et al., and apply the charge distribution ρ(|r|) over a three-dimensional mesh,19, 21, 26 q −2|r|/λ ρ(|r|) = e , (16) π λ3 with λ = 0.7Rc and |r| standing for the distance between the particle center and each grid point (see Figure 1(a)). Thus, the reduced interaction potential uij between two charge distributions in DPD, separated by a distance rij from center to center is given by21

 ρi ρj 4π uij (rij ) Rc −2Rc /λrij = r e 1− 1+ , (17)

rij λ ij with

=

e2 . k B T  0  r Rc

(18)

In this paper we introduce a water model using classical Drude oscillators similar to previous works on polarizable water,3, 7, 8 as shown in Figure 1(b). The central DPD particle is connected by stiff bonds to 2 particles representing charges with opposite sign, while both charges are connected over a harmonic bond potential with an equilibrium angle of 0◦ . This model allows to describe hydrophobic and hydrophilic solvation effects more realistically than the single-particle coarsegraining DPD approach for the solvent, as we show in the Sec. III of this paper. In our simulations, we use the numerical density of central water particles Np = 3, a dissipative force coefficient γ = 4, σ 2 = 2γ kB T as random force coefficient, and a DPD temperature kB T = 1. The Drude particles carry a charge of q = ±0.2e. The equilibrium length of the bonds between

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-3

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

FIG. 1. (a) Long-range electrostatics model using smeared-out charges in the P3M-Ewald electrostatics method combined with a Slater type charge distribution over a three-dimensional grid.19, 21 (b) Polarizable water model using the Drude oscillator model. (c) Dielectric constant of pure polarizable water as function of the harmonic angular constant kθ in comparison with experimental data. (d) Compressibility of pure polarizable water with different values for the conservative force parameter aij . We obtained a realistic value for the compressibility of water (κ −1 = 15.621) with a conservative force parameter of aij = 35. (e) Dipole moment distribution for one individual polarizable coarse grained water averaged over time in equilibrium simulations (black curve) and the average distribution over all water molecules (red curve). (f) Site-site radial distribution function of pure water (distance is in DPD units). Models used for simulations of systems with polarizable water, including (g) DPPC, (h) cholesterol, and (i) DNA segment.

central and Drude particles is d = 0.2Rc . We initially set the conservative force parameter aij = 25 for the interaction of the central water particle with the rest of the system. The 2 charged Drude particles interact with other particles only through Coulombic forces, while the intramolecular non-

bonded interactions are excluded. An harmonic angle potential with an equilibrium angle θ 0 = 0◦ and a force constant kθ are added to control the rotation and the dipole moment of each polarizable water molecule. The average dipole moment is expected to be exactly zero in an apolar environment, which

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-4

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

leads to the equilibrium angle of 0◦ between the charges.3 We adjusted the key parameter kθ for obtaining the correct dipolar distribution, and used the dielectric constant of pure water as measure for the correctness of our model. The dipole moment μi is calculated using μi = ρi |rij |, where |rij | is the distance between the centers of mass of two charge densities, while both centers have an equal absolute charge density of ρ. The dielectric constant  is computed from the correlation  of the total dipole moment M = N i=1 μi of N particles, expressed by the Kirkwood factor gk 25, 27 as gk =

1 (M 2  − M2 ), N μ2

(19)

where μ is the molecular dipole moment and N is the number of particles. Using25, 27 4 ( − 1)(2r + 1) = π Np gk μ2 /(kB T ), (2r + ) 3

(20)

and conducting boundary conditions with  r = ∞, we calculate the dielectric constant  of water as 4  = 1 + π Np gk μ2 /(kB T ). 3

(21)

Each DPD particle in our simulations of pure water corresponds to Nm = 3 water molecules,16, 28 defining the unit of length as  Rc =

Np Nm Vm NA

13

= 6.46 Å.

kθ

θ0

kD

|q|

d0

1 kB T/degrees2

0◦

1 × 105 kB T

0.2e

0.2Rc

less compressibility (κ −1 ) varies from 6.7 to 15.62 when the conservative force parameter aij for central water particles changes from 15.0 to 35.0 (see Figure 1(d)). The dimensionless compressibility is lower in comparison to the standard DPD fluid due to the additional attraction between the polarizable waters in our system, which effectively reduces purely repulsive conservative DPD forces. We chose a conservative force coefficient of aij = 35 for interactions between central water particles to match the dimensionless compressibility of water and use γ = 4.0 for the dissipative force. We use these parameters in all simulations presented in the rest of the paper. As mentioned earlier, one coarse-grained water in our simulations corresponds to 3 effective water molecules with Rc = 6.46 Å. We used this scaling between the DPD and physical length scales in all of our analyses and for comparison with the experiment. The time scale is fixed by matching the self diffusion constant of water.16 For a conservative force parameter aij = 35 we obtain a DPD unit of time equal to16, 20 τ=

(22)

Hence, the dipole moment μ is rescaled according to μ = μDP D Rc .

TABLE I. Parameters of the Drude oscillator DPD water model, i.e., angular constant kθ , equilibrium angle θ 0 , bond force constant kD , charge of Drude particles |q|, and equilibrium bond length d0 .

(23)

We employ the same degree of coarse-graining of water as in the work of Groot and Warren.16 We note, that several coarse-grained water models with different levels of coarsegraining have been developed, with the degree of coarsegraining dependent on the length scale described in the underlying models.29, 30 In Figure 1(c), the dielectric constant  of pure polarizable water is shown as a function of the angular harmonic constant kθ . We consider three different cases for assigning masses to the central and Drude particles. Specifically, in the first case, all particles have an equal mass of 1. In the second case, all particles have an equal mass of 1/3. And in the third case, the central particle has a mass 0.98, while Drude particles have equal mass of 0.01. We find that these three different mass distributions lead to an equivalent dielectric constant of approximately 80 for kθ = 1 kT/(degrees2 ), which equals the experimental value. For our final model of polarizable water, we chose a mass distribution of 1 for all three particles, because this model allows comparably large DPD timestep dt = 0.001 for integrating the equations of motion. This timestep corresponds to 25 ps in real units according to Eq. (24) below. The parameters of the Drude oscillator DPD water model are summarized in Table I. The polarizable water model behaves as a typical DPD fluid with a quadratic equation of state, while the dimension-

Nm Dsim Rc2 = 25 [ns]. Dwater

(24)

The bulk properties of the model are summarized in Table II. The numerical density of Np = 3 corresponds to the normal density of water,16 while the heat capacity is higher by one order of magnitude due to the larger fluctuations in the DPD system (experimental value is 18 cal/(mol K)31, 32 ). The dipole distribution of the model system has a maximum at 5.2 D (shown in Figure 1(e)), which is the maximum value of its dipole moment (see Table II). This value is in good agreement with the maximum dipole moment of an imaginary cluster of 3 real water molecules (5.52 D). In the radial distribution function (RDF) of bulk DPD water, we see a maximum at about r = 0.8Rc (see Figure 1(f)), which shows that our model resembles a typical DPD fluid.16 III. RESULTS A. Bulk electrolyte

We first consider bulk properties of sodium chloride solutions, using the same concentration of electrolyte as in previous studies.20, 21, 34 The system with dimensions 10 × 10 × 10 in DPD units contains 2804 (+2 × 2804 Drude particles) DPD waters and 186 ions (98 positive and 98 negative TABLE II. Bulk properties of Drude oscillator DPD water model. ρ [g/l]

cv [cal/(mol K)]

μmax [ D ]

997.87

184.3

5.2

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-5

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

B. Equilibrium simulations and electroporation of a dipalmitoylphosphatidylcholine (DPPC) membrane

FIG. 2. (a) Site-site radial distribution functions of 0.6 M solution in polarizable water. (b) Dielectric constant of sodium chloride solution as a function of concentration in comparison with experimental values.33

ions with net charge of ±1e). In physical units, the system corresponds to an aqueous solution with a salt concentration of 0.6 M.20 The DPD force parameters used in simulations are aij = 35, γ = 4.0, σ 2 = 2kB Tγ , and kB T = 1. The real forces in the P3ME scheme are truncated at Rcreal = 2.0Rc with α = 1.4. We calculate the radial distribution functions between equal and unequal ion pairs, as well as between the central particles of waters. The results (see Figure 2(a)) are in good agreement with prior works.20, 21, 34 We observe that direct ion pairs occur at a distance r = 0.7. A second peak at r = 1.7 indicates that to a minor extent separated ion pairs occur, where one single water molecule is located between two ions. To further validate our approach, we simulate sodium chloride solutions at concentrations ranging from 0.15 M up to 0.75 M and calculate their dielectric properties. A comparison of the results with the experimental data33 is shown in Figure 2(b). We note that the model slightly overestimates the dielectric constant at concentrations below 0.3M, while the dielectric constant is slightly underestimated at concentrations above 0.4 M.

We continue with simulations of a model of a DPPC membrane immersed in polarizable water. The DPD method has been successfully used to study variety of physical and biological systems.35–42 However, when modeling coarse-grained molecular systems, in most cases rigorous quantitative mapping from atomistic to effective coarsegrained interactions is still lacking.43 Using the typical DPD methodology approach here, we represent the DPPC molecule by a collection of beads. Specifically, we coarsegrain the hydrophilic head group as a charged dipolar group attached to 2 hydrophobic tails represented by 6 beads (see Figure 1(e)). All beads are connected by a harmonic bonds with equilibrium length of 0.2 Rc and a harmonic angle with equilibrium angle of 180◦ . The bond force constant is equal to 1 × 105 kB T, while the angular force constant is 1.0 kB T/degrees2 . For the parameterization of the membrane system, we scale the DPD interaction parameters in order to obtain a realistic phase behavior of the membrane under field free conditions ( Eext = 0). The anisotropy in the motion of embedded cholesterol into membranes has been observed, indicating that the tails of lipids are in a liquid ordered phase.44 It is also well known that DPPC forms stable vesicles with spontaneous curvatures in monolayers.45, 46 In particular, we systematically vary the conservative force parameter aij of the solvent-membrane interaction and the interaction between the hydrophobic tails from a very low repulsion (e.g., aij = 2.0 between the DPPC head groups and solvent) to a higher level of repulsion (e.g., aij = 15.0 between the DPPC head groups and solvent) for adjusting the relative hydrophobicity between the tails and water, as well as the interactions between the tails. This procedure is adapted from prior models which, however, used a different description of the potentials. These models aim on an inclusion of characteristic phase transition of lipid monolayers and bilayers.47–51 We note that other systematic coarse-graining methods exist,52–54 which in turn are not easily transferable to the formulation of the DPD forces.55 Throughout the parametrization procedure, we observe different phases of the DPPC membrane, i.e., solvation of DPPC in water and partial domain formation (Figures 3(a) and 3(b)), until we reach an empirical convergence to a membrane with elastic properties, partial density fluctuations, and gaussian curvature (Figure 3(c)). A further increase of the interaction parameters leads to a fully crystalline (and frozen) membrane (Figure 3(d)). We note, that we keep the electrostatic interaction as well as the water-water interaction constant throughout this procedure. The interaction parameters for the DPPC membrane system are summarized in Table III. We simulate the membrane consisting of 722 DPPC molecules in a periodic box of size 10 × 10 × 10 in DPD units, filled with 2760 waters. Tails of DPPC are mostly aligned perpendicular to the z-axis of the box, the membrane is parallel to x-y plane and located in the middle of the box. The width of the resulting bilayer is about 26 Å, which can be seen in the density profile shown in Figure 3(e). We find that water forms layers along the membrane, which is in qualitative agreement with previous

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-6

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

FIG. 3. (a) Solvation of DPPC particles with conservative force parameter of aij = 25 between DPPC and solvent. (b) Partial domain formation but still no stable membrane with a larger degree of attraction between the lipid tails. (c) Ideal membrane model, parameters are summarized in Table III. (d) Crystallization of DPPC with too high level of repulsion between hydrophobic tails. (e) Density profile of DPPC and water in equilibrium simulation. (f) Electric potential of water and DPPC with z = 0 as reference point. (g)–(l) Electroporation simulation of water-DPPC system. (h) and (i) Start of the poration after few ns. (j) Migration and partial poration of membrane. (k) and (l) Rupture and poration of the membrane.

TABLE III. Conservative force parameter aij used in the simulations. See Figure 1 for the models chosen in our coarse-grained DPD approach.

Hydrophobic Strongly hydrophobic Hydrophilic I Hydrophilic II Hydrophilic III Water Ion

Hydrophobic

Strongly hydrophobic

Hydrophilic I

Hydrophilic II

Hydrophilic III

Water

Ion

3.0 7.0 10.0 7.0 7.0 45.0 45.0

7.0 3.0 10.0 15.0 25.0 75.0 75.0

10.0 10.0 10.0 7.0 7.0 25.0 25.0

7.0 15.0 7.0 7.0 7.0 25.0 25.0

7.0 25.0 7.0 7.0 7.0 25.0 25.0

45.0 75.0 25.0 25.0 25.0 35.0 75.0

45.0 75.0 25.0 25.0 25.0 75.0 35.0

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-7

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

findings by Kasson et al.56 There is an overlap of density profiles of water and coarse-grained membrane in a range of 2.5–4 Å, which is caused by the fluctuations of the membrane. Regions with different densities along the x-y plane are present in the membrane, which leads to slight curvature formation. Water and membrane electrostatic potentials are aligned along the z-direction perpendicular to the membrane (see Figure 3(f)) and are correlated with density profiles. These effects are reminiscent of long range ordering of water which has been observed by Valley et al.57 We also find good qualitative agreement of our results with findings of Kasson et al.56 and IR-spectroscopy experiments.58 They observed anisotropic “layering” behavior of water near membranes due to the interactions between water, hydrophilic head groups and ions. Next, we apply an electric field with a magnitude of 0.6 V/cm in the direction perpendicular to the DPPC membrane, which acts on each charge density ρ i of particle i as FE i = ρi Eext ,

(25)

where Eext = (Ex , Ey , Ez ) stands for the applied electric field. We observe that water permeates a small pore after a simulation time of few ns (see Figures 3(g)–3(i)). At a later stage, the membrane deforms partially, while the water still traffics through the pore (see Figure 3(j)) at a timescale of hundreds of ns. Finally, the pore enlarges and the membrane tilts (see Figures 3(k) and 3(l)) after a total of approximately 5 μs. Pore formation on a nanosecond timescale agrees well with findings of Vasilkoski et al.,59 who investigated the mechanism of pore formation induced by nanosecond electric field pulses. On a molecular basis, the mechanism of pore formation also agrees well with prior molecular dynamics investigations.60, 61 In addition Knorr et al. observed large scale wrinkling after induction of an electric field pulse

at wavelengths in the micrometer range, which is in agreement with simulation results on longer timescales.62

C. Pure cholesterol slab immersed in polarizable water, equilibrium, and electroporation simulations

Pure phases of cholesterol have been observed in experiment,63 which is the reason for considering cholesterol slab in our simulations. Cholesterol rich regions in lipid bilayers may play a role as accumulation spots for amyloid fibrils,64, 65 while cholesterol also might play a protective role against protein-induced membrane disruptions.66 Mainali et al. observed that cholesterol slabs have crystalline behavior.63 In simulations, we approximated the individual cholesterol molecule as a linear chain with 3 bead types – hydrophilic beads carrying a negative charge, hydrophobic and strongly hydrophobic beads (see Figure 1(f)). All beads are connected by a harmonic bonds with equilibrium length of 0.2Rc and force coefficient 1.0 × 105 kB T, as well as a harmonic angle with equilibrium angle of 180◦ and angular force coefficient of 1.0 kB T/degrees2 . As in case of DPPC, this DPD representation of an individual cholesterol molecule is not rigorously linked to its atomistic representation. Here, we chose parameters for the model to obtain a stable membrane with crystalline properties. Resulting interaction parameters for the beads are summarized in Table III. We model the cholesterol bilayer (shown in Figure 4(a)) in a box with dimensions 10 × 10 × 10 in DPD units together with 2592 polarizable water molecules and 57 sodium ions. The density profiles of polarizable water, membrane and ions are shown in Figure 4(b). Profiles partially overlap due to the fluctuations of the membrane resulting in regions with slight curvature. We find that ions agglomerate in a second layer above the high density layer of water. Interestingly, the positively charged ions do not interact with the hydrophilic

FIG. 4. (a) Coarse-grained DPD model system of a cholesterol bilayer in its final configuration. (b) Density distribution of components of membrane water system as function of position along z-axis perpendicular to the membrane. (c) Electric potential in the simulation box calculated from the reference point z = 0. (d) Strong bending of cholesterol bilayer under the influence of an external electric field in the electroporation simulations. (e) Partial domain formation of cholesterol after the rupture of the bilayer has occurred. (f) Vesicle formation on the timescale of μs.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-8

E. K. Peter and I. V. Pivkin

head group of cholesterol. This is due to the electrostatic potential in the direction perpendicular to the membrane, which is mainly caused by the central slab of cholesterol molecules, as shown in Figure 4(c). To simulate an electroporation process of this system, we add an external electric field force 0.6 V/cm in z-direction perpendicular to the slab. The time-evolution of the cholesterol water system under influence of the electric field is shown in Figure 4. The cholesterol bilayer is at first bending under the influence of the electric field, and later is porated by positively charged ions (see Figure 4(d)) at a timescale of few ns. Following the poration by ions, the membrane forms a deep well and finally ruptures after the curvature of the bilayer has become too high (see Figure 4(e)) after approximately 200 ns. After the rupture, cholesterol forms vesicles which are different in shape and size at a timescale of 1 μs. These vesicles are separated by water and aligning with the electric field (shown in Figure 4(f)). We note that with the strength of the external electric field, which we applied to this system, the resistance against poration is too large and the membrane is rather ruptured through an electrophoretic migration of water on a long timescale.

J. Chem. Phys. 141, 164506 (2014)

the influence of an external electric field, the mixed bilayer shows very strong resistance against rupture or poration by water. Specifically, cholesterol starts to accumulate in the region with the highest tension at the peak of the bent bilayer and stabilizes the membrane against rupture. When the bilayer reaches a point of extremely high curvature, the membrane ruptures and vesicles form, while each of the vesicles incorporates single cholesterol molecules (see Figures 5(d)–5(f)). These newly formed vesicles migrate with the electric field and remain stable. Our findings are in agreement with experiments of Kakorin et al.,67 who found that cholesterol reduces the extent of membrane electroporation and electroelongation of vesicles. Experimental results by Koronkiewicz et al. also suggested that domain formation of cholesterol might stabilize membranes against electroporation.68 Resistance against vesicle deformation and poration of cholesterol rich vesicles has been observed in experiments by Sadik et al.69 In experiments by Zischka et al., it has been demonstrated that the application of a long-time electric field onto Mitochondria results in an irreversible electrophoretic perfusion damage,70 which might be also an indication that our findings of vesicle formation over long timescales are realistic.

D. Equilibrium simulations and electroporation of a mixed cholesterol-dipalmitoylphosphatidylcholine (DPPC) membrane

E. Translocation of an anionic DNA fragment through a model DPPC membrane driven by an electric field

We consider a mixed DPPC/cholesterol model system, where we distribute 22 cholesterol molecules across DPPC membrane consisting of 490 molecules. In equilibrium, cholesterol aligns itself between the lipid tails of DPPC (see Figure 5), and accumulates in regions where DPPC membrane has curvature (see Figures 5(a) and 5(b). Cholesterol has a stabilizing effect on the bilayer in our simulation model. Under

In simulations of the translocation process of a small DNA fragment through the DPPC membrane driven by an external electric field, we use the same parameters for DPPC as described above. We model the DNA fragment as a strongly anionic polymer consisting of 4 beads with q = −e as point charge at each bead (see Figure 1(g)). In Figure 6, transport of the DNA fragment across the DPPC membrane driven by

FIG. 5. Simulations of a mixed DPPC and cholestrol bilayer model. Cholesterol molecules are shown with green, DPPC heads and tails beads are shown with red and cyan color, respectively. (a) and (b) Equilibrium bilayer configuration with cholesterol accumulating in bending regions. (c)–(f) Simulations of cholesterol and DPPC mixture under the influence of an external electric field.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-9

E. K. Peter and I. V. Pivkin

J. Chem. Phys. 141, 164506 (2014)

FIG. 6. (a)–(f) Translocation of a small anionic DNA fragment driven by an electric field over a total simulation time of 425 ns.

an external electric field is shown. In an equilibrium simulation of this system, the DNA fragment remains in the bulk and does not adsorb to the membrane due to its strong affinity for the solvent. When an electric field is applied, the DNA approaches the membrane with its solvation shell and starts to penetrate the hydrophilic tail region of the DPPC membrane (see Figures 6(a) and 6(b)), while the water layer near the membrane is perturbed. Subsequently, the DNA fragment enters the hydrophobic core of the membrane and is completely desolvated at this stage (see Figure 6(c)). When the nucleotide enters the other side of the membrane in the hydrophilic tail region, water starts to follow the fragment through the newly formed pore (see Figures 6(d) and 6(e)). Finally, after the transport is completed, the membrane closes the pore and a new solvation shell is formed around the DNA fragment. However, the membrane structure is still perturbed and forms a well in the region where the DNA has penetrated the bilayer (see Figure 6(f)). It has been widely accepted that electric field induced DNA transport is governed by an electrophoresis mechanism of DNA through the membrane as reported previously by Xie et al.71 The formation of pores in the membrane through which DNA entries after surface binding along the transfection process also has been reported.72 IV. SUMMARY

We presented a polarizable water model based on classical Drude oscillators for Dissipative Particle Dynamics method. The model was parameterized using the dimensionless compressibility and dielectric constant of water. We validated the model using simulations of sodium chloride solutions at different concentrations. Good quantitative agreement of the dielectric constant with experimental data33 was obtained. To illustrate possible applications of the new model,

we simulated the electroporation of 3 different bilayers and also a translocation process of a small DNA segment through a membrane. The obtained results are consistent with experimental observations. We observed two different mechanisms of poration in case of the pure DPPC bilayer and mixed cholesterol-DPPC membrane. While pores form in the first case with minimal bending of the DPPC bilayer, the mixed membrane is highly stabilized due to an accumulation of cholesterol in the regions with curvature. The surface tension of cholesterolDPPC membrane is apparently much higher than in case of the pure DPPC bilayer, because it resists longer (approximately tens of ns) against rupture or poration. The rupture mechanism resembles an electrophoretic movement of water on large timescales. We speculate that for a fast poration mechanism, the strength of the external electric field has to be much higher than 0.6 V/cm, which we used in our simulations. We note that our modeling approach may overestimate the stability of this mixed bilayer system, due to our empirical parametrization approach. Although systematic mapping of atomistic interactions to effective coarse-grained potentials is lacking in our DPD simulations, our modeling approach is in good qualitative agreement with experimental data. Finally, we mention that the new model is computationally inexpensive comparing to standard molecular dynamics simulations on an all-atom level. For the electroporation calculations on the DPPC system we used 32 cores and 24 h of computation time on Swiss National Supercomputer Center (CSCS) Cray XE6 Monte Rosa supercomputer, reaching 5 μs. In comparison, simulations of 100 ns of an electroporation of all atom system (33 000 atoms) needed approximately 48 h using 72 cores on Texas Advanced Computing Center Ranger supercomputer.73 The speedup of the DPD simulations in relation to all-atom molecular dynamics simulations is about 225.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30

164506-10

E. K. Peter and I. V. Pivkin

We believe that the presented model can be used to study the effects such as water ordering and long time nonequilibrium electric field conditions, which to the best of our knowledge has not been done with Dissipative Particle Dynamics before. ACKNOWLEDGMENTS

This work was supported by the Swiss National Science Foundation (Grant No. 200021_138231) and a grant from Swiss Platform for Advanced Scientific Computing. Simulations were carried out at Swiss National Supercomputer Center (CSCS) under the project u4. Exceptional support of Maria Grazia Giuffreda from CSCS is gratefully acknowledged. We acknowledge Zachary Levine for providing information regarding the computation time for all atom electroporation simulations of bilayer systems. 1 M.

Levitt and A. Warshel, Nature (London) 253, 694 (1975). J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. de Vries, J. Phys. Chem. B 111, 7812 (2007). 3 S. O. Yesylevskyy, L. V. Schaefer, D. Sengupta, and S. J. Marrink, PLoS Comput. Biol. 6, e1000810 (2010). 4 T. A. Halgren and W. Damm, Curr. Opin. Struct. Biol. 11, 236 (2001). 5 A. Warshel, M. Kato, and A. V. Pisliakov, J. Chem. Theory Comput. 3, 2034 (2007). 6 D. van Belle, M. Froeyen, G. Lippens, and S. J. Wodak, Mol. Phys. 77, 239 (1992). 7 G. Lamoureux, A. D. MacKerell, and B. Roux, J. Chem. Phys. 119, 5185 (2003). 8 G. Lamoureux and B. Roux, J. Chem. Phys. 119, 3025 (2003). 9 S. W. Rick, S. J. Stuart, and B. J. Berne, J. Chem. Phys. 101, 6141 (1994). 10 T. Ha-Duong, N. Basdevant, and D. Borgis, Chem. Phys. Lett. 468, 79 (2009). 11 S. Riniker and W. F. van Gunsteren, J. Chem. Phys. 134, 084110 (2011). 12 A. Warshel and M. Levitt, J. Mol. Biol. 103, 227 (1976). 13 A. Warshel, J. Phys. Chem. 83, 1640 (1979). 14 T. Yan, C. J. Burnham, M. G. D. Popolo, and G. A. Voth, J. Phys. Chem. B 108, 11877 (2004). 15 W. Jiang, Y. Wang, and G. A. Voth, J. Phys. Chem. B 111, 4812 (2007). 16 R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997). 17 S. Jury, P. Bladon, M. Cates, S. Krishna, M. Hagen, N. Ruddock, and P. Warren, Phys. Chem. Chem. Phys. 1, 2051 (1999). 18 P. Espanol, Phys. Rev. E 52, 1734 (1995). 19 C. Sagui and T. Darden, J. Chem. Phys. 114, 6578 (2001). 20 R. D. Groot, J. Chem. Phys. 118, 11265 (2003). 21 M. Gonzalez-Melchor, E. Mayoral, M. E. Velazquez, and J. Alejandre, J. Chem. Phys. 125, 224107 (2006). 22 M. Lisal, J. K. Brennan, and J. B. Avalos, J. Chem. Phys. 135, 204105 (2011). 23 P. Espanol and P. Warren, Europhys Lett. 30, 191 (1995). 24 P. P. Ewald, Ann. Phys. 369, 253 (1921). 25 D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, New York 1996). 26 H. Saint-Martin, J. Hernandez-Cobos, M. I. Bernal-Uruchurtu, I. OrtegaBlake, and H. J. C. Berendsen, J. Chem. Phys. 113, 10899 (2000). 27 G. B. Litinskii, J. Struct. Chem. 39, 687 (1998). 28 E. E. Keaveny, I. V. Pivkin, M. Maxey, and G. E. Karniadakis, J. Chem. Phys. 123, 104107 (2005). 29 K. R. Hadley and C. McCabe, Mol. Simul. 38, 671 (2012). 30 L. Darre, M. R. Machado, and S. Pantano, Wiley Interdisciplinary Rev.: Comput. Mol. Sci. 2, 921 (2012). 31 C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys. 13, 19663 (2011). 2 S.

J. Chem. Phys. 141, 164506 (2014) 32 D.

Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. (Academic Press, 2002). 33 R. Buchner, G. T. Hefter, and P. M. May, J. Phys. Chem. A 103, 1 (1999). 34 P. B. Warren, A. Vlasov, L. Anton, and A. J. Masters, J. Chem. Phys. 138, 204907 (2013). 35 E. S. Boek, P. V. Coveney, and H. N. W. Lekkerkerker, J. Phys. Condens. Matter 8, 9509 (1996). 36 K. E. Novik and P. V. Coveney, Int. J. Mod. Phys. C 8, 909 (1997). 37 P. Malfreyt and D. J. Tildesley, Langmuir 16, 4732 (2000). 38 J. C. Shillcock and R. Lipowsky, J. Chem. Phys. 117, 5048 (2002). 39 V. Symeonidis, G. Em Karniadakis, and B. Caswell, Phys. Rev. Lett. 95, 076001 (2005). 40 I. V. Pivkin and G. E. Karniadakis, Phys. Rev. Lett. 101, 118105 (2008). 41 X. Li, I. V. Pivkin, and H. Liang, Polymer 54, 4309 (2013). 42 Z. L. Peng, X. J. Li, I. V. Pivkin, M. Dao, G. E. Karniadakis, and S. Suresh, Proc. Natl. Acad. Sci. U.S.A. 110, 13356 (2013). 43 J. K. Brennan and M. Lisal, Mol. Simul. 35, 766 (2009). 44 C. Gliss, O. Randel, H. Casalta, E. Sackmann, R. Zorn, and T. Bayerl, Biophys. J. 77, 331 (1999). 45 L. T. Boni, S. R. Minchey, W. R. Perkins, P. L. Ahl, J. L. Slater, M. W. Tate, S. M. Gruner, and A. S. Janoff, Biochim. Biophys. Acta. 1146, 247 (1993). 46 B. Kollmitzer, P. Heftberger, M. Rappolt, and G. Pabst, Soft Matter 9, 10877 (2013). 47 G. Brannigan and F. L. H. Brown, Biophys. J. 92, 864 (2007). 48 B. West, F. L. H. Brown, and F. Schmid, Biophys. J. 96, 101 (2009). 49 F. Schmid, D. Duechs, O. Lenz, and B. West, Comput. Phys. Commun. 177, 168 (2007). 50 C. Stadler, H. Lange, and F. Schmid, Phys. Rev. E 59, 4248 (1999). 51 D. Duechs and F. Schmid, J. Phys. Condens. Matter 13, 4853 (2001). 52 S. Izvekov and G. A. Voth, J. Phys. Chem. B 109, 2469 (2005). 53 V. Ruehle, C. Junghans, A. Lukyanov, K. Kremer, and D. Andrienko, J. Chem. Theory Comput. 5, 3211 (2009). 54 H. Wang, C. Junghans, and K. Kremer, Eur. Phys. J. E., Soft Matter 28, 221 (2009). 55 A. P. Lyubartsev, M. Karttunen, I. Vattulainen, and A. Laaksonen, Soft Mater. 1, 121 (2002). 56 P. M. Kasson, E. Lindahl, and V. S. Pande, J. Am. Chem. Soc. 133, 3812 (2011). 57 C. C. Valley, J. D. Perlmutter, A. R. Braun, and J. N. Sachs, J. Membr. Biol. 244, 35 (2011). 58 H. Binder, Eur. Biophys. J. 36, 265 (2006). 59 Z. Vasilkoski, A. T. Esser, T. R. Gowrishankar, and J. C. Weaver, Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 74, 021904 (2006). 60 D. P. Tieleman, BMC Biochem. 5, 10 (2004). 61 M. Tarek, Biophys. J. 88, 4045 (2005). 62 R. L. Knorr, M. Staykova, R. S. Gracia, and R. Dimova, Soft Matter 6, 1990 (2010). 63 L. Mainali, M. Raguz, and W. K. Subczynski, J. Chem. Phys. B 117, 8994 (2013). 64 E. Drolle, R. M. Gaikwad, and Z. Leonenko, Biophys. J. 103, L27 (2012). 65 A. Morriss-Andrews, F. L. H. Brown, and J.-E. Shea, J. Phys. Chem. B 118, 8420 (2014). 66 L. Qiu, C. Buie, A. Reay, M. W. Vaughn, and K. H. Cheng, J. Phys. Chem. B 115, 9795 (2011). 67 S. Kakorin, U. Brinkmann, and E. Neumann, Biophys. Chem. 117, 155 (2005). 68 S. Koronkiewicz and S. Kalinowski, Biochim. Biophys. Acta. 1661, 196 (2004). 69 M. M. Sadik, J. Li, J. W. Shan, D. I. Shreiber, and H. Lin, Phys. Rev. E 83, 066316 (2011). 70 H. Zischka, N. Larochette, F. Hoffmann, D. Hamoeller, N. Jaegemann, J. Lichtmannsegger, L. Jennen, J. Mueller-Hoecker, F. Roggel, M. Goettlicher, A. M. Vollmar, and G. Kroemer, Anal. Chem. 80, 5051 (2008). 71 T. D. Xie, L. Sun, H. G. Zhao, J. A. Fuchs, and T. Y. Tsong, Biophys. J. 63, 1026 (1992). 72 T. D. Xie, L. Sun, and T. Y. Tsong, Biophys. J. 58, 13 (1990). 73 Z. A. Levine and P. T. Vernier, J. Membr. Biol. 236, 27 (2010).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 199.17.206.1 On: Fri, 21 Nov 2014 19:25:30