Basic Features of a Cell Electroporation Model: Illustrative Behavior for Two Very Different Pulses

Science increasingly involves complex modeling. Here we describe a model for cell electroporation in which membrane properties are dynamically modified by poration. Spatial scales range from cell membrane thickness (5 nm) to a typical mammalian cell radius (10 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu$$\end{document}m), and can be used with idealized and experimental pulse waveforms. The model consists of traditional passive components and additional active components representing nonequilibrium processes. Model responses include measurable quantities: transmembrane voltage, membrane electrical conductance, and solute transport rates and amounts for the representative “long” and “short” pulses. The long pulse—1.5 kV/cm, 100 μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu$$\end{document}s—evolves two pore subpopulations with a valley at ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim}$$\end{document}5 nm, which separates the subpopulations that have peaks at ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim}$$\end{document}1.5 and ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim}$$\end{document}12 nm radius. Such pulses are widely used in biological research, biotechnology, and medicine, including cancer therapy by drug delivery and nonthermal physical tumor ablation by causing necrosis. The short pulse—40 kV/cm, 10 ns—creates 80-fold more pores, all small (<\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$<$$\end{document}3 nm; ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}1 nm peak). These nanosecond pulses ablate tumors by apoptosis. We demonstrate the model’s responses by illustrative electrical and poration behavior, and transport of calcein and propidium. We then identify extensions for expanding modeling capability. Structure-function results from MD can allow extrapolations that bring response specificity to cell membranes based on their lipid composition. After a pulse, changes in pore energy landscape can be included over seconds to minutes, by mechanisms such as cell swelling and pulse-induced chemical reactions that slowly alter pore behavior.

partial body levels using voxels with continuum properties that are either fixed or responsive to local fields (Benkler et al. 2008;Nagaoka et al. 2008;Christ et al. 2010). At the other extreme, there are tremendous advances in modeling voltage-gated protein channels in lipid membranes, based on molecular dynamics (MD) (Shaw et al. 2010;Jensen et al. 2013).

Continuum Planar Membrane and Cell Electroporation Modeling
The first electroporation (EP) models treated artificial planar bilayer membranes, which were also the subject of experimental studies that revealed basic behavior such as reversible electrical breakdown (REB). It began with simple, 1D descriptions of planar bilayer membranes (Litster 1975;Taupin et al. 1975;Abidor et al. 1979;Sugar 1981;Weaver and Mintzer 1981; Barnett and Weaver 1991;Weaver and Chizmadzhev 1996). Cell EP modeling was preceded by passive cell models that established basic features of responses to small fields (Pauly and Schwan 1959;Kotnik and Miklavčič 2000;Foster 2000;Yao et al. 2005;Kotnik and Miklavčič 2006). Modeling has now progressed to provide approximate descriptions of singlebilayer outer plasma membranes (PMs) (DeBruin and Krassowska 1999a, b;Stewart et al. 2004; Krassowska and Filev 2007;Talele et al. 2010;Smith et al. 2013). Some continuum models of cell membrane EP involve explicit representations of additional, intracellular organelle membranes Smith and Weaver 2008;Esser et al. 2010;Retelj et al. 2013). Other multicellular membrane models represent several irregular shaped cells in close proximity that approximate solid tissue Weaver 2003, 2006;Esser et al. 2007Esser et al. , 2009).

Molecular Dynamics (MD) Electroporation Modeling
Starting about a decade ago MD has also been applied to the formation of transient defects (''pores'') in lipid bilayer membranes, due to both mechanical stress and to elevated transmembrane voltages (Tieleman et al. 2003;Tieleman 2004;Tarek 2005;Vernier et al. 2006b;Vernier and Ziegler 2007;Gurtovenko and Vattulainen 2007;Böckmann et al. 2008;Levine and Vernier 2010;Romeo et al. 2013;Bennett et al. 2014). Most MD modeling has not only emphasized electrical creation of pores in small planar membranes, but also treat vesicular and other ''out of plane'' structural rearrangements.

Background for Present Model
Our approach uses continuum modeling. It features simple to complex membrane geometries, fixed (passive), and dynamic (active; transmembrane voltage-responsive) properties. We assign each type of membrane a resting potential source, such that as pores are created, their combined conductance shunts the source, and the transmembrane voltage, D/ m , decreases. We also assign a complex, dynamic local membrane EP model that describes pore creation, evolution (expansion/contraction), and eventual stochastic destruction (Pastushenko et al. 1979;Levine and Vernier 2010;Smith et al. 2013).
Cell level models can be based on relatively simple single-closed curved membranes that represent the outer PM. More complex and realistic cell models include explicit intracellular membranes that represent organelles. With the exception of spherical models that include a concentric inner membrane, others populate the cell interior with one or more organelle membranes at asymmetric locations. Additional realism is found in the use of double membranes that represent the nuclear envelope or the inner and outer mitochondrial membranes Esser et al. 2010;Gowrishankar et al. 2011).

Models of Biological Systems Require Approximations
Modeling is pursued to gain useful insights, including confirmation of expectations but also discovering new, unexpected behavior. We accept the view that a model is a hypothesis (Slepchenko et al. 2003), subject to testing and refutation, and can also lead to improvement of the hypothesis (see ''Future Directions'' section). Given the complexity of real cells and their membranes, we emphasize understanding model responses, related to some aspects of experiments, such as uptake or release of ions and molecules. The complexity of cell modeling also means that we must make approximations, particularly in the view of computational limitations.

Importance of Electrodiffusive Solute Transport Rates
The traditional emphasis of EP is introduction of solutes to the cytoplasm. Now EP is also of interest for removal of solutes from the cytoplasm. Obtaining samples for external analysis is one goal. Killing cells by accidental necrosis (Galluzzi et al. 2012) is another, the basis of irreversible electroporation (IRE) for nonthermal tumor ablation (Davalos et al. 2004;Miller et al. 2005;Rubinsky et al. 2008;Garcia et al. 2011;Golberg and Yarmush 2013). Release of important biomolecules and ions is basic to necrosis.
Further, intrinsic apoptosis (Galluzzi et al. 2012) is the basis of nonthermal cell death by nanosecond electric field (nsPEF) (Beebe, et al. 2002;Nuccitelli et al. 2006;Garon et al. 2007;Nuccitelli et al. 2014). In this case, release of solutes from organelles such as the endoplasmic reticulum and mitochondria is important, but at the same time molecules needed for apoptosis should remain in the cytoplasm, a quantitative issue. Both uptake and release of solutes are therefore relevant to current research topics and applications in EP.

