A Multiscale Computational Model of Skeletal Muscle Electroporation Validated Using In Situ Porcine Experiments

Objective: The goal of our study was to determine the importance of electric field orientation in an anisotropic muscle tissue for the extent of irreversible electroporation damage by means of an experimentally validated mathematical model. Methods: Electrical pulses were delivered to porcine skeletal muscle in vivo by inserting needle electrodes so that the electric field was applied in direction either parallel or perpendicular to the direction of the muscle fibres. Triphenyl tetrazolium chloride staining was used to determine the shape of the lesions. Next, we used a single cell model to determine the cell-level conductivity during electroporation, and then generalised the calculated conductivity changes to the bulk tissue. Finally, we compared the experimental lesions with the calculated field strength distributions using the Sørensen-Dice similarity coefficient to find the contours of the electric field strength threshold beyond which irreversible damage is thought to occur. Results: Lesions in the parallel group were consistently smaller and narrower than lesions in the perpendicular group. The determined irreversible threshold of electroporation for the selected pulse protocol was 193.4 V/cm with a standard deviation of 42.1 V/cm, and was not dependent on field orientation. Conclusion: Muscle anisotropy is of significant importance when considering electric field distribution in electroporation applications. Significance: The paper presents an important advancement in building up from the current understanding of single cell electroporation to an in silico multiscale model of bulk muscle tissue. The model accounts for anisotropic electrical conductivity and has been validated through experiments in vivo.


