A Reaction-Diffusion Heart Model for the Closed-Loop Evaluation of Heart-Pacemaker Interaction

The purpose of this manuscript is to develop a reaction-diffusion heart model for closed-loop evaluation of heart-pacemaker interaction, and to provide a hardware setup for the implementation of the closed-loop system. The heart model, implemented on a workstation, is based on the cardiac monodomain formulation and a phenomenological model of cardiac cells, which we fitted to the electrophysiological properties of the different cardiac tissues. We modelled the pacemaker as a timed automaton, deployed on an Arduino 2 board. The Arduino and the workstation communicate through a PCI acquisition board. Additionally, we developed a graphical user interface for easy handling of the framework. The myocyte model resembles the electrophysiological properties of atrial and ventricular tissue. The heart model reproduces healthy activation sequence and proved to be computationally efficient (i.e., 1 s simulation requires about 5 s). Furthermore, we successfully simulated the interaction between heart and pacemaker models in three well-known pathological contexts. Our results showed that the PDE formulation is appropriate for the simulation in closed-loop. While computationally more expensive, a PDE model is more flexible and allows to represent more complex scenarios than timed or hybrid automata. Furthermore, users can interact more easily with the framework thanks to the graphical representation of the spatiotemporal evolution of the membrane potentials. By representing the heart as a reaction-diffusion model, the proposed closed-loop system provides a novel and promising framework for the assessment of cardiac pacemakers.


I. INTRODUCTION
Pacemakers are implantable medical devices commonly used to resolve cardiac arrhythmias when pharmacological interventions are not effective. Over the past decade, the use of Cardiac implantable electronic devices (CIEDs) has continued to increase, particularly for cardiac resynchronization The associate editor coordinating the review of this manuscript and approving it for publication was Sung-Min Park . therapy [1]. Recent studies have highlighted a non-negligible percentage of reports of malfunctioning of CIEDs [2]. Pacemakers are programmed with complex software able to cope with a great variability of pathological cardiac situations. In general, it is not always possible to fully predict the type of interaction that a particular device configuration may have on a patient's specific cardiac activity [3], [4]. The combination of abnormal cardiac events, such as premature ventricular contraction (PVC) or ventriculo-atrial conduction, and VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ incorrect device programming can cause serious problems for the patient, e.g., endless loop tachycardia (ELT). An openloop procedure for pacemaker testing consists of the evaluation of the device behaviour upon receipt of pre-recorded cardiac signals. The open-loop test is non exhaustive as it does not consider the heart response to external stimulation. Clinical experimentation, on the other hand, while constituting an essential phase for the validation of implantable devices, can only be carried out on a limited sample of cardiac pathologies, as well as entailing significant risks and costs. In recent years, numerous studies have addressed the development of closed-loop systems for the validation of pacemakers under physiologically relevant conditions [5], [6], [7], [8]. A closed-loop system consists in the interaction between a cardiac model, which simulates the generation and propagation of the action potential in the heart tissue, and the activity of a pacemaker. A closed-loop system can be useful for model-based design when implementing the pacemaker software and algorithms in simulation or hardware emulation. In addition, a closed-loop system can be used for the validation of the physical device. A closed-loop system for the design and validation of pacemakers consists of the pacemaker (system model, emulated device, or physical device), the physiological model of the heart and the necessary interfaces [5]. The pacemaker model receives and processes the electrograms generated by the cardiac model and generates pacing signals that act on the cardiac model to induce atrial or ventricular activation. In turn, the cardiac model reacts to the pacing stimuli with a physiologically significant response and, in accordance, generates the electrograms detected by the pacemaker. The most critical component of closed-loop systems is the cardiac model. The model must be able to interact with the pacemaker (i.e., react to stimuli and generate meaningful electrograms) guaranteeing a physiologically realistic response capable of reproducing a wide range of cardiac conditions (e.g., normal sinus rhythm, tachycardia, bradycardia, atrio-ventricular (AV) blocks, etc.). At the same time, the model must be as simple as possible with a small number of parameters that are easily adaptable to different heart conditions. Furthermore, the cardiac model must be computationally efficient to allow pacemaker verification in a reasonable time. Moreover, in the case of physical device validation, real-time operation is required.
In this work, we describe a two-dimensional reactiondiffusion model (i.e., based on partial differential equations) of the electrical activity of the heart specifically designed for closed-loop pacemaker design and validation. The model is based on a phenomenological myocyte model coupled with the cardiac monodomain equations. The phenomenological myocyte model was developed in our previous works originally describing epicardial tissue [9], [10], [11]. In this work, we adapted our myocyte model to describe atrial tissue and transmural heterogeneity of ventricular tissue. The heart model includes two stimulation and sensing electrodes to simulate the interaction with a dual-chamber cardiac pacemaker. We developed a closed-loop system composed by our heart model and a DDD pacemaker model. We deployed the system on a specific hardware setup consisting of an Arduino board emulating the cardiac pacemaker, a workstation for the cardiac simulation, and a PCI acquisition board for communication between the heart and pacemaker models. A DDD pacemaker has pacing and sensing capabilities in both the atrium and the ventricle. Indeed, DDD mode is widely used in patients with combined sinus node dysfunction and AV node dysfunction. In our framework the cardiac model delivers atrial and ventricular signals to the pacemaker and receives as inputs the atrial and ventricular pacing sent from the device. We tested the behaviour of the closed-loop system in three pathological cases, to prove its functionality and potentiality. First, we considered two wellknown conditions in which cardiac pacing is a standard treatment: sick sinus syndrome and AV block. Then, we employed our closed-loop system to simulate ELT. ELT can be defined as a reentrant tachycardia and represents a common complication of dual chamber pacing systems [12], [13]. Most ELT are initiated by PVC which can lead to a retrograde atrial activation in pathological conditions. Currently, pacemakers are provided with anti-ELT algorithms to detect and terminate the ELT [14], [15].
The main novelty of this work consists of the employment of a reaction-diffusion heart model in a closed-loop system aimed at the evaluation of heart-pacemaker interaction. Indeed, previous works on closed-loop models represented the heart as a network of timed automata (TA) [5], and later hybrid automata (HA) [7], [8].
The main advantages of our implementation are the possibility to model spatial inhomogeneities, such as myocardial fibrosis and transmural heterogeneity, and the generation of more accurate electrograms (EGMs). Moreover, our closedloop environment can be exploited for the design and optimization of electrodes shape and position. Additionally, our model offers a graphical user interface (GUI) that allows the user to interact with the heart model (e.g., to simulate pathologies) and visualize the electrophysiological activity of the heart, so that even non-expert users can interact with the framework.

