Dynamic Electroporation Model Evaluation on Rabbit Tissues.

Biological electroporation is a process of opening pores in the cell membrane when exposed to intense electric fields. This work provides results for validation of a dynamic model of electroporation on biological tissues. Computational simulations were carried out and results for the electrical current through the tissue and increase of the tissue temperature were compared to experimental results. Two calculation methods were used: Equivalent Circuit Method and Finite Element Method. With Equivalent Circuit Method the dielectric dispersion present in biological tissues was included. Liver, kidney and heart of rabbit were used in the experiments. Voltage pulse protocols and voltage ramps were applied using stainless steel needles electrodes. There is good agreement between the simulated and experimental results with mean errors below 15%, with the simulated results within the experimental standard deviation. Only for the protocol with fundamental frequency of 50 kHz, the simulation performed by the Finite Element Method using a commercial software did not correctly represent the current, with errors reaching 50%. The justification for the error found is due to the dielectric dispersion that was not included in this simulator.

Introduction

The application of intense electric fields on biological tissues or in cells in suspension causes the opening of hydrophilic pores in the cell membrane, allowing the exchange of substances between the intracellular and extracellular media,12 this phenomenon is called electroporation. The applied electric field causes the accumulation of ions close to the membrane, increasing the transmembrane potential.

Electroporation has applications in several areas such as medicine, biology and biochemistry.15 Some of their important uses are: in the extraction and transfer of DNA and proteins,25 in facilitating the entry of chemotherapeutic agents into the cell and thus, increasing the efficiency of drugs9; in the vaccine development30 and the in the destruction of microorganisms with irreversible electroporation.1

Biological electroporation has been known for several decades, where many experimental aspects for the characterization of the sample have been explored such as: the use of fast-freezing microscopy in human erythrocytes,2 fluorescent markers sensitive to voltage,14 transfer of ionic markers and macromolecules to the intracellular medium,19,22 protocols for clinical applications,21 cell destruction with irreversible electroporation.18 The application of extreme electric fields far above the electroporation threshold causes cell destruction.20 This technique has been researched for the treatment of cancer in situations where the use of chemotherapy is not possible.10,34 Despite obtaining several experimental results so far, the mathematical and physical modeling of the phenomena involved still proves to be somewhat complex. This is due to the opening time of the pores being in the order of microseconds and the size of the nanometer pores, in addition to the fact that biological tissues have intense dielectric dispersion.27

The mathematical models for numerical calculation of the conductivity variation of biological tissues during electroporation are static or dynamic. The static models use the relation between electric field and conductivity value only at the end of the pulse, ignoring variations during the pulse.3,4 The dynamic models, on the other hand, provide a means for obtaining the conductivity variation over all the time domain. Langus et al.16 proposed a dynamic model for the simulation with the Finite Element Method (FEM) and obtained simulated results showing an error of only 5% in comparison with the experimental current curves. However, the proposed model is composed of fourteen parameters that must be adjusted. In addition, the authors point out that they have no theoretical explanation about the adequacy of such a set of equations for the dynamic modeling of the biological electroporation process. Ramos and Weinert27 presented a dynamic electroporation model for biological tissues based on the asymptotic model proposed by Debruim and Krassowska.7 The model contains six parameters and obtained a good representation for the electric current in the sample for one pulse of the applied field. Weinert et al.33 proposed dynamic changes in the parameters of the model by Ramos and Weinert,27 obtaining good agreement between numerical and experimental results for the electric current in the sequence of ten pulses. However, the model has become complex, increasing the number of parameters to be adjusted. Voyer et al.31 also presented a dynamic model for biological tissue based on the work of Leguebe.17 This model showed simulated response similar to the experimental values of the electric current, however it works only for one pulse and only in the pore opening process, that is, the recovery of the membrane was not modeled.

