Electroporation Using Dissipative Particle Dynamics with a Novel

Subscriber access provided by YORK UNIV

Condensed Matter, Interfaces, and Materials

Electroporation using Dissipative Particle Dynamics with a novel protocol for applying Electric field Rakesh Vaiwala, Sameer Jadhav, and Rochish Thaokar J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00911 • Publication Date (Web): 11 Dec 2018 Downloaded from http://pubs.acs.org on December 19, 2018

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 26 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

Journal of Chemical Theory and Computation

Electroporation using Dissipative Particle Dynamics with a novel protocol for applying Electric field 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 In molecular dynamics simulations of membrane electroporation, the bilayer is subjected to an electric field E either by direct addition of a force f = qE on the chargebearing species or by imposing an ion imbalance in the salt solutions on the two sides of the bilayer. The former is believed to mimic electroporation with high fields over nanosecond pulse period, during which the membrane is almost uncharged especially in the low salt limit. Conversely, the ion imbalance method elucidates a low electric fieldinduced poration over a longer period of micro-to-milli seconds with a fully charged membrane. Both these methods of applying electric field have disadvantages while investigating electroporation using dissipative particle dynamics (DPD) simulations. The method involving direct addition of force fails to address the presence of a non-uniform dielectric background for ions embedded in non-polarizable DPD water and those found in the core of bilayer. The ion imbalance method in DPD simulations suffers from its unavoidable use of a wall potential to prevent the movement of ions across the periodic boundaries. To address the above issues, we propose a simple method for imposing a

1

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

desired transmembrane potential (TMV) by placing oppositely but uniformly charged plates on either side of the bilayer. Our DPD simulations demonstrate that the profiles for bead density, mechanical stress, electrical potential as well as the transient responses in dipole moment and species fluxes obtained from the proposed method utilizing charged plates are quite similar to those obtained using the ion imbalance method. The proposed protocol is free from the aforementioned drawbacks of the direct force addition and ion imbalance methods.

Introduction Biological cells have their cytosolic contents segregated from the extracellular environment by the plasma membrane. Major constituents of the cell membrane comprise of amphiphilic phospholipids, sphingolipids and steroids. The membrane offers a barrier of around 3-5 nm thickness to cellular transport. Protein molecules also associate with the cell membrane and the extent of their association ranges from adsorption on the bilayer surface to complete penetration across the bilayer. Organization of lipids and proteins within the membrane facilitates transport of several species including ions, polar and non-polar molecules and saccharides. Electroporation is a technique that enhances membrane permeability by causing a defect in a membrane when subjected to electric fields. The imposed electric field disrupts the integrity of the membrane, whereby a pore is induced in it. Pore formation facilitates transport of small molecules, such as genes, antibodies, drugs etc. through the membrane. Electroporation has therefore been employed for biomedical applications, for instance, in vitro delivery of DNA plasmids 1 and siRNA into cells, 2 as well as in therapeutic cancer treatment 3 to enhance the uptake of anti-cancer drugs such as bleomycin and cisplatin to cancer cells. The technique helps to fuse the membranes of adjacent cells in electrofusion, and is also useful in transdermal applications to enhance the transport of ions and molecules across a skin by voltage pulses. 4 Experimental evidences of aqueous pore formation in lipid bilayers under applied elec2


Page 2 of 26

Page 3 of 26 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


tric fields have been reported. 5–9 An electric field across a membrane induces a transmembrane voltage (TMV). A single or multiple pores can be observed when the TMV exceeds a threshold, which depends on the membrane constituents. Molecular simulations have been employed to unravel a mechanism of pore formation. 10,11 Accordingly, a defect in membrane is initiated by intrusion of water molecules into the bilayer core, leading to water wires spanning across the bilayer. The lipid head groups eventually try to align alongside the water wires to stabilize the pore. In experiments, a desired TMV can be set across the membrane by applying either electric pulses, linearly rising voltage signals 12 or constant current signals. 13 It has been demonstrated by Koronkiewicz and Kalinowski 13 that application of a constant current across the bilayer enables survival of long-lived pores of fluctuating size. Applying a linearly increasing current is another experimental protocol employed in the work of Kramar et al. 14 wherein they argued that the technique facilitates more realistic estimates of size, conductance and life-time of electropores. In molecular simulations, TMV has been induced by adding a force Fi = qi E to all the atoms that have charge qi . The vector E is regarded as an applied electric field. This direct way of applying electric fields is believed to model the effect of nanosecond pulses of megavolt-per-meter fields. An issue with this protocol is an artefact arising from a combined use of periodic boundary conditions and Particle-Mesh-Ewald method to solve electrostatics. 15 Moreover, the method assumes all charges to be present in a medium of uniform dielectric constant, which is not true for bilayer-in-water systems. An alternate method for applying electric field in simulations has been proposed by Sachs et al. 16 , which supposedly sets up a more realistic TMV across the membrane and mimics the effect of low voltage pulses for milli-to-micro seconds. The simulated system in their study consists of three aqueous salt baths separated by two model membranes. By setting a charge imbalance in the salt-water baths across the membranes, electric field gradients are induced as a result of explicit ion dynamics. However, the technique becomes

