Numerical Analysis of Subcritical Hopf Bifurcations in the Two Dimensional FitzHugh-Nagumo Model

It had been shown that the transition from a rigidly rotating spiral wave to a meandering spiral wave is via a Hopf bifurcation. Many studies have shown that these bifurcations are supercritical, but, by using simulations in a comoving frame of reference, we present numerical results which show that subcritical bifurcations are also present within FitzHugh-Nagumo. We show that a hysteresis region is present at the boundary of the rigidly rotating spiral waves and the meandering spiral waves for a particular set of parameters, a feature of FitzHugh-Nagumo that has previously not been reported. Furthermore, we present a evidence that this bifurcation is highly sensitive to initial conditions, and it is possible to convert one solution in the hysteresis loop to the other.


I. INTRODUCTION
Spiral waves occur naturally in many physical, chemical and biological systems [1][2][3][4][5][6][7][8][9][10]. The motion and behavior of such waves can enlighten certain characteristics of the systems in which they occur. For instance, in cardiac tissue, the presence of these rotating spiral waves (also known as autowaves, and rotors) indicates that there is an abnormality in the hearts natural rhythm (an arrhythmia) [11][12][13][14][15]. In most electro-chemical systems, such as cardiac tissue or neurological systems, excitable properties are an essential part of creating and sustaining spiral waves. The cells ability to be stimulated in response to external energy is critical in the life cycle of spiral waves [16].
Excitable systems, such as the propagation of electrical energy along nerves, have been studied mathematically since 1940 using parameter dependent mathematical models. Spiral waves were first observed by Wiener and Rossenblueth, who created the first finite automata model to simulate spiral wave activity [2]. Since then, many models of cell excitation have been developed and studied. A classic system in modelling cellular excitation is the Hodgkin-Huxley model of nerve excitation developed in the early 1950's [17]. This model simulates the electrical energy passing through a single cell [18]. It does not however simulate how each cell reacts with other cells that are part of the whole excitable medium. In order to explore this, spatial variables must be implemented and this is usually done using reaction-diffusion models. Although the development of this type of model was initiated in the early part of the 20th century [19], it was Turing in 1952 that used them to study interactions between chemical compounds [1]. Around the same time, Richard FitzHugh was developing a mathematical model to mimic threshold phenomena in the nerve membrane * sehgals@hope.ac.uk † foulkea@hope.ac.uk [20]. In 1961, FitzHugh published a paper suggesting a model of nerve cell excitation, a simplified version of the Hodgkin-Huxley model, influenced by the Van der Pol oscillator equations [7]. In 1962, an equivalent model was published by Nagumo et. al. [7]. The system of equations is now known as the FitzHugh-Nagumo model (FHN). To this day, many new models are developed using FHN as a base [21]. Mathematically, spiral waves are parameter dependent, spatio-temporal solutions to reaction-diffusion equations. The motion of spiral waves is extremely important in understanding the dynamical behavior of the wave, with differing types of motion of the waves representing differing types of physical phenomena. There are several main categories of motion of spiral waves, which are classified by considering how the tip of the spiral wave behaves and moves around the medium. The tip of the spiral wave can be defined as the intersection of two isolines in the excitation and inhibitor fields. If the tip traces out a perfect circle around a fixed center of rotation, then it is known as rigidly rotating. A property of rigid rotation is that the shape of the arm of the spiral is fixed and the motion is periodic. Another type of motion is known as meander. This is quasi-periodic with the arm of the spiral periodically changing shape and the tip of the spiral tracing out epicycloid type patterns [22][23][24][25].
Other types of motion include hypermeander (chaotic motion) and drift (motion around a moving center of rotation) [26][27][28][29]. We are not considering these types of motion in this publication.
The transition between one type of motion to another via a parameter change has been of great interest since the study of spiral waves began. In particular, the transition from rigid rotation to meander has been shown numerically to be via a Hopf bifurcation [30][31][32]. In particular, several authors have noted that this transition is via a supercritical Hopf bifurcation, in which a stable fixed point changes to unstable and a stable limit cycle is formed. A key feature of this transition is that the growth of the limit cycle from the bifurcation point is proportional to the square root of the varying parame-ter. It must be noted that a rigidly rotating spiral wave in the fixed (laboratory) frame of reference, is represented as a stationary wave in the comoving frame of reference. In the quotient space, a rigidly rotating spiral wave is simply a fixed point.
Here we present work on a numerical approach to study the transition from rigid rotation (RW) to meander (MRW), and show that in the FHN system of equations, there are regions within the parameter space in which the Hopf bifurcations are in fact subcritical. We analyze these results to confirm their validity and show that within the hysteresis regions where there are two solutions relating to the same set of parameters, it is possible to convert one of these solutions to the other.

