Evaluation of the Current in the Cell Membrane for Numerical Simulations of Electroporation

Abstract

This study presents a rigorous numerical analysis of the current in the cell membrane, subjected to a uniform electric field, and its impact in the pore formation. The numerical model considers a single cell composed of uniform membrane and cytoplasm, in a suspension medium. The current in the cell membrane is calculated using two different approaches. The first uses a lumped parameters approach based on the geometry of the pore, while the second describes the flow of ions through the pore considering the interaction with the pore walls as an energy barrier. The formation and growth of the pores is solved using an asymptotic approximation of the Smoluchowski’s equation. The electrical potential induced in the cell membrane, which is coupled with the current in the membrane, is resolved in spherical coordinates using the finite difference method. The two approaches have qualitatively similar results but significant quantitative differences in the number and radii of pores. The ionic flow approach has resulted in the formation of fewer pores and reduced pore growth. Approximately 38,000 fewer pores are created, a 21% difference, and the largest pores are approximately 8nm smaller, a 24% difference. Thus, this approach results in a less conductive membrane and smaller electroporated area.

Index Terms
Biological cell; Electromagnetics; Finite Difference Methods; Numerical simulation

I. INTRODUCTION

A cell subjected to electric fields of sufficient strength but short duration exhibits the formation of pores, transient aqueous pathways, in its membrane, a phenomenon known as electroporation. This process can be reversible, with subsequent pore resealing, or irreversible with cell death [¹]. Both versions of the phenomenon have several applications in medicine and biology [²]. While there is vast evidence of the process, the small size and transient nature of pores have made direct observation impossible so far. Thus, theoretical models are necessary to further our understanding of the electroporation process. A variety of numerical models that describe pore formation and growth in single cells [³]-[⁶] and in cell arrangements [⁷], [⁸] has been implemented.

Electroporation models are directly dependent on calculation of the transmembrane potential, V_m, and the current in the cell membrane, which are dependent of each other. While it is possible to obtain analytical solutions for an unchanging membrane, the formation of conducting pores significantly changes the membrane properties. Thus, it is necessary to numerically solve the transmembrane potential simultaneously with the calculation of the pore density in order to accurately model the effects of electroporation.

The usual way of modeling this decrease in transmembrane potential, resulting from the increasing permeability, is by calculating an electroporation current through the pores [³], [⁴], [⁷], [⁹], [¹⁰] and treating intact areas of the membrane as uniform. This approach allows the modeling of the entire electroporation process in the scale of cells. Another possible way is molecular dynamics simulation, as seen in [¹¹], [¹²], that models in detail the movement and energy of molecules in the membrane. However, the latter is limited to simulating small areas of the membrane, in the order of nanometers, and very simple bilayers (composed of only a few types of lipids) that don’t correspond to real cell membranes [¹³] and the method has considerable computational cost even for relatively short simulations. This makes the method insufficient for simulations of the entire electroporation process.

It is possible to indirectly monitor the transport through the membrane by analyzing the intake of a marker substance, such as fluorescent dyes that bind to genetic material, membrane non-permeable dyes, magnetic nanoparticles, or cytotoxic substances [¹⁴]. This allows to discern whether electroporation has occurred, and gives a rough estimate of the size of the pores based on the dimensions of the transported molecules, but usually does not allow the observation of the process. Using more advanced methods such as rapid videoimaging systems, it is even possible to observe the temporal variation and spatial distribution of markers [¹⁵], [¹⁶]. Advanced optical techniques can also be used to detect changes in the membrane structure directly [¹⁷]. However, these methods have a low temporal resolution (in the order of seconds) that is not sufficient to model the fast stages of electroporation, particularly the formation and growth of the pores.

Other techniques are used to observe evidence of electroporation, such as using micropipettes to directly observe current flow through the membrane. This method allows for one of the best temporal resolutions (in the order of microseconds) and more direct observation, and has been used to obtain many electric and electroporation parameters used in theoretical and numerical models. However, the method is of difficult execution and, while it can be applied to small patches of the membrane, is still not able to observe individual micropores [¹⁴]. Thus, for most applications numerical simulation is the only available method for modeling the phenomenon with sufficient detail. In general, the most important parameters to model are the number of pores, their radii and distribution in the membrane.

The number of pores in the cell is especially important when using shorter pulses, with duration of few microseconds or less, as pore growth is a slower process and short pulses usually result in a large number of small pores with little variation in size. In some cases such as [⁷], the simulation is only run for a few nanoseconds and only the number of pores is considered. For longer pulses, a small number of very large pores is responsible for most of the transport of substances, but the fast creation of small pores is still responsible for a sudden increase in the cell membrane permeability in the early stages and the number of relatively large pores is important in general [³].

The size of the pore is a limiting factor for determining what substances can be introduced or removed from the cell [¹⁸]. Smaller molecules, like propidium iodide (2nm) that is used as an indicator of the creation of pores [⁶], can enter through relatively small pores. Conversely, bigger molecules like proteins and genetic code require much larger pores, though for those cases additional, non-passive transport mechanisms are also an important factor [²]. Thus, the pore radius plays a crucial role in the success of many applications. The number and size of the pores, in the form of the total electroporated area, are also an important information to determine cell survival for applications that use reversible electroporation. As more and especially bigger pores are formed, cell leakage may cause cell death [²], [⁴]. For some applications, such as microbial deactivation in food [¹⁹] and tumor treatment, achieving this threshold earlier is desired so that ohmic heating is minimized.

Accurate modeling of the current in the cell membrane is necessary for obtaining a good estimation of the size and number of pores, which is necessary in many applications. It is also important to choose the characteristics of the pulse, to build an appropriate circuit and to establish experimental protocols. However, due to the impossibility of measuring the current directly in this scale, the modeling must rely on rigorous physical modeling of the process.