3



computationally more demanding due to inclusion of an additional bilayer. Another issue with this charge imbalance protocol is the difficulty in estimating the lateral tension in the resulting asymmetrical system. 14 With an aim to avoid the over-cost of a large double bilayer system in simulations, Delemotte et al. 17 have devised a variant of the ion imbalance protocol just described above. According to their technique, a single bilayer is sandwiched between two electrolytic baths having vacuum slabs at their distal faces. The vacuum slabs prevent the movement of ions and water molecules between the aqueous baths, while maintaining periodic boundary conditions in the direction normal to the bilayer. Similar to the technique by Sachs et al. 16 , the voltage gradients across the bilayer are generated by concentration difference of electrolytes between the two aqueous baths. By showing a linear relationship between the ion imbalance and the TMV thereof, Delemotte et al. 17 have demonstrated that the bilayer behaves as a capacitor. In another work by Delemotte and Tarek 18 , it has been demonstrated by molecular dynamics simulations that at the scale of the membrane, the performances of the two protocols, namely the direct way of applying electric field by adding a force Fi = qi E to all charge bearing species and the ion imbalance protocol of Delemotte et al. 17 , are very similar. Nevertheless, once a pore is formed, the ion imbalance method fails to maintain the transmembrane voltage, because the ions move across the membrane through the pore and thus diminish the concentration difference. This leads to a pore collapse or resealing. However, this can be avoided by adopting a swapping algorithm developed by Kutzner et al. 19 to maintain a constant TMV. According to the algorithm, the number of ions of each kind in two solution baths is frequently computed, and if the numbers differ from their initial counts, an ion from one solution bath is exchanged with a water molecule from the other solution bath. Using the ion imbalance protocol augmented with the swapping algorithm, in their recent work of molecular dynamics simulations, Casciola et al. 20 have shown that the pores relax to stable conformations and enable one to make reliable estimates of pore size,

4


Page 4 of 26

Page 5 of 26 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


ionic currents, conductance, and ion selectivity. In the present work, we propose a novel way for applying electric fields across a model membrane in molecular simulations. The charge plates with opposite surface charge densities are placed on either side of the membrane, near two opposite extremes of the simulation box. Such a simulation setup thus closely resembles a typical experimental setup for electroporation wherein the electric field is applied by inserting two electrodes across the membrane. We have investigated the pore formation phenomenon on a model membrane using dissipative particle (DPD) dynamics simulations with normal DPD water, which is non-polarizable and emulates the electrostatics via its dielectric constant. The long-range electrostatic interactions are computed by solving the Poisson equation in real space by successive over relaxation method. 21 The suggested protocol of imposing an electric field is compatible with the way electrostatic interactions are solved. For comparison, we have also performed DPD simulations, employing the ion imbalance protocol to set up the electric field gradients.

Dissipative particle dynamics A system in dissipative particle dynamics is a set of coarse-grained (CG) beads called dissipative particles, wherein each bead is a cluster of atoms. These dissipative particles interact with each other through three kinds of forces within a cutoff radius Rc , namely a conservative force (f C ), a dissipative force (f D ) and a stochastic force (f R ). The conservative force arises from excluded volume interactions, and is usually modelled by a quadratic soft repulsive potential. The dissipative force tries to minimize the differences in radial component of velocities between the particles, while the stochastic force adds the kinetic energy into the system to compensate for the loss in kinetic energy due to dissipation. These forces act pairwise on the particles and are directed along the vectors joining their centers. By the interplay between the dissipative and stochastic forces, which respectively act as sink and source of momentum, the system’s energy fluctuates, and thus the temperature of system remains constant. Therefore, the dissipative and random forces are the thermostating

5



Page 6 of 26

forces, for they setup a thermostat. The conservative force between two DPD particles i and j, respectively positioned at ri and rj , is fijC

rij eij . = aij 1 − Rc

(1)

Here rij is |ri − rj |, and eij represents a unit vector joining the position vectors of particles i and j. The equations for dissipative and random forces are

fijD = −γ wD (eij · vij ) eij

(2a)

fijR = σwR ζij eij .

(2b)

The strength of the dissipative force is governed by the friction parameter γ, while the parameter σ controls the noise amplitude. The vector vi denotes the instantaneous velocity, and the relative velocity vij = vi − vj . The random number ζij is drawn from a Gaussian r

distribution. The customary choice for the weight function wR (r) is 1 − Rijc . The fluctuationdissipation theorem enunciated by Español and Warren 22 within the framework of DPD ensures the equilibrium phase space distribution of the system to be the Gibbs-Boltzmann distribution. Accordingly, the weight functions and the strengths of thermostating forces obey 2 σ 2 = 2γkB T and wD = wR .

We have used γ value 4.5

√

(3)

kB T m/Rc . There are additional connector forces for the

bonded particles due to the bond-length and the angle-bending (three-body interaction) constraints. We have used the harmonic potential function with spring constant ks for maintaining the bond-lengths (l0 ) within the lipids, 1 2 uB ij = ks (rij − l0 ) , 2

