Jian Gao(高見), Changgui Gu(顧長(zhǎng)貴), and Huijie Yang(楊會(huì)杰)

Business School,University of Shanghai for Science and Technology,Shanghai 200093,China

Keywords: ventricular fibrillation,spiral wave,pattern formation,elimination

1. Introduction

Sudden cardiac death is a major cause accounting for over 300000 deaths per year in the USA alone.[1]Under normal conditions,cardiac action potentials originated from the sinoatrial node periodically propagate at a rate of about 75 bpm in tissues.[2,3]In abnormal situations, such as during the initiation of ventricular tachycardia (VT), reentry of the propagation occurs in that a spiral wave emerges. The wave may stay as a stable spiral wave or degenerate into small fragments of waves, i.e., ventricular fibrillation (VF) resulting from a turbulent cardiac electrical activity, which is a major reason for sudden cardiac death.[4,5]During the VF event,the heart is incapable of any concerted pumping action. These VT and VF events do not often terminate spontaneously.[6-8]The termination depends on the dynamic properties of the spiral wave.[9,10]Therefore,spiral waves in cardiac muscle play an essential role in VT and VF events. Consequently,removal of spiral waves is of much practical interest.[11-17]

At present, a common method of clinical treatment for cardiac tremor is electroshock, which can stimulate the thoracic cavity with high voltage electricity(about 4-7 kV).However, the electroshock method has many disadvantages, such as damage to heart tissue, pain to the patient, and low cure rate.[18]In the research of treatment of VT and VF,one of the research directions is spiral wave dynamics,attempting to find effective methods to remove and to suppress the spiral wave by studying the dynamic properties of the spiral wave. For instance,Gaoet al.[19]proposed that the external forcing can suppress spiral waves and turbulence in excitable media, because the external forcing changes the excitability of the media. In Refs.[20,21],the authors reported that phase compression can suppress the spiral wave in excitable and oscillatory media. In Refs.[22,23]the authors studied the effects of electromagnetic induction on spiral waves and target waves in the model of cardiac tissue. Zhanget al.[24]proposed a method to eliminate spiral waves and spatiotemporal chaos using the synchronization transmission technology of network signals.Moreover,suppressing and controlling the pinned spiral wave are more difficult than the meandering spiral wave, because pinned spiral waves are often attracted to a local area (obstacle) such as heterogeneity. In the issue, in Refs. [25-28]researchers reported many effective schemes to control the pinned spiral waves. For instance,Chenet al.[28]have studied the behavior of spiral waves interacting with obstacles and the removal of pinned spiral waves by periodic mechanical deformation. In addition, some methods suppress spiral waves through competition of waves.[29-34]According to the wave competition rule,[35,36]it is required to impose a higher frequency wave to eliminate the spiral waves in excitable media.The appearance of such high-frequency waves will accelerate the beating of the cardiac muscle,[37]therefore, results in the aggravation of VT event. Accordingly, it is urgent to find a better alternative method for removal of spiral waves.

In this article, we propose a method to clear off spiral waves by exerting a global pulsating disturbance on the medium based on two models of cardiac muscle, i.e., the Fitzhugh-Nagumo (FHN) model[38]and the B¨ar model.[39]We also examine our method in the complex Ginzburg-Landau equation (CGLE) model.[40,41]The rest of the article is organized as follows. In Section 2,we introduce the models and the method. In Section 3 the numerical results are presented. The analytical results are given in Section 4. Lastly,conclusion and discussion are presented in Section 5.

2. Models and method

The two-dimensional (2D) cardiac tissue with transmembrane current is often described by the FHN model expressed as follows:[38,42]

where variableAis the trans-membrane potential, variableBis the slow variable for current,whose time ratio isε,Iis the external periodic forcing current,δArepresents an increment ofAof all local points obtained by an external one-off global pulse current,t0represents the moment when the pulse is applied,t+0represents the moment when the pulse ends,t-0indicates the infinitesimal time before the pulse is applied, ?2is?2/?x2+?2/?y2, and the diffusion coefficientDAis set to 1 for simplicity. The kinetics of this model are given byf(A,B)=A-A3/3-Bandg(A,B)=A+β-γB. Hereβandγare parameters of the medium.

