A cellular automaton model for the ventricular myocardium considering the layer structure
Deng Min-Yi†, Dai Jing-Yu, Zhang Xue-Liang
College of Physical Science and Technology, Guangxi Normal University, Guilin 541004, China

Corresponding author. E-mail: dengminyi@mailbox.gxnu.edu.cn

*Project supported by the National Natural Science Foundation of China (Grant Nos. 11365003 and 11165004).

Abstract

A cellular automaton model for the ventricular myocardium considering the layer structure has been established. The three types of cells in this model differ principally in the repolarization characteristics. For the normal travelling waves in this model, the computer simulation results show the R, S, and T waves and they are qualitatively in agreement with the standard electrocardiograph. Phenomena such as the potential decline of point J and segment ST and the rise of the potential line after the T wave appear when the ischemia occurs in the endocardium. The spiral wave has also been simulated, and the corresponding potential has a lower amplitude, higher frequency, and wider R wave, which accords with the distinguishing feature of the clinical electrocardiograph. Mechanisms underlying the above phenomena are analyzed briefly.

PACS: 05.45.–a; 87.18.Hf
Keyword: cellular automaton; electrocardiograph; ischemia; spiral wave
1. Introduction

The ventricle is one of the most important parts of the heart. The ventricle beats with high speed when ventricular fibrillation happens, which may result in sudden death, [1] so the dynamics of the transmembrane potential of the ventricle has received much attention.[25] The QRS complex and T wave of electrocardiograph (ECG) are the collective reflection of the ventricular electronic activity (Fig.  1).[6] The depolarization and repolarization of the ventricle are responsible for the QRS complex and T wave, respectively. On the temporal morphology of the standard ECG, the T wave is the first up wave after the QRS complex, where the potential is positive.[6] Myocardial diseases can usually be diagnosed with ECG; especially, the arrhythmia cannot be diagnosed without ECG. So the relation between the ECG and the electric activity of the myocardium cells should be made clear. Experiment is the main method used to study the electric activity of the myocardium cells on the macroscopic level. For example, based on a quasi-two-dimensional layer of the left ventricular epicardium of a white rabbit, Veniamin et al. studied the transition between the damped waves and the steadily propagating waves in the cardiac tissue.[2] Using the data of a patient who suffered from regular and sudden-onset fast palpitation, Sorgente et al. researched the contribution of the T wave vector in localizing the site of origin of the ventricular tachycardia.[4] But until now, no experimental reports about the relation between the myocardial electric activation and the ECG on the cellular level can be found. The reason for this may be the huge number of myocardial cells. The numerical simulation models play an important assistant role in exploring this relation. The numerical simulation models of myocardial electrical activation are classified into two major groups: one is the reaction– diffusion model, and the other is the cellular automaton (CA) model. Studies based on the former have obtained remarkable results about the dynamics of the spirals in excitable media.[710] In this paper, a CA model is established to recovery the relation between the cellular electric activation and the ECG.

Fig.  1. The standard ECG.[6]

CA is a simple and effective numerical method.[11] The CA approach usually oversimplifies the real system, but as long as the rules keep the basic properties of the system, we can still learn many interesting facts and gain insights into the mechanism underlying the phenomena.[12] Because of the easy realization of the program and the high speed of computation, CA has been used in the study of cardiac tissue for decades. For example, CA models for atrial tissue, [13] ventricular tissue, [3] and the sinoatrial node[14] have been established. Some researchers have also studied the electric signal of cardiac cells with CA. For example, Mitchell et al. performed a theoretical analysis of ventricular fibrillation and the requirements for fibrillation.[15] Zhu et al. studied the electric activity of the heart with a CA system.[16, 17] However, the CA model in Ref.  [3] did not consider the layer structure of the real ventricle, and the CA system in Refs.  [16] and [17] contained the complex ionic currents, which violates the principle of simplification of CA. Consequently, t is necessary to establish a new CA model for the ventricle tissue.

