Dynamic model of tissue electroporation on the basis of biological dispersion and Joule heating

Abstract

Electroporation is a complex, iterative, and nonlinear phenomenon that is often studied by numerical simulations. In recent years, tissue electroporation simulations have been performed using static models. However, the results of a static model simulation are restricted to a fixed protocol signature of the pulsed electric field. This paper describes a novel dynamic model of tissue electroporation that also includes tissue dispersion and temperature to allow time-domain simulations. We implemented the biological dispersion of potato tubers and thermal analysis in a commercial finite element method software. A cell electroporation model was adapted to account for the increase in tissue conductivity. The model yielded twelve parameters, divided into three dynamic states of electroporation. Thermal analysis describes the dependence of tissue conductivity on temperature. The model parameters were evaluated using experiments with vegetal tissue (Solanum tuberosum) under electrochemotherapy protocols. The proposed model can accurately predict the conductivity of tissue under electroporation from 10 to 100. A negligible thermal effect was observed at 100, with a 0.89 increase. We believe that the proposed model is suitable for describing the electroporation current on a tissue scale and also for providing a hint on the effects on the cell membrane.

I Introduction

Electroporation occurs when the transmembrane potential (TMP) of a biological cell exceeds a supraphysiological threshold by stimulation of short but intense pulsed electric fields (PEF). Excessive TMP causes local disturbances in the membrane structure. Current electroporation theory suggests that prepores are formed [1]. The prepores then expand and stabilise as hydrophilic pores. Pore formation increases cell membrane permeability, allowing non-permeable substances to cross the cell barrier [2]. Extracellular content can access the cytosol, or intracellular content can leak and even trigger apoptosis by losing homeostasis [3, 4]. There are medical and industrial applications that use electroporation to improve or replace traditional processes [5, 6].

The occurrence of electroporation depends on the distribution of the electric field, which in turn depends on the nonlinear electrical properties of the tissue. Opening the pores in the membrane changes the structure of the material and its electrical properties [7]. The system is interdependent and leads to a complex model. For this reason, electroporation is often studied in detail by computer simulations.

Electrochemotherapy is a well-known application of electroporation to catalyse the membrane transport of chemotherapeutic drugs [8]. The technique relies on standard protocols to ensure that the entire tumour is exposed above a minimum electric field threshold for electroporation to occur. The European Standard Operating Procedures for Electrochemotherapy (ESOPE) recommend a burst of eight rectangular pulses (monopolar, bipolar, or alternating) 100 long with a repetition rate between 1 and 5 [9, 10]. Electrochemotherapy studies usually do not focus on the dynamics of pore formation, but rather on the final electric field distribution and outcomes. Static models are often used to reduce computational costs [11, 12]. However, a static electroporation model is developed specifically for a PEF signature and cannot be directly used to study different signatures. If the PEF signature is changed, the static model needs to be adjusted. The reason is that a PEF has a specific energy spectrum density that affects the dynamics of electroporation, Joule heating, and the dispersive aspect of biological media [13]. To overcome this limitation, dynamic models should be used.

A dielectric dispersive medium is characterised by the dependence of its electrical properties on frequency. Biological tissue exhibits a strong dispersion from DC to hundreds of $\mathrm{GHz}\text{/}$. There are four main dispersion bands in biological tissue: $\alpha$ (from DC to about 10), $\beta$ (between 100 and 10), $\gamma$ (at 20), and $\delta$ (between $\beta$ and $\gamma$). Tissue electroporation uses square-wave PEF with a broad spectral distribution. For this reason, including biological dispersion is essential to accurately simulate the biological medium. Then a dynamic model of tissue electroporation should be developed based on biological dispersion [13, 14].

Proposals for tissue dynamic models have already been published [15, 16, 17, 18, 19, 20, 21, 22]. Some are for the analysis of irreversible electroporation, considering the dynamics of the temperature rise and its effects on the tissue conductivity. However, they do not consider time-dependent changes in electrical properties due to electroporation [15, 16, 17]. On the other hand, three models consider time-dependent changes in electrical properties due to electroporation [18, 19, 20, 21, 22]. Yet, there are some modelling simplifications. One did not consider the biological dispersion [18], while the other two consider the tissue dispersion to be $\beta$-based [19, 20, 21, 22]. Nevertheless, it is known that for most electroporation PEF protocols, more than 95% of the spectral energy is below 100 [13], which is in the $\alpha$-dispersion band. Although the three models [19, 20, 21, 22] can accurately describe the shape of the electric current during tissue electroporation, some of their parameters are adjusted according to the experimental setup, such as voltage and shape of the electrode. These adjustments hamper the use of models if the experimental setup is changed, which is common practice in PEF research. Thus, the input parameters should be redefined. Furthermore, the three models were developed using custom numerical solution software and not an immediate implementation.