The present work uses the electroporation model proposed by Ramos and Weinert27 for the simulation and comparison with experimental data of electroporation in three rabbit tissues (heart, kidney and liver) with different protocols. The numerical simulations are performed with two calculation methods: Equivalent Circuit Method—ECM27 and Finite Element Method—FEM. In the ECM, the dielectric dispersion of the tissue was included in the simulation. In addition to the electrical data, the increase of the tissue temperature was also calculated and compared with experimental measurements.

Materials and Methods

Biological Tissues

Samples of liver, kidney and heart of rabbits of the Leporidae strain Oryctolagus cuniculus species, were used. All procedures performed with the animals were approved by the Research Ethics Committee on the Use of Animals - CEUA/UNIVILLE, (process 02/2019). The liver and kidney samples were obtained with 20 mm × 20 mm × 10 mm edges. The heart does not have a homogeneous structure and its interior is hollow due to the atria and ventricles. Thus, such an organ was fragmented in half, having dimensions of 10 mm x 10 mm and 5 mm thick.

Electroporation Protocols

The electroporation protocol used in chemotherapy treatments consists of the application of 8 pulses, with a duration of 100 μs and an interval of 1 s between the pulses.21 However, the objective of this work is to evaluate whether the electroporation model reliably represents the tissue conductivity variations under some different protocols as follow:

Protocol (A) 10 voltage pulses, time at high and low level 99.5 μs, rise and fall time 0.5 μs. Voltage values: 200, 500, and 800 V.

Protocol (B) 10 voltage pulses, time at high and low level 9.5 μs, rise and fall time 0.5 μs. Voltage values: 200, 500, and 800 V.

Protocol (C) Voltage ramp, with linear increase from 0 to 400 V within the following time intervals: 100, 200, 300, 400, 500, 600, 700, 800, and 2400 μs.

For the application of protocols A and B a system of parallel stainless steel needles separated by a distance of 5.4 mm was used. The needles have a diameter of 0.6 mm and a height of 15 mm. Thus, the maximum fields obtained are: 370, 926, and 1481 V/cm for voltages of 200, 500 and 800 V, respectively. In protocol C parallel needle electrodes with a spacing of 2.8 mm, needle diameter of 0.7 mm and a height of 10 mm were also used. Protocol A and B were applied to five tissue samples, while for protocol C three samples were used.

Electroporator and Thermal Measurements

Two electroporators were used, one used for protocols A and B with the capacity to supply 800 V and 10 A. Another for protocol C which consists of a linear amplifier with a capacity of 400 V and 10 A. Waveforms voltage and current were measured and digitized and sent to the computer for analysis. To read the temperature on the sample surface, a FLIR® thermal camera, model FLIR C2, was used. The camera was positioned 5 cm above the sample. Data acquisition was performed using the FLIR tools plus software.

Electrical Properties of Tissues

Dielectric dispersion in biological tissues occurs mainly due to the accumulation of hydrated ions on the cell membranes, but also due to dipole relaxation of water and polar macromolecules in tissue fluids.26

The Cole–Cole empirical model is widely used to describe the dielectric dispersion of biological tissues due to the use of a distribution of relaxation times to represent each dispersion band11:where ε∞ is the electrical permittivity at high frequencies, σs is the static conductivity, Δε is the magnitude of the electrical permittivity dispersion, εo is the vacuum permittivity, τ is the central relaxation time of the distribution, n is the index referring to the dispersion and α is a parameter to be adjusted empirically with a value less than the unit.

The conductivity and dielectric constant spectra of the intact tissue were obtained with the aid of the Agilent® impedance analyzer model 4294A. The samples of rabbit tissues were placed in a circular parallel plate system, with a thickness of 3 mm and a diameter of 15 mm. With the aid of a genetic algorithm from the Matlab® software, the parameters of the Cole–Cole model for the beta band dispersion were obtained. The genetic algorithm was configured with the following parameters: population size = 100, crossover = 0.85 and generations = 50. Such parameters of the Cole–Cole model are used in the computational simulation of tissue electroporation with ECM.

