Saddle-node of limit cycles in planar piecewise linear systems and applications

In this article, we prove the existence of a saddle-node bifurcation of limit cycles in continuous piecewise linear systems with three zones. The bifurcation arises from the perturbation of a non-generic situation, where there exists a linear center in the middle zone. We obtain an approximation of the relation between the parameters of the system, such that the saddle-node bifurcation takes place, as well as of the period and amplitude of the non-hyperbolic limit cycle that bifurcates. We consider two applications, first a piecewise linear version of the FitzHugh-Nagumo neuron model of spike generation and second an electronic circuit, the memristor oscillator.


1.
Introduction. Bifurcations are qualitative changes in the phase portrait of families of differential equations as the parameter varies. The simplest cases are those involving equilibria. The next level of complexity would be bifurcations where limit cycles are involved. In particular, in a planar saddle-node bifurcation of limit cycles, also called fold of limit cycles, the phase portrait passes from exhibiting two hyperbolic limit cycles of different stability surrounding an equilibrium point, to the disappearance of such cycles through their collision in a non-hyperbolic semistable limit cycle.
This bifurcation is very common in applications. Specifically, the bistable configuration, where the equilibrium point in the interior of both limit cycles and the outer limit cycle are stable, whereas the inner limit cycle is unstable. In this configuration, the basin of attraction of both attractors are bounded by the unstable limit cycle, which can be considered as a threshold between the resting regime and the oscillatory one. These bistable systems are ubiquitous in biology [45], and saddle-node bifurcation is then a route for explaining some phenomena, for instance annihilation and single-pulse triggering [26]. In such phenomenon, an oscillatory behavior is ceased by injecting a sub-threshold pulse, and then, the activity is restarted by injecting a supra-threshold pulse. Annihilation has been described in several biological oscillators, such as the activity of the sinoatrial node [27], the eclosion rhythm of fruit flies, the circadian rhythm of bioluminescence in marine algae and biochemical oscillators, see [45]. Moreover, saddle-node bifurcation is also involved in the building of the elliptic bursting [30], a bursting mechanism (a type of oscillatory behavior of excitable systems whose main characteristic is an alternation of quiescent phases, and rapid oscillations phases [29]), which takes place in rodent trigeminal neurons [15]. Also in electronic circuits, phenomena involving saddle-node bifurcations can appear, see [1,43]. In particular, they are observed in the well-known Chua's circuit [23,24]. This phenomenon is of a great interest for the control system designer.
We can find very few works about the proof of the existence of saddle-node bifurcation of limit cycles in general nonlinear systems. In [12], it was shown that this bifurcation occurs in a generic two parameter unfolding of a homoclinic orbit with resonant eigenvalues. In [31,Th 3.6], authors prove the existence of a curve of saddle-node of limit cycles in canard regime in a family of slow-fast systems. However, it is not an easy task to give an expression, in the parameter space, where this bifurcation takes place. Alternatively, numerical methods can be used to analyze the saddle-node bifurcation of limit cycles. In fact, new methods and continuation packages have been developed with this goal [25].
From their appearance in the book of Andronov, Vitt and Khaikin [2], piecewise linear (PWL) systems have shown their capability not only to capture different behaviors coming from a wide class of applications [17,32], but also to reproduce a large amount of aspects of nonlinear dynamics. Furthermore, these systems show new behaviors, impossible to obtain under differentiability hypothesis [16,17,32]. Together with the property of mimicking the richness of nonlinear dynamics, PWL systems exhibit a simpler analytical treatment. This property allows, in many cases, to obtain quantitative information of the analyzed dynamical objects (for instance, the period and amplitude of limit cycles, [18,22,28]) and to explain the way some bifurcations take place [6].
Saddle-node bifurcation of limit cycles in the PWL context has been reported in different publications. In [35], the authors prove the existence of a codimension-1 manifold of saddle-node of limit cycles in a family of planar continuous PWL systems with three zones and symmetry coming from a heteroclinic connection. In [34], the authors study the number of limit cycles and prove the existence of two hyperbolic limit cycles with different stability surrounding one equilibrium point in a non-symmetric family of PWL systems. Even when this configuration is close to a saddle-node bifurcation as, for instance, we can observe in Figure 1 c) of [34], where two limit cycles are close one another, the saddle-node bifurcation is not reported in this paper. In [40], boundary equilibrium bifurcations in planar continuous PWL systems with two and three zones are studied. They find situations with two limit cycles and in this case the authors conjecture the existence of a saddlenode bifurcation of limit cycles. Regarding three dimensional PWL systems with two zones, in [23,24], the authors analyze the existence of a saddle-node of limit cycles as a degeneration of a focus-center-limit cycle bifurcation. Furthermore, they apply these results to the Chua's circuit. In [6,10], the authors describe a noose bifurcation in a PWL version of the Michelson system; such structure involves a saddle-node bifurcation, which is also analyzed. We finish this literature revision with [7,8], where a generalization of the Melnikov theory to non-smooth systems was developed and applied to prove the existence of saddle-node bifurcations of limit cycles in discontinuous and hybrid PWL systems.
In this article, we focus our attention on the proof of the existence of a saddlenode bifurcation of limit cycles in continuous PWL systems with three zones. A branch of saddle-node of limit cycles is proved to start at a degeneration of a focuscenter-limit cycle bifurcation. This conclusion is established through the application of the Implicit Function Theorem to the closing equations together with a nonhyperbolicity condition. This technique has been previously used to prove the existence of limit cycles in [3,5,22,23,36] and global connections in [9], but we remark that here it is applied not only to the set of closing equations, but also to the non-hyperbolicity condition. The starting point for applying the Implicit Function Theorem is, often, in the perturbation of a linear center, either in the plane [3,7,21], or in the space [4,6,22]. Moreover, the piecewise linear perturbation can be continuous [6,21,22], or discontinuous [3,7].
We consider two applications of our theoretical result. First, the McKean model, a PWL version of the FitzHugh-Nagumo system [20,38]. As far as we are concerned, the existence of the saddle-node bifurcation of limit cycles in the original differentiable FitzHugh-Nagumo system has not been proved yet, although there exist numerical evidences of its existence [41,42]. In [44], the author consider a discontinuous version of the McKean model with two zones and proves the existence of a saddle-node bifurcation of limit cycles. In the present paper, we consider a continuous version of the McKean model with four zones of linearity, three zones to mimic the cubic, and one small extra zone in one of the folds in order to capture the passing of the solutions through the fold, see [19]. From our main result we conclude the existence in this model of a saddle-node bifurcation and derive explicit expressions for both the curve of saddle-nodes and period of the saddle-node limit cycle. The obtained result is compatible with those in [44]. As a second application, we consider the memristor oscillator [14]. For this electronic circuit, we analyze the version of the model that was considered in [34]. Although in [34] they find cases where two limit cycles exist, they do not focus their attention on the existence of the saddle-node bifurcation of limit cycles.
The paper is outlined as follows. In Section 2, we introduce the target system and the main result. After that, Section 3 is devoted to the application of our result to first a PWL version of the FitzHugh-Nagumo neuron model and second the memristor oscillator. Subsequently, in Section 4 we include the proofs of the result established in Section 2. In Section 5 we state some conclusions and perspectives. Finally, we include two appendices. Appendix A is devoted to the most technical details of the proof of the main result and in Appendix B, we include an algorithm for fine tuning of the external impulse in order to facilitate annihilation/regeneration in a voltage trace of the McKean model. 2. PWL system with three zones. Main result. We focus our attention on the continuous planar piecewise linear (PWL) system with three zones of linearity (1) where u = (u 1 , u 2 ) T , the dot denotes the derivative with respect to the variable s and the piecewise linear vector field F is given by being v, w ∈ R, v < w, n L , n C , n R ∈ R 2 and M L , M C , M R 2 × 2 real matrices. The existence and uniqueness of solutions for the initial value problem associated to system (1) comes from the fact that F is a Lipschitz function. Note that, inter alia, matrices M L , M C and M R have to share their second columns due to the hypothesis of continuity of the vector field F . This shared column will be denoted by (m 12 , m 22 ) T , where the superscript T denotes the transposed. When m 12 = 0, the variable u 1 is decoupled and the dynamical behavior of the system (1) is not strictly bidimensional. We say, following [4], that system (1) is not observable. When the system is observable, an adequate change of variable allows to write the system into the canonical Lienard form. We enunciate this fact in the following result where we also transform the values v and w into −1 and 1, respectively. The proof of this result is direct. Proposition 1. Suppose that m 12 = 0. Then, the change of variables where Our interest begins under the assumption that system (3) possesses a center configuration in the central zone |X| ≤ 1. This fact implies D C > 0 and the existence of a unique equilibrium point (X C ,Ȳ C ) for the linear system (Ẋ,Ẏ ) T = N C (X, Y ) T + d T C . Now, we establish a result, whose proof is straightforward, where one can get, with a suitable change of variables, D C = 1 andȲ C = 0.
Proposition 2. Assume that D C > 0. Then, the change of variables where the prime denotes the derivative with respect to t, being We remark that system (4) can be written in the form where To begin with, in the following result, whose proof is direct, we impose the conditions for the existence of a continuum of periodic orbits of system (7) with two tangency points. Proposition 3. If a C = m = 0, then system (7) has a continuum of periodic orbits in the central zone. Moreover, the most external orbit Γ 0 of this continuum has two tangency points with the separation lines x = −1 and x = 1, at the points, q 0 = (−1, 0) and q 1 = (1, 0), respectively.
By perturbing this non-generic situation, we will look for non-hyperbolic periodic orbits Γ, which lives in the three regions L, C and R, coming from the most external orbit of the continuum Γ 0 . In particular, in Section 4 we prove the following result, which is the main result of this work, about the existence of a saddle-node bifurcation of limit cycles.
Theorem 2.1. Consider system (7) with fixed values of the parameters t L , t R , d L , and d R . If t L · t R < 0 and t L + t R = 0, then there exists a function a * C (m), analytic as function of m 1/3 , and defined in the neighborhood of the origin satisfying with a * C (0) = 0, and such that, for 0 < |m| 1, system (7) has a saddle-node bifurcation of limit cycles when a C = a * C (m). Specifically, if t L t R (t 2 L −t 2 R )(a C −a * C (m)) < 0, then the system has two three-zonal limit cycles with opposite stability and close to the periodic orbit C (m)) > 0, then the system has no limit cycles close to Γ 0 ; and if a C = a * C (m), then the system has a unique three-zonal semi-stable limit cycle Γ close to Γ 0 .
Moreover, approximations of the function a * C (m), the period and amplitude of limit cycle Γ are given by, and where the amplitude A has been measured as the difference between the intersection of the orbit with the separation line x = −1 with a positive y − coordinate minus the intersection of the orbit with the separation line x = 1 with a negative y−coordinate.
Remark 1. Note that, with the definition of the amplitude considered in Theorem 2.1, the tangent orbit Γ 0 of the unperturbed case of Proposition 3 has amplitude equal to zero.
Remark 2. Condition t L · t R < 0 in Theorem 2.1 is a necessary condition to prove the existence of a saddle-node bifurcation of limit cycles perturbing from the tangent orbit Γ 0 , see expression ofτ L in (22). Hence, we conclude that this bifurcation is only possible under the condition of having f (x; a C , t L , m, t R ) a quadratic shape. This conclusion does not contradict the existence of saddle-node limit cycles when f (x; a C , t L , m, t R ) has a cubic shape, (see, for instance, [35]), provided that in such case the saddle-node limit cycle perturbs from a heteroclinic connection. [34] provide conditions about the traces and determinants of the coefficient matrices of the system, including the inequality t L · t R < 0, such that there exists at most one limit cycle, which precludes the presence of a saddle-node bifurcation of limit cycles. Note that conditions in [34] are not satisfied under the hypotheses of Theorem 2.1, as it is assumed that the trace of the coefficients matrix in central zone t C = −m is small enough in absolute value.