A. Theoretical review
Our approach to studying the transition from RW to MRW is motivated by the restrictions of earlier numerical studies in studying large core spirals [30]. The problems with large core spirals is that if we study them in a numerical box that is fixed in space, then that box will need to be very large for spiral waves with large cores. This will result in very computationally expensive simulations which will take a relatively long time to complete.
The solution to this is to use the technique of simulating spiral waves in a frame of reference (FOR) that is moving with the tip of the wave. Foulkes and Biktashev [34] published a method that could achieve not only simulations for RW but also MRW, something that other authors were not able to achieve [35,36]. This means that we could afford a much smaller numerical box in which to conduct the simulations as the tip of the wave never reaches the boundaries of the box and, given a large enough box, the boundaries will not have an influence on the resulting spiral wave.
We review below the main results from Foulkes and Biktashev. Further details relating to these methods can be found in Foulkes and Biktashev [34].
Let us consider the reaction-diffusion system (RDS) of equations, where u = u(r, t) = (u 1 , u 2 , . . . , u l ) ⊤ ∈ R l , l ≥ 2, r = (x, y) ⊤ ∈ R 2 , f ∈ R l , and D ∈ R l×l is the matrix of diffusion coefficients. Foulkes and Biktashev considered a RDS that contained symmetry breaking perturbations, which forced the spiral wave solution to drift. Since drift is of no concern in this work, we consider the unperturbed RDS.
This system is invariant under the Euclidean group SE(2), the group of the isometric transformations of the plane R 2 → R 2 . This means that if u(r, t) is a solution to equation (1), thenũ(r, t) is another solution to (1) which is given byũ (r, t) = T (g)u(r, t), ∀g ∈ SE(2), where action T (g) of g ∈ SE(2) on the function u is defined as The action of the group element, g ∈ SE(2), on a spiral wave solution is such that it translates and rotates the spiral wave. It has been shown that spiral waves are equivariant under Euclidean symmetry, and that if we apply group action to a spiral wave solution, then we still have that same solution only it is now in a different position and has a different orientation [23,34]. We can therefore choose group elements, g, such that the tip of the spiral wave is always in a fixed position and has a fixed orientation. Hence, for each rigidly rotating spiral wave there is a corresponding g such that it will be stationary in the FOR, i.e. they are independent of time. Hence the isotropy group will be trivial. For meandering spiral waves, the shape of the wave changes periodically over time and therefore, we seek the set of group elements preserve the position and orientation of the tip of the meandering spiral wave, and this set of elements will be used to create the quotient data.
By considering the spiral wave solution to (1) in an appropriate Banach space, and splitting out the motion of the spiral wave across and perpendicular to an appropriate representative manifold, defined such that the tip of the spiral wave always remains on it, then we obtain the following system, where l 1 , l 2 , l 3 ∈ {1, 2, . . . , n} with l 1 = l 2 , v = v(r, t) = (v 1 , v 2 , . . . , v l ) T ∈ R l is the spiral wave solution in a frame of reference that is moving with the tip of the spiral wave, c(t) = (c x (t), c y (t)) is the translational velocity of the spiral wave and ω(t) is its rotational velocity. The position and orientation of the tip are given by R and θ respectively, hence equations (5) are equations of motion of the tip of the spiral wave. The fixed parameter, τ , is the matrix τ = 0 −1 1 0 , meaning that exp(τ θ) is the rotation matrix where the rotation is by angle θ. We note that equation (2) is a reaction-diffusionadvection system of equations [34,37], whose spiral wave solutions v(r, t) are such that their tip remains on the manifold. The conditions (3) & (4) are the tip pinning conditions. The tip can be pinned at any point within the numerical box, but the definition here is such that they are pinned at the origin, which we place at the center of the numerical box, at an orientation determined by (4).

