Dynamics Study of Constant DI and Constant TR Control for Cardiac Alternans Based on a Two Dimensional Cellular Automata Model

Min Xiong University of Tennessee Knoxville College of Engineering: The University of Tennessee Knoxville Tickle College of Engineering Kai Sun University of Tennessee Knoxville College of Engineering: The University of Tennessee Knoxville Tickle College of Engineering Xiaowen Su University of Tennessee Knoxville College of Engineering: The University of Tennessee Knoxville Tickle College of Engineering Elena G. Tolkacheva University of Minnesota Twin Cities Campus: University of Minnesota Twin Cities Xiaopeng Zhao (  xzhao9@utk.edu ) University of Tennessee Knoxville College of Engineering: The University of Tennessee Knoxville Tickle College of Engineering https://orcid.org/0000-0003-1207-5379


INTRODUCTION
Cardiac alternans is a disturbed heart rhythm in which the action potential duration (APD) of cells alternates between long and short values, and manifests on electrocardiograms (ECG) as alternantion in amplitude between two consecutive T-waves [1] . In experimental and clinical practice, cardiac alternans could be a sensitive indicator of arrhthmias leading to life-threatening sudden cardiac death (SCD), which accounts for the death of around 325,000 Americans per year. Thus, understanding and preventing cardiac alternans has been a research topic of great significance [2] .
Based on the restitution theory, which has been widely used in cardiac alternans research, the APD of a cell in one cycle is a function of its diastolic interval (DI) in previous cycle, as shown in Eq. (1), where DI is the difference between the basic cycle length (BCL) and APD in the same cycle. 1 () n n APD f DI   (1) where f is the restitution curve.
Then, as shown in Fig. 1(a), under a constant BCL or so-called periodical pacing protocol, when the BCL is long (e.g., BCL=BCL1), the cardiac response is normal and APD is locked as 1:1 rhythm, and at steady state 11 1 1 1 ( ) APD f DI BCL APD DI   (2) Also, as shown in Fig. 1(b) and Fig. 1(c), APD alternates between APD2 and APD3 when BCL=BCL2 is small and then we have: Besides, Eq. (1) indicates that alternans occurs when slope of the restitution curve is smaller than 1 [3], [4] . As shown in Fig. 1(c) and Fig. 1(d), APD alternates when BCL is smaller than BCL BCL ONSET , and stable normal response occurs when BCL is larger than BCL BCL ONSET , so BCL BCL ONSET is a bifurcation point.
To suppress cardiac alternans, different control methods have been proposed and tested in the past decades using both numerical models and/or experiments on isolated animal hearts. An adaptive proportional feedback control approach was proposed and tested on rabbit hearts in [4].
In [5], an alternating-period feedback control mechanism was tested through numerical simulation on a single cell and a 1-dimensional (1-D) tissue. Besides, the effect of a feedback control method was studied on isolated canine Purkinje fibers in [6]. However, feedback control methods are generally ineffective for tissue with a distance longer than 2 cm due to complex spatiotemporal electrophysiological dynamics, and stimulation from multiple sites are required for large tissues, also, it is difficult to determine optimal control parameters [7]- [9] . Due to the limited success of the feedback control, which utilises a beat-to-beat regulation of the BCL, a constant DI control method that decouples the relation between APD and DI was adopted to better suppress alternans [10] . As shown in Fig. 2, constant DI pacing protocol keeps the DI of the cell constant based on measured transmembrane voltage of the 1 st cell in a tissue. In this way, the instability caused by change of APD can be eliminated, thus the pacing rhythm could be easier stabilized. Constant DI control was employed and tested in single-cell models [10], [11] and ventricular fibers models [12], [13] based on numerical simulation and/or experiments. However, these studies have validated efficacy of the constant DI control only on 0-D or 1-D tissues, and demonstrated that control may fail on 1-D cables and 2-D tissues with large size due to complex spatiotemporal evolution of cardiac activity [14], [15] . In addition, it is challenging to apply constant DI control in experiments due to the difficulty of performing real-time APD and DI measurements [15] . So, it is important to develop other control approaches with higher efficacy.
As is known, the APD alternans can be observed on ECG as T-wave alternans which is easier to be detected in real-time [16], [17] , and the TR interval can be regarded as a global DI. Thus, similar to the concept of constant DI control, a constant TR control method was proposed and employed in recent research [15], [17], [18] . As shown in Fig. 3, under constant TR pacing, based on the ECG of the whole tissue, stimulations are added after detection of T peaks to keep TR intrvels on the ECG constant, and therefore the constant TR protocol can be regarded as a "globalized" version of the constant DI protocol. In [15], the antiarrhythmic effects of constant TR control were validated using experiments on ex vivo whole rabbit hearts, yet there are some limitations such as accurate detection of a distinct T-wave and size of the heart. In [17], based on a 1-D numerical model of human ventricular tissue (<7 cm), the efficacy of constant TR control of alternans was studied and compared with that of constant DI control. The results show that the constant TR control is superior for alternans control in longer fibers when compared to constant DI pacing.
Nevertheless, the simulation efficiency was low, and the study has not been performed on a higher dimensional (e.g., 2-D) model with a larger size.
Thus, it is of great significance to efficiently study cardiac alternans control on a high dimensional model with large sizes, and a cellular automata model (CAM) is a good choice to achieve this goal. As a powerful model, which has successful applications in biology, physics, and other engineering areas, a CAM uses simple rules to compute complex phenomena with fewer calculations, and the process can be easily visualized to provide an intuitive illustration of system dynamics [19] . Although CAMs are less accurate than detailed models and vivo animal hearts experiments, they can represent the main features and key dynamical behaviors of the cardiac tissue. Compared with regular full physiological models, which require extensive computation, CAMs are much more efficient. Also, compared with ex vivo hearts, CAMs are more convenient and economical for tests on tissues of different sizes.
In this paper, we tested the efficacy of both the constant DI and constant TR control on alternans prevention using a 2-D CAM.

