A Many-Shells Model for Cell Transmembrane Potentials for Pulsed Electric Field Applications

This article advances the existing theoretical analysis of pulsed electric field (PEF) treatment, an application of pulsed power technology. PEF treatment has attracted significant attention due to its potential to be used for non-thermal biological sterilization and inactivation of microorganisms, bio-extraction, and for the possible role that it may play in the treatment of tumors and cancers such as in electrochemotherapy. However, the bioelectric effects of impulsive electric fields on different biological matter are not yet completely understood. Further advances in this direction would aid in the optimization of the necessary pulse waveforms to achieve the desired PEF effects in, for example, various biomedical and food processing industries. In this work, the commonly used multishell model of microorganisms utilized for the analysis of cell transmembrane potentials (TMPs) has been generalized to include an arbitrary number of layers. This incorporates also a spheroidal geometry of arbitrary eccentricity. Analysis has been conducted on the novel mathematical model, which demonstrated the ability to estimate TMPs and relaxation times for an n-shell topology using a mesh-free approach. This allows for complex many-shelled cell models to be analyzed, free from the limitations of spatial or temporal discretization, that is, those present when using finite-element (FEM) or finite-volume-based methods. Using this model, the effects that the pulse rise time and pulse duration have on the developed TMPs in a single six-layer Saccharomyces cerevisiae, and under the microsecond-PEF regime, have been investigated. It has been found that the pulse rise time has a far lesser effect on the PEF action compared to the modulation of the pulse duration, supporting past experimental observations. The role of the electrical conductivity of the extracellular medium (and considering induced changes in the cell components) has additionally been studied, under low ( $\sigma =1$ mS/m), medium ( $\sigma =50$ mS/m), and high ( $\sigma =100$ mS/m) conductivity values. Estimations provided by this model may support the optimization of pulse waveforms for current and future PEF applications, and for the analysis of complex multilayered cell structures under pulsed electrical stress.

Abstract-This article advances the existing theoretical analysis of pulsed electric field (PEF) treatment, an application of pulsed power technology.PEF treatment has attracted significant attention due to its potential to be used for non-thermal biological sterilization and inactivation of microorganisms, bio-extraction, and for the possible role that it may play in the treatment of tumors and cancers such as in electrochemotherapy.However, the bioelectric effects of impulsive electric fields on different biological matter are not yet completely understood.Further advances in this direction would aid in the optimization of the necessary pulse waveforms to achieve the desired PEF effects in, for example, various biomedical and food processing industries.In this work, the commonly used multishell model of microorganisms utilized for the analysis of cell transmembrane potentials (TMPs) has been generalized to include an arbitrary number of layers.This incorporates also a spheroidal geometry of arbitrary eccentricity.Analysis has been conducted on the novel mathematical model, which demonstrated the ability to estimate TMPs and relaxation times for an n-shell topology using a mesh-free approach.This allows for complex many-shelled cell models to be analyzed, free from the limitations of spatial or temporal discretization, that is, those present when using finiteelement (FEM) or finite-volume-based methods.Using this model, the effects that the pulse rise time and pulse duration have on the developed TMPs in a single six-layer Saccharomyces cerevisiae, and under the microsecond-PEF regime, have been investigated.It has been found that the pulse rise time has a far lesser effect on the PEF action compared to the modulation of the pulse duration, supporting past experimental observations.The role of the electrical conductivity of the extracellular medium (and considering induced changes in the cell components) has additionally been studied, under low (σ = 1 mS/m), medium (σ = 50 mS/m), and high (σ = 100 mS/m) conductivity values.Estimations provided by this model may support the optimization of pulse waveforms for current and future PEF applications, and for the analysis of complex multilayered cell structures under pulsed electrical stress.