Remark 4.
With respect to the cases not included in Theorem 2.1, we can provide the following information. The symmetric case (t L = t R and a C = 0), has been studied in [22] and only the existence of a limit cycle that perturbs from the linear center appears. This does not mean that the saddle-node bifurcation can not be obtained from other configurations, (see, for instance, [8,35]) and, as far as we are concerned, the analysis of the saddle-node bifurcation in the symmetric case is not closed. The reversible case (t L = −t R , m = a C = 0) implies an unbounded center. If t L = −t R = 0, propositions 4 and 5, provide us the existence of a family of non-hyperbolic periodic orbits with m = a C = 0, corresponding to that of the unbounded center. In the case t L = −t R , m = 0 and/or a C = 0, there could be a saddle-node that do not emerge from the tangent, but instead arise from any threezonal periodic orbit. This study requires a different analysis, such as the Melnikov theory used in [3,8].
In Fig. 1 we plot an schematic representation of the bifurcation diagram in m, a C parameter plane in case t L t R (t 2 L −t 2 R ) < 0 for m and a C sufficiently small. Solid line represents the saddle-node curve a C = a * C (m). In the region a C > a * C (m), two limit cycles with opposite stability and close to Γ 0 exist and in the region a C < a * C (m) no limit cycles close to Γ 0 exist.
< 0 for m and a C sufficiently small. Solid line represents the saddlenode curve a C = a * C (m). In the region a C > a * C (m), two limit cycles with opposite stability and close to Γ 0 exist and in the region a C < a * C (m) no limit cycles close to Γ 0 exist.

