Electroporation Modelling of Irregular Nucleated Cells Including Pore Radius Dynamics

: When high-amplitude, short-duration electric pulses are applied to cells the permeability of their membranes is increased. From the biological point of view, the phenomenon is quite well understood, however, it is important to develop accurate numerical models to investigate the electroporation effectiveness in terms of electrical, geometrical and physical parameters. To this aim, in this paper, we illustrate a spatio–temporal, non-linear, and dispersive multiphysics approach to study the electroporation in irregularly nucleated shaped cells. The model couples the Maxwell equations with the partial differential equation describing the creation and closure of pores as well as the evolution of the pore size. The dispersive properties of biological media and the irregular geometries of the membranes have been described using the multi-relaxation Debye-based relationship and the Gielis superformula, respectively. Numerical simulations highlight the importance to include in the model the spatial and temporal evolution of the pore radius. In fact, the obtained numerical results show signiﬁcant discrepancies between our model and the one in which the pore radius dynamics is negligible.


Introduction
Electroporation (EP) is a non-thermal electromagnetic phenomenon in which the cell membrane permeability increases when exposed to a high voltage pulsed electric field (PEF). As a result, exogenous molecules that usually cannot cross the membrane barrier can enter the cytoplasm and nucleoplasm. EP is used in clinical and medical practice for various applications such as electrochemotherapy, gene therapy, DNA vaccination and nonthermal ablation [1,2]. Electroporation combined with chemotherapeutic drugs was found to be a very effective treatment both for killing tumor cells in vitro and for treating cancer in humans [3]. The effectiveness of the electroporation process depends on various physical and chemical parameters, such as the molecular composition of the membrane, but above all, on the parameters of electric pulses such as duration, frequency, rise and fall time, amplitude and number of pulses [4,5]. Conventional methods of electroporation use pulses of micro-millisecond duration to electroporate the cells. In this case, only the plasma membrane is electroporated since the cytoplasm and the extracellular medium act like good conductors. PEF with duration from the sub-microsecond to the nanosecond induce the electroporation of the plasma membrane as well as the permeabilization of the intracellular membranes. Furthermore, the electroporation process efficiency also depends on the type of cell [6], the mutual interaction with other cells and the cell positioning with respect to the electrodes [7]. Although electroporation is widely used in medical practice, EP protocols are very empirical as well as the electroporation process is not completely known. However, the accuracy of the mathematical model adopted to study the electroporation process has a considerable impact on the experimental analysis and on the choice of the EP parameters. Several numerical models have been proposed in literature to study the EP. Most of them use simplified geometric shapes to describe biological membranes as well as consider cellular structures as non-dispersive media. Moreover, in order to simplify the numerical calculation and to reduce the computational efforts, other methods do not take into account the thinness of the membranes. Among the different models of EP, studies implementing the dynamics of the pore radius have been also presented even if they were based on a simple cell geometry [8,9]. In view of these weaknesses, the present study is aimed at the evaluation of the transmembrane voltage (TMV) and pore density in a nucleated irregular-shaped cell. To this aim, a non-linear and dispersive numerical algorithm was developed. The dielectric response of cellular structures was modeled by incorporating the multi-relaxation Debye-based equations. Furthermore, the irregular shape of the cell was modeled by using the Gielis superformula. The Maxwell equations were solved in conjunction with the asymptotic Smoluchowski equation and the partial differential equations describing the spatial and time evolution of the size of the pores. Using the developed numerical algorithm, the electroporation process has been evaluated exposing the biological cell to short and long nanosecond pulses. Thus a comparison among models in which the pore radius dynamics are taken into account and the one in which it is neglected has been performed.

