A statistical model describing combined irreversible electroporation and electroporation-induced blood-brain barrier disruption.

BACKGROUND: Electroporation-based therapies such as electrochemotherapy (ECT) and irreversible electroporation (IRE) are emerging as promising tools for treatment of tumors. When applied to the brain, electroporation can also induce transient blood-brain-barrier (BBB) disruption in volumes extending beyond IRE, thus enabling efficient drug penetration. The main objective of this study was to develop a statistical model predicting cell death and BBB disruption induced by electroporation. This model can be used for individual treatment planning.
MATERIAL AND METHODS: Cell death and BBB disruption models were developed based on the Peleg-Fermi model in combination with numerical models of the electric field. The model calculates the electric field thresholds for cell kill and BBB disruption and describes the dependence on the number of treatment pulses. The model was validated using in vivo experimental data consisting of rats brains MRIs post electroporation treatments.
RESULTS: Linear regression analysis confirmed that the model described the IRE and BBB disruption volumes as a function of treatment pulses number (r(2) = 0.79; p < 0.008, r(2) = 0.91; p < 0.001). The results presented a strong plateau effect as the pulse number increased. The ratio between complete cell death and no cell death thresholds was relatively narrow (between 0.88-0.91) even for small numbers of pulses and depended weakly on the number of pulses. For BBB disruption, the ratio increased with the number of pulses. BBB disruption radii were on average 67% ± 11% larger than IRE volumes.
CONCLUSIONS: The statistical model can be used to describe the dependence of treatment-effects on the number of pulses independent of the experimental setup.

Introduction

Electroporation (EP) is a physical phenomenon in which electric fields make cell membranes transiently permeable to ions and macromolecules which are otherwise deprived of or have limited trans-membrane transport mechanisms.1–3 Electric pulses applied to the tissue induce an electric field which in turn induces a change in cell membrane potential. This change depends on various tissue related parameters such as tissue type and cell size as well as pulse parameters including pulse amplitude, shape, duration, number of pulses, and pulse repetition frequency. As a function of the induced electrical field, electric pulses can either: reversibly permeabilize the cell membrane (reversible EP) or permeabilize the cell membrane in a manner that leads to cell death (irreversible EP).4 It was recently demonstrated that when applying EP to brain tissue it also induces reversible disruption of the blood-brain barrier (BBB).5–7

Both irreversible EP (IRE), and reversible EP combined with chemotherapy, also known as electrochemotherapy (ECT), are emerging as new treatment techniques for solid tumors.3,8–16 ECT uses EP to allow increased uptake of chemotherapeutic drugs into tumor cells12 and IRE is a method aimed at inducing tumor ablation without thermal damage.17,18 Brain tumors are excellent candidates for local EP treatment. Glioblastoma multiforme (GBM) is the most frequent and most aggressive primary brain tumor with an average survival of 14 months from diagnosis. Existing treatments offer poor prognosis for GBM mainly due to tumor infiltration into the surrounding brain, high resistance to therapeutic apoptotic stimuli and poor BBB penetration of most therapeutic agents.19,20 A combined approach, consisting of inducing significant/rapid necrosis in the tumor mass and simultaneous delivery of high chemotherapy doses to the tumor and surrounding infiltrating zone is suggested as a treatment strategy. EP-induced tissue necrosis within the massive region of the tumor and surrounding BBB disruption, enabling efficient local delivery of systemically administered chemotherapy was recently demonstrated.7,21

Individual treatment planning is an important key for EP-based treatment success.22 Treatment parameters should be chosen in such a manner that will induce maximal damage to the tumor while sparing surrounding healthy tissue. This is usually done by numerical models. Several numerical models describing the electric filed distribution in the tissue have been introduced, and are applied for predicting treatment outcome and planning the electrodes placement to ensure full tumor coverage by electrical fields higher than the EP threshold.23–28 These models are usually based on experimental data. Treatment volumes calculated from MRI29 or histological data30,31 are incorporated into computerized models together with the organ characteristics and the electrodes configuration. These calculations traditionally use deterministic models, i.e all the cells exposed to electrical fields higher than a specific threshold, known in the literature, will be irreversibly/reversibly electroporated. Nevertheless, live tissues are more complex, especially malignant tissues which are inherently inhomogeneous, and therefore assuming a statistical effect of EP parameters maybe more appropriate.32,33 For this reason we chose to apply a statistical model to describe reversible/irreversible effects in vivo.

