Background and Introduction

Cancer treatment often fails at the cellular and tumor levels. At the cellular level, cancer cells frequently do not respond to pharmaceutical treatments because of acquired multiple drug resistance (1), inhibitors of apoptosis (2, 3), and inhibition of signaling molecules (4). Molecular mechanisms inhibiting apoptosis pathways can also hinder physical therapies such as localized ionizing radiation as well as systemic pharmaceutical interventions (5, 6, 7). At the multicellular level there are additional barriers within solid tumors to drug therapies that arise from the inability of drugs to fully penetrate abnormal, heterogeneous, and irregularly vascularized tumor tissue and thereby not reach all of the cancer cells at therapeutic levels (8, 9, 10, 11). Accordingly, local physical therapies that universally kill all cells within a selected tissue volume are of great interest. With this as motivation, we consider some fundamental mechanistic aspects of tissue ablation by IRE (12, 13, 14, 15) that are relevant to solid tumor treatment.

At the outset, it is important to recognize that electroporation is not simply punching small holes in a cell membrane. Instead, a widely accepted hypothesis of electroporation involves dynamic transient pores (16, 17). These pores are created electrically, expand and contract agilely in response to elevated values of the time-dependent PM transmembrane voltage, Δψ_PM, and are destroyed randomly with a mean lifetime estimated to range from milliseconds to even minutes or hours. Electroporation involves highly non-linear local electrical interactions with membranes. This results in rapid redistribution of electrical currents and fields by feedback through heterogeneous, time-dependent local populations of conducting pores, mediated through pathways within the fixed intra- and extra-cellular electrolytes. Thus, the entire cell system interacts and exhibits emergent behavior. Examples of emerging pore sites include the polar regions of spherical cells (18, 19, 20), the analogs to polar regions of cylindrical cells (21, 22, 23), and local sites in closely spaced cells within the irregular geometry of a solid tissue (24). These membrane models assume neither where electroporation will occur nor the magnitude of Δψ_PM at which significant numbers of pores will be created. Instead electroporation emerges on the basis of the underlying physics-based electroporation model and as the result of the electric interaction throughout the whole system model. Previously, dynamic pore distributions could be described only for the simple geometry of artificial planar bilayer membranes (25, 26, 27), two faces of a cubic cell model (28, 29), or the membrane of an isolated spherical cell (20). As demonstrated here, conventional electroporation of irregular cell membranes in a multicellular environment can now also be described, giving new insights in microdosimetric conditions, and leading towards a mechanistic understanding of tissue electroporation properties.

Pores may form in virtually any cell for strong but brief (nanoseconds to milliseconds) exposures to electric fields ranging from about 100 V/cm (longest pulses) to 100,000 V/cm (shortest pulses). PM electroporation itself is essentially universal because electroporation within the phospholipid regions of cell membranes appears inescapable once Δψ_PM reaches several hundred millivolts. The Δψ_PM value at which pore creation increases depends not only on the rate at which Δψ_PM rises, but primarily on a complex kinetic behavior within the energetic landscape of pore states leading to hydrophilic pores (27). Contrary to many literature statements, there is, therefore, no fixed voltage threshold for electroporation (27). The resulting pore densities can provide a membrane conductivity well beyond that of all ion channels. This high-conductivity state protects the membrane, and although this is often termed reversible electrical breakdown (REB), it actually reflects an agile but gentle structural rearrangement of the phospholipids (16, 17), as demonstrated in molecular dynamics simulations (30, 31). Significantly, both reversible and irreversible electroporation have been observed in bilayer membranes (32). Irreversible behavior is attributed to bilayer rupture by uncontrolled pore growth and the outcome is governed by the local Δψ_PM behavior (magnitude, shape, and duration), and the relatively large membrane tension (25).