The contribution of this work is an evaluation of the current in the cell membrane, using a numerical model based on an asymptotic approximation of the Smoluchowski’s equation [²⁰] and the Laplace’s equation for the electric potential. Both systems of equations are solved using the finite difference method in spherical coordinates. Two approaches for calculating the current through the pores, i.e. the lumped parameter and the ionic flow, are implemented and compared. The main results are reported in terms of the transmembrane potential, the electroporation current, the number of pores in the membrane and the radii of the pores. The two approaches result in large differences in the number of pores and their radii.

The paper is organized as follows. The models and methods section describes the biological context of electroporation, the geometry of the problem, the fundamental equations employed to model the electroporation phenomenon, the theoretical basis for the two approximations of current in the cell membrane and the implementation details of the numerical simulation. The results section presents a comparison between the two currents in terms of the temporal evolution of the number, size and distribution of the pores in the cell membrane. It also analyses the implication of those results for practical applications of electroporation process. Finally, the conclusion summarizes the findings drawn from this work. In ^{Appendix A}, further details are given on the derivation of the equations.

II. MODELS AND METHODS

A. History and theory of electroporation

Evidence of electroporation was first observed in 1958 [²¹] as a sudden and unexplainable increase in conductivity of some membranes. While many hypotheses were formulated in the following decades (weaver), it is now generally a consensus that the increase in membrane permeability is the result of the formation of transient aqueous pathways (pores) in the membrane [¹], [¹³].

The current theoretical models assume that initial pore formation happens due to energy fluctuations in the membrane. Without an applied electric field, pores are formed spontaneously, but due to their small size (smaller than 1nm) and instability (lasting less than 1ns) no significant transport occurs. The applied electric field lowers the energy threshold necessary for the creation of larger and more stable pores. A minimum transmembrane potential of approximately 1V is necessary in order to observe electroporation, with the fast creation of a large number of pores (the pore creation stage lasts a few nanoseconds). Once created, the pores continue to grow for as long as the external electric field is maintained. Thus, if the electric field is applied for a long enough time, uncontrolled pore growth can lead to membrane rupture and cell death. Conversely, if the electric pulse is sufficiently short, the pores can spontaneously reseal and the cell can survive. After the electric field is interrupted, the pores will shrink quickly (on the scale of microseconds), but will not fully reseal for a much longer time (in the scale of seconds or even minutes) [³].

Pores are initially created in a hydrophobic state, in which the polar molecules in the membrane restrict the passage of water and soluble substances. As the pore grows, reorientation of the molecules causes the pore to become hydrophilic (conducting. Most models ignore the hydrophobic stage, assuming that pores are created with a minimum size of 0.51 nm, the size at which pores become hydrophilic [¹], [³].

B. Geometry of the problem

The numerical model considers the problem of a single cell in a sparse suspension, such that there is no effect from nearby cells. The cell is assumed perfectly spherical; while that is an obvious simplification and is inaccurate for cells in tissues or clusters, it is a sufficient approximation for cells in suspensions, which are usually approximately spherical [²²], [²³].

The model is divided into three regions, as can be seen in Fig. 1: cytoplasm, membrane, and external medium. Each of these regions is assumed to have uniform properties and the membrane to have uniform thickness. While these are also simplifications, they are necessary due to the difficulty of obtaining accurate measurements for those properties in greater detail.

Fig. 1
Geometry of the problem (not to scale).

The dimensions and properties used in the model are given by Table I. The cytoplasm and external medium are described by their electrical conductivities σ_c and σ_o respectively; the electrical permittivity is ignored due to the mostly conductive nature of those media. The membrane, however, is described by a surface conductance σ_m, capacitance C_m, and resting potential V_rest. The use of a conductance and capacitance instead of conductivity and permittivity follows [³]. While the formation of pores alters the area of the membrane, and consequently C_m, this effect is insignificant for the time scales considered as seen in [²⁴].

Table I
Dimensions and properties of the cell as given by [³]

Parameter	Physical meaning	Value
R1	External radius	50*10−6 m
R2	Internal radius	49.995*10−6 m
σo	External medium conductivity	5 S/m
σc	Cytoplasm conductivity	0.455 S/m
Cm	Membrane capacitance	10−2 F/m2
σm	Membrane conductivity	2 S/m2
Vrest	Membrane resting potential	-0.08 V

The electric pulse is applied between two parallel plates, as seen in Fig. 1. It consists of a single square pulse with duration of 9μs and intensity of 40kV/m, following the methodology of [³]. This pulse is convenient as it allows the observation of the membrane charging process before the electroporation threshold is reached, the short pore formation period, and also the beginning of the much longer pore growth period. The plates are much bigger than the size of the cell and are treated as infinite planes, resulting in a uniform electric field in the ẑ direction.

The problem can be reduced to two dimensions, considering the rotational symmetry of the electric field around the ẑ axis. The problem can be further reduced to only a half-sphere, though not to a quarter sphere because the superposition of the applied electrical field with the rest potential of the membrane acts differently on each pole of the cell. This is discussed next.

C. Electroporation equations

The electroporation process is generally described by the Smoluchowski’s equation that describes pore formation and growth in terms of the pore density function n(r_p, t), as seen in [¹], [³]:

(1) $$\frac{\partial n\left(r_p,t\right)}{\partial t}+D\frac{\partial }{\partial r_p}\left(\frac{-\partial \phi \left(r_p,t\right)}{\partial r_p}\frac{1}{kT}-\frac{\partial n\left(r_p,t\right)}{\partial r_p}\right)=S\left(r_p\right)$$

where r_p is the pore radius, D is the pore diffusion constant, k is the Boltzmann constant and T is the absolute temperature. S(r_p) is the source term representing the creation and destruction of pores:

(2) $$S\left(r_p\right)=\frac{v_c}{kT}\frac{\partial \phi \left(r_p,t\right)}{\partial r_p}{e}^{-\phi \left(r_p\right)/kT}dr$$

where v_c is the fluctuation rate. φ(r_p,t) is the function that describes pore energy, given by:

(3) $$\phi \left(r_p,t\right)=U\left(r_p\right)-\pi a_p{r}_p^2{V}_m^2\left(t\right)$$

where U(r_p) is the energy value with no applied potential, a_p is a property of the membrane and V_m is the transmembrane potential.

In this work, an asymptotic approximation of the Smoluchowski’s equation is used [²⁰]. Following [³], equations (1)–(3) are rearranged into a form that directly outputs the number of pores N and the radius of the pores r_j, which are convenient for implementing the simulation directly using finite differences.

The rate of pore formation in the cell membrane, based on the equilibrium pore number N_eq, is given by [³]:

(4) $$\frac{dN}{dt}=\alpha e^{{\left(V_m/V_{ep}\right)}^2}\left(1-\frac{N}{N_{eq}\left(V_m\right)}\right)$$

where N is the number of pores, α is the creation rate coefficient, V_m is the transmembrane potential, V_ep is the characteristic voltage of electroporation. The equilibrium pore number N_eq is given by:

(5) $$N_{eq}\left(V_m\right)=N_0{e}^{q{\left(V_m/V_{ep}\right)}^2}$$

where N₀ is the equilibrium pore density for V_m = 0 and q is a constant for the pore creation rate.

The rate of growth for the radius of each of the K individual pores already created is given by:

(6) $$\frac{d{r}_j}{dt}=U\left(r_j,V_m,S_{eff}\right),j=1,2,\mathrm{\dots },K$$

in which the source term U is defined by

(7) $$U\left(r,V_m,S_{eff}\right)=\frac{D}{kT}\left[\frac{V_m^2{F}_{max}}{1+r_h/\left(r+r_t\right)}+4\beta {\left[\frac{r}{r_x}\right]}^4\frac{1}{r}\right]+\frac{D}{kT}\left[-2\pi \gamma +2\pi S_{eff}r\right]$$

with D being the pore radius diffusion coefficient, k the Boltzmann constant, T the absolute temperature, F_max the maximum electrical force for V_m = 1V, r_h and r_t constants for the advection velocity, β the energy associated with steric repulsion, γ the edge energy, r_x is the minimum pore radius for conducting pores and S_eff is the effective membrane tension given by:

(8) $$S_{eff}\left(A_p\right)=2S\text{'}-\frac{2S\text{'}-S_0}{1-A_p/A}$$

