Mixed mode oscillations in the Bonhoeffer-van der Pol oscillator with weak periodic perturbation

Following the paper of K. Shimizu et al. (2011) we consider the Bonhoeffer-van der Pol oscillator with non-autonomous periodic perturbation. We show that the presence of mixed mode oscillations reported in that paper can be explained using the geometrical theory of singular perturbations. The considered model can be re-written as a 4-dimensional (locally 3-dimensional) autonomous system, which under certain conditions has a folded saddle-node singularity and additionally can be treated as a three time scale one.


Introduction
There exist a lot of different mechanisms which can create a mixed mode oscillations (MMOs) in multiple time scale system of differential equations. In the paper of Desroches et al. (2012) one may find an information about MMOs near a singular Hopf bifurcation, near a folded node etc. analyzed from the perspective of singular perturbation theory and blow up methods (see the book written by F. Dumortier (1993) for theoretical backgrounds). Most of the research in this area is made for autonomous systems, but there exist some results for non-autonomous systems as well (e.g. Szmolyan and Wechselberger 2004). In what follows we focus on the system presented in the paper of Shimizu et al. (2011) This system could be seen as one of the versions of the forced van der Pol oscillator. The great advantage of this model is that it represents a simple electric circuit, so experimental results are easily available. The small parameter > 0 is responsible for a singular perturbation, while B 1 is also assumed to be small and positive and controls a regular non-autonomous periodic perturbation. For B 1 = 0 system (1) is a two dimensional autonomous system with classical canard explosion (see Krupa and Szmolyan 2001) controlled by the parameter B 0 . The second parameter k 1 makes enough room for a generalized Hopf bifurcation (see the book of Yu. Kuznetsov 2004) in the planar system. This bifurcation (also known as Bautin bifurcation) takes place when classical Hopf bifurcation degenerates due to the vanishing first Lyapunov index. Using MATCONT, a continuation package for Matlab (see Dhooge et al. 2003) we may easily compute that Bautin bifurcation happens at k 1 ≈ 0.513 and B 0 ≈ 0.365. The diagram presented by Shimizu et al. (2011) shows that for k 1 = 0.9 a subcritical Hopf bifurcation is observed for B 0 ≈ 0.20543. Additional computations made for k 1 = 0.2 detect a supercritical bifurcation at B 0 ≈ 0.4945. These samples will be useful for the numerical computations presented in the next sections. Systems very similar to (1) were widely studied since almost the beginning of the previous century (see van der Pol 1920). In the papers by Guckenheimer et al. (2003) and Bold et al. (2003) the results of the very extended bifurcation analysis are presented together with the references to many works on the topic.
The full system (1) reveals a very rich dynamics, for instance in the paper by Shimizu et al. (2011) authors present a sequence of stable mixed mode oscillations, which seems to be different from standard Farey sequence reported e.g. by Petrov et al. (1992). The key difference is the fact that in the system (1) we can easily find a periodic solution with more than two large oscillations in a row. In fact, it seems to be very difficult to find a MMO with only one large oscillation. Moreover the regions of parameters where the periodic MMOs exist are relatively small. We may also observe chaotic oscillations and bursting phenomena.
The theoretical background for multi-scale non-autonomous systems is not yet very developed. One of the typical approaches is to consider a time variable as an independent one (see Szmolyan and Wechselberger 2004). Then we can apply the theory valid for autonomous system and get better intuition about the dynamics. In this paper we also propose to rewrite the system as an autonomous one, but (it may seem to be artificial on the beginning) we decide to introduce two new variables z = sin ω t and p = cos ω t. This approach allows us to reveal better analogy to the known tree-dimensional systems with a fold singularity and explain the presence of MMOs using the methods known from geometric singular perturbation theory.
We write down a four dimensional system, keeping in mind an additional equation z 2 + p 2 = 1, which makes it in fact a locally three dimensional: Remark 1 For the further analysis it is convenient to see the system (2) as a three dimensional system in the space (x, y, z), but due to the trigonometrical origin of z there are two different p values in each fixed point, thus two different values ofż (difference is only in the sign). That means that at each fixed (x, y, z) there intersect two trajectories -one following increasing values of z and another one following decreasing values. On the figures it is sometimes better to show the trajectories in (x, y, p) space, the small oscillations are better visible this way. For = 0 system (2) has a slow manifold given by a surface y = −x + x 3 which has folds at x = ± The generalized Bonhoeffer-van der Pol system investigated by Sekikawa et al.
is in fact even closer to (1)  In the following sections we start with a change of coordinates which move the right fold of the slow manifold to the origin and brings the system (2) as close as possible to the form of (3). Then we project the flow on a slow manifold and check what kind of folded singularities are present. We discover that for certain values of bifurcation parameter µ the system can have a folded saddle-node. Rescaling of the coordinates makes the dynamics near the fold more clear. We show, that analogously to the results of De Maesschalck et al. (preprint) it is possible to consider the flow near the fold as a perturbation of the certain integrable system. Skipping the other blow-up charts (or exit and entrance areas, according to Krupa et al. (2008)) we make some calculations to approximate a return mechanism. Our analysis is essentially based on the theory built by Krupa et al. (2008), so in order to keep our paper short we do not repeat the information given there.

