Experimental and theoretical Brownian Dynamics analysis of ion transport during cellular electroporation of E. coli bacteria

Escherichia coli bacteria is a rod-shaped organism composed of a complex double membrane structure. Knowledge of electric field driven ion transport through both membranes and the evolution of their induced permeabilization has important applications in biomedical engineering, delivery of genes and antibacterial agents. However, few studies have been conducted on Gram negative bacteria in this regard considering the contribution of all ion types. To address this gap in knowledge, we have developed a deterministic and stochastic Brownian dynamics model to simulate in 3-D space the motion of ions through pores formed in the plasma membranes of Escherichia coli cells during electroporation. The diffusion coefficient, mobility and translation time of Ca 2+ , Mg 2+ , Na + , K + and Cl - ions within the pore region are estimated from the numerical model. Calculations of pore’s conductance have be en validated with experiments conducted at Gustave Roussy. From the simulations, it was found that the main driving force of ionic uptake during the pulse is the one due to the externally applied electric field. The results from this work provide a better understanding of ion transport during electroporation, aiding in the design of electrical pulses for maximizing ion throughput, primarily for application in cancer treatment.


INTRODUCTION
Electroporation consists in the application of an electric field to biological cells that creates aqueous pathways or pores in the membrane with an estimated minimum radius of 1 nm.The applied electric field causes an induced transmembrane voltage superimposed on the resting membrane potential present under physiological conditions, reaching voltages in the 0.2 -1 V range 68 26 .Electroporation is a subset of electropermeabilization, which ascribes the electrically induced increase in the membrane permeability to a broader range of physical and chemical changes to the membrane lipid bilayer and modulation of its proteins' function 30 .
Most biomedical applications of electroporation aim to transport extracellular ions or molecules into cells, such as calcium electroporation, wherein supraphysiological doses of calcium are internalized causing cell death 27 , irreversible electroporation (IRE) 44 , which is a non-thermal tissue ablation therapy, electrochemotherapy (ECT), consisting of established anticancer drugs with added charge groups to prevent spontaneous entry 21 , and non-viral electrogene transfer, where DNA is first introduced into a tissue, possibly by needle injection, and then the tissue's cells are electroporated within the desired treatment volume 54 62 .Thus, improving the fundamental understanding of ionic transport across the permeabilized membrane, such as diffusion, mobility, and translation times remains a priority given the importance of these phenomena in cell biology and medical applications.
Pulse duration for conventional electroporation is on the order of or longer than the membrane charging time, which also has some impact on the transmembrane potential.Higher-field (> 1 kV/cm), longer pulses ( > > 1 ms) are involved in IRE, necrosis, and electrocution injury 69 .Lower-field (< 1kV/cm), shorter pulses (between 100 μs and 1 ms) are used in ECT and gene delivery 69 .High-field short duration pulses have a negligible thermal effect on the cells and tissues and can preferentially target intracellular components.However, if the field is too high in nanosecond pulses ( > 20 kV/cm) apoptosis occurs 69 .Consequently, it is possible to tailor the operating pulse width to achieve preferential energy deposition into organelles based on their geometry, size, and location.
Manipulation of cells via electroporation can be done in vivo, ex vivo and in vitro.The present work is focused on modeling ion transport during electroporation using in vitro data of Gramnegative Escherichia coli (E.coli) bacteria for validation.The DH5a strain of E. coli bacteria is used in our research because it is a non-pathogenic facultative anaerobe which is easy and inexpensive to grow in large quantities in simple culture media.It is a well-characterized bacteria that is easy to obtain and isolate from humans and animals.It provides insights into the physiology of other pathogenic bacteria of the same family, which can aid in the development of treatments and pharmaceutical products 30 .Due to the difference in anatomy between bacterial and mammalian cells, in particular the double membrane structure with a periplasm region in E coli, vastly higher electric fields need to be applied for electroporation, on the order of tens of kV/cm.While the membrane structure is more complex in E. coli than in mammalian cells, E. coli is a simpler organism overall since it is a prokaryote which lacks a nucleus and other organelles 30 .
Molecular Dynamics (MD) is a microscopic simulation scheme which is highly representative of the molecular-level structure of the membrane.It Is well suited to simulate pore formation considering the details of the chemical components of a small region of the lipid bilayer, as performed by Piggot et.al. 50 for S. Aureus and the outer membrane of E. coli.However, long time scales are required to corroborate experimental conductivity measurements, which cannot be obtained with microscopic MD simulations due to their long computation time 11 .Besides, the atomistic structure of pores is not fully known experimentally, and inaccurate models of water molecules may be used in the calculations 51 39 .On the other hand, macroscopic analyses that treat the cellular matter as a continuum fluid with voltage barriers are more suitable for longer time scales but are less representative of the biological structure of the membrane's pore 56 55 .
Therefore, simulations of ionic transport of electroporated E. coli bacteria were carried out in the present work using a semi-microscopic Brownian Dynamics (BD) scheme, which combines the advantages of both microscopic and macroscopic models, but also inherits some of their disadvantages, namely time step dependence, exclusion of chemical reactions, loss of information due to simplified dynamics, and accumulation of errors over time 5 .Besides, BD models are Markov processes that are limited to steady state dilute systems 9 .However, the advantages outweigh the disadvantages considering the characteristics of the system being studied.For instance, the BD approach is computationally more efficient and faster than MD simulations, since it is a coarse-grained model that groups collectively multiple membrane and water molecules 51 39 .This allows for the study at longer time scales of a larger system consisting of both the inner and outer membranes of E. coli, while keeping a higher resolution than macroscopic continuum models 51 39 .The BD scheme incorporates hydrodynamic interactions and thermal fluctuations directly through stochastic forces, allowing for longer time steps ΔT and requiring fewer samples to capture the essential statistical behavior of the system 51 39 .The intermediate resolution level of BD simulations facilitates multiscale bridging between atomistic and macroscopic phenomena 39 , as has been done in the present work for conductivity validation.Thus, MD simulations are suitable to study biomolecular dynamics in solutions, diffusion transport processes and ionic conductance of channels and pores 51 .

