Article Information

• PDF (185 kb)
• Full Text with Large Figures

PubMed

• Articles by Katherine A. DeBruin
• Articles by Wanda Krassowska

Related Articles

• Modeling Electroporation in a Single Cell

  Modeling Electroporation in a Single Cell Biophysical Journal, Volume 92, Issue 2, 15 January 2007, Pages 404-417 Wanda Krassowska and Petar D. Filev Abstract Electroporation uses electric pulses to promote delivery of DNA and drugs into cells. This study presents a model of electroporation in a spherical cell exposed to an electric field. The model determines transmembrane potential, number of pores, and distribution of pore radii as functions of time and position on the cell surface. For a 1-ms, 40kV/m pulse, electroporation consists of three stages: charging of the cell membrane (0–0.51s), creation of pores (0.51–1.43s), and evolution of pore radii (1.43s to 1ms). This pulse creates ∼341,000 pores, of which 97.8% are small (≈1nm radius) and 2.2% are large. The average radius of large pores is 22.8±18.7nm, although some pores grow to 419nm. The highest pore density occurs on the depolarized and hyperpolarized poles but the largest pores are on the border of the electroporated regions of the cell. Despite their much smaller number, large pores comprise 95.3% of the total pore area and contribute 66% to the increased cell conductance. For stronger pulses, pore area and cell conductance increase, but these increases are due to the creation of small pores; the number and size of large pores do not increase. Abstract | Full Text | PDF (238 kb)

• Modeling Electroporation in a Single Cell. II. Effects of Io...

  Modeling Electroporation in a Single Cell. II. Effects of Ionic Concentrations Biophysical Journal, Volume 77, Issue 3, 1 September 1999, Pages 1225-1233 Katherine A. DeBruin and Wanda Krassowska Abstract This study expands a previously developed model of a single cell electroporated by an external electric field by explicitly accounting for the ionic composition of the electroporation current. The previous model with non-specific electroporation current predicts that both the transmembrane potential and the pore density are symmetric about the equator, with the same values at either end of the cell. The new, ion-specific case predicts that is symmetric and almost identical to the profile from the non-specific case, but has a profound asymmetry with the pore density at the hyperpolarized end of the cell twice the value at the depolarized end. These modeling results agree with the experimentally observed preferential uptake of marker molecules at the hyperpolarized end of the cell as reported in the literature. This study also investigates the changes in intracellular ionic concentrations induced around an electroporated single cell. For all ion species, the concentrations near the membrane vary significantly, which may explain the electrical disturbances observed experimentally after large electric shocks are delivered to excitable cells and tissues. Abstract | Full Text | PDF (232 kb)

• Model of Creation and Evolution of Stable Electropores for D...

  Model of Creation and Evolution of Stable Electropores for DNA Delivery Biophysical Journal, Volume 86, Issue 5, 1 May 2004, Pages 2813-2826 Kyle C. Smith, John C. Neu and Wanda Krassowska Abstract Electroporation, in which electric pulses create transient pores in the cell membrane, is becoming an important technique for gene therapy. To enable entry of supercoiled DNA into cells, the pores should have sufficiently large radii (>10nm), remain open long enough for the DNA chain to enter the cell (milliseconds), and should not cause membrane rupture. This study presents a model that can predict such macropores. The distinctive features of this model are the coupling of individual pores through membrane tension and the electrical force on the pores, which is applicable to pores of any size. The model is used to explore the process of pore creation and evolution and to determine the number and size of pores as a function of the pulse magnitude and duration. Next, our electroporation model is combined with a heuristic model of DNA uptake and used to predict the dependence of DNA uptake on pulsing parameters. Finally, the model is used to examine the mechanism of a two-pulse protocol, which was proposed specifically for gene delivery. The comparison between experimental results and the model suggests that this model is well-suited for the investigation of electroporation-mediated DNA delivery. Abstract | Full Text | PDF (359 kb)

• …more

Copyright © 1999 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 77, Issue 3, 1213-1224, 1 September 1999

doi:10.1016/S0006-3495(99)76973-0

Biophysical Theory and Modeling

Modeling Electroporation in a Single Cell. I. Effects of Field Strength and Rest Potential

Katherine A. DeBruin^,   and Wanda Krassowska

Department of Biomedical Engineering and Center for Emerging Cardiovascular Technologies, Duke University, Durham, North Carolina, USA

Address reprint requests to Katherine A. DeBruin, Ph.D., Department of Biomedical Engineering, Box 90281, Duke University, Durham, North Carolina 27708-0281. Tel.: 919-660-5131; Fax: 919-684-4488.

Electroporation is the formation of microscopic, current-carrying pores in a lipid bilayer exposed to a large transmembrane potential V_m. The pores are long lived, often surviving in the membrane for up to several minutes and providing pathways for the movement of ions, drugs, and even DNA fragments into the cell. These properties have made electroporation a common tool in biotechnology (Chang et al,Neumann et al), and the medical applications of electroporation are now being realized (River et al,Tsong, 1991,Tung et al,Zhang et al).