I. INTRODUCTION
P ULSED electric field (PEF) technology is a well-established application of pulsed power that has gained significant research attention during the past several decades.The process of PEF treatment involves the application of short and intense electrical impulses to a target medium, which can produce different bioelectric effects on various biological matter [1].Namely, electroporation can take place, where nanopores are developed inside cell membranes [2], [3], [4].Two modes of electroporation are known to be possible: reversible electroporation (RE) and irreversible electroporation (IRE).The former allows the membrane to self-repair upon removal of the pulse, which has implications for drug delivery and for intracellular manipulation.In the latter case, the membrane is unable to be repaired and may thus induce cell death, enabling its effective application to bacterial inactivation, sterilization, and decontamination [5], [6], or biomass processing technologies [7], [8].In other work, nanosecond-PEF (ns-PEF) showed promise for the treatment of tumors and cancers [9], [10], [11], while emerging studies on picosecond-PEF (ps-PEF) have suggested that precise and noninvasive pulse delivery is possible [12], which may have impact on future treatment methods for neurodegenerative diseases [13].It is believed that, to effectively modulate the strength of the electroporation process, the transmembrane potential (TMP) developed across the cell membrane is of critical importance.Studies describe a threshold TMP in the range of ∼0.5-1.5 V for IRE to take place [14], [15], [16], and results from the detailed numerical modeling of membrane pore formation generally corroborate these findings, for example, in [17] and [18].It is known that the characteristics of the pulse can have major effects on the electroporation process [19], [20], [21], [22].As such, there have been multiple studies focused on modeling the developed TMPs under a variety of pulsed waveforms, and using a variety of methods, see for instance [23], [24], [25], [26].For a comprehensive review on the field of bioelectrics, which includes the above applications and other related phenomena, readers are referred to [27].
A class of multishell analytical dielectric models of biological cells or microorganisms have been commonly employed, 0093-3813 © 2024 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
with authors incorporating a varying number of layers, for example, in [19], [23], [26], and [28], representing various combinations (where applicable) of the nucleoplasm, nuclear membrane, cytoplasm, plasma membrane, cell wall, and extracellular medium.Maxwell-Wagner theory is applied to solve the time-dependent potential field in each dielectric layer, thereby estimating TMPs for single, isolated, cells.Past works including [25], [29] have extended the traditional spherical model to more complex geometries, such as spheroids, and in different orientations.These geometries may be more representative of certain types of microbiological cell which would be poorly approximated by an ideal sphere.In other literature, closed-form solutions have been demonstrated for up to five unique regions for spherically symmetric geometries [26].Often, mesh-based methods such as the finite-element (FEM) method have been used for models which must incorporate additional layers, for instance, when modeling non-mammalian cells incorporating the cell wall [30].However, due to the thickness of membrane layers relative to the greater thicknesses of the cytoplasm or extracellular space, very fine meshes are typically necessary to attain sufficiently accurate results.Even with adaptive methods, this may necessitate significant computational resources and time.Mesh convergence studies would typically be required to ensure convergent results.When also considering the inherent time-dependency of the problem, depending on the relaxation times and nature of the voltage pulses under consideration, the time-step requirement may become prohibitively small to ensure minimal truncation error.In combination with the limitations to the mesh above, these methods may take considerable time and require many problem iterations to ensure numerical accuracy.
With particular significance to applications such as bioextraction and food processing that may deal with microorganisms with more complex cell structures than simple bacteria, this work generalizes the analytical multishell approach to include an arbitrary number of shells.The geometry, which has been considered, is that of a spheroid, with its major axis parallel to an external, uniform, electric field.It should be emphasized, however, that the effects of eccentricity are not the main focus of this approach.Rather, the main contribution is the extension of the multilayered approach to an arbitrary number of layers.The effects of eccentricity have previously been studied, for example, in [31].It is first demonstrated that the model used in [26] represents the maximum theoretical limit to the number of layers where closed-form time-domain solutions are possible using the standard method.However, a new approach to estimate the TMPs in any layer in an n-layer geometry has been developed, which relies only upon the numerical computation of polynomial roots, bypassing some limitations of traditional mesh-based transient field analysis.This "many-shells" model was first validated against a FEM solver, before a systematic study on the developed TMPs across a complex multilayered cell geometry was conducted, which investigated the role of rise time, pulse duration, and effects of the extracellular medium conductivity.