Objective of This Paper
Here we describe a cell model, an extension by Smith (Smith 2011) of the cylindrical version of the traditional Schwan passive cell model (Pauly and Schwan 1959). The simplest forms of the Schwan equation for the spherical and cylindrical cell membranes are given in ''Appendix. '' In the present model, local models for EP have been added, along with a resting potential source. Coupled electrical and poration behavior is used to predict the solute transport by electrodiffusion. Two charged solutes are considered, calcein (green fluorescence) and propidium (red fluorescence), as these fluorescent probes are often used in EP research.
We present the model's response to a relatively small, long pulse with 1.5 kV/cm strength and 100 ls duration that is typical of widely used pulses in cancer tumor treatments to deliver drugs such as bleomycin into cells (electrochemotherapy or ECT) (Poddevin et al. 1991;Dev et al. 2000;Gothelf et al. 2003;Breton and Mir 2012), and in IRE that ablates by causing accidental necrosis by nonthermal pulses alone (Davalos et al. 2004;Miller et al. 2005;Rubinsky et al. 2008;Garcia et al. 2011). Frandsen et al. describe a further extension, calcium EP (Frandsen et al. 2012).
In addition, we present model responses to a relatively large, short pulse of 40 kV/cm amplitude, and 10 ns duration (Silve et al. 2012). Pulses of this general type (submicrosecond, megavolt per m; also nsPEF or nanosecond pulsed electric fields) are of more recent interest (Schoenbach et al. 2001;Müller et al. 2001;Beebe, et al. 2002;Craviso et al. 2012;Nuccitelli et al. 2006;Garon et al. 2007;Napotnik et al. 2012). Overall, the cell system model is complicated, as required to account for observed behavior in cell experiments (Smith 2011;Canatella et al. 2001;Puc et al. 2003).

Methods
We describe our methods in both general and specific terms, with subsections that identify major features. Realistic models are a worthwhile approach to creating hypotheses for cell EP, which involves a very large parameter space that cannot be realistically addressed by experiments alone .

Model Geometry
While some prior continuum models have used toroidal pore geometry (Smith et al. 2004;Smith 2011) or even simple cylindrical geometry (Weaver and Chizmadzhev 1996), we now use a trapezoidal pore (Fig. 1a), which has a cylindrical pore interior. This region of a pore is often rate limiting. We use this interior cylinder with spreading resistance assigned outside, to approximate both a pore's vestibules and the nearby regions outside the membrane adjacent the entrance and exit to the pore. The trapezoidal pore geometry is created by rotating a trapezoid around an axis perpendicular to the local membrane plane (Smith 2011;Smith et al. 2013). Here d m % 5 nm is the average cell membrane thickness (Kotnik and Miklavčič 2000). This is larger than most bilayer membranes because it includes 50 % membrane proteins in the PM. The half membrane thickness, d m =2 is the radius of the rotated circle of the toroidal pore, and r p is the pore radius (size), which can vary significantly. The length of the inner, cylindrical pore interior is d p , here taken to be half the average membrane thickness for simplicity. Figure 1b shows the idealized (cylindrical) geometry of an elongated molecule near a pore. Cylindrical geometry is compatible with well-established engineering approximate descriptions of steric hindrance and also charged solute partitioning based on Born energy changes (Parsegian 1969(Parsegian , 1975Chernomordik et al. 1987;Glaser et al. 1988;Kakorin and Neumann 2002). Figure 1c shows the cell model geometry, with the PM interfaced between the light blue extracellular and dark blue intracellular regions, all within the simulation box, with idealized planar electrodes at the top and bottom. We use the simplest representation of a PM, a closed, curved cylindrical membrane that separates the intracellular region from the external medium. Figure 1d, e shows the full simulation box and the dense meshing near the PM. The cell radius is r m ¼ 10 lm, with cell height is h m ¼ 13:3 lm, so that the cell volume is the same as for a spherical cell with the same radius. Figure 1f shows the pore energy landscape (Smith et al. 2013), similar to previous landscapes. The landscape is populated with pores via a creation rate that injects pores at r Ã p ¼ 0:65 nm. Once within the landscape, pore radii change (expand and contract) in response to drift and diffusion contributions as the local slope of the landscape changes with time as the local D/ m changes. Pores can leave the landscape if they reach a bin just to the left of r Ã p , where stochastic destruction occurs. Both creation and destruction are governed by absolute rate equations, with creation having a strong and nonlinear dependence on the local D/ m .

Creation and Evolution of Pore Populations During a Pulse
Dynamic pore populations (pore size distributions) are an underlying feature of cell membrane EP. An applied field results in complex distributions of membrane fields, which change D/ m ðtÞ for each transmembrane node-pair. The pore populations have a complex dependence on the applied field which couples to the various local membrane sites by field amplification (Fig. SI-2, Supplemental Information). These populations are also the source of the inherent complexity of cell membrane EP. The populations are computed self-consistently with transmembrane charge transport through  Fig. 1 a Top traditional toroidal pore geometry. Bottom present, simplified trapezoidal pore geometry. Here d m is the average cell membrane thickness (5 nm), r p is the pore radius at the most restrictive (narrow) interior pore region. b The trapezoidal pore approximation is used to calculate the electrodiffusive transport of a solute with length, l s and radius, r s through a pore of radius r p and thickness d p . c The circular cell model consists of intracellular (dark blue) and extracellular (light blue) electrolyte, separated by a 5 nm thick membrane with 10 lm radius.
2 Electrical and poration behavior for a 1.5 kV/cm, 100 ls pulse (a-f) and a 40 kV/cm, 10 ns pulse (g-l). Average transmembrane voltage, hD/ m ðtÞi; is spatially averaged over the polar quadrants. Angular transmembrane voltage, D/ m ðhÞ at the end of pulse maximum shows the spatial heterogeneity in the electrical response across the membrane. Light red and black shaded regions represent electroporated regions of the membrane. Similarly, pore density, nðhÞ, at the end of pulse maximum shows the spatial extent of electroporation. The total number of pores, NðtÞ, rises after REB, and decreases exponentially after the pulse in accordance with pore lifetime, s p ¼ 4 s. Similarly, the total membrane conductance, G m ðtÞ, increases with pore creation and expansion (Color figure online) several hundred different dynamic pore populations during a pulse. Examples of pore distributions on the anodic and cathodic sides of the cell are displayed in Fig. 3.

Meshed Network Equivalent to Very Large Circuit
The meshed transport network model (Fig. 1e) is generated by assigning local model to 600 membrane sites and the surrounding electrolyte. These interconnected elements can be regarded as a very large circuit that contains both linear, fixed components (9,267 nodes within intra-and extracellular aqueous regions), and highly nonlinear active components (transmembrane node-pairs). Each of the local areas associated with a transmembrane node-pair is regarded as a very small planar membrane patch endowed with a resting potential source and a complete dynamic EP model. Solutions of the circuit yield potentials at the many extraand intracellular nodes, which allow equipotentials to be computed and displayed throughout all or portions of the system at any time. The pairs of potentials at the transmembrane nodes allow the transmembrane voltage, D/ m , to be similarly computed. These are spatially averaged and displayed as a function of time ( Fig. 2a, g).

System Model has Passive and Active Components
The model is constructed from established passive and active (nonlinear, hysteretic) components. These components range from simple conductive and dielectric (capacitive) models to a complete membrane EP model for a small planar membrane patch (Smith 2011). The most complicated local models are located at the membrane, here inserted between each of 600 node-pairs. These models account for changing pore energy landscapes and two rate equations (pore creation and destruction rates) that govern pore population evolution based on pore creation, expansion and contraction, and destruction. Dynamic pore behavior is computed at each transmembrane node-pair, based on the local D/ m and the pore population. During a pulse the pore energy changes at each of the 600 node-pairs. This means that the pore energy gradients change, leading to changing drift and diffusion rates in pore radius space at each of these transmembrane sites (Smith et al. 2013). Overall, the interconnected system of local models is complex, needed to describe the many interactions that together define cell membrane EP. The different parameters and their values used in the system model are listed in Table 1.

