Coexisting Behavior and Status Transition of the Hodgkin-Huxley Cardiac Purkinje Fiber Model Under External AC Injection

The membrane current Im\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_m$$\end{document} of the Hodgkin-Huxley (HH) cardiac Purkinje fiber (CPF) model is usually calculated as direct current (DC). In this paper, a conventional alternating current (AC) is used, namely IAC=Asin(2πft)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_{AC} = \textrm{Asin}(2 \pi ft)$$\end{document}, as the impressed injection. Dynamic characters of the model with different AC parameters and various initial conditions are observed through phase plane orbits, waveforms, bifurcation diagrams, and Lyapunov exponent spectra, which reveal multiple coexisting membrane action potential patterns, corresponding to the coexisting attractors of the same or different periods in the phase diagram. Meanwhile, the model undergoes period, quasi-period and local forward or inverse period-doubling bifurcations with changes in amplitude A or frequency f, which further proves the complex nonlinear property of the AC-injected model. In addition, by changing the external current Im\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_m$$\end{document}, the sodium-ion and potassium-ion equilibrium potentials, i.e. ENa\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{Na}$$\end{document} and EK\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{K}$$\end{document}, respectively, the regularity of the CPF heartbeat frequency is observed. The state transformations of CPF are found between normal, abnormal and sudden cardiac arrest, and the method adjusting from the dangerous state to the normal heartbeat frequency range is investigated. This study may provide a reference for exploring the evolution of nonlinear dynamics in HH CPF model and protecting the health of life and heart.