B. Reaction Kinetics
We will be using the FitzHugh-Nagumo (FHN) model for the simulations within the work presented here. This is a two variable, parameter dependent reaction-diffusion type model. Since this means that l = 2, then we let u l1 = u l3 = u 1 and u l2 = u 2 . In the equation (1), f (u) defines the model kinetics, which, for the FHN model, are given by, We see that there are two variables -u 1 (r, t), the excitation variable, and u 2 (r, t), the inhibitory variabletogether with three parameter, β, γ, and ǫ.
The parameters are varied to give a variety of solutions. Winfree [22] illustrated the spiral wave solutions in a parametric portrait, based on fixing the parameter γ = 0.5, and varying the remaining two parameters, β and ǫ, to get a plethora of solutions within a section of the parameter space, including regions of hypermeander and plane waves. In general, we usually have |ǫ| ≪ 1.

C. Numerical Implementation
Numerical implementation of this system is also detailed in [34] and resulted in software called EZRide [38]. Operator splitting was utilized to simplify the otherwise complicated equations. We can rewrite equation (2) as, such that, For discretization, we have a constant time step, ∆ t , and space step, ∆ x , covering the square spatial domain (x, y) ∈ [−L/2, L/2] 2 , where L is the length of the box in space units. The domain is divided into smaller squares by dividing the x and y axes into N x and N y intervals respectively. For our purposes, we let This means that there will be N + 1 points along each axis.
LetF andÂ be discretizations of F and A respectively. Our numerical computations are as follows, The time step, ∆ t , is given by where t s is the ratio of the time step to the diffusion stability limit, usually taken to be t s = 0.1 [39].

D. Reaction-diffusion step
Foulkes et al [34] and Barkley [21] used the initial steps for computation just as the same as used in the Barkley's EZ-SPIRAL software. Further, they added more numerical computational steps to it. So, for the reaction diffusion part, a first order accurate forward Euler method was used to calculate the temporal derivatives, and the Five Point Laplacian method for the Laplacian. E. Advection step; to calculate cx, cy and ω An upwind second-order accurate approximation of the spatial derivatives is used inÂ. In this step, the discretization ofV k+1 at the tip pinning points is used to calculateĉ k+1 x ,ĉ k+1 y andω k+1 so that after every step, the tip pinning conditions are correctly satisfied.

F. Tip pinning conditions
Pinning the tip of a meandering spiral wave was achieved by choosing two isolines, one for excitation and one for inhibitory, whose values are located within the range of values of both excitation and inhibitory variables. The values can be determined by considering the phase portrait relating to the kinetics used. The full details of the choice of tip pinning conditions are given in Foulkes & Biktashev [34], and the reader should refer these full explanations.
In summary, discretization of the tip pinning condi- where the subscripts are the variable identifier, and the superscripts represent the spatial and temporal variables respectively. The choice of tip pinning condition ensures that the tip of the spiral wave is fixed at a certain orientation and position regardless of whether the spiral wave is rigidly rotating or meandering.

G. Reconstruction of tip trajectory
The EZRide software has an in-built algorithm for reconstructing a tip in the comoving FOR [34]. Further, to check our calculations and draw the tip trajectories, we considered equation (5) and solved it using the numerical scheme which is given as Therefore, the equations (7), (8) and (9) represents the reconstructed tip trajectory.

H. Other details
Throughout the studies conducted here, we used Neumann boundary conditions. Several authors have noted the effect of the boundaries upon the behavior of the spiral waves, and different boundary conditions can lead to different solutions. However, it has been noted [40][41][42] that provided the tip of the spiral is sufficiently far from the boundary, the boundary will not affect the overall dynamics of the spiral wave determined by its tip. This in turn is related to the response functions of the spiral wave, which are localized at the tip [40]. Therefore, we can use either Neumann or Dirichlet boundary conditions, or a mixture of both, within our simulations.
As a guide to what sufficiently far from the boundary means, we usually take this to mean that there is a full wavelength between the tip and the boundary, i.e. the distance measured from the tip of the wave to the part of the arm of spiral that has the same u 1 and u 2 values as at the tip and has rotated around by 2π.
Another consideration of the numerical implementation of the system (2-4, 6) is that of stability. It was shown [34] via a von Neumann analysis that we require, for stability in the calculation of c x , c y and ω. If, for a particular timestep, the values of c x and c y went beyond these limits, then their values were restricted to these limits, i.e. they were not allowed to go any higher than those calculated in the stability limits above. One of the techniques employed in order to reduce instability was that as the advection terms are "switched on", then the tip of the spiral was moved immediately to the tip fixation point. It is also important to note that here ω was not calculated initially until the orientation of the spiral wave was met, at which point ω was given an initial value of zero, and then calculated using the method described earlier.
Allocating ω its limiting value eventually led to instabilities within the system. Therefore, unlike c x and c y , if ω went beyond its limiting value, then it was allocated a value of zero rather than its limiting value.