I. INTRODUCTION
E LECTROPORATION, electropermeabilisation, pulsed electric field (PEF) treatment, and pulsed field ablation (PFA) are terms used to describe the same phenomenon in different fields of research and applications. When a biological cell is exposed to an external electric field of sufficient amplitude, its membrane becomes transiently permeable to ions and molecules for which it is otherwise poorly permeable or not permeable at all (reversible electroporation). The exposure may lead to cell death, in which case the phenomenon is called irreversible electroporation [1], [2], [3]. Electroporation is used in various fields of biomedicine, biotechnology, as well as in food processing [4], [5], [6], [7], [8], [9], [10].
One of the important applications of electroporation in biomedicine is gene electrotransfer (GET) [11]. In GET, the phenomenon of electroporation is used to introduce genetic material (usually plasmid DNA) into cells to achieve a desired therapeutic effect. Skeletal muscle is the most commonly used target tissue for gene therapy and DNA vaccination due to its ability to express genes and secrete proteins into the bloodstream. Since skeletal muscle cells do not divide, gene expression is sustained for several months after GET [12], [13], [14]. Another increasingly important application of electroporation in biomedicine is treatment of cardiac arrhythmias, particularly of atrial fibrillation, by ablation of pulmonary vein tissue using irreversible electroporation [10], [15], [16], [17], [18]. The specific structure of skeletal and cardiac muscles dictates an anisotropic electrical conductivity, leading to a fibre-orientation dependent electric field distribution upon pulse application [19], [20], [21]. The electrical conductivity of skeletal muscle in the direction of its fibres was found to be higher than the electrical conductivity in the direction perpendicular to the fibres [22], [23], [24], [25], [26], [27]. Anisotropic electrical conductivity has also been observed in cardiac muscle [28]. It has been reported previously that the (reversible) threshold of electroporation depends on the direction of the applied electric field with respect to the direction of the muscle fibres -lower pulse amplitudes were required for parallel orientation: 80 V/cm versus perpendicular orientation: 200 V/cm [29]. Interestingly, a tissue that is not intrinsically anisotropic, such as the liver tissue, can become anisotropic after electroporation [30]. Anisotropy can also have a significant influence in PEF applications in the food industry, such as meat and fish processing [31]. Namely, many food plant tissues exhibit anisotropic properties, e.g. asparagus [32].
It is also known that electroporation at the cell level in vitro depends on the cell geometry and the orientation of the cell in the electric field [33], [34], [35], [36], whereas the preferential orientation for electroporation depends also on pulse duration and amplitude [37]. Electroporation leads to an increase in membrane conductivity, which results in bulk tissue conductivity increase [38], [39], [40], [41]. The orientation of muscle fibres is not constant throughout the muscle, which is most pronounced within the cardiac muscle, where the orientation of cardiomyocytes varies not only spatially through the heart, but also from the epicardial to the endocardial side [42]. Therefore, muscle is considered a heterogeneous tissue, which can considerably affect the success of muscle gene transfection and cardiac ablation. The increase in bulk tissue conductivity due to or during electroporation has been described previously and was used to control pulse delivery, thus limiting the damage induced due to electroporation [43]. It has been described by various functional dependencies of conductivity as a function of local electric field [44], [45], but until recently was not explicitly linked to changes in membrane conductivity. The conductivity of the cell membrane and consequently the conductivity of the cell is increased in the direction of the electric field of the delivered electroporative pulses, i.e., it can be said that electroporation induces or increases the anisotropy of the tissue [30]. Previously, a skin structure with different cell shapes and packing densities was presented and numerically compared with skin impedance measurements [46], [47], [48], which also provided the basis to connect membrane electroporation to the increase in tissue conductivity due to electroporation [49].
Hitherto, the increase in anisotropy of conductivity due to electroporation has not been well described. Namely, electroporation in one direction contributes to the change of conductivity in all directions, not only in the direction of the applied electric field / electroporation. We have therefore constructed skeletal muscle tissue from cells, similar to how it was demonstrated for skin keratinocytes in dense suspensions in the models of Huclova et al. This approach allows to connect the micro-to the macro-scale of electroporation, and to deduce about membrane changes backwards, i.e. from measurements made at the tissue scale, we can in principle infer the degree of membrane permeabilisation required to achieve the observed changes in the bulk tissue conductivity [50].
The anatomic features and physiological functioning of skeletal muscle has been well described in the literature [51], [52], [53]. When targeting skeletal muscle for electroporation, we need to consider the given architecture of the muscle cells and their relative lengths. Muscle cells are multinucleated and can have lengths ranging from millimetres to over 10 cm [54]. Each muscle fibre itself is surrounded by a connective tissue known as endomysium. These fibres, or myofibres, then form clusters The skeletal muscle is organised into various substructures, with connective tissue layers separating those structures. The major muscle is surrounded by a connective tissue layer called epimysium. The muscle contains a large number of fascicles (shown as a cross-section) which is surrounded by a tissue layer called perimysium. Each fascicle contains a number of individual muscle cells surrounded by a tissue layer called endomysium. Each muscle cell is multinucleated, and contains myofibrills, composed of the contractile muscle proteins actin and myosin. known as fascicles, which are again wrapped in additional connective tissue, the perimysium. Finally, the fascicles are then enveloped in another layer of connective tissue known as the epimysium, and the skeletal muscle such as the vastus lateralis is surrounded by a final layer of connective tissue, the perimysium (Fig. 1). These structures were described and imaged with high-resolution MRI with histology images [55]. Each connective tissue layer presents a complex system with differing densities and compositions [56]. To compound to the complexity of the skeletal muscle, there are nerve fibres and blood vessels intertwined throughout these connective tissue layers.
Our study consisted of two parts, an experimental and a numerical part, initially separated but then intertwined (see Fig. S1 in the Supplementary Materials). In the experimental part of the study, we performed in vivo experiments in the skeletal muscle tissue of pigs. We delivered electrical pulses to the tissue, inserting needle electrodes in two different orientations with respect to the muscle fibres. In the first group, the direction of the applied electric field was parallel to the direction of the muscle fibres (later referred to as parallel orientation), and in the second group, the direction of the applied electric field was perpendicular to the direction of the muscle fibres (perpendicular orientation). We then used triphenyl tetrazolium chloride (TTC) staining to determine the shape of the lesion for each application of pulses. TTC stains living cells in bright red while dead tissue appears white/pale. The white compound is enzymatically reduced to red TPF (1,3,5-triphenylformazan) in living tissues due to the activity of various dehydrogenases. This stain enables determination of the irreversibly ablated region, i.e. the lesion [57].
In the numerical part of the study, we built a multiscale numerical model of skeletal muscle. We first used a single-cell model to determine cell-level conductivity during electroporation. Next, we implemented these results into a bulk model of the tissue.
We then combined the experimental and numerical parts of the study in such a way that we used our numerical model to simulate the in vivo experiments and determined irreversible threshold of electroporation for each experiment by comparing the numerical, i.e. model results to the shape of the experimental lesion.
The purpose and main objective of the study is two-faceted. Firstly, we aimed to show how bulk tissue properties, specifically anisotropy, can be deduced from a single-cell model, i.e., to demonstrate an approach of building muscle tissue from its basic constituents, the muscle fibres. This fibre to tissue scaleup consists specifically of observing electrical properties of bulk muscle based on electrical characteristics of single cells both before treatment, and then as altered due to electroporation. Secondly, we aimed at both validating the model by comparing TTC-stained lesions in shape and size to numerically calculated electric field distributions, as well as quantifying the irreversible electroporation threshold of muscle tissue in vivo for specific pulse parameters. We made no a priori assumptions regarding the threshold in terms of its dependence on field-to-fibre orientation, i.e., we built no assumptions into the model that the threshold is dependent on field-to-fibre orientation. From this point of view the study is also an attempt at determining whether anisotropic electrical properties of muscle tissue will lead to an orientation-dependent irreversible field strength threshold on the level of a single cell, or, are the orientation-dependent lesions shapes and sizes a result of electric field distribution in an electrically anisotropic medium.