The Peleg-Fermi model is the most widely used mathematical model for describing cell death as a consequence of IRE in medicine.32–34 Although several other models have been proposed35 the Peleg-Fermi model seems the most adequate since it includes dependency on the number of pulses as well as in electrical field. For this reason we decided to apply it on our experimental data and further extend it to irreversible and reversible EP effects in vivo.

The Peleg-Fermi statistical model was first introduced as a model describing the survival of bacteria after exposure to pulsed electrical fields.36 Later on it was suggested that this model can be adapted to describe the effects of IRE.32,34 Goldberg and Rubinsky32 extrapolated experimental data obtained using prostate cancer cells and demonstrated the feasibility of applying this model to describe the effects of IRE for up to 10 treatment pulses. Garcia et al. extended the model up to 90 pulses, by theoretical analysis that is yet to be confirmed with experimental data.34

Treatment parameters such as pulse shape, amplitude, frequency, duration and number of pulses37,38 affect treatment outcome. Here, we chose to study and model the effect of number of pulses while other pulse parameters remain fixed.

A numerical model describing electric field distribution in the brain tissue based on the applied voltage, tissue and electrodes electrical properties and electrodes configuration was constructed. The calculated electrical field was then implemented in the statistical model that was estimating the effect of the number of pulses on the outcome- irreversible damage and BBB disruption.

The first goal of our study presented below was to extend the Peleg-Fermi model to describe a wider range of the number of treatment pulses in vivo and to validate the extended model using experimental data obtained from naïve rats treated with EP in the brain.

The second goal was to adapt the statistical Peleg-Fermi model to describe the effects of pulse parameters on BBB disruption. BBB disruption is a vital key in treating brain tumors since it is important to disrupt a large enough volume surrounding the tumor mass to enable efficient drug penetration into the infiltrating zone. Once established, models describing both IRE and BBB disruption can be implemented to provide a complete treatment planning for brain tumors with EP.

Materials and methods

Animal experiments

The study was approved by and performed in accordance with the guidelines of The Animal Care and Use Committee of the Sheba Medical Center, which is approved by the Israeli authorities for animal experimentation.

We have recently presented the results of an animal experiment designed to study both IRE and BBB disruption using the same experimental setup.6,7 Here we describe in detail the aspects relevant to our statistical model which are based on that experimental data. Our unique electrode setup employs a single insulated intracranial needle electrode with an exposed tip placed in the target tissue and an external surface electrode pressed against the skin. The electric field produced by this electrode configuration is highest at the exposed tip of the intracranial electrode tissue interface and then decays with the square of the distance. Therefore, the electric fields surrounding the needle electrode tip induce nearly spherical IRE effects at the target tissue and gradually decrease further away to reversible EP effects which induce BBB disruption. Regions of interest (ROIs) plotted on MR images acquired post EP treatments with various pulse parameters were used for calculating the tissue damage and BBB disruption radii. We then studied the correlation between the experimental radii and the extended statistical model.

Animal model and procedure

The study was performed by treating 46 male Spring Dawly rats with 50 μ s monopolar electric pulses at 1 Hz and 600 V, as previously described.7 The rats were divided into seven groups of 5–7 rats each, treated with varying number of pulses (N = 10, 45, 90, 180, 270, 450 and 540).

MR imaging

