# Sensitivity Analysis in a Rescaled 0D Pacemaker-Bimembrane Model

Valentin Pannetier<fnref n="1" />, Michael Leguèbe<fnref n="1" />, Yves Coudière<fnref n="1" />, Gwladys Ravon<fnref n="1" />, Delphine Feuerstein<fnref n="3" />, Guilhem Fauré<fnref n="3" /> and Richard Walton<fnref n="2" />

## Abstract

We present a new 0D model of a pacemaker device coupled to a surrogate model of cardiac tissue, based on the work by Pannetier et al. [1]. This model utilizes a simplified internal circuit of an implantable cardiac pacemaker device and includes equivalent models for electrode interfaces and a rescaled minimal myocyte model. Our model integrates the polarization effects commonly observed at electrodes-tissue interface, and allows for the simulation of pacemaker interactions with cardiac tissue, providing insights into energy consumption and device functionality.

## 1. Introduction

Pacemakers are critical devices that regulate heart rhythm by delivering electrical pulses through electrodes implanted into the cardiac tissue. These pulses are characterized by a duration and an amplitude which, when appropriately selected, provoke a cardiac beat. This is called capture by the manufacturers.

As the lifespan of pacemaker batteries is limited, understanding capture and its impact on energy consumption is essential to maximize the durability of the device. We propose an efficient computational model of a pacemaker coupled to a surrogate model of cardiac tissue, based on our previous work [1]. The model is built by simplifying the electronics of a real device into a circuit which allows us to study energy usage during a stimulation cycle. We use the model to compute the so-called Lapicque curve in the stimulation (duration, amplitude) space. The curve splits this space into two regions of functional and non-functional stimulations, respectively.

By exploring the parameters of our mathematical model, we aim to provide valuable insights that can help calibrate this model to match experimental data. This calibration is crucial for assessing the credibility of the model and ensuring its accuracy.

## 2. Model of a pacemaker coupled to a 2-membranes myocyte model

Our new 0D model is built upon an internal simplified circuit of a pacemaker and includes two equivalent models describing the electrode interfaces and a rescaled minimal myocyte model used as a surrogate for cardiac tissue. The circuit for this model is shown in Figure 1.

As in [1], we consider a stimulation signal divided into four parts, and based on a commercial device (Borea™, MicroPort CRM). First, an initial pulse is delivered that is aiming to trigger cellular excitation by discharging a large capacitor. Then, a reverse signal is applied to counteract the polarization which occurs on the surface of the electrodes. This is made possible by two switches, called $S_p$ and $S_o$, which swap the sub-circuits inside the device. To account for polarization [2], the electrodes are represented by R/C equivalent circuits (a resistor in parallel with a capacitor). We chose this model because it is the most descriptive model for which the parameters are still identifiable by calibration against experimental data [1].

The tissue is considered through a myocyte model incorporating two rescaled equivalent membranes based on the Beeler-Reuter ionic model [3], and connected by two conductances at cellular level, for the intracellular and extracellular medium.

### 2.1. Mathematical model

Let $I$ be the current flowing through the pacemaker circuit and $V_±$ the two electrode interfaces voltages. The usual Kirchhoff's laws on current yields the following equations [1] at electrodes,

$$C_{\pm} \frac{dV_{\pm}}{dt} + \frac{V_{\pm}}{R_{\pm}} = I, \quad (1)$$

with $R_±$ and $C_±$ the resistances and capacitances of equivalent circuits. We define the time constants $\tau_± := R_±C_±$.

Additionally, the surrogate tissue part of the circuit verifies the following equations on transmembrane voltages

<footnote n="1">¹ Univ. Bordeaux, CNRS, Inria, Bordeaux INP, IMB, UMR 5251, IHU Liryc, 33400 Talence, France</footnote>

<footnote n="2">² Univ. Bordeaux, Inserm, CRCTB, U 1045, IHU Liryc, 33000 Bordeaux, France</footnote>

<footnote n="3">³ Microport CRM, 92140 Clamart, France</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1. Electrical circuit describing the model. Left side (black) the device internal circuitry, center (green) the equivalent R//C model of contacts, and right side (red) the two membrane rescaled myocyte equivalent model containing extracellular and intracellular equivalent conductances.</figcaption></figure>

(TMV) $V_{m,1}$, $V_{m,2}$ and ionic states $h_1, h_2$,

$$C_m \frac{dV_{m,1}}{dt} + I_{ion,1} = \frac{G_m}{S_m} \left( V_{m,1} - V_{m,2} - \frac{I}{G_e} \right), \quad (2)$$

$$\frac{dh_1}{dt} + f_1 (h_1, V_{m,1}) = 0, \quad (3)$$

$$C_m \frac{dV_{m,2}}{dt} + I_{ion,2} = \frac{G_m}{S_m} \left( V_{m,1} - V_{m,2} + \frac{I}{G_e} \right), \quad (4)$$

$$\frac{dh_2}{dt} + f_2 (h_2, V_{m,2}) = 0, \quad (5)$$