Cell Geometry Model
As shown in Figure 1, the biological cell considered in the proposed study is the nucleated cancer-like cell. The cell shape is modelled by using the Gielis superformula [4,6,7,10,11]. The relevance of the Gielis superformula to model the cell membranes perimeter was previously highlighted in [4]. Contrary to the canonical spherical cell-based model, an unconventional behavior of irregularly shaped cells was observed. In fact, the TMV does not have a conventional cosine type behavior and moreover, non-homogeneous EP effects occur along the cell membrane perimeter. In particular, the radius vectors r 1 and r 2 describing the plasma and nuclear membranes, respectively, are given by: where θ ∈ [−π/2, π/2], m i , n i , i = 1, . . . , 4 and b 1,2 are positive real numbers, d j , j = 1, . . . , 4 are strictly positive real numbers, A 1,2 and B 1,2 are appropriate scaling factors. The plasma membrane boundary of the selected cell can be well-approximated by using the superformula parameters A 1 = B 1 = 16.8 µm, m 1 = m 2 = 6, d 1 = d 2 = 1, n 1 = n 2 = 1 and b 1 = −2 as well as the nuclear membrane can be modelled by using the following Gielis parameters A 1 = B 1 = 3 µm, m 1 = m 2 = 2, d 1 = d 2 = 1, n 1 = n 2 = 4 and b 2 = 1. Moreover, to generate the cell membrane thinness, a particular extrusion operator transformation has been suitably implemented. In particular, the following equations was implemented where (x i , y i ) are the coordinates of the internal side of the membrane, (x e , y e ) are the coordinates of the external side of the membrane, h is the membrane thickness.

Pores Radius Model
The EP has been attributed to the creation of transient pores within the lipid bi-layer of the cell membrane when an external PEF is applied. This dynamic process can be modelled by Krassowska's asymptotic equation [12]. In particular, the temporal evolutions of pore densities for plasma N Pm and nuclear N Nm membranes are calculated by solving the following first-order partial differential equations: where N 0 is the pore density in the resting membrane, α and q are EP parameters, V ep is the characteristic voltage of EP, TMV Pm and TMV Nm are the TMV evaluated on plasma and nuclear membrane, respectively. Pore density is proportional to the amplitude of the applied electric field as well as the pores creation enhances the membrane conductivity. In particular, the membrane conductivity increase that accompanies electroporation can be computed as follows: σ Pm (x, y, t) = σ 0,Pm + K Pm N Pm (x, y, t)σ p,Pm πr 2 σ Nm (x, y, t) = σ 0,Nm + K Nm N Nm (x, y, t)σ p,Nm πr 2 p,Nm , where σ 0,Pm and σ 0,Nm are the membrane conductivity at rest, σ p,Pm and σ p,Nm are the conductivities of the medium inside the pore and r p,Pm and r p,Nm the pore radius, respectively, for plasma and nuclear membrane, and where w 0 is the energy barrier inside the pore and η is the relative length of pore entrance area with where q e is the electron electric charge and k the Boltzmann constant.
The density and size of the pores depend on the intensity of the applied electric field. In stationary conditions, a dynamic equilibrium among the steric repulsion of lipid heads, the line tension acting on the pore perimeter and the surface tension of the cell membrane occurs [13]. When an external electric field is applied, this balance is lost and the pores radius start to increase. The pores initially created with radius r p,Pm (for plasma membrane) and r p,Nm (for nuclear membrane), change size to minimize the energy of the entire lipid bilayer. So, the temporal evolution of pore radius for the plasma membrane is calculated by solving the following differential equation [8]: with where T e,Pm is the effective tension of the membrane: Similarly for the nuclear membrane, the temporal evolution of pore radius is determined as follow: where T e,Nm is the effective tension of the membrane: In Equations (16) and (22), S 1,Pm and S 1,Nm account for the electric force induced by the local TMV, S 2,Pm and S 2,Nm for the steric repulsion of lipid heads, S 3,Pm and S 3,Nm for the line tension acting on the pore perimeter, S 4,Pm and S 4,Nm for the surface tension of the cell membranes [8]. This study assumes that changes of cell area, volume, and shape can be ignored for pulses duration considered here. In particular, the presented model evaluating the temporal evolutions of pore densities and pore radius for plasma and nuclear membranes makes possible a quantitative evaluation of the membranes permeability to specific molecules and ions.