III. NUMERICAL BIFURCATION APPROACH
We aim to show the nature of the bifurcation responsible for the transition of spiral waves from rigid rotation to meander by generating solutions in a frame of reference that is comoving with the tip of the spiral wave. This means that even for large core spiral waves, we can still afford a relatively small computational space. Furthermore, Foulkes & Biktashev [34] showed that within the solutions to (2)(3)(4)6), the quotient data, consisting of the dynamic variables c, ω, form limit cycle solutions. We therefore study the growth of these limit cycle solution from the onset of meander, and the nature of the growth of this data will indicate the type of bifurcation taking place.

A. Methodology: general overview
Consider the parameter portrait from Winfree for FHN, as shown in FIG.1. We see that there are different regimes of types of motion of spiral waves, according to values of β and ǫ. We decided to study the growth of the limit cycles relating to meandering spiral waves, by analyzing the quotient data for a range of spirals which, on varying one of the parameters, go from rigid rotation to meander and then back to rigid rotation. From FIG.1, we decided to fix ǫ = 0.2, and starting at β = 0.570 within the upper rigidly rotation space, we varied the βparameter in steps of ∆ β = 0.001 to get a spiral wave solution for each value of β. Our initial thoughts were that this choice of ∆ β was sufficient to generate a range of solutions which show the nature of the bifurcation. As we will see in later sections, the choice of ∆ β will lead to qualitatively similar results, but quantitatively different ones, showing sensitive dependence on initial conditions.
Each simulation records not only the quotient data, but also the final conditions from the end of the simulation. The initial conditions for each new simulation is taken as the same as the final conditions of the previous simulation. When we do this, there is always a "settling down" period from the initial conditions to the current solution. Therefore, this initial transient period needs to be eliminated. Although this transient period can sometimes vary in length, only the data taken after at least five complete periods had occurred. In some, but not many, cases, it was obvious that the transient period oc-FIG. 1: Parametric Portrait for FHN for γ = 0.5 [22].
The boundaries labeled with a preceding ∂ are boundaries between the following types of waves: ∂Pno wave -propagating wave; ∂R -propagating waverigidly rotating spiral; ∂M -rigidly rotating spiralmeandering spiral; ∂C -meandering spiralhypermeander spiral (chaotic region). The vertical red line represents the set of parameters that was taken in the main study of this publication. The parameter, λ 0 , is the scaling factor used to draw the tip trajectories on the parametric portrait (see Winfree, 1991, for further details) curred for more than five periods and therefore the data for those simulations was analyzed to see where the transient period had ceased. The remaining data was then analyzed in the usual way. Once all the simulations were completed, the quotient data was then analyzed and the "size" of the limit cycles were then plotted against the parameter, β, which we shall also call the bifurcation parameter. To do this, we need to define this size and so introduce below the quotient size, Q s , of the limit cycle of a spiral wave.
Furthermore, we shall consider the bifurcation points which are the values of the β parameter at which the bifurcation takes place, and relate directly to the last simulation for which there is no limit cycle.

B. Quotient size
The shape of the limit cycle in FIG.2 is irregular. In previous studies [43], the radius of the limit cycles were measured, but this technique cannot be applied here, due to the irregularity of the shape of the limit cycle. We therefore calculated the "distance" around the limit cycle from one point on the cycle all the way around back to that point by calculating the arc length. We know that arc length of a function q(t) in the interval t ∈ [t 1 , t 2 ] is given by, For the spiral wave limit cycles, we let q(t) = (c x (t), c y (t), ω(t)), and if we take the integral over one whole period of the spiral wave, T , then Since we do not know the exact form of c x , c y and ω, then we need to use the discretized form of the arc length formula, where N = T /∆ t , j is a starting point on the limit cycle, andĉ i x is the discretized value of c x at the i-th step, with similar notations for c y and ω.
As noted earlier, due to the transient period of the spiral wave when the simulation first starts, we neglect the first five periods of the simulation. If there are n full periods left of the simulation, the quotient size of those n periods were calculated separately using equation (10), and then Q s will be the average of those n quotient sizes.
In the case where there is rigid rotation, then Q s = 0 since once the transient period has passed, the values of c x (t), c y (t), and ω(t) all remain constant. Constant quotient is indeed an indication of rigid rotation, and once this constancy is detected, then the next simulation is started after 50 time steps. This ensured that the simulation had settled to a solution that was rigidly rotating.