A. In Vivo Experiments
In situ preclinical porcine models have been shown to be a reliable and widely used for translational ablation studies, including those for assessing the effect of electroporation therapies [58]. Further, the porcine vastus lateralis muscle has been shown to elicit reproducible lesion sizes and depths and thus was considered as a viable model for preclinical ablation studies [59].
Adult male Yorkshire pigs (75-85 kg) were used. All animals received humane care in compliance with the 'Principles of Laboratory Animal Care', formulated by the National Society for Medical Research, and The Guide for the Care of Laboratory Animals, published by the National Institutes of Health. This research protocol was approved by the University of Minnesota's Institutional Animal Care and Use Committee (IACUC protocol number: #2006-38201 A; last approval date: 28 January 2022).
The pigs were sedated with methohexital (20-50 mg/kg as needed), then intubated and mechanically ventilated. Anaesthesia was maintained using isoflurane (>2 %) and a 2:1 airto-oxygen mixture. The heart rate and blood pressure of each animal were monitored continuously (SpaceLabs, WA, USA) and properly maintained. Their core temperatures (rectal temperatures, YSI thermocouples, City, State) were monitored and maintained between 38 + 0.5 • C, external convective warming (BairHugger, 3 M, USA) was utilised as needed. Next, the vastus lateralis muscles of each pig were carefully exposed. To do so, the skin was cut using electric cautery adjacent to the muscle so not to impair function and the remaining connective tissue and fascia was blunt dissected to leave the muscle fibres exposed (see Fig. S2a in the Supplementary Materials). Thus there were no tissues adjacent to the treated muscle present that could affect the electric field distribution within the treated muscle tissue and its conductivity. The microstructure of the tissue also has no significant influence on the distribution of the electric field on a macroscopic level and can therefore be neglected, but it does influence the distribution of the electric field on a microscopic level / locally [45], [60]. As can be observed in Fig. S2a in the Supplementary Materials, the fascicles and muscle fibres of the vastus lateralis run parallel to the long axis of this muscle (origin to insertion). The muscle and fibres are on average 10 cm long and 4.5 cm wide, and the fascicles are about 4 mm in diameter.
Electroporation was performed by delivering 48 pulses with a pulse width of 100 µs in 6 trains of 8 pulses each, with a pulse repetition rate of 1000 s −1 in each train, and a pause of 2 seconds between trains. Custom built pulse generator was used. The pulses of four different amplitudes (600 V, 800 V, 1000 V, and 1200 V) were delivered through two needle electrodes with a diameter of 0.70 mm and an centre-to-centre distance of 8.0 mm. The total length of the electrodes was 8 mm, with the upper 1 mm insulated, i.e. 7 mm was the length of the non-insulated part of the electrodes, to ensure that the active part of the electrodes inserted into the tissue was the same for each experiment performed. The decision to use only 2 needle electrodes was based on the literature [50] showing that 2 needle electrodes provide an adequate gradient of electric field, have a well-defined geometry and at the same time minimise trauma to the tissue. This electrode geometry also allows comparison with other relevant studies, e.g. [61], and allows easy determination of orientation with respect to the muscle fibres. The 2 needle electrode geometry also maximises the number of possible lesions per muscle and thus minimises the number of animals needed for experiments. The electrodes were placed so that the direction of the electric field was as parallel or perpendicular as possible to the direction of the muscle fibres. Each electrode position was at least 5 cm apart and the types of ablation deliveries were randomised before each day's experiment. All voltage and current measurements were acquired with an oscilloscope (Keysight MSOS104 A, California, USA) recording at 1 GHz. High voltage differential probes (Keysight N2891 A) with an attenuation ratio of 1000:1 were fed into the oscilloscope. The probes were clamped directly onto a specialised space within the needle electrode holder. Muscle contractions were visually observed with all pulse deliveries and for the electrodes placed with the varied fibre orientations. From these observations, we could not detect any differences in contraction response between muscle stimulations with the varied electrode placements. As expected, the greater applied pulse amplitude induced visibly larger contractions: i.e., more muscle tissue was activated.
Following pulse delivery, all lesions were allowed to mature for 1.5 hours before the animal was euthanised via cardioplegia delivery. Immediately post-mortem the vastus lateralis muscles were isolated and removed and taken to a dissection table where the most superficial layer was sectioned off. This was done to allow the staining agent to penetrate the muscle cells, as