In this paper, we present a novel dynamic model of tissue dielectric properties during PEF that accounts for the three individual physical effects: electroporation, tissue dispersion, and temperature. We describe tissue dispersion using a multipole Debye function implemented in the time domain by the method of auxiliary differential equations. The electroporation effect was described using an electric-field-to-TMP relation. The time-dependent increase in tissue conductivity was evaluated by extrapolating a kinect model of cell electroporation. In addition, the conductivity dependence on temperature was included. We used in vitro potato tuber (Solanum tuberosum) to collect data and implement the model using commercial finite element method (FEM) software.

II Methods

II.1 Model Development

We solve the models using numerical simulations. We used the commercial software COMSOL Multiphysics (COMSOL Inc., Stockholm, Sweden). COMSOL is both a FEM and a computational fluid dynamic (CFD) solver with a bundle of built-in physical equations. We used COMSOL’s electric current library, which solves electrophysics for low-frequency signals using the FEM. In low-frequency electrophysics, COMSOL use the principle of charge conservation and solve the equation of continuity shown in Eq. 1 and the electric current density given by Eq. 2.

$$\vec{\nabla}\cdot\vec{J}(t)=Q_{j}$$ (1)

$$\vec{J}(t)=\sigma\vec{E}(t)+\frac{\partial\vec{D}(t)}{\partial t}+\vec{J}_{e}(t)$$ (2)

where $Q_{j}$ is the total charge density ($\mathrm{C}\text{/}{\mathrm{m}}^{3}$) and $\vec{J}(t)$ is the electric current density ($\mathrm{A}\text{/}{\mathrm{m}}^{2}$) whose components are given by the material conductivity $\sigma$ ($\mathrm{S}\text{/}\mathrm{m}$), the electric field $\vec{E}$ ($\mathrm{V}\text{/}\mathrm{m}$) and the displacement field $\vec{D}$ ($\mathrm{C}\text{/}{\mathrm{m}}^{2}$). $\vec{J}_{e}$ is an arbitrary external electric current density assumed by COMSOL to allow inclusion of external effects. The Maxwell equations define the displacement field as shown in Eq. 3, where $\epsilon_{0}$ is the vacuum permittivity and $\varepsilon_{r}$ is the relative permittivity of the material.

$$\vec{D}(t)=\epsilon_{0}\varepsilon_{r}\vec{E}(t)$$ (3)

II.1.1 Biological Dispersion

A biological medium is a dielectric that has regions susceptible to polarisation and consequently to charge relaxation times. In the frequency domain, the relaxation effect is called dispersion and can be represented as a complex relative permittivity function. There are several models for describing the dispersion behaviour in biological materials. The Debye dispersion can be implemented in the time domain using a set of auxiliary differential equations [14]. Eq. 4 represents the multipole Debye model in the frequency domain.

$${\varepsilon}_{r}^{*}(\omega)=\varepsilon_{\infty}+\frac{\sigma_{s}}{j\omega\epsilon_{0}}+\sum_{k=1}^{N}\frac{\Delta\varepsilon_{k}}{1+\left(j\omega\tau_{k}\right)}$$ (4)

where $\sigma_{s}$ is the static conductivity, $\omega$ is the angular frequency ($\mathrm{rad}\text{/}\mathrm{s}$), $\varepsilon_{\infty}$ is the permittivity of the material at high frequency, $\Delta\varepsilon_{k}$ is the permittivity variation of the pole, and $\tau_{k}$ is the relaxation time of the pole ($\mathrm{s}\text{/}$). $k$ and $N$ are the current pole number and the total number of poles, respectively. $j$ is the imaginary unit.

The implementation of the multipole Debye model in the time domain followed our previously proposed method [14], in which the dispersive effect is contained in an external electric current density for each Debye pole ($\vec{J}_{e_{k}}$), as shown in Eq. 5. Each current density is solved using an auxiliary electric field ($\vec{e_{k}}$) as shown in Eq. 6, which is a delayed value of the input electric field ($\vec{E}$) with a time constant corresponding to the relaxation of the Debye pole as shown in Eq. 7. This set of equations is implemented in COMSOL Multiphysics with the domain ordinary differential equation (DODE) physics.