Experimental Setup
Experimental studies were conducted at the Vectorology Department of the Gustave Roussy Institute using DH5a E. Coli bacteria in order to obtain measurements of conductivity and permeability of cell membranes after electroporation.In this regard, a Bio-Rad Gene Pulser Xcell™ electroporator set to voltage V0 = 1800 V, resistance R = 150 Ω and capacitance C = 50 F was used to apply a single exponential decay electrical pulse with time constant  =  = 7.5  on a 1mL capacity Bio-Rad plastic cuvette with gap width w = 1 mm, creating over time t an electric field  ⃗ =  0  −/ / = 18 kV/cm at t = 0. Figure 1 a) shows the experimental arrangement with a simplified circuitry of the pulse generator.The Bio-Rad cuvette incorporates aluminum electrodes on two opposite internal walls which get in contact with the electrical circuit.Switch 1 is first closed to charge the capacitor and then opened to isolate it from the source.The pulse is discharged on the cuvette suspension by closing switch 2. Figure 1 b) represents the waveform of the applied pulse.The amount of current released by the capacitor is related to the time at which the energy stored in it is allowed to dissipate, which depends on the resistance across the electrode assembly and sample.Therefore, the higher the resistance, the longer the pulse.Since the E. coli cell charging time is usually 100 ns, long pulses must be used if the membrane potential is to exceed 0.5-1V for permeabilization.Conductivity assays were conducted in three different physiological states: fresh bacteria in growth phase, overnight bacteria coming from a saturated culture, and electrocompetent bacteria.E. coli DH5a was grown in Luria broth (LB) media.Cultures were kept at 37 ºC, shaking at 200 rpm.For the electroporation selection, LB agar plates were prepared using 50 μg/mL kanamycin.Electrocompetent E. coli were prepared by growing the bacteria in 500 mL of LB media at 200 rpm until reaching an optical density of 0.4 -0.6.The whole volume was then centrifuged at 3000 rpm at 4 ºC for 15 min, discarding the supernatant and resuspending the pellet in 500 mL 10% glycerol in DI water.This process was repeated thrice by resuspending the cells into 200 mL, 20 mL and 2 mL, aliquoting the final volume into 100 uL fractions.
Before conducting each conductivity assay, the FiveEasy FiveGo FE30 conductometer probe was washed with deionized (DI) water, then 70% ethanol in DI water, then washed thrice with DI water.This step was necessary to reduce conductivity to background levels (1.5 μS/cm or less).
A special column for measuring a minimum amount of sample was crafted by cutting a disposable 25 mL plastic pipette of enough diameter to fit the conductometer's probe, which was then attached to a flat acrylic base.This column was also washed as described.In order to perform triplicate assays, 33 uL of each sample, or 33 uL of sample plus 1 uL of 1 μM YO-PRO-1 were taken into a final volume of 1500 DI water.
Flow cytometry assays were performed on a BD Accuri™ C6 Flow Cytometer in order to determine the cell count and differences in cell membrane permeabilization in the cited states.
A 488 nm laser was used for YO-PRO-1 excitation.The capillary was cleaned by flowing DI water before starting.The device was blanked with 233 μL DI water, 10% glycerol in DI water or 10% LB in DI water, depending on the sample.233 μL of sample dilution were then analyzed on the Fast delivery rate at 66 μL/min, by diluting 33 uL of each sample with 200 uL DI water or the corresponding diluent.For the samples tested with the nucleic acid stain, 1 μL of 1 mM YO-PRO-1 was added to 33 μL of sample plus 199 μL of DI water or the corresponding diluent.
Prior to the electroporation experiment the bacterial samples were thawed on ice.Additionally, the 1 mm electroporation cuvettes were kept on ice.33 μL of cells were transferred to the electroporation cuvette, removing any bubbles by gently tapping, drying excess external humidity and placing the cuvette into the electroporator chamber.For the samples treated with YO-PRO-1, 1.5 mL Eppendorf tubes were kept on ice before starting.Samples were then pipetted into the Eppendorf tubes by duplicate, adding 1 μL of 1 mM YO-PRO-1 to the treatment and 1 μL H2O to the mock, pipetting the volume to mix.

Basic Equation and Algorithm
From a semi-microscopic point of view, the movement of an i th ion with mass mi suspended in water with velocity vi during cellular electroporation follows Newton's equations of motion and is best described for long simulation periods by a modified Langevin eqn. 37: where  , represents an stochastic force arising from the interaction of ions with water molecules,  , is a friction force of ions with water,  , is a force due to an externally applied electric field, the force  , arises from the ion's interaction with the membrane and pores,  , is an interionic Coulomb force, and  , is the short range repulsion force from the overlap of the ions' electron clouds.
The velocity-Verlet algorithm 63  In order to speed up computation time of the forces, besides using Brownian dynamics as opposed to Molecular dynamics, the workload for ions is distributed evenly amongst various core processors using parallel processing.Since the Coulomb, repulsive and image charge forces have a computational time complexity of O(n 2 ) and they are dependent on other ions' positions, the message passing interface (MPI) protocol is used to pass the ions' positions among processors.In particular, the mpi_allgatherv function is used, which passes messages using an efficient tree algorithm.The speed up curve reaches a plateau at 50 cores, after which computation time for each time step decreases due saturation of the bus with messages.For computing the diffusion coefficient from the Green-Kubo relation (Diffusion coefficient section), ion's velocities are passed as well to the root process using mpi broadcast.For additional improvement in computational time, the simulations are conducted considering finite subregions of the cytoplasm, inner and outer membranes, periplasm and extracellular environment.The program has a modular design and includes a Makefile for compiling.