II. MANY-SHELLS MODEL A. Nomenclature and Geometry
In this section, the development of the mathematical model is described.An isolated multilayered spheroidal inclusion was considered, which was assumed to have an arbitrary number of "shells," each with a unique pair of values for its relative permittivity, ε, and electrical conductivity, σ , in S/m, as illustrated in Fig. 1.This work focused on spheroidal geometries that have their major axis aligned to that of the direction of a uniform external electric field, E 0 (t).It should be noted that, in general, the term spheroidal refers to 3-D geometries formed by the revolution of an ellipse around one of its principal axes, and that, therefore, has a circular cross section in the plane of rotation.In some texts, the term ellipsoid of revolution is equally valid.In this article, the term spherical is used to explicitly refer to spheroids, which have all three principal axes of equal length.It is remarked that ellipsoids-shapes that may have differing lengths on all three principal axesdo not feature in this study.Furthermore, the use of the term shell, in this article, refers to a bounded region between two interfaces that identify a single unique layer in the composite layered structure, including neither the innermost bounded region nor the unbounded external medium.The term layer is instead used when referring to any bounded or unbounded region, including the innermost and outermost subdomains.This is made clear here as there is no accepted meaning of the term shell, see for instance, the differing definitions in [26], [28], and [32].
In this case, the geometry of interest was represented using the prolate spheroidal coordinate system (µ, ν, θ ), which relates to the Cartesian system (x, y, z) with x = a 0 cosh µ cos ν y = a 0 sinh µ sin ν sin θ z = a 0 sinh µ sin ν cos θ where a 0 is the distance from the origin to the foci.Surfaces traced by constant µ, therefore, form spheroids, and surfaces of constant ν are confocal hyperboloids.Analysis was limited to the case of azimuthal symmetry, where the θ coordinate can, therefore, be neglected.As such, let i denote the layer number, up to a total of n layers, with i = 1, therefore, denoting the innermost region, and i = n denoting the external medium.The coordinate µ = µ i , hence, represents the boundary between i and (i + 1)th layers, for a total of n − 1 unique µ i values for n layers.Furthermore, let e i be defined as the eccentricity, given by where K i is the ratio of the minor to major axis of the layer i, which may further be shown to satisfy Therefore, as e i tends toward zero, the geometry approaches the spherical limit.It should be noted that, because of ( 1) and ( 2), a limitation to this geometry is the inability to have layers of uniform thickness for all values of ν.As the number of layers n increases, so too does K i , such that for large n (or large distance between adjacent µ i relative to a 0 ), K n approaches unity (e = 0) regardless of the eccentricity of the innermost region.It is, therefore, important to provide a definition for the shell thickness, which was taken to be the difference between the Cartesian x coordinates of neighboring layers at an angle of ν = 0, and may, therefore, be written as for i from 1 to (n − 2).Consecutive values for µ i+1 can be found recursively using (3) when the dimension of the inner or outermost region is known.Each layer is also uniquely characterized by a constant value of relative permittivity, ε i , and electrical conductivity, σ i , in S/m.

B. Governing Equations and Boundary Conditions
Following Maxwell-Wagner theory, and under the assumption that charge exists only at the layer interfaces, the potential field ϕ in the domain is governed by the Laplace equation, which in the azimuthally symmetric prolate spheroidal coordinate system reads as follows: where h is the scale factor h = a 0 (cosh 2 µ − cos 2 ν) 1/2 .Equation ( 4) is separable, and under a uniform external field can be shown to admit solutions with the general form where subscript i refers to the potential field in layer i, and with the conditions that B 1 (t) = 0 and A n (t) = −a 0 E 0 (t), which arise from the necessity for the potential to be non-singular within the innermost region, and must become self-consistent with the external potential for µ ≫ µ n−1 .The function F(µ) depends only on geometry, and is given by The time-dependent coefficients A i (t) and B i (t) can be obtained by applying an appropriate set of boundary conditions.In this case, the continuity of the potential and of current must be prescribed, such that and where ∂ t denotes the time derivative, n is the unit normal, and ⃗ E i is the electric field developed in layer i, satisfying ⃗ E i = −∇ϕ i , and results in the expression (9) when combined with ( 5) The function G(µ) follows from F(µ), satisfying the relation and therefore, has the definition

C. Matrix Representation of the System
For an n-layer system, a total of 2(n − 1) first-order differential equations arise from (7) to (9).The system of equations can be reduced to a linear algebraic system by application of the Laplace transform, replacing the time derivative operator ∂ t by multiplication with the complex frequency variable s.For brevity, functional notation is, hereafter, omitted, and the substitution is made to maintain compactness of expressions.Similarly, F(µ i ), G(µ i ), and E 0 (s) become F i , G i , and E 0 , respectively.Here, τ 1 , . . ., τ n are the intrinsic relaxation time constants for the respective layer, defined as the ratio The set of equations for i from 1 to (n − 1) then reads Equation ( 13) was transformed into a matrix representation of the form M⃗ u = ⃗ b, which facilitated the solution process by revealing the underlying mathematical structure of the system.Here, ⃗ u is the vector of unknown coefficients, It was found possible to transform the system matrix into tridiagonal form with a regular structure, shown in ( 14) at the bottom of the next page, and the vector ⃗ b in (15), as shown at the bottom of the next page, which significantly simplifies the solution process, detailed in the next section.