Rats were scanned 30 minutes post treatment and periodically thereafter up to 2 weeks post treatment, using a 1.5 T GE Optima MR system (Optima MR450w, General Electric, Milwaukee). The MR sequences included contrast-enhanced T1-weighted MRI for depiction of BBB disruption and T2-weighted MRI for depiction of tissue response. Gradient echo (GE) MRI was acquired to assess possible procedure-related bleeding.

The damage radius induced by IRE (rd) (in mm) for each rat was calculated from the hyper-intense regions on T2-weighted MR images acquired two weeks post treatment. This time point was previously determined by histology as adequate to describe IRE.7 BBB disruption radius (rb), referring to the maximal radius of tissue in which the BBB was breached, was calculated from enhancing regions on contrast-enhanced T1-weighted MR images acquired 30 minutes post EP treatment.

In both cases the radii were calculated by delineating ROIs over the entire enhancing region in each slice (excluding the ventricles). The number of pixels in the ROIs was then counted and multiplied by the volume of a single pixel to receive the ROI volume. The slice thickness was 2 mm and in-plane pixel size was 0.3 X 0.3 mm. Radii of each slice was then extracted by calculating the biggest radius based on the Euclidean distance transform of the corresponding slice. The biggest radii computed over all slices were chosen as IRE radius and BBB disruption radius.

The radii rd and rb where then plotted as a function of the number of treatment pulses (N) to determine the dependence of the radii on the number of treatment pulses.

Numerical modeling

The mathematical models were based on a two-dimensional finite element model (assuming spherical symmetry of the produced IRE lesions and BBB disruption) (Figure 1) that was implemented in the COMSOL software package (Comsol Multiphysics, v.4.2a; Stockholm, Sweden) as previously described.7,34

Figure 1.

Simulation results. (A) Electric field distribution in the numerical model. The shape of the field assumes a nearly spherical shape. (B) Temperature distribution after 540 pulses. (C) Model geometry including the location of the electrodes (red arrows)

The rat head and chest were modeled as a 30 × 15 mm ellipse (Figure 1C) with an initial conductivity of 0.258 S/m to match the conductivity used by Sel et al.37 The electric field was described by the Laplace equation for electric potential distribution in a volume conductor: where σ is the electric conductivity of the tissue, E is the applied electric field and φ is the potential. The σ(E) dependence of brain tissue was described by an smoothed Heaviside function using 500 V/cm and 700 V/cm as reversible and irreversible thresholds.37,39 These values where used traditionally for 90 pulses and are recalculated using the Peleg-Fermi model for the different number of pulses in this paper.

Dirichlet boundary condition was applied to the surface of the electrode: and to the ground where φ0 is the applied potential on the intracranial electrode.

The boundaries where the analyzed domain was not in contact with an electrode were treated as electrically isolative and Neumann boundary condition was set to zero on the outer border of the model: where n denotes the normal to the boundary.

Thermal modeling

Control of the temperature during EP treatments is important in order to avoid damage to unwanted regions. The goal is to achieve complete coverage of the targeted tissue with sufficiently high electric field while ensuring that the temperature increase during the procedure does not generate thermal damage. The thermal effects of EP were determined from solution of the modified Pennes’ bioheat equation (equation [5]) in the 2D numerical model with the inclusion of the Joule heating source term. A duty-cycle approach was used, in which a time dependent solver for the duration of the treatment was applied and the thermal dissipation was multiplied by the pulse length. where k is the thermal conductivity of the tissue, T is the temperature, wb is the blood perfusion, cb is the heat capacity of the blood, Ta is the arterial temperature, q’’’ is the metabolic heat generation, ρ is the tissue density, cp is the heat capacity of the tissue and is the local voltage amplitude. Qmet accounts for Joule heating, where φ is the electrical potential and σ is the electrical conductivity of the tissue. The initial brain temperature was set to 37°C to match human temperature although anesthesized rat temperature is around 32°C. The parameter values utilized in the bioheat equation were taken from the literature40 and were used by others to follow/measure temperature increase and determine possible thermal damage due to EP treatments.41,42 All parameters used in the simulations are summarized in Table 1. The thermal properties of the silver plating and copper were taken from the Comsol Multiphysics database.