Details on the Calculation of Various Forces
Friction Force (inelastic collision with the surrounding water) The total friction force acting on an ion from water molecules it hits at a particular instance is opposite to the direction of motion.It is represented by where 1/γi is the relaxation time constant 24 28 .Fr is considered directly proportional to the ion's velocity, as non-linear restoring forces only occur when the ions move extremely fast.The coefficient of friction miγi is given by Stokes law: where RS is the Stokes radius of the hydrated molecule and nvis is the viscosity of liquid water, which as a function of temperature T can be expressed empirically by 33   = 10 where A = 2.414x10 -5 N.s/m 2 , B = 247.8K , and C = 140 K .Considering the following relation from 66 : where kB is Boltzmann constant and Di is the bulk diffusion coefficient of an i th ion, the Stokes radius is obtained by substituting Eqn. ( 5) in (3) to obtain At T = 298 K, nvis = 8.9x10 -4 N.s/m 2 , and Di is given in Table 3.

Random Force (elastic collision with water molecules)
The interaction of ions with water molecules is represented with isotropic random collision forces which are Markovian in nature with Gaussian distribution.The Gaussian random collision force FR has a sample space SFR = (−∞, ∞) 49 and its probability density function (pdf) is given by where the mean me is zero and the standard deviation σ is expressed as 24 66  = √ Since the pseudo-random floating point numbers xi generated by a computer algorithm have a sample space Sxi = (0,1), for a random variable X with uniform pdf, f(FR) is normalized.Hence, for every xi there is a corresponding FRi such that Peforming the change of variable  =   /√2 and multiplying both sides of Eqn. ( 9) by 2/√ : The numerator is the Laplace error function erf(t) for  ≥ 0. The integral in the denominator is √/2 .Hence erf(∞) = 1 and, considering the aproximation in 25 , Eqn. (10) becomes where the coefficients a1=0.3480242,a2=0.0958798,a3=0.7478556, and  = 1/(1 + ) with p=0.47047.
To save computational time, a look-up table implemented as a 2D array is used to store 400 error function values between 0 and 1 in one column, which are matched to the generated pseudo-random numbers, and the corresponding Gaussian random forces in the other column.

Coulomb Force (other ions in the reservoir)
The Coulomb force considered for the n-ion system is given by where   is the net force on an i th ion due to all other ions within the system,   ⃗ ⃗ and   ⃗⃗⃗ are the position vectors of the ions and q is the ion electric charge.For simulation purposes, the vector   is broken down into its Cartesian projections.For the x axis it is given by , and  , have similar expressions to Eqn. (13), except that  , −  , in the numerator is replaced by  , −  , or  , −  , , respectively.

Image Charge Force (charges induced at the membrane)
Considering the effects arising from the difference in dielectric properties between the membranes and surrounding water, due to dielectric polarization when an ion gets close to a non-conducting membrane wall, it induces bound charges within the membrane's surface, which in turn induce a force  , on the ion itself.For simple geometries, this phenomenon can best be treated by the method of image charges 59 .
The membrane surfaces outside and inside the cell in the area surrounding the pore are approximated with planar boundaries, so that  , is taken into account by a Coulomb force (Eqn.12) due to an image charge q' 59 : where εm and εw are the permittivities of the membranes and the water solvent, respectively.In terms of  0 ,   =    0 and   =    0 where   and   are the relative permittivities of the membranes and water.
The charges induced close to the surface of the pore formed by electroporation are approximated considering the Green function of the interior of a cylinder of radius ρ0 and length 2z0 64 65 .
The z component on a charge q can be found by evaluating Green's function G at the cylindrical coordinates ′ = , ∅′ = ∅, ′ =  and taking the gradient following Smythe 60 : The radial force Fρ is calculated considering the charge distribution σ on the cylinder surface S, which can be found from an alternate form of Green's function 70 : where the gradient is taken with respect to , ∅, .Instead of using Coulomb's law to calculate the force on the ion due to σ, it is easier to calculate the total force on the conducting membrane pore, which is the negative of the force on the ion (Newton's third law).
Thus, the radial force on the ion is obtained from The total  , acting on an i th ion is the vector addition of forces with components Fz and Fρ arising from its own induced image charge as well as all other image charges induced by other ions in the pore (Details of derivations provided in the Appendix).