A. s-Domain Solutions
Tridiagonal systems are characterized by their ability to be solved using a set of repeated elimination row operations, followed by backward substitution.This is also to say that should either A 1 or B n become known, all other coefficients can be found through (13) by propagating the known solutions backward (or forward) through each equation in sequence.One direct and efficient method that can be used for tridiagonal systems is the Thomas algorithm [33], which when applied to ( 14) and ( 15) yields a solution to A 1 for a structure with n layers, here indicated by an additional (n) superscript where the denominator P n is the characteristic polynomial of the system, the roots of which define the relaxation characteristics of the multilayered structure, and is a critically important property when considering the application of timedependent electric fields.From the row operations, it can be shown that P n satisfies the following second-order recurrence relation: for n ≥ 3 and with the conditions that The polynomial P n can be shown to be of order n − 1, and is of the form where c i are constant coefficients.The general s-domain solution may then be written as where the additional n − 1 time constants are found from the roots of ( 19), and hence Equation ( 20) allows closed-form solutions to be found for A (n) 1 in the s-domain, and consequently allows all coefficients in ⃗ u to be determined via backward substitution into (13), with the knowledge that B (n) 1 is always zero, as per the non-singular condition described in Section II-B.As for time-domain solutions, the ability for these to be derived analytically or in closed-form is discussed in Section III-B.

B. Time-Domain Solutions
The existence of analytical time-domain solutions to ( 4) is dependent on the ability to perform the inverse Laplace transform on the resulting s-domain coefficient expressions as given by (20).To do so would require the factorization of the characteristic polynomial, P n , of order n − 1 in s as given by (21), and to define an appropriate choice of external field, E 0 (s).As shown similarly and explained in [34] for a 1-D multilayered planar model, this prohibits closed-form time-domain solutions to be sought for n > 5 using this Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
method (since, in general, roots of polynomials with order n > 5 cannot be found in radicals, in accordance with the Abel-Ruffini theorem), for which numerical root finding techniques can be employed as an alternative (e.g., using the eigenvalue theorem).It is remarked, however, that even if (19) could be symbolically factored for n > 5, the size of the resulting expressions would become increasingly impractical and unwieldy, as the number of terms in P n according to (17) can be shown to grow by a factor of four for each additional layer.Equations ( 13) and ( 20) nonetheless provide a full analytical solution to the problem, which can be processed by common computer algebra systems, provide an exact solution, and are dependent neither on a full computational mesh nor on a discrete time-stepping.Provided that ( 21) can be solved accurately for the relaxation time constants, ( 13) and ( 20) can be used to compute the electric potential and electric field anywhere in the domain, and at any time, free from the limitations of a spatial and temporal discretization.

IV. SUMMARY OF MODELING PROCESS
For clarity, this brief section provides a step-by-step description of the modeling process based on the derived general solution.The region with available closed-form solutions (n ≤ 5) where numerical root-finding is not strictly necessary is also discussed.It is remarked that numerical methods for computing roots may also be applied for n ≤ 5 despite the existence of closed-form solutions, if desired.Regardless of n, the process can be described in the following steps.
2) Compute the n − 1 layer relaxation times based on the polynomial roots using (21).This may be done analytically or numerically (see below after this list).3) Evaluate (20) using the obtained relaxation times and the desired external field expression, E 0 (s).This results in the s-domain solution for the coefficient A (n) 1 .4) Perform an inverse Laplace transform on the obtained solution from step 3), yielding the time-domain solution A (n) 1 (t).5) Repeated backpropagation of A (n)  i (t) using the linear system (13) allows the determination of all required unknown coefficients.6) With the obtained coefficients, evaluation of (5) yields the potential field, while (9) yields the electric field developed in the ith layer.
Step 2) involves the computation of the relaxation times, which as discussed in Section III-B is generally limited to n ≤ 5 if closed-form solutions are necessary, as limited by the Abel-Ruffini theorem.In those cases, numerical root-finding schemes should be used in Step 2) to determine the polynomial roots.Otherwise, closed-form solutions may be sought analytically as demonstrated in the following.
For n = 1, this corresponds to the absence of any cell or particle and represents solely the application of a uniform Fig. 2. Comparison of inner and outer TMPs developed across a spherical cell with data from [31], [36] with the present many-shells model with n = 5.Note that all models used are analytical and produce three essentially identical curves.Discrete datapoints have, therefore, been selected from [31], [36] and plotted as an overlay for visibility.
external field to a single bulk dielectric region.By definition of ( 18), P 1 is equal to unity, and the evaluation of ( 20) is, therefore, A (n) 1 = −a 0 E 0 (s).Substitution into (5) recovers the expression for the uniform external field in prolate spheroidal coordinates, as expected.
For n = 2, the polynomial is P 2 = F 1 − G 1 λ 1 as per (18), such that (21) can be solved through direct algebraic rearrangement.For the cases n = 3, 4, 5, the characteristic polynomial has order 2, 3, and 4, respectively.The roots may, therefore, be computed following the quadratic, cubic, or quartic formulas without the necessity of numerical rootfinding techniques.