The earliest reports known to the present authors of irreversible (33) and reversible (34) electroporation in cells were observed in the nodes of Ranvier of myelinated nerves. Subsequent studies of microorganisms showed cell killing attributable to irreversible membrane changes (35, 36, 37). These early reports are consistent with experimental observations of a non-thermal contribution to electrical injury. The latter led to the hypothesis of electroporation-induced necrosis within living tissue by Lee and co-workers (38, 39, 40). Their motivation was the recognition that electrical injury is often characterized by the preferential death of large mammalian cells (skeletal muscle, nerves) in tissue regions where insignificant temperature rise occurs.

PM permeabilization and the resulting pore sizes in cells are fundamental to understanding necrotic cell death by IRE. However, a basic distinction needs to be made as there appear two distinct mechanisms that can lead to necrotic cell death and, thus, irreversibility: (i) evolution of a pore population such that a significant number or size of pores become trapped or held open for sufficiently long times such that a lethal biochemical imbalance is created, and (ii) evolution of pores such that a lethal biochemical change happens within the cellular compartment, even if pores vanish in a relatively short time (17) and, thus, electroporation itself is reversible. It is presently unknown which mechanism is most relevant for necrotic cell death by IRE.

With this background and motivation it is gratifying that purposeful tissue ablation by IRE is being pursued for clinical use as a promising approach to solid tumor therapy (13). Significantly, the electric field pulse waveforms (magnitude, shape, and timescale) may be different from those generally used in tissue ablation by electrochemotherapy (ECT) (41, 42, 43). For example, Miller et al. reported cancer cell ablation by IRE using three repetitions of 1,500 V/cm pulses with 300 ?s duration (12). A further description of IRE was provided by Edd et al. (14) by using a single 20 ms square pulse of 1,000 V/cm. Clinical implications of IRE and relevant waveforms, as described by Rubinsky et al. (15), use a train of up to eight 100 μs trapezoidal pulses, temporal conditions similar to those applied in ECT.

Methods

We use the transport lattice (TL) method, which allows for a convenient description of electrical, chemical, and thermal behavior in a complex biological geometry that may contain inhomogeneities and anisotropies. Basic features of the TL method have been presented elsewhere (21, 22, 23, 24, 44, 45). Here we describe two system models, a multicellular model of irregular cells and a tissue model, and consider their electrical and thermal responses to some characteristic IRE pulses. Each system model represents rabbit liver tissue but on a different spatial scale. The modeling methods described here might also be useful in describing conventional electroporation in cell pellets (46, 47) and in dense cell suspension (48).

Multicellular Model

The system geometry, shown in Figure 1A, is based on a drawing motivated by a tissue section image. It features a layer of 20 liver cells with 14% interstitial fluid volume (49). The hepatocytes have an average cell size diameter of 21.7 μm (42). A corresponding TL (101 nodes × 101 nodes) was constructed as a large electric circuit comprising ∼ 10⁴ local models for passive charge transport and storage (resistors, capacitors) within electrolytes and active elements (pumps, electroporation) at the membrane (Fig. 1B). The linked local membrane and electrolyte models are distributed spatially and are connected to their nearest neighbors on a Cartesian lattice (21, 44). The local membrane models are interconnected at the regularly spaced nodes, with submodels that represent the PM and two contacting regions of electrolyte (Fig. 1B) (21, 44). There is no transport in the z-direction in this two-dimensional model. The lattice spacing, as well as the depth of the system model, is 1 μm, leading to a spatial scale of 100 μm × 100 μm. Voltages applied along the top and bottom boundary of the system model provide the applied uniform field.

Multiscale Tissue Model