In the simulation,the chosen parameters for this excitable medium areε=0.2,β=0.5,andγ=0.8. Equation(1)with those parameters was also discussed in Ref.[21],showing that the medium can support a rigidly rotating spiral wave and the radius of the trajectory of spiral tip is zero. Our model is integrated on a 100×100 grid points with no-flux boundary conditions. A nine-point finite difference scheme is also applied to compute the Laplacian term ?2A,then the discrete equation is calculated via the four-order Runge-Kutta method. The space steps Δxand Δyare both 1,and the time step Δtis 0.01.

To simulate the electrical activity in a normal cardiac muscle, a pacemaker (external periodic forcing current) is placed in the center of the medium to form a stable target wave,and it reads

In order to facilitate the description of the effect of pulse,δAcan be expressed askA0, i.e.,δA=kA0, wherekrepresents a coefficient andA0(≈1.83)represents the amplitude of local oscillation. In this article, the moment of applying the pulse is assumed to be at the origin of time, i.e.,t0=0. In the following sections,we will examine the effects of the key parameterk, which indicates the degree of influence on the membrane potential of different ions.[38]Snapshots plotted in white and black describe the spatial distribution of the variableA. Because the gradient between the maximum value(white)and the minimum value(black)in the medium is so large,it is difficult to display the gray scale.

3. Numerical results

The general dynamics of the medium is shown in the supplementary files,including a target wave state and the appearance of spiral wave. We use the spiral waves in Fig.S1(b)in the supplementary files as the initial situation of the medium.To examine the effects ofδAon the spiral, the value ofkis gradually increased. Figure 1 shows the effects of the parameterkon the spiral. With the increase ofk,although the stable spiral maintains in the medium, the position of the spiral tip differs(Figs.1(a)-1(c)). In particular,the position of the center fork=0 is indicated by a red solid point. Whenk=0.30,the position of the center is slightly altered and close to the red solid point(Fig.1(a)). Whenkis increased to 0.45(Fig.1(b))or 0.60(Fig.1(c)),the position of the center is farther from the red point.

Fig.1. The stabilized states of a spiral wave with selected values of k,(a)k=0.30,(b)k=0.45,(c)k=0.60,respectively.The red point in the center of the medium indicates the location of the spiral center before the pulse is applied. Each snapshot is a part of the spirals in Fig.S1(b)in the supplementary files.

Interestingly,whenkis further increased and exceeds the threshold ˉk1=0.63,the spiral waves break into many smaller spiral waves (Fig. 2(a)). Whenkis further increased and exceeds(or is equal to)the threshold ˉk2=0.76,the spiral wave is eliminated instantaneously, and a target wave is generated in the medium(Fig.2(b)). Therefore, the parameterkaffects the dynamic patterns of the medium.

Fig. 2. The stabilized states of the medium with selected values of k,(a)k=0.70,(b)k=0.77,respectively.

In order to describe the effect of the parameterkon spiral waves in more detail, we show the relationship between the parameterkand the number of spiral tips in the medium(see Fig. 3). One can find that whenkis in the interval [0,0.63],the number of spiral tip is 2, which means that the spiral waves are not broken or eliminated. Whenkis in the interval(0.63,0.76),the number of spiral tip is greater than 2,which means that the spiral waves are broken. Whenkis greater than or equal to 0.76, the number of spiral tip is 0, which means that the spiral waves are eliminated.

Further,we examine the temporal evolutions for the spiral after the disturbance is applied, ifkis less than the thresholdˉk2=0.76,and ifkis greater than or equal to the threshold ˉk2,respectively, in Fig. 4. Initially (t=0-), i.e., before the disturbance is applied,the stabilized spiral wave is observed. In Fig. 4(a), after the disturbance (k=0.70) is applied, the spiral arms are divided into many small segments (t=5), and then a number of smaller spiral waves are generated with each breakpoint as a new spiral tip (t=10 and 25). In Fig. 4(b),after the disturbance(k=0.77)is applied,the spirals(t=0-)suddenly disappear(t=5 and 15),and a target wave is generated(t=25).The temporal evolution in Figs.4(a)and 4(b)are shown in Movie 1 and Movie 2 included in the supplementary files,respectively.

Fig. 3. The relationship between the parameter k and the number of spiral tip.

Fig.4. The temporal evolution of spirals after the disturbance is applied,when k=0.70 in(a)and k=0.77 in(b),respectively. The moment when the disturbance is applied is the origin of time(t=0).

Fig.5.Alternation of the state for a local oscillation before and after the disturbance(k=0.77)is applied. (a)Temporal evolution of membrane potential A for a typical grid point (25, 25) in the medium. (b) Plots of ν versus t for the grid point(25, 25)in the medium, where ν is the frequency of this grid point.