Solute Transport by Electrodiffusion Within Bulk and Through Dynamic Pores
The important topic of specific ions and molecules (solutes) transport is addressed by electrodiffusion in the transport number approximation, which utilizes already-determined electric fields to compute electrodiffusion (Smith 2011;Smith and Weaver 2012). Two types of regions are treated simultaneously, via interconnections as in circuit theory. The first case is bulk electrodiffusive transport within the extra-and intracellular media with electrodiffusion through the individual pores of different pore populations that are located between the transmembrane node-pairs. The detailed description of individual pore solute electrodiffusion includes partitioning into pores and sterically hindered movement through the interior pore region (Smith 2011). Solute transport through many complex dynamic pores in the membrane is both rate limiting and computationally nontrivial, and is intimately connected to and within the bulk aqueous regions.
Here we consider two fluorescent probe molecules used in EP research, calcein and propidium (the cation of dissociated propidium iodide, PI). Calcein (623 Da) has r cal ¼ 0:58 nm and l cal ¼ 1:89 nm and a charge of À4. Propidium (415 Da) has charge number z = ?2, radius r pro ¼ 0:69 nm; and length l pro ¼ 1:55 nm (Smith 2011). The number of molecules within the cell is treated as free (dissolved). For simplicity, we do not represent the binding of propidium to intracellular nucleic acids (Sadik et al. 2013).

Idealized and Actual Pulse Waveforms
To aid interpretation of results, we used idealized trapezoidal pulses, with linear ramps that define the rise and fall Propidium radius (r pro ) 0.69 nm (Smith 2011) Propidium length (l pro ) 1.55 nm (Smith 2011) Propidium initial concentration (c pro;ext;0 ) 1mM Temperature (T) 298.15 K times, and a flat peak field strength ( Fig. SI-1). The model can also be used with digitized experimental waveforms (Gowrishankar et al. , 2011.

Conventional Electroporation Pulse
For mammalian cells most pulses have strengths of 0.1 to $ 10 kV/cm and durations longer than a microsecond. As an example, we use an idealized trapezoidal pulse with strength 1.5 kV/cm, duration 100 ls; and a rise/fall time of 1 ls, an approximation to a widely used pulse. Pulse trains that employ eight 100 ls pulses were first used in research (Mir et al. 1988), and subsequently in electrochemotherapy (ECT) at a repetition rate of 1 Hz (Poddevin et al. 1991;Dev et al. 2000;Miklavčič et al. 2005), in IRE (Davalos et al. 2005;Miller et al. 2005;Rubinsky et al. 2008;Rubinsky 2010;Garcia et al. 2011;Breton and Mir 2012;Kennedy et al. 2014) and in calcium electroporation (Frandsen et al. 2012).

Submicrosecond, Megavolt per Meter Pulse for Supra-Electroporation
The second pulse example is also an idealized, trapezoidal pulse, with strength 40 kV/cm, duration 10 ns, and linear ramps that define the 1 ns rise and fall times. Although terminology varies, a common designation is nsPEF (nanosecond pulsed electric field), indicating pulses with durations less than 1 ls. Such pulses have received considerable attention since 2001, when two papers appeared several months apart (Schoenbach et al. 2001;Müller et al. 2001). These pulses are close to the conditions studied by MD (Tieleman et al. 2003;Vernier et al. 2006b;Vernier and Ziegler 2007;Levine and Vernier 2010;Fernández et al. 2012;Romeo et al. 2013;Delemotte and Tarek 2012), and involve pore densities expected for supra-EP (Stewart et al. 2004;Gowrishankar et al. 2006;Esser et al. 2009). Potential applications include nonthermal cell death by apoptosis (Beebe, et al. 2002;Beebe, et al. 2003;Nuccitelli et al. 2006;Garon et al. 2007;Nuccitelli et al. 2014), and a variety of other intracellular effects. Even the introduction of bleomycin is reported for nsPEF (Silve et al. 2012).

Results and Discussion
Electrical and Poration Behavior Figure 2 shows examples of coupled electrical and poration behavior. Movement of small ions through pores and toward or away from the membrane is used to compute currents, fields, and transmembrane voltages, with the latter two quantities emphasized here. Charged solute transport contributes negligibly to electrical behavior. Accordingly, solute transport by electrodiffusion is computed using the electric fields determined solely by small, ubiquitous ions.

Transmembrane Voltage Temporal Behavior
Electroporation is characterized by electrical creation of transient aqueous pores, usually described by a rate equation that is strongly dependent on the transmembrane voltage, D/ m (Weaver and Chizmadzhev 1996;Neu and Krassowska 1999;Smith 2011;Smith et al. 2013). At one extreme, a pore can be created in about one nanosecond for D/ m % 1:2 V in a planar membrane patch (Vasilkoski et al. 2006), with similar times reported by MD studies using pure lipid bilayers (Gurtovenko and Vattulainen 2007). For longer times much smaller voltages are sufficient (Melikov et al. 2001). Cell models show behavior consistent with similar rate equations (DeBruin and Krassowska 1999a, b;Stewart et al. 2004;Krassowska and Filev 2007) and include the present model (Smith 2011;Smith et al. 2013). The transmembrane voltage, D/ m , is a key quantity for pore creation, and for dynamic pore behavior involving expansion and contraction, with the latter important to pores shrinking to a small enough size that they are destroyed by thermal fluctuations (Pastushenko et al. 1979;Levine and Vernier 2010;Smith et al. 2013). Figure 2a, g shows D/ m averaged over the cathodic and anodic quadrants for the two pulses, with D/ m (t) behavior that is quite different due to the very different durations and rise/fall times. Figure 2a shows the response to the conventional EP pulse (1.5 kV/cm, 100 ls). The display emphasizes the transient behavior at the beginning of the pulse (first 3 ls) and at the pulse end (from 99 to 101 ls). The membrane charges passively until t ¼ 0:7 ls; when the rise in D/ m is interrupted by REB (Benz et al. 1979;Barnett and Weaver 1991). This transmembrane voltage peak of 1.2 V is associated with a burst of pore creation. The accompanying large conductance increase interrupts the D/ m rise, and forces a change in direction such that D/ m quickly decreases, even though the applied field is still increasing. Within a microsecond after the REB peak the rate of D/ m decreases slowly, and reaches a quasi-plateau of about 0.5 V. Figure 2b, k shows the angular profiles of D/ m at the end of the maximum values of the two pulses. For the 100 ls pulse this is 99 ls, which is much longer than five charging times ðs m ¼ 370 ns). But the idealized 10 ns pulse peak value ends at 9 ns, which is much shorter than s m . This means that pore creation begins rapidly, at about 5 ns (Fig. 2j). At D/ m % 1:7 V (Fig. 2g) pores are made rapidly, even though the membrane never charges fully. Instead, a large local conductance exists when the pulse maximum value ends at 9 ns, and D/ m begins to relatively slowly decrease to zero over about 40 ns (Fig. 2g).