V. MODEL VALIDATION
To ensure the accuracy of the developed model, the electric fields for a ten-layer (eight total shells) spheroidal structure (parameters and geometry of which were arbitrarily chosen for the purposes of validation), energized with a step voltage corresponding to an external field magnitude of 20 kV/cm, were calculated using (20) and compared to the same solution computed using the FEM solver QuickField Professional v6.4 [35].The parameters used and resulting comparison have been included as Appendix.The time-dependent solution from the many-shells model was found to agree very well with FEM, with any discrepancies likely arising from the quantization and truncation errors from the QuickField model, or minimal errors from computing the roots of (21).
To compare the developed approach with previous models established in the literature, the reader is referred to the spherical model of Kotnik and Miklavčič [36] and to the spheroidal model of Nath et al. [31], who also validated their model by means of comparison with [36].Fig. 2 compares the TMP estimations from these works under a trapezoidal voltage pulse (10 ns duration, 1 ns rise and fall time) on a spherical cell, to those estimated by the present many-shells model with n = 5.The exact conditions and parameters used are listed in [31].In the present work and in the work by Nath et al. [31], prolate spheroidal coordinates were used, and in following [31], an aspect ratio value of K = 1.01 was assumed to represent the spherical limit.Excellent agreement with both works was found.It is further remarked that the approach in [31] is identical to the present derivation, such that the equations of [31] can be recovered from the manyshells model.The exception is that this work extends the approach with general solutions, provides analytical limits to n, and provides a solution strategy that can extend beyond this limit.
As further validation of (20), one may consider the field enhancement factor inside and outside of a single-shelled (two layered) particle of arbitrary eccentricity, e. Known solutions exist for the maximum field enhancement induced by conducting inclusions where ε 1 , σ 1 ≫ ε 2 , σ 2 , as commonly applied to the analysis of metallic particles suspended in an insulating medium.By applying this assumption to (20) and normalizing the resulting field magnitude by the external field, the field enhancement factor found at the surface of the particle at ν = 0 was found to be Since K = tanh µ, the maximum level of field enhancement outside a conducting spheroid under a uniform field is a function only of its eccentricity.Fig. 3 plots ( 22) as a function of the eccentricity, which has also been verified to be identical to an expression derived by Lekner [37], and was further confirmed to tend toward a value of f = 3 when e = 0, which is an elementary solution for the maximum enhancement on the surface of a conducting sphere.Similar analysis can be conducted for a second limiting case: a single-shelled particle with ε 1 , σ 1 ≪ ε 2 , σ 2 , often used to represent a gas void embedded within a more conductive bulk in the field of electrical insulation.Equation (23) gives the derived expression for the field enhancement internal to the void, in this case, and has also been plotted on Fig. 3 Once again, the enhancement is purely a function of the eccentricity, and as e tends toward zero, ( 23) approaches the well-established limit of 3/2 for spherical geometries [38].