3.
Applications. In this section we consider two applications of our theoretical result. First subsection is devoted to a PWL version of the FitzHugh-Nagumo model, the McKean model. Second subsection is focused on the memristor oscillator.
3.1. The McKean model. The McKean model is a simplified piecewise linear model of neuronal activity with regular firing, [37]. This model can be derived from the FitzHugh-Nagumo model just by considering a piecewise linear approximation of the cubic nullcline of the voltage V . Let us consider the differential system where C > 0 is the capacitance and the cubic nullcline of V is given here by the piecewise linear function 1/2, t r > 0 and β 2 < 1/C. Some differences with respect to the standard McKean model are remarkable. In first place, the cubic nullcline of the FitzHugh-Nagumo model is here approximated by a four-linear-segments polygonal curve. The idea is to mimic, locally, one of the quadratic folds by three linear pieces (with a close to flat central slope, namely t c ) instead of a corner. Occasionally, this approximation gives results fitting better with the ones exhibited by the FitzHugh-Nagumo model, see [19]. In second place, in the classical McKean model, the slope t r of the segment defined in the strip V ∈ (a/2 + δ, (a + 1)/2) is set equal to one, while in our version of the model it is an independent parameter, that can be changed.
By setting the values of the parameters C, β, w 0 , a and δ, in the following result we prove the existence, in the parameter space (t c , I), of a saddle-node bifurcation of limit cycles in the McKean model. Additionally, approximated expressions for the bifurcation curve and for the period and amplitude of the saddle-node limit cycle are provided. We recall that in [44] the author proves a similar result, but in a discontinuous version of the McKean model with two zones of linearity.
In Figure 2, we illustrate the passing through the saddle-node bifurcation in the case (βC + 1)(t r + 1)(t r − 2βC − 1) < 0. When the parameter I is smaller than the bifurcation value I * , then no limit cycles exist near the equilibrium point which, locally, is an attractor, see panel (a). In panel (b), the parameter is just on the bifurcation value, I = I * , and a non-hyperbolic limit cycle appears, so-called saddle-node limit cycle. This limit cycle is stable from the outside and unstable from the inside. After passing through the bifurcation value, two concentric limit cycles perturb from the non-hyperbolic one. The outer limit cycle is stable whereas the inner one is unstable, see panel (c). Even when the bifurcation takes place involving three zones, some configurations exhibiting by the system are originated at the bifurcation. In particular, when the parameter increases far from the bifurcation value, the size of the inner limit cycle decreases and it becomes a two zones limit cycle, whereas the size of the outer limit cycle increase and it becomes a four zones limit cycle. Panel (d) shows this configuration for the parameter I larger enough than I * .
Corollary 1. Consider system (11) with fixed values of the parameters C, β, w 0 , t r , a and δ. If (1 + βC)(t r − βC) > 0 and t r − 2βC − 1 = 0, then there exists a function I * (t c ), defined for t c in a neighborhood of the value t * c = βC and such that (t c − βC)(1 + t r )(t r − βC)(t r − 2βC − 1) > 0, which corresponds to a saddle-node of limit cycles of the system. Specifically, if (1 + βC)(1 + t r )(t r − βC)(t r − 2βC − 1)(I −I * (t c )) < 0, then the system exhibits two three-zonal limit cycles with opposite stability; if (1 + βC)(1 + t r )(t r − βC)(t r − 2βC − 1)(I − I * (t c )) > 0, then the system has no limit cycles; and if I = I * (t c ), then the system has a unique three-zonal semi-stable limit cycle. Approximations for function I * (t c ) and the period of the saddle-node limit cycle are given by and Proof. Following propositions 1 and 2, the change of coordinates and time where r = (1 − βt c )/C > 0 and  Figure 2. Saddle-node bifurcation in the McKean model (11) with C = 0.25, β = 0.5, w 0 = 0, a = 1, δ = 0.25, t c = 0.1 and t r = 0.8. According to Corollary 1, the bifurcation takes place at I * = 1.263 . . . In panel (a), the parameter I = 1.2 is smaller than the bifurcation value, then no limit cycles exit near the equilibrium point which is a local attractor. In panel (b) the parameter is just on the bifurcation, I = I * , and a non-hyperbolic limit cycle appears. This limit cycle is stable from the outside and unstable from the inside. In panel (c) the parameter I = 1.265 is greater than the bifurcation value and then two concentric limit cycles perturb from the non-hyperbolic one. The outer limit cycle is stable whereas the inner one is unstable. The limit cycles perturbing from the saddle-node limit cycle move away one each other as the parameter increases far from the perturbation value I * . In panel (d), for I = 1.3, the inner limit cycle becomes a two zonal limit cycle whereas the outer one becomes a four zones limit cycle.
transforms system (11) into the following one wheref Note that, restricted to the bands {x < −1}, {|x| ≤ 1} and {1 < x ≤ 1/2δ}, this system coincides with system (7), Moreover, as δ is fixed, the dynamics of the fourth zone {x ≥ 1/2δ} does not influence the analysis around the tangent periodic orbit that exists for I = 0 and t c = βC. Therefore, the rest of the proof is a consequence of Theorem 2.1.
Saddle-node of limit cycles exhibited by the McKean model can be used to explain the switching between resting and spiking activity by injecting an external impulse, see Figure 3. In fact, from the approximation of the bifurcation value I * (t c ), obtained in Corollary 1, we propose a first proof of concept to address the tuning of the external impulse I, in order to facilitate annihilation/regeneration in a voltage trace. The starting point is the assumption that the system is close to a saddle-node limit cycle. The proposed solution switches between the resting and spiking activity by injecting an external impulse during a fixed time window, depending on the parameters of the system. The algorithm is stated in Appendix B. In Figure 3, we illustrate the result of this algorithm, the oscillatory behavior is annihilated and restarted again just by injecting a single pulse.
3.2. The memristor oscillator. Memristors are two-terminal electronic passive devices where charge and electric flux are related through a nonlinear function, called the characteristic of the memristor. This device was firstly introduced by Chua [13]. These class of new generation oscillators have potential to model the behavior of synapse connections in neurons. In [34], the authors propose the following modification for the nonlinear flux-charge characteristic of the memristor oscillator appearing in [14], namely, The state equations of the system of the mathematical model of the memristor oscillator are given by where the constants a, b L , b R , G, u and v depend on the components of the circuit. Under the hypothesis v(G−a) > 0 and by means of the application of the changes of the variables stated in propositions 1 and 2, system (14) can be transformed into the form (4), where The next result is a direct consequence of Theorem 2.1.
Corollary 2. Consider system (14) with fixed values of the parameters G, b L , b R , and a. If (14) has a three-zonal saddle-node limit cycle when ) < 0, the system exhibits two three-zonal limit cycles with opposite stability; ) > 0, then the system has no limit cycles; and if u = u * (v), the system has a unique three-zonal semi-stable limit cycle. Approximations for the function u * and the period of the saddle-node limit cycle are given by and