Dynamic Model of Electroporation

The dynamic model of electroporation was proposed by Ramos and Weinert27 and is based on the asymptotic model proposed by Debruim and Krassowska.7 According to the dynamical model, Eqs. (2) and (3) below show how to calculate the rate of variation of the tissue conductivity as a function of the electric field intensity:where σp is the conductivity of electroporation, σo is the initial conductivity, E1 and E2 are electroporation thresholds and E(t) is the applied electric field. τmin and Δτ are relaxation times. The σ parameter has important functions for the model. It limits the value of σp and permits to model the electric field division between cell membrane and surround environment resulting in the reduction of the transmembrane potential when the electroporation becomes intense.

Thermal Model

In the analysis of the biological electroporation, the evaluation of the electric field distribution has central importance to define the extension of the electroporated area. Besides, side effects must be considered in the planning of electroporation protocols, the most important being thermal heating.6 In our experiments, the calculation and measurement of the temperature distribution in the electroporated tissue served as an additional check for the quality of the electroporation model.

The mathematical model used for the transfer of heat in biological tissues is known as the heat transfer equation:where k is the thermal conductivity of the tissue, T is the temperature in Kelvin, wb is the perfusion rate of the blood, cb is the thermal capacity in the blood, Ta is the temperature of the artery, λ is the generation rate of metabolic heat, ρ is the density of the tissue and cp is the specific heat of the tissue.23

The thermal model implemented in the simulators is based on Eq. (4) with some adaptations. The terms of blood perfusion and metabolic generation were neglected since the experiments were carried out in ex vivo samples. Besides, a thermal dissipation term due to the electric stimulus was included in the tissue volume represented by ρdiss:

Simulation

The simulations of biological electroporation were performed in two software: The ECM program developed in our research group and the commercial program COMSOL Multiphysics®.

The ECM was proposed by Ramos et al.,28 as a calculation method of the electric field and electric current distribution in a medium with low conductivity and high dielectric constant. It is based on lumped circuit elements, conductances and capacitances, that model the electrical conduction and polarization processes in the materials involved.26–28 To take into account the effects of dielectric dispersion in the ECM, the dispersion parameters of the β band were included.

In COMSOL Multiphysics® the following set of equations from the electromagnetic theory was solved in the time domain:

The continuity equation is shown in (6). Equation (7) are the current density components with being the conducting current density, the displacement current density. Conductivity σ is the sum of two components σ = σs + σp, where σs is the static conductivity of the tissue and σp the conductivity of electroporation. Equation (8) is the electrostatic relationship between the electric field and the electric potential.

In both simulators, the following boundary conditions are defined: electric potential defined at the border with a metallic electrode and null perpendicular electric field at the other borders of the analysis domain.

For the ECM simulation, the following rectangular discretization mesh was used: 100 divisions on the x and y axis, 10 divisions on the z axis, making a total of 100,000 elements. Thus, the edges of an element in the x and y direction are 9.9 × 10−5 m long and in the z direction 9 × 10−4 m. The COMSOL Multiphysics® software generated meshes with the following number of elements: Protocols (A) and (B)—48,969 elements; Protocol (C) 39,773 elements.

Results

Conductivity of rabbit tissues

Figure 1 shows the conductivity spectra for the heart, kidney and liver of the rabbit in the frequency range from 1 kHz to 10 MHz. These curves are averages obtained from measurements of five samples. The temperature of the samples was 22 °C. Table 1 presents the parameters of the Cole–Cole model obtained with the aid of the genetic algorithm. In this frequency range, the main dispersion process is the interfacial polarization on cell membranes that gives rise to the beta band.

Figure 1

Measured conductivity of the heart (o), kidney (.) and liver (*) of rabbit. M ± SD stands for average plus or minus standard deviation.