$$\vec{J}(t)=\sigma_{s}\vec{E}(t)+\epsilon_{0}\varepsilon_{\infty}\frac{\partial\vec{E}(t)}{\partial t}+\sum_{k=1}^{N}\vec{J}_{e_{k}}(t)$$ (5)

$$\vec{J}_{e_{k}}(t)=\frac{\epsilon_{0}\Delta\varepsilon_{k}}{\tau_{k}}\left(\vec{E}(t)-\vec{e_{k}}(t)\right)$$ (6)

$$\tau_{k}\frac{\partial\vec{e_{k}}(t)}{\partial t}=\vec{E}(t)-\vec{e_{k}}(t)$$ (7)

We have previously described the dielectric spectrum of Solanum tuberosum tissue with different numbers of Debye poles [14]. The dielectric dispersion from 40 to 10 can be parameterised with a 4-pole Debye dispersion model. The parameters are shown in Table 1.

Table 1: Parameterisation of potato tissue dispersion with the 4-pole Debye dispersion model [14].

Parameter	Value	Parameter	Value
$\epsilon_{\infty}$	$1.746\,734\text{\times}{10}^{02}$	$\tau_{2}$ ($\mathrm{s}\text{/}$)	$2.309\,457\text{\times}{10}^{-05}$
$\sigma_{s}$	$2.158\,834\text{\times}{10}^{-02}$	$\Delta\varepsilon_{3}$	$1.836\,403\text{\times}{10}^{04}$
$\Delta\varepsilon_{1}$	$2.251\,149\text{\times}{10}^{06}$	$\tau_{3}$	$1.005\,246\text{\times}{10}^{-06}$
$\tau_{1}$ ($\mathrm{s}\text{/}$)	$3.783\,432\text{\times}{10}^{-03}$	$\Delta\varepsilon_{4}$	$1.052\,989\text{\times}{10}^{04}$
$\Delta\varepsilon_{2}$	$2.917\,574\text{\times}{10}^{04}$	$\tau_{4}$ ($\mathrm{s}\text{/}$)	$1.658\,463\text{\times}{10}^{-07}$

II.1.2 Electroporation

Electroporation dynamics was described on the basis of a modified version of the Leguèbe et al. [23] cell model. We modelled membrane electroporation using three states: prepore ($P_{0}$), initial pore ($P_{1}$), and expanded pore ($P_{2}$). Each state contributes in a specific way to the increase in membrane conductivity. In the cell model, the pore state is related to the TMP. In a macro-tissue model, we cannot access TMP directly because we do not have the membrane geometric model. For this reason, we proposed to calculate the TMP based on the magnitude of the local electric field ($E$). The concentrations of the pore states grow and decay exponentially, $P_{0}$ and $P_{1}$ as a function of the electric field and $P_{2}$ as a function of $P_{1}$. The proposed pore formation system was described using a set of differential Eqs. 8 – 10. The notation in brackets indicates the concentration of the state.

$$\frac{\mathrm{d}[P_{0}]}{\mathrm{d}t}=\frac{\beta_{0}(E)-[P_{0}]}{\tau_{0}}$$ (8)

$$\frac{\mathrm{d}[P_{1}]}{\mathrm{d}t}=\frac{\beta_{1}(E)-[P_{1}]}{\tau_{1}}$$ (9)

$$\frac{\mathrm{d}[P_{2}]}{\mathrm{d}t}=\frac{[P_{1}]-[P_{2}]}{\tau_{2}}$$ (10)

$\tau_{0}$ is a constant time. $\tau_{1}$ and $\tau_{2}$ depend on whether the function is growing or declining, as represented in Eqs. 11 and 12. $\beta_{0}$ and $\beta_{1}$ describe the maximum concentration for the states $P_{0}$ and $P_{1}$ as a function of the magnitude of the electric field and are expressed in Eqs. 13 and 14, respectively. This means that both functions can vary between zero and one. Zero means that no prepore ($\beta_{0}$) or pore ($\beta_{1}$) is formed. One means that the tissue has reached saturation for each phenomenon. Note that saturation does not mean that no new prepores or pores are formed, but that their increasing number above a certain threshold no longer has significant influence on the conductivity of the tissue. The maximum number of pores is usually not considered in cell electroporation models. The asymptotic model [24], for example, does not define a limit value for the pore density.