IV. RESULTS
For all the initial simulations, we have used L = 30, ∆ x = 1 5 , ∆ t = 0.001 and t s = 0.1. These numerical parameter values were carefully taken from Foulkes et al, 2009, so that we get not only accurate simulations but computational fast generation of these simulations [34]. In the comoving FOR, we observed that for rigidly rotating spiral waves, our solution becomes stationary and that the quotient data stabilizes to constant values for c x , c y , and ω. As β varies, quotient system is no longer constant and is in fact periodic. This corresponds to meandering spiral wave solutions, exhibiting complicated quasiperiodic motion. C are MRW with outward facing "petals"; D is near the 1:1 resonance line and has an extremely large core radius; and E is MRW with inward facing "petals".
The initial results are shown in FIG.3, where the parameter, β, is plotted against the quotient size, Q s . The simulations started from β = 0.570 and β increased with a step size ∆ β = 0.001. The quotient size was zero for 0.570 ≤ β ≤ 0.601, indicating rigid rotation. For 0.602 ≤ β ≤ 0.967, we had Q s > 0, meaning that the limit cycles were present and solutions are classed as meandering spiral waves. Increasing β further, we found that for 0.968 ≤ β ≤ 0.990, Q s = 0 meaning rigidly rotating spiral waves were present.
The change in the dynamic behavior from RW to MRW is due to the Hopf bifurcations arising from the stable steady states, RW [43]. These Hopf bifurcation take place at β = 0.601 and β = 0.968. Furthermore, as we increase β from the critical point β = 0.601, there is a sudden change in the qualitative behavior of the system from a stable steady state (RW) to an oscillatory state (MRW).
It is clear that as we vary β within the MRW region, i.e. for 0.602 ≤ β ≤ 0.967, then the growth of the limit cycles, as measured by Q s , initially increases from zero up to a maximum, and then decreases again to zero. In order to see whether the arc formed by plotting β against Q s contained any special features, such as the maximum of the arc relating to a specific type of meander, we chose some values of β and plotted their reconstructed tip trajectories from the quotient data. We illustrate the various types of solutions in FIG.3. There were no particular surprises in this analysis and the results tied in with what Winfree observed in his parametric portrait of FHN that we reproduce in FIG.1. A feature of this plot is that as we increase β from the left hand Hopf bifurcation point, there is a significant gap between values of Q s which is not in line with the rest of the plot. We see that once β has increased beyond 0.601, Q s grows slowly and then growth suddenly accelerates very quickly, before decelerating near the peak of the plot. It then accelerates (this time in the negative direction of Q s ) down to the right hand bifurcation point. We therefore decided to look closer at this gap near the left hand bifurcation point.
According to the previous studies which were conducted using Barkley's model [43] or the Belousov-Zhabotinksy reaction [44], it was observed that a supercritical Hopf bifurcation is responsible for the transition from rigid rotation to meander. However, we observed a discontinuous jump in the growth of quotient size near the bifurcation points, being more prominent on the left hand side. Thus, the initial observation is that the result does not tie with the analysis of supercritical Hopf bifurcations due to the absence of square root characteristic.
However, the discontinuous jump observed in the bifurcation diagram depicted in the FIG.3, signals a subcritical Hopf bifurcation [45]. It is also known that if we vary our bifurcation parameter back and forth across the Hopf bifurcation point, we wouldn't expect to jump back to the same value of β, where it lost its stability. We therefore decided to run the simulations again, but this time starting at β = 0.990. If there is a subcritical Hopf bifurcation present, then we should observe hysteresis, which is associated with the bistable region [45].