Table 1

Material properties used for numerical model

Brain	σ - basic conductivity	0.258[S/m]
	k - Thermal conductivity	0.0565[W/(m*K)]
	Cp - Heat capacity	3680 [J/(kg*K)]
	- density	1039 [kg/m^3]
	Q’’’- metabolic heat generation	10437 [W/m^3]
	T - temperature	37°C
Blood	Cp-heat capacity	3840 [J/(kg*K)]
	density	1060 [kg/m^3]
	Wb-Perfusion rate	7.15E-3 [1/s]
copper	σ - basic conductivity	5.998E7 [S/m]
	k - thermal conductivity	400 [W/(m*K)]
	Cp heat capacity	385 [J/(kg*K)]
	- Density	8700 [kg/m^3]
Silver	σ - basic conductivity	6.273E7 [W/m^3]
k - thermal conductivity	429 [W/(m*K)]
	Cp heat capacity	234 [J/(kg*K)]
	- Density	10500 [kg/m^3]

Statistical modeling

The original Peleg-Fermi model computes the ratio (S) of surviving bacteria after EP. Here we extended this model to describe the effects of EP on brain tissue as following:

First, the model was adapted to predict tissue damage (cell death) probability induced by EP. In the Peleg-Fermi model the probability for cells survival is given by: where E is the electrical field, N is the number of pulses, Ec is the critical electric field in which 50% of the cells are killed and A is a kinetic constant which defines the slope of the curve.

The electric field calculated using the numerical model was exported to Matlab (R2011a, Mathworks, USA) and was implemented in the Peleg-Fermi model.

We have previously shown that the hyper-intense regions on T2-weighted MRI obtained 14 days post treatment were significantly correlated with rarified regions in histology, confirming that these regions represent damaged tissues.7 Based on this, rd was set as S(E,N) = 0, assuming over 99.99% of the cells were irreversibly electroporated.

An optimization based on A Nelder-Mead simplex algorithm43 with added constrains was applied to Equation [7] for each group treated with N pulses, calculating a map of S(E), until r(S = 0) matched rd. The coefficients Ec and A for the different number of pulses were extracted and behavior equations were fitted to Ec(N) and A(N).

For each group treated with N pulses, Electric field distribution (E) was calculated using COMSOL Multiphisics and extracted to Matlab. The map E, along with the equation [7] allows to associate a map of S with any pair of the Fermi distribution (Ec,A). From the S map, the two isocontours of S = 0.9999 and S = 0.0001 are fitted to circles. An optimization based on A Nelder-Mead simplex algorithm on Ec and A as variables is used to find the (Ec,A) pair of parameters best corresponding to rd/rb.

The process of extracting r from S is nonlinear as it is based on fitting S = 1 iso-contour to a circle. Therefor the global dependency between Ec, A and rb/rd is noisy. This noisiness could potentially cause problems with computation of derivatives. Additionally, since S is monotonous there is no risk of the simplex finding a local minimum.

In order to assess whether the goodness of the Ecd(N) and Ad(N) (Ec and A for IRE) fits to the experimental data, r(S = 0) for different number of treatment pulses was calculated and compared to rd.

Although the Peleg-Fermi model was originally used to describe cell death, here we adapted it to describe BBB disruption as well and calculated the relevant coefficients. For this purpose the model was fitted to the radii calculated from contrast-enhanced T1-weighted MR Images. This time rb was set as BBB(E,N) = 1, meaning less then 0.001% of the BBB was disrupted in radii larger than rb.

After determining Ecb and Ab (Ec and A for BBB disruption) for each N and the behavior equations Ec(N) and A(N), the goodness of the fit to the experimental data was evaluated by recalculating r(BBB = 1) and fitting it to rb.

The electrical field threshold for cell kill, i.e. IRE extent and for BBB disruption, i.e. reversible EP extent, for different N values were then extracted from the results of the model and compared with thresholds reported in the literature.