Transition to a standard form
It is convenient to transform the variables of the system (2) in the way that the right fold of the slow manifold is moved to the origin and is placed in a local minimum (see De Maesschalck et al. (preprint)). Thus we make a following change of coordinates: If we additionally introduce a new parameter µ = − 3 The parameter B 1 seems to disappear from the system, but in fact it is hidden in the additional relation (Z − µ) 2 + P 2 = B 2 1 . The first two equations of (6) are very similar to the ones from the system investigated by De Maesschalck et al. (after the transformation). The parameter ω in the paper of Shimizu et al. (2011) is a main bifurcation parameter and based on that results we may assume that it belongs to the interval (0, 1). Numerical results indicate that fixing ω and taking µ as a bifurcation parameter is possible and moreover we can take ω small enough so the system (6) can be treated as a three time scale system. To be more precise we assume ω = O( ).
What is clearly completely different from the systems (3) and (4) is the return mechanism: in the system (3) w grows as long as µ is positive and z is small enough. In (6) both Z and P are in fact periodic functions. The domain of P is the interval [−B 1 , B 1 ], while the domain of the variable Z depends on the parameter µ: Z ∈ [−B 1 +µ, B 1 +µ]. Geometrically speaking we observe the trajectories "reflecting" from the ends of the domain. No change of direction with respect to Z or P is possible outside the "walls" of the domain. Another difference is a lack of stationary points for (6). If we solve the equations for a stationary points we obtain that P = 0 together with Z = µ which imply B 1 = 0 and no perturbation. Instead of a singular point we obtain a periodic solution. For the same reason there is no Hopf bifurcation as such.
On the figure 1 we present some typical solutions of (6) for = ω = 0.1, B 1 = 0.01, k 1 = 0.9 and changing B 0 . Note that µ = 0 corresponds to B 0 ≈ 0.231. On the figure 1(a) we see the starting periodic solution with a small amplitude in X and Y , it corresponds to a stable stationary point in (3). With the growing |µ| this attracting periodic orbit disappears and we observe attracting MMOs. On the next pictures we observe how the shape of mixed mode oscillations are changing with B 0 . First small oscillations after reflection make almost the same way back (figure 1(c)), then the change from small oscillations to large happen earlier and earlier until on the figure 1(e) we observe no more reflection, just several small oscillations before the jump. This situations correspond to the one studied by Krupa et al. (2008). In this paper there is an assumption, that the number of small oscillation should be not too big, otherwise the approximations are not valid anymore. On the last figure we can see a trajectory corresponding to a relaxation oscillation in the paper of Krupa et al. (2008). It was shown in the paper of Szmolyan and Wechselberger (2004) that this trajectory lies on the torus.
These pictures raise a natural question what is happening between the small cycle and the first observed MMO. In (3) there exist a parameter range where we may observe a small stable cycle born in the singular Hopf bifurcation (see Guckenheimer 2008), this cycle grows and undergoes multiple series of the period doubling cascades which lead to MMOs. In the paper by Petrov et al. (1992) one can find a bifurcation diagram with so-called isolas corresponding to different patterns of MMOs. In case of the system (6) instead of a small cycle we would expect a small torus. Near k 1 = 0.9 such a solution is difficult to detect. The reason is that as we stated before the planar Hopf bifurcation is subcritical for this value of k 1 . Thus planar small cycles are unstable and in the full system the possible tori most probably have saddle properties. So to prove our hypothesis numerically we change the value of k 1 to 0.2 (which corresponds to a supercritical Hopf bifurcation) and indeed for B 0 = 0.4945 (µ ≈ −0.0058 ) we observe the small torus presented on figure2.
Let us now try to understand the mechanisms which form mixed mode oscillations with small number of SAOs. On the figure 3 we present a planar projection of a slow manifold for the unperturbed system (6) ( = 0). It consists of left and right attracting parts (S a− 0 and S a+ 0 ) and in the middle there is a repelling part (S r 0 ). When some perturbation is present ( > 0) these three parts perturb respectively into S a− , S a+ and S r . Outside the fold regions the perturbed manifolds "inherits" the dynamics from the unperturbed system following the rules of Fenichel's theory. The neighborhood of the fold  requires special analysis: we will show the presence of fold singularities and use a rescaling of the variables to understand the dynamics near Z = 0.
In the paper by Krupa et al. (2008) (among other results) there is presented a very detailed analysis of the MMO's of the type 1 k (one large oscillation followed by k small). The intervals of the parameter such as 1 k solution exist and is attracting are approximated. It is done by constructing a two-dimensional Poincaré map defined on the plane v = 0. This plane contains the fold and intersects an unperturbed slow manifold transversally. Poincareé map is divided into four parts: fold area, return mechanism and transitions between fold area and return mechanism. The key simplification of the problem is the fact that the two dimensional map can be successfully approximated by a one dimensional map. The reason is that each trajectory jumping off the left fold is getting attracted to the right attracting part of the slow manifold (S a+ ), then following the perturbed flow reaches the right fold region, then jumps off again to follow a fast flow back to the left attracting branch (S a− ). Again, following the slow flow the trajectory gets exponentially close to the straight line C − which denote the intersection of S a+ and the plane v = 0. As a result it is possible to study a one dimensional Poincaré map C − → C − , which is a much easier task. The line C − contains a critical canard point which corresponds to the transversal intersection of the left attracting manifold S a− and the repelling manifold S r . Every trajectory starting from the point on C − with w coordinate smaller that w cr undergoes small oscillations, while if it starts on the other side it jumps to S a+ as described before following the fast flow first and then creating a large loop (we call it return mechanism). It is possible to define sectors of rotation RS k -intervals of C − such as if the trajectory starts in RS k it undergoes exactly k small oscillations before leaving the fold region. Precise analysis of the dynamics near the fold is only possible for w small enough. In that case the solutions can be approximated as perturbations of a certain integrable two dimensional system.
Our numerical and analytical observations indicate that it should be possible to perform an analogous reduction for the system (6). As it was mentioned before our current goal is just to indicate some similarities and differences as well as prepare a background for a further detailed analysis.

