Establishing an electrostatics paradigm for membrane

Subscriber access provided by Macquarie University

Biomolecular Systems

Establishing an electrostatics paradigm for membrane electroporation in Framework of Dissipative Particle Dynamics Rakesh Vaiwala, Sameer Jadhav, and Rochish M. Thaokar J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00573 • Publication Date (Web): 20 Aug 2019 Downloaded from pubs.acs.org on August 21, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Establishing an electrostatics paradigm for membrane electroporation in Framework of Dissipative Particle Dynamics Rakesh Vaiwala, Sameer Jadhav, and Rochish Thaokar∗ Department of Chemical Engineering, Mumbai-400 076, India E-mail: [email protected] Phone: +91 (22) 2576 7241. Fax: +91 (22) 2572 6895

Abstract With an exclusive aim to looking into a mechanism of membrane electroporation on mesoscopic length and time scales, we report the dissipative particle dynamics (DPD) simulation results for systems with and without electrolytes. A polarizable DPD model of water is employed for accurate modelling of long range electrostatics near the waterlipid interfaces. A great deal of discussion on field induced change in dipole moments of water and lipids together with the special variation of electric field is made in order to understand the dielectrophoretic movement of water, initiating a pore formation via an intrusion through the bilayer core. The presence of salt alters the dipolar arrangement of lipids and water, and thereby it reduces the external field required to create a pore in the membrane. The species fluxes through the pore, distributions for bead density, electrostatic potential and stresses across the membrane etc. are used to answer some of the key questions pertaining to mechanism of electroporation. The findings are compared with the molecular dynamics simulation results found in literature, and the

1


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


comparison successfully establishes an electrostatics paradigm for biomembrane studies using DPD simulations.

1

Introduction

In biotechnology and biomedicine, there is a growing interest in techniques that enhance cell membrane permeability by disruption of its integrity either in a transient or an irreversible manner. One such method is electropermeabilization or electroporation, wherein a pore in the membrane is induced by application of an electric field in order to facilitate a transport of genes, antibodies, drugs etc. into a cell. Electroporation has been employed for in vitro DNA plasmids 1 and siRNA cell delivery, 2 and therapeutic cancer treatment 3 - whereby the uptake of anti cancer drugs such as bleomycin and cisplatin to cancer cells is enhanced. Electroporation also finds its application in electrofusion, 4,5 in that the membranes of adjacent cells fuse in the presence of an electric field. The transport of ions and molecules across the skin can also be enhanced by the application of the voltage pulses. 6 Experimental evidences report the aqueous pore formation in lipid bilayers under applied electric fields. 7–10 The molecular mechanism of pore formation was first put forward using atomistic simulations. 11,12 An electric field across a membrane induces the transmembrane voltage (TMV) due to Maxwell-Wagner polarization. 13 When the TMV exceeds a threshold value, which is a function of membrane composition; type of lipid chains; presence of various sterols etc., single or multiple pores can be witnessed. The usual sequence of events in electroporation discernible in simulations comprise of initiation of water fingers into the core of bilayer, leading to water wires spanning across the two leaflets, eventually followed by re-orientation of lipid head-groups which tend to align alongside the water wire to stabilize the hydrophilic pore. Although there is no dearth of MD simulations reporting the mechanism of electric field induced pore formation and membrane rupture; 11,14 the kinetics and energetics of pore growth; 15,16 pore characteristics; 17–19 electroporation in asymmetric lipid membranes; 20 2



charged membranes 21 and archaeal lipid membranes 22 etc., to best of our knowledge, there is no study using DPD method, except our previous studies- Vaiwala et al., 23,24 specifically devoted to investigating the membrane electroporation. Dissipative particle dynamics has been extensively employed virtually in all areas of research, namely flow past rigid objects, 25 polymer melts and solutions, 26–28 colloidal suspensions, 29–31 surfactant systems, 32,33 and biomembranes. 34–36 The method has advantage of simulating relatively larger system sizes on account of coarse-grained models of lipids and water, and allows a larger time stepping by merit of the soft repulsive potential between coarse-grained entities. Resealing of a pore in a reversible pore dynamics is a slow process, which occurs over a few milliseconds time scale. The long-lived pores may even take many seconds for resealing to complete. 37 In such a scenario, DPD is an ideal tool for inspecting larger membrane patch for a longer period of time. However, for the purpose of the present work which is to establish an electrostatics paradigm in DPD formalism for electroporation studies we have restricted the time scale of simulation to tens of order of nanoseconds. In Ref.- Vaiwala et al. 23 , we introduced a protocol that imposes an electric field across a model membrane using charged plates, which are permeable for DPD particles, and demonstrated that the field-induced response in regards to the structural profiles, the electrostatic potential and mechanical stress etc. are quite similar to the results obtained when the membrane is porated using the ion imbalance method. 38 The DPD model of water employed in the work of Vaiwala et al. 23 has no explicit charges to account for local polarization. Instead, the local permeativity was computed by using polarizability assigned to DPD particles. However, such a mean-field approach for polarization may not be as accurate as the water model having explicit charges especially for biomembrane studies, which involve highly inhomogeneous dielectric interfaces of lipid and water. For instance, the onset of poration reported in the work of Vaiwala et al. 23 is few tens of nanoseconds, while it is of order of nanosecond in a typical MD simulation. 12 Moreover, the process of pore creation using the water model without explicit charges is lipid initiative, while MD simulations 12 reported the

3


Page 4 of 36

Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


process to be water initiative. Seeing a need for polarizable version of water, in our previous study 24 we devised a Drude water for DPD simulations. As a part of testing and validation, a brief study on electroporation was also presented. 24 Electrostatics interactions play a crucial role in electroporation, and therefore it is imperative to establish an electrostatics paradigm that shows an accurate response of DPD models of lipid and water during electroporation in DPD formalism. The present work successfully does this by offering a detailed picture of dipolar reorientation of lipid and water using the polarizable DPD water model, and it forms a base case report for studying energetics of membrane using DPD model. Effects of salt are also investigated in the present study, and a rich set of data is provided to support the conclusions emerging from the present study. We report DPD simulations results for electroporation and address the following questions surrounding the pore formation. (1) What factor triggers the membrane to form a defect under applied electrical field? and why? (2) How do the lipids and water respond to an applied electric field? (3) What is the time scale of pre-pore formation? (4) Do the inflows of lipid head-groups and water into the pore more from anodic or cathodic interface? Or are they comparable at two leaflets of the bilayer? (5) How do the electrostatic potential and stress distributions across the membrane vary with applied fields? (6) How does membrane tension change in external electric fields? (7) Does a salt have any effect on electroporation?

2

Dissipative Particle Dynamics

Dissipative particle dynamics is a mesoscopic simulation technique. The particles in DPD are coarse-grained entities, each carrying a few atoms or molecules. The DPD particles evolve with time according to Newton’s equations of motion, and interact with each other within a range of Rc via the pairwise conservative, dissipative and stochastic forces. These forces are short ranged and centrally directed along a line joining the centers of interacting particles. The conservative force (f c ) originates from the excluded volume interactions between the coarse-grained particles with their internal atomistic degrees of freedom averaged out. 39–42 To

4



Page 6 of 36

take into account the rapid atomistic fluctuations within the DPD particles and thereby to correctly describe the dynamics, the dissipative and stochastic forces are included. The conservative interaction is modelled by a soft repulsive potential on account of coarse-graining, and the corresponding force between the particles, say i and j respectively positioned at ri and rj , is fijC

rij eij , = aij 1 − Rc

(1)

within the cutoff radius Rc . The force constant aij regulates the strength of the conservative force. The separation distance rij is |ri − rj |, and eij represents the unit vector joining the position vectors of i and j. Within the interaction cutoff Rc , the dissipative and the random forces follow fijD = −γ wD (eij · vij ) eij ,

(2)

fijR = σ wR ζij eij .

(3)

The vector vi is the instantaneous velocity, and vij = vi − vj . ζij is a random number drawn from a Gaussian distribution. A usual choice for the weight function wR (r) is 1 −

rij . Rc

The

friction parameter γ controls the dissipative strength while the amplitude of the random force is governed by σ. The dissipative and random forces are thermostating forces, because the temperature of DPD system is maintained constant by interplay of these forces according to the fluctuation-dissipation theorem, 43

2 . σ 2 = 2γkB T and wD = wR

(4)

The bonded particles have bond-length (l0 ) and angle-bending (θ0 ) constraints. The bond-stretching is usually modelled by a harmonic spring with spring constant ks , and thus the interaction potential energy is 1 2 uB ij = ks (rij − l0 ) . 2 5


(5)

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The angle-bending or three body interaction can be modelled either by a quadratic bending potential or by a cosine function. 1 2 uA ijk = ka (θijk − θ0 ) or kθ [1 − cos(θijk − θ0 )]. 2