4.
Proof of main result. This section is devoted to the proof of Theorem 2.1. By perturbing the non-generic situation described in Proposition 3, in a first subsection we will look for a non-hyperbolic periodic orbit Γ, which intersects the three regions L, C and R, coming from the most external orbit of the continuum Γ 0 . In a second subsection we will prove that it corresponds to a saddle-node bifurcation.

4.1.
Existence of non-hyperbolic periodic orbits. Let us begin by introducing some notation. For chosen parameters η = (a C , t L , m, t R , d L , d R ) and a point p ∈ R 2 , we denote by ϕ(t; p, η) = (x(t; p, η), y(t; p, η)) the solution of system (7) with initial condition ϕ(0; p, η) = p. The coordinates of ϕ(t; p, η) will be referred to as x i (t; p, η) and y i (t; p, η), with i ∈ {L, C, R}, depending on the region where the solution belongs to, for that value of t. Consider a point p 0 = (−1, y 0 ). Assume that it exists a flight time τ Cu > 0 such that x C (τ Cu ; p 0 , η) = 1 and x C (s; p 0 , η) ∈ (−1, 1) for all s ∈ (0, τ Cu ). In such a case, we can define the Poincaré half-map between the switching lines x = −1 and x = 1 at the point y 0 as P Cu (y 0 , η) = y C (τ Cu ; p 0 , η). Similarly, we can define the Poincaré half-map between the switching lines x = 1 and x = −1 at the point y 2 as P C d (y 2 , η) = y C (τ C d ; p 2 , η), where τ C d > 0 is the flight time and p 2 = (1, y 2 ). On the other hand, consider a point p 1 = (1, y 1 ). Assume that it exists a flight time τ R > 0 such that x R (τ R ; p 1 , η) = 1 and x R (s; p 1 , η) > 1 for all s ∈ (0, τ R ). In such a case, we can define the Poincaré half-map between the switching line x = 1 and itself at the point y 1 as P R (y 1 , η) = y R (τ R ; p 1 , η). Similarly, we can define the Poincaré half-map between the switching line x = −1 and itself at the point y 3 as P L (y 3 , η) = y L (τ L ; p 3 , η), where τ L > 0 is the flight time and p 3 = (−1, y 3 ).
At this point, the Poincaré map for system (7) can be defined.
A periodic orbit which visits the three regions must satisfy P (y 0 , η) = y 0 , (see Fig. 4), or equivalently, the following eight equations, Figure 4. Representation of a three-zonal periodic orbit of system (7).
x C (τ Cu ; (−1, y 0 ), η) = 1, y C (τ Cu ; (−1, y 0 ), η) = y 1 , together with the constrains τ L , τ R , τ Cu , τ Cd > 0 and x C (s; (−1, y 0 ), η) ∈ (−1, 1), for all s ∈ (0, τ Cu ), x R (s; (1, y 1 ), η) > 1, for all s ∈ (0, τ R ), x C (s; (1, y 2 ), η) ∈ (−1, 1), for all s ∈ (0, τ Cd ), must be satisfied. Equations (16) and inequalities (17) are the conditions of the existence of a periodic orbit living in the three regions. To take into account the non-hyperbolicity of the periodic orbit, we consider the derivative of the Poincaré map, which corresponds to the exponential of the integral of the divergence of the system along such a periodic orbit, see [11]. In the particular case of PWL systems, the integral of the divergence can be explicitly computed as the sum of the products of the traces and the flight times in each region of linearity, see [21]. Hence, for a fixed point y 0 of Poincaré map P . Therefore, the condition of non-hyperbolicity in our case reads, Now, consider fixed values for the parameters t L , t R , d L , and d R . The existence of a non-hyperbolic limit cycle arising from the last periodic orbit of the linear center reduces to the existence of a set of parameters a C , m ∈ R, τ L , τ R , τ Cu , τ Cd > 0 and real values y 0 , y 1 , y 2 , y 3 , so that they satisfy equations (16) and (18) and inequalities (17). We begin by looking for a solution for the equations and later we will check that the solution satisfies the required inequalities.
We would like to point out that we have used the symbolic manipulators Mathematica and Maxima in order to check the truthfulness of expressions in (22), obtained first by hand.
Remark 5. Note that we have chosen parameter τ R in order to apply the Implicit Function Theorem, but parameter τ L could be equally chosen to obtain a dual result.
Once we have proven the existence of solution of system (16) and we have found an approximation of the solution (22), we proceed to prove in the following result that these solutions correspond to non-hyperbolic periodic orbits of system (7), checking that they satisfy inequalities (17).
Proposition 5. If t L · t R < 0, the functions given in Proposition 4 correspond to a non-hyperbolic three-zonal periodic orbit of system (7).
Proof. To ensure that the functions given in Proposition 4 correspond to a nonhyperbolic three-zonal periodic orbit, we need to check that functionsτ L (τ R ),τ Cu (τ R ) andτ Cu (τ R ) satisfyτ L (τ R ) > 0,τ Cu (τ R ) > 0 andτ Cd (τ R ) > 0 for τ R > 0 sufficiently small and that these functions and the rest of the functions verify inequalities (17) for τ R > 0 and small. For this purpose, we will assume in what follows that τ R is strictly positive and sufficiently small.
A similar reasoning allows to prove that the function x R verifies that x R (s; (−1,ȳ 3 ),η) > 1 for s ∈ (0, τ R ) and the second inequality of (17) is true.