Two-Dimensional CAM for Cardiac Simulation
Both 2-D and 3-D CAMs have been developed to simulate cardiac electrical wave propagation and fibrillation [20]- [26] ; however, dynamical phenomena related to cardiac alternans have never been studied in these publications. Recently, a 2-D CAM capable of simulating cardiac dynamics including ectopic beat, wave break, spiral wave, conduction block, and alternans was proposed by our group [27] . Here, we adopted the same model with better practical parameter settings for the study of cardiac alternans control. Detailed mechanism of the employed 2-D CAM is described below.

Wave Propagation Rule in the CAM
To simulate dynamics of a 2-D cardiac tissue, an n×n CAM is constructed by MATLAB, and each small square in the CAM represents an individual cell of the tissue. Then, based on their transmembrane voltages, the cell states can be divided into four categories as listed in Table I [20] , where the transmembrane voltage of a cell becomes 1 when it is excited or stimulated.
A 3×3 square, which contains 9 cells in the bottom left corner of the CAM, acts as the paced cells, and they excite according to some protocols such as periodic pacing. Besides, as shown in Fig. 4, a single cell (black square) is excited at next time step, if its transmembrane voltage is less than 0.1 V and the voltages of at least three adjacent cells (gray squares) are bigger than 0.9 V [26] . Otherwise, its transmembrane voltage just evolves with time according to the voltage waveform function, which will be introduced below.

Transmembrane Voltage Waveform Function
The transmembrane voltage of a cell will suddenly change to 1.0 V the moment it is excited, otherwise, its voltage will progress according to the function below: where t is the time elapsed since last stimulation; c is a constant number which ensures V(APD,t=0)≈1.0 V and is set as 0.001 in this model; APD is determined by a restitution curve which will be given in next section. Thus, when the APD is determined, the parameter T can be calculated by Eq. (5), and then transmembrane voltage V in Eq.(4) will decrease from 1 to nearly 0 with the increase of t.
In addition, it can be easily proved or calculated that V(APD,t=APD)=0.1 V, which is exactly the refractory threshold. This also helps to define that the APD in this CAM is the time duration for the voltage to decrease from 1 V to refractory threshold 0.1 V.

Restitution Curve
To provide the value of APD for calculation of transmembrane voltage in Eq. (4), the restitution curve used in this paper is given below according to cardiac restitution theory: where the subscript n+1 and n denote the (n+1) th and n th stimulation cycle, respectively.
From the restitution curve, it can be observed that APDn+1 increases with DIn in a nonlinear way.
With the CAM introduced above, cardiac scenarios such as periodic normal conduction, normal conduction considering scar tissue, spiral wave, and especially alternans can be simulated efficiently in MATLAB, and the case illustrations can be found in [27].

Pacing Protocols
Downsweep protocol is a widely used pacing method to study cardiac dynamics. Tpically, pacing of the tissue starts from a long BCL=BCL0 (around 1000 ms), and then BCL is progressively decreased (equivalent BCL is used for constant DI and constant TR pacing). For each specific BCL value, a certain number of stimulations (typically 50) are applied to reach steady state [28] [29] .

Constant BCL downsweep pacing protocol
Under constant BCL downsweep pacing, the paced cells are stimulated each time after a constant period. Based on the simulation model introduced above, the values of BCL are set as a series of numbers decreasing from 350 ms to 300 ms as shown in Eq. (7).
For each specific constant BCL value in the series, 50 stimulations are applied to ensure the wave reaches steady state.