B. Single-Cell Numerical Model
Calculations of cell-level conductivity during electroporation were performed in Comsol Multiphysics v6.0 (Comsol Inc., Sweden). We used a unit-cell approach as described in [46] and applied for cell electroporation previously [49]. First, we constructed a 3D cuboidal biological cell in Matlab R2022a (MathWorks, USA) using the superformula [62] ((1)-(9) in [46]) and parameters for a cuboidal cell (Table I in [46]). The constructed cell was saved in a drawing interchange format (.dxf file) and imported to Comsol using the Import function under Geometry with the CAD module. The imported biological cell was then scaled to obtain the desired geometry (r long and r short , Table S-I in the Supplementary Materials, and Fig. 2) of a single muscle cell. A unit cell around the muscle cell was constructed as a block. Its geometry was adapted so that the volume fraction of the biological cell was 76 %. This value is lower than that of the actual muscle tissue [63], but it is the highest value that can be achieved with the geometry used without impractically increasing the computational cost. The effect of the volume fraction value used in the model on the calculated cell-level conductivity results is further addressed in the Results and discussion section. Because of the symmetry of the model, only 1/8 of the model was calculated to decrease the computational cost.
For the calculation of electroporation, we coupled the Weak Form Boundary PDE (partial differential equation) and the Electric Current interface in the same way as described in [49] with the parameters [47], [64], [65] summarised in Table S-I in the Supplementary Materials. The Weak Form Boundary PDE was used to calculate the change in pore density, while the Electric Curents were used to calculate the Laplace equation. Eight square 100 µs long pulses with a rise-and fall-time of 1 µs were applied either in parallel or perpendicular direction with respect to the direction of the long axis of the muscle cell to two opposite boundaries with a pulse repetition rate of 1000 s −1 . The other boundaries were set as electrically insulated. The muscle cell membrane was modelled as a contact impedance boundary condition and its conductivity (σ m ) increased due to the change in pore density (N ): (parameters described in Table S-I in the Supplementary Materials). Pore formation was calculated with the asymptotic pore equation with fixed pore radii, which then increased the membrane conductivity, as described in [49], [65]. The transmembrane voltage (U m ) was calculated as the difference between the potential on the inside and the outside of the cell membrane. For calculation of electrical conductivity of the unit cell, two additional Laplace equations were added to the model via the Electric Current interfaces, one for calculation of the electrical conductivity in the direction of the long axis of the muscle cell (later referred to as parallel conductivity), and the other for calculation of the electrical conductivity in the direction perpendicular to the long axis of the muscle cell (perpendicular conductivity). In both, a fixed voltage of 1 V (V test ) was applied throughout the whole course of simulation. The voltage of 1 V used had no influence on the results of the electroporation part of the model, as these additional interfaces were used exclusively for the calculation of the two conductivities and were not coupled to the primary Electric Current interface. First, a stationary study was run to initialise the parameters of the two physical interfaces, used for calculations of the conductivities. Then, a time-domain study was performed for all included physical interfaces to model the pore formation and calculate the corresponding conductivity change. The simulation time was 30 seconds with denser time points during the pulse application, which is sufficient to observe the closing of the pores, as well as the faster dynamics during the pulse application.
Unit cell conductivity was obtained by using the currents caused by V test and calculated within the two corresponding Electric Current interfaces. These two currents were obtained by calculating a boundary integral of a normal current density across a cross-section (y-z plane for the parallel current and x-z plane for the perpendicular current, see Fig. S3 in the Supplementary Materials for the coordinate system orientation). Taking into account the geometry of the cross-section, i.e. of the unit cell (d x , d y , d z ) as well as the applied voltage (V test ), the conductivity of the muscle cell during and after pulse application was determined in S/m and used in further calculations. The results of these calculations are shown in Fig. 3 for the two cases where the direction of the electric field is either parallel or perpendicular with respect to the direction of the long axis of the muscle cell, and for both cases, the conductivity both in direction parallel or perpendicular with respect to the direction of the long axis of the muscle cell is given.
To summarise: In the numerical model of a single cell that serves as a means of bulk tissue properties determination, the pulses are applied to the unit cell directly by prescribing the appropriate initial and boundary conditions to the surfaces of the unit cell, and the geometry of the electrodes is not involved. The results of this model, specifically the electric conductivity distribution in and around the cell, is then fed into a bulk model of tissue where the electrodes are modelled (as needle electrodes) and the bulk tissue properties used to calculate the electric field distribution in bulk tissue.