However, the process of electroporation is not well understood. Numerous experimental studies have been aimed at revealing the mechanism of electroporation in various types of membranes ranging from artificial lipid bilayers (Chernomordik and Chizmadzhev, 1989,Glaser et al) to red blood cells (Chang, 1992,Kinosita and Tsong, 1979) to chick myocyte monolayers (Jones et al,Jones et al). These studies investigated the properties of pore formation and resealing using pulse charge techniques (Benz et al,Zimmermann, 1982), measured the kinetics of electroporation in voltage-clamped membranes (Chernomordik and Chizmadzhev, 1989,Tovar and Tung, 1992), tracked the movement of ions and fluorescent dyes across electroporated membranes (Kinosita et al,Mehrle et al,Rossignol et al), imaged the transmembrane potential using voltage-sensitive dyes (Hibino et al,Knisley, 1994), and visualized large pores using freeze-fracture electron microscopy (Chang, 1992). With the wide variety in membrane composition and experimental techniques, the literature on electroporation is difficult to compare and often conflicting. A model is needed to help understand the experimental results and draw qualitative, universal conclusions about the electroporation process and the behavior of electroporated cells.

Until recently, the development of theoretical models of electroporation has lagged behind the experimental research, with the available models unable to fully replicate or explain the experimental observations. The first model described the basic biophysics of electroporation using the Smoluchowski equation, which governs the evolution of the pore distribution function in the space of the pore radii (Pastushenko et al). Weaver and coworkers derived the equations of Pastushenko et al. from statistical mechanics and expanded the biophysical description into a numerical model (Barnett and Weaver, 1991,Freeman et al). However, these formulations are mathematically and computationally complex and therefore only suitable for use in space-clamped membranes. Recognizing the need to model electroporation in spatially extended systems, Weaver suggested a cubic cell model for electroporation that consists of two space-clamped membrane patches connected by a resistor (Weaver and Barnett, 1992). This representation captures some features of cellular electroporation, but it does not allow for spatial variation in the transmembrane potential or pore density.

The need for a model that provides a closer relationship between theory and experiments can be fulfilled by the macroscopic model of electroporation recently developed by DeBruin and Krassowska, 1998, Krassowska, 1995, and Neu and Krassowska, 1999, which provides a means for investigating the mechanisms and effects of electroporation in a variety of tissue geometries. To date, the model has been used successfully to reproduce experimental results involving guinea pig papillary muscle fibers exposed to large electric fields (DeBruin and Krassowska, 1998) and to investigate the influence of electroporation on the shock-induced transmembrane potential in a two-dimensional sheet of cardiac tissue (Aguel et al).

Part I of the present study uses the macroscopic model of electroporation as a basis for the development of a model of a single cell electroporated by an external electric field. This model is used to investigate the process of electroporation in a spherical cell, including the time evolution and spatial distribution of the transmembrane potential and pore density as well as the effects of the rest potential and shock strength. The modeling results are compared to experimental data reported in the literature.

The transmembrane potential on the surface of an isolated single cell exposed to an external electric field can be computed using Laplace’s equation, because both the intracellular and extracellular domains are source-free:

(1)

(2)

(3)

(4)

(5)

(6)

(7)

Display large version of this figure

Figure 1

Schematic of a spherical single cell with radius a immersed in a spherical shell of extracellular space with thickness 2a. The electric field E is oriented such that the depolarized pole is at θ=0 and the hyperpolarized pole is at θ=π. All profiles of the transmembrane potential V_m or the pore density N around the cell are plotted from −π/2 to 3π/2.

The pore density N is governed by a first-order differential equation (DeBruin and Krassowska, 1998,Neu and Krassowska, 1999),

(8)

For any cell shape, Eqs. (1) must be solved numerically. The intracellular and extracellular space is discretized using a finite difference method, and the resulting linear system of equations is transformed using LU decomposition. In each time step, the intracellular and extracellular potentials are computed using forward and backward substitution. The results are used to find V_m, I_ion, I_ep, and N at the present time step to be used in the calculation of Φ_i and Φ_e at the next time step. This approach solves the original problem described by Eqs. (1). An alternative approach is to use singular perturbation to derive an asymptotic approximation to those equations. By expanding the potentials in powers of a small parameter ɛ and using only the leading order terms, the time dependence of the boundary conditions disappears and Eqs. (1) become a quasi-stationary system that may be solved at widely spaced time intervals. In this case, changes in the transmembrane potential are driven by the time dependence in Eq. (8) for the pore density N. To achieve both accuracy and computational efficiency, this study developed a combined solution method in which Eqs. (1) are solved during the shock and the singular perturbation approximation is used during resealing. Additional details are given in Appendix B .

As an example, this study uses a spherical single cell with radius a=50μm immersed in a spherical shell of extracellular space with thickness 2a=100μm (Fig. 1). Whenever possible, the cell parameters (diameter, passive kinetics), stimulus protocol (electric field strength, duration), material constants (intracellular, extracellular conductivities), and electroporation characteristics (significant effects at V_m ≈ 1V) are matched to the values reported by Kinosita and coworkers for unfertilized sea urchin eggs (Hibino et al,Hibino et al,Kinosita et al,Kinosita et al,Kinosita et al). When a rest potential of −80mV is required (Chambers and de Armendi, 1979), the reversal potential of the ionic current E_l is set to −83.75mV. For a rest potential of 0mV, E_l is set to 0mV. The shock protocol consists of a 400-V/cm field applied for a duration of 1ms. The electroporation parameters α, V_ep, N_o, and w_o appearing in Eqs. (7) depend on the type of membrane. For this study, the values of α, N_o, and w_o are based on experimental results from artificial lipid bilayers (Glaser et al), whereas the parameter V_ep was altered such that the critical transmembrane potential V_cr at which electroporation becomes significant is approximately 1V. Values for all parameters are given in Table 1.