In this paper, a CA model for the ventricle tissue considering the three-layer structure is reported. There are three types of cells in this CA model. The computer simulation is carried out. The results show the R, S, and T waves when the system evolves with the electric travelling waves under the normal conditions. In the case of an ischemia happening in the endocardium, the point J and the segment ST decline, while the potential line after the T wave rises. The spiral-wave state of the electric impulse and the corresponding electrogram are also simulated. In the following sections, the model, the simulation results, and the conclusion will be presented.

2. The CA model for the ventricle

The experimental results have demonstrated three types of cells in the ventricular wall: endocardial, subepicardial, and epicardial cells.[18] These three types of cells have different repolarization characteristics: the action potential duration (APD) of the subepicardial cells is the longest; however, the deference between the epicardial cells and the endocardial cells is not worthy of mention.[19] On the other hand, the endocardium is the thickest, then the subepicardium, and the epicardium is the thinnest.[18] Considering the above factors, we construct the CA model for the ventricle as follows: the ventricle slice is put on the xoy coordinate surface, where the y axis direction is from endocarium to epicardium (Fig.  2(a)). This two-dimensional ventricle slice is divided into 128 × 32 sites. The ventricular cells are distributed on these sites. According to the thicknesses of the cells, the endocardium, subepicardium, and epicardium occupy 25 sites, 5 sites, and 2 sites, respectively. The transmembrane potential of a cell can be in m states. The state of cell (x, y) at discrete time t is represented by ux, y (t) ∈ (0, 1, 2, … , m − 1) with m ≥ 3. The ux, y (t) = 0, 1 denotes the rest and the excited states, respectively, while ux, y (t) = 2, … , m − 1 represents the refractory states. Obviously, the APD of a cell equals m − 1. Referring to Ref.  [19], the different types of cell have different m; for example, m in this paper is 14, 15, and 12 for endocardial, subepicardial, and epicardial cells, respectively (Figs.  2(b)– 2(d)).

Fig.  2. CA model for a ventricle slice. (a) Ventricle slice containing three layers. Transmembrane potentials of (b) endocardial cells, (c) subepicardial cells, (d) epicardial cells. Symbols ∘ , • , and ★ represent the rest, the excited, and the refractory states, respectively.

The no-flux boundary condition and the Moore neighborhood are used throughout this paper. For radius r, there are (2r + 1)2 − 1 neighbors in a cell’ s neighborhood. If the number of excited neighbors of cell (x, y) at time t is n, then the time evolution of the cell is given as follows:

where K is the excitable threshold.

3. Simulation results

When the computer simulation is done, the potential at the field point is calculated using the formula proposed by Zhu et al.[17]

where σ in = 1.28  S/m is the intracellular conductivity; [20]S0 is a constant that just affects the baseliner, and here S0 = 32; N is the total number of the cells; there are eight nearest neighbors of each cell; Vi is the transmembrane potential of the i-th cell in the system; Vj is the transmembrane potential of the j-th neighbor of the i-th cell; θ j is the included angle between two directions: one points from the i-th cell to its j-th neighbor, the other points from the j-th neighbor of the i-th cell to the field point; σ i is the average conductivity of the cells along the direction from the i-th cell to the field point; and, rj is the distance from the j-th neighbor of the i-th cell to the field point. For simplicity, we let the conductivity of the other medium besides the cardiac cells be 0.74  S/m.[20]

3.1. Travelling wave and the corresponding electrogram