where $G_m := G_i G_e (G_e + G_i)^{-1}$ denotes the equivalent conductance of extracellular and intracellular media in mS, $C_m$ is the membrane capacitance in $\mu F cm^{-2}$, $S_m$ is a scaling coefficient in $cm^2$, and with Beeler-Reuter ionic models ($I_{ion,1}, f_1$) and ($I_{ion,2}, f_2$).

Finally, the system is closed by writing Kirchhoff's law on voltage on the device part of the circuit,

$$R_k \frac{dI}{dt} + \frac{I}{C_k} = \frac{V_-}{\tau_-} + \frac{V_+}{\tau_+} + \frac{G_m}{G_e} \left( \frac{dV_{m,1}}{dt} - \frac{dV_{m,2}}{dt} \right), \quad (6)$$

where $R_k$ and $C_k$ are the resistance and capacitances of the device circuit depending on states of switches $S_p$ and $S_o$.

The complete system is solved in our CirCE<fnref n="1" /> software.

## 2.2. Experimental data

Mathematical model design and parametrisation was guided by experiments using ex vivo sheep hearts (N=7, aged 1-2 years old and body weights of 43-52 kg).

Dual coronary-perfused ventricular wedge preparations from hearts opened along the posterior left and right ventricles and posterior septum were stretched on support frames

for ex vivo assessment into a bath of saline solution, as described in [4] (Fig. 2).

<figure><img src="image_2.png" /><figcaption>Figure 2. Tissue wedge preparation with the pacemaker electrode in place into the second experiment protocol.</figcaption></figure>

Controlled electrical pulses were applied initially with fixed pulse durations but variable amplitudes that decreased every 8 pulses, starting above the capture threshold. The amplitude decreased until absence of capture is detected by the device. Repeating this process for several durations yields the experimental Lapicque curve.

However, we observed that the pacemaker device is unsuitable due to its limited temporal and voltage resolutions, which are insufficient for an accurate localization of the Lapicque curve. Instead, using a PSA (Pulse System Analyzer) combined with electrical monitoring of the tissue and the lead (PowerLab systems), allowed to test of a larger band of amplitudes and durations with more steps. This setup also allows reaching 100 kHz temporal data acquisition for precise observation of events, such as the rapid pacemaker pulse events.

<footnote n="1">¹https://carmen.gitlabpages.inria.fr/circe/</footnote>
### 3. Global sensitivity analysis

To effectively compare experimental and numerical Lapicque curves, it is essential to establish a good criterion to characterize the capture phenomenon. Whereas experimentally we only have the voltage measured at electrodes, the model gives direct access to the TMVs to define more criteria. Consequently, we introduced several specific quantities derived from both TMVs and voltage $V_{meas}$, as measured during the experiment at the pacemaker entries (MV). These quantities should be able to provide insights into the membrane's behavior, and local/global response to stimulations.

After studying several quantities on voltages (TMVs and MV), we kept the following observables:

* baseline-adjusted integrals; the baseline is defined as a first value of predefined window starting 5 ms before each stimulation,

* values at end of measurement window; these values typically correspond to the state of the model just before the next stimulation,

* duration between the peak of TMVs and the average distance between it and TMVs baseline (D50), corresponding to the APD50 for capturing stimulations,

* extrema of voltages and of their time-derivatives,

* relative squared differences on pulse, switch and discharge parts of the cycle compared with experimental data obtained on sheep hearts.

We conducted a Sobol sensitivity analysis [5] using the SALib<fnref n="2" /> Python library, to measure the effects of the most uncertain parameters on the observables. The sampling technique and bounds of tested parameters are chosen carefully, as we have no significant a priori estimates for them. The bounds (see table 1) were extended compared to the previous ones [1] since the model was partially modified. The model was evaluated for 786, 432 samples of parameters, but the simulation failed for 2% of them due to physically incompatible values. The evaluation of these failed samples is replaced by the average evaluation of those remaining during the analysis of the Sobol indices. The results are shown on figure 3.

<table><caption>Table 1. Bounds of parameters used in the Sobol sensitivity analysis of the 0D model.</caption><thead><tr><th>Name</th><th>Min</th><th>Max</th><th>Unit</th></tr></thead><tbody><tr><td>C<sub>ets</sub>, C<sub>ts</sub></td><td>10<sup>-2</sup></td><td>10<sup>3</sup></td><td>μF</td></tr><tr><td>R<sub>p</sub>, R<sub>int</sub></td><td>10<sup>-4</sup></td><td>10<sup>1</sup></td><td>kΩ</td></tr><tr><td>R<sub>±</sub></td><td>10<sup>-4</sup></td><td>10<sup>3</sup></td><td>kΩ</td></tr><tr><td>τ<sub>±</sub></td><td>10<sup>-3</sup></td><td>10<sup>2</sup></td><td>ms</td></tr><tr><td>G<sub>e</sub>, G<sub>i</sub></td><td>10<sup>-4</sup></td><td>10<sup>2</sup></td><td>mS</td></tr><tr><td>S<sub>m</sub></td><td>10<sup>-1</sup></td><td>10<sup>3</sup></td><td>cm<sup>2</sup></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 3. Results of total order indices of Sobol analysis on several outputs features on 786, 432 samples. The three top blocks report the indices for observables derived from $V_{m,1}$, $V_{m,2}$ and $V_{meas}$ respectively. The lower block is dedicated to $L^2$ relative difference between the model and experimental data.</figcaption></figure>