The constants V_ep and V_cr are related, but they are not equivalent. V_ep is a parameter in Eq. (8) indicating that the change in V_m causes an e-fold increase in the pore creation rate. Hence, V_ep is analogous to a time or length constant. If V_m=V_ep, then the pore creation rate changes by only a factor of e¹=2.7, too small to be detected experimentally. However, if V_m=V_cr ≈ 4V_ep, then the pore creation rate changes by e⁴=55, a factor large enough to cause an experimentally detectable change in the membrane conductance.

Throughout the shock, Eqs. (1) are solved with a combined solution method using varying time steps. When the transmembrane potential of the cell is changing quickly (i.e., charging, discharging), the original equations are solved with a time step of τ_c/32=0.034μs, where τ_c is the time constant of cellular polarization (Hibino et al),

(9)

According to the literature, a spherical cell with a passive resistive-capacitance (RC) membrane exposed to an external electric field will polarize such that the maximum and minimum transmembrane potentials V_m occur at the poles of the cell, and V_m at the equator is equal to the rest potential V_rest (Schwann, 1989). The polarization arises with a time constant of τ_c=1.1μs (Eq. (9)), and the time course of V_m is consistent with the exponential charging expected of an RC membrane (Figure 2A, dashed line). Once charging is complete, the transmembrane potential varies cosinusoidally around the circumference of the cell according to the relationship

(10)

Display large version of this figure

Figure 2

(A) Time course of the transmembrane potential V_m at the depolarizing pole of a spherical cell exposed to a 400-V/cm field. The nonelectroporating single cell (dashed line) charges to its steady-state value with a time constant τ_c=1.1μs, but the charging of the electroporating cell (solid line) is interrupted within 1μs of shock application and V_m settles into a constant value of approximately 1V. (B) V_m around the cell at the end of a 400-V/cm, 1-ms shock. The transmembrane potential for the nonelectroporating cell (dashed line) shows the cosinusoidal shape predicted by RC theory. V_m around the electroporated cell (solid line) is lower and the profile is flattened in the polar regions.

However, when a cell is exposed to a shock that induces larger transmembrane potentials, electroporation occurs and the RC theory fails. The V_m charging transient is interrupted, and the transmembrane potential settles into a nearly constant value of approximately 1V, the critical value of transmembrane potential V_cr required to produce significant electroporation in this preparation (Figure 2A). At the end of a 1-ms shock, the transmembrane potential profile around an electroporated cell is smaller than the profile predicted for an RC cell (Figure 2B). The largest decrease in V_m occurs at the poles where the induced potential is largest, and the smallest decrease is near the equator. The profile has also lost its cosinusoidal shape, appearing flattened as approximately two-thirds of the cell’s circumference has a nearly uniform V_m magnitude of V_cr ≈ 1V.

This dramatic change in the electrical behavior of the cell can be explained with a detailed examination of the time courses of V_m and the pore density N. When exposed to an electric field, the cell initially polarizes with a cellular time constant τ_c=1.1μs. Near the equator, where the induced potential is less than the critical value for electroporation V_cr, the membrane will polarize to the steady-state potential predicted for an RC cell (Figure 3A, θ=3π/8, π/2). The membrane in this region will contain a small, baseline number of pores (e.g., N_o at V_m=0mV), but N does not change significantly during charging (Figure 3B) and the current through these pores does not influence V_m. Near the poles, the transmembrane potential quickly exceeds the critical value for electroporation (V_cr ≈ 1V) and creates a very fast increase in the pore density N (Fig. 3, θ=0, π/8). A portion of the stimulus current is shunted across the membrane through these pores, interrupting the V_m charging transient within about 1μs of exposure to the electric field. In the region between the pole and the equator (θ=π/4), the increase in N is more gradual because the induced potential is smaller.

Display large version of this figure

Figure 3

Time course of (A) V_m and (B) N at five locations around a spherical cell. Near the poles (θ=0, π/8), the V_m transient is initially steep but quickly truncated, and N experiences almost a step increase once V_m exceeds the critical value for electroporation V_cr ≈ 1V. Near the equator, where the induced potentials are smaller (θ=3π/8, π/2), V_m follows the time course for an RC cell. Since V_m≪V_cr, there is no significant increase in N in that region.