Table 1

Parameters of the Cole–Cole model of heart, kidney and liver of rabbit.

Parameters	Heart	Kidney	Liver
σs (S/m)	0.127 (0.0327)	0.0589 (0.003)	0.0475 (0.0166)
Δσ (S/m)	0.076 (0.0079)	0.0479 (0.0076)	0.0957 (0.023)
τ3 (μs)	7.93 (0.61)	11.1 (1.88)	7.17 (0.71)
β	0.41 (0.032)	0.39 (0.027)	0.37 (0.040)

Values are given as average (Standard Deviation)

Waveforms of Electric Current and Conductivity in Rabbit Liver for Voltage Pulses

In the following, the waveforms of electric current during electroporation are presented only for rabbit liver. The graphics for the other tissues are similar, changing only the amplitude of the waveform. Figure 2 show the experimental and numerical waveforms of electric current for rabbit liver obtained with protocol A. Figure 3 show the results with protocol B. In both protocols the nominal electric fields, calculated as the applied Voltage amplitude divided by the distance between needles, are: 370 V/cm for voltage of 200 V, 926 V/cm for 500 V and 1481 V/cm for 800 V.

Figure 2

Waveform of experimental and numerical electric currents for rabbit liver samples stimulated with 5 kHz (Protocol A). The ECM no EP curve is the simulation result without electroporation. Exp is the average of 5 samples. SD is the experimental standard deviation. Nominal electric fields: (a) 370 V/cm; (b) 926 V/cm; (c) 1481 V/cm.

Figure 3

Waveform of experimental and numerical electric currents for rabbit liver samples stimulated with 50 kHz (Protocol B). The ECM no EP curve is the simulation result without electroporation. Exp is the average of 5 samples. SD is the experimental standard deviation. Nominal electric fields: (a) 370 V/cm; (b) 926 V/cm; (c) 1481 V/cm.

In Fig. 2 there is an increase in current during the first pulse, which is attributed to the creation and expansion of hydrophilic pores in the cell membranes. In the subsequent pulses there are additional increases in current. Such memory effect occurs because the pores created during the previous pulse are unable to close in the time interval between the pulses. It is proposed that this effect is due to lipid peroxidation of the cell membranes.29 However, the increase in current reaches a limit. Saturation occurs when the number of pores created only compensates for the closed pores in the interval between pulses. This tends to happen as the initial conductivity in each pulse is enough high.

With the increase in the applied voltage, it is observed that the electric current in each pulse increases and the electroporation becomes more intense. For the smallest nominal electric field of 370 V/m, the current variation in the pulses is much smaller than in the other cases, suggesting that this field value is close to a lower limit (threshold field) for the occurrence of the phenomenon. However, when comparing this curve with that obtained numerically for the condition without electroporation (ECM no EP), it can be concluded that electroporation is in fact occurring.

In Protocol B the fundamental frequency of the pulse train was increased to 50 kHz. The experimental current curves in Fig. 3 have practically no evidence of electroporation due to the intense displacement current that overlaps the current resulting from the opening of pores in the tissue. However, when comparing with the simulated curve in the ECM without electroporation (ECM no ep), it is evident that the phenomenon is occurring. For comparison, the frequency response of electroporation for cells in suspension was assessed by Gárcia-Sánchez et al.,8 where different sinusoidal protocols were used to assess the percentage of permeabilized cells, using fluorescence spectroscopy. Up to 50 kHz the percentage of permeabilized cells is close to 90%, while at the frequency of 100 kHz the percentage of electroporated cells was less than 20%.8 This effect may be a consequence of the transmembrane potential reduction with the increase of frequency of the applied field.

The values of the electroporation model parameters used in the simulations in both protocols for liver, kidney and heart of rabbit are shown in Table 2. The parameters σo = 1 × 10−8 S/m, τmin = 1 × 10−6 s, Δτ = 1 × 10−6 s and E1 = 2300 V/m are the same in all simulations.