Results

MR images of 46 rats that were previously treated with EP as described above were included in the current analysis. Treatment parameters were 600 V, 50 μ s pulses at 1 Hz with varying number of treatment pulses from 10 to 540 pulses. The extent of tissue damage and BBB disruption, i.e. rd- the irreversible damage radius and rb – the BBB disruption radius were calculated from the MR images acquired 30 minutes post treatment and 2 weeks post treatment as described in the Methods section.

Dependence on the number of treatment pulses

The average radius of each treatment group as calculated from the MR images is presented in Table 2.

Table 2.

Average radii of IRE and BBB disruption for each treatment group. Each group of 5–7 rats was treated with different number of pulses (10–540) at 600V, 50μ s pulses at 1Hz

# of pulses	10	45	90	180	270	450	540
IRE radius (mm)	0.62 ± 0.15	1.35 ± 0.18	0.89 ± 0.20	1.42 ± 0.15	1.37 ± 0.16	1.92 ± 0.07	1.80 ± 0.21
BBB disruption radius (mm)	1.25 ± 0.06	1.74 ± 0.04	1.84 ± 0.07	2.54 ± 0.15	2.19 ± 0.14	2.84 ± 0.04	2.69 ± 0.12

The dependence of r on N has been previously described by both logarithmic and power functions.44 Here, by fitting the mean rd of each treatment group to the number of electric pulses - N, we found the logarithmic function to provide a better fit to the data, resulting in the following dependence of rd on N: r

Similarly, by fitting the mean rb of each treatment group to N, the dependence of rb on N was found to be: r

rd and rb can be seen in Figure 3. The average ratio between rb(N)and rd(N) was found to be 1.67 ± 0.11 (s.e.m), confirming the coverage of significant volumes surrounding the IRE with BBB disruption. The small error suggests that the ratio between rd(N) and rb(N) is not affected by the number of applied pulses. The ratio between rd(N) and rb(N) plotted as a function of the number of treatment pulses supports this observation (Figure 3B). The coefficients of the empirical function for the BBB disruption are higher, because the BBB is disrupted by electric fields lower than those required for IRE ablation.

Figure 3. (A)

radii of irreversible damage and BBB disruption calculated from the MRIs, as a function of the number of treatment pulses, and the logarithmic equations fits (B) ratio between rb(N) and rd(N) as a function of number of the number of treatment pulses.

Irreversible damage model

The coefficients Ecd and Ad of equation [7] were calculated for each value of N as shown in Figure 4A–B. In order to find Ecd and Ad we used rd values obtained from equation [8] rather than using the average values obtained from the experiments, as this equation describes the dependence of rd on the number of treatment pulses based on the experimental data. Although Ec(N) is traditionally described with an exponential function we chose to describe it here using a power function as it fitted the data considerably better (r2 was considerably larger: 0.89 for the power function versus 0.5 for the exponential function), especially in the high N range. Still, when fitting the optimization results of Ec(N) of only the first 90 pulses to an exponential function, r2 increased to 0.83 (Figure 4C).

Figure 4.

Dependence of Ecd (A) and Ad (B) on the number of treatment pulses. (C) Exponential dependence of Ecd on the number of treatment pulses with N limited to 90 pulses. (D) Correlation between radii obtained from experimental data and radii obtained from the statistical model for IRE. Error bars represent 95% confidence level.

MRI example. (A) 3 slices of contrast-enhanced T1-weighted MR images of a rat treated with 45 electroporation pulses. The MR images was obtained 30 min post treatment. Each slice is 2mm thick. The enhancing region represents BBB disruption. (B) ROI (green) plotted in the MR image to mark the enhancing region.

Linear regression analysis confirmed that r(S = 0), calculated from the Peleg-Fermi equation with Ecd(N) and Ad(N)described well the rd obtained from the experimental data: F(1,5) = 45, p < 0.008, r2= 0.79. The resulting regression equation was: rd = 0.19 + 0.87 x, (x = r(S = 0)).