The tissue system comprises a large tissue region and two ideal cylindrical electrodes, each with radius r_e = 0.25 mm, separated by L_e = 10 mm (Fig. 1C) (50). The nominal applied electric field, E_app is defined here as the voltage difference between the electrodes, V_app, divided by the electrode spacing, L_e. Although needle electrodes actually have spatially varying fields, the term nominal electric field may be used for convenience, and it is also in line with the original definition in the limit of infinite needle radii (i.e., planar electrodes) (51). The tissue system with scale 200 mm × 200 mm is symmetric about y=0. Thus, a no-flux boundary was placed at y=0 and only the region y ≥ 0 was actually meshed and simulated. Nodes were optimally distributed throughout the tissue system using a meshing algorithm developed by Persson and Strang (52). Figure 1D shows the mesh near the electrodes. The multiscale tissue model accounts for the electrical response at both the microscopic (e.g., PM electroporation) and macroscopic (e.g., needle geometry) scales and the interplay between the two (50). The scale of the tissue system is orders of magnitude larger than the scale of the cells and therefore a discretization of the system could not realistically resolve individual cells and membranes. The multiscale model uses representative simple cell models distributed throughout the system model to calculate the local cell and membrane response, and the macroscopic electrical transport properties are determined by the distributed models.



Figure 1: (A) Geometry of the multicellular model of a region of liver tissue with 14% interstitial space. The extracellular electrolyte is shown in blue, the cells in red. (B) Circuits: Functional local models that represent electrolyte (M_el), membrane (M_m), and the PM-electrolyte interface in the 101 × 101 TL; details described elsewhere (21, 22, 44). The equivalent circuit representing the asymptotic electroporation model (left gray box) (21) and its extension to the SE (right gray box) (54) is solved at every local membrane site. The pore distribution is discretized, and pore drift and diffusion in pore-radius space determine the associated non-Ohmic pore conductance that is input for the membrane current I_m in module M_m. (C) The tissue system model contains two circular electrodes in a 200 mm × 100 mm tissue region (only a small subregion of the entire simulation region is shown). The electrodes have radii r_e = 0.25 mm and separation L_e = 10 mm (edge-to-edge). (D) The tissue system mesh (only a small subregion of the mesh close to the needle electrodes is shown). (E) The tissue model cell unit comprises extracellular (e), membrane (m), and intracellular (i) regions in series and a parallel shunt region. The spatial and electrical parameters of the regions are labeled. (F) The passive tissue conductivity of the tissue model, obtained from an ac frequency sweep in SPICE, is compared with experimental measurements on rat liver by Raicu et al. (53).

Tissue Level: The impedance of a region of tissue is equal to the impedance of a cell scaled to have the same relative dimensions, assuming that tissue comprises a uniform grid of such cells. Thus, we use the simple cell model shown in Figure 1E, which has a membrane enclosed region of intracellular electrolyte surrounded by extracellular electrolyte (50). This simple model can be straightforwardly translated into an equivalent circuit. The membrane and each region of electrolyte have an associated conductivity σ and permittivity ε (Table I). Additionally, each electrolyte region has a tortuosity ν to account for the structural complexity of tissue not otherwise represented by the model. The relative sizes were chosen such that 14% of the total volume was extracellular. The tortuosities ν were used as free parameters in fitting the frequency-dependent rat liver tissue conductivity to that measured experimentally by Raicu et al. (53). By this choice of electrical and geometrical parameters, the model reproduces the static conductivity value and approximately the same trend in the frequency dependence (Fig. 1F). The equivalent circuit for the simple cell model (Fig. 1E) is placed between each pair of adjacent nodes in the mesh with electrical components scaled to the local mesh geometry. The effective conductivity of the membrane changes in accordance with the local degree of membrane electroporation as determined by the distributed cell models.


Cell Level: An equivalent circuit for a single cell is created for each node in the tissue level mesh to determine the cellular response to the local electric field. Each of these circuits is distinct from the primary, macroscopic tissue-level circuit network, but all of the circuits are solved simultaneously. The asymptotic model of electroporation is used as described below. The voltage across the cell unit is equal to the local electric field magnitude, as determined from the electric potential of the nodes in the mesh, multiplied by the cell unit length l_u. The distributed cell models determine the transmembrane voltage and pore density throughout the tissue domain as functions of time. The pore density then determines the membrane conductivity used in the macroscopic transport network. Thus, there is continual feedback between the macroscopic (tissue level) and microscopic (cellular level) models: the macroscopic behavior determines the local electric field in the microscopic model and the subsequent behavior at the microscopic scale (e.g., electroporation) then determines the local electrical properties at the macroscopic scale.