The electric impulse is produced at the sino-atrial node, and then transmits to the atria, the atrioventricular node, and finally to the ventricle.[21] Because of the properties of the conduction system in the ventricle, the impulse can spread throughout the ventricle in a short time.[21] To imitate this physiological characteristic, all of the endocardial cells are set at the excited state initially (t = 0), while the other cells are set at the rest state. Let r = 5, which makes sure that the subepicardium can be excited in one step, and then the epicardiaum. After three steps, the system returns to the refractory states. According to the difference of APD between these three types of cells, the epicardial cells first come back to the refractory states, and then the endocardial cells, and finally the subepicardial cells. At last, the system returns to the rest state. The system will repeat the above state transformation again as long as the endocardial cells accept other impulses. The excitable threshold is K = 2, the period of impulse is Te = 33, and the field point is (130, 40), which is close to the V4 lead of the electrocardiographic leads.[6] The simulated potential trajectory is shown in Fig.  3(a). To analyze the mechanism underlying the trajectory in Fig.  3(a), the simultaneously recorded transmembrane potentials from endocardial, subepicardial, and epicardial cells are shown in Fig.  3(b).

Fig.  3. (a) Electrogram recorded at (130, 40) and (b) the simultaneously recorded transmembrane potentials from endocardium, subepicardium, and epicardium cells, where K = 2, r = 5, and Te = 33.

By comparing Fig.  3(a) with Fig.  1, we find that the potential trajectory shows the typical R, S, and T waves when the electric impulse evolves as a travelling wave in the ventricle slice. The occurrence mechanism of these waves can be found in Fig.  3(b). According to Eq.  (1), the potential comes from the differences of the transmembrane potentials; in our CA model, it comes from the differences of the transmembrane potentials of the three types of cells. For example, at the first step (t = 0), only the endocardial cells are excited, so the transmembrane potential of the endocardial cells is higher than that of the subepicardial cells, which contributes a positive potential to the field point (see the point at t = 0 in Fig.  3(a)). The other points in Fig.  3(a) can be analyzed similarly. It is necessary to point out the significance of the T wave. According to the definition of the T wave, one can find in Fig.  3(a) that the first T wave starts at t = 7 and ends at t = 15. In Fig.  3(b), the potential differences between the three layers hardly exist between these two moments without the existence of the subepicardium. So the existence of the subepicardium is the foundation of the T wave, which is a well-known physiological conclusion.[19] The ventricle slice considered here does not include the ventricular septum, which is responsible for the Q wave according to the physiological version, [6] so the Q wave does not appear here.

3.2. Ischemia and the corresponding electrogram

The segment ST is an iso-electric level which starts at the end of the S wave and ends before the start of the T wave; [6] for example, the segment in Fig.  3(a) starts at t = 3 and ends at t = 7. The point J is the first point of the segment ST, [6] the point is at t = 3 in Fig.  3(a). The potentials of segment ST and point J will decline when the ischemia occurs in endocardium.[6] To simulate this phenomenon, we suppose that the region within x = 120 to x = 128 and y = 0 to y = 25 is ischemic. The cells in the ischemic region act with higher resting potential and shorter APD[22] compared with that of the normal cells (Fig.  4(a)). The simulated potential trajectory is shown in Fig.  4(b).

Fig.  4. (a) Simultaneously recorded transmembrane potentials from endocardium, subepicardium, and epicardium; and (b) electrogram recorded at (130, 40). Parameters are the same as those in Fig.  3 except the existence of the ischemia in endocardium.

The results in Fig.  4(b) show that: (i) the potentials of point J and segment ST really decline. For example, the potential at t = 3 (the first point J) is no longer 0 in Fig.  3(a) but − 0.0036, and the potential reduces continuously afterwards until it reaches the minimum − 0.1035 at t = 8. Because of the shorter APD, the transmembrane potential of the ischemia endocardial cells is much smaller than that under the normal conditions between these two moments, which results in the even more negative potential for the field point (130, 40). (ii) The potential line after the T wave raises; for example, the potential between t = 15 and t = 32 in Fig.  4(b) is 0.0722, but it is 0 in Fig.  3(a). These results can be explained with Fig.  4(a): between t = 15 and t = 32, because of the higher resting potential, the transmembrane potential of the ischemia endocardial cells is no longer equal to but larger than that of the subepicardial cells, which contributes a positive potential to the field point (130, 40), and then results in the rise of that potential line.

3.3. Spiral wave state and the corresponding electrogram