VI. TRANSMEMBRANE POTENTIALS
The developed many-shells model of a biological cell has been used to analyze a complex, four-shelled (six-layer) cell structure, as shown in Fig. 4.This geometry has particular significance to biomass processing technologies, which may deal with cells that incorporate a cell wall (e.g., algae or yeast).The parameters chosen in this work reflect typical values for the commonly studied yeast, Saccharomyces cerevisiae.In addition, since conventional closed-form solutions for geometries of greater than five layers do not exist (as discussed in Section III-B), this example, thus, serves to demonstrate how the novel many-shells approach can be applied despite this limit.A systematic study was conducted to investigate the effects of the pulse rise time,  pulse duration, and extracellular conductivity on the developed TMPs.The waveform of choice follows that of a doubleexponential (DE) impulse, analytically represented in the time-domain as where Ē0 is the peak electric field magnitude in V/m, α and β are wave-shaping parameters with units s -1 , and A 0 is a dimensionless scaling constant, unique to each impulse.The selection of a DE waveform was informed by the output of capacitive-based impulse generators, which have been used in PEF treatment, which, in practice, produces the characteristic Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I GEOMETRIC PARAMETERS FOR THE SIX-LAYERED SHELL
MODEL OF FIG. 4, BASED ON DATA FROM [28], [44] DE waveform with defined rising and falling edges.Unless otherwise stated, the geometric parameters, which were used to model the cell, are tabulated in Table I.Since the external field was assumed to be uniform and directed in the positive Cartesian x-direction, the maximum TMPs exist at ν = 0 following: It is widely believed that a TMP in the range of ∼0.5-1.5 V is necessary to induce electroporation, while TMPs in excess of this range almost guarantee cell necrosis.Once electroporation takes place, the formation of membrane pores leads to a significant increase in the membrane conductivity, local to the region where pores were formed, and causes the TMP to collapse.As the present model does not consider time-varying membrane conductivity, calculated TMPs the electroporation threshold are unlikely to be representative.However, the observed characteristics of the TMPs beyond ∼1 V are nevertheless discussed, for several reasons.Firstly, there exists some variation to the electroporation threshold, since the parameter is heavily dependent on cell type, geometry, pulse type, and other conditions.Developed threshold TMPs upward of 2 V have been previously suggested in some cases [16], [39], [40].Additionally, in other work [41], [42], [43], authors have elucidated details regarding the finite delay time for pore formation, which may induce a further overshoot of the TMP.Considering this, the results presented here do not limit the estimated TMPs to a particular voltage.However, the moment when the TMP reaches 1 V has been demarcated within all figures and discussed, as it represents a useful point of reference.It is further remarked that the present work has been restricted to µs-PEF, and has solely focused on waveforms which lie within a range typically used for microbial inactivation and bio-extraction purposes.The use of this model in the ns-PEF and ps-PEF regime has been left as a subject for future investigation.Sections VI-A-VI-C also assume a spherical cell, as the focus here is not on the analysis of eccentricity.Nevertheless, Section VI-D briefly discusses the cell eccentricity as a separate matter.The shell parameters were chosen to represent the commonly studied yeast, Saccharomyces cerevisiae, for which there is a wealth of existing and available literature, and have often been studied using spherical or spheroidal dielectric models.These parameters are tabulated in Table II, and unless otherwise stated, have been used.