Table 2

Parameters of the electroporation model for liver, kidney and heart of rabbit and errors in the calculated electric current for the first pulse.

Rabbit tissue	σt (S/m) A	σt (S/m) B	E2 (V/m)	Error (%) A	Error (%) B
Liver 370 V/cm	0.135	0.040	7000	7.02 (10.53)	15.3 (39.9)
Liver 926 V/cm	0.350	0.150	20,000	0.08 (4.43)	14.4 (18.2)
Liver 1481 V/cm	0.465	0.200	30,000	1.48 (3.93)	10.2 (9.72)
Heart 926 V/cm	0.600	0.700	17,000	0.5 (2.5)	9.3 (11)
Heart 1481 V/cm	0.855	0.855	27,000	0.3 (2.7)	8.2 (8.6)
Kidney 370 V/cm	0.290	0.130	5000	0.8 (4.8)	4.8 (51)
Kidney 926 V/cm	0.345	0.270	17,000	0.2 (3.2)	6.0 (24)
Kidney 1481 V/cm	0.550	0.410	28000	0.3 (3.5)	2.8 (8.6)

The indices A and B refer to the 5 kHz (protocol A) and 50 kHz (protocol B), respectively. Averaged errors are informed for ECM (FEM)

The starting point for the parameter values was obtained from the work of Weinert et al.33 Small adjustments in the parameters were necessary to decrease the error between the experimental and simulated data. The parameters of the model depend on the type of tissue and the applied field strength. The dependence on the type of tissue can be justified by the different electrical properties of the same. The applied field strength influences the amount and size of pores created in the membrane, hence there is also variation in the model parameters, mainly on σ.

The errors shown in the last two columns of Table 2 were calculated taking the average difference between numerical and experimental values of the electrical current in the first pulse and using the average value of the experimental current as reference. It is observed that the simulations performed with the ECM show minor errors when compared to the FEM, mainly in the B protocol. This is because the dielectric dispersion was not included in the FEM simulations.

Figure 4 shows the mean conductivity of electroporation calculated with ECM in the volume with dimension of 9.9 × 4.8 × 9.0 mm in x, y and z axis, respectively, between electrodes. Indices a and b refer to the results for Protocols A and B, respectively. These curves corroborate the characteristics that have already been known of the conductivity during electroporation process: the increase during the pulse, the decrease between pulses, the effect of memory and the saturation of the process.

Figure 4

Mean conductivity of rabbit liver calculated with ECM as a function of time in the sequence of 10 pulses for 200, 500, and 800 V. Protocols A and B presented in figures a and b respectively.

Temperature profiles on the tissue surface

The temperature rise curves on the surface of the samples are presented along two axes: Axis (1) that cuts both electrodes; Axis (2) perpendicular to the previous one that is located at the midpoint between the two electrodes. Figure 5 shows the temperature distribution along these axes for Protocol (A) stimulated with the voltage of 800 V. The parameters of the thermal model are shown in Table 3, where ht is the convective coefficient at the surface of the sample.

Figure 5

Temperature profiles on the tissue surface for rabbit liver samples stimulated with 800 V and 5 kHz (protocol A). (a) Along Axis 1; (b) along Axis 2.

Table 3

Parameters of the thermal model.

Parameters	Tissues	Electrodes (stainless steel)
cp (J/kg K)6	3540	502
ρ (kg/m3)6	1079	7850
k (W/mK)6	0.52	15.9
ht (W/m2 K)6	500	–

Around the electrodes there is an intense heating, a fact expected due to the greater intensity of the electric field and the greater power dissipation in this region. The ECM and FEM simulation methods solve the heat transfer equation by different approaches, but the results obtained are practically identical. The purpose of calculating and measuring the temperature in the tissue is to have a countercheck to verify that the proposed electroporation model and the adjusted parameters result in good estimates of the tissue conductivity. Note that the numerical and experimental results are very similar.

