A general, self-consistent scheme for analyzing cellular electroporation for bio-medical applications is developed to probe realistic biological shapes and different length scales ranging from nanometers to hundreds of micrometers. The COMSOL Multiphysics suite is used with suitable embellishments to incorporate the details of the electroporation (EP) process and the inherent internal physics. The results are obtained for the voltage pulse driven electroporation for a Jurkat cell with mitochondria (as an example organelle) where spatial dimensions on the order of a few nanometers become important, to hundreds of cells (with Bacillus as an example) where collective effects and mutual interactions can dominate. Thus, scalable computing to generalized geometries with the ability to include complex organelles is made available. The results obtained for mitochondrial EP in Jurkat cells compare well with available data. In addition, quantitative predictions of field attenuation and shielding in Bacillus clusters are made, which point to highly nonuniform field distributions and a strong need to engineer novel electrode designs.
I. INTRODUCTION
It is well known that electric pulses of sufficient duration and intensity can induce a variety of effects in cells and their organelles, such as membrane electroporation,1–12 increase in reactive oxygen species,13–17 capacitive intracellular entry of calcium, which can disrupt intracellular calcium homeostasis and trigger a cascade of other events,18–21 changes in the mitochondrial transmembrane potential (TMP) with the likely opening of the mitochondrial permeability transition pore (mPTP),22,23 and adverse effects on adenosine triphosphate (ATP) dynamics. A large influx of Ca2+ after electroporation (EP) can also lead to intracellular ATP depletion and the inhibition of ATP production in mitochondria.24,25 Furthermore, ATP can leak through the permeabilized plasma membrane from the cell,26,27 which then threatens cell survival.
Electroporation is probably the most studied of the various bioeffects given its various applications that include electrochemotherapy28–30 for tumor treatments, gene delivery,31,32 cell fusion,33 micro-organism inactivation,34 and immunotherapy,35 just to name a few. It is an effective non-viral approach to transfection that is independent of cell surface receptors. More recent uses include tissue ablation based on irreversible electroporation (IRE) that mostly employs 20–200 pulses over 70–100 ms with fields from 1000 to 2000 V/cm, or milli-second excitations consisting of single pulses of 1–20 ms duration with fields from 1000 to 2000 V/cm.36 High-frequency irreversible electroporation (H-FIRE) is yet another related modality that has shown great promise for tumor treatment while mitigating pain and muscle contractions.37–39 The H-FIRE protocol typically uses 50–200 bursts of bipolar microsecond electric pulses, each burst at a frequency of 1 Hz containing about eight pulses having a 100 μs energized duration with electric fields in the 500–4000 V/cm range.40 Pulsed field ablation has also been used to treat cardiac arrhythmias using different combinations of pulses with varying electrical characteristics.41,42
Despite the relatively simple conceptual path toward various electrically driven bioengineering applications, the biological system itself is very complicated, making it difficult to discern and unentangle direct effects from indirect triggers and pathways. For example, pulsed electric field effects on mitochondria are well known and have been studied by several groups.20,23,43,44 In this regard, changes in the permeability and dissipation of the mitochondrial membrane potential (MMP) have been observed. The increased outer mitochondrial membrane (OMM) permeability has been discerned based on the observations of cytochrome c release.45,46 However, it is not yet clear whether electric pulse-induced EP occurs predominantly at the OMM or the inner mitochondrial membrane (IMM), or if both processes are operative to varying degrees. Furthermore, it has not yet been determined if the observed outcomes are a direct effect of mitochondrial membrane electroporation, or if it is triggered indirectly via the opening of the mPTP, or other apoptosis-related pathways. For example, it may be possible that the dissipation of mitochondrial membrane potential (MMP) first causes swelling and rupture of the OMM, leading to eventual cytochrome c release. In any event, many of the processes are most likely very interconnected and driven by electric fields, making it germane to quantitatively study and first understand the degree to which voltages could be modulated across the membranes of such organelles.
In the above context, the IMM is of particular interest as it is known to be a major site of the electron transport chain and ATP synthase. Cristae are involuted structures that form tubules or lamellae that are part of the IMM. It is hypothesized that the role of crista morphology is to increase the surface area of the IMM to enable greater oxidative phosphorylation. It is also conjectured that the proton concentration and proton motive force (PMF) are higher at the crista portions of the IMM.47 Such proton concentration and gradients are essential for ATP synthesis. Furthermore, the sharper membrane curvatures at the sites of the various cristae suggest that the induced electric field distributions at some cristae (at least those whose axis of symmetry roughly lines up with the direction of the applied electric field) might favor enhanced EP over certain crista sites and segments. If so, the localized EP at the IMM could be the trigger that leads to the unfolding of multiple downstream bioeffects. Hence, investigating the electric field distribution and its impact on electroporative outcomes in these pleomorphic IMM structures becomes crucial in understanding the progression of electric field driven dysfunction associated with the mitochondria.
Apart from an appropriate simulation study of sub-cellular organelles such as mitochondria with their intricate geometry and small length scales (in the nanometer range), it is also important to simulate and understand other aspects that could be driven by collective effects. This realm could include biological entities that might be spread out over a larger spatial span. In such cases, correlation, mutual coupling, and interactions, as well as spatial changes in the field distribution, may become important. An example in this regard would be tissues consisting of a high number of closely packed individual cells of varying geometries. Cells lying at the periphery could then shield interior regions from the externally applied fields. Such an outcome might require a clever engineering design of multi-pronged electrode applicators, or even possible development of phased-array or asymmetric cathode–anode systems48 to appropriately tailor the electric field drivers. Tissue conductivity is typically anisotropic and varies between different parts of the body; therefore, specific electro-transfer protocols or electrode geometries and configurations would be needed to enhance electro-transfer mediated delivery. As a simple example of complex tissues and clusters, here we choose to look at different geometric arrangements of Bacillus bacteria to probe polarization and the extent of associated spatial nonuniformities. By simulating three-dimensional arrangements, the goal is to fold in field distortion and screening effects, and evaluate the net electric field operative at the various locations and sites. This presents a first step toward including realistic spatial variations within heterogeneous tissues and organs for predicting bio-responses driven by external electric pulse trains.
Most simulation of electric field driven bioeffects in bio-medical engineering has relied on spherical cellular geometries for simplicity. Yet most cells tend to have complex shapes, and the geometries of organelles (e.g., the endoplasmic reticulum or mitochondria) are even more complex. In biology, spherical shapes are uncommon with protoplasts, murine myeloma cells,49 and some bacteria such as Streptococcus50 being rare examples of near-symmetric structures. However, most other cells are nonspherical and complex, although some cells such as mammalian red blood cells, retina photoreceptor cells,51 and many bacteria such as E. coli, Pseudomonas,50 and yeasts52 do exhibit rotational symmetry. This contribution models the buildup of the mitochondria with its two membranes, including geometric features of the cristae, within a spherical Jurkat as an example of a complex biological shape.
Toward the abovementioned goals, COMSOL-based, three-dimensional, time-dependent simulations have been carried out to analyze the spatiotemporal development of the electric fields and resulting membrane EP. Changes in the membrane conductivity and the dynamic evolution of pore radii are taken into account. The predicted response has been obtained and discussed for two different trapezoidal excitation pulses: one with fast rise- and fall-times and short overall duration, and the other for a longer-lasting trapezoidal shape. Two separate durations were chosen because transient features have been known to differentially modulate mitochondria and their viability.53 For example, previous studies evaluated mitochondrion membrane effects using tetramethylrhodamine ethyl ester (TMRE) to determine mitochondrion membrane potentials (ΔΨm). Single pulses with short rise and fall-times caused an electric field-dependent increase in calcium influx, dissipation of ΔΨm, and cell death. Pulses with long rise and fall-times exhibited electric field-dependent increases in calcium influx but diminished effects on the dissipation of ΔΨm (due to membrane poration and discharging) and cell viability. Their data seemed to suggest that high-frequency components would have a significant differential impact on mitochondrion membranes.
Here, we probe and discuss two separate scenarios: (i) single cell electroporation with particular emphasis on the mitochondria and their characteristic structures including the OMM, IMM, and cristae. (ii) Clusters of Bacillus bacteria for the quantitative predictions of EP in various arrangements, with attention to possible screening of fields and differential outcomes within the volume distribution based on the location and proximity of the periphery. An outline of the simulation scheme and other pertinent details are provided in Sec. II, which details the coupling between the time-varying electric fields at each discretized point and evolution of the membrane electrical parameters, such as the conductivity, upon EP. A Smoluchowski-based approach is utilized for assessing the pore generation and their growth in radial-space. Section III reports on the results obtained and pertinent discussions on the predicted behaviors. Finally, the conclusions and future research directions are provided in Sec. IV.
II. MODELING DETAILS
Simulations of electroporation as a way to theoretically assess directed energy bioeffects have been carried out extensively. The simplest starting point for understanding electric field interactions with biological cells is to calculate the transmembrane potential (TMP). Most work in this regard has been based on the simplifying assumption of spherical cells. Initial reports on TMP evaluation were based on a plasma membrane whose conductivity was much less than that of the surrounding medium.54 Subsequent work generalized these TMP calculations to account for general membrane conductivity,55 concentric spheres to represent the cell with an internal organelle (e.g., the nucleus),56,57 spheroidal and ellipsoidal cells,58,59 and irregular shapes.60,61 Applying electric pulses of sufficient duration and intensity to biological cells induces a sufficiently high TMP to permeabilize the cell membrane through the process of electroporation.62 The formation of pores dynamically changes the conductivity of the cell membrane both during63 and long after64 stimulation. This requires self-consistent solutions for both the membrane conductivity and electroporation dynamics.65 Hu and Joshi subsequently applied this self-consistent approach for linking TMP with electroporation dynamics for assessing spheroidal cells exposed to nanosecond duration pulsed electric fields (PEFs).66 Distributed circuit approaches that represent a cell, its membrane, and surrounding aqueous regions, with the inclusion of the formation and growth of membrane pores based on the Smoluchowski equation, have been used within a time-domain nodal method.67–69
The challenges in using home-grown codes or non-commercial software involve the computational burdens arising from complex geometries inherent in biology, or the need to adequately treat the features of disparate length scales (e.g., accurately accounting for complicated organelles, such as the endoplasmic reticulum, to tissues and organs), and variability in shapes and sizes due to heterogeneous environments. As an example of the latter aspect, bacteria70,71 can have complicated geometries, with rod-shaped bacteria growing in different directions.72 For instance, E. coli divides perpendicularly and grows in length while Thiosymbion, which is similar genetically, divides asynchronously while growing in width.
A general goal in the present context is to build a general scheme that is flexible and easily amenable to changes in geometric shapes, electrode positions, and cell populations with acceptable accuracy, while allowing user-specified additions in the cell number and degree of heterogeneity. In this context, using a scripted or home-grown solution would invariably require in-depth numerical derivations and/or re-formulations to resolve meshing and field solutions for each individual change. The result would be a loss in scalability and reproducibility for the system and require a substantial effort from changes in simulation scenarios. As a potential solution, we use the COMSOL Multiphysics suite in the present case, as a simple and efficient way to alleviate the above issues. Taking advantage of the generalized field solvers and dynamic meshing that is available in COMSOL allows for quick and easy changes in the cell geometry or size, or variations in multi-electrode positions, without placing heavy computational re-adjustment costs and burden on the user. The generalized field solver would also allow for better definition of the geometry in future modeling, including the use of piecewise continuous data gathered from the images of biological entities. This unshackles the restraint of adhering to simple shapes such as spheres or ellipsoids, which usually fall short of modeling true biological entities. The dynamic meshing also further removes computational burden from the user as the derivation of mesh sizes and spacings is done automatically within COMSOL.
In the above, N0 represents the equilibrium pore density at a TMP [ = V(t)] value of zero, while D denotes the pore radius diffusion coefficient taken to be 5 × 10−14 m2/s, while rH and rt are constants for advection velocity defined in Ref. 79. In addition, α (= 109 m−2 s−1) is the creation rate coefficient, q = (rm/r*)2 with r* = 0.51 × 10−9 m, rm = 0.8 × 10−9 m, β (=1.4 × 10−19 J) is the steric repulsion energy, Fmaximum is the maximum electric force for V(t) = 1 V, and γ is the energy per unit length of the pore perimeter previously mentioned in Eq. (2), while the effective membrane tension σeff is given as σeff = 2σ′−(2σ′−σ0)/(1−Ap/A)2. In the last term, σ′ (=2 × 10−2 J/m2) denotes the tension of the hydrocarbon–water interface and σ0 (=10−6 J/m2) is the tension of the lipid bilayer without pores. Thus, by using the coupled set of Eq. (3), the pore density N(t) at a single effective pore of radius [=ρ(θ,φ,t)] can be computed at every time step and angular location θ and φ. Using the above formulation based on the asymptotic model significantly reduces the calculations since the full pore radius distribution is no longer taken into account.
It may also be mentioned for completeness that other works utilizing COMSOL Multiphysics tools for EP have assumed a fixed pore radius of 0.76 nm. Such a fixed pore radius has then been used to directly evaluate the membrane conductivity based on an analytic expression.81,82 This radius corresponds to the minima of the pore formation energy and should ideally be the average pore radius after an exposure as the pore formation energy will pull the distribution back to the local energy minima. This method, however, would perhaps be more appropriate for quasi-static conditions, or to probe the responses to small signal excitations containing a small DC component, or if heating is well controlled. In general, having a fixed pore radius ignores the full physics and would discount segments or regions having larger electric field exposures, which may be critical for the stability of the biosystem. In any case, tests of the three methods {i.e., a full pore radius dependent distribution [n(ρ,t)], a single effective time-dependent variable pore radius, and a fixed pore radius} for a single spherical cell based on invoking azimuthal symmetry to essentially yield variability in the radial distance (r), polar angle (θ), and time (t) were carried out as an illustrative exercise. The goal was to compare the time-evolution of the TMP at the pole (θ = 0) of a spherical cell. The result is shown in Fig. 1 where a single 100 ns trapezoidal pulse with a 10 ns rise- and the fall-times is applied. All three models achieve similar peak TMPs, but in the presence of the full pore distribution, the fall in the TMP is much more drastic as the exposure continues and pores grow toward larger radii. The asymptotic model with the variable pore radius is shown to be able to approximate this result closely. With a fixed effective radius, the TMP is predicted to be unphysically locked into a high value exceeding 1.2 V.
The implementation of the variable pore radius with the simpler N(t) asymptotic formulation was adopted here. Calculations of the temporal evolution of the membrane conductivity over each angular segment involved a multi-step process. First, the electric fields at each time step for all the discretized angular positions on the membrane were computed from the COMSOL-based fine-element method. The electric fields yielded the TMPs [ = V(θ,t)] at the various membrane segments. These values were used to obtain time-dependent changes in the pore density N(θ,t) and the effective pore radius ρ(θ,t) as given by Eq. (3). The pore density and pore size were then used to obtain information on the fractional poration at each segment F(θ,t) [ = π ρ2(θ,t) N(θ,t)]/[2π rc2 sin(θ)Δθ] at every time step. This enabled computations of the dynamic membrane conductivity σ(θ,t) as σ(θ,t) = σmem [1−F(θ,t)] + σfl F(θ,t), with σmem denoting the conductivity of the unporated membrane and σfl denoting the conductivity of the extracellular fluid that could fill the pore space within the membrane at each segment. The starting values of the membrane and extracellular fluid used here are given in Tables I and II. Dynamic evolution of the membrane permittivity was obtained in an analogous manner.
Geometry and material properties for COMSOL Bacillus model.
Description/symbol . | Value . | Reference . |
---|---|---|
Cylindrical radius | 0.5 μm | 83 |
Cylindrical height | 3.0 μm | 83 |
Semi-spherical cap radius | 0.5 μm | 83 |
Bacillus dimensions | 4.5 × 0.5 × 0.5 μm3 | 83 |
Cell spacing (multi-cell models) | 10 μm (multidimensional) | Based on simulated cases |
Outer membrane thickness | 7 nm | 84 |
CM conductivity | 1.1 × 10−7 S/m | 84 |
CM permittivity | 5 ε0 | 85 |
Extracellular and intercellular permittivity | 72.3 ε0 | 85 |
Intracellular permittivity | 67 ε0 | 86 |
Extracellular conductivity | 1.2 S/m | 87 |
Intracellular conductivity | 0.3 S/m | 87 |
Electrode conductivity | 4.1 × 107 S/m | 88 |
Electrode permittivity | 6.9 ε0 | 88 |
Description/symbol . | Value . | Reference . |
---|---|---|
Cylindrical radius | 0.5 μm | 83 |
Cylindrical height | 3.0 μm | 83 |
Semi-spherical cap radius | 0.5 μm | 83 |
Bacillus dimensions | 4.5 × 0.5 × 0.5 μm3 | 83 |
Cell spacing (multi-cell models) | 10 μm (multidimensional) | Based on simulated cases |
Outer membrane thickness | 7 nm | 84 |
CM conductivity | 1.1 × 10−7 S/m | 84 |
CM permittivity | 5 ε0 | 85 |
Extracellular and intercellular permittivity | 72.3 ε0 | 85 |
Intracellular permittivity | 67 ε0 | 86 |
Extracellular conductivity | 1.2 S/m | 87 |
Intracellular conductivity | 0.3 S/m | 87 |
Electrode conductivity | 4.1 × 107 S/m | 88 |
Electrode permittivity | 6.9 ε0 | 88 |
Jurkat cell model parameters with the inclusion of the mitochondria.
Description/symbol . | Value . | Reference . |
---|---|---|
Cell radius | 5 μm | 89 |
Mitochondrion dimensions | 2 × 0.75 × 0.75 μm3 | 90 |
OMM–IMM separation width | 8 nm | 90 |
OMM thickness | 7 nm | 91 |
IMM thickness | 5 nm | 91 |
OMM conductivity | 10−4 S/m | 91 |
IMM conductivity | 5 × 10−6 S/m | 91 |
CM conductivity | 1.1 × 10−7 S/m | 84 |
OMM permittivity | 5 ε0 | 91 |
IMM permittivity | 5 ε0 | 91 |
CM permittivity | 5 ε0 | 85 |
Extracellular and intercellular permittivity | 72.3 ε0 | 91 |
Extracellular conductivity | 1.2 S/m | 91 |
Intracellular conductivity | 0.5 S/m | 91 |
Crista radius | 14 nm | 90 |
Description/symbol . | Value . | Reference . |
---|---|---|
Cell radius | 5 μm | 89 |
Mitochondrion dimensions | 2 × 0.75 × 0.75 μm3 | 90 |
OMM–IMM separation width | 8 nm | 90 |
OMM thickness | 7 nm | 91 |
IMM thickness | 5 nm | 91 |
OMM conductivity | 10−4 S/m | 91 |
IMM conductivity | 5 × 10−6 S/m | 91 |
CM conductivity | 1.1 × 10−7 S/m | 84 |
OMM permittivity | 5 ε0 | 91 |
IMM permittivity | 5 ε0 | 91 |
CM permittivity | 5 ε0 | 85 |
Extracellular and intercellular permittivity | 72.3 ε0 | 91 |
Extracellular conductivity | 1.2 S/m | 91 |
Intracellular conductivity | 0.5 S/m | 91 |
Crista radius | 14 nm | 90 |
With regard to the geometry, a separate route was taken for each of the two examples considered here. The rod-shaped bacterium in our cluster modeling was taken as a three-dimensional (3D) cylinder with spherical caps on either end. The second example simulated in this contribution required the geometric construction of mitochondria with its many cristae at the IMM. This geometry was chosen to illustrate and underscore the ability to model more complex structures that previously had not been studied numerically. Initially, a ring of eight cristae was generated around the small axis of the ellipsoidal organelle. After this, a random generation of identical cristae are moved throughout the organelle to simulate the geometries present in real mitochondria. A list of the material properties and parameters for Bacillus and Jurkat cells containing mitochondria used in our simulations is provided in Tables I and II.
The simulations were based on a three-dimensional, time-dependent analysis with the electric field distributions being solved on the basis of the finite-element method inherently available within the COMSOL Multiphysics software tool. All the domains and cell geometries were meshed using triangular elements with a much smaller mesh refinement of 2.5 nm around the cell membrane. The biological entities under study (whether a single Jurkat or Bacillus cell, or such cell clusters) were taken to consist of cytoplasm, a cell plasma membrane having a time-dependent conductivity, with a surrounding extracellular medium in a rectangular box with six faces. Voltage values, as appropriate for the chosen waveform, were applied to the top and bottom boundaries, with Neumann boundary conditions operative at the four remaining surfaces. The COMSOL tool allows for flexibility of the boundary conditions, thus making it possible to implement a variety of electrode geometries, including needle-like arrays. The dimensions of the overall rectangular simulation volume were changed based on the number of cells used in each simulation case as needed. Instead of constructing a finite membrane thickness for the bacillus or for the plasma membranes or the OMM and IMMs, these regions were treated as distributed impedances within the electric current solver within the COMSOL package. This allowed the definition of a thin surface without having to specifically construct a membrane. This would significantly reduce the computational weight of the system as it no longer required meshing down to the nanometer-levels corresponding to the membrane thickness. Our procedure and implementation were, thus, similar to other COMSOL-based approaches that have been reported in the literature.61,92–96
Finally, in our all-COMSOL implementation for EP, the decision to choose the electric field was left to the user. Thus, either plate or needle electrodes could be selected, although other multi-array options could have also easily been constructed. Details of the pulse shape and applied fields are easily handled via COMSOL's functional interface. To solve for the pore density and pore radius [as given in Eq. (3)], the Mathematics Toolbox in COMSOL was used. This methodology allows one to apply weak form boundary partial differential equations (PDEs) to the cells61 to match Eqs. (3a) and (3b). In so doing, the treatment of the cell membrane could be handled by obtaining the required TMP based on the voltage drop across the thin membrane surface. In this process, a self-consistent analysis was developed and implemented to simulate both electrical and porative developments at all cell membranes. After calculating a pore radius and pore density, the local conductivities were easily updated dynamically for each mesh point.
Finally, it may be mentioned that in the present analysis, the role of extracellular conductivity or imbalances between the internal and extracellular media on cell electroporation has not been probed specifically. We do not expect a large effect, and so this aspect was ignored. In addition, reports in the literature seem to present a slightly murky picture. For example, early work97–99 seemed to point to ease in poration (or its occurrence at lower voltage levels) with increasing conductivity of the extracellular medium. Molecular dynamics simulations some years ago seemed to suggest easier or faster poration for a higher conductivity medium.100 However, a report by Li et al.101 that probed the dependence of ion transport (and the inflow of propidium iodide) through porated membranes on the extracellular conductivity, based on a coupled Nernst–Planck, Smoluchowski approach, showed an inverse relation. Thus, their theory and results indicated that the transport of the ions or molecular delivery into cells would be lower at higher extracellular conductivity. More recent work by Novickij et al.102 seemed to indicate that small differences in conductivity would not alter cell viability significantly but could impact permeabilization to some extent. Another experimental study seemed to indicate an inverse relationship between extracellular conductivity and the efficiency of the transfer of molecules such as propidium iodide (PI) and bleomycin.103 The difference was hypothesized to arise from Maxwell stress leading to different deformations and cell shapes as conductivity affected the electric fields at cell membranes. Recent work by Sabri et al.104 showed that EP efficiency would be maximized as the conductivity ratio between the intra- and extracellular media approached unity. Their work also seemed to agree with the predictions by Smith et al.,105 which had indicated a current (i) throughput that depended on the extracellular conductivity g with the mathematical form: i ∼ g/(A + Bg). Thus, the model of Smith et al. would lead to larger current at higher conductivity. In any event, given the somewhat divergent results in the literature and minor effects for a dissimilar conductivity ratio, it was assumed that the effect of extracellular conductivity would be ignored.
III. SIMULATION RESULTS AND DISCUSSION
Simulations for a single Jurkat cell with its mitochondria were performed. The shape and geometric structure used in the calculations are shown in Fig. 2. Thus, Fig. 2 illustrates a spherical cell mimicking a Jurkat cell encompassing a mitochondrion model with its double-membrane structure, together with the cristae on the IMM. The geometry was not obtained from any digital image and so is not exact, but it does roughly mimic the general mitochondrial shapes, sizes, and features. These simulations with the mitochondria serve to highlight the implementation of our COMSOL-based model with the use of the asymptotic electroporation model for computational simplicity. The simulation runs were carried out for two input trapezoidal pulses with the long axis of the mitochondria taken to be normal to the axis of the electric field. This was done to ensure that there would be a crista subject to a direct normal exposure to the applied electric field. This orientation is likely to lend itself to high field enhancements at the local regions of the cristae, leading to EP. Two trapezoidal pulses are used: one having a shorter exposure (SE) and the other with a longer exposure (LE). The primary goal of these simulations for the two pulsed cases was to probe how changes in the rise time of the external field might affect the organelles inside the cell. As already mentioned, previous experimental reports by Beebe et al.53 showed a pulse rise-time differential in the modulation of the mitochondrial response.
(a) Mitochondrion geometry for a single Jurkat cell containing mitochondria and (b) a closeup of the mitochondria themselves with their cristae.
(a) Mitochondrion geometry for a single Jurkat cell containing mitochondria and (b) a closeup of the mitochondria themselves with their cristae.
The SE pulse was taken to consist of a total duration of 100 ns, made up of a 20 ns rise time, a 20 ns fall-time, and a 60 ns ON time at a constant field of 15 kV/cm. The LE pulse, on the other hand, was taken to have a total duration of 1 μs, with 100 ns rise- and fall-times. Figures 3(a)–3(d) illustrate the results obtained from the 100 ns short exposure (SE). Both the pore density and TMP distributions at the cristae over a region directly normal to the electric field have been shown. Thus, Figs. 3(a) and 3(b) are snapshots of the transmembrane voltage profiles at the 10 and 80 ns instants, respectively, and capture the conditions over the rising and falling edges of the applied pulse. As evident from Fig. 3(a), a strong TMP buildup occurs over times as short as 10 ns, with voltage values at the crista tip almost reaching 2 V. This voltage at the crista tip is predicted to fall to about 0.8 V at the 80 ns instant due to the creation of high density of pores leading to an electrical shorting effect. A comparison between Figs. 3(a) and 3(b) also shows the TMP gradually moving away from the tip of the cristae in a downward direction toward the base. The quick TMP rise at the tip is shown to quickly form pores at 10 ns with a peak density of 6 × 1014 m−2 in Fig. 3(c). Despite the decrease in the TMP at the tip in Fig. 3(b), the pore density, however, continues to remain strong and reaches a higher value of ∼3 × 1017 m−2 in Fig. 3(d) at 80 ns. This is simply the result of a continuous and cumulative ongoing electroporation process in the presence of the applied electrical pulse.
Simulation results obtained for the TMP and the pore density distributions across a region of the cristae normal to the applied electric field for a short exposure (SE) having a total duration of 100 ns. Snapshots of the TMP at (a) 10 and (b) 80 ns. The predicted pore density distributions (in m−2) at the (c) 10 and (d) 80 ns instants.
Simulation results obtained for the TMP and the pore density distributions across a region of the cristae normal to the applied electric field for a short exposure (SE) having a total duration of 100 ns. Snapshots of the TMP at (a) 10 and (b) 80 ns. The predicted pore density distributions (in m−2) at the (c) 10 and (d) 80 ns instants.
The resulting rise and fall of the TMPs seen in the SE results do not carry over directly to the longer exposure. In the case of the longer pulse, as is shown in Fig. 4, the average TMP across the cristae is predicted to be lower than for the SE case. This is a consequence of the longer rise time. The resulting peak pore density formation is then lower for the LE case. For example, in Fig. 4(d) at 800 ns, which is during the pulse termination phase, pore densities of at most 1.5 × 1016 m−2 are predicted, as compared to about ∼3 × 1017 m−2 in Fig. 3(d) at the falling edge of the SE case. This result is qualitatively in agreement with the observations by Beebe et al.,53 who noted stronger EP effects in mitochondria upon quicker stimulation with shorter pulses.
Simulation results obtained for the TMP and the pore density distributions across a region of the cristae normal to the applied electric field for a longer exposure (LE) having a total duration of 1000 ns. Snapshots of the TMP at (a) 100 and (b) 800 ns. The predicted pore density distributions (in m−2) at two different time instants of (c) 100 and (d) 800 ns.
Simulation results obtained for the TMP and the pore density distributions across a region of the cristae normal to the applied electric field for a longer exposure (LE) having a total duration of 1000 ns. Snapshots of the TMP at (a) 100 and (b) 800 ns. The predicted pore density distributions (in m−2) at two different time instants of (c) 100 and (d) 800 ns.
Figure 5 shows an alternate depiction and illustration on how changes in pulse rise-time and duration affect the mitochondrial response. This figure illustrates the transmembrane potential at different points normal to the electric field at the cristae, the OMM, and the cell membrane (CM). The results show a much slower growth in TMP for the OMM compared to other membrane structures for both pulsing scenarios. This behavior is directly attributed to the OMM conductivity of this membrane being orders of magnitude larger than that of the cristae and the CM. The higher conductivity mitigates growth in potential across the membrane barrier and slows the development of TMP. For shorter exposures, Fig. 5(a) shows that the cristae undergo a significantly quicker TMP rise compared to both the OMM and CM. This is caused by their low conductivity and the field enhancement associated with their curved geometry. As a result, a very sharp rise and an immediate drop occur in the cristae for shorter pulses. Also, the TMP is predicted to reach levels higher than those of the longer pulses. The immediate drop in TMP after around 120 ns in Fig. 5(a) is due to EP, which creates local electrical shorting across the crista membrane. In contrast, the TMP at the OMM is very slow to rise, not crossing 0.65 V for the shorter duration, faster rise-time pulse. Since several hundreds of milli-volts is a typically accepted value for EP initiation, our results imply that it would be difficult to porate the OMM with electrical pulses having short duration and fast rise- and fall-times. This illustrates a situation where although one might expect a buildup in TMP, there may be negligible effects in terms of induced permeability or even transport through the OMM.
TMP trends for a short exposure (SE) of total duration 100 ns and a long exposure (LE) of 1 μs duration. The panels show details of the TMP growth at the tip of the cristae, the outer mitochondrial membrane (OMM), and the cell membrane (CM).
TMP trends for a short exposure (SE) of total duration 100 ns and a long exposure (LE) of 1 μs duration. The panels show details of the TMP growth at the tip of the cristae, the outer mitochondrial membrane (OMM), and the cell membrane (CM).
While the mitochondrion modeling presented here suggests that the EP will most likely occur at the cristae and be driven quickly by short (nanosecond range) electrical pulses, the resulting transport may or may not be appreciable for two reasons: (a) first, since the tips of the cristae where most of the EP is predicted have small area, this route would present much smaller pathways (or doorways) for transport and (b) furthermore, since the EP depends on the normal component of the local fields, the effect will not be uniform and can be expected to be based on the orientation of the cristae with respect to the applied field. It is even possible that not all crista sites would be porated. In this connection then, it may practically be more useful to have multi-pronged electrodes or even a phased-array system with applied electric fields rotating across the set of electrical prongs for a better sweep of the entire spatial volume. Such novel electrode engineering was previously suggested in the context of enhanced plasmid DNA delivery.48 For completeness, it may be mentioned that the composition of the cristae and the inner mitochondrial membrane differs from the outer plasma membrane. For instance, the IMM has fewer lipids and more proteins. Hence, a more thorough analysis with regard to the parameters (such as the conductivity), possible spatial variation, and dynamic response to externally applied electric fields would be necessary. Molecular dynamics simulations are one approach for incorporating such microscopic details.106,107 However, to the best of our knowledge, such techniques have not been applied to this area of electrical stimulation and external voltage-driven bioeffects. It remains an interesting area of multiscale work in the future.
Having presented and discussed results for a single cell, with an emphasis on organelles and length scales as short at a few nanometers, we next turn to the analysis of cell clusters where collective effects may become important and begin to dominate bio-outcomes. For simplicity, this study examines the simple Bacillus (a Gram-positive, rod-shaped bacterium), which is a relevant example to study electrically induced sterilization and de-contamination in various settings based on directed energy. In addition, the rod-like geometry presents a convenient opportunity to move beyond single cell studies and probe effects in clusters, while also including a different nonspherical shape as a test of the present COMSOL-based approach.
Figure 6 shows the predicted EP behavior in a single cell. The ability to porate a cylindrical entity with changes in orientation is an important measure for gauging the overall success of an electric field driven bacteria-killing. Toward this goal, simulations were carried out for different orientations, and four representative results of the ensuing EP are given in Fig. 6. A single trapezoidal pulse with a total duration of 560 ns was used, with rise- and fall-times of 5 ns and an ON time of 550 ns. For convenience, the pore density is shown on logarithmic plots to scale the large disparities in the population range. Four different angles of 0°, 30°, 60°, and 90°, between the long axis of the Bacillus and the applied field direction, were chosen in Figs. 6(a)–6(d). As might be expected, in each case, the strongest EP is predicted to occur at the regions subjected to large normal field components, specifically at the poles for the 0° orientation, and along the sides of the cylindrical surface for the 90° angle. The maximum pore density is predicted to be roughly 1015 m−2. The simulations predict much lower membrane poration for the Bacillus compared to the Jurkat cell discussed above.
Logarithmic pore density formations of Bacillus at orientations of (a) 0°, (b) 30°, (c) 60°, and (d) 90° for the angle of the long axis relative to the direction of the applied longitudinal field. A trapezoidal 560 ns pulse, as discussed in the text, was used.
Logarithmic pore density formations of Bacillus at orientations of (a) 0°, (b) 30°, (c) 60°, and (d) 90° for the angle of the long axis relative to the direction of the applied longitudinal field. A trapezoidal 560 ns pulse, as discussed in the text, was used.
The results of Fig. 6 also reveal that as the angle between the long axis of the bacteria and the electric field increases, a definite relocation of the poration to different points along the surface of the Bacillus begins to occur. Second, the resulting pore population then [e.g., as illustrated in Fig. 6(d)] becomes smaller. This points to a potential issue when beginning to look at cell ensembles and clusters. For example, in clusters or tissues, certain cells that did not have the optimal orientation might not be porated as effectively. As a step in this direction of multicellular analysis, a collection of 125 cells with random orientations were simulated in a cubic volume with a 15 μm spacing between electrodes, as shown in Fig. 7. The ensemble was exposed to the 560 ns trapezoidal pulse discussed previously. Looking only at a single layer containing a 5 × 5 bacteria array, the resulting plot showing EP at the various cells for the 150 ns time instant is given in Fig. 8 with a logarithmic scale used again to conveniently scale the disparate population range. Figure 8 reveals varying degrees of electroporation within the bacteria population. Some of this arises from the orientation of the rod-like cell, while the distance from the electrode vicinity also plays a role. Thus, for example, a cell with a large angular separation from the applied field axis and located at the periphery has electroporation that is orders of magnitude smaller than the highest values.
Transparent down view of the 5 × 5 × 5 array of randomly oriented cells with electrodes separated by 15 μm.
Transparent down view of the 5 × 5 × 5 array of randomly oriented cells with electrodes separated by 15 μm.
Logarithm of the predicted pore density distributions (in m−2) for a 5 × 5 arrangement of bacillus cells with random orientation. A 150 ns snapshot is shown. The logarithmic scale used to conveniently scale the disparate population range.
Logarithm of the predicted pore density distributions (in m−2) for a 5 × 5 arrangement of bacillus cells with random orientation. A 150 ns snapshot is shown. The logarithmic scale used to conveniently scale the disparate population range.
The consequences of spatial distribution, orientation, and electrode separation were again probed by scaling the system to an even higher population of 1000 bacteria cells. Figure 9 shows the simulation results for a 450 ns snapshot demonstrating an EP of a 10 × 10 × 10 array with all cells oriented horizontally and exposed to the same electric pulse as chosen in Fig. 8. As evident from Fig. 9, the greatest EP takes place at the topmost and bottommost layers. Cells within these layers with positions close to the electrode are subject to the strongest EP. Although flat electrodes were used here, it is very likely that with pin/needle electrodes, the variation and spatial drop-off in electric fields with distance from the sources (i.e., an attenuation) would be even stronger. This is a point that would need to be addressed through careful electrode engineering. Furthermore, Fig. 9 shows that a large fraction of the cells in the middle region and at the periphery remain nearly intact with no EP. Finally, to better clarify this aspect of a nonuniform EP in a cluster, the 450-ns snapshot of the top three layers of the previous arrangement is shown in Fig. 10. Again, only the topmost layer exhibits significant EP, with a quick quenching going deeper into the cluster.
A 450 ns snapshot of a 10 × 10 × 10 cell arrays of bacillus illustrating pore densities at the various Bacillus cells in the cluster. For convenience, the logarithm of the pore population was taken to suitably scale the wide variation in population over the volume.
A 450 ns snapshot of a 10 × 10 × 10 cell arrays of bacillus illustrating pore densities at the various Bacillus cells in the cluster. For convenience, the logarithm of the pore population was taken to suitably scale the wide variation in population over the volume.
A 450 ns snapshot showing results for the pore density (on the logarithmic scale) across the top three layers of a 10 × 10 × 10 array of Bacillus cells with horizontal orientation. The applied field was along the vertical direction and, hence, orthogonal to the long axis of the individual Bacillus cells.
A 450 ns snapshot showing results for the pore density (on the logarithmic scale) across the top three layers of a 10 × 10 × 10 array of Bacillus cells with horizontal orientation. The applied field was along the vertical direction and, hence, orthogonal to the long axis of the individual Bacillus cells.
The primary goals of this contribution have, thus, included (a) presentation of a numerical framework with capabilities for analyzing electrical drivers that could trigger bio-responses, (b) development of numerical procedures to probe realistic biological entities of interest including cells of various shapes, clusters and colonies, and even tissues, and (c) laying the groundwork for predictions over a large spatial range of sizes starting from nanoscale dimensions such as membranes, to orders of magnitude larger entities such as clusters of cells and even heterogeneous tissues.
The framework presented here could, for example, be extended to other complicated organelles such as the endoplasmic reticulum (ER). The spatiotemporal evolution of the transmembrane potential (TMP) at the ER membrane could then be analyzed. This capability could be quite useful, especially in its application to predictions of calcium, a known secondary messenger. It could then lead to the numerical evaluation of calcium dynamics, which are known to play a role in cell death following electrical pulsing. Furthermore, it must be stated that the above goals are modest starting points and represent the first steps toward more comprehensive capabilities that could include mass transport (such as protons, hydroxides, sodium, potassium, calcium, or chloride ions), field-assisted changes in reaction rates, analysis of downstream events such as cytochrome c release and its outflow triggered by field effects or membrane peroxidation initiated by formation of reactive oxygen species, or even mechanical stresses developed from Maxwell tensor arising from inhomogeneous electric fields and/or permittivities. Some discussion of these aspects has already begun in the literature,108 although simple shapes or single entities had been assumed. Such analysis could even facilitate comparative assessments of parallel pathways.109
However, it must be acknowledged that the subject of mutual interactions within the mitochondria and/or other organelles is very complex and involves as yet unknown processes. Furthermore, quantitative data on possible rate coefficients and parameters are not fully known.
IV. CONCLUSIONS
A self-consistent analysis of cellular electroporation for bio-medical applications has been presented. The focus here was on two levels. The first objective was to devise a general scheme to transcend the limitation of simple cellular geometries and provide a simulation tool capable of incorporating realistic biological shapes and sizes. This was accomplished by using the COMSOL Multiphysics suite, with suitable embellishments to incorporate the details of the electroporation process and its many inherent internal physics. Second, we assessed the electric field driven biophysics for both microscopic organelles (such as the mitochondria), where spatial dimensions on the order of a few nanometers become important, and larger length scales encompassing multiple cells, where collective effects and mutual interactions can dominate. Thus, this approach provides scalable computing to generalized geometries with the ability to analyze different complex microscopic organelles. Although we chose EP here as a simple, representative biological effect, the technique is general and could be used to analyze any electrically triggered biological phenomenon, such as nerve signaling, electrically stimulated chemical release, changes in reaction kinetics due to the variation in the electrochemical potential, or proton transport at membranes. Also, although this study features a simple trapezoidal excitation as an illustrative example, different waveforms (e.g., monophasic, biphasic, combinations of different frequency components, alternating current) could also be assessed. Finally, though not discussed here, the present simulation scheme would also allow for the inclusion of thermal effects (as was reported in a separate context by our group),63 pressure changes, fluidic flows, and other stimuli.
Both Bacillus bacteria and mitochondrion organelles were chosen as test examples with widely differing length scales for the present approach. Our results for a single Jurkat cell containing a mitochondrial structure suggest that the EP at the cristae within the IMM will be likely and quickly driven by short (nanosecond range) electrical pulses. The buildup of TMP at the OMM, on the other hand, will most likely not be very rapid. Hence, either longer pulses or a wavetrain consisting of multiple pulses might be required for poration (as was used by Ref. 23 with consequent transport through the OMM). For the IMM, on the other hand, despite fast EP mainly at the crista sites, the overall transport volume may or may not necessarily be appreciable for two reasons: (a) since the tips of the cristae where most of the EP is predicted to have a small area, this route would present much smaller doorways for transport and (b) because EP depends on the normal component of the local fields, the effect will not be uniform and can be expected to be based on the orientation of the cristae with respect to the applied field. It is even possible that not all cristae would be porated. In this connection then, it may practically be more useful to have multi-pronged electrodes or even a phased-array system with applied electric fields rotating across the set of electrical prongs for a better sweep of the entire spatial volume.
Our simulations of Bacillus as simple examples of cell clusters spanning several micrometers show a distinct screening effect and the difficulty for electric field penetration into a cluster. This is a point that would need to be addressed through careful electrode engineering. The development of phased-array or asymmetric cathode–anode systems is a possibility for overcoming this shortcoming. In any event, the numerical technique described and applied here could be used for such analysis toward devising optimized electrode-based pulsing strategies.
Finally for completeness, some limitations and constraints of the present technique that could be addressed to attain even better simulation capability are briefly outlined. For example, the focus of the current treatment is on membrane electroporation based on developing strong TMPs. However, this does not necessarily guarantee tumor cell death for a variety of reasons: (a) the electric fields that could be generated in a tumor mass consisting of a collection of cells would have spatial variation. Hence, not all regions might be subjected to the necessary field values above the threshold for irreversible electroporation or field-assisted chemotherapeutic drug uptake. In this connection, the role of the proximity effect has already been discussed on this contribution. (b) Second, the cell parameters within a cluster may have a variation due to various factors including cell lifecycle and natural variations in shape and size. (c) Besides differences in the spatial distance from the applied excitation source site, heterogeneous details of the biological tissues would also need to be carefully addressed. For example, the presence of vascular structures with blood has been shown to lead to the redistribution of electric fields.110 The reductions in the electric field strength were found to be in the proximity to large blood vessels and clustered vessel structures. Consequently, such cells would likely be less affected by externally driven field treatments due to their location near larger vessel structures. (d) From a practical standpoint, it would also be useful to acquire better knowledge of possible cell repair109 and various other mechanisms and factors that might reduce repair for greater tumor cell killing. (e) Finally, from a pure simulation standpoint, it may be possible to improve the efficiency, run time, and memory requirements of the present numerical scheme that then allow the inclusion of realistic anatomical details for actual treatment planning. Such approaches have already begun to emerge,111 and the current effort should move toward similar goals.
ACKNOWLEDGMENTS
The support from the U.S. Office of Naval Research under Grant No. N00014-21-2055 is gratefully acknowledged. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for government purposes notwithstanding any copyright notation herein. In addition, R.P.J. acknowledges useful discussion with T. Batista Napotnik (University of Ljubljana).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
W. Milestone: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). C. Baker: Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal). A. L. Garner: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). R. P. Joshi: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.