V. HYSTERESIS
The simulations were now run backwards from β starting from 0.990 and decreasing in step of ∆ β = 0.001 again to 0.570. We performed the same calculations as in the previous section to calculate Q s for all the values of β and observed its growth against β, as shown in FIG ever, for these simulations, the Hopf bifurcation point on the left-hand side in the reverse case has now shifted to β = 0.592, whereas on the right-hand side the Hopf bifurcation point shifts very slightly to left at β = 0.967.
If we combine the data for the forward and backward runs, we see that the data coincide exactly except for only a small range of values. We show this in FIG.4. This depicts the presence of the hysteresis region in which steady and oscillating states coexist. Hence, this region is associated with bistability. The presence of the hysteresis zone is an important characteristic for subcritical Hopf bifurcations, where the system can be in more than one state [45][46][47]. In our study, the hysteresis region exhibits both RW and MRW solutions for a small range of β. Therefore, the existence of both the solutions and change in the β-value corresponding to the equilibrium clearly corresponds to the case of subcritical Hopf bifurcation.
Consider only the hysteresis region on the left hand side of the FIG.4. We note that as the lower branch "jumps" to the higher branch, there are only a few data points on this right hand jump. These data points represents simulations in which both the solutions are meandering spiral waves. As these data points are of low magnitude, i.e they are very few, then we decided to see if we could increase their magnitude by decreasing the step size in β, ∆ β . We re-ran the simulation in exactly the same way as before, but this time we took ∆ β to be 0.0005 and 0.0002. The effect of this was not as expected, and rather than getting more data points in the same hysteresis loop, we had fewer points in the meandermeander region and a change in the size of the hysteresis region.
One such explanation for this is that for larger β-steps, the step is large enough to induce a perturbation such that there is sensitivity to initial conditions causing a premature switch from RW to MRW motion. This de-scribes the effect of sensitivity to initial conditions [48]. Furthermore, the bifurcation points for the forward and the reverse run across the β values are highlighted in  Further, analyzing the behavior across the hysteresis region, we decided to run the simulations with the same range of β but now for a different value of epsilon. We chose ǫ = 0.25 and varied β both forward and backward across the range. From FIG.6, we can see that both the graphs coincide with each other apart for the few values of β near the bifurcation point, which seems more obvious on the left hand side, similar to our original result. The width of the hysteresis region represents the parameter range where the system is bistable. It can be clearly seen that the width of the hysteresis zone is also affected with the change in the value of ǫ. Here, we discovered that the width of the hysteresis zone decreases with the increase in ǫ. The presence of hysteresis also confirms that the transition to instability is subcritical.

VI. CONVERSION
Our next concern was to consider only the hysteresis region and to check if it is possible, for a pair of solutions corresponding to the same β value, to convert from one solution to another. We only tested for the region on one side of the original bifurcation diagram (∆ β =0.001), being more prominent on the left hand side, shown in FIG.4. We described the conversion of one type of solution to another under perturbation such as single shock defibrillation [42,49,50].
The main focus here is on the bistable region where there exists both steady and periodic states. In our case, we investigate a parameter region in an excitable media, where RW and MRW solutions coexist for the same parameter values. They differ significantly as rigid motion depicts steady states whereas meander depicts oscillatory states.
It can be seen that for a particular value of β within the hysteresis region given by 0.592 < β < 0.602, depicted in FIG.4, we have two types of solutions. Here, we consider the transition between the two solutions and chose β = 0.595, which exhibits steady state (RW) while we run the simulations forward across the β-range whereas it shows the periodic behavior (MRW) with the reverse run. Now, to test if the spiral wave solutions can be converted from one type to another, there are many ways in which this can be achieved. We decided to use the method of single shock defibrillation by adding a perturbation to the system. Foulkes et al. [42] have shown that we can convert one type of spiral into another by means of a single shock, and so we utilize techniques from their work.
We added a perturbation to the u-field, uniformly throughout space at a specified time T [42].
where u = (u, v) and h(t) = (h u (t), h v (t)) ⊤ is a timedependent perturbation, defined as where A is a constant. Since we consider solutions in a comoving frame of reference, we consider solutions to the reaction-diffusionadvection system of equations [34] ∂v ∂t whereh(t) is defined as Since h is dependent on time, thenh = h. We initially applied a shock of minimum amplitudes of A = 0.1, 0.5, and observed that it was not sufficient to convert from a rigidly rotating to a meandering spiral wave. Further, increasing it to A = 1.0, we could see the conversion from rigidly rotating wave to a meandering wave pattern. The conversion was also possible with the shock amplitude of A = 1.2 whereas an increase in the amplitude to A = 1.3 resulted in the elimination of spiral wave activity, thereby causing defibrillation.
All these shocks were applied after 40,000 steps to check for the conversion. After applying the shock, we can see the change in the steady spiral state changes to a periodic state. It concludes that the shock successfully converted a rigidly rotating spiral wave solution into a meandering solution, as shown in FIG.7(top).
We also tested if we could shock from a periodic solution to a stable solution. We now considered the meandering spiral wave solution for β = 0.595 and applied a shock of same amplitude to it. In this case, we did not observe the conversion from meandering to rigid rotation by a shock of any amplitude. This shows that conversion is not possible from meander to rigid state.
Alonso et al. [51] noted that the meandering spiral waves with the inward facing petals in two-dimensions always have negative filament tension in three-dimensional case, and meandering solutions with outward facing petals have positive filament tension. Hence, this gives us an explanation for the conversion not being possible in our case, as shown in FIG.7. For β = 0.595, we can see that these are the outward meandering spiral waves corresponding to positive tension in 3D. In addition, Foulkes et al. [42] have shown that it is only possible to convert a wave with negative filament tension to a wave with positive filament tension. Therefore, for β = 0.595, it is possible to convert from rigid rotation to meander but not conversely.