A. HEART MODEL
We built a two-dimensional cardiac model in which the cardiac tissue is divided into 6 distinct regions: sinoatrial (SA) node, atrium, AV node, endocardium, midmyocardium and epicardium (Fig. 1). For the interaction with a dual chamber (DDD) pacemaker we introduced two stimulation and sensing sites: an atrial electrode placed near the right atrial auricle and a ventricular electrode placed in the apical area of the right ventricle (light grey areas of Fig. 1). The cardiac model presents spontaneous activation in the SA node, takes in input FIGURE 1. Geometry of the 2D heart model made of 6 distinct regions: SA node, atrium, AV node, endocardium, midmyocardium and epicardium. The light grey regions show the location of the pacing electrodes.
activation signals from the pacemaker, and gives in output the atrial and ventricular electrograms.

1) MYOCYTE MODEL
To describe the electrophysiological behaviour of atrial tissue and ventricular endocardium, midmyocardium and epicardium, we adapted the phenomenological model of epicardial tissue described in our previous works, which proved to be computationally efficient and easy to parameterize [9], [10], [11]. To obtain the steeper restitution curves that are experimentally observed in the endocardium and midmyocardium compared to the epicardium, we modified e 0 1 and e 2 by introducing a linear dependence on the state variable u: We also introduced a quadratic dependence between A and u for an accurate representation of the action potential dome in the four types of cells: Finally, we changed the d w parameter to prevent it from getting too high, thus facilitating numerical integration into the heterogeneous myocardium: The model parameters of the four types of tissue are reported in Table 1. Thanks to the low number of parameters and the easy understanding of their role in determining the dynamic behaviour of the system, we obtained the specific value of each parameter for the four tissue types by performing a twostep optimization procedure aimed at reproducing available experimental data on human heart tissue. In the first step, we determine the model parameters that reproduce the main characteristics of experimental action potential morphology, such as action potential amplitude, upstroke velocity, notch amplitude (if present), plateau voltage, and resting membrane potential [16], [17], [18], [19], [20], [21], [22], [23]. The first optimization step is straightforward since the mentioned action potential features are directly related to one or two parameters of the model. For the complete description of the  role of each parameter we refer the reader to [9]. In the second optimization step we adjusted the model parameters by fitting experimental action potential duration and conduction velocity (CV) steady-state restitution curves [19], [24], [25], [26], [27]. In particular, we used the Gauss-Newton method backtracking line search to find a constrained minimum of the least square deviation between simulated and experimental restitution curves. For each cell type, we calculated the restitution curves on a cable model as described in [9].