Transmembrane Voltage Angular Profiles
Flattening of the transmembrane voltage profile around a cell is a classic indicator of EP, usually shown as the angular profile of D/ m . Flattening is observed experimentally (Kinosita et al. 1988;Frey et al. 2006;Flickinger et al. 2010), and is also exhibited by cell membrane models and has been found in cell model responses for both spherical (DeBruin and Krassowska 1999a; Smith and Weaver 2008;Talele et al. 2010;Sadik et al. 2013), and cylindrical membranes (Stewart et al. 2004;Gowrishankar et al. 2006;Smith et al. 2013). Figure 2b, h shows the model's angular response for the two pulses. In addition to profile flattening around the poles, the local peaks in Fig.  2b indicate transition regions that separate porated (shaded regions) from unporated membrane (DeBruin and Krassowska 1999a; Krassowska and Filev 2007;Esser et al. 2010).

Pore Densities
The model's poration response is shown as pore density plots (Fig. 2c, i) for both pulses, and shows where pores are present, independent of their size at the time indicated.
Here the pore density (number of pores per square lm) is indicated by white bands superposed on the cell membrane. The width of the white bands is proportional to pore density. Figure 2c shows two thin white bands that occupy about a third of the membrane, with 2:7 Â 10 4 pores (Fig. 3a-d) created by a 1.5 kV/cm, 100 ls pulse. The time is the end of the pulse peak, 99 ls, when pores have evolved to form two subpopulations (Fig. 3c). Similarly, Fig. 2i shows two broad white bands that occupy more of the membrane, with 2:1 Â 10 6 pores ( Fig. 3f, h) created by a 40 kV/cm, 10 ns pulse. The time is also at the end of the pulse peak, 9 ns. For this much larger and shorter pulse, there is very little pore expansion ( Fig. 3f-h). All the pores are small, but there are about 100-fold more pores.

Total Pore Number
The total number of pores in the cell membrane, N, is another partial indicator of EP. A general expectation is that more pores provide greater transport capability for small ions that dominate electrical behavior and for small solutes. According to the model, NðtÞ rapidly increases during a pulse, but after a pulse NðtÞ decreases exponentially with a 4 s pore lifetime (recovery or resealing time). This general behavior is a consequence of the importance of electrical interactions during a pulse. However, membrane depolarization results for most pulses, resulting in  Fig. 3 Pore histograms show the distribution of pore radii, Nðr p Þ, at various times in response to the 1.5 kV/cm, 100 ls pulse (a-d) and the 40 kV/cm, 10 ns pulse (e-h). The maximum pore radius r max ¼ 60 nm, but the display has been truncated to 40 nm in consideration of the maximum pore radii observed D/ m % 0 post pulse, with pore destruction proceeding spontaneously, due to local fluctuations. This is consistent with MD simulations, but MD lifetimes are much shorter, for reasons not yet understood.
For the conventional EP pulse, Fig. 2d shows the initial increase in total pores, which occurs in less than a microsecond, well before the pulse has reached its peak value. On a much longer time scale Fig. 2e shows N both during the pulse and for five pore lifetimes (20 s total) of stochastic pore decay post pulse. During essentially all of the post pulse decay, there are only thermalized pore distributions (Fig. 3d, h). The shape of the thermalized distribution is narrow, collapsed to only small pores. This narrow distribution is on average retained, but the magnitude, NðtÞ, of the population decays exponentially (Fig. 2e,  k).
The response to the nsPEF pulse ( Fig. 2j) also shows a rapid initial increase in total pores. In this case, the pulse ends (at 9 ns) before REB occurs, but already about two million pores have appeared, all small ( Fig. 3f-h). Both the small pores and their large number are features of supra-EP (Stewart et al. 2004;Smith and Weaver 2008;Joshi and Schoenbach 2010), which is also predicted to porate many organelles. The longer time scale result again shows a sharp increase in N, followed by the same exponential decay (Fig. 2k), but the magnitude of N is two orders of magnitude larger.

Membrane Conductance
The cell membrane conductance, G m , is closely related to N for the simple reason that the predominant small ions (Na þ ; K þ and Cl À ) are so small that they easily pass through even small pores. A D/ m -dependent partitioning of these ions into pores is expected which gives G m an obvious dependence on N, but also on D/ m (Chernomordik et al. 1987;Glaser et al. 1988;Kakorin and Neumann 2002). The total conductance is important, because almost independent of how it is spatially distributed a large G m is effective in shunting the resting potential source (DeBruin and Krassowska 1999a;Smith et al. 2013). For this reason most porating pulses lead to post pulse membrane depolarization. The shunting of that source by the 1,000 or more pores created by most pulses leads to D/ m % zero after a pulse.

Pore Population Life Cycles: Creation, Evolution and Destruction
Pore populations, also termed pore size distributions, yield more detailed insights beyond NðtÞ. Here we show examples of pore populations created by the two pulses. For both pulses, the pre-pulse condition is an average of two pores total, one on each side (Fig. 3e). The most interesting phase of a population life cycle is from the beginning of additional pore creation by an applied pulse, to the time where the pulse ends. After that there is a general pore contraction that leads to a thermalized pore population in about 10 ls (not shown). The pore contraction depends on discharge of the membrane, such that D/ m goes to zero somewhat after the pulse end (Fig. 2a, g).
The thermalized populations have the same shape on average (Fig. 3d, g, h). Such distributions arise from thermal equilibrium and the shape of the pore energy, W, around the minimum (the 0 V curve of (Fig. 1f)) (Neu and Krassowska 1999;Smith 2011;Smith et al. 2013). This collapsed, post pulse population then decays with a relatively long pore lifetime (here 4 s) that is believed to depend on cell membrane composition and temperature. While these are common features before and after a porating pulse, the detailed population behavior is quite different for the two pulses. Long porating pulses lead to populations with a wide range of sizes ( Fig. 3a-d), whereas short pulses do not (Fig. 3e, h). According to essentially all theories and MD simulations, once created, pores tend to evolve to complex size distributions, while D/ m is elevated, and only vanish when they contract to a minimum size. A detailed discussion of pore dynamics follows.