After the charging transient, N in the electroporated regions of the cell settles into a slow upward drift because V_m is still greater than V_cr. This continued creation of pores provides additional pathways to shunt current across the membrane, and V_m throughout the electroporated region slowly decreases toward V_cr. This feedback between V_m and N occurs at different rates in different locations. Near the poles, the V_m transients are steeper and create larger pore densities than in the surrounding regions (Fig. 3). With more current pathways across the membrane, the post-transient transmembrane potential is smaller. This situation can be observed at the end of the 5-μs shock shown in Figure 3A, where the transmembrane potential at the pole (θ=0) is smaller than V_m at θ=π/8, which is, in turn, smaller than V_m at θ=π/4. At the θ=3π/8 location, the transmembrane potential is subcritical and no electroporation occurs.

The post-transient differences in V_m create concavities, or dips, in the transmembrane potential distribution around the cell (Figure 4A). The magnitude of these concavities decreases over time as N compensates by increasing nonuniformly (Figure 4B). At the end of a 1-ms shock, V_m is nearly constant at 1V throughout the electroporated regions (Figure 4A, heavy solid line). The transmembrane potential may be considered symmetric about the equator, because the magnitude of V_m is the same at the depolarized and hyperpolarized ends of the cell. The pore density N is also symmetric.

Display large version of this figure

Figure 4

(A) V_m and (B) N around a spherical cell for four time instants during a 1-ms shock. The concavity in V_m near the poles disappears over time, while the pore density distribution gradually widens and increases.

Theory applicable to RC cells predicts that the intrinsic rest potential V_rest of the cell will alter the transmembrane potential profile by shifting it in the direction of V_rest. V_m at the poles would still be symmetric about the equator, but the transmembrane potential at that location would be equal to the rest potential. For example, if V_rest=−80mV, then the cell in this study would have V_m equal to +2.92V at the depolarized pole and −3.08V at the hyperpolarized pole, both 3V from the rest potential. However, in this scenario, the transmembrane potentials still far exceed the critical potential for electroporation, V_cr ≈ 1V, so significant electroporation will occur at both ends of the cell. Intuitively, one would expect that the negative bias of the rest potential would cause the hyperpolarized end to electroporate earlier than the depolarized end, and the correspondingly steeper slope of the V_m transient would produce a larger pore density. Since the cell is a source-free system (Eq. (1)), the net transmembrane current must be zero under all conditions (Krassowska and Neu, 1994). After the first 1–2μs, the capacitive transient is complete, and the electroporation current is much larger than the ionic current. Therefore, one would also expect that V_m at the hyperpolarized end would be smaller than at the depolarized end.

Comparing V_m and N at points around the cell at the end of a 400-V/cm, 1ms shock confirms that intuitive scenario qualitatively (Table 2). With V_rest=−80mV, the pore density at the hyperpolarized pole is larger than N at the depolarized pole, while the opposite is true for the magnitude of V_m. Quantitatively, the asymmetry in V_m and N is very small, and it is unlikely that this minor variation would be detectable experimentally. However, Table 2 also shows a surprising result that can be measured experimentally: V_m at the equator (θ=π/2) is approximately equal to 0mV even if the intrinsic rest potential of the cell is −80mV. This negative offset disappears during the initial charging transient because the nearly step increase in N increases the electroporation current I_ep by four orders of magnitude, making I_ep≫I_ion, the ionic current that supports the rest potential. The electrical behavior of the cell is governed by I_ep, even in regions which are not electroporated. The intrinsic rest potential of the cell plays only a minor role, producing the slight asymmetry in V_m and N observed in Table 2.

The bimodal shape of the pore density distribution around the cell (Figure 4B) is directly related to the cosinusoidally varying magnitude of the transmembrane potential initially induced by the electric field. Larger potentials, such as those near the poles, produce more pores to shunt the extra stimulus current across the membrane. Near the equator, the subcritical V_m does not significantly influence N. Increasing the electric field strength will increase the number of pores throughout the cell in an effort to dissipate the extra stimulus current, and a larger fraction of the cell membrane will attain the critical transmembrane potential V_cr and electroporate. However, the shape of the pore density distribution and the transmembrane potential are qualitatively unchanged (Fig. 5).

Display large version of this figure

Figure 5

(A) V_m and (B) N around a spherical cell at the end of a 1-ms exposure to three electric field strengths. Larger fields did not alter the maximum magnitude of V_m, but did increase the height and width of the pore density profile. As a result, the fraction of the cell membrane with V_m ≈ 1V also increased.

The exception to this behavior occurs when the cell is exposed to a shock that induces transmembrane potentials just over V_cr. For example, if the cell in this study is exposed to a 150-V/cm field, the maximum RC transmembrane potential equals 1.125V instead of the usual 3V. With the smaller field strength, the initial transient in V_m is less steep, fewer pores are formed, and the initial bias in pore density produced by the rest potential becomes important. Figure 6 shows the time course of V_m and N at both poles during a 150-V/cm shock. The hyperpolarized pole electroporates first because of the negative value of V_rest, and N in that region increases by approximately three orders of magnitude (Figure 6B, dashed line). The pore density at the depolarized pole is still small because V_m<V_cr. If the shock ended during this time frame (duration less than 80μs), there would be a significant asymmetry in the pore density profile. As V_m at the hyperpolarized pole becomes less negative, the balance of current increases V_m at the depolarized pole. For shock durations between 100μs and 240μs, N at the depolarized pole is larger as the increasing V_m causes that end of the cell to electroporate. For shocks longer than about 240μs, the difference in N at the two poles become less significant, and, after 500μs, both the transmembrane potential and the pore density are almost symmetric. These results imply that shock strength and rest potential may be important, but only when the cell is polarized to just over the critical V_m with shocks of very short duration. Larger or longer shocks eliminate the effects of V_rest.