(6)

Here, θ0 is the angle between three particles i, j and k which have the angle constraint, while ka and kθ are the stiffness parameters. The charge species interact via the long-range electrostatic force f E . Equations (1),(2),(3), (5) and (6), taken together with the long-range electrostatic force f E , sum up to a net force fi acting on a charged DPD particle, fi = fiE +

X

fijC + fijD + fijR + fijB +

A fijk .

(7)

j6=k k6=i

j6=i

3

XX

Electrostatics in DPD formalism

As mentioned earlier, the DPD particles interact by the soft repulsive potential, and this allows the non-bonded particles to overlap. Therefore, it is highly likely for two oppositely point-charge species to form an artificial ion-pair. The ion-pairing can be avoided by making the potential well (electrostatics + soft repulsion) shallower than the thermal energy (−kB T ). One way of rendering the potential well shallower than (−kB T ) is to smear out the charges of the ions by some kind of distribution. A complete recipe to solve long range electrostatic interactions with a linear charge smearing was suggested by Groot 44 . González-Melchor et al. 45 adopted an exponential charge distribution (Slater type density), and demonstrated that the resulting pairwise electrostatic force can be written down as the sum of real and reciprocal contributions, which can be computed using Ewald summation. The detail derivation for electrostatics of Slater type distribution and the error estimates, which are used to optimize the electrostatic parameters including the charge smearing length, can be found in our earlier work, Vaiwala et al. 46 . Bore et al. 47 recently addressed an issue of solving electrostatics in non-homogeneous dielectrics, and proposed a density functional-based for6



Page 8 of 36

malism to compute electrostatic forces for mesoscopic systems. In line with Groot’s method of solving electrostatics, the derivatives in generalized Poisson equation are discretized by a finite difference scheme, and the resulting system of linear equations is solved via successive over relaxation. The approach by Bore et al. 47 reproduces the concentration-dependent partitioning of an ideal salt in oil/water systems, and correctly captures the polarization effects attributed to dipoles and charged particles. For solving electrostatics in the present work, we followed Groot’s method with a linear charge smearing function, as delineated below. The distribution of electrostatic potential is governed by the Poisson equation,

∇ · (s (r)∇φ(r)) = −ρc (r)/0 .

(8)

Here ρc (r) is the charge density, 0 denotes the dielectric permittivity of vacuum, while s is regarded as a local dielectric screening constant. We have considered the screening constant s to be uniform, because the DPD water model used in the present study is the polarizable water, which was parametrized to reproduce the dielectric constant of water with the background screening constant s = 2. 24 The Poisson equation is made dimensionless by scaling the electrostatic potential by kB T /e, the length by Rc and the charge by e. Here e is the elementary charge of electron, kB denotes the Boltzmann constant, and T is the temperature of system. The equation (8) in its dimensionless form then scales to

∇2 φ(r) = −Γρc (r).

(9)

Note that φ, ρc and the differential operator ∇ are all dimensionless now. Γ is the nondimensional Bjerrum length, e2 /(0 s kB T Rc ). For smearing of the charges, the computational domain is discretized into cells of size Rc . The point-charges are then smeared out over all lattice nodes within the spheres of Re = 1.6Rc centred around the charges using a

7


Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


linear charge distribution function, 44 3 f (r) = πRe3

r 1− . Re

(10)

By employing the successive over relaxation technique, 48 the field Equation (9) is solved in the real space. Accordingly, the nodal values of potential are updated iteratively via

φ

new

=φ

old

2

− ζ Γ¯ ρe + ∇ φ .

(11)

Here ρ¯e is the total charge on a lattice node, and ζ is the relaxation parameter, which is taken to be 0.15. The iterative loop in equation 11 terminates when the relative error defined in equation (12) reaches a value within a tolerance limit (usually 0.03). sP err =

(Γ¯ ρe + ∇2 φ)2 P 2 2 nodes (∇ φ)

nodes

(12)

With the use of the same distribution function f (r), the grid-based gradients of electrostatic potential are back interpolated to compute the electrostatic force on the charged particles,

fiE (ri ) = −qi

X

f (|r − ri |)∇φ(r).

(13)

nodes

The summation runs over the lattice nodes within the distance Re from the particle location ri . The field gradient ∇φ(r) at lattices is calculated by central difference formula,

∇µ φ(rj ) = [φ(rj+1 ) − φ(rj−1 )]/2h.

The grid size h is chosen to be equal to Rc .

8


(14)


CH3 H3C

CH3

N

Choline

1

CH2

0.47

H2C O

O P O

0.47

Phosphate

2

O

Carbonyl - 1

CH2

H

0.31

C

3 O

0.59

O

CH2

C

Carbonyl - 2

7

O C

O

0.59

4 8

0.59

Hydrophobic tails

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 36

0.59

5

9

0.59

0.59 6 10

(a)

(b)

Figure 1: (a) A schematic of polarizable water model. The Drude particle (DP) carries a charge |q| = 0.83e, the equilibrium bond-length l0 is 0.25 Rc , and the equilibrium angle is θ0 = 0. The mass of central W site is 2/3 m, while Drudes have mass 1/6 m each. (b) A DPD model of DMPC lipid. The bond lengths are specified in unit of Rc , while the numbers within coarse-grained beads are the bead labels.

4

Simulation Details

A bilayer comprised of 614 zwitterionic model lipids, 49 shown in Fig. 1(b), is simulated for electroporation at room temperature with and without electrolytes in a box of 20×20×26 Rc3 size. The bilayer is hydrated with a sufficient quantity of polarizable DPD water 24 (see Fig. 1(a)). The overall bead density is 3 Rc−3 . For lipid and water being four-to-one coarse-grained (Nm = 4), the length scale Rc = 0.711 nm is deduced by mapping the volume of a coarsegrained (CG) water with that of four real water molecules. The time scale of simulation τc is chosen to be the diffusional time scale of the CG water. Mapping the diffusivity of the water, which is obtained using a mean-squared displacement data from a bulk water simulation, the time scale turns out to be τc = 57.034 ± 0.24 ps. The bond-lengths in water and lipids are kept constrained by harmonic springs, with stiffness ks = 512 kB T /Rc2 , while the angle between DP-W-DP is loosely constrained through an angle bending potential u(w) = 21 ka (θ − θ0 )2 to its equilibrium value θ0 = 0. The anglebending constant ka is 1 kB T /rad2 . Referring to Fig. 1(b), the angles in lipids are maintained 9


Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Repulsion force parameter for lipid and water in DPD units aij (kB T /Rc ) Choline Phosphate Carbonyl Tail W site Ions

Choline

Phosphate

Carbonyl

Tail

W site

Ions

100 100 102 130 98 98

100 100 102 130 98 98

102 102 100 110 102 102

130 130 110 100 130 130

98 98 102 130 100 100

98 98 102 130 100 100

Table 2: Three sets of simulated systems studied Item lipid CG water 1:1 electrolyte

System-A 614 25060 0 0-2

System-B 614 24036 1024 ions (∼ 0.25 M) 0-2

System-C 614 20964 4096 ions (∼ 1 M) 0-2

Applied Electric field |Ez | (kB T /eRc ) Simulation runs (repeated) at poration field

R-A1, R-A2 R-A3, R-A4

R-B1, R-B2 R-B3, R-B4

R-C1, R-C2 R-C3, R-C4

by the angle-bending potential u(l) = kθ [1 − cos(θ − θ0 )] to an equilibrium angle θ0 = 120o for the triplet of CG beads 2−3−7, and to an angle θ0 =180o for the triplets of CG beads 2−3−4, 3−4−5, 4−5−6, 7−8−9 and 8−9−10. The stiffness constant kθ is 20 kB T . The soft repulsion force parameters among various interacting sites are tabulated in Table 1. All charge bearing species, namely choline and phosphate beads of lipids, Drude particles of water, and ions of electrolytes, interact via the electrostatics. The long range electrostatic interactions are solved in real-space by successive over relaxations method, 44 as discussed in the previous section. The noise amplitude of random force is chosen to be σ =3 p 3 1 √ (kB T ) 4 m 4 / Rc , while the strength of the dissipative force is kept γ =4.5 (kB T m)/Rc in accordance with the fluctuation-dissipation theorem (Equation 4). Within the same CG water, the particles do not interact at all. As detailed in Table 2, there are three sets of systems simulated under external electric fields. The charges are considered to be in a uniform dielectric background of screening constant s = 2. 24 A constant electrical field E is applied in −z direction across the bilayer

10



3.5

3.5

R-A1

-3

2 1.5

0.05

1.5 0.025

1

0.5

0.5

-2

-1

0

1

2

0.075

2

1

-3

0