Basal Rate and Spontaneous Pore Creation
The basal rate of pore creation is not zero at the resting potential (here -50 mV). The result of a field pulse is to increase an existing rate, not to start a new process. Accordingly, even a fully depolarized membrane ðD/ m ¼ 0) has a very small but nonzero spontaneous pore creation rate, due to thermal fluctuations. Present theories and models use an absolute rate equation that is very nonlinear in D/ m , which results in large numbers ([1,000) of pores for typical pulse strengths and durations. From this it also follows that there is no critical field or transmembrane voltage. Instead, both an elevated D/ m and a finite time are involved in producing a detectable number of pores. Accordingly, the model's pre-EP condition of N [ 0 is consistent with the idea that nonionizing fields only modulate ongoing processes in biological systems (Vasilkoski et al. 2006).

Local, Spontaneous Pore Destruction
In the present model, stochastic pore destruction occurs at r Ã p ¼ 0:65 nm (Fig. 1f). Most pores are destroyed after a pulse, removed stochastically from thermalized distributions at D/ m % 0 (Fig. 3d, g, h). In pore radius space this means that some pore radii move over the landscape (Fig.  1f) on the 0 V curve to reach r Ã p , where they become candidates for random removal by local fluctuations, according to an absolute rate equation. Pores expand and contract within the changing (D/ m -dependent) landscape by both drift (force) and diffusion (random walk). Mathematically, this is analogous to electrodiffusion, and is solved using the descretized method employed for solute electrodiffusion (Smith and Weaver 2012).
One interpretation is that ''pore lifetime'' may be explained by involving two distinguishable processes: (1) the analog of electrodiffusion within the landscape, and (2) stochastic destruction that is local, here at slightly below r Ã p . That is, pore loss is due to both the local stochastic destruction rate and transport in r p space to the destruction site slightly below, r Ã p . In one case, the destruction process could be mostly governed by diffusion, analogous to diffusion-limited chemical reactions. Alternatively, MD studies may show that the final stages of contraction and destruction comprise a single, multi-molecule event.

Supra-Electroporation Pore Populations by nsPEF
Because of its relative simplicity, we begin a brief detailed discussion of the pore population examples by considering the 40 kV/cm, 10 ns pulse response (Fig. 3e-h). The pores are small and the population behavior is limited. Figure 3e is early (t = 1 ns), and can be regarded as the pre-EP condition (Fig. 2j). Both the top (anodic) and bottom (cathodic) display the total pore number summed over the respective cell side. In Fig. 3e this is 1.03 pores on each side. This is the equilibrium number of spontaneous pores generated by equilibrium between two different, nonzero rates: spontaneous pore creation and stochastic pore destruction. The occasional spontaneous pore creation is followed by random destruction with a mean lifetime of 4 s. Together, these rates of creation and destruction are the basis of N þ and N À ¼ 1:03 on the anodic and cathodic sides. The average total number of pores is thus N % 2, a small value which guarantees major fluctuations.
The nsPEF example shows that many pores exist at 10 ns, when the field pulse drops to zero. D/ m peaks at a magnitude of about 1.7 V on both the cathodic and anodic sides, at 9 ns, the end of the field pulse maximum. This abrupt interruption of membrane charging nevertheless results in many pores, with slightly different numbers on the anodic (N þ ¼ 1:08 Â 10 6 pores) and cathodic sides ðN À ¼ 9:99 Â 10 5 pores), so the total pore number is N ¼ 2:08 Â 10 6 . Here we show three significant figures, not because the model is that precise, but to make comparisons easier.
The main point of Fig. 3f-h is that there are many pores, all small due to the short pulse duration, only during which pores can expand beyond r min . In addition, the highly nonlinear dependence of pore creation on D/ m leads to many pores. The resulting large membrane conductance holds down D/ m (Frey et al. 2006), to values far less than predicted for passive (nonporated) membranes (Stewart et al. 2004;Kotnik and Miklavčič 2006). Figure 3f-h shows narrow pore size distributions as expected for supra-EP (Stewart et al. 2004;Smith and Weaver 2008), which is a logical consequence of the model's construction (Methods), and are consistent with experimental electrical behavior (Frey et al. 2006) and observed fluorescent probe uptake (Vernier et al. 2006a;Pakhomov et al. 2007aPakhomov et al. , b, 2007a and also bleomycin (Silve et al. 2012). Existing experimental and MD modeling results are consistent with the present model.
At 10 ns (pulse end), the pore population is represented by histograms with only two occupied bins (0.5 nm width), with a total of about 2:1 Â 10 6 pores (Fig. 3f). There is a slight asymmetry, with about 8 % more pores on the anodic than cathodic side. There are also about 1,000-fold fewer pores in the larger bin. At t ¼ 1 ls (Fig. 3g) the pulse has long ended, and the population has relaxed into the thermalized post pulse distribution (Fig. 3h), which but for the total pore number has the same distribution shape as Fig. 3d. A comparison of the populations in Fig. 3d, g, h shows the same pore size distribution shape, which means that each of these populations is in thermal equilibrium around r p;min for 0 V in the landscape of Fig. 1f (Neu and Krassowska 1999;Smith 2011;Smith et al. 2013). The similar shape of these thermalized distributions occurs because D/ m ¼ 0, and for our model only the 0 V curve of the landscape matters. The properties of such post pulse distributions are important to the post pulse diffusion of moderate size solutes.

Conventional Electroporation Pore Populations
Figure 3a-d shows the more complex pore population behavior of the longer 1.5 kV/cm, 100 ls conventional EP pulse. Pores are created early in the pulse, and subsequent interactions within the model allow the emergence of a subpopulation of large pores. Emergence of such subpopulations has now been reported for three related cell EP models (Krassowska and Filev 2007;Esser et al. 2010;Smith et al. 2013). The time for emergence appears to be several microseconds. Pore creation essentially ceases (Fig.  2d) near the REB peak at $ 0:7 ls, with average D/ m falling to a plateau of about 0.5 V (Fig. 2a). For this transmembrane voltage insignificant additional pore creation occurs on a 100 ls time scale.
When the pulse reaches its peak value at 1 ls, the distribution of pore radii has already broadened (Fig. 3a). The anodic distribution extends to 10 nm radius pores, and the cathodic population to slightly less at 9 nm. At 9 ls, about 10 % of the pulse duration, the pore size distribution has become more complex, with Fig. 3b showing the pore distribution beginning to divide into two subpopulations. Both the anodic and cathodic populations have become complex. There is a large peak at $ 1.5 nm, and a barely discernible local maximum at $ 8 nm for both the anodic and cathodic sides. The pore size distribution at 1 ls is slightly after the REB peak, when pore creation is complete (Fig. 2a, d). The total pore number is then N ¼ 2:7 Â 10 4 , unchanged as the population evolves during the pulse (Fig.  3a-c). Some asymmetry occurs due to the resting potential source, with somewhat more pores on the cathodic (black histogram) than anodic (red histogram) side. The largest pores ( $ 11 nm radius) appear on the anodic side, while the cathodic side has a maximum pore size of $ 9 nm.
At 99 ls, the end of the pulse maximum, two distinct subpopulations have emerged on each side of the cell (Fig. 3c), but the number of pores on each side has not changed. This is consistent with insignificant pore creation for D/ m % 0:5 V (Fig. 2a) on this time scale. Instead pore expansion and contraction occur, due to gradients within the pore energy landscape ( Fig. 1f) (Smith et al. 2013). Finally, well after the pulse end at 110 ls the population has collapsed to form a thermalized pore distribution. This occurs at 10 ls after the pulse end (Fig. 3d). When the field pulse reaches zero approximately ten microseconds are needed for the cell membrane to discharge and for the pore populations to collapse by individual pore contraction until the thermalized distribution is achieved. With the exception of scaling by total pore number, N, the distributions in Fig. 3d, h are the same.
Overall, Fig. 3 illustrates a general feature of the model's response. Complex pore populations can develop during a pulse, but after the pulse the membrane conductance shunts the resting potential source, and the membrane is depolarized ðD/ m % 0Þ. In pore radius space pores contract toward r p;min % 1 nm, by drift, due to the slope in W at D/ m ¼ 0 (Fig. 1f). This collapse into a quasi-equilibrium distribution takes about 10 ls. If the scaling factor of N is taken into account, all of the post pulse thermalized populations are the same, within fluctuations. Depending on pore lifetime (resealing or recovery time), post pulse solute transport occurs only through this type of pore distribution. The initial number of post pulse pores, and their lifetime, determines post pulse transport.