Constant DI downsweep pacing protocol
Unlike the constant BCL downsweep pacing protocol, under a constant DI downsweep pacing protocol, the paced cells are stimulated each time after a constant period (i.e., DI) from the instant when their transmembrane voltage drops to the refractory threshold 0.1 V, which is the threshold to distinguish between APD and DI intervals.
In this paper, to test the efficacy of constant DI control on the CAM, the values of DI in control are set as a series of numbers decreasing from 80 ms to 1 ms as shown below: For each specific constant DI value in the series, 50 stimulations are applied to test the control efficacy [10] [15] .

Constant TR downsweep pacing protocol
Under a constant TR pacing protocol, the paced cells are stimulated each time after a constant period (i.e., TR) from the instant when a T-wave peak on ECG is detected. A peak detection algorithm from [30] is applied in our work to conduct online T-wave peak detection.
The values of TR in control are set as a series of numbers decreasing from 145 ms to 1 ms, as shown below: For each specific constant TR value in the series, 50 stimulations are applied to test the efficacy of the control [15] .
To conduct constant TR pacing, a ECG of the tissue should be recorded. In our CAM-based tissue, with electrode A and electrode B placed out of the tissue (see Fig. 3), the pseudo-ECG can be calculated as: where  denotes gradient operation; (xA, yA) is the location of electrode A; (xB, yB) is the location of electrode A; and the transmembrane potential Φe can be obtained by: where the distances rA and rB in Eq. (11) are provided by: where (x, y) is the location of any cell in the tissue.
Pseudo-ECG takes the action potential of all cells in the tissue into consideration; thus, it is a global representation of cardiac dynamics [31] .

SIMULATION RESULTS ON A 50×50 CAM
This section presents the simulation results under different pacing protocols using a 50×50 CAM. The central cell located at (25,25) is selected as a representative example to illustrate the pacing stability, and the straight diagonal cable from the paced cell at (0,0) to the corner cell at (50,50) is chosen to show the spatial dynamics. The simulation time step step is 0.1 ms, and the threshold for alternans is set as 1 ms.

Cardiac dynamics under constant BCL downsweep pacing protocol
Firstly, under constant BCL pacing protocol, the bifurcation diagram of the central cell (25,25) under different BCL values and the spatial alternans along the diagonal cells at a constant BCL=320 ms are shown in Fig. 7.
It can be observed from Fig. 7(a) that alternans occurs at BCL BCL ONSET =335 ms, which is the same as the value calculated from the restitution curve in Eq. (6), and conduction block happens at BCL BCL CB =311 ms. Spatial concordant alternans (SCA) occurs along the diagonal cable as shown in Fig.  7(b) when BCL=320 ms, which also indicates SCA in the 2-D tissue.
Illustrations of normal response occurs at BCL=340 ms and alternans occurs at BCL=320 ms generated by a 200×200 CAM are shown in Fig. 8. From Fig.8 one can see that the waves have the same length (yellow color) under normal response but oscillate between small and large values under alternans.
Thus, the CAM provides an intuitive visual illustration of cardiac dynamics. Also, transmembrane voltage of the cell (25,25) under constant BCL=340 ms and BCL=320 ms are shown in Fig. 9. From Fig. 9(a), we can also see that when BCL=340 ms, the APD remains constant, thus it is normal response. However, the APD under BCL=320 ms shown in Fig. 9(b) alternates.

Cardiac dynamics under constant DI downsweep pacing protocol
Secondly, under the constant DI pacing protocol, the results of the same central cell and diagonal cable as that in the constant BCL test are shown in Fig. 10. Here, the equivalent BCL is calculated by fresituttion(DI)+DI. Also, because BCL= fresituttion(1)+1≈117 ms when DI=1 ms, the limit of BCL under constant DI control pacing is given as BCL DI LIMIT =117 ms.
From the results in Fig. 10, there is no bifurcation, and the wave propagation can always be stable under constant DI protocol even when DI is 1 ms with a corresponding equivalent BCL less than 120 ms.
From Fig. 11, it can be observed that no matter under relative large or small equivalent BCL, e.g. under DI=60 ms or DI=20 ms, pacing of the cell (25,25) can always be stablilized.
Thus, simulation results validate the high efficacy of constant DI pacing on the 2-D CAM-based tissue.