C. Bulk Tissue Numerical Model
The bulk tissue model of the skeletal muscle, used for the calculations of the distribution of the electric field during the delivery of the electroporation pulses was built in Comsol Multiphysics v6.0 (Comsol Inc., Sweden). The geometric representation of the skeletal muscle tissue in the model was a cuboid with dimensions of 40 mm × 40 mm × 10 mm, with the muscle fibres assumed to be aligned along the x-axis in the model. Two needle electrodes with the same dimensions as used in the experiments (i.e., diameter of 0.70 mm, insertion length in the tissue of 7 mm, and the centre-to-centre distance between the two electrodes of 8.0 mm) were added to the model. In each simulation, the direction of the applied electric field with respect to the muscle fibres was the same as it was in each experiment. In the model we used the Electric Current physics interface with a time domain study.
Based on the results of the single-cell model, shown in Fig. 3, we defined interpolation functions of electrical conductivity for use in our bulk tissue model. We used sequential linear interpolation. We defined four interpolation functions: two for the case when the electric field is applied in direction parallel with respect to the direction of the muscle fibres (parallel orientation), and two for the case when the electric field is applied in direction perpendicular with respect to the direction of the muscle fibres (perpendicular orientation). For each of these two cases, we defined one function for the conductivity in direction parallel with respect to the direction of the muscle fibres (parallel conductivity), and another function for the conductivity in direction perpendicular with respect to the direction of the muscle fibres (perpendicular conductivity) -all four functions are shown in Fig. 4 for a single 100 µs pulse. In this way, we have included electrical conductivity in the model as a function of time and also as a function of the distribution of the magnitude of the electric field. However, since the distribution of the electric field in the bulk tissue model is also a function of space, this means the conductivity in our bulk tissue model is defined as a function of time and also as a function of space, separately for the value of the conductivity in direction parallel with respect to the direction of the muscle fibres, and for the value of the conductivity in direction perpendicular with respect to the direction of the muscle fibres.
The four conductivity functions defined above can be used for the case when the electric field is applied in direction either exactly parallel or exactly perpendicular with respect to the direction of the muscle fibres. However, for an arbritary orientation of the applied electric field ϕ with respect to the direction of the muscle fibres (see Fig. 4(e) for a schematic representation of ϕ with respect to the muscle fibres), we need to define additional functions by combining the above defined four. To define the electrical conductivity in direction parallel with respect to the Fig. 4. Conductivity as a function of time and the magnitude of the electric field, in different directions with respect to the muscle fibres, used in the bulk-tissue model. The functions were defined by using the results of the unit-cell model (Fig. 3) and sequential linear interpolation. For clarity, the functions are only shown for the first pulse. In (a) and (b) the direction of the applied electric field was parallel with respect to the direction of the muscle fibres (parallel orientation), and in (b) and (d) perpendicular with respect to the direction of the muscle fibres (perpendicular orientation). In (a) and (c) the functions for the conductivity in direction parallel with respect to the direction of the muscle fibres are shown (parallel conductivity), and in (b) and (d) the functions for the conductivity in direction perpendicular with respect to the direction of the muscle fibres (perpendicular conductivity). A schematic representation of the orientation of the electrodes with respect to the direction of the muscle fibres is given in (e). direction of the muscle fibres, we used (2): and to define the electrical conductivity in direction perpendicular with respect to the direction of the muscle fibres, we used (3): wherein the subscript index indicates the direction of the applied electric field with respect to the direction of the muscle fibres, the superscript index indicates the direction of conductivity with respect to the direction of the muscle fibres, and the superscript ϕ indicates the case when the electric field is applied in an arbitrary orientation ϕ [ Fig. 4(e)] with respect to the direction of the muscle fibres. We used the bulk tissue model described above to calculate the spatio-temporal distribution of the electric field in the tissue for each experiment. To reduce the computational cost, we used only a single train of eight 100 µs long pulses in our simulations, as the results in Fig. 3 show that the calculated conductivity reaches a plateau during the first pulse train for all values of the electric field.

D. Irreversible Electroporation Threshold Determination
On each image of the lesion obtained from the in vivo experiments, we marked four pixels: two representing the location of the two needle electrodes and the other two defining the direction of the muscle fibres. We used these markings to determine the actual orientation of the electrodes with respect to the muscle fibres ϕ, and to manipulate the images so that they could later be compared with each other and also be compared with the calculated electric field distributions. We magnified and rotated each image so that the electrode locations were always at the same coordinates: (−4.0 mm, 0) and (4.0 mm, 0), since the centre-to-centre distance between the two electrodes is 8.0 mm. We then used our bulk tissue model to calculate the time-dependent electric field distribution in the muscle tissue for a single train of 8 pulses with a pulse width of 100 µs and a pulse delivery rate of 1000 s −1 . We performed the same number of in silico simulations as was the number of the in vivo experiments: in each of the simulations, we used the actual orientation of the electrodes with respect to the muscle fibres ϕ as determined from the lesion images, and used the same voltage as in the experiments. For each simulation, we generated contour images of the electric field distribution at the end of the 8 th pulse in the same manner as the lesion images (i.e., with the electrode positions at the same coordinates). We generated a seperate image for each value of the electric field. We chose range of values from 100 V/cm to 350 V/cm to ensure that the irreversible threshold value was within expected range and used a step size of 1 V/cm. To determine the value of the irreversible threshold of electroporation for each of the experiments performed, we calculated the Sørensen-Dice similarity coefficient for each of the generated contour images using the following equation [66]: where X is the actual lesion image and Y i is the i th calculated contour image. We chose the value of the electric field from the contour image that had the highest Sørensen-Dice similarity coefficient with the actual lesion image from one of the experiments as the value of the irreversible threshold of electroporation for that particular experiment. For calculations of Sørensen-Dice similarity coefficients we used Matlab R2022a (MathWorks, USA). The flowchart of the methodology of the present study is shown in Fig. S1 in the Supplementary Materials.