Nasha et al. indicated that many dangerous cardiac arrhythmias may come from the spiral wave state of the cardiac electric impulse, and fibrillation appears once the spiral wave state turns into the turbulence state.[23] The spiral wave is a self-sustained wave state. When the excitation spreading along the longside of an infinitely inhomogeneous stripe is cut, the stationary propagation of a wave segment can be observed.[24] But in the ventricle CA model used here, the impulse spreads across the longside of the slice, which makes it impossible to produce a spiral wave under the normal conditions: the electric impulse spreads very fast under the normal conditions, so even the travelling wave is cut, the cut-down point goes out off the system quickly, and then the spiral wave state cannot come into being. Pathological changes will reduce the conduction ability of the conduction system in the ventricle.[21] On the other hand, the decrease of the excitation ability can also reduce the wave speed. To imitate these two situations, only the endocardial cells within y ≤ 5 are set at the excited state initially (t = 0), and the value of K is no longer 2 but is now 14. The other parameters are the same as those in Fig.  3. When the computer simulation is done, the travelling wave is cut down at t = 6 (Fig.  5(a)). The results show that when the cut-down point returns to the endocardium, the endocardial cells have come back to the rest state, so the impulse can spread continually, and the spiral wave finally comes into being (Figs.  5(b) and 5(c)). The further computer simulation results indicate that if there is no more impulse coming from the endocardium, then the spiral wave will be maintained for ever. The potential at field point (130, 40) is recorded again, and the results are shown in Fig.  5(d).

Fig.  5. Distributions of the excited state at (a) t = 6, (b) t = 14, (c) t = 139, the black represents the excited state; (d) the electrogram recorded at (130, 40). The parameters are the same as those in Fig.  3, except K = 14.

Figure  5(d) indicates that when the impulse evolves with spiral state, the potential trajectory appears as follows: (i) the amplitude of the potential is 0.5447, which is much smaller then 2.0039 in Fig.  3(a) under the normal conditions; (ii) the period of the potential is 15 steps, not 33 steps in Fig.  3(a) under the normal conditions; and, (iii) the wide R wave appears, which corresponds to the clinic ECG as the ventricular tachycardia occurs.[6] The last two phenomena have been observed in the clinic.[4] So the above results support the conclusion that the spiral wave state may be the reason for the ventricular tachycardia.

4. Conclusion

In this paper, a CA model considering the three-layer structure is proposed for a ventricle slice. Three types of cells with characteristic repolarization properties are distributed on the three layers. The computer simulation results exhibit the following phenomena. (i) Under normal conditions, the R, S, and T waves appear in the potential trajectory, which corresponds to the standard ECG. This result agrees with the physiological conclusion about the existence of the subepicardium. (ii) When the ischemia occurs in the endocardium, the potential decline of point J and segment ST and the raise of the potential line after the T wave can be observed, which is close to the clinic conclusion. This result comes from the shorter APD and higher resting potential of the ischemia cells in endocardium. (iii) When the impulse evolves with the spiral-wave state, the potential trajectory has a smaller amplitude, shorter period, and wide R wave, which supports the conclusion that the spiral-wave state maybe the reason for the ventricular tachycardia. The above phenomena indicate that the CA model in this paper is appropriate for the ventricle slice, and the work has reference value for further study in the dynamics of the electric impulse of the cardiac.