Force from Externally Applied Electric Field
The Lorentz force acting on an i th ion is given by where  ⃗ and  ⃗ are the electric and magnetic fields acting on the system. ⃗ has two components: one due to the externally applied pulse throughout the cuvette, as explained in the experimental setup section, and another one at the membrane related to the resting potential arising from differences in intracellular and extracellular ionic concentrations.The resting potentials Vrest at the inner or outer membranes of E. coli (Table 2) enhanced by the externally induced voltage Vm are related to the  ⃗ field as follows: where points a and b are at the edges of the pores as shown in Figure 2, and  is an infinitesimally small section of its total length.In the simulated system,  ⃗ is assumed to remain constant in space over the entire region of integration (i.e. the pore's area) for each time step, while exponentially decaying over consecutive time steps due to the nature of the applied pulse.The magnetic field  ⃗ is induced by the change in the electric field   due to the discharge of the capacitor, with magnitude where µ0 is the is the permeability of free space,   = √ 2 +  2 is the position of the ions in a plane parallel to the electrodes with area Ae = 200 mm 2 , and  =  0  −/ / is the current flowing from the capacitor through the wires connected to the electrodes.
The induced transmembrane potential in relation to the shape and size of the long prolate spheroidal E. coli cell and the externally applied electric field  ⃗ is given by 19   =   (), where the form factor f = 1.12 if the long axis of E. coli (Laxis in Table 2) and the field are parallel (angle θ = 0), and f = 1.8 if they are perpendicular 19 , the physiological factor g = 1 assuming the membrane is a pure dielectric, and rl is half the axis of the cell in the direction of the field, The induced potential is highest when the rod is parallel to the field since the shorter axes are three times smaller than the long axis (1:1:3 axial ratio).To calculate Vm at the inner membrane, the periplasm and outer membrane widths are subtracted from the cell's axis.
Thus, when ions are in the pores, the electric field in Eqn. ( 18) is calculated by substituting Vm from Eqn. ( 21) and Vrest from Table 2 into Eqn.(19), divided by the membranes' thickness.This approach is in accordance with standard RC cell theory for small electric fields 15 .However, experiments showed a flattening of Vm to approximately 1V in the polar regions for larger  ⃗ fields 15 .Thus, a cutoff has been implemented for voltages above 1V.The  ⃗ field is different at the inner and outer membranes, since the membranes' thickness, resting and induced potentials are different.In regions surrounding the membrane,  ⃗ in Eqn. ( 18) is the externally applied electric field.The induced magnetic field  ⃗ is calculated using Eqn.( 20) throughout the simulation regions, since only the external  ⃗ field component changes with time, assuming Vrest remains constant.

Short Range Repulsion Force from the Overlap of the Ions' Electron Clouds
For ions having different charges with opposite polarities, a short range repulsion force arising from the overlap of the electron clouds is considered.In this regard, an interatomic potential U, based on the Thomas-Fermi-Dirac approximation can be expressed in the so-called Born-Mayer form within an interval of internuclear separation R 1 10 : where Rl and Rii are the lower and upper limits of the interval, typically ~1.5a0 (a0=0.52917Å) and ~3.5a0, respectively, and parameters As and b are constants for a given pair of same-type atoms.For pairs of unlike atoms 1 and 2, the following empirical "combination rule" is used: The interatomic force is  , = −∇ 12 .

Pore Creation and Evolution
Under a typical physical scenario, pores are formed in the cell membrane during electroporation, reach a maximum radius, and subsequently lead to plasma membrane recovery through shrinkage, or membrane disintegration.The rate of pore creation and destruction is governed by the induced potential Vm (Eqn.21) according to the following first order ordinary differential equation 15 46 : where N is the pore density, and α, β, Vep, and N0 are constants provided in Table 2.The total number of pores K at a given time is the surface integral of the pore density:  = ∬   .Since pores concentrate on both ends of the E. coli rod after it rotates parallel to the field (angle θ goes to zero in Vm), their surface areas combined can be approximated as A=2πrw 2 , where rw is half the shortest axis (Laxis,short in Table 2), i.e.K=2πrw 2 N. Eqn. ( 24) is solved with the initial condition  = 0 (no pores).
Due to the differences in transmembrane voltage, thickness, position and composition of both membranes, the sizes of the inner and outer pores evolve differently over time to minimize the energy of the lipid bilayers.The change in the radius ri of an i th pore is determined by the discretized equation The first term accounts for the electric force induced by the local transmembrane potential Vm and applies to pores with arbitrary size and a cylindrical inner surface 48 ; the second term accounts for the steric repulsion of lipid heads 35 46 ; the third, for the line tension acting on the pore perimeter 35 47 46 ; and the fourth, for the surface tension of the cell membrane, where the effective tension σeff is given by 47 46 where σ0 is the tension of the membrane without pores, σ' is the energy per area of the hydrocarbon-water interface, and is the combined area of pores at a given time 47 35 , which are initially created with minimum radius   =  * 35 .The average pore radius is 24), ( 25) and ( 26) are defined in Table 2, assumed to be the same for both inner and outer membranes except Vm, which has the greatest effect.

Initial Conditions Given to the Simulated Ions
The Maxwell-Boltzmann probability density function () describing the distribution of molecular speeds in terms of the translational kinetic energy εi of a an i th ion is given by Since the pseudo-random floating-point numbers xi has a sample space Sxi = (0,1), for a random variable x with uniform pdf, f(FR) is normalized.Hence, for every xi there is a corresponding εi such that The integrals in the above equation are solved directly giving where   ′ = 1 −   , it follows Substituting the kinetic energy Eqn.  =     2 /2 in the above, we get the initial speed of an i th ion: noting that ln(x') is a negative number.The ion's initial velocity vector with magnitude equal to its initial speed is in a direction specified by the spherical angles  and , which are randomly calculated based on a uniform pdf.The velocity spherical coordinates are transformed to Cartesian when incorporated in the model.
The initial positions of ions have been randomly assigned with uniform pdf within the cytoplasm subregion under simulation.

Geometry and Boundary Conditions of the Simulated Regions
The E. coli rod is subjected to a torque and rotates because of the induced dipole, reorienting itself with the longest axis parallel to the applied field direction.Only at the end of orientation permeabilization can be efficient 19 .In the simulation, pore formation is considered in the inner and outer membranes (Figure 2).Both pores are aligned with each other in a membrane region at the hemispherical end of the rod along its larger axis in the z direction parallel to the field.In this position, the highest transmembrane potential is induced (Eqn.21), and thus there is a higher probability of pore formation.
For practical purposes, the pores are modeled with a cylindrical shape with sharp boundaries and evolve over time according to Eqns. ( 24) to (26).The inner and outer sides of the cell are modeled with finite cubic volumes with imaginary or soft boundaries (See Figure 2).
When an ion leaves the system, another ion is symmetrically inserted back from the opposite boundary, thereby keeping the number of particles in the system constant in order to approximate a canonical ensemble in thermal equilibrium with a fixed temperature and volume.For the sake of simplicity and computational efficiency, this approach has been chosen over stochastic Monte Carlo techniques 31 , which are more complex but do not necessarily improve the accuracy of the predictions 14 .The initial specific ion concentrations in the intracellular and extracellular subregions are provided in Table 3. Ions leaving the extracellular subregion across the boundary opposite the outer membrane are reinserted in the cytoplasm.
In order to avoid short range interactions when applying periodic boundary conditions, the shortest side of the cytosol, periplasm and extracellular boxes have been chosen to be at least twice as big as the cut-off radius for short range electrostatic forces 16 .The box sizes are sufficiently small to keep computational efficiency while having a proportional ionic concentration, with 10 nm sides for 1 and 4 nm pore radius, and 25 nm to contain the pores with 20 nm radius.Bigger simulation boxes have a negligible effect on both the thermodynamics and kinematics of ions at the studied regime 20 .On the other hand, water diffusion in MD simulations is affected by the box size 20 .In the case of the present BD work, however, water is treated as a continuum and its diffusion is hence unaffected.
Interactions between the ions and the pore or the cell membrane surrounding it are elastic in nature, i.e., the total momentum and kinetic energy are conserved, and the collisions are handled with specular scattering at hard boundaries, with the angle of incidence equal to the angle of reflection/deflection.

Calculation of Transport Parameters Diffusion Coefficient
The diffusion coefficient D is calculated by integrating the velocity autocorrelation function C(t) for 3 dimensions (Green-Kubo relation) 23 36 : C(t) is given by where   ( 0 ) is the ion's initial velocity and the notation <> signifies the average.When evaluated over N ions, the above Eqn. is expressed as follows: The integration of C(t) is approximated with the trapezoidal rule.

Mobility
The diffusion coefficient D of ions is related to their mobility by Einstein's relation 17 :

Translation (Membrane crossing) Time
Typically, only a few ions flow through a pore simultaneously because of the small pore size and repulsive forces among ions.Thus, simulations were run for only one ion crossing the pore at a time, thereby saving computational time (approximately seven hours for each ion crossing on one 2.1 Ghz core, different crossings were sent to different cores with MPI in parallel).In the case of bigger pores with a 20 nm radius or greater, there is a possibility of multiple ions being present simultaneously.Thus, a multi-ion treatment was utilized for calculating the membrane crossing time in bigger pores, as well as for determining the conductivity validated with experimental data.

Current and Conductivity
The conductivity for an n concentration of ions is given by The current in the pore can be calculated by The current density J is the current per cross sectional area A: The first term on the right-hand side is the drift current (Ohm's law), the second term is the diffusion current (Fick's 1 st law) Also from Eqns ( 36) to (38), the drift current is

Diffusion Calculation of Diffusion Coefficient in Bulk for Model Validation
The bulk diffusion coefficient is known from experimental observations, and hence it is possible to compare with a computer simulation to verify the accuracy of the calculations.Figure 3 shows the velocity autocorrelation function (Eqn.( 33)) as a function of time for 1000 ions of each type suspended in a bulk system.The calculated bulk diffusion coefficients at 298 K (provided in Table 4) are the areas under the curves (Eqn.( 32)) which closely match the values reported by other authors (provided in Table 3), with an average root mean square error equal to 2.8132x10 -10 m 2 /s for ten simulation runs of each ion type.The difference between the values does not have a significant effect on the results during the pulse since the main driving force is due to the applied electric field, as explained later in the discussion section.However, the small difference in the coefficients might have a minor effect on the calculations after the pulse when transport is primarily due to diffusion.The dependence of the diffusion coefficient on temperature has also been studied.Simulations were run for ions within the pore under experimental conditions (T = 298 K, E = 18 kV/cm, ravg = 18 nm from Eqn. ( 25) after 1 ms) and at 278 K and 318 K.These temperatures were chosen so they are not far apart from each other, are within the temperature range for E. coli growth (277 K to 322 K -ideally 310 K) and its survival limits, typically greater than T = 193 K (below freezing point) and less than water's boiling point 3 .This temperature range can be easily considered in future experiments.Care is taken to avoid phase change, since at temperatures below 273 K water is frozen so friction increases sharply and Brownian Motion is no longer observed nor can be simulated with the model.
It has been observed that the diffusion coefficient increases as the temperature gets higher, as shown in Figure 4. Error bars represent one standard deviation of uncertainty of ten simulation runs.Smaller ions have less friction with water molecules, thus having higher diffusion with steeper dependence on temperature, which is related to the decrease of water viscosity with heating proportional to the ion radius (Eqns.2-6) The findings agree with the kinetic theory of gases.

FIGURE 4.
Dependence of the diffusion coefficient on temperature.Values are shown with standard deviation error bars for Ca 2+ , Mg 2+ , K + , Na + and Cl -ions within the pore.

Diffusion in Pore vs. E-field
Simulations were run for ions in the pore under the influence of an external force due to an electric field.The results agree with the theory which states that if a field is applied, the particles' diffusion will remain the same because no energy-dependent frictional force or nonuniform scattering process has been incorporated (other than its relation to temperature).In the case of an energy dependent scattering, e.g., the monotonic increase with energy observed in semiconductor transport, the diffusion coefficient would have exhibited a corresponding decrease.Figure 5 depicts the diffusion coefficient values in the pore with one standard deviation error bars of ten simulation runs at the experimental room temperature condition (T = 298 K) for different magnitudes of an externally applied field.The pore radius is related to the applied electric field by Eqns.( 21) and (25).The slight variation of the diffusion coefficient seen in the figure is statistical in origin arising from the Brownian motion of ions due to the surrounding water molecules, accounted for by the random Gaussian force FR.

FIGURE 5.
Diffusion coefficient values with standard deviation error bars for Ca 2+ , Mg 2+ , K + , Na + and Cl -ions in the pore for various magnitudes of an externally applied electric field at temperature T = 298 K.

Diffusion in Pore vs. Pore Size
The diffusion coefficients of ions in pores with different radius sizes are shown in Figure 6.The effect of small and large pores are explored, with radius values differing by similar factors (x4 or x5).Other conditions correspond to the experiment (T = 298 K, E = 18 kV/cm).It can be concluded from the results in the figure that the diffusion of ions does not vary much according to the pore's radius.This might be due to the fact that random forces do not move the ion far away from the straight trajectory it would otherwise have if no random forces were present.Hence the ion hardly ever touches the boundary, if at all, when initially placed at the center of the pore.Some model results in the literature, however, point at an increase in ionic uptake with increasing pore size, due to a crossover from cation to anion specific permeation 38 , or a monotonically decreasing energy barrier with increasing pore radius 32 .In particular, MD simulations found a Cl -flux increase from 0.02 ns -1 nm -2 for a medium sized pore (Area = 4 nm 2 ) to 0.1 ns -1 nm -2 for a larger pore (Area = 9 nm 2 ), while Na + showed no significant increase in the permeation rate 38 .

FIGURE 6.
Diffusion coefficient values with standard deviation error bars for Ca 2+ , Mg 2+ , K + , Na + and Cl - ions in the pore for various radius sizes at T = 298 K.

Mobility and Drift Velocity
Table 4 below lists calculated mean mobility values of ions moving through the pore using Einstein's Eqn.(35) for the corresponding mean diffusion coefficient of ten simulation runs calculated with the velocity autocorrelation function (Eqns.32-34) at temperature T = 298 K, applied electric field = 18 kV/cm, pore radius r = 1 nm.The diffusion coefficient values in Table 4 correspond to those shown in Figure 4 at room temperature, and Figures 5 and 6 irrespective of the applied field and pore radius.The drift velocity is the average velocity in the z direction parallel to the pore's axis and is obtained directly from model output.

Table 4
Calculated mean mobility values for the corresponding ionic diffusion.

Translation (Membrane crossing) time Histograms of Transport Time of Ions vs. Pore Size
In order to evaluate ionic throughput dependence on pore size, simulations were run 100 times for each ion type to have sufficient statistics with a low fractional sampling error of 0.142 8 .Figures 7, 8 and 9 show histograms of Ca 2+ , Mg 2+ , K + , Na + and Cl -ions' transport time through pores with radius = 1, 5 and 20 nm under experimental conditions which are standard for bacterial elecgtroporation, i.e., a high 18 kV/cm applied electric field at temperature T = 298 K to permeabilize both the inner and outer membranes of E. coli.The Cl -anion flows in the opposite direction of the cations in pores created on the other side of the cell.Table 5 tabulates the mean and standard deviation of the calculated inner membrane ion crossing times through pores with different radii (shown in Figures 7, 8 and 9).Applying a T-test statistical analysis to these results, it is determined that there is no statistically significant difference in transport time.This might be due to the fact that the random forces arising from water molecules and induced charges at the membrane do not move the ions far away from a straight trajectory propelled by the electric field.

Conductivity -Model Validation with Experiments and Comparative Evaluation for Different Pulses
In order to estimate the average number of permeated E Coli cells in the samples, measurements were taken with a flow cytometer in reference to an unpulsed control, indicating 1.57x10 10 electroporated cells in the fresh sample, and 1.47x10 11 in the overnight electrocompetent sample.The percentage of electroporated cells is 87.65 % and 85.87 % for the fresh and overnight samples, respectively.Representative graphs of the flow cytometry data are provided in Figure 10.The computed contributions of each ion type to the increase in the conductivity (Eqn.36) of electrocompetent E. coli are shown in Figure 11 as a function of time during the first nanoseconds of the exponential decay pulse with V0 = 1800 V and  =7.5 ms.The electric pulse starts 1 ms before the results shown, which is the time it takes for all E. coli cells to realign themselves and trigger an efficient permeabilization at their hemispherical ends, with their longer axes parallel to the applied electric field 19 .The conductivities shown in Figure 11 are calculated multiplying the ionic efflux of an inner/outer membrane pore obtained from the BD model times the number of pores per cell considering their average radius, and the total number of cells in the sample from the flow cytometry data.Applying Eqn.(24), 11 inner/outer membrane pores start to be created per cell at time t = 1 ms with an average radius of 35 nm.The discrete jumps in conductivity seen in Figure 11 are related to the nature of the simulation, since ion efflux is stored every 100 time steps.These results are comparable to the temporal growth in ionic concentration depicted in the continuum model work by Q. Hu and R. P Joshi 29 , although they focus on calcium ion transport in mammalian cells subjected to monopolar and bipolar square pulses with higher amplitude.Cations flow in the direction of the applied electric field and the anion in the opposite direction in pores created on the other side of the cell.The corresponding drift current density is the first term in Eqn.(38), i.e., the conductivity multiplied by the field.

FIGURE 11.
Computed electrophoretic conductivity of electrocompetent E. coli as a function of time for each ion type in the nanosecond regime during the pulse.Time t = 0 corresponds to 11 ms after the pulse rising time.
After application of the pulse, ions continue to diffuse through long-lived pores or permeable structural defects until homeostasis is reached, pores are fully resealed, and the ion membrane permeation barrier is restored.The increase in the reservoir's extracellular concentration of ions cions with time is characterized quantitatively by the Nernst-Planck equation, which has the following analytical solution for passive diffusion after the pulse 52 : where , , C and τ1 have the same meanings as in reference 52 .In the above equation, the exponential decay accounts for the resealing of pores with time.While the BD model simulates a canonical ensemble with a constant number of ions in nanoscale subregions, it keeps track of ions exiting the pores which contribute to an overall increase in the macroscopic reservoir's ionic concentration.
Considering Eqn.(36), the passive diffusion conductivity   is related to the ionic concentration as follows: where the ionic charges q are provided in Table 3 and the mobilities  in Table 4.The equation terms correspond to the conductivity change from individual ions, with the decimal numbers being proportions of the total ion concentration determined by the BD model in nanosecond resolution with no applied field.Under the same external conditions, this proportionality trend remains constant over longer periods of time.The number of pores per cell npores is calculated with Eqn.(24), ncells is the number of cells in the sample derived from flow cytometry measurements.Thus, it is possible to scale down to an individual pore conductivity for better comparison with the model.During the pulse, the solution to the Nernst-Planck takes the linear form   () =  52 .
Figure 12 shows the conductivity values obtained using Eqn.(41) curve fitted using linear regression to the average experimental data of electroporated triplicate samples of fresh and overnight electrocompetent DH5a E coli bacteria.The latter sample has a higher conductivity due to its higher permeability and number of cells.Fresh DH5a conductivity measurements were taken every 30 s for 7 mins, starting 45 s after the pulse.Electrocompetent DH5a conductivity measurements were taken every 30 s for 10.5 mins, starting 25 s after the pulse.The average coefficients of variation are 2.3181 % and 4.9292 % for the fresh and overnight triplicate data, respectively.The FiveEasy FiveGo FE30 conductivity probe employed has an experimental uncertainty equal to 0.5 % of the measured values.Considering human reaction time, the time intervals measured with a stopwatch have an absolute experimental uncertainty of 0.2 s.The experimental uncertainty of the time intervals is not depicted with horizontal error bars since it is unnoticeable in the figure's time scale.The model closely fits the measured conductivity data, with root mean square error equal to 2.0839 µS/m for fresh bacteria, and 2.898 µS/m for electrocompetent bacteria.The increase in the conductivity of the suspension over several seconds is in accordance with an experimental report 58 which indicates that longer electric pulses create larger pores which take longer to reseal, e.g., 6-7 mins for a 100 us pulse, 15-20 mins for a 2 ms pulse.In the present study the pulse duration time constant is 7.5 ms.According to Eqn. (24), the pore density increases steadily and reaches a maximum value, after which it decreases slowly over several milliseconds.The maximum pore density is related to the magnitude of the applied electric field; the higher the field, the longer it takes to reach a higher maximum value.The average radius decreases over time due to the exponential decay of the electric field and related induced membrane voltage in the first term in Eqn.(25).
We explore the effect on the conductivity of various electric fields ranging from 12 to 21 kV/cm after 11 ms, when pore density has reached a maximum steady value for all fields considered, thus allowing a better comparison.Considering exponential decay pulses, for E(t=0) = 12 kV/cm, 38 pores are formed per cell after 11 ms with average radius ravg = 14 nm; for E(t=0) = 15 kV/cm 47 pores are formed with ravg = 16 nm; for E(t=0) = 18 kV/cm 54 pores are created with ravg = 18 nm; for E(t=0) = 21 kV/cm 60 pores with ravg = 20 nm.Considering a square wave pulse with a constant field E(t) = 18 kV/cm for 11 ns, 58 pores are created with ravg = 20 nm.
Figure 13 shows the conductivity of electrocompetent E. coli for different applied electric fields for a nanosecond time period starting 11 ms after the pulse rise time (at t = 0).At t = 11 ms, the electrophoretic conductivity can be determined from the Nerst-Planck equation 52 .For better comparison, only the increase in conductivity is shown, i.e., the conductivity is zero at t = 11 ms.The difference in conductivity with respect to the applied field is mainly due to the variation in the number of pores formed per cell and their average radius, with higher conductivity for higher pore number and radius.The difference in the ionic efflux with different E-fields of a single pair of inner/outer pores is less significant in the nanosecond regime.As can be seen in Figure 13, the conductivity is higher with the square wave pulse than the corresponding exponential decay pulse, due to the initial E(0) = 18 kV/cm field being sustained for the entire 11 ns period.Due to cytoplasmic leakage present after the pulse, E coli electroporation is not a fully reversible process.Thus, two successive pulses under crossed polarization are not possible 19 .
Figure 13.Simulation results probing the effect of changes in the pulse electric field amplitude on the temporal increase in the conductivity of electrocompetent E. coli.Time t = 0 corresponds to 11 ms after the pulse rising time.The starting conductivity of zero corresponds to the conductivity of the sample at 11 ms.

DISCUSSION
The Brownian dynamics approach has been chosen over the MD method (microscopic) and the continuum model (macroscopic) because it is computationally fast while at the same time the level of abstraction is low enough so that effects can be reasonably measured.From the BD model calculations, it was found that the main driving force of ionic uptake during the electric pulse is the one due to the externally applied electric field (electrophoretically mediated transfer), while other effects such as random collisions of ions with the medium or induced charges at the membrane only decrease ion throughput slightly.In Eqn.(1), the order of magnitude for the Coulomb and repulsive forces depends on the distance amongst the ions; the larger the distance the smaller the force.On average, the short-range repulsive force is several orders of magnitude smaller than the other forces ( < 10 -17 N), while Coulomb, induced image, random and friction forces are in the order of 10 -10 N. The applied electric field produces a force outside the membrane ~2.8x10 -1 N.However, inside the membrane it becomes stronger and comparable to the Coulomb force.The main driving force is due to the applied electric field since it is oriented in one direction, while other forces point in all directions and cancel out or have their overall magnitude greatly diminished.
When applying the PEF the main contributors to the conductivity (Figure 11) are the K + cations due to their higher intra-cellular concentration, small size, and high diffusivity.Mg 2+ cations also contribute significantly to the conductivity due to their double charge and relatively high concentration.After the pulse, however, Mg 2+ contribute less due to their bigger size and thus lower diffusion.Despite having positive charges, few Ca 2+ and Na + ions are transported during the pulse because of low concentration and diffusivity (or mobility).The negatively charged Cl - anions move slowly in the opposite direction to the field, with strong attraction by the positive ions.
The bulk diffusion coefficient and conductivity calculations were compared to experimental observations, constituting a framework for calibrating the variables and validating the model in order to obtain realistic values for the unknown transport parameters, i.e., the mobility, diffusion, and translation time of ions within the pores.With the external electric field set to correspond to the applied pulse, the BD model calculates the conductivities due to electrophoretic transport of ions.With the electric field set to zero, the model calculates the increase in conductivity due to passive diffusion.These results present a tendency which is integrated into an analytical solution of the Nernst-Planck equation (Eqn.40) for comparison with the experimentally determined extracellular conductivity.Measurements were taken over several seconds after the pulse when ion transport is due to facilitated and simple passive diffusion through long lived pores and protein channels 52 .The model fits reasonably well to the experiments considering its limitations (Figure 12).Further validations could be carried out by comparison of results with other microscopic and macroscopic models.A physical approximation introduced by the BD simulation is the neglect of the hydration (or solvation) spheres.Typically, the presence of polar water molecules leads to a shielded arrangement around the ions, which can best be captured by full MD simulations.This could become an issue for transport of ions across narrow pores, as the radius of the hydration shell approaches the channel opening.In order to incorporate in the BD model the average energetic contribution of the solute-solvent hydrogen bonds in the first solvation shell, implicit solvation models may be utilized 42 43 .Certain calculations, such as those involved in Green's function for the image charge forces, require summations to infinity.For practical purposes, a cut-off length has been introduced (.e., m,n = 4).Since the numbers used on a computer have a finite amount of digits, only a small subset of the rational numbers can be represented exactly with the IEEE binary floating-point arithmetic standard.Both rounding up and down introduce several errors in the calculations, such as the cancelation of significant digits due to subtraction of nearly equal number, division by a number with small magnitude, or multiplication by a large magnitude.Whenever possible, polynomials and equations have been written in forms that minimize the number of these calculations in order to make the occurrence of round-off errors less probable.
In order to better simulate pore formation and closure, more membrane details could be included in the model.Although the exact structure of pores is unknown experimentally, adding the effects of hydrophobic and hydrophilic molecules that form the membrane's lipid bilayer might help get more refined simulation results.To this aim, a hybrid MD-BD simulation scheme may be used which includes dipole shifts and configurational changes.Poisson solvers may be utilized to obtain the potential energy within the membrane and pore for a given configuration of ions.Changes in the resting membrane potential due to the ionic gradient may also be calculated using the Goldman-Hodkin Katz equation 7 .
Potential future applications of the BD numerical simulation scheme include the study of cell death caused by ion dysregulation during electroporation.The model could be extended to simulate transport through proteins, e.g., voltage-gated ion channels as modeled in work reviewed by Sansom et.al. 57 , as well as the transport of molecular chains.This will be an important step towards predicting the electric-field assisted entry and insertion of macromolecules such as polymers and DNA into cells.

FIGURE 1 .
FIGURE 1. a) Simplified electrical circuit diagram of the electroporator device.b) Applied unipolar exponential decay waveform with time constant τ = RxC