where A_p is the combined electroporated area (adding all the pores) and A is the total membrane area, S’ is the tension of the hydrocarbon-water interface and S₀ is the tension of the membrane without pores. The values of all the constants and parameters are given by Table II.

Table II
Parameters of the asymptotic smoluchowski's equation as given by [³]

Parameter	Physical meaning	Value
k	Boltzman constant	1.38*10−23J/K
T	Absolute temperature	300 K
N0	Equilibrium pore density	1.5*109 1/m2
D	Pore radius diffusion coefficient	5*10−14m2/s
q	Constant for pore creation rate	2.4606
Vep	Characteristic voltage of electroporation	0.258V
rx	Minimum pore radius	0.51*10−9 m
rh	Constant for advection velocity	0.97*10−9 m
rt	Constant for advection velocity	0.31*10−9 m
α	Creation rate coefficient	1*109 1/(m2s)
β	Steric repulsion energy	1.4*10−19 J
γ	Edge energy	1.8*10−11 J/m
Fmax	Max electric force for Vm = 1 V	0.7*10−9 N/V2
S’	Tension of hydrocarbon-water interface	2*10−2 J/m2
S0	Tension of membrane without pores	1*10−6 J/m2

D. Calculation of the transmembrane potential

The transmembrane potential, V_m, is the difference between the potential on both sides of the membrane:

(9) $$V_m\left(\theta \right)={V_c\left(r,\theta \right)|}_{r=R_2}-{V_o\left(r,\theta \right)|}_{r=R_1}$$

in which V_c is the potential in the cytoplasm and V_o is the potential in the external medium.

The potential in the external medium and inside the cell is obtained by solving the Laplace’s equation at each time instant, according to:

(10) $${{\nabla }^2{V}_o=0|}_{{\Omega }_o}$$

(11) $${{\nabla }^2{V}_c=0|}_{{\Omega }_c}$$

The problem can be reduced to a half-sphere due to the symmetry of the problem, with a simple Neumann boundary condition in the symmetry axis Γ₂. The problem has also a Dirichlet boundary condition in the limits of the domain Γ₁ (set as a semi-spherical shell with three times the pore radius). This boundary condition indicates that the electric field far from the cell converges to the electrical field between the parallel plates with value of E_0(t), as given by:

(12) $$V_o\left(t,r,\theta \right)=-E_0\left(t\right)rcos\left(\theta \right),r\to \infty \text{in} {\Gamma }_1$$

The interface conditions are given by the continuity of potential and current in the interface. The potential continuity is given by:

(13) $${V_o|}_{r=R_1}={V_{mem}|}_{r=R_1}$$

(14) $${V_c|}_{r=R_2}={V_{mem}|}_{r=R_2}$$

and the current by:

(15) $$-\widehat{n}\cdot \left({\sigma }_c\nabla V_c\right)=-\widehat{n}\cdot \left({\sigma }_o\nabla V_o\right)=C_m\frac{\partial V_{mem}}{\partial t}+{\sigma }_m\left(V_m-V_{rest}\right)+I_p$$

in which I_p is the electroporation current, defined as the sum of the currents through each individual pore:

(16) $$I_p\left(\theta \right)={\displaystyle {\sum }_{j=1}^{K(\theta )}\frac{i_p\left(r_{j,V_m}\right)}{△A}}$$