3

0

0.3

2

0.24 0.18

1.5

0.12 0.06

-3 -2 -1 0

1

2

1

3

0

-3 -2 -1 0

1 2

3

0.5

-3

-2

-1

0

1

2

0

3

z [nm]

z [nm]

(a)

2.5

0.1 -3

2.5

R-C1 3

ρ(z) [Rc ]

2.5 ρ(z) [Rc ]

3

0

3.5

R-B1

3

-3

ρ(z) [Rc ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 36

(b)

-3

-2

-1

0

1

2

3

z [nm]

(c)

Figure 2: Bead density profiles for (a) the simulation run R-A1, (b) run R-B1 and(c) run R-C1, averaged over simulation time ∼11.41 ns. For comparison, solid lines correspond to density profiles under no applied electric field, while dotted lines refer to the density distributions at applied field Ez = -1.8 kB T /eRc for R-A1 and Ez = -1.4 kB T /eRc for R-B1 and R-C1. Black color indicates choline, red phosphate, green carbonyl-1, orange carbonyl-2, brown lipid, cyan water, and violet total density. The inset depicts ion concentrations for (blue) cation and (pink) anion. by a convenient way of adding a force fi = qi E to all charge bearing species, mimicking the nanosecond intense pulsed electric field (nsPEF) induced electroporation. 50 For each system with its poration field, the simulation is repeated for multiple runs to get better insights and averaging. After equilibration without external electrical field for ∼20 ns, the field is applied in most cases for ∼11.41 ns time. The particles are evolved by integrating Newton’s equations of motion using the modified velocity Verlet scheme 51 with a time step size 0.01 τc .

5

Results and discussion

In order to understand the mechanism of pore formation, below we present our exclusive analysis on the membrane’s response to electroporation in terms of structure of lipid and water; fluxes of choline and water through the core region of bilayer; the field induced dipole re-orientation; a transient change in the head group density; profiles for electrical potential and stresses; and a transient behavior for ionic conductance.

11


Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5.1

Density distribution

Figure 2 shows the distributions of bead densities for lipid, its constituents and water for the three sets of systems studied (see Figs. 1, 2 and 3 in Supporting Information (SI) for all repeated simulation runs). The profiles for head groups (choline, phosphate, carbonyl-1 and -2) become broader with applied electric field, implying the field induced structural redistribution of the particles. The water particles penetrate deeper into the bilayer core. These observations are in agreement with the density profiles reported in the work of Tieleman. 12 The insets in Figs. 2(b) and 2(c) depict the density distribution of free ions in the salt systems B and C, respectively. The concentration of cation near the interfaces is higher than that of anion. The depth of penetration of cations into the core of bilayer is also relatively high compared to anions. The phosphate group repels anions and thus restricts them from further spreading close to the carbonyl groups. It is seen that the screening of the phosphate group at high salt concentration leads to its diminishing influence on the cation distribution around the hydrophilic region of bilayer.

5.2

Species fluxes with time

The detail information about the pore initiation can be known by analyzing a transient change in the species fluxes shown in Fig. 3 (see Figs. 4, 5 and 6 in SI for all repeated runs). As apparent in Fig. 3(P), a sudden rise with time in the number of water particles inside the bilayer core, within a zone −0.6Rc < z < 0.6Rc , signifies the time for prepore (hydrophobic pore) formation. This occurs at times ∼ 7 − 10 ns, ∼ 3 − 4 ns and ∼ 5 − 6 ns for the systems A, B and C, respectively. That a prepore triggers at these times is confirmed by visiting the simulation snapshots given in Figs. 16, 17 and 18 in SI, and also shown in Fig. 4. These images, which are created by visual molecular dynamics software, 52 show that the pore creation process initiates with an intrusion of water, followed by movements of lipid head-groups reaching the hydrophobic pore. This followup of lipids making the hydrophobic pore hydrophilic is implied by a rise in the lipid counts after a short delay (especially for 12



40

R-A1

P

50

20 10

R-B1

R-C1

P 90

40

counts

counts

counts

120

60

P 30

30 20

60 30

10

0

0

0

50

100

200

Q

20 10

60 40 20

0

100 50 0

0

0

1

2

3

4

5

6

7

8

9

10

11

12

0

time [ns]

Q

150 counts

30

(a)

Q

80

counts

40 counts

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 36

1

2

3

4

5

6

7

8

9

10

11

0

12

time [ns]

(b)

1

2

3

4

5

6

7

8

9

10

11

12

time [ns]

(c)

Figure 3: Species flux through a pore in (a) simulation R-A1, (b) run R-B1 and (c) run R-C1, showing (P) number of water (pink) and choline (green) in the core of bilayer, and (Q) cumulative flow of water and choline into the core through anodic (pink and green) and cathodic (blue and brown) interfaces. Also shown in (Q) are the cumulative fluxes for cation (violet and black) and anion (cyan and orange) through anodic (violet and cyan) and cathodic (black and orange) interfaces. For visual clarity, the counts for the choline are shifted by 20 in (a), by 40 in (b) and by 80 in (c). The vertical dotted line indicates onset of pre-pore creation.

Figure 4: Simulation snapshots from run R-A1 showing side views at specified times. Water is indicated by red color with yellow Drudes. Blue color represents the choline bead of lipids. Other head beads and tails are not shown for clarity.

13


Page 15 of 36

250 System-A water

150

Ez = 0 R-A1 R-A2 R-A3 R-A4

lipid

100

250

200 System-B (0.25 M) 150 water

µz [D]

50 0

200

100

50

50

0

-5

-4

-3

-2

-1

0

1

water

-200

2

-250 -6

3

4

5

6

lipid

-5

-4

-3

-2

-1

water

-150 -200

0

1

2

3

4

5

6

-250 -6

z [nm]

z [nm]

(a)

0

-100 water

-150 -150 lipid

water

-50

-100

-100

Ez = 0 R-C1 R-C2 R-C3 R-C4

System-C (1 M)

150

100

-50

-50

-200 -6

lipid

Ez = 0 R-B1 R-B2 R-B3 R-B4

lipid

µz [D]

200

µz [D]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(b)

lipid

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

z [nm]

(c)

Figure 5: Distributions of z-component of dipole moments of lipid and water for (a) system-A, (b) system-B and (c) system-C. The profiles are averaged over simulation times. Solid lines refer to lipid dipoles while dotted lines correspond to water dipoles. Black color indicates µz values at Ez = 0. system A) in Fig. 3(P). As evident in Fig. 3(Q), which shows the cumulative fluxes of various species- water, choline and free ions- into the pore through anodic (+ve electric potential) and cathodic (-ve electric potential) leaflets, water is more likely to run through the anodic leaflet. A reason for more inflow of water from the anodic side is the dielectrophoretic force on water dipoles. In order for lipids to stabilize the pore, more lipids with their choline enter the pore from the anodic interface. The unfavourable dipolar arrangement of lipids under electric field can be a reason for this. Being under electrophoretic drifts, ions also migrate through the pore. However, having been repelled by negatively charged phosphate, anions admit a flux which is quite less than the flux for cations.

5.3

Dipole distribution

When an electric field is applied in a system containing dipolar particles, they try to re-orient themselves to oppose the very cause of their re-orientation- the applied field. Following this notion, the lipid and water dipoles in electroporation tend to align their dipole direction (negative charge to positive charge within a dipole) with the direction of the applied field, which is from the higher electrostatic potential (anodic leaflet) to the lower potential (cathodic 14



0.4 0.2 0

Ez [V/nm]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

-0.02 V/nm

-0.2 -0.4

Ez = 0

Polarization = (εr-1)ε0 Ez

Ez = -1.8 kBT/eRc (R-A1)

µz = Polarization * Volume

-0.6

Ez = -1.8 kBT/eRc (R-A2)

3

µz = -58 D for Vol. = 20x20x0.1 Rc

Ez = -1.8 kBT/eRc (R-A3) Ez = -1.8 kBT/eRc (R-A4)

-0.8 -8 -7 -6 -5 -4 -3 -2 -1 0

1

2

3

4

5

6

7

8

z [nm]

Figure 6: Electrostatic field distribution across the membrane for four simulation runs at Ez = −1.8 kB T /eRc . The symbols refer to () R-A1 run, (4) R-A2, (5) R-A3, and (+) run R-A4. Also shown is the field distribution at ( ) Ez = 0.The profiles are averaged over the simulation times.

leaflet). The dipole values for water and lipids are computed for each slice of thickness 0.1 Rc parallel to the bilayer plane. As expected according to the notion just stated and evident in Fig. 5(a), the peak in z-component dipole moment of lipids at anodic leaflet reduces from its no field value µz = 197 ± 2 D to an average µz = 140 ± 2 D over the simulation run for the system A. Similarly, there is a significant decrease in the peak value of z-dipole moment of water near the cathodic interface from µz = 142 ± 2 D to an average value 98 ± 2 D. The water dipoles facing toward −z direction at anodic leaflet are only slightly enhanced from −µz = 142 ± 2 D to 160 ± 2 D under the applied field. The lipid dipoles at cathodic leaflet, which are also pointing in −z direction, shows only a little change from −µz = 197 ± 2 D to 188 ± 2 D, probably due to interaction with water dipoles. For the systems containing salt, as evidenced in Figs. 5(b) and 5(c), the changes in the peak values of z-dipole moments of interfacial water and lipids show similar qualitative trends as observed in the system without salt. It should be noted that the changes in the peak values of z-dipole moments of lipids have dual causes; a change brought about by local density of lipids at interfaces and a change by field induced re-orientation of the lipid 15


Page 16 of 36

Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Dipole moment of lipids in z-direction within 0.1 Rc thick slabs at interfaces, averaged over four independent simulation runs for systems-A, B and C. All values are within the standard error 2 D. Systems→ Field (kB T /eRc ) ↓ Ez = 0 Ez = −1.4 ∆4

System-A (leaflets) (anode) (cathode) (µz D) (−µz D) 197 197 135.6 177.9 -31.2 -9.7

System-B (leaflets) (anode ) (cathode) (µz D) (−µz D) 221 221 159.2 201.2 -28 -9

System-C (leaflets) (anode) (cathode) (µz D) (-µz D) 235 235 185.6 215.4 -21 -8.3

δ1 1 (leaflets) (anode) (cathode) (%) (%) 12.2 12.2 17.4 13.1 -

δ2 2 (leaflets) (anode) (cathode) (%) (%) 19.3 19.3 36.9 21.1 -

molecules. Although it seems that the applied field causes even lipids at cathodic leaflet to reduce their negative dipole moments, subtracting out the contribution by density reduction in z-dipole moment of these lipids makes it clear that the applied field in fact favours the lipids present in cathodic leaflet to re-orient themselves, enhancing their negative z-dipole moment. With no electrolytes being present in the system A, the water dipoles in the bulk regions are (partially) aligned in the direction of applied field. The water particles in the bulk therefore have the dipole moment as high as µz = −54 ± 2 D in each slice of z, as indicated in Fig. 5(a). This value of µz in bulk water within a slice of volume dV is consistent with µz obtained from the local electric field in bulk water, µz = (r − 1)0 Ez dV = −58 D, as confirmed by calculations from electric field distribution in Fig. 6. The high dielectric constant of water leads to a nearly iso-potential condition in bulk water and thereby vanishingly small net electric field there. However, at the anodic leaflet, where the lipid dipoles are oriented upward, the dipoles create a negative field (i.e. pointing in −z direction), while the lipid dipoles which are oriented downward at the cathodic leaflet create a positive electric field (i.e. pointing in +z direction). Therefore the external electric field, which is applied in negative z direction), is bolstered at the anodic side and diminished at the cathodic side. 5.3.1