The properties for a selected local oscillation,before and afterδAis applied,are shown in Fig.5,including the temporal evolution(a)and the corresponding frequency(b). Before the application (t ＜0), the frequency is about 0.95×10-1. After the disturbance is applied(t ≥0),the frequency is quickly reduced to 0.60×10-1. Interestingly,the duration for this reduction is about ?t=40,which is about 2.5 times of the period of the target wave.

In addition to common spiral waves above,we also examine whether our method can suppress the pinned spiral waves,i.e., spiral waves, which are pinned (anchored) to localized heterogeneities,are more difficult to suppress and control than the common spiral wave. In order to obtain the pinned spiral wave on the basis of Fig.1(b),we operate on the medium: fix variablesAin a small area near the spiral tip at-1.03 (potential of resting state). Figure 6(a) shows the pinned spiral waves. The temporal evolution of the pinned spiral waves is shown in Movie 3 in the supplementary files. After the pulse withk=0.77 is applied, the pinned spiral waves are eliminated instantaneously, and a target wave is generated by the pacemaker in the medium of Fig. 6(b). The target wave can cross the obstacle and continue to spread. Movie 4 in the supplementary files shows the temporal evolution.

Fig.6.(a)The pinned spiral waves,(b)a target wave with two localized heterogeneities. The red points are the localized heterogeneities.

In the supplementary files,we examine our method based on the B¨ar model and the CGLE model. We find that in these two models, through our method, the spiral tip can also be moved and the spiral wave can be eliminated.

4. Analytical results

There are several questions: Why can this method eliminate spiral waves? Why can the position of the spiral tip be moved whenkis less than the threshold? In this section we try to explain.

To facilitate the analysis,we answer these questions based on the CGLE model (the model and numerical results are in the supplementary files). In Fig. 7(a), the pointO, called the point defect, is the spiral center which is at the fixed point(Re(A),Im(A))=(0,0), the lineΩrotates counterclockwise withOas the center. Except the center pointO, oscillations of all the other points are simple harmonic oscillation,and the closer to the center,the smaller the amplitude.The points(e.g.,P3andP′3) that have the same distance from theOhave the same phase trajectory. All local points (e.g.,P1,P2andP3)in the radiusrhave the same phase,points inr′(e.g.,P′3)and points inr(e.g.,P1,P2orP3) exhibit an anti-phase synchronization,that is,the difference between their phases isπ.

The spiral center is a self-sustaining wave source that does not require external action. To facilitate the analysis,we discretize the medium. Figure 8 shows nine oscillators near the spiral center in a discrete system, whereA(0,0)is the point defect at the center of the spiral wave, and the oscillators symmetric about pointOare anti-phase synchronous,i.e.,A(-1,1)=-A(1,-1),A(-1,0)=-A(1,0),A(-1,-1)=-A(1,1),A(0,1)=-A(0,-1).A(0,0)is at the fixed point (0,0). SinceA(0,0)=0, their dynamic termsA-(1+iβ)|A|2A=0. Now consider the diffusion term ofA(0,0),

Consequently, pointOis always at the fixed point(0, 0), and the point defect can be self-sustaining,resulting in that the area nearOis always in an uneven state, and this uneven state is what the self-sustaining wave source should satisfy.Moreover,disrupting the anti-phase synchronization around the point defect will break the existing balance, and completely breaking the anti-phase synchronization around the center of the spiral wave will cause the point defect to automatically collapse.Our method is to eliminate the spiral wave by breaking the anti-phase synchronization.

Fig.7. The dynamic behavior near the spiral center:(a)A snapshot of a circular region near the center of the spiral. Point O is the center of the spiral, the black solid line Ω is a set of local points where Re(A)is 0,the black arrows indicate the direction of rotation of Ω. The blue dotted line r is any radius of this circular region,and the yellow dotted line r′is another radius on the same diameter. Points P1, P2 and P3 are local points on r,P′3 is the symmetry point of P3 relative to point O.P1,P2 and P3 have the same phase,P3 and P′3 are in anti-phase,i.e.,the difference between their phase is π. (b)A phase diagram. Trajectories C1 and C2 in (b) are the phase trajectories of P1 and P2 in (a), respectively. C3 is the phase trajectory of P3 and P′3.

Fig.8. Nine adjacent oscillators near the point defect in a discrete system. The circles and the square are adjacent local oscillators,the square is the point O in Fig.5(a),Δx and Δy are the space steps.