with i_p being the current for an individual pore, ΔA is the area for that section of the membrane and j being the index of individual pores within the group K of all pores in that division.

Fig. 3 illustrates the continuity of current in the problem. The presence of the I_p term in (15) causes V_m to depend strongly on the current through the pores, and consequently on the pore formation process. As a result, V_m must be calculated numerically, coupled with the pore formation simulation, and an approximation to calculate i_p must be used.

Fig. 2
Discretization of the problem (not to scale).

Fig. 3
Current continuity and cell properties (not to scale).

The finite difference method is chosen to solve the system of equations (10)–(11) due to its simplicity and the ability to easily model the domain. A spherical coordinate system is adopted to allow for more accurate representation of the spherical cell and the Laplace’s equation is discretized using central differences.

The domain is discretized in 60 divisions on the r direction, 20 in the cytoplasm and 40 in the external medium, and 64 divisions in the θ direction [³], as represented in Fig. 2. The membrane, however, is treated differently, as it is much thinner than the r division (5nm and 2.5μm respectively), and using a sufficiently fine discretization in the entire domain would greatly increase the computational cost. The membrane is treated as a single division, with length equal to membrane thickness h, between segments 20 and 21.

A system of equations is created using the coefficients of the discretized Laplace’s equation, which is solved at each instant (with a time step of 1.5ns) for updated values of I_p. Equations (4) and (6) are discretized using forward differences and solved at each instant for the newly updated value of V_m, resulting in the updated number of pores and their respective radii. That, in turn, allows the calculation of I_p using the approximations described in sections II.D and II.E. This is detailed in ^{Appendix A}.

E. Lumped parameters approximation for the current through a pore

This approach is a simplification first proposed in [²⁵] for the modeling of large pores. Due to the large radius of the pores, the non-ohmic interactions between the pore walls and the flowing charges were ignored. However, in later works such as [³] this assumption was not explicitly mentioned and later works by different groups, such as [⁷], used this simplification uncritically for the modeling of populations consisting entirely of small pores.

In this approximation, the current through each pore is calculated by the series association of two resistances:

(17) $$i_p=\frac{V_m}{R_p+R_i}$$

where R_p is the resistance of a cylinder representing the pore. Using Ohm’s law:

(18) $$R_p=\frac{h}{\pi \sigma r_p^2}$$

and R_i is a non-linear input resistance that models the non-uniform current in the interface between a large conducting medium and a relatively narrow conducting cylinder embedded in a non-conducting surface, following the methodology found in [²⁶]:

(19) $$R_i=\frac{1}{2\sigma r_p}$$

where h is the membrane thickness, r_p is the pore radius and σ = 2 S/m is the conductivity of the solution filling the pore.

The main advantage of this approach is its simplicity. It uses a simple geometric model that does not require extensive knowledge of the pore formation mechanism. However, it has the disadvantage of not representing more complex effects, like interaction between charges and the pore walls.

F. Ionic flow approximation for the current through a pore

In this approximation, the current is calculated as the flow of ions through the electropores, as described in [²⁷]. This approach uses a variation of the Nernst-Planck equation, that describes active (potential-driven) and passive (osmosis and diffusion) transport of substances through the membrane. By modeling the transport of ions across the pore caused by the potential difference, it derives a relation between voltage and current for each individual pore.

This approach was used in [¹⁰], but the model failed to consider the spatial distribution of the potential in the cell. The same approach was used for a more complete simulation in [⁴], [²²], using a commercial finite elements software. The formulation used in this paper follows that presented in [²⁸].

The current on a pore of radius r_p is given by:

(20) $$i_p=V_m\frac{\pi \sigma \frac{r_p^2}{h}\left(1-e^{-V_m^x}\right)}{e^{-V_m^x}\frac{w_0{e}^{\left(w_0+n{V}_m^x\right)}+n{V}_m^x}{w_0+n{V}_m^x}-\frac{w_0{e}^{\left(w_0-n{V}_m^x\right)}-n{V}_m^x}{w_0-n{V}_m^x}}$$

where w₀ is the energy barrier inside the pore that models the interaction between the pore walls and the ions, q_e is the electron charge, n is the relative size of the pore entrance region and V_m^x is a normalized potential given by:

(21) $$V_m^x=V_m\frac{q_e}{kT}$$

The parameters are given in Table II and Table III. This approach has the advantage of being more complete, with the inclusion of an energy barrier and a more detailed physical model for the conduction through the pore. Its main disadvantage is the increased complexity.

Table III
Parameters for ion flow calculation as given by [²⁸]

Parameter	Physical meaning	Value
qe	Electron charge	1.602*10−19 C
w0	Energy barrier	8.411*10−19 J
n	Relative size of the entrance	0.15

G. Numerical implementation

The simulation was run for the duration of a single electrical pulse (9μs), simulating the pore formation and the start of the pore growth period but not the pore resealing – which happens on a much slower timescale, of milliseconds or even seconds, making simulation less practical. The pores created at each time instant and in each angular division of the membrane are treated as a single group for purposes of radius increase calculations.

III. RESULTS

A. General characteristics of the electroporation process

Electroporation happens in three distinct stages, as can be seen in Fig. 4 and Fig. 5: membrane charging, pore creation and pore growth. In the first stage, the membrane is intact and behaves mostly like a dielectric with very small losses, being charged by the applied pulse until the transmembrane potential reaches the electroporation threshold value of 1V. The first pore is created at 49.45ns, marking the start of the very fast pore formation stage in which a large number of new pores are created very quickly. All pores are created with a minimum radius of 0.51 nm, which is the size at which pores become conductive. As current starts being conducted through the pores, the potential falls and pore creation slows down. After 1.75μs the potential falls below the threshold and pore creation stops completely, starting the pore growth stage in which no new pores are created but the existing pores continue to increase in size for the remaining duration of the simulation.