Folded singularities
Let us start our analysis with the examination of the folded singularities of the system (6). In order to do so we project system (6) on the slow manifold Y = √ 3X 2 − X 3 = f (X). De-singularized (see Desroches et al. 2012) system has the following form: Here () in the RHS denotes the rescaled time-derivative ()f (X). Time flow changes the directory for X such as f (X) < 0, so on the attracting branches of the manifold.
Solving the equations for the stationary points we obtain (0, ± B 2 1 − µ 2 , 0) on the left fold and ( 2 1)) on the right fold. P 0 denotes a value of P corresponding to Z = 2 3 √ 3 ( 2 3 k 1 − 1). For the considered values of the parameters: k 1 ∈ (0, 1), B 1 and µ small there are no singularities on the right fold. Assuming |µ| < B 1 > 0 we have two stationary points A 1 = (0, B 2 1 − µ 2 , 0) and A 2 = (0, − B 2 1 − µ 2 , 0) on the left fold. Corresponding eigenvalues are: In both points zero eigenvalue corresponds to P -direction. Assuming B 1 and µ small enough the first point is a node while the second one is a saddle in the (X, Z) space. For µ = ±B 1 or ω → 0 two singular points melt into a single one placed in the origin and it is a folded saddle-node (in (X, Z) space again) with the eigenvector corresponding to zero eigenvalue transversal to a fold. In case of the system (3) we have only one singular point and the saddle-node appears as folded saddle bifurcating into a folded node. However in the analyzed range of parameters singularity moves outside the assumed domain (w ∼ O( )). The same situation is present in the system (6) for |µ| > B 1 (which seems to be satisfied for the parameter range where we observe MMOs with a small number of SAOs). Krupa and Wechselberger (2010) investigate the folded saddle-node singularity in the domain wide enough to take into account the singularities but in that case different phenomena play a key role.

Fold region in the rescaled coordinates
In this section we investigate the fold area using a classical approach of geometrical theory singular perturbations. We make a rescaling which is in fact identical to considering a family chart of blown-up coordinates: After an additional time rescaling (we keep the() notation) the system (6) takes form For = 0 andZ = 0 (which is satisfied for B 1 and µ = 0) first two equations of the system (9) form an integrable system with a constant of motion given by The special solution corresponding to H(X,Ȳ ) = 0 separates the area of periodic solutions and open level curves. This solution also corresponds to a singular canard (singular limit of intersection of the manifolds S a− and S r ).
The integrable system is exactly the same as the one examined by De Maesschalck et al. (preprint). Moreover, if we assume ω = O( ) as mentioned before, we can treat the full four-dimensional system as a perturbation of the integrable two dimensional one. The singular canard trajectory perturbs into a strong canard lying on the intersection of perturbed attractive and repelling manifolds. To compute the corresponding value ofZ (we call itZ cr ) we may follow the calculations from the paper by De Maesschalck et al. (the only difference is a plus sign atZ in the second equation) and conclude that the critical value ofZ is given bȳ The idea behind this formula is as following:Z cr divides the area of almostperiodic orbits from non-periodic ones. Thus the orbit passing this critical point can be treated as a periodic one, but with an infinite period time. Small ω allows us to treatP andZ as constants, while small guarantees that the orbits could be approximated by the orbits of the integrable planar system. The computations are almost identical to those made by De Maesschalck et al., but their presentation would involve a lot of additional notations.
If we recall the original variables we obtain Z cr = − . Since there is a relationship between P and Z we can compute P cr = ± B 2 1 − (Z cr − µ) 2 . That means that there exist two canard trajectories (obviously if and µ are small enough both Z cr and P cr belong to the domain). We can observe this phenomena for instance on the figure 1(e): the area of jumps is followed by the small oscillations and then by jumps again. However in the system (6) this double canard occur due to the symmetry of the trigonometric functions, it could be interesting to construct a three-dimensional polynomial system with similar geometry.