2) NUMERICAL METHODS
To simulate the action potential propagation during pacemaker-heart interaction, and to generate the atrial and ventricular electrograms, we incorporated the myocyte model described in the previous section in the monodomain formulation of cardiac tissue: where D is the diffusion coefficient. We reported the values of D for the six regions of our cardiac model in Table 2. The values of diffusivity were selected to replicate time intervals corresponding to complete atrial activation [28], atrioventricular conduction [29], and complete ventricular activation [30]. VOLUME 10, 2022 To reproduce the auto-rhythmic property of the SA node and the consequent depolarization of the surrounding tissue, we introduced an additional current component (I stim ) to the second member of Eq. 5, in the SA region of the model, characterized by non-null values at the times of SA node spontaneous activation. We enforced the following boundary condition at the domain boundaries: where DB N is the imposed flux. DB N changes in space and time and it is non-zero at the electrode-atrium and electrodeendocardium boundaries when the stimulation current is on. In all the other boundaries DB N is always zero. We spatially discretized Equation 7 by employing a centered finite difference scheme with a resolution of 0.02 cm, whereas we performed time integration by applying the forward Euler method with a time resolution of 0.02 ms. Moreover, we employed the smoothed boundary method, described by Fenton et al. in [31] and Yu et al. in [32], to implicitly solve for the boundary condition of Equation 6 on the non trivial 2D geometry of Fig.1. According to the smoothed boundary method, the partial differential equation (PDE) incorporating the boundary condition is: where ψ is a scalar field that varies between 0 and 1 on a thin band spanning the external boundaries of the heart domain. Equation 7 can be defined on a square computational domain containing the cardiac domain, so that it can be easily represented with a standard finite difference approximation, thus avoiding complex meshing on the boundaries. Note that we updated the variables of the model only in the regions of the computational domain where ψ is greater than a certain threshold (i.e., on the cardiac domain and just outside the boundaries).
To depolarize the tissue surrounding the electrodes due to the external pacing events, the DB N parameter is set to a value equal to twice the diastolic threshold for of a duration of 1 ms, when the heart model receives the atrial or ventricular activation signal. In addition, we evaluate the atrial and ventricular electrograms (i.e., the outputs of our cardiac model) by calculating the potential at the center of the electrode interfaces, in the hypothesis of homogeneous volume conductor: where β indicates a dimensionless constant, and r defines the distance vector between the source and the center of the atrial and ventricular electrodes.
We implemented the heart model in Matlab (R2022a, MathWorks, Natick, Massachusetts) and we used the ''GPU coder'' toolbox of Matlab to optimize the speed performance. By using ''GPU coder'' we obtained an optimized CUDA code in the form of a Matlab MEX file. To test the speed performance of our heart model we executed a Matlab script containing the heart model MEX file on a workstation with the following specifications: CPU AMD Ryzen Threadripper 3960X 24-Core 3.79 GHz, graphics card NVIDIA GeForce RTX 3090.