Fig. 4
The three stages of electroporation, compared with the number of pores and the transmembrane potential in the division θ = 0°. Data obtained using the lumped parameters method.

Fig. 5
The three stages of electroporation, compared with the number of pores and the maximum pore radius in the division θ = 0°. Data obtained using the lumped parameters method.

Spatially, the electroporation process differs as one moves along the membrane, as can be seen in Fig. 6. Due to the alignment with the direction of the applied pulse, the potential increases very quickly and reaches the threshold early in the regions around the poles θ=0°and θ=180°. This causes the creation of a very large number of pores in these regions, as can be seen in Fig. 6. The largest pores appear on θ=129.9° and its vicinity, and on the other half of the cell the pore radius peaks on θ=53.44°. These regions have fewer pores than the poles, which causes the potential to fall more slowly during the pore creation stage and stabilize at a higher value during the pore growth stage and consequentially causes faster pore growth. In the central region of the membrane, from approximately θ=60° to approximately θ=120°, the potential never increases above the threshold necessary for pore creation and as a result no pores are present in this region.

Fig. 6
Spatial distribution of the number of pores and maximum radius of pores after 9μs. Data obtained using the lumped parameters method.

B. Transmembrane potential

The simulation for both approximations resulted in qualitatively similar but quantitatively different transmenbrane potentials and electroporation currents. The temporal variation of those quantities can be seen on Fig. 7 and Fig. 8., taken for a single angular division at θ=0° and truncated at 2μs for easier visualization of the relevant phenomena. The spatial variation of the transmembrane potential can be seen on Fig. 9, taken at the final time instant simulated.

Fig. 7
Comparison of V_m with different methods for the division θ=0°.

Fig. 8
Comparison of I_p with different methods for the division θ=0°.

Fig. 9
Transmembrane potential across the membrane after 9μs.

Fig. 7 shows that the transmembrane potential reaches approximately the same peak value. The lumped parameters approximation results in a bigger potential by less than a 0.5% difference, but the ionic flow approximation results in a noticeably faster drop. This is reflected by the similarly earlier fall in the current in Fig. 8, and impacts the rate of pore creation and eventual end of the pore creation stage. For the θ=0° division, the potential stabilizes at approximately 0.474V for the lumped parameters approximation and 0.468V (at 6.622μs) for the ionic flow, a difference of approximately 1.1% that, as will be seen, impacts the pore growth stage. As can be seen on Fig. 9, the equilibrium value varies significantly for different divisions but the lumped parameter approach results in larger values across the entire membrane.

As can be easily seen on Fig. 8, the ionic flow approximation results in a much larger current peak (77.91kA/m²) than the lumped parameters approximation (72.11 kA/m²). In the earlier stage the difference between the two is proportionally much larger, up to a 23.1% difference, decreasing after approximately 600ns. The ionic flow current also reaches its peak and starts decreasing earlier than the lumped parameters current. For both cases, after a brief oscillation the current mostly stabilizes with less than a 0.1% difference between the two approximations after 4.143μs.

C. Number of pores

The number of pores created over time for the θ=0° division can be seen on Fig. 10, while the number of pores created for the entire simulation can be seen on Fig. 11.

Fig. 10
Number of pores of all sizes for the division θ=0°.

Fig. 11
Number of pores across the membrane after 9μs.

It can be easily seen that the lumped parameters approximation results in the creation of a larger number of pores. For both approximations, the fist pore appears at the same time instant, but the rate of pore creation becomes significantly different after a few hundred pores have been created. When the pore creation stops, the lumped parameters approach has resulted in the creation of 2,034 more pores than the ionic flow approach in the θ=0° division alone. The total difference is of 37,904 pores in the simulated half of the cell.

D. Pore radius

The pore growth can be seen on Fig. 12 and Fig. 13. Fig. 12 shows only the maximum pore radius for the division θ=0° and is truncated at 5μs for better visualization of the early stages. Fig. 13 compares the pore growth for two divisions, θ=180° and θ=45°, for the entire 9μs simulated. The pore radius distribution across the membrane can be seen on Fig. 14, that plots the maximum pore radius for each division after 9μs.

Fig. 12
Radius of the biggest pore for the division θ=0°.

Fig. 13
Radius of the biggest pore for the divisions θ=0° and θ=45°.

Fig. 14
Radius of the biggest pore across the membrane after 9μs.

From those results, it can be seen that the lumped parameters approach results in not only more pores, but also faster pore growth. This can be easily attributed to a higher transmembrane potential during two different stages: in the pore creation stage the decrease in the potential happens later for the lumped parameters approach, leading to a period where the voltage is higher right after the formation of many pores. During the pore growth stage the equilibrium potential is slightly higher, causing faster pore growth over an extended period of time. Likewise, pore growth is faster for pores near θ=129.9° and θ=53.44°, with an increasing difference between regions as the simulation progresses.

It should be noted that for the division θ=0° the pore radius has a local maximum during the pore creation stage, while pores in division θ=45° grow monotonically. In the former, potential increases very quickly resulting in fast pore growth, but the sudden drop in the potential due to the creation of a large number of pores results in a small shrinkage.

E. Analysis of the differences between the approaches

It can be deduced that the current obtained through the ionic flow approximation is significantly smaller than that obtained through the lumped parameters approximation for the early stages of the electroporation. This causes a slower drop in the transmembrane potential, with the membrane subjected to a higher potential for a longer time. This in turn results in higher pore creation and pore growth rates during the short pore formation stage.