Electrolyte and Membrane Models

The passive electric components for the electrolyte are resistors and capacitors (44) (Fig. 1B). The membrane circuits (Fig. 1B) include components for charge storage and conduction, resting potential, and the asymptotic or SE model of electroporation (21). This provides a convenient means for combining the d_m = 5 nm thick membrane with a TL of much larger scale. The dielectric constant of the extra- and intracellular electrolytes, ε_e, is 72. Close to a pore the membrane dielectric is treated as pure lipid and assigned a dielectric constant, ε_l, of 2.1. This choice recognizes that local membrane properties are relevant to pore formation. In contrast, the PM capacitance involves a spatial average over membrane lipid and protein regions resulting in a relative permittivity, ε_m, of 5. The extracellular electrolyte has a conductivity, σ_e, of 1.2 S/m while the medium inside the cell has a conductivity, σ_i, of 0.4 S/m. These conductivity values are identical to the validated macroscopic model. Following Läuger (54) we use a simplified, single resting potential source model (Fig. 1B) comprised of an active voltage source, V_ip, and source series resistance, R_ip (21). Here, the fixed quantities V_ip and R_ip, together with a negligible conductivity of the equilibrium pores (27), determine the membrane resting potential, ψ_PM,rest, in the absence of applied fields.

Electroporation Models

The transient aqueous pore hypothesis of electroporation is based on continuum models of membrane pores, mechanical and electrostatic energy differences W, and thermal fluctuations, usually in the form of the Smoluchowski equation (SE) (25, 26, 27, 55, 56, 57):


The SE describes the evolution of the local PM pore distribution n in terms of the number of hydrophilic pores and pore radius r_p, given by the diffusion constant D_p. The asymptotic model (56) is an approximation to the SE-based models that disregards pore size change, and can be used to describe cell and tissue responses to nanosecond time scale pulses (22, 24). For longer ECT and IRE pulses, the asymptotic model approximates the system electrical response. These pores readily transport Na⁺, Cl^-, and K⁺ ions that dominate extracellular and cytosolic conductivity, but not significantly larger molecules. The implementation of the asymptotic electroporation model in terms of an equivalent circuit is described in detail elsewhere (21). It has been recently generalized to include dynamic pore expansion and contraction based on the SE electroporation model (57). The extended equivalent circuit is given in Figure 1B and represents drift and diffusion in pore radius space from the minimum pore radius (r_min = 0.8 nm) to a maximal pore radius (considered here as r_max = 3 nm). From the pore distribution, we find the local non-Ohmic membrane conductivity σ_m (27), that, together with the local Δψ_PM, determines the membrane current I_m as input into the membrane circuit M_m. This approach describes the PM response to conventional electroporation, for which the pulse duration may be much longer than the PM charging time, typically 0.2 to 1 μs for isolated cells. Irregular shaped cells close together in a tissue have a heterogeneous and generally longer PM charging time due to long and often tortuous aqueous pathways within the interstitial space. Pore lifetimes reported in the literature vary over many orders of time (from milliseconds to hours) (17) and there is presently no mechanistic understanding of this large range; we use an illustrative experimental value, 3 ms (58), but note that this paper only considers electric conditions during the pulse for which the value of the pore lifetime is not relevant. All electroporation parameters are given in Ref. (27).

TL Solution

The system model circuits are solved for the electric potential φ by means of Kirchhoff?s laws using Berkeley SPICE 3f5. SPICE generates solutions that are processed and displayed in Matlab as equipotentials and distributions of electroporated regions (21, 22, 23, 24).

Specific Absorption Rate