Display large version of this figure

Figure 6

Time course of (A) V_m and (B) N at the poles of a spherical cell exposed to an electric field of 150 V/cm. This field induces potentials that barely exceed the critical value of electroporation, V_cr ≈ 1V, at the poles of the cell. V_m experiences some minor fluctuations around V_cr, but the time course of N shows that electroporation at the depolarized pole is delayed by 80μs with respect to the hyperpolarized pole. If the shock were terminated during this time period, a very asymmetric pore density profile could be obtained. After 500μs, N is almost identical at each pole.

When the shock ceases, the cell discharges the potential induced by the electrical field. This process is faster than cellular polarization because electroporation increases the total conductance of the membrane. If V_rest=0mV, the transmembrane potential around the cell discharges to zero within a few microseconds. If V_rest=−80mV, V_m follows a similar time course, discharging to a value very close to 0mV within 1–2μs. The cell requires about 20s to return to its preshock conditions (Fig. 7).

Display large version of this figure

Figure 7

Time course of (A) V_m and (B) N at the depolarized pole of the cell after a 400-V/cm, 1-ms shock. The vertical line in panel A is the V_m trace during the shock, which appears very short on a time scale of many seconds. V_m slowly repolarizes to its rest potential of −80mV over a period of 20s, the time required for the pores to completely reseal and N to return to its preshock value throughout the cell.

The cell’s prolonged recovery period is due to the slow rate of pore resealing, whose time constant can be evaluated from Eq. (8),

(11)

This study developed a computationally efficient model of a spherical single cell with an electroporating membrane. The modeling results demonstrate that electroporation substantially alters the transmembrane potential around the cell. As compared to V_m predicted for an RC cell, electroporation decreases the transmembrane potential throughout the cell and flattens the predicted cosinusoidal profile near the poles. The pore density increases with shock strength such that V_m in the electroporated regions remains nearly constant at 1V regardless of the strength of the applied electric field. After the shock, the pores reseal with a time constant of 1.5s, and complete recovery of the cell to preshock conditions requires approximately 20s. The intrinsic rest potential of the cell was found to have essentially no effect on either V_m or N.

The majority of results reported in this study are similar to experimental observations made by Kinosita and coworkers, who used a voltage-sensitive fluorescent dye to investigate the transmembrane potential induced in unfertilized sea urchin eggs exposed to large electric fields (Hibino et al,Hibino et al,Kinosita et al,Kinosita et al,Kinosita et al). First, the researchers observed that the transmembrane potential throughout the cell was much lower than predicted for an RC cell, similar to the modeling results shown in Figure 2B of this study. The experimental profile of V_m showed a flattening in the polar regions, and a concavity existed at both poles. However, in contrast to Figure 4A of this modeling study, the degree of concavity did not appear to decrease over time. Kinosita and coworkers also found that significant electroporation occurred during the first microsecond of the shock, followed by a slower increase in the electroporation conductance throughout the duration of the shock. This qualitative description of the time course of N agrees with the model’s predictions (Figure 4B), but V_m has more complicated behavior experimentally. This discrepancy may be due to changes in the radii of pores during the shock, a feature not presently included in the model of electroporation. In experiments, the maximum electroporation conductance G decreased by an order of magnitude within the first millisecond postshock, and resealing was not complete after 2s (the longest time interval measured). These results are also consistent with the predictions of this modeling study, in which G at the poles decreased by 85% in the first postshock millisecond because of the non-ohmic nature of the pores, and complete resealing to preshock conditions requires 20s. A similar time course for electroporation was reported for green algae cells (Neumann et al).

Second, Kinosita’s group tested the saturation of the transmembrane potential with shock strength and found that, for sufficiently large shocks, increasing the field strength did not increase V_m. This observation is consistent with the modeling results indicating that larger shocks create more pores, shunting the excess stimulus current across the membrane and limiting V_m to approximately 1V throughout the electroporated region (Fig. 5). Knisley found a similar relationship in his study of rabbit myocytes (Knisley, 1994), in which larger shocks produced a more pronounced decay in V_m such that the transmembrane potential at the end of a 20-ms shock was approximately equal for all field strengths. This saturation of V_m appears to be a phenomenon that is independent of tissue geometry, because it was also observed in both experimental and modeling studies of voltage-clamped lipid bilayers (Freeman et al), one-dimensional fibers (DeBruin and Krassowska, 1998,Krassowska, 1995,Zhou et al), and two-dimensional sheets (Aguel et al).

Third, Kinosita and coworkers observed a disappearance of the rest potential consistent with the modeling results (Table 2). Other researchers (Knisley and Grant, 1995,Teruel and Meyer, 1997) also found that the intrinsic rest potential did not play an important role in the electroporation process. These two experimental studies eliminated the intrinsic V_rest by altering the extracellular ionic concentrations, but the results were the same as those observed with a negative rest potential.