4.2.
Correspondence to a saddle-node bifurcation. To finish with the proof of Theorem 2.1, in this second subsection we are going to prove that the nonhyperbolic limit cycle whose existence has been proved in the previous subsection corresponds, indeed, to a saddle-node bifurcation when t L + t R = 0. To obtain this, we will prove that the nondegeneracy conditions on the Poincaré map hold, that is, the second derivative of the Poincaré map defined in (15) with respect to the initial condition and the derivative of the Poincaré map with respect to parameter a C are both different from zero, see [33,39].
The section is divided into two parts. In the first one, we compute the second derivative of the Poincaré map with respect to the initial condition and we will see that this derivative is nonzero. In the second one, we compute the derivative of the Poincaré map with respect to parameter a C and we test that this derivative does not vanish.
Computation of the second derivative of the Poincaré map with respect to the initial condition. Consider the non-hyperbolic periodic orbit ϕ(t; (−1,ȳ 0 ),η) given in Proposition 5. From the study in [11], it is easy to see that the second derivative of Poincaré map with respect to y 0 in the non-hyperbolic fixed pointȳ 0 is given by where the derivative of the flight times with respect to y 0 have to be evaluated in Denote M ij the component ij of a matrix M and F k i the component i of the vector field defined by system (7) in the corresponding zone k ∈ {L, C, R}. From [11], the derivative of the flight times with respect to the initial condition are given by , and dτ L dy 0 = dτ L dy 3 dy 3 dy 2 dy 2 dy 1 The substitution of expressions (26)- (27) in (25) allows us, after some straighforward but tedious computations, to ensure that the second derivative of the Poincaré map is different from zero if and only if Finally, the substitution of variables τ Cu , y 0 , y 1 , y 2 , τ Cd , y 3 , τ L , a C , m, τ R by functions developed in expressions (22), provides us the following condition Thus, since t L · t R < 0 and t L + t R = 0, the first term in (28) is different from zero. Therefore, the first nondegeneracy condition on the Poincaré map holds.
Computation of the derivative of the Poincaré map with respect to parameter a C .
To compute the derivative of the Poincaré map with respect to parameter a C in the neighborhood of the non-hyperbolic periodic orbit γη, we will use the following expression from [39], where being w 0 is the orientation of the curve, if x, y, ∈ R 2 , we define the wedge product x ∧ y = x 1 y 2 − y 1 x 2 , div(F ) is the divergence of the vector field F defined by system (4) and F a C denotes the derivative with respect to parameter a C . Note that, although this identity is originally only valid for smooth systems, it can be generalized to continuous piecewise smooth systems with similar reasoning as those done in [8]. In our case, w 0 = −1 and F (γη(0),η)) = 0, so we will compute only the factor Q(ȳ 0 ,η) of expression (29) given in (30). The integral along a periodic orbit will be divided into four parts, depending on the zone of linearity that the orbit is located. By using Lemma A.5 of Appendix, these four addends are given by: 1. Part of the orbit located between y 0 and y 1 , in central zone, 2. Part of the orbit located between y 1 and y 2 , in right zone, 3. Part of the orbit located between y 2 and y 3 , in central zone, 4. Part of the orbit located between y 3 and y 4 , in left zone,

