Pacemaking function of two simplified cell models

Simplified nonlinear models of biological cells are widely used in computational electrophysiology. The models reproduce qualitatively many of the characteristics of various organs, such as the heart, brain, and intestine. In contrast to complex cellular ion-channel models, the simplified models usually contain a small number of variables and parameters, which facilitates nonlinear analysis and reduces computational load. In this paper, we consider pacemaking variants of the Aliev-Panfilov and Corrado two-variable excitable cell models. We conducted a numerical simulation study of these models and investigated the main nonlinear dynamic features of both isolated cells and 1D coupled pacemaker-excitable systems. Simulations of the 2D sinoatrial node and 3D intestine tissue as application examples of combined pacemaker-excitable systems demonstrated results similar to obtained previously. The uniform formulation for the conventional excitable cell models and proposed pacemaker models allows a convenient and easy implementation for the construction of personalized physiological models, inverse tissue modeling, and development of real-time simulation systems for various organs that contain both pacemaker and excitable cells.

Nowadays computer modeling of various organs and tissues is indispensable part of the 2 physiology research. Computational models of different levels of complexity are being 3 utilised, with complex ion-channel multi-variable models on the top of the list. A 4 rigorous review of the cardiac models was presented in [1]. The number of such models 5 and their updates increases yearly following new physiological findings and 6 measurements. Consisting of many differential equations for ion channels and their gate 7 formulations, the models provide detailed description of cell behavior under various 8 normal and pathological conditions, including, for example, influence of drugs [2]. The 9 number of variables and parameters in such models can reach several dozens, leading to 10 the necessity to utilize significant computational resources, even despite of current Timeline of simplified physiological cell models development. Solid shapes correspond to excitable cell models, dashed shapes -to pacemaking cell models, respectively.
(FHN), and Hodgkin-Huxley [7] (HH) equations. These models are based on small 20 number (usually 2 -3) of variables, i.e. ordinary differential equations (ODEs). Later 21 additions to this class include modifications of the HH and FHN models, namely Van 22 Capelle-Durrer (VCD) [8], Aliev-Panfilov [9] (AP), Morris-Lecar [10] and its 23 pacemaking variants [11,12], Fenton-Karma [13], Mitchell-Schaeffer [14] (MS), and its 24 modification by Corrado and Niederer [15] (CN). Most of the above mentioned models 25 are included into the modeling software packages and repositories, such as 26 openCARP [16] and Physiome Project [1,17]. 27 Fig 1 demonstrates timeline of the historical development of main simple 28 physiological models. Solid shapes correspond to the excitable models, dashed -to the 29 pacemaking models, and the models capable to exhibit both types of behavior are 30 represented by both shape types. Even though the majority of the models shown in 31 Fig 1 is capable to provide pacemaking operation, their utilization is limited. The VDP 32 and FHN models and their modifications are being predominantly used as simple models 33 of natural pacemakers in physiological simulations of different levels of complexity (see, 34 for example, [18][19][20][21][22]). The VDP model of relaxation oscillator was proposed to describe 35 general heartbeat dynamics and by its nature does not have the quiescent excitable 36 form. The two-variable FHN model, as a reduction of the four-variable HH model of the 37 squid giant axon action potential, was developed to model neuronal excitability. It has 38 the properties of producing spike trains (tonic spiking) at sufficiently large stimulating 39 constant current, and the apparent absence of a firing threshold [23], which are 40 undesirable for simulation of other organs, such as the heart. 41 The disadvantage of other simplified models simulating both pacemaker and 42 non-pacemaker action potentials of cardiac cells like VCD, Morris-Lecar, and their 43 modifications is the significant number of parameters that limits the areas of their 44 utilization. On the other hand, computationally lightweight models such as AP, MS, 45 and CN are being successfully used for solution of cardiac inverse problem and building 46 patient-specific models [24][25][26][27]. The MS model, however, is known to have some 47 stability problems [15]. 48 For multi-scale and real-time simulations the models with modest number of 49 parameters are preferable. Such models have been implemented in many studies of 50 electrophysiology of the heart and cardiac tissue [28][29][30][31][32] as well as intestine [18,33], 51 stomach [34], uterus [35], and bladder [36]. Recently, the AP model was used in deep 52 learning-based reduced order modeling [37], allowing to boost the solution of 53 parameterized problems in cardiac electrophysiology. 54 Another recently proposed type of pacemaker models developed to satisfy real-time 55 requirements is represented by parametric [38] and resonant [39] (8-24 variables) models. 56 They quantitatively reproduce action potential shapes and some cellular behavior, but 57 do not include several important physiological properties such as interactions due to 58 electrotonic coupling. Although the models provide relatively low computational load, 59 the latter may significantly rise with inclusion of the detailed features. The above 60 mentioned models, however, underline the need for simple computationally efficient 61 models, in particular, having unified description for pacemaking and excitable cells.