III. RESULTS AND DISCUSSION
Electric pulses delivered in vivo in parallel or perpendicular orientation with respect to the direction of the muscle fibres yielded markedly different lesions both in size and shape (Fig. 5, Table I). The lesions from the parallel orientation group were consistently smaller than the lesions from the perpendicular orientation group at each of the voltages used [ Fig. 5(c)]. This also applies to the lesion width [dimension n in Fig. 5(e)], but not to the lesion length [dimension m in Fig. 5(d)]. We were able to distinguish between two very characteristic and morphologically disparate lesion types when the field is delivered in direction parallel to the muscle fibres [ Figs. 5(a), 6(a-d)] as opposed to   [61] and is the typical, expected shape of a two-needle electrode delivery reversible/irreversible electroporation area in electrically homogeneous (i.e. isotropic) material [50] [ Fig. 6(a-d)]. Quite differently to the parallel, in the perpendicular direction of electric field with respect to the direction of muscle fibres, the irreversibly damaged area extends far out and away from the line connecting the electrodes, and this extension seems to be most prominent along the plane bisecting the line connecting the two electrodes [ Fig. 6(e-h)]. This may seem counter-intuitive, since the electric field right in the middle between the two electrodes is normally the lowest (for parallel orientation or in homogeneous, isotropic tissue [50]), however, the strong anisotropy in tissue electrical conductivity seems to result in the electric field strongly extending into this central plane between the electrodes.
As outlined in the Materials and Methods section, the multiscale approach allows for calculating both the initial as well as electric field strength-dependent muscle tissue conductivity in both the parallel to fibres and perpendicular to fibres orientation (both of which then need to be further analysed with respect to the direction in which the electric field is being generated). This is given in brief by (2) and (3), and by interpolating the electric field-dependent conductivity evolution curves as given in Fig. 4, we can determine (simulate) the electric field strength distribution (and current density, if desired) for an arbitrary placement of electrodes in tissue and an arbitrary voltage applied to the electrodes. This is possible as long as the pulse protocol Fig. 6. Comparison of the experimentally-obtained lesions with the results of the numerical model. One sample is shown for each experimental group; parallel orientation (i.e., the electric field was applied in direction parallel with respect to the direction of the muscle fibres) in (a) -(d), and perpendicular orientation (i.e., the electric field was applied in direction perpendicular with respect to the direction of the muscle fibres) in (e) -(h); the voltage used in the experiments was 600 V in (a) and (e), 800 V in (b) and (f), 1000 V in (c) and (g), and 1200 V in (d) and (h). The value of determined irreversible threshold is given for each sample presented. The arrows mark the discrepancy in the morphology of the calculated as compared to the experimentally determined lesion shape.
that is simulated matches that for which the evolution of conductivity was calculated (see Fig. 3). The resulting calculated electric field distribution can be compared to the appropriately oriented and magnified lesion image, and using an appropriate surface matching algorithm, we can determine the field strength value for which the calculated lesion shape best matches the experimental one. We refer to this process as thresholding in continuation. In Fig. 6 and Table I, we present the results of thresholding for eight representative lesions with their optimally matched calculated electric field thresholds. We can observe that the shape of the calculated distribution closely matches experimentally obtained lesion shape.
There is however some discrepancy in the morphology of the calculated as compared to the experimentally determined lesion shape that our model cannot explain through merely modelling electric field distribution, and we attribute this discrepancy to a yet unidentified mechanism or phenomenon related to either complex muscle anatomy, or cell death, or the damaged cell staining mechanism. This discrepancy is most obvious in Fig. 6(e) and (g) (marked with red arrows), where we can observe a rather sharp change in continuity of the lesion edge in several places along the lesion's circumference.
The comparison between experimentally determined lesions and calculated electric field distributions allows us to determine the effective irreversible electroporation field strength threshold for this particular muscle and pulse protocol. Results are given in Fig. 7 and Table I. According to the theory of electroporation [3], this irreversible threshold should not depend on the voltage applied to the electrodes. Indeed, Fig. 7(c) clearly demonstrates that the determined field strength threshold is independent of the applied voltage on the electrodes be it either 600, 800, 1000 or 1200 V, with an average value of 193.4 V/cm and a standard deviation of 42.1 V/cm [ Fig. 7(a)]. This scatter may seem large, however, in view of the enormous number of variables difficult to control for in the in vivo experiments, this is perhaps surprisingly small (the coefficient of variation is 21.8 %). Moreover, there is no statistically significant difference in the average threshold if comparing parallel and perpendicular orientation of applied electric field with respect to the direction of the muscle fibres [ Fig. 7(b)]. This is in line with theoretical expectation, as we do not expect the orientation to influence the sensitivity of the cell to electroporation-induced damage and its capability of recovery following electroporation. The different lesion shape (and size!) that can be observed between the parallel and perpendicular orientation can thus be explained largely by the difference in distribution of the electric field (and current density) during pulse delivery as resulting from anisotropy in electrical conductivity, without resort to introduction of additional complexity such as orientation-dependent irreversible electroporation threshold. The determined value of the irreversible electroporation threshold of 193.4 V/cm is much lower than the one reported previously (i.e. 450 V/cm in [29]). This is not surprising since we used a much higher number of pulses (48 compared to 8), knowing that membrane electroporation depends on the number of pulses applied [67]. In addition, the experimental conditions differed considerably (e.g. type of electrodes, pulse repetition rate).
In Fig. 8 we present the anisotropy ratio, that is the ratio between the conductivity in direction parallel with respect to the direction of the muscle fibres and the conductivity in direction perpendicular with respect to the direction of the muscle fibres, as a function of the elongation of the cell, or more precisely, of the cell diameter (fixed) as opposed to its length (varied) at the end of the 8 th pulse of the first train. By presenting this analysis/dependence we wish to justify the particular r short -to-r long cell dimension ratio that we chose and used for all the simulations of conductivity evolution, i.e. 1:6.7, as representing a sufficiently anisotropic cell geometry (in other words, a sufficiently long cell to represent a unit of muscle fibre from the electrical conductivity perspective) that does result in an appropriate degree of anisotropy of bulk muscle, facilitating the subsequent lesion shape analysis and thresholding processes. Fig. 8 illustrates that elongating the cell further, to ratios of 1:10 or 1:20, does not significantly impact the conductivity anisotropy ratio at the end of the 8 th (i.e. the last) pulse in the train.
We also analysed the impact of varying the volume fraction and the elongation of the unit-cell on the initial, preelectroporation conductivity of tissue in both parallel and perpendicular direction with respect to the direction of the long axis of the muscle fibre. Results of this analysis are given in Fig. S4 in the Supplementary Materials. While the dependence on volume Fig. 8. The effect of the geometry of the cell (r short -to-r long ratio), used in the unit-cell model, on the calculated anisotropy ratio (defined as AR = σ /σ ⊥ ). The results are shown for a train of eight 100 µs long pulses with a pulse delivery rate of 1000 s −1 . The magnitude of the applied electric field was 1000 V/cm. In (a) the electric field was applied in direction parallel with respect to the direction of the long axis of the muscle cell (parallel orientation), and in (b) in direction perpendicular with respect to the direction of the long axis of the muscle cell (perpendicular orientation). The insets show the anisotropy ratio during the first pulse. In (c), the unit-cell model results of anisotropy ratio at the end of the 8 th pulse for the parallel and the perpendicular orientation are given, together with the result of fitting an exponential function (f (r long ) = a (1 − exp(−r long /b)) + c) of r long to the unit-cell model results of anisotropy ratio. fraction is of negligible importance (Fig. S4b), the same cannot be said of the geometrical diameter-to-length ratio of the unit cell (Fig. S4a). Fig. 9 shows the number of pores formed on 1/8 of the muscle cell membrane as calculated with the unit cell model. At electric field values below 400 V/cm, more pores are formed when the electric field is applied in direction parallel to the muscle fibres than when it is applied in direction perpendicular to them, whereas at field values above 400 V/cm, more pores are formed in perpendicular orientation, which may explain different orientation sensitivity reported in the literature [36], [37].
The calculated values of electric current ranged from 1.97 A to 4.58 A, while the measured currents (of the first pulse train) ranged from 3.09 A to 11.70 A (all values are listed in Table I). The measured currents were consistently higher than the calculated ones. This may be due to unaccounted for phenomena related to high-voltage pulse delivery to in vivo animal tissue, which is far more complex than a mere assembly of unit cells [68]. The mean values of the measured electric current for each applied voltage are shown graphically in Fig. S5 in the Supplementary Materials. In spite of large scatter of measured currents -which is not surprising -higher voltages applied resulted in larger currents measured. When comparing the current of the first pulse in different trains, there is no statistically significant difference between the current values at all voltages, which allowed us to model only one train of pulses, rather than six delivered experimentally, to reduce computational cost. Through comparison of model results and in vivo experimental results, we believe to have successfully demonstrated the power of this approach in explaining especially the particular lesion shape that results from the highly geometrically asymmetrical muscle cell and consequential anisotropic properties of muscle tissue. Our study's main results include the determination of the field strength threshold required for irreversible electroporation of porcine skeletal muscle cells using the specific pulse protocol that was used, and the finding that the field-to-fibre orientation does not seem to affect the irreversible field strength threshold. We have thus successfully demonstrated that it is not the susceptibility of individual cells to damage caused by the magnitude of the electric field strength, but rather the anisotropic electric properties of tissue that result in the disparately shaped lesions whose shape strongly depends on field-to-fibre orientation.
The reader should beware our study concerns a single type of tissue, i.e., the porcine vastus lateralis skeletal muscle, and was performed in both the numerical as well as experimental part employing a specific electroporation (pulse) protocol. The obtained threshold in both magnitude and its non-dependence on field-to-fibre orientation must be interpreted in this light. Changing the muscle tissue type (e.g. cardiac muscle) or pulse protocol (e.g. pulse length, repetition rate, number, etc.) will most likely result in different thresholds and necessitates a repeat of the numerical study, with, ideally, a validation in form of experimental in vivo work. Furthermore, the calculated temporal evolution of conductivity during electroporation, calculated using the unit cell approach (Fig. 3), does not show the immediate drop in conductivity after the pulse that was observed in the experiments [69]. The reason for this shortcoming is that our unit-cell model is based on the models of Huclova et al. [46], and does not account for all the dynamics of recovery in conductivity. We were therefore unable to account for the immediate drop in conductivity after the pulse in the current developmental phase of our multiscale model.
We would also like to emphasise that in our study we used the values for intracellular (0.55 S/m) and extracellular (1.2 S/m) conductivity at the upper end of the ranges available in the literature. We chose these values to achieve comparable currents in the numerical electrical model to those observed in vivo. However, we also performed a parametric study in which we modelled with lower conductivity values; as low as 0.3 S/m for intracellular and 1.0 S/m for extracellular conductivity. We found that the impact on the main result of the study (i.e., the determined irreversible electric field strength threshold) was negligible. The threshold values for lower conductivity values differed no more than 5 % from the values obtained using the higher conductivity (i.e., 0.55 and 1.2 S/m). In absolute numbers, this means the average value when using lower conductivity resulted in a change in the irreversible threshold from 193.4 to 197.7 V/cm. Moreover, we can offer some theoretical support for the use of the higher-conductivity values for the extracellular space with negligible consequences on the calculated field strength thresholds by reasoning that the extracellular conductivity is of significant importance only when close to the tissue conductivity, i.e., when the tissue / cell suspension consists of only a few cells. This was extensively studied and demonstrated, for instance, in [70], where the authors compared the results of the FEM simulations with the analytical results of Maxwell, Rayleigh, Tobias, and Bruggeman. Their results of normalised effective conductivity versus the volume fraction (Fig. 5 in the aforementioned study) show that the effective conductivity at a cell volume fraction of 0 (where there are no cells in the medium) is equal to the extracellular conductivity, and then steadily decreases with increasing cell volume fraction. The effective (bulk tissue) conductivity amounts to only about 20 % of the extracellular conductivity at a cell volume fraction of cells of 0.76, which is the volume fraction used in our study. Thus, at a relatively high-volume fraction of 0.76, the intracellular conductivity has a much stronger influence of the bulk tissue conductivity as does the extracellular; meaning that even if we had overestimated the extracellular conductivity, the impact of such an overestimation should not be critical to the validity of the model-based data and derived conclusions.

IV. CONCLUSION
This paper presents a fresh attempt at implementing the very first steps in building a comprehensive model of electroporated skeletal muscle tissue by accounting for its anisotropy. We chose the ground-up multiscale approach to the study of macroscopic conductivity and its changes due to electroporation, whereby we try to deduce electrical properties and irreversible damage observed at the macroscale based on properties of single elongated muscle cells and their behaviour under exposure to electric field.
While we cannot -at this stage of development -claim the model is all-embracing and complete, as there are indications that unaccounted for phenomena and mechanisms exists that need to be either further elucidated and explained, and then included as co-factors into the model (e.g. the dynamics of current, mechanisms of cell death, etc.), we nevertheless believe it has the capability of serving as the basic proof of concept for this particular multiscale approach, and thus strongly warrants and invites further validation, exploration, elaboration, and consequentially exploitation.