Electrodiffusion of Solutes Within the System Model
The model describes transport of ions and molecules by a descretized representation of electrodiffusion (Smith and Weaver 2012). Figures 4, 5 give examples of solute concentration distributions at different times for calcein (green) and propidium (red) for the conventional and supra-EP pulses, respectively. In both cases, the initial extracellular concentration is 1 mM. During a pulse, the drift contribution usually dominates. After a pulse the membrane is quickly depolarized and diffusion dominates. For a solute to enter a pore its radius r s should be less than r p , the radius of the restrictive pore interior. This leads to a ''soft cutoff'' for transport as a function of size. It should also be energetically favorable for a charged solute to enter, which in the model involves a partition coefficient estimate based on electrostatic energy (Born energy) (Parsegian 1969(Parsegian , 1975. Post pulse diffusion should therefore favor propidium over calcein, in agreement with the model (Figs. 4, 5). Inside the pore solute movement is sterically hindered (Smith 2011).
To understand post pulse transport, a basic question is: What size and charge solutes pass through the thermalized distributions? We address this by providing a known amount of a particular solute within the extracellular region, and then applying a pulse. We know that small, biologically important ions readily pass through the decaying pore population. Water and small sugars can enter and diffuse through these pores, as can charged biologically active solutes such as ATP (507 Da, mainly charge number z = -4 and -3) (Smith 2011). Many membrane integrity probes can also pass through these pores, based on the cylindrical representation of their size (solute radius, r s , is most important; charge number also affects partitioning into a pore) (Smith 2011). Colorimetric and fluorescent dyes are already widely used in cell membrane electroporation research. Descriptions of transport of two medium size charged probes, calcein and propidium, are shown in Figs. 4,5,6,7. Calcein Uptake for the 1.5 kV/cm, 100 ls Pulse Figure 4a-d shows calcein concentration (spatially distributed green) and the number of internalized calcein molecules (white number inside cell) at the indicated times for the 1.5 kV/cm, 100 ls pulse. The difference in intracellular calcein numbers between Fig. 4a (pulse start; t = 0 s) and Fig. 4a (pulse end; 100 ls) shows that electrodiffusive drift dominates, with 3:2 Â 10 7 molecules internalized at 100 ls. There is only $ 3 % additional transport during a post pulse interval of five pore lifetimes (20 s; Fig.  4d) with 3:3 Â 10 7 molecules. This is consistent with inhibition of calcein partitioning into the pore interior, due to its larger average charge number.
The intracellular green color distribution changes significantly between Fig. 4a, b. At the pulse start the intracellular space is black, indicating concentration,  1 nM), i.e., essentially zero (Fig. 4a). At pulse end ð100 lsÞ, calcein concentration appears as a green fuzzy curved band near the membrane on the cathodic side of the cell (Fig. 4b). This is a diffusively broadened distribution of calcein that has entered mainly by drift from the cathodic side, expected for negatively charged calcein. Close examination reveals a much thinner green band near the membrane around the anodic pole, with two small green areas away from the anodic polethese areas are transition regions that favor large pores (DeBruin and Krassowska 1999a; Krassowska and Filev 2007;Smith 2011). Small amounts of calcein can enter pores and diffuse through the membrane as D/ m ramps down before the pulse end. This type of detailed behavior is expected from the model. Internalized calcein is redistributed by diffusion within the cell, which smooths the intracellular concentration distribution. This is shown by the progression of fuzzy green within the cell at 10 ms (Fig. 4c). The difference between Fig. 4b at 100 ls (intracellular green near membrane, particularly on the cathodic side) and Fig. 4c (intracellular green spread out by diffusion to clearly mark two anodic locations associated with transition regions) is striking. The visual result of transition regions is a fuzzy upside down black ''T'' that occupies the upper half of the cell. By the end of the simulation at 20 s (Fig. 4d), diffusion has smoothed the calcein distributions, so the extraand intracellular regions are uniformly green. But the total uptake of n cal ¼ 3:3 Â 10 7 calcein molecules corresponds to a concentration somewhat smaller than the extracellular concentration, which is why the cell interior is darker green than exterior (Fig. 4d).
Overall, calcein enters by electrodiffusion through the dynamic pore populations (Fig. 3a-d). No significant calcein concentration gradients are seen near the cell membrane on either side. This means that the calcein gradient is mostly within the membrane, consistent with even the porated membrane presenting the dominant barrier to solute transport. Also, consistent with a significant porated membrane barrier, calcein uptake results show that only slight additional transport occurs by diffusion after the pulse (Figs. 4c, d, 6).
Propidium Uptake for the 1.5 kV/cm, 100 ls Pulse  Fig. 6 Electrodiffusive transport rates of propidium and calcein, dn s =dt, into the cell for a 1.5 kV/cm, 100 ls pulse on two different timescales (a, b). Cumulative solute influx, n s ðtÞ (c, d) shows the total amount of solute inside the cell as a function of time. During the 1.5 kV/cm, 100 ls pulse, 3:2 Â 10 7 molecules of calcein and 1:5 Â 10 7 molecules of propidium are transported into the cell. After the 100 ls pulse, 3:0 Â 10 5 molecules of calcein and 2:1 Â 10 7 molecules of propidium are transported into the cell band on the anodic side (top), and only a small hint of entry on the cathodic side (bottom). Here 1:5 Â 10 7 molecules are internalized, about half that for calcein (Fig. 4b). This is consistent with drift-dominated electrodiffusion during the pulse. While nearly the same molecular mass, propidium has a charge number magnitude, jzj of 2, half that of calcein, consistent with accumulating about half the amount as calcein during the pulse.
Strikingly, comparison of Fig. 4f at 100 ls , the pulse end, with Fig. 4g at 10 ms, show no change during this small portion of the post pulse time. This interval is only 5 Â 10 À4 of the 20 s post pulse interval. Figure 4 shows that during the 100 ls pulse 1:5 Â 10 7 molecules are internalized, but that during the entire post pulse period an additional 2:1 Â 10 7 molecules are taken up. For the model's 4 s pore lifetime, less than 1 % of the pores remain at 20 s, a time at which post pulse transport is essentially complete. The cumulative transport of 3:6 Â 10 7 molecules at 20 s (Fig. 4h), is a significant increase from 1:5 Â 10 4 molecules at 10 ms (Fig. 4g). Figure 4d shows about the same number ð3:3 Â 10 7 moleculesÞ for calcein.
The model shows that transport mechanisms for calcein and propidium differ significantly for the 1.5 kV/cm, 100 ls conventional EP pulse. Almost identical total uptake occurs for calcein ( Fig. 4a-d) and propidium (Fig. 4e-h), but is accomplished differently. Calcein is internalized almost entirely by drift during the pulse, with only slightly more (1 Â 10 6 molecules) after 100 ls. In contrast, propidium is internalized about half by drift during the pulse and half by diffusion after the pulse. See also Fig. 6.
Calcein Uptake for the 40 kV/cm, 10 ns Pulse The nsPEF pulse molecular uptake response is quite different, revealed by both uptake rates (Fig. 7a, b) and cumulative amounts (Fig. 7c, d). A general prediction is that nsPEF creates many more but small pores (Fig. 3e-h) than conventional EP pulses (Fig. 3a-d). Insignificant transport occurs during the pulse due to its very short duration, and post pulse diffusion accounts for essentially all transport (Smith and Weaver 2011).  Fig. 7 Electrodiffusive transport rates of propidium and calcein, dn s =dt, into the cell for a 40 kV/cm, 10 ns pulse (a, b). Cumulative solute influx, n s ðtÞ (c, d) shows the total amount of solute inside the cell as a function of time. During the 40 kV/cm, 10 ns pulse, 1.5 molecules of calcein and 0.17 molecules of propidium are transported into the cell. After the 10 ns pulse, 1:3 Â 10 7 molecules of calcein and 7:9 Â 10 8 molecules of propidium are transported into the cell 1969Smith 2011). Nevertheless some calcein does enter by partitioning, and then moves through pores by sterically hindered diffusion (Smith 2011;Smith and Weaver 2011). As the thermalized pore population (Fig.  3g, h) decays there is a diffusive calcein influx that also decays with a 4 s lifetime, leading to cumulative intracellular calcein uptake of 4:3 Â 10 3 ; 4:5 Â 10 5 ; and 1:3 Â 10 7 molecules at the post pulse times shown in Fig. 5a-d.
Propidium Uptake for the 40 kV/cm, 10 ns Pulse Uptake of propidium is similar, but larger, owing mainly to its twofold smaller charge number and the dominance of post pulse diffusion. Figure 5e shows zero internalized molecules at the pulse end (10 ns). The subsequent entry is similar qualitatively but not quantitatively. The times shown are the same as for calcein (Fig. 5a-d), but the cumulative propidium uptake numbers are larger: 0, 3:3 Â 10 5 ; 3:4 Â 10 7 ; and 7:9 Â 10 8 molecules. By 20 s inspection shows that the intra-and extracellular concentrations of propidium are almost identical (Fig. 5h). The model, however, does not take into account binding of propidium to numerous binding sites within the cell (Sadik et al. 2013). Our present model therefore underestimates propidium uptake, because if binding was included then the intracellular concentrations (free solute) would be smaller, and the gradient across the membrane would be larger. In our present model, the uptake of both calcein and propidium is not that different. The more interesting observation is the relative contribution during a pulse (D/ m is elevated, supporting drift) and after (D/ m ¼ 0, supporting diffusion over a much longer time).