B. CLOSED-LOOP SYSTEM
We developed a closed-loop system composed by our heart model and a DDD pacemaker model. We included the heart model in a Level-2 MATLAB S-function thus enabling to use it in Simulink. A level-2 MATLAB S-function is a Simulink block with multiple input and output ports which allows to include a MATLAB function in a Simulink model. It comprises a set of callback methods that the Simulink engine invokes when updating or simulating the model. We developed the pacemaker model in Simulink (R2022a, MathWorks, Natick, Massachusetts) by using the Stateflow toolbox of Matlab. According to previous works [33], [34], we set five different timers to define the DDD pacemaker operation: • AVI (Atrio-Ventricular Interval) represents the time period from an atrial event to a ventricular pace. It starts with an atrial activation (sensed or paced) and triggers ventricular pacing (VP) when the timer runs out if no ventricular event has been sensed (VS). It is used to synchronize the ventricular pacing with the atrial activity.
• LRI (Lowest Rate Interval) defines the longest interval between two consecutive ventricular activation (VP or VS). The LRI timer is reset at each ventricular event and delivers an atrial pacing (AP), after the time defined by Atrial Escape Interval (AEI), if no atrial event has been sensed. AEI is defined as a difference between LRI and AVI.
• URI (Upper Rate Interval) defines the maximum ventricular pacing rate.
• PVARP (Post-Ventricular Atrial Refractory period) starts at each ventricular event and serves to block unexpected atrial signals. If an atrial activation is detected during PVARP, it is marked as AR and does not impact the pacing schedule.
• VRP (Ventricular Refractory Period) defines a blocking interval for ventricular events. It is used to prevent the detection of unwanted ventricular activation. In each closed-loop simulation, we set the DDD pacemaker timing parameters with standard values: AVI = 150ms, LRI = 1000ms, URI = 400ms, PVARP = 185ms, VRP = 200ms [35]. The pacemaker model also replicates an algorithm aimed at terminating ELT, similar to the algorithms currently employed in cardiac pacemakers [14], [15]. During ELT, V-A conduction works as a retrograde reentrant circuit while AVI works as ''virtual'' A-V conduction pathway. Thus, ELT consists of a sequence of VP followed by atrial sensed events (AS), generally maintained at the maximum ventricular rate (URI). The algorithm detects ELT by recognizing 8 consecutive VP-AS cycles. In addition, ELT is confirmed only if the difference between each VP-AS interval and the first VP-AS interval detected is within ±32 ms. Thus, the tachycardia is suppressed by increasing the PVARP to a fixed value of 500 ms for the next cardiac cycle so that atrial activity induced by retrograde conduction is not sensed.
As described in Fig. 2, we developed an hardware setup for our closed-loop system (Fig.2). We deployed the Simulink model of the pacemaker on an Arduino Due board using the Matlab toolbox ''Simulink Support Package for Arduino Hardware'', whereas the heart model runs on the Workstation, as described in the previous section. The Workstation was also equipped with a data acquisition board (DAQ PCI-6023E, National Instruments), which allows the communication between the heart model and the pacemaker model. In particular, the pacemaker model sends atrial pacing (AP) and VP stimuli to the heart model by using two analog outputs of the Arduino board. In turn, the heart model sends to the Arduino board the atrial and ventricular activation signals (Ain, Vin) through two digital output channels. Ain and Vin are obtained by processing atrial and ventricular electrograms. The processing procedure consists of squaring, moving average filtering, and thresholding of the electrograms.
Lastly, we developed two GUI in Simulink to make the system usable by non-experts in cardiac modelling (see also Supplementary material). The first one, on the heart model side (cardiac GUI, Fig. S1), permits to regulate the heart rate and to simulate cardiac rhythm dysfunction (e.g., AV block, partial AV block, bradycardia) by using a dropdown menu. Additionally, the cardiac GUI enables to change the position of the electrodes of the pacemaker. By using two dedicated buttons of the GUI, it is also possible to generate additional SA node activation or PVC, elicited in a region in proximity to the ventricular electrode.
The second GUI (pacemaker GUI, Fig. S2) allows to set the timing parameters of the device (AVI, URI, LRI, PVARP, VRP) by using appropriate sliders in order to search for the best configuration for the cardiac situation analyzed. Fig. 3 shows the action potentials simulated at a pacing frequency of 1 Hz, with the parameter sets corresponding to the different cardiac regions. The myocyte model accurately reproduces the experimental action potential morphologies corresponding to the different cardiac regions. The simulated epicardial action potential shows a spike-and-dome morphology resembling the action potential measured by Nabauer et al. [18]. On the contrary, endocardial and midmyocardial action potentials do not possess a prominent notch, coherently with experimental recordings [18], [19]. The atrial model reproduces a spike-and-dome morphology comparable with experimental recordings [21], [22], [23]. However, other action potential patterns (e.g., triangular shaped) can be observed in the human atrium, depending on the relative intensity of delayed rectifier and transient outward currents [21]. Indeed, spatial heterogeneity of the atrial action potential has been observed in animal experimental models [36], [37]. Note that in this work we neglected atrial heterogeneity. Similarly, we did not consider ventricular apico-basal [38] and interventricular heterogeneity [39]. In our model, the maximum upstroke velocity in tissue is 236 V /s for the endocardium, 326 V /s for the midmyocardium, and 209 V /s for the epicardium, matching the experimental values of upstroke velocities reported in the literature for the three type of ventricular cells [16], [17]. For the atrial model the maximum upstroke velocity in a single cell is 181 V /s, which is comparable to experimental recordings in isolated atrial myocytes [40]. The resting membrane potential is −85 mV in the ventricular regions, and −75 mV in the atrial model, matching experimental recordings [16], [17], [18], [19], [21], [22]. Fig. 3 shows action potential duration (APD) restitution curves for the different cardiac regions, together with experimental data points from [19], [24], [25], [26], and [27] used also in the parameter fitting procedure. Our model formulation captures the APD rate dependence of all the four types of cells. The midmyocardium has the longest APD, up to approximately 590 ms at very large pacing cycle length (i.e., greater than 5 s). At short diastolic intervals (DIs), the APD restitution curve of midmyocardial cells approaches those of endocardial and epicardial cells. At long DIs, the APD in the endocardium is longer than in epicardium, whereas at short DIs the endocardial action potential becomes shorter than the epicardial AP. Furthermore, epicardial cells have a relatively flat APD restitution curve at long DIs compared to endocardial and midmyocardial cells. The atrial action potential is slightly longer than in the epicardium at rest, but it becomes shorter at short DIs approaching the endocardial restitution curve. The epicardial and atrial models also match the experimentally measured conduction velocity in epicardial (about 70 cm/s) [41] and atrial (about 55 cm/s) [42] tissue, respectively. The conduction velocity restitution curve for the epicardial model was already FIGURE 3. APD restitution curves for midmyocardial, endocardial, epicardial, and atrial cells compared to experimental data. Experimental data for midmyocardial cells are taken from [19], for endocardial cells from [24] and [25], for epicardial cells from [26], and for atrial cells from [27]. In the bottom of the figure the APs for the four types of cells are shown at a pacing frequency of 1 Hz. described and compared with experimental data in our previous work [9]. The atrial model reproduces the flat conduction velocity restitution curve experimentally observed by Feld et al. [42]. Fig. 4 shows the results of a 5 seconds simulation of the 2D whole heart model in open loop (see also Video S1) with the set of parameters describing a healthy activation sequence ( Table 2). The cardiac tissue is stimulated in the SA node at a heart rate of 1 Hz. The electrical activation in the sinoatrial node spreads out through the atria (Fig. 4A, top left) until it reaches the AV node approximately after 70 ms, as reported experimentally [28]. Complete atrial depolarization requires about 125 ms, coherently with experimental data (116 ± 18 ms) [28]. The depolarization of the atrial tissue generates two wide deflections in the atrial electrogram (Fig. 4B, top), whereas it does not affect the ventricular electrogram (Fig. 4B, bottom). The AV node introduces a significant delay in the AV conduction, slowing down the action potential propagation towards the cardiac septum (Fig. 4A, top right). Atrio-ventricular conduction time is about 130 ms, which is inside the range for healthy subjects (120-200 ms) [29]. The septum triggers the ventricular depolarization (Fig. 4A, bottom left), which generates a deflection in the ventricular electrogram (Fig. 4B, bottom). Complete ventricular depolarization requires about 100 ms, as reported experimentally [30]. In our model, the delay between atrial and ventricular activation depends on the diffusivity values, and mostly on the value set in the AV node. Thus, it is possible to modify such value to simulate different cardiac conditions such as partial or complete AV block. Note that the conduction in the endocardium is faster than in the other ventricular regions due to its higher diffusivity value ( Table 2). The atrial repolarization occurs almost simultaneously with the ventricular repolarization. Finally, the ventricular repolarization terminates in the midmyocardium, which indeed presents the longest APD (Fig. 4A, bottom right). The ventricular repolarization can be also appreciated in the ventricular electrogram as a small T-wave (Fig. 4B, bottom). In terms of simulation speed, we achieved a ratio of about 5, meaning that 1 second of cardiac simulation requires 5 seconds to run. Instead, for the code including the computational operations for the electrograms, we obtained a ratio approximately 10.