FIGURE 2 .
FIGURE 2. Geometry and dimensions of the simulated region of an electroporated E. coli bacteria, containing the cytoplasmic and external membranes.a1, b1 are the edges on the z axis of the inner membrane pore with radius r1; a2, b2 are the edges of the outer membrane pore with radius r2.Ions are pushed through the pores by the electric field towards the periplasm between the membranes and the extracellular region.Hard boundaries are highlighted in blue.Other boundaries are soft in the cytosol or opened in the periplasm and extracellular region.(Dimensions as in Table2).

FIGURE 7 .
FIGURE 7. Histograms of Ca 2+ , Mg 2+ , K + , Na + and Cl -ions transport time through a 5nm long pore with radius = 1 nm, under an 18 kV/cm applied electric field at temperature T = 298 K.

FIGURE 8 .
FIGURE 8. Histograms of Ca 2+ , Mg 2+ , K + , Na + and Cl -ions transport time through a 5nm long pore with radius = 5 nm, under an 18 kV/cm applied electric field at temperature T = 298 K.

FIGURE 9 .
FIGURE 9. Histograms of Ca 2+ , Mg 2+ , K + , Na + and Cl -ions transport time through a 5nm long pore with radius = 20 nm, under a 18 kV/cm applied electric field at temperature T = 298 K.