Total Cell Transmembrane Solute Transport Rates
During a pulse, both D/ m ðtÞ and electric fields within the bulk media are elevated. At each membrane node-pair the local D/ m creates heterogeneous local fields near the entrance/exits of individual pores, and also within the pore interior. Both are included in the model. These pore-associated fields mostly determine electrodiffusive solute transport, as the fields in bulk away from the membrane are much smaller. That is, electrodiffusion through the individual pores of each population is rate limited, with bulk media fields serving mainly to supply and removing solute from the membrane regions. Even for the highly mobile small ions that govern D/ m during supra-EP the porated membrane is the main barrier.
After a pulse the membrane becomes depolarized in about a microsecond. Post pulse solute transport occurs by diffusion that is sterically hindered within pores of a thermalized pore population. Entry into pores is governed by partitioning, which increases as solute charge decreases. Within simple bulk media away from the membrane partitioning is a nonissue, and both diffusion and drift are unhindered.
Calcein and Propidium Influx Rates for the 1.5 kV/cm, 100 ls Pulse Figure 6a shows the calcein influx rate (green dot-dash curve) jumping to $ 2 Â 10 8 molecules/s in a few microseconds. But at 100 ls (pulse end), the rate plummets by more than four orders of magnitude (Fig. 6a), as D/ m drops from $ 0.5 V to zero (Fig. 2a). The entire cell membrane is then depolarized. With D/ m ¼ 0 solute drift ceases, and post pulse solute transport is by diffusion through the thermalized pore population (Fig. 3d). Similar behavior occurs for propidium (red dot-dash curve). The driftdominated propidium transport during the pulse is slightly smaller owing to its smaller charge number, but is larger post pulse due to more partitioning into pores. Figure 6a, b both show a rapid rate decrease at the pulse end at 100 ls. After the pulse there is transient, complicated behavior while the membrane is discharging, and the solute transport is transitioning from drift-dominated to diffusion dominated. This transition interval involves general pore contraction (Fig. 3c, d) as D/ m decreases in two phases from 99 to $ 101 ls (Fig. 2a). After this brief transition, both calcein and propidium are transported over the next 20 s by an exponentially decaying diffusion rate (Fig. 6a, b). The very different rates during and after the pulse are readily seen in Fig. 6b, where the red/green spike (drift) near t = 0 is followed by two parallel dot-dash lines (diffusion) on this semi-log plot. The details of partitioning and steric hindrance are used in all cases of electrodiffusion through pores of different sizes in each population, and these rates appropriately match those in the nearby bulk media, a basis for the model's solutions being computationally nontrivial.
Calcein and Propidium Influx Rates for the 40 kV/cm, 10 ns Pulse The 40 kV/cm, 10 ns pulse has very different solute transport rates (Fig. 7a, b). The calcein and propidium rates rapidly increase from 1 molecule/s (log scale constraint) to a peak values of order 10 9 molecules/s at about 9 ns (peak pulse value; when total pore number, N, is reached Fig. 2j), and then falling by several orders of magnitude by 400 ns. As shown by the behavior in Fig. 7b-d, these huge rates for very short times are unimportant, given the overwhelming contribution of post pulse diffusion (Fig. 7c, d).
Calcein and Propidium Cumulative Uptake for the 40 kV/ cm, 10 ns Pulse Pore creation is complete within only $ 9 ns (Fig. 2j). In this very short time drift dominates for both calcein and propidium (Fig. 7c), with diffusion taking over after a transition involving pore contraction and membrane depolarization. Although there are noticeable transport changes between 9 and about 300 ns (Fig. 7b), the contribution to cumulative transport in this interval is negligible, being overwhelmed by the post pulse diffusive contribution (Fig. 6e-h). The importance of post pulse diffusion is consistent with cumulative calcein uptake of 1:3 Â 10 7 molecules, about 70-fold less than the 7:9 Â 10 8 molecules for propidium, which can partition more readily into the interiors of the small pores that comprise the post pulse thermalized pore populations. The predicted dominance of diffusion for total uptake by nsPEF (Smith and Weaver 2011) is consistent with the model's behavior.
The detailed discussion of solute uptake rates and amounts show that cell membrane EP behavior is complex. This complexity is a result of many interacting membrane sites, using a model for dynamic EP at each of the 600 sites. Once established cell EP models with irregular and multiple membranes can be used with confidence, if comparison to experiments are also pursued.

Future Directions
Modeling should confront three phases of EP phenomena, which involve different time scales and also different types of membrane interactions. The first is pore creation, which can be approached by both continuum models and MD simulations. The second is pore population evolution based on interactions that tend to expand or contract individual pores within a population. The third is post pulse behavior, when the membrane is believed fully depolarized shortly after a pulse. In the absence of an applied field, pore behavior may be affected by other interactions such as intracellular field sources that extend through the porated cell membrane, time-dependent membrane tension from cell swelling, and lipid alteration by interaction with products of stimulated reactions (e.g., reactive oxygen species; ROS) due to EP.
Modeling during the first and second phases involves large fields, and causes changes in D/ m (t). Modeling during the third phase is much more challenging, and has received little attention. Fields and transmembrane voltages are small to zero. Candidate interactions are of different types (mechanical, chemical, and thermal), and may involve specificity based on pore constituents and also the transported ions and molecules (solutes). For this reason, and outlined below, we identify new approaches that can extend the present model's capability to address these problems.
The present model is based on well known passive electrical properties of aqueous media and the fixed portions of cell membranes. Typical values from the literature are used. More composition-specific properties could be used, but the present state of quantitative information from cell EP experiments is not good enough to support detailed testing and validation. Nevertheless, as noted above this and other members of a family of cell models can account for several key experimental observations. Below we briefly note extensions of the present model which should bring improved realism and quantification.