Voltage Ramps

Figure 6 shows the experimental and simulated current curves for the voltage ramps. Results are presented for the rise times of 100, 500, and 800 μs. The waveforms of the experimental current have two types of behavior: at the beginning of the ramp there is a linear behavior, in which electroporation is not occurring; when the electric field reaches some threshold, the time derivative of the current increases, changing the slope of the waveform.

Figure 6

Waveforms of electric current for rabbit liver samples stimulated with Protocol C): (a) rise time of 100 μs; (b) 500 μs; (c) 800 μs. The ECM no EP curve is the simulation without electroporation. Exp is the average of 5 samples. SD is the experimental standard deviation.

The model was able to adequately represent all the ramps since the average errors between the curves of experimental and simulated currents fluctuate around 6%. Pliquett et al.24 used trapezoidal waveform voltage pulses with a rate of 150 kV/s and obtained similar results. Such an effect is also seen in patch clamp experiments at the cellular level.5,13 Voltage ramps can be used to obtain the electroporation thresholds as highlighted by the work of Weinert and Ramos.32

Only σ and E2 model parameters were adjusted for minimizing the averaged errors. The values obtained are: σ = 0.3, 0.18 and 0.23 S/m and E2 = 35, 40.8 and 40.8 kV/m for rise times of 100, 500, and 800 μs, respectively.

Discussion

The dynamical model of tissue electroporation proposed by Ramos and Weinert27 describes the variation of the tissue electrical conductivity as a relaxation process in which the parameters such as long term conductivity and relaxation time are functions of the electric field applied in the tissue. This model can be used only in numerical simulators in which the electric fields distribution is calculated in time domain.

The waveforms in Figs. 2 and 3 show that the numerical estimates of electrical current in the tissue for a train of voltage pulses applied between metallic needles is in good agreement with the experimental results for ECM method. In Table 2 the maximum error for the first pulse In Protocol A is 7% after parameter adjustments with the ECM. With FEM method the error reach 11%. In Protocol B the maximum error with ECM is 15% while with FEM is 50%. The difference between the numerical methods is the tissue dielectric dispersion, which was included only in ECM. In Protocol C the errors with both methods were about 8%. In Fig. 1 it can see that dielectric dispersion in liver, kidney and heart of rabbit is relevant in the frequency range from 10 kHz to 10 MHz. Thus, the dielectric dispersion influences electroporation dynamic and should not be ignored when high energy harmonics of the applied field localize inside the dispersion frequency band.

The parameters σ and E2 were adjusted for different levels of voltage and tissues and between different protocols. The other parameters did not change according to the protocol or the amplitude of the applied electric field. The dependence of σ on the intensity of the applied electric field is justified due to the level of pore creation in the cell membrane. This parameter is closely related to membrane and electrolyte conductivity. σt establishes a limiting value for σp and because the term Ememb = E σ/(σ + σp), it also limits the electric field in the membrane. The term Ememb/E2 defines the relaxation time. It is obtained that E2 increases almost linearly with the increase of the electric field. It can be concluded that the characteristic time of opening and closing the pores is little dependent on the applied electric field.

The experimental and computational results of the increase in temperature on the surface of the samples are similar and the differences can be justified in part by errors in the estimates of the parameters of the thermal model as well as by the low frequency of acquisition of the thermal camera (2.5 frames per second). Despite this, the similarities in the experimental and theoretical temperature distribution curves suggest that the electroporation model allows the realistic calculation of the increase in the conductivity.

The results presented in this study suggest that the empirical model of electroporation by Ramos and Weinert27 presents the appropriate characteristics for describing the variation in the electrical conductivity of biological tissues when stimulated with high intensity electric fields. Thus, it can be useful as an analysis tool in the computational simulation of this phenomenon with a view to planning experiments and therapies based on biological electroporation.