C. CLOSED-LOOP TESTS
We performed three different trials of the closed-loop system, corresponding to three pathological cases, to prove its functionality and potentiality.
Firstly, we simulated an AV block (Video S2) by setting to zero the diffusivity of the AV node (see Table 2), through an apposite drop-down menu of the cardiac GUI. Fig.5(a) shows five snapshots of the membrane potential map captured during the test. After an initial atrial depolarization, the action potential cannot propagate towards the ventricles due to the complete AV block. Coherently, when the AVI timer runs out, the pacemaker model delivers a ventricular stimulation (Fig.5(d)). In addition, Fig.5 shows both atrial (b) and ventricular (c) electrograms, in order to provide a complete overview about the test.
In second instance, we simulated a bradycardic condition (Video S3). From Fig.6(a) can be observed that after a first SA node activation, which leads to a normal ventricular depolarization, SA node fails to auto-depolarize again. As a result, consistently with the pacemaker expected behaviour, the pacemaker model delivers an atrial stimulation (Fig.6(d)) when the AEI runs out.
Lastly, we simulated a condition of partial AV block and slowing of conduction, in which a PVC leads to the onset of ELT (Video S4). Partial AV block and conduction slowing were simulated by reducing the diffusivity values according to Table 2. In Fig.7 are shown the atrial (a) and ventricular (b) electrograms and the pacemaker signals (c) obtained during the simulation. After three initial cycles in sinus rhythm, we generated an ectopic ventricular activation that propagates retrogradely towards the atria. Thus, after the time required for the V-A conduction, the device detects the atrial depolarization and marks the event as AS, which in turn triggers a ventricular stimulation, as imposed by AVI and URI. Thereby, VP generates a second V-A conduction, establishing ELT. The pacemaker paces the ventricle for every AS, increasing the ventricular rate inappropriately. Afterwards, following the recognition of eight consecutive VP-AS patterns combined with the aforementioned requirements about anti-ELT algorithm (reference in Sec.II-B), the device detects ELT and extends PVARP to 500ms to terminate the undesired tachycardia.