$$\tau_{1}=\left\{\begin{matrix}\tau_{1G}(1-0.5[P_{0}])&\text{if}\ \beta_{1}(E)-[P_{1}]\geq 0\\ \tau_{1D}&\text{otherwise}\end{matrix}\right.$$ (11)

$$\tau_{2}=\left\{\begin{matrix}\tau_{2G}&\text{if}\ [P_{1}]-[P_{2}]\geq 0\\ \tau_{2D}&\text{otherwise}\end{matrix}\right.$$ (12)

$$\beta_{0}(E)=\frac{1}{1+e^{-(|E|-E_{0})/\Delta E_{0}}}$$ (13)

$$\beta_{1}(E)=\frac{1}{1+e^{-(|E|-E_{1})/\Delta E_{1}}}$$ (14)

where $E_{0}$ and $E_{1}$ are the central values of each logistic function, and $\Delta E_{0}$ and $\Delta E_{1}$ shape the slope. $\tau_{1G}$ and $\tau_{2G}$ are the characteristic times of $[P_{1}]$ and $[P_{2}]$. $\tau_{1D}$ and $\tau_{2D}$ are the relaxation times of $[P_{1}]$ and $[P_{2}]$.

The increase in the number of pores ($P_{1}$) would probably decrease the number of prepores ($P_{0}$); the same is true for increasing the size of the pores ($P_{2}$) and their initial shape ($P_{1}$). We have not included this state interaction in our model definition because we are working with concentrations and not absolute numbers. It is not expected that all prepores lead to a pore and that all pores expand. According to [25], about only 2.2% of the pores will expand. We believe that its implementation would have little effect on the final result of the model while increasing the number of parameters.

Since we cannot directly deal with membrane conductivity in tissue simulation, the increase in conductivity was implemented in tissue dispersion. There are works that deal with the effects of electroporation on biological dispersion [26, 27]. Impedance analysis before and after PEF stimulation showed an increase in tissue conductivity throughout the spectrum. Permittivity is also affected, but to a lesser extent. To increase conductivity across the spectrum, we increased the static conductivity of the biological dispersion ($\sigma_{s}$) according to Eq. 15.

$$\sigma_{P}=\sigma_{s}+\sigma_{P_{0}}[P_{0}]+\sigma_{P_{1}}[P_{1}]+\sigma_{P_{2}}[P_{2}]$$ (15)

where $\sigma_{P_{0}}$, $\sigma_{P_{1}}$, and $\sigma_{P_{2}}$ are the increasing coefficients of the states $P_{0}$, $P_{1}$, and $P_{2}$, respectively. $[P_{0}]$, $[P_{1}]$, $[P_{2}]$ are evaluated using Eqs. 8 – 10.

II.1.3 Thermal Dependence

The electrical conductivity of biological tissue increases with temperature [28]. The electric current flowing through the tissue generates Joule heating. We simulated the temperature development in the sample during the pulse burst. The equation for heat diffusion with Joule heating term is presented in Eq. 16.

$$\rho c_{p}\frac{\partial T}{\partial t}-\nabla\cdot\left(k\nabla T\right)=\vec{J}\cdot\vec{E}$$ (16)

where $T$ is the temperature ($\mathrm{K}\text{/}$), $\rho$ is the density of the material ($\mathrm{kg}\text{/}{\mathrm{m}}^{3}$), $c_{p}$ is the heat capacity of the material at constant pressure ($\mathrm{J}\text{/}\text{(}\mathrm{kg}\text{\,}\mathrm{K}\text{)}$), and $k$ is the thermal conductivity ($\mathrm{W}\text{/}\text{(}\mathrm{m}\text{\,}\mathrm{K}\text{)}$). $\vec{J}$ and $\vec{E}$ are the electric current density and electric field calculated with Eqs. 1 and 2. The thermophysical properties of Solanum tuberosum [29] and the electrode used during the experiments (316L Stainless Steel [30]) are given in Table 2.

Table 2: Thermophysical properties at 20 for the materials used in the experiments.

Material	$\rho$	$c_{p}$	$k$
Solanum tuberosum [29]	1053	4410	0.56
316L Stainless Steel [30]	8000	480	13.50

Temperature influences the increase in electroporation conductivity (Eq. 15) because it affects all components of the tissue. The conductivity temperature coefficient was defined as $\chi$ = $1.7\text{\times}{10}^{-2}$ $\text{/}\mathrm{K}$ [15]. Thus, the conductivity is adapted as follows.

$$\sigma_{T}=\sigma_{P}\left(1+\chi(T-T_{0})\right)$$ (17)

where $T_{0}$ is the initial temperature.

After including thermal and electroporation effects, Eq. 5 is adapted to Eq. 18 and used to describe the apparent conductivity in the simulator.

$$\vec{J}(t)=\sigma_{T}\vec{E}(t)+\epsilon_{0}\varepsilon_{\infty}\frac{\partial\vec{E}(t)}{\partial t}+\sum_{k=1}^{N}\vec{J}_{e_{k}}(t)$$ (18)

II.2 In vitro experiment

Potato tubers (Solanum tuberosum) were brought from local growers. Growers were certified by the Brazilian Ministry of Agriculture, Livestock and Food Supply (MAPA) for organically grown products. Cylindrical fragments were cut using a 18.50 diameter stainless steel cutter. The cylindrical fragments were then cut into 5 tall samples. The samples were wrapped in paper towels to reduce denaturation and oxidation. Cutting was performed immediately before each experiment. The time between cutting and the end of the experiment was less than 30 minutes. The laboratory temperature was 20 (293.15).

Figure 1: The experimental sample and the rotated 2D axisymmetric geometry created in the simulator. (a) Experimental setup. (b) 3D representation. (c) Transversal view. (d) Frontal view. The geometry is centred at the origin.

The samples were placed between two 30 mm diameter circular 316L stainless steel plates as shown in Fig. 1a and carefully fixed with a spring clamp. We subjected the samples to a PEF protocol according to the ESOPE guidelines [31]. The repetition rate was fixed at 5. Each sample was subjected to a PEF protocol and then replaced. The voltage was swept from 50 to 250 (50 steps) and from 300 to 500 (100 steps). Because the samples were 5 high, the equivalent electric field ranged from 10 to 50 (10 steps) and from 60 to 100 (20 steps).

Data were collected using a Tektronix DPO2012B oscilloscope (Tektronix Inc, Oregon, USA) with a Tektronix TPP0100 voltage probe and a Tektronix A622 current probe. We post-processed the data using a Python script to determine the average, standard deviation, and confidence intervals for each protocol.

II.3 In silico experiment

Iterative simulations were used to evaluate the parameters. We used a 2D axisymmetric geometry to replicate the experimental setup. The sample geometry was a rectangle 9.25 long and 5 high. Rectangles 15.00 long and 1 high were placed on the top and bottom of the sample geometry to form the plate electrodes. All geometries were rotated 360 to form the cylindrical shape. Fig. 1b shows the final geometry. We implemented the biological dispersion (Eqs. 5 – 7) with electroporation and thermal conductivities (Eqs. 17 and 18) in the sample material. The conductivity and relative permittivity of the electrode were 1.74 and 1, respectively.

The boundary conditions were defined as follows. A lateral boundary of the entire geometry was used for the axisymmetric rotation. For electrical analysis, the boundaries of one electrode were defined as terminal and those of the other as ground (Dirichlet boundary condition). For electrical and thermal analysis, the external boundaries of the geometry were defined as electrical and thermal insulating, respectively (Neumann boundary condition). The COMSOL’s multiphysics module automatically provided the electrical information as input for thermal analysis, all domains were considered as Joule heating source. The initial temperature was 20.

The input voltage followed the experimental signal. The mesh was created with the COMSOL mesh creation tool using finer resolution. The final mesh resulted in 736 domain elements. We used the intermediate generalised-$\alpha$ method to adjust the time step to improve convergence. The intermediate generalised-$\alpha$ method allows the solver to strictly decrease the time step. The maximum time step was established at 0.1 in the transition regions and at 1 otherwise.

III Results

The parameters of the dynamic model of Solanum tuberosum electroporation are shown in Table 3. Fig. 2 presents the plot of the functions $\beta_{0}$ and $\beta_{1}$. Fig. 3 shows the experimental and simulated electric currents for the used PEF ($10$ to 100). The experimental results were summarised as average, confidence interval (95%), and standard deviation. If PEF is higher than 20, the electric current has a nonlinear increase during PEF. Also, if PEF increases in magnitude, the maximum current values increase nonlinearly. The overshoot of the in vitro current in the opposite direction in the PEF transitions is not due to dispersion, temperature rise, or electroporation, but is a common parasitic effect in the experimental setup.

Table 3: Electroporation parameters for potato tubers. Initial temperature was 20.

Parameter	Value	Parameter	Value
$E_{0}$	43 $\mathrm{kV}\text{/}\mathrm{m}$	$\tau_{1G}$	40 $\mathrm{\SIUnitSymbolMicro s}\text{/}$
$\Delta E_{0}$	5.5 $\mathrm{kV}\text{/}\mathrm{m}$	$\tau_{1D}$	150 $\mathrm{\SIUnitSymbolMicro s}\text{/}$
$E_{1}$	22 $\mathrm{kV}\text{/}\mathrm{m}$	$\sigma_{P_{1}}$	0.11 $\mathrm{S}\text{/}\mathrm{m}$
$\Delta E_{1}$	2.7 $\mathrm{kV}\text{/}\mathrm{m}$	$\tau_{2G}$	500 $\mathrm{\SIUnitSymbolMicro s}\text{/}$
$\tau_{0}$	0.5 $\mathrm{\SIUnitSymbolMicro s}\text{/}$	$\tau_{2D}$	1 $\mathrm{s}\text{/}$
$\sigma_{P_{0}}$	0.375 $\mathrm{S}\text{/}\mathrm{m}$	$\sigma_{P_{2}}$	0.04 $\mathrm{S}\text{/}\mathrm{m}$

Figure 2: Shape of sigmoidal functions $\beta_{0}$ and $\beta_{1}$ for potato tissue with parameters of Table 3.

Figure 3: Experimental and simulation results for the electroporation experiments following the ESOPE guidelines (eight pulses 100 long at 5). Only the first, second, and eighth pulses are shown; see Supplementary Information for the full set of pulses. Magnitude of the PEF protocol at (a) 10, (b) 20, (c) 30, (d) 40, (e) 50, (f) 60, (g) 80, and (h) 100. The overshoot of the in vitro current in the opposite direction in the PEF transitions is a common parasitic effect in the experimental setup. The circled numbers indicate the pulse number. Exp Avg is the experimental average, CI is the confidence interval (95%) and Std Dev is the standard deviation.

We evaluated the evolution of dynamic states and the thermal increase in the centre of the sample (coordinates (0;0) in Fig. 1c and d). The dynamic evolution of the concentrations of prepore ($P_{0}$), initial pore ($P_{1}$), and expanded pore ($P_{2}$) using 20 and 100 is shown in Fig. 4 (see Supplementary Information for the complete set of input electric fields). $P_{0}$ are created and decimate during the first 500 after the pulse rise and fall times. $P_{1}$ are created at a rate faster than its decimation. A higher magnitude of the PEF leads to more pore formation. $P_{2}$ has the slower dynamics and accumulates over pulses.

Figure 4: Concentration of electroporation dynamic states $P_{0}$, $P_{1}$, and $P_{2}$ for 20 and 100 at the centre of the sample.

The increase in temperature is shown in Fig. 5a. The total increase in temperature ($\Delta T=T-T_{0}$) after the eight pulses was $0$, $0.004$, $0.025$, $0.074$, $0.1822$, $0.3046$, $0.5647$, and 0.8862 for $10$, $20$, $30$, $40$, $50$, $60$, $80$, and 100, respectively. The maximum temperature variation (0.8862) represents an increase in conductivity of 1.5%.

Figure 5: (a) Temperature rise due to Joule heating at the centre of the sample. (b) Apparent conductivity with thermal and electroporation influences at the centre of the sample. For the apparent conductivity, only the first, second, and eighth pulses are shown; see Supplementary Information for the full set of pulses.

Fig. 5b presents the evolution of the apparent conductivity for all input PEF. Measurement was taken in the centre of the sample. This curve is calculated through Eq. 17, where the contribution of electroporation and thermal effects to tissue conductivity is included.

IV Discussion

Modelling complex biological systems is a challenging undertaking due to complex interactions, nonlinear dynamics, computational demand, and uncertainties inherent in these systems. Another difficulty is addressing parameters that correlate with the microphysical processes. Studies in single cells provide an opportunity to study the effects directly at the cell membrane [32]. On the other hand, tissue studies can only evaluate a macroscopic effect [33]. In electroporation, the increase in membrane conductivity (and thus the tissue conductivity) is one of the primary observable effects. Our tissue model subsumed the complex dynamics of electroporation into three main states $P_{0}$, $P_{1}$, and $P_{2}$, which correlate with the membrane-level hypothesis of pore creation and expansion. These states differ in characteristic and relaxation times, contribution to increase in conductivity, and dependence with the applied electric field [34, 35]. Although the model was built on the basis of the cell electroporation model of Leguèbe et al. [23] we proposed a different approach. Our model has an extra pore-detailing state (three states instead of two). In our definition, $P_{0}$ should explain the instantaneous increase in conductivity associated with rapid opening of hydrophobic pores (so-called prepores), $P_{1}$ should explain the initial formation of hydrophilic pores (initial pores), and $P_{2}$ the final expansion of the pores (final pores). Because of that, the relation between each state is also slightly different, $P_{0}$ and $P_{1}$ depend on the electric field, and $P_{2}$ depends directly on $P_{1}$.

We must note that the dynamics of electroporation is not yet fully understood [36], and there are several hypotheses on how the phenomenon occurs [25, 37, 38, 39, 40]. On the tissue scale, we can only assess a combined effect, which limits our conclusions about the physical changes on the membrane. Although we assume that each state is explained mainly by the aforementioned reasons, one would expect pores to form during $P_{2}$ and pores to expand during $P_{1}$, for example. Further work should be done to validate the correlations between tissue and membrane effects for each state.

The differences between the curves in Fig. 2 show that the saturation of $\beta_{1}$ is reached at electric field magnitudes of about half the saturation threshold of $\beta_{0}$. We believe that $\beta_{1}$ is saturated at lower electric field magnitudes due to the spatial distribution in the cell membrane for prepore and pore formation. Prepores are usually less stable and require less energy to form. Pores, on the other hand, will require more energy to form. For this reason, pores are likely to occur in the poles of the cell, where the electric field perpendicularly strikes its structure [25, 41]. In low electric field stimulation, prepores and pores would likely form on the cell poles. Under these conditions, the rate of pore formation would be higher and a small amount of prepores would almost saturate the number of initial pores. Increasing the electric field magnitude would increase the number of prepores along the cell structure, but those would not form an initial pore. The difference in energy threshold would hamper the formation of pores throughout the cell structure. This could explain why increasing the magnitude of the electric field after a certain threshold keeps increasing the conductivity at the beginning of the pulse but does not change the increase amount over the pulse. As mentioned, the mechanics of single cell electroporation is not yet fully understood, and conclusions about its causes on the tissue scale are speculative.

The final pores cause larger flaws in the membrane than pre-pores or initial pores. Since the membrane is an insulating material, the flaw size leads the pre-pores to be less conductive than the initial and final pores. In fact, studies on cell suspensions under electroporation indicate an increase in conductivity over time during PEF stimulation, consistent with the higher conductivity of the final pores [42]. There is a lack of information on conductivity during initial pore formation. Pre-pores are expected to form spontaneously even if the cell is at resting potential [34]. However, our proposed model does not consider the absolute values, but the concentration of the individual states. This consideration explains why the conductivity of the pre-pores state ($\sigma_{P_{0}}$) is greater than that of the initial pores ($\sigma_{P_{1}}$) and the final pores ($\sigma_{P_{2}}$). In absolute terms, we expect the number of pre-pores to be greater than the number of pores and final pores ($P_{0}\gg P_{1}\gg P_{2}$), but the concentration analysis normalises these values. Therefore, the average increase in conductivity is reflected in a higher value for the conductivity coefficient of the pre-pores than for the initial and final pores ($\sigma_{P_{0}}\gg\sigma_{P_{1}}\gg\sigma_{P_{2}}$).

The higher number of prepores formed in the first moments of the pulse accelerates the formation of the initial pores. This effect could be linked to theories that pores may be formed and expanded by coalescence of prepores [43, 44]. Therefore, a higher concentration of $P_{0}$ would lead to a faster increase in $P_{1}$, as introduced in Eq. 11. $P_{1}$ and $P_{2}$ are also related. Here, we consider that the expansion of pores is proportional to the occurrence of initial pores. The difference is that the pores would take longer to expand and then longer to close. The timing of opening and closing of the pores should also be taken into account. There are differences in the mechanisms of opening and closing pores [45]. The mechanism of closing the pores is slower than that of opening them. Eqs. 11 and 12 adjust the characteristic or relaxation time of states $P_{1}$ and $P_{2}$ depending on whether it is a growing or a decaying curve.

Fig. 3 shows that our proposed model can describe the dynamics of the electric current during electroporation for tested input voltages. The solid foundation of the biological dispersion can be observed when 10 (50) is applied (Fig. 3a). At 10, the influences of electroporation are small, so the electric current is explained mainly by biological dispersion (see that there is no increase in apparent conductivity in Fig. 5b). Electroporation phenomena start to visually occur above 20 (Fig. 3b), when the electric current deviates from the natural waveform of the biological dispersion. This threshold for the occurrence of electroporation is assessed in the curves of Fig. 2, especially in the first dynamic dependence $\beta_{0}$. The $\beta_{0}$ shape is similar to the static model of Solanum tuberosum proposed by Ivorra et al. [46]. The authors have developed a static model applying a single pulse 400 long and evaluating instantaneous conductivity at 100. The similarity between $\beta_{0}$ and the static model is consistent, as $P_{0}$ has the greatest influence on the increase in tissue conductivity for the most magnitudes.

The results of the thermal analysis in Fig. 5a show that the ESOPE protocol has a minimal temperature effect even at the higher electric field. The maximum increase was only 0.88, which is reflected in a change in conductivity of approximately 1.5. Although we did not perform an experimental analysis of the thermal rise, our thermal simulation results are similar to those of other studies for the first eight pulses [15, 16, 17].

Fig. 4 explains the electroporation dynamics for 20 and 100. We can see that $P_{1}$ and $P_{2}$ concentrations increase from 20. As we defined that $P_{1}$ and $P_{2}$ are the states of hidrophylic pore formation, these states are expected to have the main influence on the increase in tissue permeability for macromolecules. Thus, cell permeability has already increased significantly at initial thresholds, where the effects of electroporation are not completely saturated. This confirms the statements in preclinical electrochemotherapy simulations that drug delivery is achieved with similar reliability after a certain threshold value of the electric field (the so-called reversible electroporation threshold) [47, 48, 49, 50, 51]. In further research, our aim is to better understand the effects of dynamic states on tissue permeability.

The increase in apparent conductivity with the influence of thermal and pore formation dynamics is shown in Fig. 5b. We can observe that $P_{1}$ and $P_{2}$ are more influential up to 30, while $P_{0}$ significantly influence the increase in conductivity for thresholds higher than 40. As previously mentioned, this phenomenon may be related to the formation of pre-pores throughout the cellular structure that do not have sufficient energy to form a final pore under higher-intensity stimuli. Thus, the states $P_{1}$ and $P_{2}$ determine the increase in conductivity for thresholds below 30, while $P_{0}$ has a greater influence on stimuli above 40.

We found characteristic times for the first two dynamics similar to Voyer et al. [19]. The authors analysed the electric current during the first pulse and used a similar set of equations to describe the first two dynamics. As only one pulse was analysed, they do not implement the relaxation between pulses or the third dynamic, which limits their model to a single-step analysis. In terms of increased conductivity, our results follow the same magnitude as those found by the static model of Solanum tubersum by Ivorra et al. [46] and the dynamic model of Weinert et al. [22] evaluated in rabbit tissues. Weinert et al. compressed all electroporation dynamics into a single differential equation, which led to distinct dynamics of the increase in conductivity. We suspect that this compress resulted in parameter adjustments for voltage variations in their model (see Table 2 in [22]). A common factor among the three dynamic models of tissue electroporation proposed to this date is the adjustment of parameters based on input variations [18, 19, 20, 21, 22]. In this sense, our model can accurately describe a wide range of input voltages with a single set of parameters.

The largest differences between the experimental average and the simulation results occur at 20 (Fig. 3b) and 50 (Fig. 3e). The difference at 20 arises because the simulated increase in conductivity is slightly faster than that observed experimentally. We could implement a new function for the characteristic time of $P_{1}$ to account for this difference. However, this would introduce new parameters to fit only the first two pulses of an input electric field. For simplicity, we prefer not to make this consideration. The difference at 50 is due to an overestimation of $\beta_{0}$ for this particular electric field. At 50, we are in the transition region of $\beta_{0}$, where small deviations in parameter definition would lead to large differences in electric current. It would be possible to better fit the result for 50 while increasing the deviation for the others. Since the result for 50 is at the upper edge of the standard deviation and all other input electric fields are close to the experimental average, we considered this set of parameters as the best choice.

V Conclusion

We have proposed a novel dynamic model that describes tissue dielectric properties during PEF while accounting for electroporation, tissue dispersion, and temperature. The model divides the electroporation phenomenon into three dynamic states: prepore, initial pore, and final pore formation. The states are associated with microscopic effects on the membrane. The model can accurately describe the electric current during PEF. We believe that our proposed model can improve the study of PEF for electroporation-based applications.

Author Declarations

Conflict of Interest

The authors have no conflicts of interest to disclose.

Author Contributions

R.G.: Conceptualisation, Methodology, Formal analysis, Investigation, Data Curation, Writing – Original Draft, Writing – Review & Editing, Visualisation. D.L.L.S.A.: Methodology, Investigation, Data Curation, Writing – Review & Editing. J.R.S.: Formal analysis, Investigation, Writing – Review & Editing. G.B.P.: Resources, Supervision, Writing – Review & Editing. D.O.H.S.: Resources, Supervision, Project administration, Writing – Review & Editing.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References