Introduction
Purkinje fiber is a reticulated terminal fiber that is part of the cardiac conduction system and plays a key role in the synchronization of the left and right ventricles. In 1952, Hodgkin and Huxley [1] and Hodgkin and Keynes [2] established a quantitative equation named Hodgkin-Huxley (HH) model to describe the membrane current of the squid giant axon, it was modified by D. Noble to form the HH cardiac Purkinje fiber (CPF) model applied to heart action potentials [3]. The parameters of the HH CPF model correspond to the intrinsic characteristics of the membrane or its physiological environment, which has attracted many scholars to carry out numerical experiments and theoretical studies. A parametric study on the ionic conductance of the original and simplified CPF model is conducted [4]. Chaotic attractors are shown when the leakage current is changed [5]. The Na + and K + channels of the model are demonstrated as memristors, and the state transitions between locally-active, edge of chaos and locally-passive domain under external DC stimulation are discussed [6].
Cells may encounter interference from various external factors in the actual environment, and respond accordingly by external stimuli. H. Reuter found that during the application of the voltage clamp technique to cardiac Purkinje fibres in a calcium-containing solution, calcium ions carry an appreciable membrane current in the inward direction when the membrane of the Purkinje fibre is depolarized, which may be involved in excitation-contraction coupling [7,8].
In sodium-free solution calcium ions which flow into the fibre activate contraction directly, in sodium-containing solution, calcium ions current serves primarily to fill some intracellular stores from which calcium can be released by moderate depolarization [9]. In order to better simulate the real world and discover richer nonlinear properties, many studies on neurons and neural networks have introduced different external factors, including current [10][11][12], magnetic field [13][14][15], electric fields [16][17][18] and noise [19][20][21], etc. Under various stimuli, complex dynamic phenomena, bursting behavior, coexisting attractors, chaotic attractors, and multistability are discovered. Coexisting asymmetric bursts are found when AC or bipolar pulsed currents are injected into Hindmarsh-Rose neurons [22,23]. Multi-scroll chaotic attractors are generated when Hopfield-type neural networks are stimulated by electromagnetic radiation and multilevel logic pulses [24]. The switching between the two steady states of the model is accomplished by voltage pulses [25]. Two different chaotic states are displayed when the neuron's external current changes [26]. Influenced by electromagnetic induction, multiple firing patterns and multistability are exhibited in neurons [27][28][29]. Mode transitions of neurons are achieved under electromagnetic radiation [30][31][32][33].
Common alternating current is widely used as external electrical stimulation, different electrical modes driven by variational injection parameters are found [11]. Bursting patterns of simplified memristive Hodgkin-Huxley neuron model with sinusoidal currents are explored [12]. The effect of sinusoidal current on the electrical activity of a four-variable Hindmarsh-Rose neuron model under a depolarizing field is investigated [34]. The firing behaviors of coupled systems connected by magnetron memristors under periodic sinusoidal signals are observed and the memory function of the model to the stimulus signal is verified [35]. The bursting behavior appears in Hopfield neural networks when there is a step difference between the external AC frequency and the natural frequency of the system [36]. The deep brain stimulation system is optimized by analyzing the state changes of neurons stimulated by AC [37]. In this paper, instead of constant DC, alternating current I AC = Asin(2 ft) is applied as a continuous injection stimulus to the HH CPF model, multi-equilibrium points, coexisting attractors and other nonlinear interweaving bifurcation are presented when amplitude A and frequency f are changed. Compared with the stimulation of DC supply, the injection of AC into HH CPF model has richer and more peculiar dynamic phenomena. Aiming at the normal heartbeat frequency of 60 − 100 beats/min proposed by the World Health Organization, this paper also studies the effect of applied alternating current I AC on the heartbeat frequency and makes it better or stable, of course, inappropriate I AC settings can also worsen heart rate. At the same time, the parameters of sodium-ion equilibrium potential E Na and potassiumion equilibrium potential E K are adjusted to find the sudden change of heart rate in CPF.
The main contributions of this paper are: • Combined with the original Hodgkin-Huxley circuit and equation, the AC-injection HH CPF model and circuit's layout are reconstructed. • The dynamics of the AC-injected HH CPF model and the stability are analyzed through bifurcation diagrams, Lya-punov exponent spectra, waveforms and phase diagrams with different initial conditions. • In the case of the injection AC I AC = Asin(2 ft) , the coexisting attractors and state changes are observed through the largest Lyapunov exponents (LLEs) with variable frequency f and amplitude A. • Compared with fixed DC injection, AC injection has rich coexisting phenomena and the characteristics of interweaving bifurcation. • The external AC injection current I AC = Asin(2 ft) is applied to the HH CPF model to study and compare the effects on the heartbeat frequency, and combing with ion equilibrium potentials E K and E Na , a method for the transformation between normal and abnormal heartbeat frequency is sought in this paper.
The Hodgkin-Huxley heart Purkinje fiber model induced with AC is proposed in Sect. 2. The evolution and types of AC equilibrium points are analyzed in Sect. 3. In Sect. 4, numerical analysis is used to observe the dynamic characteristics under the change of I AC parameters and initial conditions. Abnormal heartbeat frequency is restored to normal state by AC injection in Sect. 5. Finally, the research work is summarized in Sect. 6.

AC-injected Hodgkin-Huxley Cardiac Purkinje Fiber Model
In the original cardiac Purkinje fiber equation proposed by Noble in Ref. [3], the total membrane current I m consists of the sum of the sodium-ion current I Na , potassium-ion current I K , chloride ion current I An , and the current into the membrane capacitance I C m . The equilibrium potentials of sodium ion ( Na + ), potassium ion ( K + ) and chloride ion are E Na , E K and E An , respectively. The multiplication term before subtracting the Nerntes potential from the membrane potential indicates the conductance of the corresponding ion. Then the specific form of the individual ionic current can be described as: The whole HH CPF model is made up of four-dimensional nonlinear differential equations, where V represents the membrane potential, m and h are used to compute the sodium-ion current I Na , and n is also related to the calculation of the potassium-ion current I K . During the integration process, the values of m, h, and n depend on the previous value of V and react to V again. Considering chloride ions and other anions with different permeability have little effect on the membrane current [3,6] and [38], it is assumed that the I An anion leak current is zero. By injecting a microampere sinusoidal alternating current I AC of the same size as the ionic current to replace the previous membrane current I m in the model, the circuit structure is shown in Fig. 1  where In the model, parameters such as sodium-ion equilibrium potential E Na = 40 mV, potassium-ion equilibrium potential E K = −100 mV remain unchanged. The I AC item is constructed as I AC = Asin(2 ft) , where A represents the amplitude in A, f means the frequency in Hz, and t is the time in ms. When I m is regarded as sinusoidal AC rather than DC with a fixed constant, more interesting dynamic phenomena will appear as amplitude A or frequency f evolves.

