Sinusoidal RF Simulations for Optimized Electroporation Protocols

Protocols surrounding electroporation have long been based on trapezoidal pulsing of biological cells. Here, we revisit cellular electroporation for bio-medical applications, including tumor treatment, based on a self-consistent electro-thermal analysis with sinusoidal RF excitation. Predictions for the evolution of pores and their surface angular distribution, as well as potential heating and temperature increases, are given. Our results show an optimum frequency range from 5–7 MHz to achieve increased mass transport without detrimental heating in Jurkat cells. Through parametrized frequency sweeps, this work establishes potential optimized regimes that could guide experimental and clinical protocols. More signi�cantly, the optimal frequency for porating healthy B-cells is predicted to be ~ 2.5 MHz, with almost no poration at 7 MHz. This opens up the exciting possibility for treating malignant tissue with a well-tuned optimal frequency range for bioeffects, while minimizing deleterious effects on healthy cells and tissues.

Such electric stimulation opens up two slightly different routes to tumor treatment. If the intensity of the externally applied electric eld is not very high, then the effects are transient, and the membrane pores can reseal and the cell can completely recover. This scenario is utilized in electrochemotherapy, which combines external electric elds with chemotherapeutic drugs (e.g., bleomycin or cisplatin) for targeted uptake leading to tumor cell death [22,23]. On the other hand, applying higher intensity or longer duration electric pulses, termed irreversible electroporation (IRE), has proved to be an effective method for solid tumor treatments. This route leads to loss of homeostasis and has been used for non-thermal ablation therapy [24][25][26][27][28] with important advantages over conventional ablation methods [29]. In IRE ablation, cell death predominantly results from electroporation rather than local temperature increases, though some degree of local heating of tissues can occur, especially when many pulses are applied [30,31]. Compared to electrochemotherapy, IRE tumor ablation uses more pulses and higher voltages.
Further research has shown IRE to be a successful noninvasive method for cell ablation with low risk to other nearby organs, including nerves and blood vessels. In this regard, high frequency irreversible electroporation (H-FIRE) was proposed as a potential strategy [32][33][34][35][36] based on the delivery of bursts of short bipolar pulses with durations from 1-5 µs rather than the longer duration 70-100 µs monopolar pulses typically used in traditional IRE. This novel approach has been successfully applied to treat solid tumors with minimal muscle contractions [32,[35][36][37][38][39]. The electric elds required to kill cells with H-FIRE are signi cantly higher than in conventional IRE protocols [40]. The biphasic H-FIRE protocol combines pulses having ON-times in the 0.5-2.0 µs range, separated by OFF-times in the 0.5-5.0 µs range [32]. These pulses are repeated for overall durations on the order of hundreds of microseconds. Tests in canines [41] achieved signi cant ablation of brain tumors with minimal muscle contraction.
While IRE protocols typically use microsecond duration electric pulses, early electroporation experiments used much longer, millisecond duration pulsing [42][43][44]. These longer duration pulses generated pores much larger than the size of the bleomycin and/or sucrose molecules with slow pore resealing kinetics.
For example, pore resealing took ~ 15-20 min after a 2-ms duration pulse [45]. On the other hand, Schoenbach et al. proposed electrically induced tumor cell killing based on much shorter nanosecond pulses [46][47][48]. This strategy used a much higher eld intensity (~ 50-100 kV/cm) compared to 0.5-5 kV/cm in the H-FIRE scheme [32][33][34][35][36]. Numerical modeling predicted much higher pore densities, but with lower nanometer sizes [49][50][51]. The presence of nanopores was also con rmed by experimental measurements [52]. These short pulse durations were chosen to be shorter than the charging time of the plasma membrane but on the order of the charging time of intracellular organelles, thereby allowing for poration of the inner membranes of organelles, such as the nucleus, mitochondria, or the endoplasmic reticulum, before fully permeabilizing the plasma membrane [53,54]. Electrically modulated cellular responses in this mode have included intracellular calcium release [55,56], shrinkage of tumors [57,58], temporary blockage of action potential in nerves [59], and activation of platelets for accelerated wound healing [60]. Finally, as a consequence of the higher elds and shorter exposure time, a sharp temperature gradient at the membrane has been predicted at the start of such nanosecond excitation [61]. This could have consequences on mass and ion transport based on thermo-diffusive ows.
The primary considerations in electroporation-based therapy are the need to maintain open pores of su cient size and number for adequate throughput of drugs, exogeneous molecules and ions as needed. Pore opening and their evolution in size, are dynamic and depend on the transmembrane potentials across the surface. In particular, for pulsed operation, the frequency at which the external stimulation is applied relative to the inherent time constants of the bio-membrane system assumes an important role.
Time durations of pulsed stimulation for the electrically assisted tumor treatments have ranged from nanoseconds to milliseconds. However, logically there could likely be an optimum operating point, an aspect that has not been studied to the best of our knowledge. While there may be many aspects to consider in gauging the success of such treatment, the primary focus might be achieving adequate poration, while minimizing excessive membrane heating to avoid collateral damage. It is reasonable to expect that an optimized domain of operation in terms of operating frequency and applied bias to achieve this condition will most likely exist.
Towards this goal of optimizing operating conditions, at least in terms of excitation frequency, here we undertake a comprehensive electro-thermal simulation analysis of single cell electroporation. The method is similar to that employed in a recent report [61]. Changes in the membrane conductivity, dynamic evolution of pore radii, and changes in the membrane temperature are all evaluated for a broad range of sinusoidal biphasic excitation with operating frequencies from 1 MHz to 50 MHz. We also include the role of delays in the input signal wavetrain, since inherent charging times of the circuit preclude continuous operation and energy delivery over long times without any OFF periods. The simulations are carried out over periods spanning up to 50 µs which exceeds the temporal span of the H-FIRE burst. Since internal heating and temperature is a slow process compared to electroporation, the longer 50 µs simulation times allow for reasonable assessment of heating effects as well. Rather than rectangular pulses, a sinusoidal input was deliberately chosen to focus on a single frequency at a time, i.e., a narrowband situation as has been used in some reports [62]. Use of rectangular or other types of waveforms would have inherently added a mixture of constituent frequencies, making it di cult to unfold the response behaviors at a particular frequency.
The following sections detail the simulations, the results obtained and the relevant discussions. An outline of the simulation is provided in Section II which details the discretized circuit model, changes in surface tension of the membrane with poration, and coupling of the time-dependent heat-ow equation. Section III reports the results and pertinent analysis, including the role of the resulting transmembrane potential (TMP) waveforms, heating along the membrane, and the optimum frequency regimes to achieve successful electroporation with minimal thermal collateral damage. We provide our conclusions and future research directions in Section IV.

Ii. Modeling Details
Modeling for the single spherical cell was carried out using a coupled numerical scheme involving the (e) Finally, the possibility of oxidative degeneration of lipids due to oxidative stress [67-69] often associated with electric stimulation that can occur as a downstream event on longer times, was ignored as well.
The simulation region consisted of a sphere encompassing the cell and its surrounding extra-cellular medium. Discretization of the simulation region with spacing dr, dθ, and dφ yielded elemental volume segments. Each segment was represented by a parallel combination of a resistor and a capacitor to form a local impedance and thus emulate the elemental electrical characteristic. The parallel resistivecapacitive (R l C l ) combination accounted for both the conduction and displacement current components in the elemental space. The local resistivity, permittivity and physical dimensions were taken into account in computing values for the R l and C l elements for each volume, with "l" denoting the l th local region (i.e., either membrane, intracellular, or extra-cellular media). This essentially divided the entire region of interest into a distributed network of interconnected impedances. By applying the Kirchchoff current continuity law at each of the n interior nodes of this discretized system, a set of n-equations coupling the voltage values V ij with its six nearest neighbors was obtained. Solving for the node voltages V(r i , θ j , φ k ) within this distributed network then yielded the local electric eld drivers at each time step. Voltages (V j ) at the boundary nodes of this spherical volume were assigned values based on their distance R from the center of the cell (R being the boundary radius) and their angular position θ j (relative to the z-axis) for a chosen external eld E 0 as: V j = -E 0 R cos(θ j ). Finally, it may be mentioned that for simplicity, the cell membrane itself was taken as one integral grid-spacing unit in the radial r-direction, i.e., this sub-region was not further discretized due to its small (~ 5 nm) spatial extent. It may also be noted that the idea of a distributed circuit for electroporation analysis had also been reported some years ago by Sweeney and Davalos [70].
Results of the time-evolving TMP were then used in a continuum Smoluchowski equation [71][72][73] to analyze pore formation, growth and decay. This Smoluchowski-equation based continuum description accounts for all electropores through a size distribution without tracking the kinetics of each individual electropore. The pore formation energy E depends on the effective bilayer membrane tension σ e and the transmembrane potential ΔV(r i , θ j , φ k ). Most analysis in the literature have tended to assume a xed value for the membrane tension. Here, however, a variable tension, in keeping with the approach by Mukherjee et al. [74] and Neu et al. [75], was used in the calculations. The effective membrane tension was taken to be: 1 where σ denotes the unperturbed membrane surface tension in the absence of any pores, σ' is the energy of the hydrocarbon-water interface, and A f represents the ratio of the porated area to non-porated area.
The above Eq. (1) correctly accounts for a reduction in the membrane tension with increasing pore density (and hence, total pore area) which works to increase the barrier for pore formation in the Smoluchowski representation. Very roughly, as pores form, the membrane loss its tension through local relaxation, and thus leads to decreases in the effective value . Not only does such a pore-area dependent approach match expected experimental results [75], but it also stabilize the simulations by restricting uncontrolled pore growth and expansion.
The membrane was discretized in\to M angular segments, and the full Smoluchowski model was solved. This yielded the complete pore density n(ρ,θ,t) with its various spatial dependencies rather than using the asymptotic approximation of Neu and Krassowska [72] for n(t) alone. The variable "ρ" in the term "n(ρ,θ,t)" represents the distance in r-space for the pore radii (r) at a membrane segment with an angular position "θ". Values of the electric elds E(r,θ,t) and local currents J(r,θ,t) at all grid points in the simulation space were used to obtain the power density S k (r,θ,t) at the corresponding grid point from their dot-product as: 2 .
In the above, "k" denotes the k th component (membrane, pore, or intra-or extra-cellular aqueous medium). Thus, at each angular position, the time-dependent local membrane conductivity σ k (r,θ,t) within each elemental volume was factored in. The time-dependent heat conduction equation for calculations of the internal temperature is given as: Here, the implicit second-order Crank-Nicholson technique [76] was utilized based on its numerical stability. All parameters, such as ρ vk , C k and k ck , were assumed to be temperature independent, given the relatively modest temperature increases expected in the present study. The parameter values used are given in Table 1 and are for Jurkat tumor cells. Iii. Simulation Results And Discussion The central goal of this work is to perform self-consistent electro-thermal analyses of cellular response to varying sinusoidal pulses to assess the possibility of an optimal operating frequency range from the standpoint of electrically assisted tumor treatments. The precise value of an optimal frequency, if found, could conceivably vary with both cell type and cell-cycle. Furthermore, identifying a desirable range of operating frequencies might help tailor the input excitation, and guide in the selection of the frequency parameter to speci cally act on tumor targets over healthy cells.
In the present simulations, the electric eld amplitude of all sinusoidal pulses was maintained 15 kV/cm. As regards the total duration, 50 µs long wavetrains were used, under two different scenarios. In the simplest implementation, the pulses were continually applied without any intermittent delays or OFF times. In an alternative implementation, a slight delay was included after each pulse to account for any practical time lag in charging or circuit timings. From a numerical standpoint, this is not a critical issue, and can easily be adjusted in accordance with actual experimental implementations that might be in use. But regardless of the details, it can be expected that without any OFF times, the predicted membrane temperature would be higher, and possibly lead to a slightly larger population of pores with greater mass transport. In their experiments, ingle frequency sinusoidal narrowband excitation was applied to Chinese Hamster Ovary (CHO) cells. The noticeably large cell death at the lower frequency excitations was attributed to greater heating.
In order to reduce temperature increases, a slight delay (or OFF time) was implemented in our pulse-train. Figure 2(a) illustrates the scheme with each single full cycle sinusoidal pulse being followed by a 500 ns OFF-time. This delay is in the range of the OFF times used in the H-FIRE protocol. Figure 2b is an expansion of Fig. 2(a) along the time-axis to better illustrate the constant 500 ns delay implemented after each cycle. Three different excitation frequencies of 1, 5 and 10 MHz were used. This selection of a common OFF-time for all three frequencies was meant to reduce the number of parameters in the system. Since the scheme is general, for precise comparisons against actual experiments, the OFF-time values could be adjusted if and as necessary. However, the qualitative trends predicted here should continue to hold.
Simulation results obtained at frequencies of 1, 5, 7, and 10 MHz are shown in Fig. 3. The values in all cases were at the polar region (i.e., an angular position of zero degrees.). For clarity, the predicted evolution of the temperature rise over the initial 0.25 µs is given in Fig. 3(a), while the response over the entire 50 µs duration is shown in Fig. 3(b). The ambient was again chosen to be 300 K. The plots of Fig. 3a shows a clear initial temperature rise (ΔT) above the ambient after about 34 ns from the start of the pulse. The temperature increase at all four operating frequencies is qualitatively similar, though the maximum value is seen to occur for 7 MHz after about 40 ns after the start of the pulse. The initial peak for ΔT at the 1 MHz excitation is the lowest at about 20 K. The maxima of the other peaks are predicted to rise above 25 K with ease. While this level of heating would be a concern if sustained for long periods of time, this rise lasts for less than 40 ns. This time scale is so small that any resulting thermal damage, which is often measured in terms of an integral associated with the Arrhenius equation [85,86], can be assumed to be negligible. More prolonged exposures occur at later time periods. During these exposures, lower temperatures occur at either optimal frequencies or frequencies that are too high. The plots of Fig. 3 An interesting feature seen in Fig. 3(b) is that at higher frequencies, the cellular system maintains a lower average temperature. This is the result of two factors. (i) One, since a constant 500 ns OFF time had been assumed in our simulations, and since the number of cycles over any given xed time duration is larger at higher frequencies, the number of OFF periods is also correspondingly more at higher operating frequencies. This allows for greater cooling within the system (ii) More importantly, there predominantly is a lack of the transmembrane potential (TMP) build up at higher frequencies. This is best understood by rst considering the situation at very high excitation frequencies. Under these conditions, the time period would be shorter than the membrane charging times. As a result, the TMP values would not rise much before a swing in the input polarity would begin to decrease the membrane voltage. Consequently, the TMP across cell walls remains poor at high frequencies, leading to poor poration, current throughput and corresponding weak thermal dissipation or temperature enhancements. However, though the membrane heating at the high frequencies is under control, the poor poration implies minimal mass throughput and weak injection of drugs to affect the tumor cell.
As already discussed, the temperature rise at lower frequencies is high since a large TMP can be created across cell membranes, though it occurs at longer times. However, since each of the individual cycles are also long (i.e., longer operating time periods), pores can be formed and remain open during the excitation. The longer the system is subjected to the input low-frequency elds, the power dissipation continues, leading to greater heating over time. In contrast, operation within a frequency regime that is in between these two extremes leads to an optimal situation. The heating though marginally stronger than that at very high frequencies, is still not large or damaging, and nor would it lead to much by way of collateral damage. However, the TMP can be sustained, leading to mass transfer and drug delivery into tumor cells.
The temperature rise at the 5, 7, and 10 MHz excitations fall off within about 0.1 µs, and thereafter, remain xed at about an incremental value below 6 K. For the 1 MHz excitation, on the other hand, much higher temperature enhancements of about 17-18 K are predicted in the long run in Fig. 3(b).
The angular dependence of the temperature increase over the cell membrane is brought into focus through the results shown in Fig. 4 for a 7 MHz excitation frequency. Though the highest electric elds are applied at and near the polar regions, the transmembrane potentials are not necessarily the highest in these areas. The large and rapid increase in electric elds over time starting time of pulse application, leads to strongest poration early on near the polar regions. This produces a large drop in the average TMP due to large increases in membrane conductivity and a consequent ″electrical shorting″ effect. Figure 4(a) illustrates the result of this behavior. An initial spike in temperature, is quicky quenched as the TMP reduces leading to sharp reductions in the local electric elds, and hence, the power dissipation density which is a product of the local current density and electric eld. A clearer picture over the longer term for the temperature rise can be seen from Fig. 4(b). The largest increments are predicted to be at the 36˚ angular position (the 0˚ position corresponds to the pole), with periodic spikes predicted to form with delays equal to the time period of 0.14 µs and inter-pulse delays of 500 ns.
Having discussed the temperature changes, the time-dependent variations in TMPs for four different excitation frequencies are discussed. The values in all cases were obtained at the polar region, i.e., an angular position of zero degrees. Over early times, shown in Fig. 5(a), the peak TMP values of almost 2 volts are predicted to be reached for the 5, 7, and 10 MHz excitations. The corresponding time periods are 0.2 µs, 0.1428 µs, and 0.1 µs, respectively. All three are lower than the timescale of Fg. 5(a), and so the peak in TMP is seen for all three cases, before the values start to fall as the polarity of the excitation signal is reversed. In the case of the highest 10 MHz case, a negative potential is also seen to be reached. However, at the slowest 1 MHz excitation, with a corresponding 1 µs time period, only the drop from the peak due to localized poration at the polar cap, is seen to occur at about 95 ns in Fig. 5(a). Thereafter, the TMP at the 1 MHz excitation begins a slow drop off and eventually assumes negative values after about 0.5 µs in Fig. 5(b). For the 5 and 7 MHz excitations, the rise and fall of the driving signal is roughly tuned to the charging time constant and pore reduction phase following the growth in poration. The spikes in the TMP are then roughly in step with the external excitation, as seen in Fig. 5(c). Furthermore, the peak TMP values at these two frequencies are predicted to be at about 0.32 Volt, which place the operating point at a local minima of the pore formation energy curve [72,73]. This is a stable operating point, and so the system can continue to maintain mass transfer through the membrane in a reliable and lasting manner without any danger of pore runaway or cell destruction.
Combining the results of the electro-thermal analysis, the frequency-dependent behavior of both the median pore radii distribution and average temperature rise are shown in Fig. 6 spanning a range from 1 to 100 MHz. As the frequency increases, a clear trend of a decreasing median pore radius is seen. This in line with the results shown and discussed in Fig. 5. At higher frequencies, lower TMP values are set up, with the operating point near a local minimum of the pore formation energy curve. On the other hand, with lower input frequencies, the cell faces prolonged exposure to higher TMPs as was shown in Fig. 5(c).
This takes the operating point beyond the minima of the pore formation energy characteristic, thus leading to bigger pores. However, the density of such pores at the lower frequencies is not large since the voltage overshoot effect is only possible at high frequencies. The dashed line shown in Fig. 6 represents the location of the pore radius corresponding to a minima in the pore energy characteristic. One can expect that successfully controlled electroporation protocols would have the largest amount of pores remaining in the stable regime. This is exactly the case when looking at the previously discussed ″ideal″ frequencies of 5 and 7 MHz. The other consequence highlighted in Fig. 6 is the decrease of average temperature with increasing excitation frequency. This trend is also in line with the results shown in Fig. 3, and also matches the data from Katsuki et al. [62] of greater thermal damage at lower operating frequency.
Next, differences in the bio-response between normal and tumor cells are presented through our simulations. Figure 7 shows a comparison of the outcomes for pore density creation between normal Bcells and Jurkat cells. Parameters for the B-cells were obtained from a report by Polevaya et al. [87] and those for Jurkats are listed in Table 1. The maximum conducting pore density values are shown in Fig. 7. The principal aspect of this result is the difference in the two cases, and the overall shift in frequency space between the healthy and tumor cells. The maximum pore formation for Jurkat cells is predicted to occur near the maximum predicted at around 7 MHz, while the maxima in Fig. 7 occurs around 2.5 MHz for the normal B-cells. This is a desirable goal and a signi cant result as it underscores the potential for selectivity. For example, applying a 7.5 MHz RF eld would be very effective in porating Jurkat cells, while the same excitation would have almost no effect on healthy cells. Obviously, the values obtained here are subject to change with cell parameters, and hence, depending the speci c cell-type, the optimal excitation frequency of operation for selective treatment would change. However, the overall existence of a desirable or favorable operating point would remain and could be used for targeted tumor killing, or cell targeting for selective membrane permeabilization. It would also be very bene cial to generate a series of data tables and plots for optimal operating frequency points for various biological cell types as a potential guide for parameter selection towards clinical protocols and treatments.

Iv. Conclusions
A self-consistent electro-thermal analysis of the cellular electroporation for bio-medical applications, including tumor treatment, has been presented. The focus has been on quantitative assessment of the cell membrane response to sinusoidal RF excitation pulse trains. Predictions for the evolution of pores and their angular distribution on the surface, as well as potential heating and temperature increases, have been given. Unlike the trapezoidal input signals studied in most previous reports which inherently fold in multiple frequency components, here the predicted responses were obtained for a range of speci c input frequencies. Such a frequency sweep allowed three potential bene ts: (i) It clearly brought out an operating window in the frequency domain from the standpoint of creating optimal cell poration and thus mass throughput, while keeping the temperature rise to within moderate levels. (ii) Additionally, by knowing the bio-response at different individual operating frequencies, one could tailor or roughly predict the expected behavior to different waveforms containing a slew of frequency components. (iii) Finally, and perhaps the most interesting and useful aspects were the predicted differences in the operating frequency between cell types at which strong pore formation could occur. Here in comparing the Jurkat tumor cell response with that of healthy B-cells, a substantial difference between frequency values of 2.5 MHz and 7.5 MHz, respectively, was obtained. Even more signi cant was the near absence of poration at 7.5 MHz in the case of B-cells, for which the Jurkats had maximal electroporation. This would, for example, portend the possibility of selectively eliminating Jurkats without harming the B-cells based on a sinusoidal RF excitation protocol at 7.5 MHz.
While other cell types might have slightly different optimal ranges for the input electrical parameters depending on their inherent properties, a maximal operating frequency for greatest poration can be expected to hold. The present simulation technique is general and so would be able to ascertain the bene cial frequency range for different cell types. The simulations also offer variability in the off-times between pulses as dictated by the driving circuit capabilities and practical consideration of the switching electronics, for assessing a variety of potential scenarios. Finally, exciting future work based on the present method could explore optimized regime and differences in bio-responses between various healthy and tumor cells, for tabulating selectivity enhancement regimes as a guide for parameter selection for actual treatments.

Con ict of Interest
The authors have no con ict to disclose.