Cardiac dynamics under constant TR downsweep pacing protocol
Thirdly, results under the constant TR pacing protocol are shown in Fig. 12. In the simulation, two electrodes are placed at (0,0) and (51,51) to generate the ECG.
From Fig. 12 we can see that the pacing is stable under constant TR pacing when TR is bigger than 26 ms with an equivalent BCL bigger than 117 ms. Also, conduction block occers at BCL TR CB =117 ms, and no alternans is observed in the whole process.
To illustrate detailed dynamics of the tissue under constant TR control, ECGs of the tissue and transmembrane voltage of the cell (25,25) under TR=118 ms (equivalent BCL=313 ms) and TR=61 ms (equivalent BCL=200 ms) are shown in Fig. 13 and Fig. 14, respectively.

DIFFERENT MODEL SIZE
In the last section, the efficacy of both constant DI and constant TR approaches were validated on a 50×50 model using a downsweep protocol where the DI and TR decrease by 1 ms each step.
The results show that both constant DI and constant TR control approaches could help to significantly decrease the equivalent BCL at which unstable dynamics occurs. More specifically, both approaches are fully effective in controlling alternans until an equivalent BCL is smaller than BCL DI LIMIT , which is the smallest BCL that can be given by restitution curve. However, the results in [17] show that both constant DI and constant TR pacing protocols have limited effect on alternans preventation when a 5 ms decrease step is adopted in the downsweep protocol, and both alternans and conduction block will occur.
Thus, to fully study efficacy of the two control methods, in this section, dynamic responses of the system under the three pacing protocols, e.g. constant BCL, constant DI, and constant TR, using downsweep protocols with different decrease steps are studied along with a consideration of model size.
Performance of the three pacing protocols was tested and results are shown in Fig. 15-17.
From simulation results, system dynamics remain the same when the model size is bigger than 20×20, thus only performance on tissues with a size less than 20 is given in BCL ONSET comparison. In addition, the BCL CB does not change with model size and is plotted in the histogram.
From the results, it can be observed that: (1) For constant BCL pacing protocol, the simulation results, including BCL ONSET and BCL CB remain the same independend of the model size and Δ.
(2) Both constant DI and constant TR pacing protocols are effective to prevent alternans compared with constant BCL pacing, and both protocols stablize cardiac response till the minimum equivalent BCL, e.g. BCL DI LIMIT , for Δ=1 ms. However, their efficacy decays with the increase of Δ in the downsweep protocol.
(3) For Δ>1 ms, i.e. 5 ms and 10 ms, the constant TR control is generally more effective than constant DI control.
In comparison, paper [17] used a fixed decrease step Δ=5 ms and presented that both constant DI and constant TR will cause alternans and conduction block when values of DI and TR are small. This is the same as what we observed in the CAM when decrease step Δ is set as 5 ms or 10 ms in the downsweep protocol. However, the effect of different decrease steps on control efficacy has not been studied and is first investigated in this paper. As observed from the numerical simulation, both constant DI and constant TR approaches are fully effective if the decrease step Δ in downsweep protocols are small enough.

CONCLUSION AND DISCUSSION
Based on a 2-D CAM, the efficacy of constant DI and constant TR pacing protocols on alternans preventation or elimination are studied in this paper by numerical simulation. Unlike the existing simulation work focusing on 0-D or 1-D tissue, the work in this paper is done on a 2-D model.
Although less accurate than real heart experiments, CAM provides an economical and efficient approach for cardiac dynamic research and establishes a basis for both higher dimensional cardiac simulation and alternans prevention in heart tissue.
Simulation experiments in this paper validate that: (1) CAM provides high efficiency and intuitive visual illustration for cardiac electrical wave propagation simulation.
(2) Both constant DI and constant TR protocols are effective to resolve cardiac alternans.
However, their efficacy depends on the adopted decrease step in downsweep protocols, and a small decrease step could improve their performance.
In the meanwhile, there are some limitations of this study: (1) The employed CAM uses simple rules to simulate wave propagation, and the accuracy is not as high as the traditional full physiological model. The rules can be improved to keep a better balance between efficiency and accuracy.
(2) The intracellular calcium cycling which is important for alternans development is not considered.
Future work can be done on the following aspects: (1) Using a 3-D CAM to simulate caidiac wave propagation and study the efficacy of different control approaches.
(2) Improving the rules of CAM so that more details can be included without a significant sacrifice of computation efficiency.
(3) Study the essential mechanism in the dependence of control efficacy and decrease step in downsweep protocols.
(4) Conduct tests on well-developed two-dimensional full physiological simulation models or ex vivo whole animal hearts to validate the results obtained from the CAM-based simulation in this paper.

CONFLICT OF INTEREST
The authors declare that they have no conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. ONSET , and when the slope is bigger than 1, APD+DI>BCL BCL ONSET and the pacing is normal, also, when slope is smaller than 1, then APD+DI<BCL BCL ONSET and alternans occurs.