FIGURE 10 .
FIGURE 10.Flow cytometry graphs of pulsed electrocompetent DH5a E. coli cells with YOPRO-1.33 uL were taken in 200 uL 10% glycerol (in Dl Water as blank).Plot 6 depicts the nonelectroporated (green) and electroporated w/ YOPRO-1 absorption (blue) populations based on their forward (FSC-A) and side scatter (SSC-A) light intensity profiles.Plot 7 shows a histogram of the permeabilized cell count after gating for various light intensities (FL1-A).

FIGURE 12 .
FIGURE 12. Calculated passive diffusion conductivity after the pulse (solid lines) compared to the average experimental data for a) fresh DH5a E. coli bacteria and b) electrocompetent DH5a E coli grown overnight.Triplicate samples are shown with standard deviation error bars.

Table 1
lists the model variables which are tuned to maximize ion throughput while keeping electroporation reversible.Table2lists the physical, chemical, geometry and electroporation constants used in the calculations of the theoretical Brownian dynamics model.Table3lists the physical parameters for the ions considered in the model.

Table 5 Mean and standard deviation of calculated inner membrane ion crossing time through pores with various radiuses r
()    [0.016  2+   2+ + 0.348  2+   2+ + 0.579  +   + (41) +0.026  +   + + 0.031  −   − ⌋ , 61e experimental patch-clamp technique can better validate the model with measurements taken in the microsecond time scale.Other experimental techniques, such as voltage sensitive dyes, Raman spectroscopy and confocal microscopy, as well as the giant unilamellar vesicle model system might also help validate the model, as in the work of Sozer et.al.forMD simulation validation61.The primary measurements for comparison are conductivity, diffusion, mobility, crossing time, electric field and voltage disturbances at different temperatures within the electroporated region.