Evolutions of Equilibrium Points with Fixed I AC Parameters
The dynamic characteristics of the model can be observed through the equilibrium point. Let the left side of Eq. 2 be zero, the equations for calculating the equilibrium point E of the HH CPF model can be derived as: Correspondingly, the Jacobian matrix of the model at the equilibrium point E is given by: There are different values of Jacobian matrices and eigenvalues at various equilibrium points. And it has been proved that the HH CPF model has at most three equilibrium points when I AC is a fixed value [6]. Considering that the equilibrium points with various initial conditions may be different where there are multiple equilibrium points, six sets of different initial conditions (V 0i , m 0i , h 0i , n 0i ) (i = 1, 2, 3, 4, 5, 6) are chosen to observe the changing trend of the equilibrium point E on the basis of a large number of comparative experiments, which are listed in Eq. 6. When the alternating current with frequency f = 7.93 Hz and amplitude A = 8.18 A evolves in one cycle with a step size of 0.001 * (1∕f ) , the V-values of the equilibrium points computed by the Fsolve function of matlab under six initial conditions are shown in Fig. 2. In the colorful picture, the blue, magenta and cyan trajectories emerge from initial conditions (V 0i , m 0i , h 0i , n 0i ) , i = 1, 2, 3 , respectively, while the yellow, orange and green trajectories correspond to (V 0i , m 0i , h 0i , n 0i ) , where i = 4, 5, 6 . In order to make the image clear, the size of the drawn points is larger, leading to the illusion of two equilibrium points corresponding to the fixed t for the same initial condition. The specific details are local enlarged at the right of the image.
When t increases, I AC = Asin(2 ft) changes periodically, and the loci of the equilibrium points with different initial conditions are distinct. From Fig. 2, we can find the following phenomena:   Fig. 3 The equilibrium-point curves for six groups of initial conditions with frequency f = 7.93 Hz, amplitude 7 A, 7.93137 A and 12 A, which labeled with different colors, respectively. All the 1 sequences on the left represent the V − t curve, and all the 2 sequences on the right represent the V − I AC curve. The evolution of all graphs is in a period of I AC = Asin(2 ft) . It is obvious from the figures that there are 3 equilibrium points in some places If the initial conditions are closer to one of these three types of equilibrium points, the obtained trajectories will continue steadily along the current type of equilibrium point. Otherwise, the solved equilibrium points will jump between different types of equilibrium points, as shown in the partial enlarged view of Fig. 2. For example, the blue, cyan and green curves correspond to the trajectory types of Q 1 , Q 2 and Q 3 , respectively. For the yellow and orange curves, the equilibrium points are distributed not only in Q 2 but also in Q 3 . And the equilibrium points of magenta curve are mainly Q 2 , and sometimes Q 1 .