Finally, Kinosita’s group estimated the electroporation conductance G based on their experimental data and reported a maximum G of 4.3×10³mS/cm². The distribution of the electroporation conductance around the cell was bimodal, with the largest value of G at either pole and a small value near the equator. The modeling study produced similar results, with a maximum G of 2.2×10³mS/cm² and a pore density distribution with a qualitatively similar shape (Figure 4B). Both model and experiment found that approximately two-thirds of the cell are significantly electroporated with an electric field strength of 400V/cm. Kinosita and coworkers also calculated the maximum fractional area of the membrane occupied by the pores to be 10⁻⁴ to 10⁻³, consistent with the experimentally and theoretically determined values for artificial lipid bilayers reported in the literature (Chernomordik et al,Freeman et al). This modeling study predicts that the fractional area of the example cell occupied by pores is 2×10⁻⁵, again in agreement with the experimental results.

The steep dependence of the pore creation rate (Eq. (8)) on the transmembrane potential is inherent to the electroporation process and cannot be avoided by choosing different electroporation parameters. For example, manipulating α and N_o within a physiologically valid range (α, 92cm⁻²ms⁻¹ [Glaser et al] to 200cm⁻²ms⁻¹ [DeBruin and Krassowska, 1998]; N_o, 1.5×10⁴ to 1.5×10⁶cm⁻² [Benz and Hancock, 1981,Chernomordik and Chizmadzhev, 1989,Rosenberg and Jendrasiak, 1968]) will not significantly alter the dependence, because the rate of change of the pore density is not strongly affected by either of these parameters. In comparison, dN/dt is exponentially dependent on V_ep, but altering that parameter will change the value of the critical transmembrane potential V_cr, which is determined by experimental data for a particular cell type.

This steep dependence of dN/dt on V_m has two consequences. First, at equilibrium, V_m is the square root of a logarithmic function of N, implying that V_m is almost insensitive to changes in N. This relationship explains the saturation phenomena observable in Fig. 5, where increasing the shock strength from 150V/cm to 400V/cm increased N by a factor of 8.2, but left V_m in the electroporated regions unchanged. Second, the increase in N is a very fast process, and the creation of pores is complete within about 1μs (Fig. 3). This feature of the model does not necessarily contradict the experimental results that show electroporation occurs on a millisecond time frame (Hibino et al) and the critical transmembrane potential V_cr decreases with shock duration (Hibino et al). Instead, it is possible that the slow (millisecond) change in membrane conductance observed experimentally is due to an increase in the radii of the pores, a feature not represented in this model of electroporation. Likewise, the decrease in V_cr for longer pulses may be due to an increase in the pore radius, which increases the current through each pore and decreases V_m below V_cr. Including the effects of pore radius requires a substantial addition to the model that will be the subject of a future study.

Although the model of an electroporating cell successfully reproduced the experimental data published by Kinosita and coworkers (Hibino et al,Hibino et al,Kinosita et al,Kinosita et al,Kinosita et al), it does have additional limitations. First, the model is a simplified description of the extremely complex processes occurring in a cell membrane. Important biophysical elements such as the stretching of cells exposed to an electric field (Isambert, 1998) are not captured. Second, the value of the electroporation parameter V_ep was chosen to give a critical V_m for electroporation of ±1V (value reported by Kinosita’s group), but the values of other parameters were estimated from experiments performed on artificial lipid bilayers and therefore may not be wholly applicable to sea urchin eggs. Third, the macroscopic model of electroporation describes only primary pores, those formed as a direct result of large transmembrane potentials. Secondary pores, which are thought to be a later stage of development that provides transport routes for macromolecules including DNA (Weaver and Chizmadzhev, 1996), are beyond the scope of this model. Finally, the pores have been shown experimentally to be cation selective (Weaver and Chizmadzhev, 1996), but that feature is not included in this model of electroporation.

Despite these limitations, the only significant difference between the experimental results from Kinosita and coworkers and the modeling results reported here concerns the asymmetry of the electroporation process. The majority of Kinosita’s studies show that the transmembrane potential is symmetric around an electroporated cell (Hibino et al,Kinosita et al,Kinosita et al), but the most recent experiments indicate that there may be a transient asymmetry in V_m when the shock is first applied (Hibino et al). More studies of the transmembrane potential are needed, but many researchers have reported an asymmetry in the uptake of marker molecules with entry predominately at the hyperpolarized end of the cell (Djuzenova et al,Gabriel and Teissie, 1997,Knisley and Grant, 1995,Mehrle et al,Mehrle et al,Rossignol et al,Tekle et al,Teruel and Meyer, 1997). Several studies attribute this asymmetry in uptake to the rest potential, because the negative value is thought to bias electroporation toward the hyperpolarized pole (Djuzenova et al,Gabriel and Teissie, 1997,Mehrle et al,Mehrle et al,Tekle et al). The findings of this modeling study imply that this hypothesis is only valid when the induced potential is very near the critical value for electroporation (Figure 5 and Figure 6). All the experimental studies quoted here use electric fields that far exceed the critical value, and, in those cases, the model predicts that V_rest will cause only a very minor asymmetry in the transmembrane potential and the pore density. Thus, these modeling results rule out V_rest as a cause of the asymmetric uptake of marker molecules.