Eventually, the system reaches a state of near equilibrium, as the system displays a “negative feedback” behavior: high potential causes the formation and growth of pores, which in turn increase the current causing the potential to decrease. As the potential falls below the pore creation threshold new pores stop being created, ending the pore creation stage. While the pores continue to grow due to the small but still significant potential in this equilibrium, it is a much slower process compared to the sudden creation of thousands of pores and fast increase in their radii and its impact on the electroporation current is not observable in the time scales of this study – after several milliseconds it will result in pore shrinkage, with complete pore resealing taking even longer, as seen for example in the similar simulations in [³]. However, it can be observed that the slightly higher equilibrium potential for the lumped parameters approximation results in an increasing difference in pore radius during the early parts of the pore growth stage.

F. Impact on applications

As mentioned, many applications are strongly dependent on the size and number of pores in the membrane. Those using shorter pulses, on the order of nanoseconds, are generally more concerned with the number of pores rather than with the pore radius. On the other hand, applications that require the transport of substances with large molecules have pore radius as a primary concern. For applications where the conductivity is important, as well as applications where cell death or viability are the main concern, both quantities are relevant.

The different approximations for the current in the cell membrane result in significant differences in those quantities, with an obvious impact in those applications. The ionic flow approach predicts the creation of more and larger pores, and consequently a more conductive membrane with a larger electroporated area. This suggests that electric pulses with lower intensity and/or shorter duration should be used to obtain the same effect, compared to the pulses predicted by the lumped parameters approach.

IV. CONCLUSION

The different approximations for the calculation of current through electropores have clear and significant quantitative effects on the simulation of the electroporation process. While they are qualitatively similar and display broadly the same overall characteristics such as the distinct stages of the process and the same spatial distribution, any application that relies on the accurate calculation of electroporation characteristics must pay special attention to accurately modeling the electroporation current.

The lumped parameters approximation is a simpler theoretical model, that relies exclusively on geometric and average electrical properties. It estimates a larger current, especially for small (recently formed) pores, which results in less pores for all divisions of the membrane, smaller pore radius through the entire process and slower pore growth resulting in increasing differences over time.

The ionic flow approximation uses a more complex and comprehensive physical model to calculate the current through electropores. This results in a lower estimate for the current and consequent formation of more pores and faster pore growth.

The tests were run using Matlab, in an Intel i5-3210M with 8GB of RAM using Kubuntu Linux. The tests using the lumped parameters approach took approximately 45 seconds and those using the ionic flow approach took approximately 53 seconds. Thus, the lumped parameters approach is significantly faster, as expected, though the computational cost of the ionic flow approach is not prohibitive.

Both models are approximations derived from physical models, due to the impossibility of directly observing the current through individual pores. However, the assumptions, such as pores being relatively large and the energy barrier being negligible, made by the lumped parameters approximation are less general. Such assumptions might not be true for the starting stages of electroporation or for very short pulses, as pore growth is significantly slower than pore formation. The ionic flow approximation does not make those assumptions, and thus it is recommended despite its higher complexity and computational cost.

ACKNOWLEDGMENT

This work has been supported by the Brazilian agencies Coordination for the Improvement of Higher Education Personnel (CAPES) and National Council for Scientific and Technological Development (CNPq).

REFERENCES

• ^[1]
  T. Kotnik, P. Kramar, G. Pucihar et al., “Cell membrane electroporation- part 1: The phenomenon,” Electrical Insulation Magazine, IEEE, vol. 28, no. 5, pp. 14–23, September 2012.
• ^[2]
  S. Haberl, D. Miklavčič, G. Sersa et al., “Cell membrane electroporation-part 2: the applications,” Electrical Insulation Magazine, IEEE, vol. 29, no. 1, pp. 29–37, January 2013.
• ^[3]
  W. Krassowska and P. D. Filev, “Modeling electroporation in a single cell,” Biophysical Journal, vol. 92, no. 2, pp. 404 – 417, 2007.
• ^[4]
  P. Lamberti, S. Romeo, A. Sannino et al., “The role of pulse repetition rate in nspef-induced electroporation: A biological and numerical investigation,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 9, pp. 2234–2243, Sept 2015.
• ^[5]
  W. Ying and C. S. Henriquez, “Hybrid finite element method for describing the electrical response of biological cells to applied fields,” IEEE Transactions on Biomedical Engineering, vol. 54, no. 4, pp. 611–620, April 2007.
• ^[6]
  J. A. Ramírez, W. P. Figueiredo, J. F. C. Vale et al., “Investigation of the electroporation effect in a single cell,” COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, vol. 32, no. 5, pp. 1692–1706, 2013.
• ^[7]
  J. Dermol–Černe and D. Miklavčič, “From cell to tissue properties - modeling skin electroporation with pore and local transport region formation,” IEEE Transactions on Biomedical Engineering, vol. 65, no. 2, pp. 458–468, 2018.
• ^[8]
  R. Susil, D. Šemrov, and D. Miklavčič, “Electric field-induced transmembrane potential depends on cell density and organization,” Electro-and magnetobiology, vol. 17, no. 3, pp. 391–399, 1998.
• ^[9]
  G. Pucihar, D. Miklavčič, and T. Kotnik, “A time-dependent numerical model of transmembrane voltage inducement and electroporation of irregularly shaped cells,” IEEE Transactions on Biomedical Engineering, vol. 56, no. 5, pp. 1491–1501, May 2009.
• ^[10]
  R. P. Joshi, Q. Hu, R. Aly et al., “Self-consistent simulations of electroporation dynamics in biological cells subjected to ultrashort electrical pulses,” Phys. Rev. E, vol. 64, p. 011913, Jun 2001.
• ^[11]
  Q. Hu, Z. Zhang, H. Qiu et al., “Physics of nanoporation and water entry driven by a high-intensity, ultrashort electrical pulse in the presence of membrane hydrophobic interactions,” Phys. Rev. E, vol. 87, p. 032704, Mar 2013.