Reference
1 Kohl P, Sachs F and Franz M R 2005 Cardiac Mechano-Electric Feedback & Arrhythmia Tianjin Tianjin Science & Technology Translation & Publishing Co. 121(in Chinese) [Cited within:1]
2 Sidorov V Y, Aliev R R, Woods M C, Baudenbacher F, Baudenbacher P and Wikswo J P 2003 Phys. Rev. Lett. 91 208104 DOI:10.1103/PhysRevLett.91.208104 [Cited within:2]
3 Pourhasanzade F and Sabzpoushan S H 2010The 3rd IEEE International Conference on Computer Science and Information TechnologyJuly 9–11, 2010Chengdu, China 502 [Cited within:2]
4 Sorgente A and Josephson M E 2012 J. C. Cases 5 e28 DOI:10.1016/j.jccase.2011.09.003 [Cited within:2]
5 Zhu J J, Jiang S Q, Wang W Y, Zhao C, Wu Y H, Luo M and Quan W W 2014 Chin. Phys. B 23 048702 DOI:10.1088/1674-1056/23/4/048702 [Cited within:1]
6 Yan Q, Guan D L, Yang H and Yan H W 2012 A Concise Hand book of Electrocardiogram Beijing People’s Medical Publishing House(in Chinese) [Cited within:8]
7 Chen J X, Mao J W, Hu B B, Xu J R, He Y F, Li Y and Yuan X P 2009 Phys. Rev. E 79 066209 DOI:10.1103/PhysRevE.79.066209 [Cited within:1]
8 Ma J, Ying H P, Pan G W and Pu Z S 2005 Chin. Phys. Lett. 22 2176 DOI:10.1088/0256-307X/22/9/010 [Cited within:1]
9 Dong L F, Bai Z G and He Y F 2012 Acta Phys. Sin. 61 120509 DOI:10.7498/aps.61.120509(in Chinese) [Cited within:1]
10 Gan Z N and Cheng X M 2010 Chin. Phys. B 19 050514 DOI:10.1088/1674-1056/19/5/050514 [Cited within:1]
11 Wolfram S 1983 Rev. Mod. Phys. 55 601 DOI:10.1103/RevModPhys.55.601 [Cited within:1]
12 Makowiec D 2010 Int. J. Mod. Phys. C 21 107 DOI:10.1142/S0129183110015002 [Cited within:1]
13 Moe G K, Rheinboldt W C and Werner Abildskov J A 1964 Am. Heart J. 67 200 DOI:10.1016/0002-8703(64)90371-0 [Cited within:1]
14 Makowiec D 2010 Acta Phys. Pol. B: Proc. Suppl. 3 377 [Cited within:1]
15 Mitchell R H, Bailey A H and Anderson J 1992 IEEE T. Bio-Med. Eng. 39 253 DOI:10.1109/10.125010 [Cited within:1]
16 Zhu H, Sun Y, Yin B S and Zhu D M 2001 Acta Bioph. Sin. 17 123(in Chinese) [Cited within:2]
17 Zhu H, Sun Y, Rajagopal G, Mondry A and Dhar P 2004 BioMed. Eng. OnLine 3 29 DOI:10.1186/1475-925X-3-29 [Cited within:3]
18 Drouin E, Charpentier F, Gauthier C, Laurent K and Marec H L 1995 J. Am. Coll. Cardiol. 26 185 DOI:10.1016/0735-1097(95)00167-X [Cited within:2]
19 Antzelevitch C 2001 Cardiovasc. Res. 50 426 DOI:10.1016/S0008-6363(01)00285-1 [Cited within:3]
20 Tinniswood A D, Furse C M and Gand hi O P 1998 Phys. Med. Biol. 43 2361 DOI:10.1088/0031-9155/43/8/026 [Cited within:2]
21 Yu C G, Bai R, Chen D L and Huang Y 2008 Cardiac Electrophysiology- Foundation and Clinic Wuhan Huazhong University of Science & Technology Press(in Chinese) [Cited within:3]
22 Lu W G, Wang K Q, Zuo W M, Li J and Zhang H G 2011 J. Biomed. Eng. 28 1200(in Chinese) DOI:10.1056/NEJMoa0903753 [Cited within:1]
23 Nasha M P and Panfilov A V 2004 Prog. Biophys. Mol. Bio. 85 501 DOI:10.1016/j.pbiomolbio.2004.01.016 [Cited within:1]
24 Gao X, Zhang H, Zykov V and Bodenschatz E 2014 New J. Phys. 16 033012 DOI:10.1088/1367-2630/16/3/033012 [Cited within:1]