Other factors must be considered to explain this experimentally observed asymmetry. First, the lipid bilayer itself may be asymmetric, in which case the polarity of the shock would affect the local creation of pores (Genco et al). Second, there may be interactions between the pores and the ionic channels, proteins, and other structures in the membrane that are not replicated by the model. Finally, electroporation may be influenced by different ionic concentrations in intracellular and extracellular space (Djuzenova et al,Knisley and Grant, 1995,Tekle et al). This last hypothesis will be investigated theoretically in Part II of this study, which focuses on the interaction between electroporation and ionic concentrations (DeBruin and Krassowska, 1999).

Appendix A

Origin of Eq. (8) governing pore density

This Appendix shows the connection between Eq. (8), used in this paper to compute the pore density, and the existing theory of electroporation. It is based on the results of a study by Neu and Krassowska, 1999, who derived Eq. (8) as an asymptotic limit of the Smoluchowski equation, generally recognized in the literature as describing the biophysical mechanisms of electroporation.

Neu and Krassowska assumed a relationship between the pore radius and the pore energy that was proposed by Chizmadzhev and colleagues (Abidor et al,Glaser et al). As shown in Fig. A1, the energy E(r) of a pore with radius r is the lesser of the two curves,

(A1)

the energy of nonconducting (hydrophobic) pores, and

(A2)

the energy of conducting (hydrophilic) pores. In Eqs. (A1), r_* and E_* are the minimum radius and energy barrier for the creation of conducting pores (Fig. A1), γ is the pore edge energy, σ is the membrane surface tension, and C is a constant. The third term in Eq. (A2) represents the steric repulsion between the lipid heads lining the pore (Israelachvili, 1992) and is responsible for the increase in pore energy with shrinking radius (Weaver and Chizmadzhev, 1996).

Display large version of this figure

Figure A1

The energy of a pore as a function of radius at the transmembrane potential V_m=0mV. The dashed and solid lines show the energy of hydrophobic and hydrophilic pores, respectively. To better illustrate the relationship between the two pore types, the plot shows the pore energies only for small pore radii.

The pore energy E(r) in Fig. A1 corresponds to the situation when there is no externally applied transmembrane potential. In the presence of a transmembrane potential V_m, the pore energy, denoted by φ(r), is given by

(A3)

where the term −πa_pV_m²r² is the capacitive contribution (Abidor et al,Weaver and Mintzer, 1981). The coefficient a_p can be estimated based on a continuum model as (Glaser et al,Powell and Weaver, 1986)

(A4)

where h is the membrane thickness, κ_w and κ_m are dielectric constants of water and membrane, and ∈_o is the permittivity of a vacuum.

Given the pore energy, electroporation is described mathematically by the Smoluchowski equation (Barnett and Weaver, 1991,Freeman et al,Pastushenko et al,Powell and Weaver, 1986,Weaver and Mintzer, 1981). If n(r, t) denotes the pore density distribution function such that at a given time t, the number of pores per unit area with radii between r and r+dr is n(r, t)dr, then n(r, t) is governed by the equation,

(A5)

where D is the diffusion coefficient of pores, k is the Boltzmann constant, T is the absolute temperature, and S(r) is the source term that represents the creation and destruction of pores. S(r) can be written as

(A6)

where ν_c is the attempt rate density (Weaver and Mintzer, 1981), ν_d is the frequency of lipid fluctuations (Glaser et al), and U denotes the pore energy φ of nonconducting pores (r<r_*). H(r), the Heavyside step function, represents the fact that only nonconducting pores are destroyed.

The Smoluchowski equation (Eq. (A5)), used with constants typical for electroporation, contains several small parameters. Their presence facilitates the use of singular perturbation to perform a rigorous simplification of Eq. (A5), and such an asymptotic reduction (Neu and Krassowska, 1999) transformed the Smoluchowski equation into an ordinary differential equation (ODE). This ODE describes the dynamics of the pore density N(t), which is related to the pore distribution function n(r, t) by

(A7)

The asymptotic ODE for N(t) has the form

(A8)

In the quasistatic case, K and N_eq are given by Eqs. 77–78 of the paper by Neu and Krassowska,

(A9)

(A10)

Substituting Eqs. (A9) into Eq. (A8) yields Eq. (8) used in the main body of this paper.

The paper of Neu and Krassowska relates the coefficients of Eqs. (A8) to constants appearing in the expressions for pore energy (Eqs. (A1)) and in the Smoluchowski equation (Eq. (A5)) (Neu and Krassowska, 1999):

(A11)

(A12)

(A13)

(A14)

Eqs. (A8) are the dimensional versions of Eqs. 68–78 from Neu and Krassowska. Energies E, φ, and U are in units of kT, and U′, φ′, and φ′_m denote derivatives with respect to r evaluated at r_* and r_m.