Complex Permittivity Model
The interaction of the PEF with the dielectric medium of cell compartment can be investigated by considering the electric polarization, i.e., the electric field induced disturbance of the charge distribution. From the macroscopic point of view, such interaction is generally due to different phenomena as the orientation of dipoles, ionic diffusion, interfacial polarization etc. Dielectric dispersion is the corresponding frequency dependence of permittivity which typically display extremely high dielectric constants at low frequencies, falling off in more or less distinct steps as the excitation frequency is increased. Cell suspensions typically exhibit a significant β-dispersion in the radio frequency range. This is due to the Maxwell-Wagner effect at the interface between the intra or extracellular solution and the phospholipid membrane. They also display a large α-dispersion at low frequency caused by the sarcoplasmic reticulum, gap junctions, and counterion relaxation on the cell surface. Moreover, an additive β-dispersion coming from the dispersion of organelles appears as a small tail to the large cell membrane β-dispersion [14]. In our numerical analysis, the following Debye-based relationship has been used to model the dispersive behavior of each cellular compartment [4,7,10].
where M is the order of Debye equation, ε ∞ is the high frequency permittivity, ∆ε i are the relaxation amplitudes and τ i are the relaxation times. As the dielectric relaxation of the plasma (Pm) and nuclear (Nm) membrane takes place in the MHz region, a second order Debye equation has been used to describe their dispersive properties. Instead, the dielectric relaxation of the extracellular medium (Ex), the cytoplasm (Cp) and the nucleoplasm (N p) occurs in the GHz zone and thus a first order Debye equation has been adopted to model the dispersive properties of these materials. The time domain relations pertaining the polarization vector P 1 and P 2 , corresponding to the first and the second order Debye dispersion model, can be formulated as: where and a 1,Nm = τ 1,Nm + τ 2,Nm (39)

Electromagnetic Model
The electromagnetic model is based on the Maxwell equations In particular, taking into account the small cell size and the microsecond/nanosecond electric pulses, the rate of the magnetic flux can be neglected. So, the resulting Maxwell equations in the differential form can be written in terms of the electric potential φ and the electric field E as Finally, the TMV for the Pm and Nm membrane is evaluated as the difference between the electric potential inside (i) and outside (o) of each membrane:

Results
Implementing the aforesaid computational algorithm, the EP process was studied in both the cases characterized by variable and constant pore radius. The simulations were carried out considering the nucleated cancer cell inserted into a 2D space domain having length L = 100 µm and width W = 100 µm (see Figure 1). Table 1 shows the geometric, electric and EP parameters used in the numerical computations. Our computations and plots have been performed using MATLAB suite software. In our calculations, the nuclear envelope consists of two membranes in near contact between them, each one with a thickness equal to half the thickness of the entire nuclear envelope. Such hypothesis can be justified taking into account that the voltage drop across the perinuclear space is negligible [6]. In particular, both the inner and outer membranes of the nuclear envelope have the same conductivity σ 0,Nm = 1 × 10 −4 S m −1 . Thus, supposing that the transmembrane voltage is equally distributed between the two lipid membranes, the TMV through the single membrane is calculated as half of the voltage on the entire nuclear envelope [15]. Moreover, in the case of the numerical model with a constant pore radius an appropriate validation was carried out in [6] using the experimental results concerning real cancer cells electroporated by a unipolar rectangular pulse with duration of 100 µs. Firstly, the biological cell has been exposed to a PEF having rectangular shape with amplitude E = 100 kV cm −1 , duration T = 10 ns, rise time t r = 0.9 ns, and fall time t f = 0.9 ns. Moreover, a time delay of 5 ns and a computational window for the time coordinate of 20 ns was considered. Figure 2a-c highlight the dynamics versus the time of the TMV, pore density and conductivity of the membranes at the top of the cancer-like cell (θ = 90 • ), respectively, for the model with the pore radius variable, r v p , and the model with the pore radius constant and equal to r p = 0.8 nm, r c p . The application of the external pulse generates a rapid increase in the transmembrane voltage on the plasma and nuclear membrane. The resulting pore creation, as well as the increase of the membrane conductivity, lead to a decrease of the TMV reaching a negative peak value. For both models, the activation of electroporation occurs at the same time instant. The subsequent decrease of TMV is more rapid in the model with a variable pore radius compared to that with a constant pore radius. In fact, the enhancement of both pore density and radius creates a rapid increment of membranes conductivity generating a reduction of the transmembrane voltage. Moreover, Figure 2d displays the dynamics versus the time of pore radius at the top of the cancer-like cell for the plasma and the nuclear membrane. As the applied pulse is of short duration, the temporal dynamics of the pore radius is similar for the plasma and the nuclear membrane, increasing from 0.51 nm (minimum radius of hydrophilic pores) to a maximum value of about 0.64 nm. Therefore, given the temporal shortness of the applied pulse, the pore radius for both the plasma and the nuclear membrane does not reach the value of r p0 = 0.8 nm which represents the minimum energy radius at TMV = 0 V. In addition, the pore radius size assumes equal values for the Pm and Nm membrane up to 5 ns, which denotes the time instant in which the applied pulse starts its activation. Moreover, an appreciable difference between the temporal evolution of the pore radius for each membrane is obtained starting by 10 ns. Furthermore, Figure 2e illustrates the pore density along the membranes circumference at t = 20 ns for both models. As highlighted by the obtained numerical results, the values of the pore density for the nuclear and the plasma membrane obtained using the model with pore radius variable are higher than that evaluated with the model having pore radius constant.  [15] For both models, the plasma and nuclear membrane are electroporated in the points facing the electrodes. In particular, the plasma membrane is electroporated between the top points and an angle of about 22 • and in the equatorial zone. Furthermore, the plasma membrane is not electroporated at θ = 0 • due to its irregular perimeter. Moreover, the irregular perimeter of the plasma membrane creates for both models peaks of pore density values at θ = 60 • . The nuclear membrane is characterized by significant electroporation between the top points and an angle of about 35 • . A more evident difference between the two models is further highlighted when the biological cell is exposed to longer pulses.   , respectively, for the model with the pore radius variable and the model with the pore radius constant. In these simulations, the biological cell has been electroporated using a rectangular pulse having amplitude E = 50 kV cm −1 , duration T = 100 ns, rise time t r = 0.1 ns and fall time t f = 0.1 ns. Moreover, the applied rectangular pulse starts at the time instant t start = 10 ns and the computational window for the time coordinate is 120 ns long. As reported in Figure 3a,b, the EP on the Nm membrane occurs about 1 ns faster than of Pm one. Furthermore, the TMV on each membrane quickly increases to about 1.6 V. The subsequent decrease of TMV is more rapid in the model with variable pore radius. In fact, for this model, the growing of the pore density and radius creates a rapid increment of membranes conductivity generating a quickly reduction of the TMV. Furthermore, as shown in Figure 3a the temporal dynamics of the TMV characterizing the two models are quite different. In fact, in the case of the model with variable pore radius, the temporal evolution of the TMV shows a significant gap between the peak value and the value that the curve assumes in the central flat area. Instead, for the model with a constant pore radius, the value assumed by the TMV in the central flat area is close to the value assumed by the peak. As shown in Figure 3d, the pore radius initially increases from r pm to a maximum value of about 1.1 nm, for the Pm membrane, and of about 1 nm, for the Nm one. In particular, for both membranes the maximum value of the pore radius is reached at the time instant of 110 ns. Moreover, the pore radius reaches the value of r p0 = 0.8 nm at the time instant of about 30 ns, for the plasma membrane, and at the time instant of about 40 ns, for the nuclear membrane. In addition, the pore radius size assumes equal values for the Pm and Nm membrane up to 10 ns, which is the pulse activation time instant. A relevant discrepancy between the temporal evolution of the pore radius of each membrane is appreciated starting by 20 ns. Finally, when the pulse is removed the pore radius starts to decrease. Moreover, Figure 3e illustrates the pore density along the membranes circumference at t = 120 ns for both models. The EP for the Pm membrane occurs in the angle range starting from 90 • and ending at about 25 • as well as in the equatorial zone. In particular, it is not electroporated at θ = 0 • due to its particular geometric shape. Moreover, the proper perimeter of the plasma membrane creates for both models peaks of pore density values at θ = 60 • . The nuclear membrane is characterized by significant electroporation between the top points and an angle of about 35 • . As highlighted by the obtained numerical results, a higher level of electroporation is obtained for the model in which the pore radius is assumed variable for both the plasma and nuclear membrane.

Conclusions
In this paper, a nonlinear dispersive numerical algorithm of electroporation for irregularly nucleated shaped cells including the variation of pore radius has been presented. The developed algorithm solves simultaneously the Maxwell equations, the Smolouchouski equation and the equation describing the dynamic of the pore radii evolution. In particular, the presented algorithm describes the irregular perimeter of the Pm and Nm membranes using the Gielis superformula and it considers the dielectric dispersion properties of cell compartments using a Debye formulation. With reference to a nucleated cancer-like cell, a number of simulations have been carried out in order to compare the model with the pore radius variable and the one characterized by a constant pore radius. The performed analysis highlights a discrepancy between the two models underlining the importance to take into account in the numerical algorithm the differential equation describing the variation of the pore size.

Conflicts of Interest:
The authors declare no conflict of interest.