6


(4)

Page 7 of 26 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


and a cosine function for retaining the stiffness of the lipid chains,

uA ijk = kθ [1 − cos(θijk − θ0 )].

(5)

Here kθ denotes the stiffness parameter, while θ0 is an equilibrium angle between three consecutive particles i, j and k. The charge species, for instance explicit ions and the charges of lipids, also interact through a long-range electrostatic force f E , which is delineated in the next section.

Electrostatics in the framework of DPD The non-bonded DPD particles overlap by dint of the soft repulsive interaction, allowing two oppositely point-charged particles to artificially form an ion-pair. To address this issue of ionpair formation, Groot 21 suggested that the ion-pairing can be avoided if the potential well is made shallower than the thermal energy (−kB T ) by smearing of the charges. The recipe to solve electrostatic forces involves spreading of charges on lattices by a linear function, determining the field gradients and divergence at the lattice nodes, and back interpolating the electrostatic forces from lattice to the charge species. González-Melchor et al. 23 employed an exponential charge distribution (Slater type density), and solved the electrostatic force using Ewald summation. A full description of the electrostatics of smeared out charge species with Slater type distribution can be found in our earlier work, Vaiwala et al. 24 . We have adopted a linear charge distribution for solving electrostatics in the present work. P Given a spatial distribution of the charges qi ( i qi = 0) by a charge density ρc (r), the distribution of electrostatic potential φ(r) will be governed by the Poisson equation

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

(6)

where 0 and r are the dielectric permittivity of vacuum and local relative dielectric constant, respectively. The local permittivity (r) can be expressed in terms of the dielectric constant 7



Page 8 of 26

of water w and the local polarizability relative to water, P (r) = r (r)/w as,

(r) = 0 r (r) = 0 w P (r).

(7)

Equation (6) in its dimensionless form reads as

∇ · (P (r)∇φ(r)) = −Γρ(r).

(8)

Henceforth, φ, ρ and the differential operator ∇ are dimensionless. The factor Γ is the non-dimensional Bjerrum length, e2 /(0 w kB T Rc ). The electrostatic potential is scaled by kB T /e, and length is made dimensionless by Rc . The symbols e, kB and T denote the elementary charge of an electron, the Boltzmann constant, and temperature of the system, respectively. The computational domain is discretized into cells of uniform linear size Rc . The point-charges are spread out over the lattice nodes within a distance of Re = 1.6Rc using a linear charge distribution function, 21 3 f (r) = πRe3

r 1− . Re

(9)

The field for electrostatic potential is solved in real space by the successive over relaxation technique, 25 whereby the nodal values of the potential are iteratively computed via

φnew = φold − ζ (Γ¯ ρe + ∇ · (P ∇φ))

(10)

till the relative error defined by equation (11) attains a value below a tolerance limit. ρ¯e indicates the total charge on a lattice node by contributions from all charged DPD particles within the distance Re . The relaxation parameter ζ is 0.15. sP err =

+ ∇ · (P ∇φ)]2 2 nodes [∇ · (P ∇φ)]