BBB disruption model

The same optimization method that was used to calculate Ecd and Ad was applied to the BBB disruption data with rb set to BBB(E,N) = 1, meaning that for radii larger than rb the BBB was not breeched(less than 0.001% ). Ecb and Ab were calculated for each value of N as can be seen in Figure 5A–B.

Figure 5.

Dependence of Ecb (A) and Ab (B) on the number of treatment pulses for BBB disruption. Error bars represent 95% confidence level.

As for the IRE models, the goodness of the fit to the experimental data was also evaluated. r(BBB = 1) and r(BBB = 0) were calculated from BBB(E,N) for each value of N using Ab(N) and Ecb(N).

Next we evaluated the correlation between rb(N) obtained from the experimental data, and r(BBB = 0) linear regression analysis confirmed that r(BBB = 1), calculated from the extended Peleg Fermi model with Ecb(N) and Ab(N) described well the behavior of rb obtained from experimental data: F(1,5) = 45, p < 0.001, r2= 0.91. The regression equation was: rb = 0.19 + 0.87 x (x = r(S = 0)).

Electric field thresholds

The electrical field thresholds for E(S = 0) and E(S = 1) were calculated from the model for cell death. E(BBB = 0) and E(BBB = 1) were calculated for BBB disruption. In the cell death model, E(S = 0) represents the threshold needed for over 99.99% cells death whereas electric field lower then E(S = 1) will cause cell death lower than 0.001%. In the BBB disruption model E(BBB = 0) represents the threshold needed for over 99.99% of the BBB to be breeched while electric field lower then E(BBB = 1) will not disrupt the BBB(BBB disruption lower than 0.001%). Thresholds are presented in Figure 6.

Figure 6.

Electrical field thresholds. (A) IRE thresholds. Dashed line represents published IRE thresholds for white matter for 80 50 μ s pulses at 4 Hz. (B) BBB disruption thresholds. Dashed line represent previously published threshold for 90 50 μ s pulses at 4 Hz.5 (C) Thresholds for E(S = 0) for the IRE and E(S = 1) for BBB disruption. (D) Ratio between E(S = 1) and E(S = 0) for IRE and E(BBB = 0) and E(BBB = 1 ) for BBB disruption. Error bars are smaller than markers.

The ratio between S(E,N) = 0 and S(E,N) = 1 thresholds, representing the transition zone between over 99.99% cell death and no cell death (S(E,N) = 0 / S(E,N) = 1) thresholds was calculated. The ratio is relatively high (between 0.88 and 0.91) even for small numbers of pulses and depends only weakly on the number of treatment pulses. This means that the transition between 99.99% cell death threshold and no cell death threshold is narrow and gets even narrower for large numbers of treatment pulses. This is not the case with BBB disruption, where the ratio between the thresholds increases with the number of treatment pulses eventually converging to one (Figure 6). The average ratio between BBB disruption ratio and damage ratio is 1.67 ± 0.11 (s.e.m).

Thermal model

The initial temperature of the rats’ brain in the simulation was set to 37°C to show that the treatment should not induce thermal damage in clinical use. The maximum temperature reached in the tissue after treatment at 600 V was 38.9°C (Figure 1B). This temperature was reached using 540 pulses. Although temperature was not measured during EP treatment, histological analysis of brains extracted 60 min post treatment revealed no signs indicative of thermal damage such as coagulation, or extensive hemorrhages.45 Connective tissue and blood vessels were preserved in the treated area suggesting damage induced only by IRE.46

Discussion

When treating tumors by EP, it is important to deliver the electric pulses so that the entire tumor volume will be treated to avoid recurrence. It is also vital to treat the infiltrating zone surrounding the tumor mass with high efficacy while preserving the healthy tissue. This is especially important in the case of brain tumors, where the infiltrative zone is relatively large47 and the preservation of healthy brain tissue is of critical importance.