A. Effect of Rise time
To investigate the effects of the rise time of the applied field, the wave-shaping parameters for five waveforms following the form of (24) were computed, using a swarm-like   5. TMPs developed across the plasma and nuclear membranes of the six-layer spherical cell of Fig. 4 at various rise times given by Table III, with the FWHM fixed at 10 µs.
parameter optimization algorithm.Each signal is characterized by its 10%-90% rise time (t f ) and its full-width at half-maximum (FWHM).In this section, the FWHM was held constant at 10 µs, while the front-time was varied from 100 to 500 ns (typical for µsPEF applications, e.g., [46], [47], [48]), in steps of 100 ns.The computed wave-shaping parameters are shown in Table III.For each wave shape, the maximum TMPs across the plasma and nuclear membranes were estimated using the many-shells model, at an angle of ν = 0 and over 10 ms.The obtained results are shown in Fig. 5, where red circles indicate the moment when the TMP value exceeds a threshold of 1 V.
In general, pulses in the present microsecond range produced significant TMPs across the plasma membrane, which exceeded the threshold value of 1 V in all cases.Although the nuclear membrane also developed a TMP, the maximum TMP over the course of the input pulse did not exceed 1 V, in agreement with the known characteristics of µsPEF, which predominantly targets the outer membrane.The rise time, however, had a clear effect on the time necessary for the plasma membrane TMPs to reach 1 V, which were prolonged with a longer signal front-time, as evidenced from Fig. 5.An increase from t f = 100 to t f = 500 ns prolonged the time to reach a TMP of 1 V from around 100 to 300 ns.Longer front-time also saw the decrease in the maximum TMP magnitude across the nuclear membrane, the peak of which was also delayed in time.Beyond 1 V, plasma and nuclear TMPs appear to converge, regardless of fronttime.This may suggest that, after the formation of pores and the stabilization of the pore density, the subsequent bioelectric effects on the electroporated membrane may not differ significantly with different rise times under the range of conditions studied here.However, the authors emphasize that this suggestion is made based on the results here where explicit pore formation (e.g., using Smoluchowski's equation) was not considered.Validation of these results with explicit pore formation taken into account is an area for future work.

B. Effect of Pulse Duration
Following the methodology described in Section VI-A, this section presents the results of a complementary study on the effects of the pulse duration, in this case, defined as the FWHM.Holding the front-time constant at 200 ns, parameters for five waveforms of 1, 5, 10, 100, and 500 µs FWHM were computed, as shown in Table IV.The resulting TMPs are shown in Fig. 6.The variation of the FWHM over this range found effects neither on the time necessary for the plasma membrane to reach the threshold TMP nor  [28], [51] on the peak TMP magnitude across the nuclear membrane.However, since the FWHM effectively controls the dynamics of the wave-tail, significant effects were observed beyond the ∼1 V threshold.It was estimated that a prolonged pulse duration resulted in higher maximum plasma membrane TMP, accompanied by a significant delay in the time for this maximum TMP to be reached, as shown in Fig. 6.If the electroporation threshold TMP is within the higher (>1.5 V) range, this may result in waveforms of shorter pulse duration being unable to induce electroporation, or produce only reversible effects.The time that the plasma membrane TMP was above the voltage threshold was also prolonged with longer pulse duration, which may suggest the intensification of any bioelectric phenomena at the membranes, during and after pore formation.This supports the findings of experimental studies such as [20], [49], which reported a positive correlation between effective pulse duration and PEF inactivation performance.

C. Effect of Extracellular Medium Conductivity
In this section, the effects of the electrical conductivity of the external medium on the TMP dynamics were investigated.PEF systems used in bio-extraction and food processing technologies will encounter a variety of target media, which may have vastly different characteristics.In general, changing the properties of the extracellular media will also cause changes to the properties of the other cell layers [27], [28], which must be incorporated when modeling TMPs to better reflect practical values.In this work, low (σ e = 1 mS/m), medium (σ e = 50 mS/m), and high (σ e = 0.1 S/m) conductivity external media were considered, following the strategy in [28].These may represent standard laboratory solutions such as deionized water [28] or on the higher end, certain liquid foodstuffs that may be subject to PEF treatment in the food processing industry [50].Alongside changing the parameters of the extracellular medium, the induced changes on the cell wall and membranes were also incorporated according to [28].It is remarked that these parameters were originally measured and reported in [51], where both living and dead (heat-treated) yeast cells were investigated.The parameters used reflect the case of living cells, and are tabulated in Table V.The waveform was fixed with a 200-ns front-time and 10 µs FWHM.It should be noted that, for simplicity, it was assumed that the applied field can be fully established regardless of the medium conductivity.Fig. 7. TMPs developed across the plasma and nuclear membranes of the six-layer spherical cell of Fig. 4 with different extracellular conductivity as given by Table V, using a 0.2/10 µs pulse.
As before, plasma and nuclear TMPs were estimated in each case and plotted together in Fig. 7. Higher extracellular conductivity resulted in more rapid establishment of threshold TMPs than for lower conductivity values, resulting from the rapid (and more intense) redistribution of the electric field to the far less conductive membranes.A corresponding increase in the nuclear membrane TMPs was also observed, such that, for the parameters used in this study, the nuclear membrane also came close to the 1-V threshold in the highest conductivity case.Beyond 1 V, lower conductivity resulted in a prolonged time for the TMP to peak and a slight reduction of the maximum TMP.For applications targeting cell lysis, a higher conductivity medium may, therefore, prove beneficial to PEF that aims to induce IRE.

D. Note on Cell Eccentricity
The model developed in this work assumed the arbitrary eccentricity of the modeled cell, by making use of prolate spheroidal coordinates.This provides a somewhat better representation of cells that deviate more substantially from the spherical case, or may be used to approximate longer, cylinder-or rod-like geometries characteristic of some bacterial cells.However, the limitations imposed from the use of the confocal spheroid approach were already discussed in Section I.When applied to model geometries, which may include many layers or have certain layers which are significantly thicker than others, one must be cautious when interpreting the estimated TMPs, as the geometry may begin to deviate significantly from reality.Since the effects of the cell eccentricity are highly dependent on the specific cell characteristics, a detailed analysis will not be presented in this work.For the Saccharomyces cerevisiae cell demonstrated in this study, microscopic images indicate that they are often ovular, but not significantly [44], [52].When the eccentricity was changed to have a minor to major axis ratio of K = 0.75, insignificant changes in the TMP responses were observed, showing that the deviation from the spherical form was not sufficient to have a large impact on this case.VII.CONCLUSION In this work, the generalization of the commonly used multishell model for estimating transmembrane voltages in PEF applications has been presented.A many-shells dielectric model for biological cells has been developed using prolate spheroidal coordinates, incorporating an arbitrary number of poorly conducting "shells," and subjected to a uniform time-dependent electric field of DE form.Mathematical analysis indicates the bounds for which closed-form solutions are able to be obtained, and a full methodology that can be applied to cases beyond this range has been described.The model, therefore, allows TMP estimates and relaxation characteristics to be computed, that are free from the limitations of mesh-based approaches such as the FEM method, for layered cell geometries of significant complexity.
For an isolated yeast cell model under the application of a 30-kV/cm (peak) DE impulse in the microsecond regime, the many-shells model indicated that the signal rise time likely has a weak effect on the developed TMPs across the plasma and nuclear membranes, despite prolonging the time required to reach threshold TMPs on the lower range (∼0.5-1V).Pulse duration, on the other hand, had little effect on the time to reach a TMP of ∼1 V or lower, but may have more pronounced impact on cells with threshold TMPs in the higher (≳1.5 V) range.It is also thought that behavior beyond the threshold voltage may indicate that the pulse duration has major effects on the subsequent efficacy of PEF bioelectric effects on the porated membrane.Under three different values of extracellular medium conductivity (and also considering the induced changes to the other layer properties), results suggest that high-conductivity media may be beneficial for PEF processes targeting lysis and inactivation, though the balance with thermal effects must also be considered.
The present model, however, does not consider dynamic pore formation and the subsequent changes induced in the membrane conductivity.It would be of high interest to couple such an approach with pore formation models, such as the Smoluchowski equation, to analyze the developed pore densities under various configurations.It would also be of significant interest to further investigate PEFs in the nanosecond and picosecond regimes, believed to be effective for the treatment of cancers and neurodegenerative diseases.The ability for the model to incorporate n layers also opens up the possibility for its application as a first-order approximation of radially nonuniform layer parameters, as a further step toward higher fidelity models.

APPENDIX TEN-LAYER MODEL VALIDATION WITH FEM
The parameters used for the comparative validation study of a ten-layer inclusion subjected to a step function are shown in Table VI, and the comparison results of the developed electric fields are shown in Figs. 8 and 9 for various timesteps.The utilized step function incorporates an instantaneous transition from zero to +20 kV/cm.The numerical simulation was performed using QuickField Professional version 6.4 [53].

Fig. 1 .
Fig.1.Prolate spheroidal coordinate system used to develop the many-shells dielectric model.

Fig. 3 .
Fig. 3. Comparison of the two limiting cases of field enhancement factor between the present model and known limits.

Fig. 4 .
Fig. 4. Six-layer cell domain considered within this work, as a model of the yeast Saccharomyces cerevisiae.Note that the cell is assumed spherical (e = 0), details on the effects of e are discussed further in Section VI-D.Pictorial representation only; relative layer thicknesses are not to scale.

Fig. 8 .
Fig.8.Electric field magnitude of (top row) QuickField FEM computed solution and (bottom row) analytical many-shells to a ten-layer spheroidal inclusion subject to +20 kV/cm step function, with parameters as defined in TableVI.

Fig. 9 .
Fig. 9. Comparison of QuickField and analytical field solution at t = 0.01 s for ten-layer spheroidal inclusion, along the line y = 0, x ≥ 0.
Graduate Student Member, IEEE, Igor Timoshkin , Senior Member, IEEE, Scott MacGregor , Senior Member, IEEE, Mark Wilson , Member, IEEE, Martin Given , Senior Member, IEEE, and Michelle Maclean

TABLE IV WAVE
-SHAPING PARAMETERS FOR THE INVESTIGATION OF THE PULSE DURATION IN SECTION VI-B Fig. 6.TMPs developed across the plasma and nuclear membranes of the six-layer spherical cell of Fig. 4 at various FWHM given by TableIV, with the rise time fixed at 0.2 µs.

TABLE V ELECTRICAL
PARAMETERS FOR LOW, MEDIUM, AND HIGH EXTRACELLULAR CONDUCTIVITY RANGES USED IN SECTION VI-C, GATHERED FROM

TABLE VI PARAMETERS
USED FOR TEN-LAYER COMPARISON STUDY WITH FEM SOLVER