[Γ¯ ρe nodes P

8


(11)

Page 9 of 26 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


Using the same distribution function f (r), the electrostatic force on the charged DPD particles is calculated by back interpolating the grid-based gradient of potential,

fiE (ri ) = −qi

X

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

(12)

nodes

The summation is restricted to all lattice nodes which are within the distance Re from the particle location ri . Since a typical lipid bilayer has a dielectric constant ∼ 2, the polarizability (P = r /w ) for DPD particles constituting the lipid tails is taken to be 0.025, while that of water is 1. The local polarizability along any edge of cells is calculated by averaging the polarizability of particles over the four cells that share the common edge. 21

Simulation details A fluid bilayer, composed of 1108 molecules of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipids and fully hydrated by 42500 DPD water, is investigated for pore formation in 27 × 27 × 24Rc3 size simulation box at room temperature. The periodic boundary conditions are applied in all the three directions. The coarse-grained DPD water effectively represents a cluster of four water molecules. The length scale Rc is derived by mapping volume of a DPD water particle with that of four water molecules at room temperature. Accordingly, the length scale turns out to be 0.711 nm. The time scale τc is determined by mapping the diffusion coefficient of water at room temperature, 26 and it computes to be 90.861 ps. A four-to-one coarse-grained model of DMPC and its interaction parameters were taken from the recent work of Li et al. 27 . A schematic representing the model DMPC dealt in Figure 1 shows the lipid constituents and the bond-lengths between various kinds of coarse-grained (CG) beads. The bonds in the lipids are modelled by harmonic springs through spring constant 512 kB T /Rc2 , and the stiffness in lipid chains is maintained using kθ = 20 kB T . Referring Figure 1, the equilibrium bond-angles (θ0 ) for CG beads 2 − 3 − 4, 3 − 4 − 5, 4 − 5 − 6, 7 − 8 − 9 and 8 − 9 − 10 are kept 180o , whereas the triplet of CG beads 2 − 3 − 7 is constrained to 120o . The repulsive 9



CH3 H3C

1 Choline

CH3

N

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 26

0.59

5

9

0.59

0.59 6 10

Figure 1: A schematic of DMPC model used for lipid bilayer simulations. The bond-lengths are in units of Rc , and the numbers within coarse-grained beads are bead labels. Choline bead is positively charged (+1e) while the phosphate carries negative (−1e) charge.

force parameters between various interacting beads are given in Table 1. In a simulation setup shown in Figure 2(A), two charged plates with oppositely charged surface densities ±σ are positioned at the two extreme ends of the simulation box at z = ±11Rc . In another setup shown in Figure 2(B), the water baths surrounding the bilayer contain a symmetric 1:1 electrolyte. An ion imbalance (∆Q) across the bilayer is created by keeping an excess cation in the upper solution bath and an excess anion in the lower solution bath. The total number of electrolytic ions (cations and anions) are 1024. This amounts to a salt concentration of ∼ 170 mM. The mixing of the two solution baths is avoided by use of phantom walls at their distal ends. The phantom walls exert a soft repulsive force Table 1: Repulsive force constant aij for DMPC lipid, CG water (W) and electrolytic ions in reduced DPD units. aij (kB T /Rc ) Choline Phosphate Carbonyl Tail W Ions

Choline 100 100 102 130 98 98

Phosphate 100 100 102 130 98 98 10

Carbonyl Tail W 102 130 98 102 130 98 100 110 102 110 100 130 102 130 100 102 130 100


Ions 98 98 102 130 100 100

Page 11 of 26 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


(a)

(b)

Figure 2: Simulation setups; (a) setup-A for the proposed charged plates protocol, and (b) setup-B for the ion imbalance method. f = 21 Awall (1 − ∆z/zc ) in z-direction within the range zc = Rc . The strength of the force Awall = 100 kB T /Rc is sufficient to avoid migration of the particles across the boundaries in z-direction. The box length in z-direction is extended to 48 Rc by creating vacuum slabs. With this setup, the swapping algorithm 19 is adopted to prevent the drop of electric potential caused by movement of the ions across the pores during electroporation. After initial equilibration without external electric field for few nanoseconds, an external electric field is applied by the charged plates in setup-A (Figure 2(A)) and by the ion imbalance in setup-B (Figure 2(B)) for a period of 109 ns. In both the setups, the applied field points in the −z direction, with higher electric potential at the upper leaflet than the lower leaflet. Therefore, the upper leaflet is now regarded as an anodic interface and the lower leaflet as a cathodic interface.

Results and discussion In this section we compare the results obtained at various surface charge densities σ = 0, 1, 2 and 2.5 e/Rc2 in the proposed protocol (setup-A) with those derived from the ion imbalance

11



protocol (setup-B) at imbalances of ∆Q = 0, 100, 200 and 300 e. The electric field set up by the charged plates depends on the surface charge density, relative position of the plates and the local polarizability (for pure water column refer Appendix-A). Pore formation is observed in the respective setups at σ = 2.5 e/Rc2 and ∆Q = 300 e. The lower field values in setup-B are due to the presence of free ions. The setup-A should be relevant for nanosecond pulsed electroporation wherein the dielectric response is expected; while the setup-B should simulate a micro-to-milli seconds pulse wherein the membrane is fully charged by free ions.

Density profiles In Figure 3 the density profiles for choline and tail beads are compared. The structure of choline in Figure 3(A) corresponding to the case of no pore formation at σ ≤ 2 e/Rc2 as well as that of pore formation at σ = 2.5 e/Rc2 is fairly identical with the respective density distributions in setup-B at various values of ∆Q shown in Figure 3(B). The distribution becomes broader with an increase in the imposed electric field, and this is in line with the density distribution reported by Tieleman 11 . The effect is more pronounced at the porating field. Similar is the case with the density profile for tail beads indicated in Figures 3(C) and 3(D), where the distribution becomes wider at σ = 2.5 e/Rc2 in setup-A and at ∆Q = 300 e in setup-B. The spread of ion concentrations for ∆Q = 0 and 300 e is depicted in Figure 4. For no imbalance of ions across the membrane, the distributions for cation and anion are symmetric across the bilayer, with more cations near the interfaces than the anions. Moreover, the cations have penetrated deeper into the bilayer core, whereas the anions are repelled by the phosphate groups, which restrict them from further spanning into the core of bilayer. At ∆Q = 300 e, it is plain to see that the gradients of cations and anions are established, which are responsible for setting up a gradient of electric potential.

12


Page 12 of 26

Page 13 of 26

(A)

-3

-3

0.5

0.4 0.3 Electroporation observed

0.2

0.4 0.3 Electroporation observed

0.2

0.1 0 -5

∆Q = 0 ∆Q = 100 ∆Q = 200 ∆Q = 300

0.6

ρ(z) [Rc ]

0.5

ρ(z) [Rc ]

(B)

σ=0 σ=1 σ=2 σ = 2.5

0.6

0.1 -4

-3

-1

-2

0

1

2

3

4

0 -5

5

-4

-3

-1

-2

Z [nm]

2

3

(D)

σ=0 σ=1 σ=2 σ = 2.5

4

5

2.5 -3

2 1.5 Electroporation observed

1

∆Q = 0 ∆Q = 100 ∆Q = 200 ∆Q = 300

3

ρ(z) [Rc ]

2.5 -3

1

3.5 (C)

3

2 1.5

Electroporation observed

1

0.5 0 -5

0

Z [nm]

3.5

ρ(z) [Rc ]

0.5 -4

-3

-1

-2

0

1

2

3

4

0 -5

5

-4

-3

-1

-2

Z [nm]

1

0

2

3

4

5

Z [nm]

Figure 3: (A) and (B) respectively refer to the density profiles for choline in setup-A and setup-B, while (C) and (D) indicate the density profiles for tail beads of lipids in setup-A and setup-B, respectively. The values for σ are in units of e/Rc2 , while those of ∆Q are in units of e. 0.05

0.04

-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


0.03

0.02 cation, ∆Q = 0 anion, ∆Q = 0 cation, ∆Q = 300 anion, ∆Q = 300

0.01

0 -7 -6

-5 -4

-3 -2

-1

0

1

2

3

4

5

6

7

Z [nm]

Figure 4: Distributions for ion concentrations across the membrane at (−) ∆Q = 0, and (−−) ∆Q = 300 e. Colors refer to (blue) cation and (pink) anion densities. The profiles are averaged over the time 72.69 ns.

13



3000

(C)

(A) counts

counts

600

400

200

2000

1000

0

0 30

40

50

60

70

80

90 100 110

30

40

50

time [ns]

80

90 100 110

(D)

(B)

300

100 50

200 100

0 30

70

400

counts

150

60

time [ns]

200

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 26

0 40

50

60

70

80

90 100 110

30

40

50

time [ns]

60

70

80

90 100 110

time [ns]

Figure 5: Transient response in species within the bilayer core and their inflows through bilayer leaflets in setup-A at σ = 2.5 e/Rc2 . (A) number of water particles (pink) in the core of bilayer, (B) number of choline (green) within the bilayer core, (C) cumulative inflow of water through anodic (pink) and cathodic (blue) leaflets, and (D) cumulative inflow of choline at anodic (green) and cathodic (brown) leaflets. The onset of poration is marked by a vertical line.

Species fluxes The transient dynamics of water and choline at the porating field (σ = 2.5 e/Rc2 ) is shown in Figure 5. A sharp increase in the water counts within a core region of the bilayer (−0.6Rc < z < 0.6Rc ) is a signature of the onset of pore creation, which occurs at ∼80 ns. This is the time when the counts for lipids within the core also rise abruptly. In absence of water charges, a role of lipids in forming a pore should be prominent. This hypothesis of lipid initiative mechanism can be confirmed by visiting the molecular images generated by visual molecular dynamics software. 28 In the image-A of Figure 7(a) it is distinctly clear that a significant amount of lipids with their head groups is present even at the onset of poration. 14


Page 15 of 26

500 3000

(A)

(C) 2000

300

counts

counts

400

200

1000

100 0

0 0

10

20

30

40

50

60

70

0

10

20

time [ns]

30

40

50

60

70

50

60

70

time [ns] 800

200

(B)

(D) 600

counts

150

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


100

400

50

200

0

0 0

10

20

30

40

50

60

70

0

10

time [ns]

20

30

40

time [ns]

Figure 6: Transient response in species within the bilayer core and their inflows through bilayer leaflets for setup-B at ∆Q = 300 e. (A) number of water particles (pink) in the core of bilayer, (B) number of choline (green), cation (violet) and anion (orange) in the core, (C) cumulative inflow of water through anodic (pink) and cathodic (blue) leaflets, and (D) cumulative inflow of choline at anodic (green) and cathodic (brown) leaflets. Also shown are the inflows of cations (violet and black) and anions (cyan and orange) with colors respectively indicating the inflows at anodic and cathodic leaflets. The onset of poration is marked by a vertical line. In (B) and (D) the lipid counts are shifted by 50 for visual clarity.

On the scale of (C) and (D) of Figure 5, the inflows of water into the pore from anode and cathodic interfaces are comparable, whereas more lipids are flowing in from the anodic leaflet at most of the instances. This is clearly due to the unfavourable dipole moment of the lipids present in the anodic leaflet. Figure 6 depicts the transient dynamics of water, choline and ions through the core region −0.6Rc < z < 0.6Rc of the bilayer in setup-B at ∆Q = 300 e. A sudden rise in water and choline counts is an indication of onset of pore formation, which occurs at ∼40 ns time. The

15



(a)

(b)

Figure 7: Simulation snapshots in (a) setup-A with σ = 2.5 e/Rc2 and (b) setup-B with ∆Q = 300 e. The beads shaded with violet and orange are head groups and tails of lipids, respectively, while cyan color beads are water particles. Cations and anions are in black and red colors, respectively. For visual clarity, water is not shown in top views, and only choline and water are shown in side views. In top view of image-A in (b), head groups are not shown for pore visualization. inflow of water through anodic leaflet is higher at all times, whereas the choline fluxes are comparable at two leaflets. Substantiated by Figure 7(b), image-A, the presence of lipids inside the core region and their forming a channel across the leaflets at the onset of pore creation itself is a proof positive that the pore formation is a lipid initiative process. For the ion fluxes, it is clear in Figure 6(D) that a net flow of cations is always from the anodic interface, while the anions migrate in the opposite way with a number quite less than that of the cations. Figure 6(D) also indicates a huge imbibition of cations at the anodic leaflet. As mentioned earlier, the anions are repelled by phosphate groups, and therefore the flux of anions is expected to be lower than the flux of cations.

16


Page 16 of 26

Page 17 of 26

0.8

0.8

(A)

0.6

0.6

0.4

0.4 φ(z) [v]

φ(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


0.2 0 -0.2

(B)

0.2 0 -0.2

σ=0 σ=1 σ=2 σ = 2.5

-0.4 -0.6

∆Q = 0 ∆Q = 100 ∆Q = 200 ∆Q = 300

-0.4 -0.6

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

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

Z [nm]

Z [nm]

Figure 8: Distribution of electrostatic potential across the membrane in (A) setup-A, and (B) setup-B. The values for σ are in units of e/Rc2 , while those of ∆Q are in units of e. The vertical dash lines indicate locations of interfaces.

Electrostatic potential Figure 8 depicts the distribution of electric potential in direction normal to the bilayer. With an increase in the surface charge density in setup-A, the potential gradient across the bilayer shows proportionate rise in TMV. This is congruous with the observation in Figure 8(B) of higher TMV at high ∆Q in setup-B. It is interesting to note that in both the setups the minimum TMV required for pore formation is ∼ 1 V. In Figure 8(A), the potential distribution within water baths at low surface charge density is nearly uniform. However, at high applied fields the potential increases in magnitude while moving away from lipid-water interfaces. This is due to an enhanced field in the ion-free water near the charged plates with increasing surface charged density. It is worth to note that the relative distribution of the field (potential drop) between the aqueous slab and the lipid bilayer is proportional to the dielectric constants of water and lipids. Due to high conductivity of the aqueous salt bath, the distribution of ions renders the potential in the bulk regions (Figure 8(B)) nearly uniform for all cases of ion imbalances studied.

17



800

800

(A)

600

600

400

400

200

200

0 t = 45.43 ns t = 54.52 ns t = 63.60 ns t = 72.69 ns t = 81.77 ns t = 90.86 ns t = 99.95 ns t = 109.03 ns

-200 -400 -600 -800

µz [D]

µ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 18 of 26

(B)

0 t = 9.08 ns t = 18.17 ns t = 27.26 ns t = 36.34 ns t = 45.43 ns t = 54.52 ns t = 63.60 ns t = 72.69 ns

-200 -400 -600 -800

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

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

Z [nm]

Z [nm]

Figure 9: Component of dipole moment of lipids in direction normal to bilayer with time, in (A) setup-A at σ = 2.5 e/Rc2 and (B) setup-B at ∆Q = 300 e.

Lipid dipoles On application of the electric field across the bilayer, the lipid dipoles reorient themselves to try and minimize the local stresses. This transient change in the z-component of dipole moment is shown in Figure 9. The data for each curve is the lipid dipole moment averaged over 100 configurations sampled over the time period 1.817 ns. Since non-polarizable water particles do not have explicit charges to shield the lipid dipoles from the imposed electric stresses, the lipid dipoles have to play a key role in sustaining their equilibrium structure at the lipid-water interfaces. The role of water is thus indirect through the electric field distribution via its dielectric constant. Since the applied field points in the −z direction, the lipid dipoles at anodic interface are unfavourably oriented, while the dipoles at the cathodic interface are favourably aligned with the direction of the applied field. Referring Figure 9(A), though the dipole moment fluctuates for an initial period of t < 72.69 ns, the imposed electrical stress eventually succeeds in lowering the dipole moments by causing structural defects, leading to creation of a pore in the bilayer with permeation of both lipids and water through it. Once the pore is formed, the dipole moment continues dropping, implying the growth of the pore with time. 18


Page 19 of 26

6

5 pore radius [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


4

3

2

1

0

30

40

50

60

70

80

90

100

110

120

time [ns]

Figure 10: Evolution of pore size with time. Circles correspond to the pore size in setup-A with σ = 2.5 e/Rc2 , while square symbol indicates pore size in setup-B at ∆Q = 300 e. Similar observation in regards to mechanism of pore formation emerges from setup-B, as depicted by the change in dipole moment dealt in Figure 9(B). The fluctuation of the dipole moment near the interface at initial times may be due to a field driven readjustment and the local screening of the head-groups by ions, which are dynamic. Once the imposed stress overcomes the lipid readjustments and thus succeeds in lowering the lipid dipole moment below a threshold value, the membrane breaks down into a pore. These observations of transient change in dipoles of lipids taken together with the simulation snapshots (Figure 7) lead us to conclude that the electroporation in absence of explicit charges assigned to water is a lipid initiated process. The poration itself seems to initiate at the anodic side.

Pore size Figure 10 shows the kinetics of pore growth. A pore is modelled as a two dimensional ellipse, with radius of the pore defined as a geometrical mean of major and minor axes of the elliptical pore. We analyze the pore radius by projecting the membrane’s configurations on the x-y plane normal to the membrane. A notable difference in performances of the two setups under comparison is that the onset of pore formation is delayed to ∼ 82 ns in the

19



400

400 (A)

γ = 14.10 mN/m γ = 14.09 mN/m γ = 12.33 mN/m γ = 10.63 mN/m

σ=0 σ=1 σ=2 σ = 2.5

300

(B)

γ = 16.04 mN/m γ = 15.51 mN/m γ = 13.53 mN/m γ = 10.67 mN/m

Σ [bar]

200

100

100

0

0

-100

-100

-200 -4

∆Q = 0 ∆Q = 100 ∆Q = 200 ∆Q = 300

300

200 Σ [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

Page 20 of 26

-3

-2

-1

0

1

2

3

-200 -4

4

-3

-2

Z [nm]

-1

0

1

2

3

4

Z [nm]

Figure 11: Stress distribution across the bilayer in (A) setup-A and (B) setup-B. The tension is indicated by γ within the standard error of 0.14 mN/m. The values for σ are in units of e/Rc2 , while those of ∆Q are in units of e. charged plates method. This may be due to charging time of the membrane, for higher the conductivity lesser is the charging time. The free ions in the charge imbalance protocol quickly accumulate at the membrane-water interfaces, enabling the poration within ∼40 ns. While the absence of free ions in the charged plates method delays electroporation. The process of pore creation can be advanced within a few nanoseconds by the use of polarizable DPD water, as demonstrated by Vaiwala et al. 29 .

Stress distribution Using the pressure tensor components Pµµ (µ = x, y, z), the distribution of stress is computed as per Σ(z) =< Pzz − (Pxx + Pyy )/2 >. Here again the stress profiles shown in Figure 11 are comparable. The local stresses decrease with applied field. In simulations with charged plates the profile shows asymmetry, which may be due to buckling of the membrane as observed during the simulation. The membrane tension is deduced by integrating the stress profile, and it shows a decreasing trend with the electric field.

20


Page 21 of 26 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


Conclusions Having summarized the different ways of imposing the electric gradients across the membrane we propose a novel method that sets up the electric field normal to the membrane. Using dissipative particle dynamics simulations, we have demonstrated that the response of the system against applied field by the charged plates is quite similar to the observations made when the membrane is subjected to electrical gradients by the ionic imbalance. The proposed protocol is compatible with the way the electrostatics is being solved by successive over relaxation method in real space. The advantages of the method with the charged plates over the charge imbalance technique are (i) the proposed protocol does not require vacuum slabs, (ii) the method also avoids the necessity of artificially imposing of wall potential, which is inevitable in the case of the ion imbalance method by merit of absence of attractive forces in DPD, (iii) it does not necessitate the need for explicit ions in the system to study electroporation, and (iv) it is easy to implement in DPD simulations. These advantages make the proposed protocol a suitable tool for investigating the membrane electroporation that involve slow formation or closure of electropores.

Appendix The electric field distribution in a column of water (p=1) due to charged plates having a surface charge density ±σ located at z = ±zplate from the center of the column is shown in Figure 12 at various separations of the plates. As confirmed by circles in Figure, the field distribution follows Ez = Γσ1 /p for z < −zplate and z > zplate , while it obeys Ez = −Γσ2 /p for the region between the plates. The σ1 and σ2 values are calculated using σ1 + σ2 = σ and σ1 /σ2 = Lplates /(Lz − Lplates ). Here the size of the column in z-dimension is denoted by Lz , and the distance between two plates is Lplates . It should be noted that the periodicity is applied in three directions. In the case of oil-water slabs, where dielectric inhomogeniety exists, the ratio σ1 /σ2 also involves the dielectric constants of oil and water.

21



40

20 10

40 ± 6Rc ± 7Rc ± 8Rc ± 9Rc ± 10Rc

30 20

± 11Rc

0 -10

10

-10 -20

-30

-30 0

2

4

6

-40 -14 -12 -10 -8 -6 -4 -2

8 10 12 14

± 1Rc

0

-20

-40 -14 -12 -10 -8 -6 -4 -2

± 6Rc ± 5Rc ± 4Rc ± 3Rc ± 2Rc

Ez [kBT/eRc]

30

Ez [kBT/eRc]

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 26

z [Rc]

0

2

4

6

8 10 12 14

z [Rc]

Figure 12: Electric field distribution set up by charged plates with σ = ±2.5e/Rc2 at various plate separations in a column of water (Γ = 13.87). The solid lines indicate the simulation results, while the circles refer to Ez = −Γσ2 /p between the plates and Ez = Γσ1 /p in outer regions. The legends show the locations of the charged plates.

Acknowledgement The authors acknowledge Prof. Ateeque Malani, Department of Science and Technology India, and Indian Institute of Technology Bombay 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. Proceedings of the National Academy of Sciences 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. Proceedings of the National Academy of Sciences 2011, 108, 10443–10447.

22


Page 23 of 26 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


(3) Heller, R.; Gilbert, R.; Jaroszeski, M. J. Clinical applications of electrochemotherapy. Advanced Drug Delivery Reviews 1999, 35, 119–129. (4) Prausnitz, M. R.; Bose, V. G.; Langer, R.; Weaver, J. C. Electroporation of mammalian skin: a mechanism to enhance transdermal drug delivery. Proceedings of the National Academy of Sciences 1993, 90, 10504–10508. (5) 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. Journal of Electroanalytical Chemistry and Interfacial Electrochemistry 1979, 104, 37–52. (6) Benz, R.; Beckers, F.; Zimmermann, U. Reversible electrical breakdown of lipid bilayer membranes: a charge-pulse relaxation study. Journal of Membrane Biology 1979, 48, 181–204. (7) Weaver, J. C.; Chizmadzhev, Y. A. Theory of electroporation: A review. Bioelectrochemistry and Bioenergetics 1996, 41, 135–160. (8) Weaver, J. C. Electroporation of biological membranes from multicellular to nano scales. IEEE Transactions on Dielectrics and Electrical Insulation 2003, 10, 754–768. (9) Chen, C.; Smye, S.; Robinson, M.; Evans, J. Membrane electroporation theories: a review. Medical and Biological Engineering and Computing 2006, 44, 5–14. (10) Tieleman, D. P.; Leontiadou, H.; Mark, A. E.; Marrink, S.-J. Simulation of pore formation in lipid bilayers by mechanical stress and electric fields. Journal of the American Chemical Society 2003, 125, 6382–6383. (11) Tieleman, D. P. The molecular basis of electroporation. BMC Biochemistry 2004, 5, 10.

23



(12) Kramar, P.; Miklavčič, D.; Lebar, A. M. Determination of the lipid bilayer breakdown voltage by means of linear rising signal. Bioelectrochemistry 2007, 70, 23–27. (13) Koronkiewicz, S.; Kalinowski, S. Influence of cholesterol on electroporation of bilayer lipid membranes: chronopotentiometric studies. Biochimica et Biophysica Acta -Biomembranes 2004, 1661, 196–203. (14) Kramar, P.; Delemotte, L.; Lebar, A. M.; Kotulska, M.; Tarek, M.; Miklavčič, D. Molecular-level characterization of lipid membrane electroporation using linearly rising current. Journal of Membrane Biology 2012, 245, 651–659. (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. Biophysical Journal 2008, 95, 1837–1850. (16) Sachs, J. N.; Crozier, P. S.; Woolf, T. B. Atomistic simulations of biologically realistic transmembrane potential gradients. Journal of Chemical Physics 2004, 121, 10847– 10851. (17) Delemotte, L.; Dehez, F.; Treptow, W.; Tarek, M. Modeling membranes under a transmembrane potential. Journal of Physical Chemistry B 2008, 112, 5547–5550. (18) Delemotte, L.; Tarek, M. Molecular dynamics simulations of lipid membrane electroporation. Journal of Membrane Biology 2012, 245, 531–543. (19) Kutzner, C.; Grubmüller, H.; De Groot, B. L.; Zachariae, U. Computational electrophysiology: the molecular dynamics of ion channel permeation and selectivity in atomistic detail. Biophysical Journal 2011, 101, 809–817. (20) 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. 24


Page 24 of 26

Page 25 of 26 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


(21) Groot, R. D. Electrostatic interactions in dissipative particle dynamics - simulation of polyelectrolytes and anionic surfactants. Journal of Chemical Physics 2003, 118, 11265–11277. (22) Español, P.; Warren, P. Statistical mechanics of dissipative particle dynamics. Europhysics Letters 1995, 30, 191. (23) González-Melchor, M.; Mayoral, E.; Velázquez, M. E.; Alejandre, J. Electrostatic interactions in dissipative particle dynamics using Ewald sums. Journal of Chemical Physics 2006, 125, 224107. (24) Vaiwala, R.; Jadhav, S.; Thaokar, R. Electrostatic interactions in dissipative particle dynamics—Ewald-like formalism, error analysis, and pressure computation. Journal of Chemical Physics 2017, 146, 124904. (25) Beckers, J.; Lowe, C.; De Leeuw, S. An iterative PPPM method for simulating Coulombic systems on distributed memory parallel computers. Molecular Simulation 1998, 20, 369–383. (26) Groot, R. D.; Rabone, K. Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophysical Journal 2001, 81, 725–736. (27) 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. (28) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. Journal of Molecular Graphics 1996, 14, 33–38. (29) Vaiwala, R.; Jadhav, S.; Thaokar, R. Four-to-one coarse-grained polarisable water model for dissipative particle dynamics. Molecular Simulation 2018, 44, 540–550.

25



Table of Contents (TOC) graphic

26


Page 26 of 26

Electroporation Using Dissipative Particle Dynamics with a Novel

Recommend Documents