VII. CONVERGENCE ANALYSIS
We also need to test that the numerical results that we observed are in fact a true representation of the analytical results. Therefore, we performed several tests to see if the solutions we originally observed converged by changing ∆ x , ∆ t and the domain size, L, to prove that it is an accurate representation of the true solution. We note that ∆ x = L/N x and ∆ t = t s (∆ x ) 2 /4, so to vary ∆ x while keeping L fixed, we varied N x . Similarly, if we vary ∆ x , and keep t s fixed, then ∆ t will also vary. If we vary L, then in order to keep ∆ x fixed, we must vary N x . A summary of the four tests we conducted is shown below.
• ∆ x varied, t s fixed implying ∆ t varies, L fixed; • ∆ x varied, t s varied so that ∆ t is fixed, L fixed; • ∆ x fixed, t s varied so that ∆ t varied, L fixed; and • ∆ x , t s and ∆ t are fixed, L varies.
The default values of the various numerical parameters are L = 30, ∆ x = 0.2, ∆ t = 0.001, t s = 0.1. FIG.8 shows a selection of the convergence tests that we conducted, plotting β against Q s .
The first test, shown in FIG.8(a)-(b) was to vary ∆ x while keeping t s fixed, meaning that ∆ t varies too. It was observed that with the decrease in the space step, the size of the left hand hysteresis region becomes more prominent. We also noted that other changes included an increase in the maximum value of Q s as a function of β, and the bifurcation point positions also changed. However, despite all these changes, the overall shape of the graph depicting the growth of quotient size against β remains the same.
The next test varied the ∆ x but this time we ensured that ∆ t was fixed by varying t s . Partial results are shown in FIG.8(c)-(d). As with the first test, certain features such as the location of the bifurcation points, the maximum value of Q s , etc, changed as we varied ∆ x , but the overall shape of the arc remained qualitatively the same.
The third test varied ∆ t while keeping all other numerical parameters fixed, FIG.8(e)-(f). While we ran the simulations for different time steps, instabilities were present for the larger time steps but these instabilities tend to decrease in size and frequency as we decreased the time step. But, yet again, while features such as the location of the bifurcation points, maximum value of Q s , etc changed slightly, the overall shape of the curve stayed the same.
The final test is the convergence in box size and some of the results are shown in FIG.8(g)-(h). We considered two variants: first with half the box size and the other one when we double the original box size. In both the cases, we got almost the same Hopf bifurcation points. Only the properties (area, β with maximum quotient size) affecting quotient size differs while we increase the box size.
These tests, show that as we vary the numerical parameters, some of the features of the results changes such as the location of the Hopf bifurcation points, the maximum value of Q s . This should be expected. With varying numerical accuracy comes variation in the error associated with the corresponding results. But one thing remains clear, and that was the overall shape of the curve remained similar throughout the studies. It is clearly evident that the original data represents a true likeness of the actual data.