Effect of Salt

In the presence of salt, the local polarization induced by electrolytes lowers down the zcomponent of dipole moment of water in bulk regions. For instance at the applied electrical field Ez = −1.4 kB T /eRc , the z-dipole moment within each slice of 0.1Rc thickness (parallel

16



Page 18 of 36

Table 4: Dipole moment of water in z-direction within 0.1 Rc thick slabs at interfaces, The values are averaged over four independent simulations for systems-A, B and C, and are within the standard error 2 D. Systems→ Field (kB T /eRc ) ↓ Ez = 0 Ez = −1.4 ∆4 (%)

System-A (leaflets) (anode) (cathode) (−µz D) (µz D) 141 141 148.2 108.6 5.1 -23

System-B (leaflets) (anode) (cathode) (−µz D) (µz D) 103 103 107.7 77.2 4.6 -25

System-C (leaflets) (anode) (cathode) (−µz D) (µz D) 68 68 78.9 37.8 16 -44.4

δ1 1 (leaflets) (anode) (cathode) (%) (%) -27 -27 -27.3 -28.9 -

δ2 2 (leaflets) (anode) (cathode) (%) (%) -51.8 -51.8 -46.8 -65.2 -

1

δ1 = (System B - System A)/System A. δ2 = (System C - System A)/System A. 3 ∆ = (at E=1.4 - at E=0)/at E=0. 2

to x-y plane) in the bulk water is −39 ± 2 D in no salt case, while it drops to −9 ± 2 D in 0.25 M salt system (see Fig. 5(b)). However, there is no further significant change when the salt is increased to 1 M; the dipole moment remains µz = −7 ± 2 D at both sides of the bilayer in bulk water, as indicated in Fig. 5(c). Thus, the salt-induced local polarization causes the water dipoles in the bulk to be configured more randomly. Referring to Table 3, which compares the z-dipole moment of lipid at varying salt concentration, the salt induces lipids to be more ordered in absence of external field, for the z-dipole moment of lipid increases from |µz | = 197 ± 2 D without salt to 221 ± 2 D (12.2% change) at 0.25 M, and further increases to 235 ± 2 D (19.3%) at 1 M. However, the z-dipole moment of interfacial water tabulated in Table 4 asserts that the water dipoles are randomized by addition of salt, as the dipole moment near the interfaces falls from |µz | = 141 ± 2 D without salt to 103 ± 2 D (−27% change) at 0.25 M, and further decays to 68 ± 2 D (−51.8% change) in 1 M case. This trend of ’alignment of lipid dipoles and dealignment of water dipoles with the bilayer normal by adding a salt’ observed in absence of applied field continues even in the case with applied field, Ez = −1.4 kB T /eRc . The effect is more pronounced on the interfacial water than the lipids though. Since the Debye layer thickness (κ−1 ), which is of order of few angstroms (∼ 6 Å for 0.25 M, and ∼ 3 Å for 1 M salts), is less than the system size Lz = 18.5 nm, and since the simulation run is longer than the membrane charging time (∼ 2 − 10 ps) and the charge relaxation time (∼ 0.01−0.05 ns), the presence of salt is expected to affect poration. With the electric field at the interfaces scaling as Ez ∼ φs κ, φs being the potential at the interfaces, the 17


Page 19 of 36

16

16

16 R-C1

R-B1

R-A1

12

8

8

8

4 0

4 0 -4

-4

0 -4

A

-8

A

-3

-2

-1

C

-8

C -4

C 0

1

2

3

4

5

-12 -5

z [nm]

(a)

4

A

-8 -12 -5

Fz [pN/D]

12

Fz [pN/D]

12

Fz [pN/D]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-4

-3

-2

-1

0

1

2

3

4

5

-12 -5

z [nm]

(b)

-4

-3

-2

-1

0

1

2

3

4

5

z [nm]

(c)

Figure 7: The dielectrophoretic force on a unit dipole along the field direction for simulation runs (a) R-A1, (b) R-B1, and (c) R-C1. The colors refer to the force per unit dipole moment along z-direction under no external fields (violet), and at poration field Ez = -1.8 kB T /eRc (pink) for R-A1, and -1.4 kB T /eRc for R-B1 and R-B2. The black line indicates the density distribution of carbonyl-2 beads just to show precisely the locations of water-carbonyl interfaces. The profiles are averaged over the simulation run times. field Ez near the interfaces enhances by addition of salt. This is a reason why the external field required to build up ∼ 1 V transmembrane potential (TMV) is lower (Ez = −1.4 kB T /eRc ) in the presence of the salt than the threshold electric field (Ez = −1.8 kB T /eRc ) in the salt-free system. Again it is due to enhanced electrical field by screening action of the ionic crowd near the lipid-water interfaces that the lipid dipoles are more ordered than the lipids in salt-free environment. However, the interaction of interfacial water with lipids is now screened by the Debye layer of ions, and therefore the interfacial water dipoles become more disordered with respect to the bilayer normal. 5.3.2

Dielectrophoretic force

We have mentioned earlier that the inflow of water through the pore is more at anodic leaflet than the cathodic one, and that the dielectrophoretic force is responsible for this. Fig. 7 shows the force per unit dipole, Fz /µz = ∂Ez /∂z (refer Figs. 10, 11 and 12 in SI for repeated runs). In absence of external field, the force on water dipoles tends to push the water particles away from the bilayer; the force is in +z direction on negative water dipoles at anodic interface, and it is in −z direction on the positive water dipoles at cathodic side. 18



R-A1

R-B1

R-C1

0.06

0.6

0.06

0.6

0.5

0.04

0.5

0.04

0.5

0.04

0.4

0.02

0.4

0.02

0.3

0

0.3

0

0.02

0.3 0

-0.8 -0.4

0

0.4

-3

-3