The electrode configuration in this experiment, which consists of a single intracranial insulated electrode with an exposed tip combined with an external surface ground electrode, provides an electric field distribution that is strongest at the intracranial needle tip tissue interface and decreases with the square of the radius. This setup provides a well-controlled region of permanent damage induced by IRE, further surrounded by significant BBB disruption zone. This combined response offers the potential of this configuration for the treatment of brain tumors combining IRE and chemotherapy. This setup induces rapid tissue damage in the tumor mass surrounded by significant BBB disruption, thus potentially enabling efficient drug delivery of systemically administered drugs to the infiltration zone surrounding the tumor. Since GBM cells are highly resistant to therapeutic apoptotic stimuli, however, they exhibit a paradoxical propensity for extensive cellular necrosis19,20, IRE may be efficient for treating the tumor mass. Disrupting the BBB in the local vicinity of the tumor can also improve drug intake since peripheral administration of therapeutic agents is inefficient due to poor penetration of most drugs across the BBB. As clinical trials using IRE or ECT for the treatment of deep seated tumor are becoming more common8–16, a model that can predict treatment outcome and enable individual treatment planning is increasingly being recognized as a need.

The statistical model of EP-induced cell death used in this manuscript was originally suggested by Golberg et al.32 who validated the model using experimental data in vitro. Here, we present for the first time an experimentally validated statistical model for tissue IRE where the Peleg-Fermi model was extended to a wide range of number of treatment pulses r and to BBB disruption.

The results of this study demonstrate the feasibility of applying the Peleg-Fermi model for describing irreversible EP in the brain and for treatment planning. Furthermore, since our model is based on experiments with up to 540 pulses we were able to extend the model beyond the up to 90 traditional pulses used for IRE. This is important since protocols outside the traditional 100 pulses are being evaluated48,49 and a tool to evaluate protocols with higher number of pulses is needed.

The results of the model indicate, as expected, that with increasing number of treatment pulses it is possible to treat larger volumes of tissue and that the IRE threshold decreases with the number of pulses. This is however true up to a limited extent since both Ec(N) and the thresholds eventually plateau. This suggests that although increasing the number of pulses while lowering the treatment voltage may represent a safe way to avoid thermal damage while still achieving large enough treatment volumes, there is an upper limit for this effect. In addition, when using higher voltages, raising the number of pulses will eventually lead to increased damage induced by Joule heating but will not increase the damage induced by IRE. This plateau phenomenon does not only result from the logarithmic behavior of r(N),as can be seen in Equations [8], [9] (Figure 6) but also in Ecd(N) (Figure 4B) and Ecb(N) (Figure 5B).

Although the behavior of the equations describing Ec(N) was previously described as exponential32,36, which supports the claim that larger number of pulses increases lethality, we found that Ec(N) is better described by a power function. As power functions plateau faster than exponential functions it further supports the limited effect of increasing the number of treatment pulses. One explanation for the difference might be that previously the model was limited to 10-90 pulses32,36, and therefore the plateau effect was not yet reached. In a paper published about evaluation of the Fermi equation as a model of dose-response curves on dose response the author described Ec as a Weibull function suggesting exponential function is just a simplification for limited range of pulses.50

Once we found that the Peleg-Fermi model can be used to describe IRE, we continued to further extend the model to describe BBB disruption induced by EP. For the model, we correlated radii calculated from contrast-enhanced T1-weighted MRI with BBB(E,N) = 1 since even a relatively small BBB disruption can be visible. We found that the extended Peleg-Fermi model describes well not only the behavior of IRE radii but also that of BBB disruption induced by EP with high statistical significance. This also indicates that there are possibly similar underlying mechanisms at play, which cause the effects.

The combination of the two models can be used for efficient treatment planning for brain tumors where IRE is used for ablating the tumor mass while BBB is disturbed in the rims and infiltrating zone thus allowing efficient access of therapeutic agents. The rims in this setup are on average 1.67 ÷ 0.11 mm wider than the damage, with no correlation to number of pulses. This suggests that the volume of BBB disruption is over 4 times larger than the volume of IRE.