SADDLE-NODE OF LIMIT CYCLES IN PLANAR PWLS AND APPLICATIONS 5293
Using Lemma A.5, some direct but tedious computations allow us to compute and And then, we obtain Q(ȳ 0 ,η) = S 1 + S 2 + S 3 + S 4 . The substitution of variables τ Cu , y 0 , y 1 , y 2 , τ Cd , y 3 , τ L , a C , m, τ R by functions developed in expressions (22) provides us the following expression Thus, if t R = 0, the first term in the development is different from zero and therefore, the proof of the existence of the saddle-node bifurcation is concluded.
Finally, we proceed to obtain relations (8)- (10). From expression ofm in (22), we can obtain the expression of τ R in terms of the parameterm, by inverting the series. Moreover, taking into account that m · t R (t 2 L − t 2 R ) > 0, τ R > 0. Then, substituting this expression into the expression ofā C in (22), we obtain the existence of function a * C (m) and the approximation given in (8). The expression of the period T γ given in (9) has been obtained by summing up the series of the first three expressions of (22) plus the expression of τ R obtained. The approximation of the amplitude has been obtained with the first order ofȳ 0 −ȳ 2 from (22), by substituting the expression of τ R in the resulting expression.

5.
Conclusions and perspectives. In this article, we have proven the existence of a saddle-node bifurcation of limit cycles perturbing from a local center in planar continuous PWL systems with three zones. Power series of the bifurcation manifold, the amplitude and the period of the saddle-node limit cycle are also provided. The obtained result is general in the sense that it characterizes the bifurcation in any continuous PWL system with three regions of linearity. In particular, this result follows up theorems 7, 8, 9 and 10 in [34] in the following sense: In statement c) of these theorems, the authors prove the existence of two limit cycles, in a certain region of the parameter space. The region is given in terms of the trace of the coefficient matrix in central zone, and of a value ε > 0. On the other hand, in Theorem 2.1 we obtain, for each m sufficiently small, a curve a C = a * C (m) that provides a boundary for the existence of two limit cycles, as we can see, for instance, in Figure 1. By using the injectivity of function a * C , for each a C small enough, it is possible to find a value m = m * (a C ) such that (m * (a C ), a C ) belongs to the curve.
The value |m * (a * C )| corresponds to the ε of theorems 7-10 (c) in [34], after the appropriate change of variables.
This main result can be applied to a wide variety of models whose oscillatory behavior is well-known, and where these oscillations are shown to born in a saddlenode bifurcation of limit cycles. In particular, we consider two different applications.
The first one is a PWL version of the FitzHugh-Nagumo system, the McKean system. For this model, the bifurcation value is written in terms of the natural parameter of the system, the applied current. Note that there are not restrictions on the conductance value C, so it can be as small as needed. Therefore, these results can be applied in the slow-fast regime. Nevertheless, the saddle-node limit cycle, and the two limit cycles which collide are far from the canard regime. The analysis of fold limit cycles in the canard regime is beyond the objective of this work and is part of an ongoing project.
The second application, is an electronic circuit, namely the memristor oscillator [14]. We consider the version of the model studied in [34]. Moreover, the bifurcation value considered here is the boundary of linearity zones of the flux-charge characteristic of the memristor.
In both applications, we provide approximated expressions for the bifurcation curve and the period of the saddle-node limit cycle.
the following variational problem with respect to the initial conditions, Notice that we have already substituted m = 0 in the corresponding matrix in the central zone A C , see (5). Solution of (31) is given by Now, it suffices to substitute τ = π in previous expressions to conclude the proof of the first statement. The proof of the second statement is completely analogous. Consider the first equality in the third statement. Function ∂y R ∂y1 (τ ; (1, 0), η 0 ) corresponds to the second component of the solution of the following variational problem with respect to the initial conditions, As we need to evaluate the solution in τ = 0, it is straightforward that ∂y R ∂y1 (0; (1, 0), η 0 ) = 1. Analogously, considering the corresponding variational problem for the solution in the left zone, it follows that ∂y L ∂y3 (0; (−1, 0), η 0 ) = 1. In the next lemma we compute the non-null components of the 7th colum of the Jacobian matrix. The components in the 3th and 7th row can be obtained from Lemma A.1.
Now, it suffices to substitute τ = π in previous expressions to conclude the proof of the first statement. The proof of the second statement is completely analogous. Finally, consider the third statement. As functions ∂y R ∂a C (τ ; (1, 0), η 0 ) and ∂y L ∂a C (τ ; (−1, 0), η 0 ,) are the second component of the solution of the corresponding variational problem with respect to parameter a C , whose initial condition is always (0, 0), and we want to evaluate them in τ = 0, it is obvious that ∂y R ∂a C (0; (1, 0), η 0 ) = ∂y L ∂a C (0; (−1, 0), η 0 ) = 0. In next lemma we compute the components of the 8th column of the Jacobian matrix. The components in the 3th and 7th row can be obtained from Lemma A.1.
Now, it suffices to substitute τ = π in previous expressions to conclude the proof of the first statement. The proof of the second statement can be analogously done. The third statement is analogous to the proof of the third statement of Lemma (A.3).
Finally, we present an auxiliary result that is used for the computation of the derivative of the Poincaré map with respect to parameter a C .
Lemma A.5. Consider the autonomous linear systeṁ where x ∈ R n , A ∈ M n , b ∈ R n , n ≥ 1 and λ ∈ R. If A + λI is regular, then, where I the identity matrix of dimension n.
Appendix B. Fine tuning algorithm. In the following steps, we provide an algorithm to fine tuning the external impulse I in the McKean model, in order to switch from an oscillatory behavior to a resting behavior (step 1 to step 7), and vice versa (step 8 to step 10).
1. Given an oscillatory voltage trace V (t), corresponding with a limit cycle of system (11), with known parameter values C, β, w 0 , a, δ, t c and t r , and assuming that this limit cycle is close to a saddle-node limit cycle, compute the impulse, I 1 = I * (C, β, w 0 , a, δ, t c , t r ), from expression (12). 2. Compute the equilibrium point (V 1 , w 1 ) of system (11) with I = I 1 . Let λ = µ ± iη be the eigenvalues of the Jacobian matrix of the vector field at (V 1 , w 1 ). 3. From the first equation of system (11), approximate the trace of w(t). 4. Look for valuesṼ = V (t) andw = w(t) such thatṼ = βw − w 0 withw < w 1 . 6. Compute the impulse I 2 = w 2 − t c V 2 − (δ − a/2)(t c + 1) such that system (11), with I = I 2 , has an equilibrium point at (V 2 , w 2 ). 7. At timet, apply impulse I = I 2 − I 1 during the time π/η, where η is the imaginary part of the eigenvalue given in the step 2. 8. Compute the values 9. Compute the impulse I 3 = w 3 − t c V 3 − (δ − a/2)(t c + 1) such that system (11), with I = I 3 , has an equilibrium point at (V 3 , w 3 ). 10. At any time after π/η, apply impulse I = I 3 − I 1 during the time π/η, where η is the imaginary part of the eigenvalue given in the step 2.