0.4

ρ(z) [Rc ]

0.6

ρ(z) [Rc ]

-3

ρ(z) [Rc ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 36

-0.8

-0.4

0

0.4

0.8

0.8

0.2

0.2

0.2

0.1

0.1

0.1

0 -3 -2.5 -2 -1.5 -1 -0.5

0

0.5

1

1.5

2

2.5

0 -3

3

(a)

-2

-1

0

1

2

0 -3

3

z [nm]

z [nm]

(b)

0.06

-2

-0.8

-1

-0.4

0

0

0.4

0.8

1

2

3

z [nm]

(c)

Figure 8: Transient change in distribution of choline density for simulation runs (a) R-A1 at Ez = −1.8 kB T /eRc , (b) R-B1 at Ez = −1.4 kB T /eRc , and (c) R-C1 Ez = −1.4 kB T /eRc . For comparison, dotted line shows the choline density average over simulation time at Ez =0. The other symbols are ( ) in time interval 0-2.28 ns, () in 2.28-4.56 ns, (4) in 4.56-6.84 ns, (5) in 6.84-9.13 ns, and (+) in 9.13-11.41 ns. So water penetration is prevented by the inherent electrical fields of the charges in absence of external field. However, the distribution of the dielectrophoretic force does not remain symmetric when an external electrical field is imposed across the bilayer. As indicated by the label A, the magnitude of the force under an applied electric field reduces at the anodic leaflet, while the force amplitude increases at the cathodic side, shown by the label C. This explains why there has been a greater extent of water succeeding in penetrating the core from the anodic leaflet, and supports the hypothesis of Ziegler and Vernier 53 that for a saturated lipids a pore is initiated by the water at the anodic interface.

5.4

Transient change in density distribution

As seen before in Fig. 3(Q), the lipids that make the pore hydrophilic are mainly from the anodic leaflet. This fact is further substantiated by a transient change in density distribution of choline in Fig. 8 (refer Figs. 7, 8, and 9 in SI for all repeated simulation runs). The peak values of the profiles in a time interval of 2.28 ns are summarized in Table 5. The percentage reduction in the peak density values (%∆) under poration field with reference to the peak densities in absence of applied field indicates that the density reduction following the poration

19


Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5: The peak density values of choline and their relative % change reported with time for simulation runs R-A1, R-B1 and R-C1. The change is with reference to the corresponding peak densities in absence of applied electric field. For all repeated runs refer Tables 1, 2 and 3 in SI. time (ns) 0 - 2.28 2.28 - 4.56 4.56 - 6.84 6.84 - 9.13 9.13 - 11.41

0 - 2.28 2.28 - 4.56 4.56 - 6.84 6.84 - 9.13 9.13 - 11.41

0 - 2.28 2.28 - 4.56 4.56 - 6.84 6.84 - 9.13 9.13 - 11.41 1

Density of choline in run R-A1 Anodic leaflet (%∆1 ) Cathodic leaflet (%∆1 ) −3 [0.625 Rc at E=0] [0.638 Rc−3 at E=0] 0.517 (-17.4) 0.519 (-18.6) 0.491 (-21.5) 0.555 (-13.0) 0.521 (-16.7) 0.526 (-17.4) 0.451 (-27.9) 0.512 (-19.8) 0.438 (-30.0) 0.550 (-13.8) Density of choline in run R-B1 [0.643 Rc−3 at E=0] [0.629 Rc−3 at E=0] 0.506 (-21.2) 0.625 (-0.6) 0.463 (-27.9) 0.555 (-11.8) 0.485 (-24.6) 0.589 (-6.4) 0.497 (-22.7) 0.556 (-11.6) 0.479 (-25.4) 0.562 (-10.6) Density of choline in run R-C1 [0.635 Rc−3 at E=0] [0.640 Rc−3 at E=0] 0.555 (-13.3) 0.619 (-2.5) 0.546 (-14.7) 0.654 (3.0) 0.483 (-24.5) 0.563 (-11.3) 0.403 (-37.0) 0.453 (-28.7) 0.408 (-36.2) 0.495 (-22.0)

relative % change, with standard error within ∼0.5.

is relatively higher at the anodic leaflet than the cathodic side. With density values higher in the pore for z > 0 than z < 0, the inset in Fig. 8 shows that these lipids that have left the anodic interface are now inside the pore.

5.5

Transient change in dipole moment

The dynamical change in the z-component dipole moments of lipids and water is depicted in Fig. 9 ( and in Figs. 13, 14 and 15 in SI for repeated runs) with the peak values tabulated in Table 6. As mentioned earlier, the change in z-dipole moment of lipids under electric field has dual causes; first is due to a change in density, and second is the field induced re-orientation of dipoles, which tends the lipids to align with the field direction. The reduction in density of lipids at interfaces lowers the peak values of +µz at the anodic interface and −µz at the cathodic interface. On the other hand, the field induced re-orientation always tries to reduce 20



200

R-A1

200

lipid

150

R-C1

R-B1

lipid

200

lipid

water water

100

100

100

water

0

µz [D]

µz [D]

50 µz [D]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 36

0

0

-50 -100

-100

-100

water

water

-150 -200 -5

water

lipid

-4

-3

-2

-1

0

1

2

3

4

-5

5

-4

-3

-2

lipid

-1

0

1

2

3

4

5

-5

z [nm]

z [nm]

(a)

-200

lipid

-200

-4

-3

-2

-1

0

1

2

3

4

5

z [nm]

(b)

(c)

Figure 9: Transient change in distribution of z-component dipole moments of lipid and water for simulation runs (a) R-A1 at Ez = −1.8 kB T /eRc , (b) R-B1 at Ez = −1.4 kB T /eRc , and (c) R-C1 at Ez = −1.4 kB T /eRc . For comparison, the black curves show the distributions for lipid (solid) and water (dotted), averaged over 11.41 ns at Ez = 0. The other symbols are ( ) in time interval 0-2.28 ns, () in 2.28-4.56 ns, (4) in 4.56-6.84 ns, (5) in 6.84-9.13 ns, and (+) in 9.13-11.41 ns. Table 6: The peak values of z-component of dipole moment of lipids, water and their relative % change reported with time. The change is with reference to the corresponding peak values at anodic and cathodic leaflets in absence of applied electric field. For all repeated runs refer Tables 4, 5 and 6 in SI. time (ns) 0 - 2.28 2.28 - 4.56 4.56 - 6.84

lipid dipole moment in run R-A1 µz anodic side (%∆1 1 )(%∆2 2 ) −µz cathodic side (%∆1 1 )(%∆2 2 ) [197.4 D at E=0] [196.7 D at E=0] 144.1 (-27.0)(-9.6) 177.9 (-9.5)(9.1) 137.2 (-30.5)(-8.9) 194.5 (-1.1)(11.9) 149.9 (-24.0)(-7.3) 181.6 (-7.7)(9.8)

water dipole moment in run R-A1 −µz anodic side (%∆3 ) µz cathodic side (%∆3 ) [140.9 D at E=0] [142.2 D at E=0] 157.7 (12.0) 101.8 (-28.4) 156.4 (11.0) 107.6 (-24.4) 161.8 (14.8) 99.5 (-30.0)

6.84 - 9.13 9.13 - 11.41

124.5 (-36.9)(-9.0) 134.0 (-32.1)(-2.1)

0 - 2.28

184.3 (-6.3)(13.5) 188.3 (-4.3)(9.5) lipid dipole moment in run R-B1 [223.7 D at E=0] [218.6 D at E=0] 165.6 (-26.0)(-4.7) 229.6 (5.0)(5.7)

159.6(13.3) 92.1 (-35.2) 162.3 (15.2) 92.3 (-35.1) water dipole moment in run R-B1 [100.1 D at E=0] [105.6 D at E=0] 120.8 (20.6) 99.6 (-5.6)

2.28 - 4.56 4.56 - 6.84 6.84 - 9.13 9.13 - 11.41

146.8 171.1 162.6 150.6

0 - 2.28 2.28 - 4.56

199.9 (-8.5)(3.2) 218.4 (-0.1)(6.2) 200.1 (-8.4)(3.1) 204.5 (-6.4)(4.2) lipid dipole moment in run R-C1 [235.0 D at E=0] [234.6 D at E=0] 203.2 (-13.5)(-0.3) 241.5 (2.9)(5.5) 206.9 (-12.0)(2.7) 262.6 (11.9)(8.9)

100.5 (0.3) 77.5 (-26.5) 104.4 (4.2) 81.1 (-23.2) 97.7 (-2.4) 73.1 (-30.8) 118.4 (18.2) 72.6 (-31.3) water dipole moment in run R-C1 [70.2 D at E=0] [66.2 D at E=0] 80.9 (15.3) 46.4 (-29.9) 81.2 (15.7) 54.9 (-17.1)

4.56 - 6.84 6.84 - 9.13 9.13 - 11.41

171.3 (-27.1)(-2.6) 137.3 (-41.6)(-4.5) 141.1 (-40.0)(-3.7)

70.0 (-0.2) 60.1 (-14.3) 69.7 (-0.7)

1 2 3

(-34.4)(-6.4) (-23.5)(1.1) (-27.3)(-4.6) (-32.7)(-7.3)

220.1 (-6.2)(5.2) 166.7 (-28.9)(-0.3) 187.0 (-20.3)(1.7)

41.8 (-36.9) 18.6 (-71.9) 24.8 (-62.6)

∆1 is total relative % change, within the standard error 2 ; this includes changes due to density reduction and lipid re-orientation. ∆2 is the relative % change, solely due to lipid re-orientation caused by the electric field. ∆ is the total relative % change, within the standard error 2.

+µz at the anodic interface and causes to enhance −µz at the cathodic interface. This fact is affirmed by %∆2 values being negative at anodic leaflet and positive at anodic side. A reason why more lipids come into the pore from the anodic leaflet is their z-dipole moment being positive at the anodic interface. The applied field tries to render their z21


Page 23 of 36

0.8 System-A

0.6

0.4

0 Ez = 0 Ez = -1.0 kBT/eRc Ez = -1.2 kBT/eRc Ez = -1.4 kBT/eRc Ez = -1.6 kBT/eRc Ez = -1.8 kBT/eRc

-0.2 -0.4 -0.6 -8 -7 -6 -5 -4 -3 -2 -1 0

1 2

3 4

5 6

7

0.4

0.2

0.2

0 -0.2

-0.6 8

-0.8 -8

-6

-4

-2

0

2

4

6

System-C (1 M)

0 -0.2

Ez = 0 Ez = -0.8 kBT/eRc Ez = -1.0 kBT/eRc Ez = -1.2 kBT/eRc Ez = -1.4 kBT/eRc

-0.4

z [nm]

(a)

0.6

0.4

φ(z) [v]

0.2

0.8 System-B (0.25 M)

φ(z) [v]

0.6

φ(z) [V]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ez = 0 Ez = -0.8 kBT/eRc Ez = -1.0 kBT/eRc Ez = -1.2 kBT/eRc Ez = -1.4 kBT/eRc

-0.4 -0.6 8

-0.8 -8 -7 -6 -5 -4 -3 -2 -1 0

z [nm]

(b)

1 2

3 4

5 6

7

8

z [nm]

(c)

Figure 10: Electrostatic potential distribution across the membrane with applied electric fields for (a) system-A, (b) system-B, and (c) system-C. Potential is averaged over time 11.41 ns. dipole moment negative by flipping them, as it successfully does on water dipoles in the bulk. However, restricted by the hydrophobic tail, the dipole moment of lipids reduces only partly. A way out of this frustrating situation is therefore to go into the pore. Thus, the flux of lipids into the pore is higher from the anodic leaflet. Note that the net percentage change in z-dipole moment of lipids at the cathodic leaflet is slightly positive or negative due to cancelling effects of the field induced orientation and density change. For water, the field induced change in dipole moment is pronounced on +µz at the cathodic interface, since the applied field (pointing into −z direction) strongly opposes the positive dipoles.

5.6

Electrostatic Potential

In atomistic molecular dynamics simulations, the electrostatic potential is computed via a double integration of charge density across the bilayer. However in DPD, the potential is readily available at grid points by solving the Poisson equation in real space. In absence of applied fields the overall potential inside the bilayer (due to lipids and water) is reported to be positive (∼ 1 V ) in all-atom molecular dynamics simulation works. 54 With their polarizable MARTINI model,Yesylevskyy et al. 55 observed negative (∼ −2 V ) potential. In our DPD simulations, the potential within the bilayer in absence of external field is ∼ −0.2 V , as seen in Fig. 10. 22



The transmembrane voltage (TMV) increases with the field E, as depicted in Fig. 10. Moving away from the interfaces on both the sides of membrane, the magnitude of potential continues to increase in the salt-free system, as shown in Fig. 10(a), while the potential remains almost constant for the systems B and C, which contain the salt. This is because, for the system without salt, there is a high non-zero dipole moment of water (µz = −39±2 D at Ez = −1.4 kB T /eRc ) in z-direction in bulk regions, which is reduced to around µz ∼ −9±2 in systems containing the salt (see Figs. 5(b) and 5(c)). The lipid dipoles near the interface at anodic leaflet are pointing upward in +z-direction, and they are surrounded by water dipoles orienting in −z-direction, followed by water in the bulk facing again in −z-direction under electric field. However, at the other side of the bilayer, the interfacial water and the water in the bulk have their dipoles facing opposite. This asymmetry across the membrane in configurations of the dipoles is a reason that explains an abrupt change in the potential distribution near the cathodic interface. At the porating field, the transmembrane voltage is close to 1 V , which is less than the minimum reported threshold values (2 − 2.6 V ) for phospholipids in studies with molecular dynamics simulations. 12,54 The reasons for over estimated TMV in molecular dynamics simulations have been stated by Ziegler and Vernier. 53 In our DPD simulations, the TMV value ∼ 1 V is consistent with φ = Ez Lz . For the size of simulation box in z-dimension as high as Lz = 26 Rc , the applied field necessary to set up 1 V TMV is as low as 0.06 V /nm compared to other molecular simulation works 12,14 where the applied field is as high as 0.6 − 1 V /nm on account of much smaller Lz in MD calculations. For systems B and C, as mentioned earlier, the electric field at the interfaces scales as Ez = φs κ, where κ and φs are the inverse Debye screening length and the electrostatic potential at the interface, respectively. By adding a salt the Debye length decreases, which in turn enhances the field at the interfaces. Therefore, the external field required to set up 1 V TMV is as low as Ez = −1.4 kB T /eRc in the systems containing salt as compared to the field Ez = −1.8 kB T /eRc required to impose 1 V TMV in the salt-free system A.

23


Page 24 of 36

Page 25 of 36

400

400

300

Ez = -1.4 kBT/eRc Ez = -1.6 kBT/eRc Ez = -1.8 kBT/eRc

Ez = 0 Ez = -1.0 kBT/eRc Ez = -1.2 kBT/eRc

300

System-C (1 M)

γ = 9.6 mN/m γ = 7.9 mN/m γ = 5.7 mN/m γ = 4.1 mN/m

Ez = -1.4 kBT/eRc

100

200

100

100

0

0

0

-100

-100

-100

-200 -3

-2

-1

0

1

2

-200 -3

3

-2

z [nm]

(a)

γ = 9.7 mN/m γ = 7.3 mN/m γ = 5.3 mN/m γ = 1.6 mN/m

Ez = 0 Ez = -1.0 kBT/eRc Ez = -1.2 kBT/eRc Ez = -1.4 kBT/eRc

300

200 Σ [bar]

200

400 System-B (0.25 M)

γ = 7.6 mN/m γ = 6.1 mN/m γ = 4.7 mN/m γ = 4.5 mN/m γ = 2.7 mN/m γ = 0.9 mN/m

Σ [bar]

Ez = 0 Ez = -1.0 kBT/eRc Ez = -1.2 kBT/eRc

System-A

Σ [bar]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-1

0

1

2

-200 -3

3

z [nm]

(b)

-2

-1

0

1

2

3

z [nm]

(c)

Figure 11: Stress distribution across the membrane with the applied electric field for (a) system-A, (b) system-B, and (c) system-C. The stress is averaged over simulation time 11.41 ns. The membrane tension is indicated by γ within the standard error 0.6 mN/m for the respective fields.

5.7

Stress distribution

Using the pressure tensor, namely normal component Pzz and lateral components Pxx , and Pyy , the stress profile Σ(z) =< Pzz − (Pxx + Pyy )/2 > is computed. To this end, the box length in z dimension is split into slabs of thickness 0.2 Rc . Following the Irving-Kirkwood ext arising from an external electrical field Ez is deduced method, 56 the normal component Pzz

for a patch of membrane with the lateral area A using

ext Pzz (z) = −

Ez X qi sign(z − zi ). 2A i

(15)

ext ext While the lateral or in-plane pressure component (Pxx + Pyy )/2 would simply be zero. Here

zi are z-coordinates of charges qi . With a uniform field pointing in negative z-direction, the lipids in the anodic leaflet are the most stressed because of having their dipoles oriented opposite to the direction of applied field. This is apparent in Fig. 11, which shows the stress distribution at varying fields, with an enhanced asymmetry at the anodic leaflet. The stress distribution becomes increasingly asymmetric with the applied field, and this fact is in agreement with a recent report in the

24



1.2 System-B (0.25 M), R-B1

cation anion cation + anion

4.0

conductance [nS]

1.0

conductance [nS]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 36

0.8 0.6 0.4

cation anion cation + anion

System-C (1 M) R-C1

3.0

2.0

1.0

0.2 0.0

0.0 2

3

4

5

6

7

8

9

10

2

3

4

time [ns]

(a)

5

6

7

8

9

10

time [ns]

(b)

Figure 12: Transient conductance in simulation runs (a) R-B1, and (b) R-C1 at Ez = -1.4 kB T /eRc . The symbols are ( ) cation, () anion, and (♦) both cation and anion. The vertical dotted line indicates the onset of electroporation. work of Bruhn et al., 57 which states that the membrane potential alters the stress profile. The membrane tension, which is obtained by integrating the stress profile, drops with an increasing field. This is again consistent with the effective tension γ ∗ = γ − 1/2CV 2 as per the capacitive model. 7 Referring to the stresses at Ez = −1.0, −1.2, and −1.4 kB T /eRc in Figs. 11(b) and 11(c) for the systems containing salt, and comparing them with the stresses in Fig. 11(a) of the salt-free system, the higher stress associated with the salt systems is due to the higher value of Pzz arising from the osmotic pressure of ions. This is a tension increasing action of ions. But the applied field per se has the tension reducing action. Therefore, the effect of electric field is due to a combined effect of tension reducing action of the applied field and tension increasing action of the osmotic pressure due to ions.

5.8

Ionic Conductance

Under electric potential gradients, ions exhibit electrophoretic motion; cations migrate toward a lower potential while anions tend to move toward a higher potential. Therefore, since the potential gradient imposed upon the membrane system is pointing in −z direction, the (net) flux of cation is found to be in direction from the anodic leaflet to the cathodic one,

25


Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


whereas anions migrate in the opposite direction. Furthermore, as observed earlier in Fig. 3(Q), the cumulative flux for cations is more than the counter flux of anions at all times. Since the negatively charged phosphate repels anions, there is a lower concentration of anions near the interfaces than cations (see the inset in Fig. 2), and therefore it is obvious to have anion flux lesser than the flux for cations. Using the fluxes for ions, the ionic currents (slopes of fluxes) and conductances (current divided by TMV) are calculated, and shown in Fig. 12. Since the transmembrane voltage is ∼ 1 V for the systems studied, the ordinate in Fig. 12 can be regarded as the ionic currents in nA. It is plain to see that the conductance (and current) contributed by cations is higher than the conductance by anions, and that the conductances rise with time, which is more clearly evident in the system with high salt concentration. Since the membrane pores grow with time (see the images given in SI), the rise in conductance with time qualitatively agrees with the increased conductance against pore radius reported in molecular dynamics simulation work by Rems et al. 58

5.9

Comparison with MD simulations

Table 7 summarizes the key findings from our DPD simulations, and compares them with other MD simulations. Here we see a fair agreement between the two simulation techniques in regards to the findings, and it thus successfully establishes the DPD model’s ability (with polarizable water) in capturing the electrostatics, and provides a strong base case for studying field induced phenomena in biological systems with DPD model.

6

Conclusions

To sum up, a model bilayer hydrated by polarizable version of DPD water is simulated for electroporation in presence of electrolytes. A great deal of discussion on dipolar rearrangement of lipids and water has been made. An external electrical field required to porate the membrane within 11.41 ns simulation run is found to be more in absence of electrolytes. Nevertheless, it is the intrusion of water that initiates the prepore creation in both the caseswith and without salt. The water gushing into the pore is more from the anodic leaflet 26



Page 28 of 36

Table 7: Findings summarized and compared with atomistic simulations. Item Bead density distribution Poration time

DPD simulation becomes broader, peak reduced order of few ns

Fluxes of water & lipid Mechanism

more from anodic leaflet water fingering water dipole becomes asymmetric with applied field decreases with applied field in NVT (1) reduces required poration field, (2) lipid and water dipoles change grows with time stronger for cations

Stress profile Membrane tension Salt effect

Conductance Ion fluxes

MD simulation Tieleman 12

Consistent? yes

Tieleman 12 Tarek 14 Ziegler and Vernier 53

yes yes yes

Tieleman 12 Ziegler and Vernier 53 Bruhn et al. 57

yes yes yes

Tarek 14 in NPT Tieleman 12 almost no salt effect

no

Rems et al. 58 Rems et al. 58 anion flux is more

yes no

no

on account of the dielectrophoretic force. The hydrophobic pore is stabilized by lipid head groups that come mainly from the anodic leaflet. Cations migrate from a higher electrostatic potential (anode) to a lower potential (cathode), and the flux of cations is always more than the counter flux of anions. By merit of higher cationic flux, the conductance due to cations is observed to be more than anions. The conductance grows with a pore size. The applied field induces asymmetry in stress distribution across the membrane, and the effect is more pronounced at the anodic leaflet. The membrane tension shows declining trend with the applied field. The salt reduces the poration field. By strengthening the lipid dipoles and weakening the water dipoles, the salt has effects on the interfacial structure of lipid and water. These findings are also compared with other molecular dynamics simulations.

27


Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ASSOCIATED CONTENT Supporting Information Figures for bead density profiles, species fluxes, transient density profiles, dielectrophoretic force, transient profiles for water and lipids dipole moments, and simulation snapshots. Tables for peak values for choline density and dipole moments of water and lipids. This information is available free of charge via the Internet at http://pubs.acs.org.

Acknowledgement The authors acknowledge Prof. Ateeque Malani from IIT Bombay, and DST India for providing the high performance computing facility.

References (1) Golzio, M.; Teissié, J.; Rols, M.-P. Direct visualization at the single-cell level of electrically mediated gene delivery. Proc. Natl. Acad. Sci. 2002, 99, 1292–1297. (2) Paganin-Gioanni, A.; Bellard, E.; Escoffre, J.; Rols, M.; Teissie, J.; Golzio, M. Direct visualization at the single-cell level of siRNA electrotransfer into cancer cells. Proc. Natl. Acad. Sci. 2011, 108, 10443–10447. (3) Heller, R.; Gilbert, R.; Jaroszeski, M. J. Clinical applications of electrochemotherapy. Adv. Drug Deliv. Rev. 1999, 35, 119–129. (4) Abidor, I.; Sowers, A. Kinetics and mechanism of cell membrane electrofusion. Biophys. J. 1992, 61, 1557–1569. (5) Mekid, H.; Mir, L. In vivo cell electrofusion. BBA-General Subjects 2000, 1524, 118– 130. 28



(6) Prausnitz, M. R.; Bose, V. G.; Langer, R.; Weaver, J. C. Electroporation of mammalian skin: a mechanism to enhance transdermal drug delivery. Proc. Natl. Acad. Sci. 1993, 90, 10504–10508. (7) Abidor, I.; Arakelyan, V.; Chernomordik, L.; Chizmadzhev, Y. A.; Pastushenko, V.; Tarasevich, M. Electric breakdown of bilayer lipid membranes: I. The main experimental facts and their qualitative discussion. J. Electroanal. Chem. Interfacial Electrochem. 1979, 104, 37–52. (8) Benz, R.; Beckers, F.; Zimmermann, U. Reversible electrical breakdown of lipid bilayer membranes: a charge-pulse relaxation study. J. Membr. Biol. 1979, 48, 181–204. (9) Weaver, J. C.; Chizmadzhev, Y. A. Theory of electroporation: A review. Bioelectrochem. Bioenerg. 1996, 41, 135–160. (10) Weaver, J. C. Electroporation of biological membranes from multicellular to nano scales. IEEE Trans. Dielectr. Electr. Insul. 2003, 10, 754–768. (11) Tieleman, D. P.; Leontiadou, H.; Mark, A. E.; Marrink, S.-J. Simulation of pore formation in lipid bilayers by mechanical stress and electric fields. J. Am. Chem. Soc. 2003, 125, 6382–6383. (12) Tieleman, D. P. The molecular basis of electroporation. BMC Biochemistry 2004, 5, 10. (13) Kotnik, T.; Bobanović, F.; Miklavčič, D. Sensitivity of transmembrane voltage induced by applied electric fieldsâĂŤa theoretical analysis. Bioelectrochem. Bioenerg. 1997, 43, 285–291. (14) Tarek, M. Membrane electroporation: a molecular dynamics simulation. Biophys. J. 2005, 88, 4045–4053.

29


Page 30 of 36

Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(15) Böckmann, R. A.; De Groot, B. L.; Kakorin, S.; Neumann, E.; Grubmüller, H. Kinetics, statistics, and energetics of lipid membrane electroporation studied by molecular dynamics simulations. Biophys. J. 2008, 95, 1837–1850. (16) Majhi, A. K.; Kanchi, S.; Venkataraman, V.; Ayappa, K.; Maiti, P. K. Estimation of activation energy for electroporation and pore growth rate in liquid crystalline and gel phases of lipid bilayers using molecular dynamics simulations. Soft Matter 2015, 11, 8632–8640. (17) Casciola, M.; Kasimova, M. A.; Rems, L.; Zullino, S.; Apollonio, F.; Tarek, M. Properties of lipid electropores I: Molecular dynamics simulations of stabilized pores by constant charge imbalance. Bioelectrochemistry 2016, 109, 108–116. (18) Marracino, P.; Castellani, F.; Vernier, P.; Liberti, M.; Apollonio, F. Geometrical characterization of an electropore from water positional fluctuations. J. Membr. Biol. 2017, 250, 11–19. (19) Marracino, P.; Liberti, M.; Vernier, P.; Apollonio, F. A statistical analytical model for hydrophilic electropore characterization: a comparison study. RSC Advances 2017, 7, 31997–32007. (20) Gurtovenko, A. A.; Lyulina, A. S. Electroporation of asymmetric phospholipid membranes. J. Phys. Chem. B 2014, 118, 9909–9918. (21) Dehez, F.; Delemotte, L.; Kramar, P.; Miklavčič, D.; Tarek, M. Evidence of conducting hydrophobic nanopores across membranes in response to an electric field. J. Phys. Chem. C 2014, 118, 6752–6757. (22) Polak, A.; Tarek, M.; Tomšič, M.; Valant, J.; Ulrih, N. P.; Jamnik, A.; Kramar, P.; Miklavčič, D. Electroporation of archaeal lipid membranes using MD simulations. Bioelectrochemistry 2014, 100, 18–26.

30



(23) Vaiwala, R.; Jadhav, S.; Thaokar, R. Electroporation using Dissipative Particle Dynamics with a novel protocol for applying Electric field. J. Chem. Theory Comput. 2019, 15, 603–612. (24) Vaiwala, R.; Jadhav, S.; Thaokar, R. Four-to-one coarse-grained polarisable water model for dissipative particle dynamics. Mol. Simul. 2018, 44, 540–550. (25) Hoogerbrugge, P.; Koelman, J. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett. 1992, 19, 155. (26) Kong, Y.; Manke, C.; Madden, W.; Schlijper, A. Simulation of a confined polymer in solution using the dissipative particle dynamics method. Int. J. Thermophys. 1994, 15, 1093–1101. (27) Groot, R. D.; Madden, T. J.; Tildesley, D. J. On the role of hydrodynamic interactions in block copolymer microphase separation. J. Chem. Phys. 1999, 110, 9739–9749. (28) Spenley, N. Scaling laws for polymers in dissipative particle dynamics. Europhys. Lett. 2000, 49, 534. (29) Boek, E.; Coveney, P. V.; Lekkerkerker, H.; van der Schoot, P. Simulating the rheology of dense colloidal suspensions using dissipative particle dynamics. Phys. Rev. E 1997, 55, 3124. (30) Pryamitsyn, V.; Ganesan, V. A coarse-grained explicit solvent simulation of rheology of colloidal suspensions. J. Chem. Phys. 2005, 122, 104906. (31) Chatterjee, A.; Wu, L.-M. Predicting rheology of suspensions of spherical and nonspherical particles using dissipative particle dynamics (DPD): methodology and experimental validation. Mol. Simul. 2008, 34, 243–250. (32) Jury, S.; Bladon, P.; Cates, M.; Krishna, S.; Hagen, M.; Ruddock, N.; Warren, P.

31


Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Simulation of amphiphilic mesophases using dissipative particle dynamics. Phys. Chem. Chem. Phys. 1999, 1, 2051–2056. (33) Anderson, R. L.; Bray, D. J.; Del Regno, A.; Seaton, M. A.; Ferrante, A. S.; Warren, P. B. Micelle formation in alkyl sulfate surfactants using dissipative particle dynamics. J. Chem. Theory Comput. 2018, 14, 2633–2643. (34) Kranenburg, M.; Venturoli, M.; Smit, B. Molecular simulations of mesoscopic bilayer phases. Phys. Rev. E 2003, 67, 060901. (35) Kranenburg, M.; Smit, B. Phase Behavior of Model Lipid Bilayers. J. Phys. Chem. B 2005, 109, 6553–6563. (36) Foram, M. T.; Ayappa, K. G. Effect of Polymer Grafting on the bilayer Gel to LiquidCrystalline Transition. J. Phys. Chem. B 2010, 114, 2738–2748. (37) Winterhalter, M. Lipid membranes in external electric fields: kinetics of large pore formation causing rupture. Advances in colloid and interface science 2014, 208, 121– 128. (38) Delemotte, L.; Dehez, F.; Treptow, W.; Tarek, M. Modeling membranes under a transmembrane potential. J. Phys. Chem. B 2008, 112, 5547–5550. (39) Li, Z.; Bian, X.; Caswell, B.; Karniadakis, G. E. Construction of dissipative particle dynamics models for complex fluids via the Mori–Zwanzig formulation. Soft Matter 2014, 10, 8659–8672. (40) Yoshimoto, Y.; Kinefuchi, I.; Mima, T.; Fukushima, A.; Tokumasu, T.; Takagi, S. Bottom-up construction of interaction models of non-Markovian dissipative particle dynamics. Phys. Rev. E 2013, 88, 043305. (41) Kinjo, T.; Hyodo, S. Equation of motion for coarse-grained simulation based on microscopic description. Phys. Rev. E 2007, 75, 051109. 32



(42) Akkermans, R. L.; Briels, W. J. Coarse-grained dynamics of one chain in a polymer melt. J. Chem. Phys. 2000, 113, 6409–6422. (43) Español, P.; Warren, P. Statistical mechanics of dissipative particle dynamics. Europhys. Lett. 1995, 30, 191. (44) Groot, R. D. Electrostatic interactions in dissipative particle dynamics - simulation of polyelectrolytes and anionic surfactants. J. Chem. Phys. 2003, 118, 11265–11277. (45) González-Melchor, M.; Mayoral, E.; Velázquez, M. E.; Alejandre, J. Electrostatic interactions in dissipative particle dynamics using Ewald sums. J. Chem. Phys. 2006, 125, 224107. (46) Vaiwala, R.; Jadhav, S.; Thaokar, R. Electrostatic interactions in dissipative particle dynamics- Ewald-like formalism, error analysis, and pressure computation. J. Chem. Phys. 2017, 146, 124904. (47) Bore, S. L.; Kolli, H. B.; Kawakatsu, T.; Milano, G.; Cascella, M. Mesoscale Electrostatics Driving Particle Dynamics in Nonhomogeneous Dielectrics. J. Chem. Theory Comput. 2019, 15, 2033–2041. (48) Beckers, J.; Lowe, C.; De Leeuw, S. An iterative PPPM method for simulating Coulombic systems on distributed memory parallel computers. Mol. Simul. 1998, 20, 369–383. (49) Li, X.; Gao, L.; Fang, W. Dissipative Particle Dynamics Simulations for Phospholipid Membranes Based on a Four-To-One Coarse-Grained Mapping Scheme. PloS One 2016, 11, e0154568. (50) Casciola, M.; Tarek, M. A molecular insight into the electro-transfer of small molecules through electropores driven by electric fields. BBA-Biomembranes 2016, 1858, 2278– 2289.

33


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(51) Groot, R. D.; Warren, P. B. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys. 1997, 107, 4423–4435. (52) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33–38. (53) Ziegler, M. J.; Vernier, P. T. Interface water dynamics and porating electric fields for phospholipid bilayers. J. Phys. Chem. B 2008, 112, 13588–13596. (54) Delemotte, L.; Tarek, M. Molecular dynamics simulations of lipid membrane electroporation. J. Membr. Biol. 2012, 245, 531–543. (55) Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable water model for the coarse-grained MARTINI force field. PLOS Comput. Biol. 2010, 6, e1000810. (56) Irving, J.; Kirkwood, J. G. The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics. J. Chem. Phys. 1950, 18, 817–829. (57) Bruhn, D. S.; Lomholt, M. A.; Khandelia, H. Quantifying the relationship between curvature and electric potential in lipid bilayers. J. Phys. Chem. B 2016, 120, 4812– 4817. (58) Rems, L.; Tarek, M.; Casciola, M.; Miklavčič, D. Properties of lipid electropores II: Comparison of continuum level-modeling of pore conductance to molecular dynamics simulations. Bioelectrochemistry 2016, 112, 112–114.

34



Table of Contents (TOC) graphic

35


Page 36 of 36

Establishing an electrostatics paradigm for membrane

Recommend Documents