Although both the IRE and the BBB disruption models were constructed separately, when planning a treatment protocol for brain, both should be used since BBB disruption with no irreversible damage only occurs in relatively low electric fields and when higher voltages or higher number of pulses are used, irreversible damage is difficult to avoid.

The electric field thresholds for IRE and reversible EP are mostly limited to the traditional treatment protocol, i.e. 90 pulses, but when using protocols that include a different number of pulses, different thresholds should be used.51 In this study we calculated the thresholds needed for the different number of pulses and fitted them to a power function. The thresholds we found for 90 pulses fit well within the thresholds previously reported in the literature for IRE in white matter52 and BBB disruption5 as can be seen in Figure 6.

The ratio between the thresholds of BBB(E,N) = 0 and BBB(E,N) = 1 were found to increase with the number of pulses suggesting that the window between 99.99% of BBB disruption and no BBB disruption narrows with the number of pulses. This suggests that while increasing the number of pulses will eventually not lead to bigger radius of BBB disruption, larger percentage of the BBB will be disrupted thus improving drug penetration to the tissue. This is not the case for IRE where the ratio between the thresholds for S(E,N) = 0 and S(E,N) = 1 seems to be nearly independent on the number of treatment pulses. This is somewhat surprising but could be explained by the fact that the ratio is relatively high to begin with (between 0.88–0.91) and that our dataset starts with 10 pulses, however the range of pulse numbers in this study covers the most commonly used IRE protocols in clinical practice. This is also consistent with previous publications saying there is a sharp delineation of IRE treated and healthy tissue53 and demonstrates that the sharp delineation is maintained even for high number of pulses.

The ratio between rb and rd was found to be nearly independent on the number of pulses. This is further supported by the relatively constant ratio between the thresholds that where calculated from the model. Thus, during treatment planning it might be sufficient to calculate one radius. It also indicates that BBB disruption may be used as a safety limit for irreversible EP. Though when using other electrode configurations, caution is needed. If thermal damage occurs, typically at high voltages or high number of pulses, it may influence the ratio between cell death and BBB disruption.

The thermal model showed only a mild increase in brain temperature. The maximal temperature at the end of 540 pulses reached 38.9°C. Since 42°C is often considered the thermal damage threshold if sustained for long durations42, it is safe to assume that the tissue damage found in our experiments was induced solely by EP and not by thermal effects. This was also confirmed by histology7 showing no signs of thermal damage, although a temperature assessment in real time is advisable.

Despite our understanding that this model may be used by physicians and researchers for the selection of treatment protocols, a model that also incorporates dependence on additional treatment parameters such as frequencies and pulse durations should be developed.28 Such all-inclusive model would enable physicians to choose the safest and most efficient protocol on a per-patient bases. Another point to bear in mind is that although the electrode configuration suggested in this paper produces very low Joule heating, using other electrode configurations with high number of pulses might induce thermal damage in addition to IRE.54 Although this study indicates that the combination of IRE and BBB disruption may be applied for the treatment of brain tumors, experimental validation using animals bearing intracranial tumors is yet to be done.

In conclusion, the results of our study indicate that it is possible to apply high voltage electric pulses in a manner that induces localized focused irreversible damage in the brain surrounded by a larger volume of BBB disruption while using a single minimally invasive intracranial electrode. We used existing statistical models of cell kill by electric pulses that were based on theoretical cases and validated them using in vivo experimental data and extended the knowledge of EP thresholds beyond the traditional 90 pulses protocol used in IRE. We further extended the model to describe BBB disruption induced by EP. These models can assist physicians and researchers in selecting optimal treatment protocols allowing them to achieve the desired outcome in treating brain tumors. Although validation of the model in tumors is yet to be done, the results confirm that treatments outside the most commonly used protocols can achieve expected outcome.