Figure 9 shows the set of local points where the real or imaginary part ofAis 0,which is called the symbol boundary curve(SBC).The SBC is symmetric about pointObefore theδAis applied. The intersection(O)of the black solid line and the dotted line in Fig.9 is the spiral center, i.e., the point defect,beforeδAis applied. After the disturbance is applied to the medium,the SBC of the real part loses its central symmetry and is close to one end of the SBC of the imaginary part,and the largerk, the closer to the SBC of the imaginary part(Fig.9). For instance,the blue solid line in Fig.9 is nearer to one end of the dotted line than the yellow and red ones. The alteration of the SBC of the real part,on the one hand,leads to the alteration of the position for the point defect,on the other hand, leads to the increase of the curvature of the SBC at the point defect(Fig.9). For example, pointsO1,O2andO3are the new point defects,the curvature atO3is greater than atO2andO1. Simultaneously, the diffusion effect has a strong influence on the place where the curvature of the SBC is large,which causes the curvature to become smaller and smaller,resulting in further alteration of the position of the point defect.

Fig. 9. The set of local points where Re(A)=0 and Im(A)=0. The black solid line and dotted line in are the SBCs of the real and imaginary parts of the initial spiral wave respectively. The red, yellow and blue solid lines,which are the SBCs of the real part after δA is applied,correspond to k=0.50,k=0.85 and k=0.92,respectively. Point O is the point defect of the initial spiral wave,each of points O1,O2 and O3 is point defect at t=0+ after δA is applied.

With our method, whenk ＜1, the SBCs of the real and imaginary parts always have an intersection point(Fig.9),that is, the point defect cannot be directly eliminated. However,whenkis large enough(about 0.92),the curvature of the SBC for the real part(at the defect pointO3)is large enough,so that the effect of the diffusion on this SBC will be strong enough.Eventually, under the effect of diffusion, point defect is removed from the medium along the SBC of the imaginary part.

5. Discussion and conclusion

In conclusion,we have presented a simple method to clear off spiral waves in three models of cardiac muscle by applying a global instantaneous disturbance. We have observed that the spiral wave can be eliminated in a short time, and then, the medium reaches the homogeneous oscillation in about one oscillation period. During this elimination,different from other methods,the oscillation frequency of the medium does not increase, but decreases directly and rapidly. In our method, we can eliminate the spiral wave by increasing the membrane potential globally. This increase is less than the amplitude of the membrane potential (about 80 mV). In addition, our method can eliminate the spiral wave within a short period of time.We believe that the method in this paper can provide a better guide to the clinical treatment of VT and VF.

The traditional method for treating the VT and VF events in the clinic is electroshock that stimulates the thoracic cavity with high voltage electricity(about 4-7 kV).This method brings great pain to the patient and cause great damage to heart tissue and other tissues of the human body.[18]It is widely believed that spirals in cardiac muscle are play a vital role in life-threatening situations(VT and VF).[4,5]Therefore,experts and scholars have conducted a lot of research on how to remove and suppress spiral waves. By most previous methods of the expulsion of spiral waves,[19,29-33]a target wave was created after applying a periodic forcing at a local point,which leads to the removal of the spiral wave from the heart.A target wave whose frequency is lower than that of the spiral wave cannot suppress the spiral wave,and a target wave with a higher frequency is required to suppress the spiral wave.[35,36]However,the emergence of the target wave with a higher frequency will cause the heartbeat to be accelerated again,which will aggravate the VT and directly endanger life. Accordingly,these methods are not clinically feasible. In this article,in the process of eliminating the spiral wave, the frequency of the medium is not increased again.

It is possible to effectively eliminate spiral waves and pinned spiral waves in cardiac tissue by applying a global pulse to increase the trans-membrane potential. Note that the value of the parameterkmust exceed the threshold ˉk2=0.76,otherwise more spiral waves will be generated (Fig. 2(a)),making the VT and VF events more serious.

About the control of spiral tip,Liuet al.[43]have studied the motion of spiral tip controlled by a local periodic forcing imposed on a region around the spiral tip in an excitable medium,and observed three types of trajectories of spiral tip.Chenet al.[44]moved the spiral tip by propagation time delay.In fact,using the method in this article can also make the spiral tip move. Therefore, the use of periodic pulses can make the spiral tip meander,and by adjusting the time interval between pulses,different trajectories of spiral tip can be obtained.