• ^[12]
  V. Sridhara and R. Joshi, “Evaluations of a mechanistic hypothesis for the influence of extracellular ions on electroporation due to high-intensity, nanosecond pulsing,” Biochimica et Biophysica Acta (BBA) – Biomembranes, vol. 1838, no. 7, pp. 1793 – 1800, 2014.
• ^[13]
  A. A. Gurtovenko, J. Anwar, and I. Vattulainen, “Defect-Mediated Trafficking across Cell Membranes: Insights from in Silico Modeling,” Chem Review, vol. 110, pp. 6077–6103, 2010.
• ^[14]
  T. B. Napotnik and D. Miklavčič, “In vitro electroporation detection methods – An overview,” Biolectrochemistry, vol. 120, pp 166 – 182, 2018
• ^[15]
  L. Rems and D. Miklavčič, “Tutorial: Electroporation of cells in complex materials and tissue,” Journal of Applied Physics, vol. 119, no. 20, 2016.
• ^[16]
  M. A. Chiapperino, L. Mescia, P. Bia et al., “Experimental and Numerical Study of Electroporation Induced by Long Monopolar and Short Bipolar Pulses on Realistic 3D Irregularly Shaped Cells,” IEEE Transactions on Biomedical Engineering, vol. 67, no. 10, pp. 2781 – 2788, Oct 2020.
• ^[17]
  E. K. Moen, B. L. Ibey, H. T. Beier et al., “Quantifying pulsed electric field-induced membrane nanoporation in single cells,” Biochimica et Biophysica Acta (BBA) – Biomembranes, vol. 1858, no. 11, pp. 2795 – 2803, 2016.
• ^[18]
  T. Kotnik, W. Frey, M. Sack et al., “Electroporation-based applications in biotechnology,” Trends in Biotechnology, vol. 33, no. 8, pp. 480 – 488, 2015.
• ^[19]
  E. J. Araújo, I. J. S. Lopes, and J. A. Ramírez, “Computational modeling approach for the optimisation of a pulsed electric field system for liquid foods,” IET Science, Measurement & Technology, vol. 13, pp. 337–345(8), May 2019.
• ^[20]
  J. C. Neu and W. Krassowska, “Asymptotic model of electroporation,” Phys. Rev. E, vol. 59, pp. 3471–3482, Mar 1999.
• ^[21]
  R Stämpfli, “Reversible electrical breakdown of the excitable membrane of a ranvier node,” Anais da Academia. Brasileira de Ciências, vol. 30, pp. 57 – 63, 1958.
• ^[22]
  W. Krassowska and J. Neu, “Response of a single cell to an external electric field,” Biophysical Journal, vol. 66, no. 6, pp. 1768 – 1776, 1994.
• ^[23]
  T. Kotnik and D. Miklavčič, “Analytical description of transmembrane voltage induce by electric fields on spheroidal cells,” Biophisical Journal, vol. 79, no. 2, pp. 670–679, 2000.
• ^[24]
  A. Garner, N. Chen, J. Yang, J. Kolb et al., “Time domain dielectric spectroscopy measurements of hl-60 cell suspensions after microsecond and nanosecond electrical pulses,” Plasma Science, IEEE Transactions on, vol. 32, no. 5, pp. 2073–2084, Oct 2004.
• ^[25]
  K.C. Smith, J.C. Neu, W. Krassowska, “Model of Creation and Evolution of Stable Electropores for DNA Delivery”, Biophysical Journal, vol 86, pp. 2813-2826, 2004.
• ^[26]
  J. Newman, “Resistance for flow of current to a disk,” Journal of the Electrochemical Society, vol. 113, no. 5, pp. 501–502, 1966.
• ^[27]
  A. Barnett, “The current-voltage relation of an aqueous pore in a lipid bilayer membrane,” Biochimica et Biophysica Acta (BBA) - Biomembranes, vol. 1025, no. 1, pp. 10 – 14, 1990.
• ^[28]
  P. Lamberti, V. Tucci, S. Romeo et al., “nspef-induced effects on cell membranes: use of electrophysical model to optimize experimental design,” IEEE Transactions on Dielectrics and Electrical Insulation, vol. 20, no. 4, pp. 1231–1238, August 2013.
• ^[29]
  O. Momoh, M. Sadiku, and C. Akujuobi, “Solution of axisymmetric potential problem in spherical coordinates using exodus method,” Session 4P3 Novel Mathematical Methods in Electromagnetics 2, p. 673, 2010.

APPENDIX

As described in II.B, the electroporation process is described by (4)–(6), which are derived from an asymptotic approximation of Smoluchowski’s equation.

In order to simulate the process, those equations are discretized using forward finite differences, with a discrete time step Δt.

The discretization of (4) results in:

(22) $$N_k^{n+1}=\left(\frac{\left(1-u_n\left(V_m\right)\right)}{1+u_n\left(V_m\right)}{N}_k^n+\frac{\alpha \Delta t{e}^{{\left(Vm/Vep\right)}^2}}{1+u_n\left(V_m\right)}\right)$$

where the superscript n is used to indicate the current time step. The auxiliary term u_n is given by:

(23) $$u_n\left(V_m\right)=\frac{\Delta t\alpha }{2N_0}{e}^{\left(1-q\right){\left(V_m/V_{ep}\right)}^2}$$

while the discretization of (6) results in:

(24) $$r_j^{n+1}=r_j^n+\frac{D\Delta t}{kT}\frac{\left(r_j^n+r_t\right)}{r_j^n+r_t+r_n}{F}_{max}{V}_m^2+\frac{D\Delta t}{kT}\left[4\beta {\left(r_x\right)}^4\frac{1}{{\left(r_j^n\right)}^5}+2\pi {\sigma }_{eff}{r}_j^n\right]-4\beta {\left(r_x\right)}^{4}2\pi \gamma$$

The potential is obtained by solving the Laplace’s equation (10). However, due to the geometry of the problem, spherical coordinates are used, and the expanded Laplace’s equation is given in this coordinate system by:

(25) $${\nabla }^{2}V\left(r,\theta \right)=\frac{{\partial }^{2}V}{\partial r^2}+\frac{2}{r}\frac{\partial V}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}V}{\partial {\theta }^2}+\frac{cos\theta }{r^{2}sin\theta }\frac{\partial V}{\partial \theta }=0$$

By discretizing this equation using central finite differences, with divisions Δθ and Δr and indexes i and j respectively, we obtain:

(26) $$\begin{array}{c}V\left(i,j\right)=\frac{1}{2\left(1+\frac{1}{{\left(i\Delta \theta \right)}^2}\right)}\times [\left(1+\frac{1}{i}\right)V\left(i+1,j\right)+\left(1-\frac{1}{i}\right)V\left(i-1,j\right)\\ +\left(\frac{1}{{\left(i\Delta \theta \right)}^2}+\frac{cosj\Delta \theta }{2i^2\Delta \theta sinj\Delta \theta }\right)V\left(i,j+1\right)+\left(\frac{1}{{\left(i\Delta \theta \right)}^2}-\frac{cosj\Delta \theta }{2i^2\Delta \theta sinj\Delta \theta }\right)V\left(i,j-1\right)]\end{array}$$

To prevent singularities, special treatment must be given to equations at the lines corresponding to θ=0° and θ=180°, as seen in [²⁹]. This is done by applying Neuman boundary conditions at the boundary Γ₂ and replacing V(i, j + 1) = V(i, j − 1) in the equations above. This results in a slightly different equation given by:

(27) $$V\left(i,0\right)=\frac{1}{2\left(1+\frac{1}{{\left(i\Delta \theta \right)}^2}\right)}\times \left[\left(1+\frac{1}{i}\right)V\left(i+1,0\right)+\left(1-\frac{1}{i}\right)V\left(i-1,0\right)\right]$$

(28) $$V\left(i,j_{max}\right)=\frac{1}{2\left(1+\frac{1}{{\left(i\Delta \theta \right)}^2}\right)}\times \left[\left(1+\frac{1}{i}\right)V\left(i+1,j_{max}\right)+\left(1-\frac{1}{i}\right)V\left(i-1,j_{max}\right)\right]$$

More importantly, the interface condition in the membrane also requires special treatment, as the potential in different media can not be used to calculate the potential on either side of the interface. Instead, the first term of the Taylor expansion for the current in the interface as a function of the potentials and properties of the same media is replaced in the interface condition equation (12). This results in a different set of equations for the points on the interface, as seen in:

(29) $$\begin{array}{c}\left[2\left(1+\frac{1}{{\left(i_{R_2}\Delta \theta \right)}^2}\right)-\frac{2C_m\Delta r}{{\sigma }_o\Delta t}\left(-1+\frac{1}{i_{R_2}}\right)\right]V\left(i_{R_2},j\right)\\ \frac{+2}{{\sigma }_o}\left(1+\frac{1}{{\left(i_{R_2}\Delta \theta \right)}^2}\right)\left(\frac{C_m}{\Delta t}{V}_m\left(j\right)+{\sigma }_m{V}_m\left(j\right)+I_{ep}\right)=2V\left(i_{R_2},j+1\right)-\frac{2C_m\Delta r}{{\sigma }_o\Delta t}\left(1-\frac{1}{i_{R_2}}\right)V\left(i_{R_2},j-1\right)\\ \frac{+1}{{\left(i_{R_2}\Delta \theta \right)}^2}\left(1-\frac{\Delta \theta cot\left(j\Delta \theta \right)}{2}\right)V\left(i_{R_2}-1,j\right)+\frac{1}{{\left(i_{R_2}\Delta \theta \right)}^2}\left(1+\frac{\Delta \theta cot\left(j\Delta \theta \right)}{2}\right)V\left(i_{R_2}+1,j\right)\end{array}$$

(30) $$\begin{array}{c}\left[-2\left(1+\frac{1}{{\left(i_{R_1}\Delta \theta \right)}^2}\right)+\frac{2C_m\Delta r}{{\sigma }_c\Delta t}\left(1+\frac{1}{i_{R_1}}\right)\right]V\left(i_{R_1},j\right)\\ \frac{+2}{{\sigma }_c}\left(1+\frac{1}{{\left(i_{R_1}\Delta \theta \right)}^2}\right)\left(\frac{C_m}{\Delta t}{V}_m\left(j\right)+{\sigma }_m{V}_m\left(j\right)+I_{ep}\right)=2V\left(i_{R_1},j-1\right)+\frac{2C_m\Delta r}{{\sigma }_c\Delta t}\left(1+\frac{1}{i_{R_1}}\right)V\left(i_{R_1},j+1\right)\\ \frac{+1}{{\left(i_{R_1}\Delta \theta \right)}^2}\left(1-\frac{\Delta \theta cot\left(j\Delta \theta \right)}{2}\right)V\left(i_{R_1}-1,j\right)+\frac{1}{{\left(i_{R_1}\Delta \theta \right)}^2}\left(1+\frac{\Delta \theta cot\left(j\Delta \theta \right)}{2}\right)V\left(i_{R_1}+1,j\right)\end{array}$$

Equations (27), (28), (29) and (30), taken for every point in the domain, create a system of equations that relates the potential at each point to the one in the nearby points. In the boundaries of the domain, a fixed value is imposed representing the current value of the applied pulse. This system is resolved at every time instant for the updated pulse value and also the updated electroporation current. The resulting potential is used to solve the electroporation equations, with the corresponding electroporation current to be used in the next iteration.

Publication Dates

• Publication in this collection
  09 Mar 2022
• Date of issue
  Mar 2022

History

• Received
  22 Dec 2020
• Reviewed
  11 Jan 2021
• Accepted
  07 Oct 2021