VIII. DISCUSSION
We have examined the type of bifurcation responsible for the transition from rigidly rotating spiral wave solutions to meandering spiral wave solutions using FHN model.
We chose the FHN model, where β was considered as a bifurcation parameter. It is also important to note that we conducted numerical simulations in the FOR moving with the tip of a spiral wave. This helped us to overcome illustrating the solutions have settled. This is illustrated by using parallel black lines on the t-ω plots.
the drawbacks of previous studies which were unable to quickly study the large core spirals [43]. We were able to successfully study the bifurcation analysis of spiral waves for the full range of β depicting the transition from rigid to meander and then back to rigid motion.
Since, our analysis was conducted in comoving FOR, we study the limit cycle solutions within meandering spiral waves described by the advection coefficients. Here, the limit cycles are not of any particular shape. Therefore, we presented a new technique to check for the type of bifurcation responsible for the transition from rigid to meander motion in FHN model with a specific set of parameters. We calculated the quotient size of limit cycles and plotted it as a function of β. With this approach, we discovered that unlike other authors [43,44], there is a subcritical Hopf bifurcation responsible for the transition from rigid rotation (RW) to meander states (MRW) in the parameter space that we studied. Near the Hopf bifurcation points, it was observed that there is a discontinuous jump, being more prominent on the left hand side of the bifurcation diagram. Due to this discontinuous growth in Q s , the analysis was conducted in a reverse direction across the β-range, which disclosed a region of hysteresis where one could find two spiral wave solutions for a single set of parameters. It made it more clear about the existence of subcritical bifurcations within FHN model.
The hysteresis region is the region of bistability where there exists two types of solutions for all the values of β within the region. We investigated this region by applying the technique of single shock defibrillation. We showed that only a rigidly rotating spiral wave can be converted into a meandering spiral wave but not conversely. Reasons for this could be related to filament tension, as noted in the conversion section above.
From theory, it is known that for subcritical Hopf bifurcations, as well as the boundary of the hysteresis region there is also a locus of unstable limit cycles. This locus extends from one end of the lower branch of the hysteresis loop to the opposite end of the upper branch. In all our studies, we were unable to simulate any solutions along this locus of unstable limit cycles. The absence of this locus should not distract from the result that a hysteresis region is present. Other studies involving FHN-like models and subcritical Hopf bifurcation have not reported this locus too, which asks the question about whether this locus is possible to simulate in FHN-like models for spiral waves [33]. Possible reasons for this are that, firstly, FHN is notoriously stable [33] and therefore simulating unstable solutions will either be extremely difficult or impossible. Also, the numerical techniques that we employed into simulating spiral waves in a comoving FOR relies on numerical stability in the calculation of the quotient data, particularly ω. If stability limits are breached, then the calculations of the quotient data is halted, and the simulation comes out of the comoving FOR. Hence, it may not be possible for our techniques to simulate unstable limit cycles. This is an area of research that we intend to look into for future work.
Furthermore, the β-step in the system was arbitrarily closely approximated by other β-steps with distinguishably different bifurcation points, resulting in wider hysteresis region. In other words, we can say that a small change to initial conditions may lead to different behavior, thereby showing the sensitive dependence on initial conditions. We have also demonstrated the numerical convergence analysis which confirms the true representation of the solutions. It highlighted that the overall shape of the bifurcation diagrams remains the same for all numerical and model parameters. It was clearly noted that the change in the numerical or model parameters does not qualitatively affect the bifurcation diagram depicting the growth of quotient size. Adding to it, these tests can be used in future to check for different numerical parameters which would enable us to run faster simulations for various studies.
In 2018, Fu et al published work relating to subcritical Hopf bifurcations in a FHN-like model. Although based on FHN, their model was modified significantly to achieve the results they needed. They noted that the FHN was very stable and we believe that their modifications enabled them to observe subcritical Hopf bifurcations by creating instability within the model. Hence, we believe that our results is novel in that it is the first report instance of subcritical Hopf bifurcation in the original FHN model. However, it is encouraging that within FHN type models, subcritical Hopf bifurcations can be found. Whether this can be extended to other FHN type models has yet to be seen.
In terms of practical uses of the discovery of subcritical Hopf bifurcations in FHN, there are many applications across the biological, physical and chemical sciences. In cardiac dynamics, the presence of two solutions for the same parameter set has been observed before but in the rigidly rotating region of FHN [42]. They showed that there were two rigidly rotating spiral wave solutions for the same parameter set within the hysteresis. In our current study, the knowledge of the presence of two solutions being two different types of spiral wave motion (one rigidly rotating and the other meandering) can assist in studying the conversion of certain arrhythmia such as monomorphic VT, as simulated by rigidly rotating spiral waves, to those such as ventricular fibrillation (VF) or polymorphic ventricular tachycardia (VT), both represented by meandering spiral waves. Subcritical Hopf bifurcations are also important in areas such fluid dynamics in which the presence of subcritical Hopf bifurcations relates to the onset of turbulence in fluid flow.