Extrapolation of MD Structure-Function Results to Cell Models
Continuum cell EP models developed over several decades feature homogeneous representations of membrane properties. Even if more than one membrane is included to represent one or more organelles, the individual membranes at most have different homogeneous properties for the unaltered membrane regions, with poration represented by adding highly nonlinear, hysteretic pore models (e.g., Fig. 1a, b). In essence, the traditional Schwan model is extended by adding relatively simple representations of dynamic pores that change in size and conductance as D/ m varies with time and position. This is the basis of the present model (Smith 2011;Smith et al. 2013), which includes an unusually complex representation of dynamic pores.
Going forward, MD can provide a further, larger increase in realism using atomistic and molecular constructions, presently of the lipid bilayer regions of a cell membrane. MD examines the response of a small cell membrane patch, initially emphasizing lipids, and within which typically a single, dynamic pore typically appears with both fluctuations and overall expansion and contraction responding to D/ m . MD studies are increasingly enabled by increasing computing power (hardware and software), and the knowledge accumulating from a decade of EP simulations (Tieleman et al. 2003;Böckmann et al. 2008;Levine and Vernier 2010;Delemotte and Tarek 2012;Bennett et al. 2014). MD simulations are relevant to several major and inter-related predictions of the cell model, such as pore creation and lifetimes, pore density, and pore radius. We expect longer simulated times ðlsÞ, larger simulated areas ðlm 2 Þ; and more complex membranes (saturated and unsaturated glycerophospholipids, sphingolipids, and cholesterol) for electric field strengths corresponding to porating transmembrane voltages of $ 0:5 V or more. The incorporation of proteins, with molecular detail, will soon be possible. This offers the prospect of including lateral and transverse mechanical constraints imposed by membrane-associated proteins on the restructuring of the lipid bilayer to form a porated complex.
Finally, one specific connection between MD and continuum models, and associated experimental observations is the transport of small solutes through porated membranes. MD simulations can generate conductance and poration results for comparison to results from both continuum cell models and experiments. Quantitative comparisons will improve cell models and will result in a deeper understanding of the existing literature of measured values, and can stimulate and guide new experiments that will support and extend measurement methods and applications.

Post Pulse Landscape-Altering Mechanical and Chemical Interactions
We hypothesize that the pore energy landscape can change after the pulse, when D/ m % 0. More specifically, we suggest that the effective membrane tension may be significantly altered for cell membranes after field pulses, for durations of seconds to minutes.
An example of membrane tension change is porationinduced osmotic swelling. Large intracellular molecules cannot cross the porated membrane, and therefore create an osmotic gradient, and water diffuses into the cell to restore osmotic equilibrium. Cells respond to this osmotic imbalance in different ways. Some swell, resulting in a volume increase, and the extent and rate of swelling can be correlated to the extent and rate of ion concentration equilibration. The small ion flux that comes at the beginning of this process can be described in detail by MD, enabling a quantitative comparison between ion transport in the models and the swelling in cells.
Similarly, a candidate for poration-induced chemical alteration is reactive species (ROS) generation, in which ROS progressively alters the lipids of the membrane both during and after a pulse. Such chemically-induced changes can also alter the effective membrane tension, changing the landscape, and possibly also the effective pore diffusion coefficient. Relatively little work has been done on either mechanical or chemical interactions that could alter pore behavior after a pulse. Both are attractive general candidates for extending cell EP modeling capability.

Summary
We present illustrative responses of our present cell EP model, which consists of a closed curved cylindrical membrane that separates two aqueous media. The model is constructed from well known passive electrical properties (conductive and dielectric) of water-based media and also the fixed properties of a cell membrane. The active constituent of the model is a 1D model for membrane EP with dynamic pores. It is a recent version of a family of EP models that feature two rate equations, one for stochastic pore creation, the other for random pore destruction. These two equations introduce or remove pores from a pore energy landscape that changes locally with transmembrane voltage, D/ m ðtÞ. The various local models are connected by a meshed network that has 1D models for charge transport, forming a system model.
The model is a hypothesis which can be tested by subjecting the model to various stimuli, here two very different electric field pulses. Although these pulses are idealized trapezoids, digitized representations of actual experimental pulses can also be used. The model is complex, but nowhere near the complexity of a real cell membrane. However, the existing capabilities and the prospect for significant extensions together suggest that cell level modeling is an important complement to the many experimental studies of electroporation. A collaborative effort involving coordinated experiments and models could significantly accelerate basic studies, and also begin to set the stage for personalized models for tumor treatments, in which various imaging modalities aid model construction and biochemical information from circulating tumor cells may allow development of cell EP models that target cancer cells and tumors to improve treatments. Both EPbased drug delivery (ECT using bleomycin or cisplatin, or calcium EP) or drug-free, nonthermal ablation by IRE or nsPEF should benefit from this approach that begins to set the stage for a more focused biomedical engineering based on electroporation.

Appendix: Two Closely Related, Simple Passive Cell Models
While referring to all other underlying equations through citations that provide mathematical details and explanations, here we present two versions of a basic passive cell model equation, often referred to as the Schwan equation in the spherical case (Pauly and Schwan 1959;Kotnik and Miklavčič 2000;Foster 2000;Kotnik and Miklavčič 2006). The above papers contain considerable detail, with different approximations, and these include the transient responses of the two models. This simplest version is for the case that the passive cell membrane becomes fully charged. Then, fields and transmembrane voltages are no longer changing, and this direct current (dc) version is often presented. The two simplest versions are U m;sph ¼ 3 2 E e r cell cos h Spherical membrane Schwan equation ð1Þ U m;cyl ¼ 2E e r cell cos h Cylindrical membrane Schwan equation Here the notation is E e for the magnitude of the uniform field within the extracellular aqueous medium that fills the simulation box (entire system model) with no cell present. E e can differ significantly from the nominal field, E n , a potential source of error if not accounted for (Pliquett et al. 1996). This simple comparison shows that the 2D cylindrical version has the same form, with the induced transmembrane voltage differing only by a factor 4/3. The two fields can be essentially the same for some pulsing media and electrode composition, but for others the electrolyte field can be smaller, e.g., 90 % of the nominal strength (Canatella et al. 2001). This difference can be significant in the ''standard model'' given the D/ 2 m dependence in the creation rate equation. Our model assumes E e (t), the actual field within the electroporating medium.
This simplest case has zero resting potential (strictly resting transmembrane voltage, D/ m;rest ). These equations are consistent with field amplification (Fig. SI-2). We use the version that yields the transmembrane voltage at each membrane site, with the resting potential included through assignment of a nonequilibrium source of the resting potential (Läuger 1991;DeBruin and Krassowska 1999a).