62
In this work, we consider variants of the AP and CN cellular models providing them 63 intrinsic pacemaker properties (hereinafter called pAP and pCN, respectively), and 64 demonstrate their main characteristics for both single pacemaking cells and coupled 65 pacemaker-excitable systems.

66
As application examples for the proposed pacemaker models, we include simulations 67 The two-variable AP model [9] proposed to describe non-oscillatory cardiac tissue that 72 supports stable propagation of excitation waves is represented by the following set of 73 reaction-diffusion type nonlinear ordinary differential equations: where u and v are normalized transmembrane potential and slow recovery variable, (HB), which is a typical mechanism for the onset of oscillations with a stable limit cycle 92 (periodic orbit) in nonlinear dynamical systems [40]. For the sake of convenience, 93 further we introduced a parameter b AP = −a in Eq (1), controlling the intrinsic 94 oscillation frequency when b AP > 0: The idea to use of the parameter b AP for cubic u-nullcline in the FHN-like systems 96 to control excitability range is straightforward. It has been used, for example, to study 97 the propagation of action potential in combined pacemaking-excitable FHN model 98 tissue [41]. However, in the case of the AP model, which was developed primarily to 99 represent excitable cardiac tissue, the intrinsic pacemaking function in a single cell and 100 coupled systems was not considered yet.

101
The resulting phase space geometry of the pAP model is shown in Fig 2A.  The two-variable CN modification [15] of the ionic MS model [14] for cardiac excitable cells is represented by the following set of nonlinear differential equations: Here h is the gate variable for the inward current (sodium ion channels), u gate > 0 is 106 the excitation threshold potential, τ in , τ out , τ open , and τ close are the time constants 107 affecting the corresponding characteristic phases of the evolution of transmembrane 108 potential u (shown after suprathreshold stimulation in Fig 2B.1). As in the pAP model, 109 the latter is also limited in the region u > 0 by the second u-nullcline u = 0.

110
To introduce pacemaking behavior into the CN model, we did the following 111 modifications. First, the piece-wise function (Eq 5) was replaced by the formulation 112 with sigmoid function [7,42] of the transmembrane potential (Eqs 7-9) with the slope 113 factor u s (in dimensionless voltage units), similar to the approach used in [43], changing 114 the shape of h-nullcline [44]. Second, we replaced u gate in Eq 4 by a parameter 115 b CN = −u gate , allowing to shift independently left branch of the u-nullcline Similar to the pAP model, b CN becomes a parameter suitable for control both the CN 118 cell excitability and the pCN intrinsic oscillation frequency (see  Though the method is similar to the previously proposed pacemaking modification 132 of the MS excitable cell model [43], the considered pCN model (Eqs 6-9) possesses 133 different nonlinear dynamic properties.  Because the load in a tissue is not limited to an integer value (n can be non-integer [46]), 152 apart from normal case (n = 1) we considered higher (n = 2) and lower (n = 0.5) loads. 153 Minimal and maximal frequencies of complete 1:1 synchronization   As an application example for the pAP-AP system we replicated simplified 2D SAN 175 simulation model used to illustrate the effect of SAN -atrium coupling on the 176 pacemaking behavior (Fig 11 in [47]). The model consists of a rectangular area 10 mm 177 by 10 mm (200 × 200 mesh, ∆x = 0.05 mm spatial step size) of atrial tissue represented 178 by AP model cells (Fig 3A). The whole 2D model is anisotropic with the ratio of  In the second application example for the proposed pacemaker models, we 195 considered a simple 3D elecrophysiological model of the small intestine. As a reference, 196 we took some of the results from the papers [18,33], in which the MS and FHN models 197 were used. In contrast to the works, we described ICC and SMC layers with the pCN 198 and CN model cells, respectively, to demonstrate the broad applicability of the pCN 199 model. It has been demonstrated that the excitable MS model is apt to spontaneous 200 excitations at some conditions, with its modification proved to be robust to such 201 pacemaker behavior [15]. Using the pCN model with improved frequency control in the 202 combination with robust CN model instead of the MS and FHN models may improve 203 the behavior of the pertinent computational tissue models.

204
Similar to [18,33], our intestine model geometry is presented by a long two-layer tube, 205 cut along its axle y on one side and stretched to a dual-layer plane (Fig 3B). The blue 206 and red surfaces in Fig 3B represent external SMC and internal ICC layers, respectively. 207 Such a simple tube geometry corresponds to anatomy of the small intestine in general. 208 The simulation domain for both layers has dimensions N x × N y = 252 × 2400, with 209 uniform spatial mesh size ∆x = 0.5 mm, which corresponds to the 1200 mm long tube 210 (one half of that in [18]) with the mean circumference of 126 mm (40 mm diameter).

211
Electrical dynamics in the intestine layers is described by the following set of ODEs September 10, 2021 6/16 for transmembrane potentials u I of ICC and u M for SMC: where I I tot and I M tot represent right-hand side of the Eq 6 for ICC and SMC layers, 212 respectively, and last terms in the right-hand side of Eq 10 and 11 describe electrotonic 213 coupling between the layers with diffusion coefficient D IM = 0.006 mm 2 ms -1 .
where i is the cell index counting from the duodenum.

224
To demonstrate appearance of intestinal dysrhythmias in a similar way as in [33], Equilibrium points for bifurcation diagrams were calculated with MatCont software [50]. 235 The parameters for the AP and pAP, CN and pCN models used in the simulations are 236 listed in Table 1 and Table 2, respectively. The initial conditions for the models were 237 chosen to be u(0) = 0.01, v(0) = 0.01 for the AP and pAP, and u(0) = 0.01, h(0) = 0.2 238 for the CN and pCN models.

239
The accuracy of the results for single cell simulations with fE method at different time steps ∆t were compared with that calculated with unconditionally stable implicit backward Euler (bE) method with small time step ∆t = 0.0001 ms. The relative norms (for a single cycle) and relative frequency error (for 60 s runs), as well as speedup of calculations are given in Table 3. For 1D coupled pacemaker-excitable systems the maximum relative frequency error with ∆t = 0.1 ms was also about 0.16%. As seen from Table 3, in most cases of the simulations with the pAP and pCN models the fE method with time step of 0.1 ms would be enough to obtain reasonable accuracy. This where N is the dimension of the simulation domain, should be taken into account as 240 well [51].  All the bifurcation parameters affect the intrinsic oscillation frequency (central 256 columns of Fig 4 and Fig 5). For pAP the highest variation in the frequency was 257 observed with b AP (with almost linear dependency, Fig 4A.2) and I ext (Fig 4E.2). The 258   (Fig 4D.2). For the pCN model significant variation of POP and MDP was observed for 263 all bifurcation parameters (left columns in Fig 5).  Fig 4C.3). For the pCN model this takes place only for 270 I ext dependence (Fig 5D.3).