The specific absorption rate (SAR), defined as Joule heating σ|∇& phi;|² divided by the tissue (electrolyte) mass density ρ, is traditionally used in the electrical characterization of tissue to time-varying electromagnetic fields (e.g., cell phones and MRI exposures). At a given time point, SAR is calculated as follows (44). In a local volume of electrolyte, represented by model M_el (Fig. 1A), a SAR contribution of V(t)²σ_el/(ρl²)/2 is added at the two nodes (in the SAR array) between which the model is connected. Here, σ_el is either the intracellular or extracellular electrolyte conductivity, l is the lattice spacing, and V=Δφ is the voltage difference between the two nodes at time t. In the case of a membrane model (M_m; Fig. 1A), the SAR contribution is (V₁(t)²σ_e/(ρl²))/2 + Δψ_PM(t)²σ_m(t)/(ρd_m²) + (V₂(t)²σ_i/(ρl²))/2 where V₁(t) and V₂(t) are the voltages across the two resistors that represent the extracellular and intracellular electrolytes contacting the membrane, and the membrane thickness d_m << l.

Results

First, we consider a single 100 μs trapezoidal pulse with 1 μs rise and fall times that represent widely used ECT pulses and some IRE pulses. We consider both 400 V/cm and 700 V/cm magnitudes, corresponding for those pulse waveforms, respectively, to the recently reported thresholds for reversible and irreversible electroporation in liver tissue (42, 59). Figure 2 shows the multicellular responses to a 700 V/cm pulse at the end of the pulse (t = 99 μs) for a passive membrane (without electroporation) and a membrane with the SE electroporation model. For the passive response, the current is entirely driven through the interstitial space, and the voltage contours are concentrated around the cells, indicating elevated transmembrane voltages at the PM. In striking contrast, electroporation leads to local high-conductivity transients that cause a strong redistribution of the current being driven partially through the cell interior. The voltage contours are still concentrated around the cells, but to a smaller extent. This reflects electroporation limiting the magnitude of the PM transmembrane voltage to values Δψ_PM < 1.5 V (not shown).



Figure 2: Passive versus active electric response for the 700 V/cm trapezoidal pulse at the end of the pulse at t=99 μs. (A) Passive model response (without electroporation): Currents are driven entirely through the interstitial space and the intracellular electric fields are negligible. Voltage contours are concentrated around cells indicating elevated transmembrane voltages at the PM. (B) The active SE electroporation model response: Currents go through cells and interstitial space, leading to significant electric fields inside the cells. The voltage contours are still concentrated around cells, but to a smaller extent, indicating lower PM transmembrane voltages than in (A).

Figure 3 shows the distributed response to the 400 V/cm pulse with the SE electroporation model at different time points during the pulse. White dots correspond to local membrane areas with at least one local pore. Prior to the end of the pulse rise time (t ≤ 1 μs), we find the electric field penetrating the entire multicellular region due to membrane displacement currents. Intracellular electric fields decline substantially after the initial charging phase, and it requires electroporation to re-establish significant intracellular fields. Initial electroporation is seen at t = 4.6 μs at distinct sites. The electroporation pattern evolves with time and includes more and more membrane sites. At the end of the pulse, not all the cells are found to be electroporated; smaller cells are spared as conventional electroporation is strongly cell-size dependent. Significantly, however, our results agree with the notion of a minimum field strength that is needed to reversibly electroporate all cells in liver tissue (42).



Figure 3: Electrical response of multicellular model to 400 V/cm pulse at different time points: (A) Charging phase at t ≤ 1 μs: Displacement currents at the PM lead to intracellular electric fields and currents through the cell. (B) After the charging phase (e.g., t = 4.6 μs), the intracellular electric fields are negligible again, initial electroporation occurs at distinct sites (indicated by white dots corresponding to more than one local pore). (C) At t = 21 μs, the electroporation pattern has spread laterally but the intracellular fields remain small. (D) At the end of the pulse at t = 99 μs, we find still more electroporated sites, and some intracellular equipotentials. Not all cells are electroporated, as the cell size is important.

Continued