IV. DISCUSSION
Previous works modelled the interaction between heart and pacemaker in a closed loop system. Whereas the representation of the pacemaker algorithm as a timed automaton is a standard in the literature, our system differs from previous works in the modelling of the heart electrical activity. Previous works modelled the heart as an interconnected network of nodes that represent specific cardiac regions. In particular, two modelling choices have been adopted for the heart model. First, it has been represented as a TA [5], and later as a HA [7], [8]. A TA is a finite state machine whose transitions are defined by external inputs and a continuous timer. Thus, TA can capture cardiac timing properties and are also theoretically feasible for model checking [15]. A HA is a generalization of a TA, its internal states are defined by the values of a set of continuous variables that evolve in time, according to the ordinary differential equations that reproduce cardiac action potentials. Thus, while the behaviour of a HA is limited by the number of its internal states and transitions like a TA, it is able to capture more complex physiological dynamics [8]. In the work of Ai et al. [7], the heart model is composed of a network of HA. The spatial organization of the cardiac regions was taken into account as a delay in the nodal connections allowing the reciprocal interaction between neighbour nodes.
Differently from previous works, we used a PDE formulation to implement the heart model in the closed loop system. The PDE formulation is well established in the cardiac electrophysiological modelling literature and has several applications [43], [44]. To the best of our knowledge, our work is the first to introduce a PDE heart model to simulate the closed-loop interaction between CIEDs and the human heart. The use of PDEs to model the electrophysiological activity of the heart represents an important advancement in closedloop systems. First, in our heart model the cardiac tissue and AP propagation are continuous in space, avoiding the use of ''lumped'' nodes to represent the electrophysiological activity of specific cardiac regions. Notably, in HA models the representation of distributed and heterogeneous regions with ''lumped'' nodes introduces additional parameters in VOLUME 10, 2022  order to reproduce the electrophysiological behaviour of a continuous geometry, both in terms of electrotonic interaction and generation of cardiac signals. Indeed, each ''lumped'' node requires a specific parameter tuning procedure. On the contrary, in a PDE formulation the number of model parameters that needs to be estimated does not increase with the geometric dimension of the problem. In fact, once the myocyte models (and respective tissues) are defined, the use of a realistic geometry allows for the computation of realistic EGMs and electrophysiological behaviour without further adjustment.
Additionally, previous works based on HA formulation did not consider the role of transmural heterogeneity neither in the heart model nor in the genesis of electrograms. On the contrary, in a PDE formulation, the properties of the myocardial tissue follow from the definition of the myocyte model and tissue layers. Thus, it is possible to represent different tissue layers by adjusting the parameters of the myocyte model. It is worth mentioning that heterogeneity of electrophysiological properties in the myocardial layers contributes to the genesis of recorded signals [45] and is involved in the insurgence of arrhythmias for some cardiac pathologies [46], [47].
In this work, we considered midmyocardial cells to have a longer APD as reported in [16] and [19]. However, recent electrophysiological studies highlighted that there is no evidence for a coherent midmyocardial region with long action potential in the human ventricular wall [45], [48], [49], thus the role of M-cells in transmural heterogeneity of repolarization is still debated. Notably, a PDE formulation allows to model arrhythmic rotors (i.e., tachycardia and fibrillation), whereas a heart model made by lumped nodes is not able to capture such behaviour. Indeed, it is possible to represent a great variety of pathological conditions that contribute to the activation pattern of the heart and are relevant when testing the interaction between the heart and a CIED. As an example, cardiac pacemakers have been used to treat atrial fibrillation and atrial flutter patients [50], [51], [52]. Thus, to model such a condition, it is important to reproduce the pathological substrate of atrial fibrillation, which has been often associated with diffuse fibrosis in atrial myocardium [53], [54]. Note that HA and TA cannot model small spatial inhomogeneities, such as cardiac fibrosis.
Moreover, reaction-diffusion models can directly generate simulated cardiac signals, such as electrograms, which could be of great interest in pathological conditions, such as those mentioned above. For example, our framework could be employed to develop and validate new algorithms that reveal rhythm dysfunction based on specific features of cardiac electrograms. Similarly, our closed-loop environment could be exploited for the design and optimization of electrodes shape and position, whereas in TA and HA representations the electrodes are not physically included in the heart model.
Furthermore, in HA models, conduction velocity depends on the interaction between two neighbouring nodes; thus, the restitution properties of the conduction velocity do not depend only on the myocyte model but are determined by the interaction between nodes. Notably, both APD restitution and CV restitution can be easily optimized in a PDE model by using a fitting procedure directly applied to literature data, from in vivo or in vitro experiments [9], [23].
Recently, with the increase in technological availability and detail of diagnostic procedures, it has become possible to accurately represent the electrophysiological activity of a patient's heart. In particular, procedures to fit to a PDEs model patient-specific data, such as activation maps and cardiac geometries, were proposed. As an example, it is possible to estimate the distribution of the Purkinje network of a patient by using inverse problem methods [55], [56], [57]. Combined with a cardiac geometry from Cardiac Magnetic Resonance or Computed Axial Tomography, it is possible to represent the propagation of the depolarization waves directly on the patient's heart. This procedure has already been carried out in the literature with accurate results and is generalized to different pathological conditions, such as ischemic lesions or right bundle branch block [44], [58], [59], [60]. Indeed, visualizing the spread of action potentials on a heart geometry makes the simulations easily understandable to non expert users. Additionally, our GUIs allow users to interact with the models during simulations to reproduce pathological conditions and/or modify pacemaker behaviour.
The main drawback in the use of PDEs to model the electrophysiology of the human heart is the computational time and memory that are needed to solve the system of coupled equations. However, recent implementations of the cardiac equations have reached almost real-time speed by using proper GPU environments [61]. Indeed, the use of efficient numerical schemes and parallel computing [61], [62] can enhance the computational performance. In this work, we employed the finite difference method, which is typically faster than the finite element method [31]. Moreover, we adapted a previously published phenomenological model for cardiac cells, which is computationally efficient (i.e., 3 state variables) and easy to parameterize [9], [11]. Finally, we exploited parallel graphical computing that already proved to significantly accelerate cardiac simulation [61]. For the previous reasons, we believe that the choice of using a PDE representation of the heart model offers a greater degree of flexibility and adaptability of the framework when compared to a HA representation.
Our results demonstrate that by using our framework it is possible to replicate scenarios that arise exclusively when pacemakers interact in closed loop with a heart model, similar to what has been observed in previous articles employing HA models [8]. Furthermore, our system can be used to optimize the pacing algorithms and electrode placement [63], [64], [65]. Indeed, our closed-loop model can replicate different cardiac conditions and include different pacemaker algorithms. However, a more complex formulation of the heart model would be necessary in order to carry out more accurate simulations, which could also provide a platform for closed-loop model validation. In particular, a 3D, anatomically realistic geometry with a Purkinje network would grant a more realistic modelling of the cardiac depolarization and its interaction with the pacing activity.
In this work, we employed the monodomain formulation, neglecting current flows between the cardiac tissue and the surrounding passive conductor. Nevertheless, the heart tissue could also be represented by the bidomain formulation and considering fiber anisotropy. The bidomain formulation accounts for the changes in extracellular potential and can capture pacing artifacts like electrode crosstalk. However, the bidomain formulation requires solving an elliptic PDE at each time step, which is typically solved with iterative methods [43]. Thus, the time required for a simulation step is not known a priori, making hardware implementation difficult. Moreover, bidomain models are computationally much more complex than monodomain models. Computational efficiency is an important aspect in closed-loop models, particularly when hardware-in the-loop validation is of interest. Indeed, physical device validation requires real-time cardiac models. Our heart model cannot run in real-time yet, but we believe that further optimizations could fill the gap towards this goal.