271
In the pCN model b CN (Fig 5A.2) and u gate (Fig 5C.2) are the parameters with the 272 highest variation of intrinsic frequency. The parameter b CN , similar to b AP , allows to 273 obtain very low oscillation rates (less than 0.1 Hz). The influence of I ext (Fig 5D.2) on 274 the frequency is relatively weak, thus the pCN model looks less sensitive to the coupling 275 strength (compare with Fig 4E.2, see also Fig 6).  characteristics is similar to that for the original and modified MS models [14,15,43].

284
Dynamics of 1D coupled pacemaker-excitable system 285 Depending on the coupling (diffusion coefficient) the pacemaker-excitable system may 286 be either in a fully synchronized (1:1) regime, when all excitable cells in a strand follow 287 the driving frequency, or in a asynchronous/chaotic regime, when not all of the 288 pacemaker action potentials are able to propagate to the end of the strand [47,53].  For the pAP-AP coupled system, the complete synchronization was confined in two 295 separate areas created by interlocks of minimal and maximal frequency curves 296 (Fig 6A.1). The areas partially merged in the case of n = 0.5. Similar corresponding 297 separate areas of the parameter b AP are seen in Fig 6A.2. This indicates that at certain 298 values of d complete 1:1 synchronization in the pAP-AP system with heavy load can not 299  frequency in the pAP-AP system (Fig 7A), while in the pCN-CN system it reaches a 329 maximum before final fall-off (Fig 7B).  Though both the original AP and CN models were introduced to describe the cardiac 366 action potential, they have been also utilized for simulations of various organs which 367 contain both pacemaking and excitable cells. We intentionally did not try to tune 368 precisely the parameters of the pAP and AP cells to the existing experimental SAN and 369 atrial tissue data (including the shape of action potentials), with the aim just to 370 demonstrate the applicability of the pAP model for simulations of complex physiological 371 configurations. Additional tuning of the parameters might be necessary for accurate 372 simulation of organs other than the heart, as seen from the case of the 3D intestine 373 model.

374
The utilization of the considered pacemaker models can be beneficial for the 375 development of simplified real-time simulation systems (e.g. personalized 376 medicine [24,25] and inverse cardiac modeling [26]), deep learning [37], as well as for 377 medical device testing platforms [38,52].   One of the important features of the considered pacemaker models is that they do 389 not include neither any additional expression for currents nor additional ODE, thus 390 having the same number of variables as the original excitable models. We believe that 391 our proposed models can be effectively used for computationally efficient 392 electrophysiological modeling of tissues that include primary and subsidiary pacemaking 393 cells, allowing to build different models of whole organs with simple uniform model 394 description, as well as for patient-specific medicine and medical device testing requiring 395 real-time simulations.