In application to a single cell, the following simplifications were made. First, the formulation given above represents a quasistatic limit, i.e., it is assumed that the pore distribution function n adjusts instantaneously to temporal variations in pore energy. As argued in the original study (Neu and Krassowska, 1999), this approximation is valid when the changes in V_m occur on a time scale of at least 5μs. Here, cellular polarization has a time constant of 1.1μs, so the quasistatic approximation introduces an error. However, since this assumption affects only the coefficient N_o, one can expect only a modest difference between solutions using the quasistatic and time dependent versions of the asymptotic ODE.

Second, the model used here suppresses the dependence of α and N_o on the transmembrane potential and treats them as constants. This simplification is acceptable because the dependence of dN/dt on V_m is dominated by the exponential exp[(V_m/V_ep)²]. In comparison, the dependence on V_m through α and N_o is much weaker and is unlikely to be detectable experimentally. The radius at the minimum pore energy r_m also depends on V_m, but it changes very little for V_m between 0mV and the critical value V_cr (Neu and Krassowska, 1999). Hence, r_m and, consequently, q are constant.

In principle, Eqs. (A11) can be used to determine values for the parameters of the model. However, this method would use several molecular-level constants whose values are known only up to an order of magnitude (Barnett and Weaver, 1991). Alternatively, the four parameters can be determined experimentally. Glaser et al performed specially designed voltage-clamp experiments on artificial lipid bilayers that yielded estimates for α and V_ep. The single cell model adopted Glaser’s value for α, but decreased V_ep from 460 to 258mV so that V_cr ≈ 1V, the value reported by Kinosita and coworkers for unfertilized sea urchin eggs (Hibino et al,Hibino et al,Kinosita et al,Kinosita et al,Kinosita et al). N_o was computed by dividing the measured background conductivity of a lipid bilayer by the conductance of a single pore (Benz and Hancock, 1981,Chernomordik and Chizmadzhev, 1989,Rosenberg and Jendrasiak, 1968). Finally, q was chosen based on the experimental estimates of Glaser and coworkers for r_* (0.3–0.5nm) and r_m (0.6–1.0nm) (Glaser et al).

Appendix B

Singular perturbation approximation to Eqs. (1)

For the pore resealing phase, this study uses singular perturbation to develop approximate, quasistationary equations governing the intracellular and extracellular potentials Φ_i and Φ_e. Once the shock has ceased and the induced potential has been discharged, V_m assumes a nearly constant value V_* everywhere around the cell,

(B1)

where G_*=89.32mS/cm² is the average conductance of the electroporated membrane as determined from simulations. G_* is due to pores remaining in the membrane after the shock, as they reseal with a time constant τ_N=1.5s (Eq. (11)). Recognizing that τ_N is 10⁶ times larger than the cellular time constant τ_c=1.1μs (Eq. (9)) motivates the use of singular perturbation. The method used here is similar to the one proposed for an excitable cell in an external electric field (Krassowska and Neu, 1994). The first step is to convert the governing equations into nondimensional form using the system of units shown in Table B1. Equations (1) remain unchanged because they are invariant under scaling. Equation (4) for the boundary conditions on the membrane S is written as

(B2)

where μ≡σ_e/σ_i, κ≡g_lN_o/αC_m, and ν≡d_cG_*/σ_i are O(1) constants and ɛ≡τ_c/τ_N=d_cαC_m/σ_iN_o=1.4×10⁻⁶ is a small parameter. The presence of this small parameter in Eq. (B2) allows the expansion of potentials in powers of ɛ. For Φ_i,

(B3)

and similar expansions are written for Φ_e and Φ_m. Substituting these expansions into Eqs. (1) gives

(B4)

(B5)

Since ɛ ∼ 10⁻⁶, only the leading order terms will be considered. The contribution of the first-order terms is less than 0.1% as determined by simulations. Collecting powers of ɛ and discarding all but the leading order terms results in a simplified system of equations governing ϕ_i⁰,

(B6)

(B7)

Analogous equations can be derived for the extracellular potential Φ_e ∼ ϕ_e⁰,

(B8)

(B9)

(B10)

With I_ion and I_ep known from the previous time step, the system of equations for ϕ_i⁰ and ϕ_e⁰ can be treated as a time-independent boundary value problem. Eqs. (B6) are converted to spherical coordinates and discretized in r and θ using the finite difference method. The resulting linear systems of equations is solved in each time step using Gaussian elimination. The transmembrane potential is computed from ϕ_i⁰ and ϕ_e⁰ and used to update N and calculate I_ion and I_ep. The time step during resealing is governed by the convergence requirements for Eq. (8) describing the rate of change of N. With τ_N=1.5s the maximum time step is 100ms, and a 20s simulation (complete resealing) can be completed in 7 minutes, 44 seconds on a Sun Ultra 1. If, instead, the original problem (Eqs. (1)) is used, the computational time is estimated to be 1800 hours. The substantial savings of the singular perturbation approximation make it feasible to conduct investigations of the resealing process in an electroporated single cell.

Table B1 Singular perturbation scaling units

	Parameter	Unit	Typical Value
	x	dc=2a	100μm
	t	No/α	1.5s
	σi, σe	σi	4.55mS/cm
	Vm	V*	−0.178mV
	Iion	gl(V*−Vrest)	15.9μA/cm2
	Iep	G*V*	−15.9μA/cm2