Evolutions of Equilibrium Points with Different Amplitude A of AC
For the equilibrium-point curve, the frequency only changes its period. In order to visually compare the changes in V-value and the stability evolutions of the equilibrium points, in addition to A = 8.18 A in Fig. 2, the amplitude A of AC is set as 7 A, 7.93137 A and 12 A, respectively. With frequency f = 7.93 Hz, the evolution step size is consistent with the one in Fig. 2. The equilibrium-point curves for different AC amplitudes are still periodic functions with a period of 1/f, and the curves within one period are shown in Fig. 3.
With amplitude A increasing, the variation range of the external current becomes larger, new equilibrium points are calculated on the new current values, then the range of the V-value at the equilibrium points is magnified. It is reflected in the equilibrium-point curves, which generally appear to be stretched in the direction of the V-axis, and the isolated circles in Fig. 3(a1) gradually intersect the sine-like curve in (b1). When I AC is less than −7.93137 A, there is only one type of equilibrium point, so a "channel" appears at (c1). Next, because there are three equilibrium points in the range of [−7.93137, −0.7951] A, our main focus is on the part where I AC ≤ −0.7951 A.
• For (V 01 , m 01 , h 01 , n 01 ) group: When there are three equilibrium points at the same time, the type of equilibrium points with (V 01 , m 01 , h 01 , n 01 ) is always Q 1 . As amplitude A grows, the blue "sinusoid" curve changes from flat to steep, and then becomes a clearly discontinuous curve segments due to the disappearance of the equilibrium points at partial I AC where the shape gradually becomes sharper as the AC amplitude evolves in Fig. 3 It can be seen from the above that the distribution of equilibrium-point types and corresponding range of I AC do not change with amplitude A under different initial conditions. But the equilibrium points jump at the fixed I AC so that the shift in amplitude significantly alter the shape of the V − t curve. With the six amplitude parameters, the evolution waveforms of the equilibrium points with (V 01 , m 01 , h 01 , n 01 ) and (V 05 , m 05 , h 05 , n 05 ) are quite different in the experiments. So, the commonly used initial condition (V 05 , m 05 , h 05 , n 05 ) and the initial condition (V 01 , m 01 , h 01 , n 01 ) furthest away from (V 05 , m 05 , h 05 , n 05 ) among the six sets of initial conditions are selected for further study of the coexisting behavior of attractors in the next section.

