ERK'2021, Portorož, 415-418 415 XXX-X-XXXX-XXXX-X/XX/$XX.00 ©20XX IEEE Time-Dependent Numerical Model of The Induced Transmembrane Voltage in real-shaped Cardiomyocyte Maria Scuderi, Damijan Miklavčič, Janja Dermol-Černe University of Ljubljana, Faculty of Electrical Engineering Ljubljana, Slovenia maria.scuderi@fe.uni-lj.si, Abstract —Pulsed-field ablation (PFA) is a nonthermal new promising method used to treat heart arrhythmias with the application of a sufficiently high pulsed electric field to the target tissue. Different orientations and elongated shapes of cardiomyocytes may hinder obtaining the homogeneous distribution of the applied electric field and homogeneous ablation of tissue. To increase the efficiency of PFA treatments, numerical models were developed to study the effects of the electric field on a cardiomyocyte. Cardiomyocyte shape had usually been approximated as a prolate spheroid cell or a cylinder but recently, Milan et al. [1] for the first time compared different cardiomyocyte approximations to the realistic shape of a cardiomyocyte at steady-state conditions. They suggested prolate spheroid geometry as a good-enough and numerically less time-consuming approximation of the real-shaped geometry of a cardiomyocyte. The aim of our research work is to further advance the existing real-shaped cardiomyocyte model by developing a time-dependent numerical model of a cardiomyocyte comparing prolate spheroid and a real-shaped geometry to evaluate if it is reasonable to use prolate spheroid geometry to represent a cardiomyocyte in time-dependent numerical studies. Keywords —Transmembrane electric potential, electric field, myocardial cell, electroporation 1 Introduction If an electric field is applied to a biological cell, transmembrane voltage (ITV) is induced across its plasma membrane and is superimposed to the physiological resting transmembrane voltage. The resting transmembrane voltage is in the range of -70 mV to - 40 mV for eukaryotic cells [2]. Up to the membrane breakdown, the ITV is proportional to the amplitude of the applied electric field and depends on the cell orientation with respect to the applied electric field [3,4]. ITV depends, among others, also on the density of cell suspension, the cell geometry, position on the cell membrane, and extracellular medium conductivity. When a sufficiently high pulsed electric field is applied to the biological cells, an increase of cell membrane conductivity and permeability is detected and this phenomenon is called electroporation (EP). Electroporation occurs on the regions of the cell membrane in which the ITV is above 0.2-1 V. [3]. Conventionally, heart arrhythmias are treated with targeted ablation of arrhythmogenic regions in the heart with radiofrequency and cryoablation in which thermal energy is used to induce t tissue necrosis. However, the application of thermal energy causes pulmonary vein stenosis, mitral valve trauma, and can cause phrenic nerve, and oesophageal injury [5]. Thus, there is the need to find a safer ablation technique. Recently, electroporation i.e. pulsed-field ablation (PFA) emerged as a non-thermal ablation method to treat heart arrhythmias [5,6]. In PFA, nonthermal energy is used to ablate cardiac tissue, promoting cell apoptosis, preserving the nearby tissue like nerves and vessels, and allowing the complete and durable isolation of the pulmonary vein and other structures in the heart depending on the aim [7]. However, different orientations and the cylindrical shape [8] (its diameter is 4 times langer than its length) of cardiomyocytes with respect to the applied electric field represents a challenge to achieve a homogeneous electric field in the area to be treated [4,9]. Thus, to increase the efficiency of the treatment, it is important to study the effects of the applied electric field on cardiomyocytes in time domain especially when short pulses, sometimes shorter than the membrane charging time, are applied. In 2019, Milan et al., for the first time, developed and validated a 3D steady-state numerical electromagnetic model (NMReal) of a realistic geometry of a single mammalian ventricular cardiomyocyte [1]. In their study, they compared and tested simplified cardiomyocyte geometries models (e.g. cylinder, ellipsoid, prolate spheroid, and parallelepiped) with the realistic shape geometry (NMReal) of a cardiomyocyte. In conclusion, they estimated an average error of < 5 %, when the electric field was oriented <30° with respect to the long axis of the cell, between the prolate spheroid model and the real- shaped model. Thus, the prolate spheroid model could be considered as a good-enough and much simpler alternative to the real-shaped NMReal model [1]. However, it was not determined if the steady-state simplifications used by Milan et al. [1] also hold if short pulses are used. Thus, the aim of this study was to develop a time-dependent numerical model of the ITV of a single cardiomyocyte comparing the prolate spheroid (simplest geometry) and the realistic geometry (NMReal, the most complex geometry). 2 Methods In the following subsections, the methods to develop and validate the models are described. 2.1 ITV Electromagnetic models 416 The finite element model of a cardiomyocyte, exposed to an external pulsed electric field, was constructed in COMSOL Multiphysics (v5.5, COMSOL AB, Sweden) for the stationary and time-dependent case. First, the numerical model was developed at the stationary condition to compare the results with Milan et al paper [1]. Then, the model was adapted for the time-dependent case when a short pulse was applied. The 3-dimensional simulation was made in the AC/DC module, Electric current physics, and Time domain study by solving the Laplace equation (1) −𝛻 (𝜎 𝑖 𝛻 𝑉 ) − 𝛻 𝜕 (𝜀 𝑖 𝛻 𝑉 ) 𝜕 𝑡 = 0 (1) where 𝜎 𝑖 is the conductivity and 𝜀 𝑖 is the dielectric permittivity of the extracellular and intracellular medium and cell membrane subdomains of the model. For the numerical model at the stationary condition, the Steady State study was used and the second term of equation 1 was neglected. Two different geometries, prolate spheroid geometry (120 µm long with a diameter of 30 µm) and the real- shaped geometry (NMReal) of a cardiomyocyte (143 µm length and 36 µm width and 22 µm height) were used to represent a single cardiomyocyte as in Milan et al. [1], Figure 1. The real-shaped geometry of a cardiomyocyte was kindly provided by Milan et al. [1]. Both geometries were modeled in the center of the simulation square cube with dimensions 400 x 400 x 400 μm. The box size did not affect the results. A finer mesh was used for the cell membrane, which is the most important part of the geometry for the ITV analysis. The normal mesh was used for the other parts of the geometry. The voltage was applied as a trapezoidal electric pulse and as a constant value on the two opposite boundaries of the square block for the time-dependent and steady-state numerical model, respectively. The remaining four faces of the block were modeled as insulating surfaces. The pulse had to be trapezoidal, to avoid convergence problems. The electric pulse was obtained by subtracting two Heaviside functions using the COMSOL function flc1hs [10]. The pulse rise time was 1 μs and the pulse length was 100 μs. The electric field was applied perpendicular and parallel to the long axis of the cardiomyocyte for both geometries and was 1 V/cm as in Milan et al. paper [1]. The ITV was calculated as the difference between the extracellular, V e, and intracellular, V i, electric potential, ITV=V e-V i, for both geometries. The maximum absolute value of ITV, ITV max, on the cell membrane is presented in the regions of the cell membrane facing the cathode and anode. For a prolate spheroid geometry, the ITV max was evaluated at the point of the cell membrane facing the cathode or anode. For the real-shaped geometry, it was difficult to identify the point at which the ITV is the highest due to the irregularity of the shape. Thus, the ITV max was evaluated by the Boundary Probe type Maximum in Comsol, which allows determining the maximum value of the ITV on the cell membrane as shown in Figure 3. The cell membrane was modeled as a conductivity surface at the interface between the cell interior and the exterior, using the Contact Impedance Boundary Condition 𝒏 · 𝑱 𝟏 = 1 𝑑 𝑚 (𝜎 𝑚 + 𝜀 0 𝜀 𝑚 𝜕 𝜕𝑡 ) (𝑉 1 − 𝑉 2 ) 𝒏 · 𝑱 𝟐 = 1 𝑑 𝑚 (𝜎 𝑚 + 𝜀 0 𝜀 𝑚 𝜕 𝜕𝑡 ) (𝑉 2 − 𝑉 1 ) (2) where n is the normal vector, J is the current density, 𝑑 𝑚 is the cell membrane thickness, 𝜎 𝑚 is the cell membrane conductivity, 𝜀 𝑚 is the cell membrane permittivity, 𝜀 0 is the permittivity of the vacuum, V 1 and V 2 are the electric potentials at the outer and inner surface of the membrane in equation (2). For the numerical model at the steady- state condition, the time-dependent permittivity term of equation 2 was neglected. The parameters used in the numerical model are the same as in Milan et al paper [1] and are given in Table 1. Table 1. Parameters used in the models Parameter Symbol Value Intracellular Permittivity ε i 80 Extracellular Permittivity ε e 80 Membrane Permittivity ε m 5 Intracellular Conductivity σ i 0.8 [S/m] Extracellular Conductivity σ e 1.4 [S/m] Membrane Conductivity σ m 1.493x10 -8 [S/m] Membrane thickness d m 10 [nm] Block length l 400 [µm] Figure 1. Geometries used to represent a cardiomyocyte. 1A) Prolate spheroid geometry, 120 µm long and 30 µm wide and high. 1B) Real-shaped geometry (NMReal model) 143 µm long, 36 µm wide, and 22 µm high. 2.2 Model Validation The estimation of the ITV relative error between the numerical model and the analytical solution using a prolate spheroid geometry was evaluated to first validate the model only for the stationary state. The relative error was determined by solving the following equation 𝑅𝑒𝑙𝑎𝑡𝑖𝑣𝑒 𝑒𝑟𝑟𝑜𝑟 = |𝑣 𝑎 − 𝑣 𝑒 | 𝑣 𝑒 𝑥 100% (3) 417 where 𝑣 𝑎 is the observed value (ITV with the numerical model) and 𝑣 𝑒 is the exact value (with the analytical equation). The analytical equation that allows calculating the ITV for a prolate spheroid geometry is : 𝐼𝑇𝑉 (𝜑 ) = 𝐸 𝑅 1 2 − 𝑅 2 2 𝑅 1 − 𝑅 2 2 √𝑅 1 2 − 𝑅 2 2 log 𝑅 1 + √𝑅 1 2 − 𝑅 2 2 𝑅 2 𝑥 𝑅 2 cos (𝜑 ) √𝑅 1 2 𝑠𝑖𝑛 2 (𝜑 ) + 𝑅 2 2 𝑐 𝑜𝑠 2 (𝜑 ) (3) where E is the external electric field, R 1 is the radius along the axis of rotational symmetry, R 2 is the radius perpendicular to this axis, φ is the polar angle measured from the center of the cell concerning the direction of the field [11]. For the time-dependent state, an estimation of the ITV relative error between the prolate spheroid, 𝑣 𝑎 , and the real-shaped numerical model, 𝑣 𝑒 was evaluated. 3 Results and discussion Induced transmembrane voltage was calculated for a prolate spheroid, and a real-shaped geometry, in stationary, Figures 2 and 3, and time-dependent Figures 4 and 5 conditions. The electric field, 1 V/cm, well beyond the threshold for electroporation, was applied parallel and perpendicular to the long axis of the cell for both geometries. The direction of the electric field is indicated by plus and minus signs. A higher value of the ITV is present when the electric field is applied parallel to the long axis of the cell (~8 mV, Figure 4A) than when it is applied perpendicularly (3 mV, Figure 4B), which corresponds to already published studies [10,13]. This shows that the induced transmembrane voltage depends on the cell orientation with respect to the applied electric field [2]. For the steady-state case, the ITV values showed in Figure 3 are similar to the ones published in Milan et al. paper [1]. The relative error of ITV max estimation between the prolate spheroid and the analytical model, when the applied electric field was parallel to the long axis of the cells, is <1%. ITV max was defined as the maximum absolute value of ITV presents in the regions of the cell membrane facing the anode and cathode. The estimated ITV max relative error <1% is in agreement with Milan et al paper [1]. For the time-dependent case, the induced transmembrane voltage was evaluated as a function of time when a single trapezoidal pulse 100 µs long was applied. In Figures 4 and 5, the blue and red curves showed the values of the induced transmembrane voltage as a function of the time using prolate spheroid and real- shaped geometry, respectively. After a few microseconds, a steady-state value of the ITV is obtained, Figure 4. The time required to reach a steady-state condition of the ITV for a real-shaped model is longer than with using the prolate geometry Figure5. This might be due to the irregularity of the real-shaped geometry of the cardiomyocyte compared to the prolate one. The solid curves represent the ITV value evaluated when θ=0° and the dashed line for θ=90°. θ was defined as the polar angle measured from the center of the cell concerning the electric field. Figure 2. Variation of the induced transmembrane voltage of a prolate spheroid is obtained applying 1 V/cm electric field parallel (Fig. 2A) and perpendicular (Fig. 2B) to the long axis of the cell. The direction of the electric field is indicated by plus and minus signs. The legend, in both Figures, indicates the values of ITV in mV. Figure 3. Variation of the induced transmembrane voltage of a mammalian cardiomyocyte is obtained applying 1 V/cm electric field parallel (Fig. 3A) and perpendicular (Fig. 3B) to the long axis of the cell. The direction of the electric field is indicated by plus and minus signs. The legend, in both Figures, indicates the values of ITV in mV. 418 The relative error of ITV max estimation comparing the prolate spheroid and the realistic shape geometry was 17% and 15.7% when the electric field was applied parallel and perpendicular to the long axis, respectively. The prolate spheroid was under- or overestimating the results obtained with the realistic geometry. The obtained estimations of the relative error of ITV max are different compared to Milan et al. paper [1]. The reason for this difference might be due to numerics and the use of a different finite element simulation software (Milan et al. used Elmer FEM) than Comsol Multiphysics used in this study. Besides, the maximum value of the ITV in Milan et al. paper is not specified in how it was evaluated [1]. Thus, different ways the evaluate the maximum value of the ITV might affect the final results. Figure 4. The induced transmembrane voltage (ITV) as a function of time when a single pulse 100 µs long was applied. The blue and red curves represent the ITV values obtained with a prolate spheroid and real-shaped geometry, respectively. The ITV values were evaluated at the pole (θ=0°) and the equator (θ=90°) of the cell. 1 V/cm electric field was applied parallel (Fig. 4A) and perpendicular (Fig. 4B) to the long axis of the cell. Figure 5. Zoom-in of Fig. 4 to the first 10 µs of the simulation. It shows the charging time of the cell membrane using the prolate spheroid and real-shaped geometries. 4 CONCLUSION The aim of this research work was to develop a time- dependent numerical model of a cardiomyocyte and establish if the prolate spheroid geometry is accurate enough to represent a single mammalian cardiomyocyte in calculating the induced transmembrane voltage. The ITV max value obtained using prolate spheroid and real- shaped geometry of a cardiomyocyte is in the same range and the relative error of ITV max is 17% and 15.7%, when an electric field is applied parallel or perpendicular to the long axis of the cell, respectively. The computing time of the prolate spheroid model is 3 or 4 times faster than the real-shaped model when the electric field is applied parallel or perpendicular to the long axis of the cell, respectively. Thus, the prolate spheroid geometry might be considered as a good-enough approximation to represent a cardiomyocyte, taking into account the decrease in the computational load in comparison with the realistic geometry in numerical models. Based on the ITV results, we need to apply a higher electric field when the cells are oriented perpendicular than when are oriented parallel to it to obtain similar results [4,9]. The ITV, without electroporation, increases proportionally the applied electric field but with electroporation, the ITV increase is non-linear due to membrane breakdown. Acknowledgment This work was supported by the research funding from Medtronic PLC and the Slovenian Research Agency (ARRS). References [1] H. F. Milan, R. A. Bassani, L. E. Santos, A. C. Almeida, and J. W. Bassani, “Accuracy of electromagnetic models to estimate cardiomyocyte membrane polarization,” Medical & biological engineering & computing, vol. 57, no. 12, pp. 2617–2627, 2019 [2] T. Kotnik, L. Rems, M. Tarek, and D. Miklavčič, “Membrane electroporation and electropermeabilization: mechanisms and models,” Annual review of biophysics, vol. 48, pp. 63–91, 2019. [3] G. Pucihar, T. Kotnik, B. Valič, and D. Miklavčič, “Numerical determination of transmembrane voltage induced on irregularly shaped cells,” Annals of biomedical engineering, vol. 34, no. 4, p. 642, 2006. [4] J. Dermol-Černe, T. B. Napotnik, M. Reberšek, and D. Miklavčič, “Short microsecond pulses achieve homogeneous electroporation of elongated biological cells irrespective of their orientation in electric field,” Scientific Reports, vol. 10, no. 1, pp. 1–17, 2020. [5] S. McBride et al., “Ablation Modalities for Therapeutic Intervention in Arrhythmia-Related Cardiovascular Disease: Focus on Electroporation,” Journal of Clinical Medicine, vol. 10, no. 12, p. 2657, 2021. [6] C. J. Bradley and D. E. Haines, “Pulsed field ablation for pulmonary vein isolation in the treatment of atrial fibrillation,” Journal of Cardiovascular Electrophysiology, 2020. [7] J. S. Koruth et al., “Pulsed Field Ablation Versus Radiofrequency Ablation: Esophageal Injury in a Novel Porcine Model,” Circulation: Arrhythmia and Electrophysiology, vol. 13, no. 3, p. e008303, 2020. [8] P. Pham, G. Cauffet, A. Bardou, J. Olivares, and E. Novakov, “Development of a linear transient model for stimulation of isolated cardiac cells,” The European Physical Journal Applied Physics, vol. 12, no. 3, pp. 217– 222, 2000. [9] R. Ranjan and N. V. Thakor, “Electrical stimulation of cardiac myocytes,” Annals of biomedical engineering, vol. 23, no. 6, pp. 812–821, 1995. [10] L. Rems, M. Ušaj, M. Kandušer, M. Reberšek, D. Miklavčič, and G. Pucihar, “Cell electrofusion using nanosecond electric pulses,” Scientific reports, vol. 3, no. 1, pp. 1–10, 2013. [11] T. Kotnik and D. Miklavčič, “Analytical description of transmembrane voltage induced by electric fields on spheroidal cells,” Biophysical Journal, vol. 79, no. 2, pp. 670–679, 2000.