Return mechanism
To obtain an approximation of the return mechanism let us go back to the projection on the slow manifold described by the system (7). For and ω small enough we may assume that the change of P during "fast time" intervals can be neglected, so the slow manifold projection provides us enough information to approximate the map. If we divide the second equation of (7) by the first one we obtain We assume that Z is small comparing to X − k 1 f (X) which is justified by the domain where Z is defined. Then (recalling the relationship (Z−µ) 2 +P 2 = B 2 1 ) we may rewrite (10) as where the sign at the square root depends on the direction the trajectory is following. It is + for increasing P and − for decreasing. Let us introduce an auxiliary variable W = ω t. To cover all possible situations it is sufficient to take W ∈ [0, 3π]. We can easily compute that dW = dP ± √ B 2 1 −P 2 , which transforms (11) into a simpler dependence It is clear from the figure 3 that for the full return loop the variable X needs to run from √ 3 to 2 3 √ 3 (right attractive manifold) and then from − √ 3 3 to 0 (left one).
If we denote by P 0 a start point on S a− close to a jump area then a corresponding W 0 is assumed to be given by arccos P 0 B 1 if in P 0 P is a decreasing function of time and 2π − arccos P B 1 otherwise. We may approximate the result of the return mechanism as Now, depending in which part of the domain we landed we have Checking all possible combinations (the trajectory can come back without or with hitting the domain boarder) we obtain the formula where the sign is − if we follow the direction of increasing P and + otherwise. Numerical experiments confirm this approximation. Let us underline that this formula covers also the case when there is a "reflexion" effect in the way. The trigonometric functions are responsible for the "compression" close to the border of the domain while in the paper by Krupa et al. (2008) the distance between initial and end points of the return mechanism is approximately constant.
The formula (13) shows that the number of LAOs per length of the domain depends only on the parameters k 1 and ω. This approximation is good enough for the Poincaré map, but unfortunately does not allow us to generate periodic orbits. It is however practically easier to catch a (numerically) periodic orbit for k 1 = 0.2 then k 1 = 0.9, one of the possible reason is a clearer picture due to less LAOs per period. An example of such an orbit is presented on the figure 4.

Discussion
In this paper we presented a locally three dimensional system which is equivalent to a Bonhoeffer-van der Pol model. Keeping the parameter region close to the examples from the paper of Shimizu et al. (2011) we discussed the similarity of that system to known models with folded saddle-node II singularity. Some additional assumption were made in order to separate the time scales and eliminate the effect of the stationary points. A brief analysis of the fold region and an approximation of the return mechanism were provided.
We observed several differences such as a presence of two canards per period and non-homogeneous distribution of the return mechanism. Basically, these effects are present due to the choice of the variables. On the other hand this choice not only draw a parallel to the known and well analyzed systems, but also could be useful for the construction of polynomial systems with interesting dynamics.
Further analysis of the system could be especially interesting if we recall that the original system (1) models a simple circuit, which can be use for experiments with the parameters. That implies, that this circuit can be a very simple physical model for the three time scale systems with a folded saddle node singularity.
In our research we concentrated on the indication of the parameter range where the particular type of dynamics takes place. It could be interesting to go outside this region and take a look on the nearby bifurcations based on the papers of Guckenheimer et al. (2003) and Bold et al. (2003). To underline how different the behavior of the trajectory could be, on the figure 5 we present the periodic MMO with so-called bursting phenomena present. This type of solution is typical for the neuronal models (see Izhikevich 2007). The parameters are like in case of the figure 1(e), except B 1 = 0.1. This change also cause that the stationary points on the fold exist even for relatively large |µ|. Another interesting phenomenon to study is presence of the delayed Hopf bifurcation (see Wechselberger 2010, Neishtadt 1987).