Analysis of Stability Types
To reveal the stability of equilibrium points with the six sets of initial conditions, the AC waveform, the V-values of equilibrium points, the first two larger real part of eigenvalues of Jacobian matrix and the imaginary part of the complex eigenvalue are combined and depicted in Fig. 4, when t varies in steps of 0.001 * (1∕f ) over one period of the AC with amplitude A = 8.18 A and frequency f = 7.93 Hz. The solution method of the eigenvalues of the Jacobian matrixs at different equilibrium points is similar to the code 2.11 in Ref. [39].
For V − t curves in the figure, black curve denotes only one equilibrium point exists, and the colored curves show the part with three equilibrium points. Among them, the red, green and blue V − t curves indicate that the solved equilibriumpoint types are the Q 1 , Q 2 , Q 3 , respectively. Interrupted blanks in the V − t curves indicate that the equilibrium points cannot be solved. And the distribution of equilibrium points in Fig. 4 is the same as that in Fig. 2.
With the imaginary and the first two larger real parts of the eigenvalues, the type of equilibrium-point stability can be divided. In Fig. 4, some representative time points are marked with arrows whose colors depend on the equilibriumpoint types. In order to numerically judge the relationship between stability and equilibrium-point types under different initial conditions and analyze the distribution and trend of stability, the equilibrium points, Jacobian eigenvalues, and the type of stability calculated at these time points are listed in Table 1.
According to the observation of the curves in Fig. 4 and the datum in Table 1 , it can be found that: Stable node point The stable equilibrium points only exist in the blue V-value evolution curve of the Q 3 equilibrium points and the black V − t curve at t ∈ [90.4161, 99.3695] ms in Fig. 4 Unstable saddle-focus Both positive and negative values of the imaginary part are drawn in yellow. Combined with the real and imaginary part of the eigenvalues, the full range of the yellow curve in the figure indicates that the stability type of the equilibrium points is unstable saddle-focus. In the fixed range [ Fig. 4(b), the saddlefocus also appears.

Distribution and trend
Under the present parameters, there are 4 types equilibrium-point stability, including stable node point, unstable saddle-focus, the rare equilibrium points with zero eigenvalues and unstable saddle point. With t < 65.0694 ms, the stability is either unstable saddle-focus or unstable saddle point, if the equilibrium point exists. The stability of the Q 1 equilibrium points is unstable saddle point, except that it is unstable saddle-focus at some fixed time or I AC range. The Q 2 equilibrium points uniquely correspond to the stability of unstable saddle point. And the Q 3 equilibrium points must be stable. Therefore, when the V-value of the initial condition is smaller and closer to the V-value of the Q 3 equilibrium points, the range of the stable equilibrium points is larger.
The injection of alternating current makes HH CPF model evolve from single unstable saddle point at t = 0 ms to unstable saddle-focus, stable node point, and even equilibrium points with eigenvalues of 0, showing a variety of dynamic characteristics.

Coexisting Behavior of AC-induced Asymmetric Attractors
By

Dynamic Characteristics of the HH CPF Model with Amplitude A as Variable
When the AC frequency f is set to 10.50 Hz, amplitude A is adjusted from 0 to 12 A with 0.01 interval. Under the chosen two sets of initial conditions, the first two Lyapunov exponent spectra and the bifurcation diagrams with fixed frequency parameters are plotted in Combining the Lyapunov exponent spectra and the bifurcation diagrams, it can be obtained the injecting alternating current allows the CPF model to exhibit diverse coexisting asymmetric behaviors, which include the coexistence of (a) V01, m01, h01, n01) (b) (V02, m02, h02, n02 ( periodic and quasi-periodic attractors, the coexistence of identical periodic attractors with different topologies and the coexistence of quasi-periodic attractors with different topologies. For nine representative amplitude parameters A, the phase diagrams projected on the V − n plane, membrane potential waveforms, and time series of n are shown in Fig. 6, where the initial conditions for the orange and blue orbits are consistent with those used in Fig. 5(b).
When frequency f = 10.50 Hz, the asymmetric bifurcations under the two initial conditions are closely intertwined. In Fig. 6(a), (c), the phase diagrams of the periodic parts in the bifurcation diagrams are shown. As the bifurcation enters the quasi-periodic part, there is the obvious coexistence of period and quasi-period as well as the coexistence of quasiperiod with different topologies, as shown in Fig. 6(d)-(f). Similarly, around the amplitude A = 5.88 A and the beginning of the inverse bifurcation, the coexistence of attractors with different topologies in the same period appear, as shown in Fig. 6(b), (g)-(h). Finally, the critical states of period 1 to period 2 are shown in Fig. 6(i). Various coexisting attractor combinations are also summarized in Table 2.
Other frequency parameters, such as 7.93 Hz and 11.00 Hz, are also fixed during the experiment, and the coexistence of attractors occurs too. However, since the coexistence types of attractors are similar, it is not elaborated here. It can be seen from the bifurcation diagrams and the phase diagrams that the amplitude A not only affects the shape and type of the attractors, but also makes the coexistence behavior appear intermittently, resulting in complex dynamic behaviors such as local inverse perioddoubling bifurcation.

Dynamic Characteristics of the HH CPF Model with Frequency f as Variable
Here, another current parameter amplitude A is fixed as 8. 18 A, frequency f varies within [0, 20] Hz with a step size of 0.01 Hz, and the initial conditions are the same as those used when the AC frequency f is fixed. When frequency f varies, the dynamic analysis of the HH CPF model under the two sets of initial conditions is conducted by the first two Lyapunov exponent spectra in Fig. 7(a) and the bifurcation diagrams in Fig. 7(b). The corresponding relationship between the color and the initial condition is consistent with Fig. 5, as shown in the identification in the figure.
In Fig. 7, it can be observed that the curves of the first two Lyapunov exponents under the initial conditions (−40, 3, 30, 5) and (−80, −0.05, 0.8, 0.463) have a high degree of coincidence. The orange and blue trajectories are gradually the same in the intertexture with amplitude A increasing in the bifurcation diagrams. It is clear that the period-doubling bifurcation appears locally. And with the increase of f, the model goes from the period-1 state through the period-2 pattern into the irregular bifurcation part, and eventually alternates between quasi-period and period.
By observing the Lyapunov exponent spectra and bifurcation diagrams, 9 parameter examples are selected to analyze the specific situation of the coexistence of attractors with different frequency f, respectively. The phase orbits in the V − n plane and the time series of V and n are plotted in Fig. 8.
From Fig. 7(b), it can be known that there are blue period-2 limit cycles coexisting with various orange attractors in a relatively large f range, and their phase diagrams are drawn in Fig. 8(a)-(e). Among them, period-2 and quasi-periodic attractors coexist in Fig. 8(a), and the orange period-8, period-4, period-2, and period-1 attractors in Fig. 8(b)-(e) illustrate the occurrence of local inverse period-doubling bifurcations. Likewise, in the quasi-periodic interleaving part of the bifurcation diagrams, the coexistence of quasi-periodic attractors with different topologies in Fig. 8(f) and the coexistence of quasi-periodic and periodic states in Fig. 8(g)-(i) appear. Therefore, during the frequency change from 0 to 20Hz, a variety of coexisting asymmetric attractors are obtained. The corresponding relationships between their current parameters, Lyapunov exponents and states are also shown in Table 3.
Different membrane potential patterns can be observed in AC-injected HH CPF model with different initial states, which suggests this asymmetric coexisting behavior is triggered by AC stimulation. The complex nonlinear dynamics with a variety of states under different initial conditions also exists in Hindmarsh-Rose model with AC [22], non-ideal active voltage-controlled memristor based Chua's circuit [40], locally active memristive neuron model [41], memristor emulator-based dynamical circuit [42], chaotic system based on two memristor [43] and so on.

Dynamic Characteristics of the HH CPF Model with Amplitude A and Frequency f as Variables
As mentioned, the dynamics of the HH CPF model for AC stimulation depend on two parameters: amplitude A and frequency f, which one or both will affect the behavior exhibited by the model. Therefore, in this section, the two parameters are simultaneously changed within a certain range, and the corresponding dynamic behaviors are analyzed by calculating LLE. By observing the range where periodic and quasi-periodic coexistence occurs in the bifurcation diagrams Figs. 5, and 7, the focus is on the region where amplitude A ∈ [7.00, 10.00] A, frequency f ∈ [9.50, 12.50] Hz, with a step size of 0.01, respectively. The calculated maps of LLE are drawn in Fig. 9. Among them, the gradient color in Fig. 9(a) indicate the development of negative LLE, and red signifies quasiperiodicity. In the selected parameters, the negative values of the LLE are slightly more than positive, and the probability of quasi-periodic occurrence is lower when amplitude A is larger, frequency f is smaller, or frequency f is larger, amplitude A is smaller. In addition, it is found that the states represented by dark blue, cyan, green, yellow, and red are all inclined at about 45 degrees in the figure. Therefore, no matter frequency f or amplitude A is fixed in this range, period and quasi-period will appear alternately in the bifurcation diagrams, that is, the transforms between stability and instability, as shown in the bifurcation diagrams and phase trajectories in the previous section. Figure 9(b) highlight the changing trend and existing range of positive LLE by coloring. In the colorful area, from the blue edge through the green and yellow areas to the red center, the LLE gradually increases. When amplitude A ∈ [8, 8.5] A, and frequency f ∈ [10.5, 11] Hz, the clustered red indicates that the LLE are the largest in the parameters, which may point to a direction for exploring chaos. By comparing the left and right graphs, it is found that in the  range where amplitude A belongs to [9.5, 10] A, and f is [10.5, 11] Hz, there is an obvious coexistence of periodic and quasi-periodic states. Combined with Fig. 9(a)-(b), it is found that the parameter space of the HH CPF model contains multiple shrimp-shaped periodic domains [44][45][46], which are mainly concentrated in the transitional areas from periodic to substantially quasi-periodic. Meanwhile, a large number of "islands" are mutually scattered in the zone of another state, forming the periodic and quasi-periodic interweaving, which means that the stability of the model is sensitive to AC parameters. To further observe the parameter ranges of the coexisting attractors, periodic and quasi-periodic states are marked in yellow and red in Fig. 9(c), respectively. These two colors clearly distinguish periodic from quasi-periodic regions. Unlike the initial conditions (−40, 3, 30, 5) , the quasi-periodic state of the model under the initial conditions (−80, −0.05, 0.8, 0.463) breaks through the diagonally upward periodic state defense line as amplitude A increases, forming a hook-shaped quasi-periodic domain. Moreover, there are some subtle differences between the two figures. These parts do not form a continuous region, but appear as "islands" in small areas. In order to clarify the obvious coexistences of different states and facilitate labeling, Fig. 9(c) is enlarged in Fig. 10.
In the parameter space, 6 parameter points where periodic and quasi-periodic coexist in the previous contents are selected and marked in the Fig. 10. The circled part in yellow means the period, and the circled part in red represents the quasiperiod. Color bars are generated due to discontinuous values, and some periodic states have LLEs greater than 0 but around 0 + , so the red-yellow boundary is above the 0 value. Among the 6 points, most of the coexisting behaviors are caused by unconspicuous periodic or quasi-periodic "islands", and only the quasi-periodic state of point F exist in the hook-shaped red domain where the state types are both the coexistence of period-4 and quasi-period. The differences of state domains under two sets of initial conditions reflect the influences of initial conditions on model's stability and state. The corresponding relationships between the points and the parameters are listed in Table 4.

Restoration of Normal Frequencies for Anomalous HH CPF Model by AC Injection
Cardiac Purkinje fibers play an important role in the initiation and maintenance of ventricular arrhythmias, such as the trigger beats preceding some reentrant tachycardias and ventricular fibrillation caused by long QT syndrome, ischemia and Brugada syndrome [47]. When the heart is abnormal, the action potential duration of the CPF is prolonged or shortened, or even the excitability is lost [48,49]. This means the waveform of the transmembrane voltage in the CPF will be anomalous.  In the absence of any external current stimulation, that is, I m = I AC = 0 A, E Na = 40 mV, E K = −100 mV, the beat rate of the membrane potential of the HH CPF model is 71.5137 beats/min, as shown by the brown curve in Fig. 11. This is in line with the normal frequency of the human heartbeat of 60−100 beats per minute promulgated by the World Health Organization. Once any of the total membrane current I m , sodium ion equilibrium potential E Na and potassium ion potential E K is changed, the heartbeat frequency of HH CPF will be too fast, too slow, or even stop.
• Effect of E Na on heartbeat frequency. When E K = −100 mV, six E Na values are used to show different representative membrane potential states in Fig. 11. According to the calculation formula in the figure, the heartbeat frequency f hb of different curves can be obtained, and the corresponding relationships between them are shown in the legend.
It can be seen from the Fig. 11 that the membrane potential waveforms are abnormal except for the one at E Na = 40 mV. The orange and blue curves represent the action potential during bradycardia at E Na = 36 mV and E Na = 37 mV, respectively. The slow repolarization plateau and electric relaxation phase of membrane potential disappear at E Na = 58.1 mV, see magenta curve. The green and black curves represent two unexcited states resulting from inability to depolarize and incomplete repolarization at E Na = 35 mV and E Na = 60 mV, respectively. In addition, the experimental results show that when E Na gradually moves away from 40 mV, the heartbeat frequency decreases or increases, and finally reaches 0 beats/min. • Effect of I m on heartbeat frequency. For the HH CPF model with different E Na , AC and DC are selected as the external stimulation current respectively. The original Table 4 The detail parameters of the 6 dots in Fig. 10   waveforms and the stimulated waveforms are combined and drawn in green, blue and orange, respectively, as shown in Fig. 12. In order to clearly show the waveform changes after DC and AC stimulation respectively, the selected current parameters make the three curves not coincide, and the stimulated curves are relatively close to the original heartbeat frequency of 71.5137 beats/min. We have also deliberately listed the state transitions of Fig. 12 in Table 5 Fig. 12(d). After AC stimulation, the shape of the waveform in Fig. 12(d) is acceptable and the frequency is still within the normal range, but DC makes the membrane potential abnormal. • Comparison of AC and DC in adjusting heartbeat efficiency. The frequency f = 1.2 Hz is necessary for AC current. If the period of AC 1∕1.2Hz = 833.33 ms is substituted into the formula in Fig. 11, the calculated f hb is 72 beats/min, which is roughly consistent with the heartbeat frequencies of the waveforms after AC stimulation. This means that the AC frequency f has a strong modulation effect on the action potential. For AC amplitude A, its value is based on adjusting the shape of the waveform. The more abnormal the waveform, the more energy it takes to recover. The external positive DC always accelerates the heartbeat frequency f hb , and has a positive correlation with it. The higher the current is, the bigger the f hb is. Therefore, DC is not applicable to the stimulation of tachycardia and cardiac arrest on the tachycardia side. As shown in Fig. 12(e)-(f), the weak DC of 0.1 A makes the situation worse. • Effect of E K on heartbeat frequency. Similarly, when E Na is fixed at 40 mV, E K is adjusted to observe the action potentials. During the experiments, when E K is less than −100 mV, the f hb decreases to 0 beats/min gradually. When E K is relaxedly greater than −100 mV, the f hb increases, and finally f hb = 0 beats/min due to the waveform into a straight line. So the five abnormal E K values are used to compare with the normal E K = −100 mV. The corresponding action potential waveform is plotted in Fig. 13. Combining the heartbeat frequencies and the shape of the waveforms, it can be found that the green and black tracks at E K = −105 mV and E K = −85 mV are unexcited, the orange track at E K = −100.8 mV and the brown waveform at E K = −97 mV represent bradycardia and tachycardia respectively, and the magenta track at E K = −90 mV indicates the abnormal action potential waveform with the disappearance of the slow repolarization plateau and the electrical diastole. The waveforms in Fig. 13 are also stimulated by AC and DC respectively, and the action potentials are redrawn in Fig. 14. In Fig. 14(a)-(b), for bradycardia and cardiac arrest caused by too small E K , the waveforms after DC and AC stimulation are very similar to those under the conditions like Fig. 12(a)-(c). Anyway, the size of DC current cannot be well controlled without experience, which may lead to the heartbeat frequency of the stimulated action potential easily exceeding or not reaching the range of 60-100 beats/ min. Since the influence of AC amplitude on waveform frequency can be ignored, the absolutely normal heartbeat frequency at the f = 1.2 Hz means the AC stimulation is effective. When facing the normal waveform in Fig. 14(c), DC is more likely to make CPF enter a dangerous state than AC. When E K is too large, tachycardia or cardiac arrest occurs, as shown in Fig. 14(d)-(e). At this time, in addition to increasing the heartbeat frequency, DC cannot produce other effective effects, while AC can restore the dangerous states to normal. The advantages of sinusoidal AC compared with DC and the guidance of CPF transformation between normal and abnormal may stimulate the research on clinical treatment of heart disease.

Conclusion
In this paper, a model of HH CPF based on AC stimuli is proposed. With different AC parameters and initial conditions, the AC equilibrium point jumps back and forth between the three equilibrium points. The stability and evolution trend of the equilibrium points are diversified, which makes the Purkinje fiber have the conditions for complex dynamic behavior. Through bifurcation analysis, coexisting asymmetric bifurcations, local forward and inverse period-doubling bifurcations are found. In further numerical analysis, the coexisting behavior of attractors generated by AC stimuli and quasi-periodic attractors with more complex dynamics are revealed. Phase diagrams and waveforms show that multiple types of coexisting attractors appear regardless of the changes in the amplitude or frequency of the AC. Furthermore, the intertexture of periodic and quasi-periodic reflects the range of periodic and quasi-periodic coexistence through successive AC parameter changes. The change of heartbeat frequency may be caused by many factors. The conversion between normal heartbeat, abnormal heartbeat and stopped heartbeat is studied by setting parameters of external currents such as I AC , I m , and parameters E K , E Na in HH CPF model. Changes in amplitude A and frequency f in AC current are observed to induce variations in cardiac depolarization patterns, somewhat consistent with the excitation-contraction coupling resulting from the calcium ion solution injection concentration in H. Reuter. We will continue to study their correlation in the future, these will provide more reference value for heart protection and disease prevention.