V. CONCLUSION
In this manuscript, we described a PDE-based heart model specifically designed for closed-loop evaluation of heartpacemaker interaction. Moreover, we developed a specific hardware setup for the proposed closed-loop framework. The results showed that our approach can accurately resemble the interaction between the heart and pacemaker models, demonstrating the feasibility of reaction-diffusion heart model for closed-loop systems. Whereas the main limitation of the employment of reaction-diffusion models is the computational complexity, we showed that the gap towards realtime closed-loop simulations is affordable in near future. Besides the computational complexity, the main limitations of our closed-loop model are the isotropy of the myocardial tissue, 2D geometry, and lack of Purkinje system. Despite the aforementioned limitations, we believe that the proposed closed-loop system provides a novel and promising framework for the assessment of cardiac pacing. Indeed, our closed loop framework offers a graphical representation (see supplementary videos) of heart-pacemaker interaction, which makes simulation easily understandable also to non expert users (like physicians). Moreover, the user-friendly GUI of the proposed system would allow easy handling of the framework to explore different combinations of pacemaker parameters and dysfunctions of the heart rhythm. By using our framework, the user can verify the efficacy of the parameters chosen for the device in a controlled environment and check if the current set of parameters gives rise to hazardous situations when heart-pacemaker interaction occurs. In our opinion, the proposed framework offers the possibility to reduce the knowledge gap between the expert that programs the pacemaker and the physician that sets the requirements for its use on the patient. ALESSANDRO TOGNETTI (Member, IEEE) received the Graduate degree in electronic engineering and the Ph.D. degree in bioengineering from the University of Pisa. He is currently an Associate Professor in bioengineering at the Department of Information Engineering and the ''E. Piaggio'' Research Center, University of Pisa, where he is also the Vice President of the Degree Course in Biomedical Engineering. He holds the courses in biosensors, bioelectric phenomena and bionics senses with the University of Pisa. His research interests include the design and development of wearable multisensory technologies for measurement and analysis-outside the laboratory and during daily activities-of posture and human movement and physiological signals, development of soft sensors (flexible and e-textiles) for monitoring and analyzing human movement (upper and lower limbs with particular attention to both physical and neurological rehabilitation applications), non-invasive biomedical devices (physiological signals, physical activity, innovative interfaces for prostheses) and sensory fusion (fusion between physiological signals and inertial sensors for an accurate classification of human activity, and fusion of soft and inertial sensors to improve movement reconstruction performance by systems wearable). His main scientific skills range from the development of innovative biomedical sensors-starting from the physical principle and enabling technology-to system integration, from signal processing and information to the development of ICT applications for health monitoring, rehabilitation, robotics, and human-machine interaction. Much of the research was carried out in the context of national and international research projects (more than 20 projects of which 12 European) and through numerous scientific collaborations with universities, research institutes, and companies both nationally and internationally. His research activity has generated more than 100 international contributions. Open Access funding provided by 'Università di Pisa' within the CRUI CARE Agreement