The analysis shows a strong influence of the tissue parameters ($S$, $G_i$, $G_e$) on the various characteristics of the TMVs, particularly for the first membrane. They also have a significant influence on the integral and on the last value of $V_{meas}$. Surprisingly, the capacitance value $C_{ets}$ has little influence on these outputs. Additionally, we can also note the parameters that are not identifiable from the experimental data, such as the capacitance $C_{ets}$ and the two resistances of the pacemaker circuit $R_p$ and $R_{int}$. These parameters should be excluded from any forthcoming calibration procedure.

Finally, this study also provides information on which observables can be used to detect capture. For both TMVs, all quantities are ranked at the same level, but only the integral and the last value of $V_{meas}$ seem to be impacted by parameter changes.

<footnote n="2">²https://salib.readthedocs.io</footnote>
4. Hysteresis of stimulation threshold detection

To further investigate the issue of capture characteriza-
tion, we used the model to reproduce the hysteresis that is
observed when searching the stimulation threshold using
either ascending or descending amplitudes. This is com-
monly observed in [6,7]. It became apparent that although
all the criteria can be used to identify the transition from
non-capturing to capturing stimulation, the criteria based
on the TMVs derivatives or the Vmeas integral do not have
good mathematical properties and do not define a clear
threshold. In figure 4, we plot the remaining quantities
defined in the previous section, depending on the stimula-
tion amplitude during a threshold search, as performed by
the pacemaker.

Different thresholds are observed depending on whether
the stimulation amplitude is increasing or decreasing. This
confirms that our model is able to reproduce the observed
hysteresis. Except for the D50 on the first membrane, all
quantities show a large variation around the same threshold
amplitude. This should tell that any of these quantities is a
good indicator of capture. On the contrary, the values at the
end of the stimulation cycle show a very small difference
compared to the amplitude of signal, less than 1% for the
TMVs (not plotted), and even lower than 10⁻⁶ for the MV.
This suggests that other observables, such as the maximum
or the D50, are safer to characterize capture.

<figure><img src="image_4.png" /><figcaption>Figure 4. Illustration of hysteresis in the 0D model of 4 quantities. The amplitude varies on x-axis.</figcaption></figure>

**5. Conclusion**

The hysteresis effect observed numerically underlines
the model’s ability to reproduce the capture phenomenon
on different observables, even for quantities absent from
the experimental data (TMVs). However, high Sobol in-
dices for these markers indicate that calibration of the pa-
rameters with respect to data must be performed with care
to ensure the validity of the model.

Acknowledgments

This work was supported by the H2020 EU SimCar-
dioTest project (Digital transformation in Health and Care
SC1-DTH-06-2020; grant agreement number 101016496).
This study received financial support from the French Gov-
ernment as part of the “Investments of the Future” program
managed by the National Research Agency (ANR), Grant
reference ANR-10-IAHU-04. Experiments presented in
this paper were partially carried out using the PlaFRIM ex-
perimental testbed, supported by Inria, CNRS (LABRI and
IMB), Université de Bordeaux, Bordeaux INP and Conseil
Régional d'Aquitaine.

**References**

[1] Pannetier V, Leguèbe M, Coudière Y, Walton R, D'Hiver P, Feuerstein D, Amaro D. Modeling cardiac stimulation by a pacemaker, with accurate tissue-electrode interface. In International Conference on Functional Imaging and Modeling of the Heart. Springer, 2023; 194–203.

[2] Martinsen OG, Heiskanen A. Bioimpedance and bioelectricity basics. Elsevier, 2023.

[3] Beeler GW, Reuter H. Reconstruction of the action potential of ventricular myocardial fibres. The Journal of physiology 1977;268(1):177-210.

[4] Walton RD, Pashaei A, Martinez ME, Constantin M, Duchateau J, Bear L, Cros C, Pascarel-Auclerc C, Guo Y, Benoist D, et al. Compartmentalized structure of the moderator band provides a unique substrate for macroreentrant ventricular tachycardia. Circulation Arrhythmia and Electrophysiology 2018;11(8):e005913.

[5] Soból I. Sensitivity estimates for nonlinear mathematical models. Math Model Comput Exp 1993;1:407.

[6] Huang X, Qian Y, Zhang X, Hu G. Hysteresis and bistability in periodically paced cardiac tissue. Physical Review E—Statistical Nonlinear and Soft Matter Physics 2010; 81(5):051903.

[7] Malinowski K, Unger P, RES J. Evaluation of different hys-
teresis functions-a study with the biotronik actros dr. Prog
Biomed Res 1999;2(4):154-159.

Address for correspondence:

Valentin PANNETIER
351, cours de la Libération - F 33 405 TALENCE
valentin.pannetier@math.u-bordeaux.fr