Mixed mode oscillations and phase locking in coupled FitzHugh-Nagumo model neurons
Elizabeth N. Davison Zahra Aminzare Biswadip Dey Naomi Ehrich Leonard







Mixed mode oscillations and phase locking 
 in coupled FitzHugh-Nagumo model neurons
    Elizabeth N. Davison
Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08540, USA.

Zahra Aminzare
 Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA.

Biswadip Dey^*

Naomi Ehrich Leonard ^*

    
=========================================================================================================================================================================================================================================================================================








We study the dynamics of a low-dimensional system of coupled model neurons as a step towards understanding the vastly complex network of neurons in the brain. We analyze the bifurcation structure of a system of two model neurons with unidirectional coupling as a function of two physiologically relevant parameters: the external current input only to the first neuron and the strength of the coupling from the first to the second neuron. Leveraging a timescale separation, we prove necessary conditions for multiple timescale phenomena observed in the coupled system, including canard solutions and mixed mode oscillations. For a larger network of model neurons, we present a sufficient condition for phase locking when external inputs are heterogeneous. Finally, we generalize our results to directed  trees of model neurons with heterogeneous inputs.



Key words. Mixed mode oscillations, coupled FitzHugh-Nagumo, geometric singular perturbation, bifurcation 



AMS subject classifications. 34A26, 34C15, 34C60, 34D15, 34E17, 37G05, 37G10, 92B20, 92C20 
Mixed mode oscillations and phase locking 
 in coupled FitzHugh-Nagumo model neurons
    Elizabeth N. Davison
Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08540, USA.

Zahra Aminzare
 Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA.

Biswadip Dey^*

Naomi Ehrich Leonard ^*

    
=========================================================================================================================================================================================================================================================================================







§ INTRODUCTION




The study of model neurons has a rich history, dating back to the pioneering work of Hodgkin and Huxley <cit.> on the action potential in the squid giant axon. A two-dimensional model that captures salient qualities of the four-dimensional Hodgkin-Huxley model was developed independently by FitzHugh <cit.> and by Nagumo<cit.>. In this model, commonly known as the FitzHugh-Nagumo (FN) model, one variable represents the membrane potential and the other represents a gating variable. A constant external input to the FN model neuron can produce quiescent behavior (a low-voltage stable equilibrium point), firing (a stable limit cycle), or saturated behavior (a high-voltage stable equilibrium point). The FN model neuron captures realistic neuronal behavior such as spike accommodation, bistability, and excitability <cit.>. 

A system of two coupled neurons can represent a larger network of neurons that cluster into two groups in which neurons within each group synchronize but neurons in different groups do not.
 A cluster synchronized network can be reduced to a quotient network <cit.> by leveraging balanced conditions on coupling and graph structure <cit.>, as well as  
bounds on coupling strength 
<cit.>.
A system of two FN model neurons with gap junction diffusive coupling (two-FN system) has been studied numerically and analytically in the symmetric case<cit.>, where both neurons receive the same external input and are coupled bidirectionally (undirected coupling). Gap junction diffusive coupling is modeled as a difference between the membrane potentials of the two model neurons multiplied by a parameter that represents the coupling strength. The two-FN system has also been studied numerically<cit.> in a context where the intrinsic properties of both models are the same but the neurons are coupled unidirectionally (directed coupling). Here, we add to the existing literature by analytically describing the bifurcation structure of the directed two-FN system in terms of two parameters, the external input to the first model neuron and the unidirectional coupling strength from the first model neuron to the second.

The FN model neuron is a classic example of a fast-slow system, and the coupled pair of FN model neurons exhibits rich dynamics characterized by the timescale separation. Under certain conditions on external input and coupling strength, the system exhibits canard solutions, which are solutions that pass from a stable to an unstable manifold in the slow system and stay near the unstable manifold for a long time relative to the slow system timescale <cit.>. Canard solutions result from the presence of folded singularities, particularly folded nodes, which give rise to robust families of canard solutions <cit.>. Canard solutions can lead to mixed mode oscillations (MMOs), which are periodic oscillations that alternate between canard-driven oscillations and a relaxation oscillation <cit.>. The existence of canards and MMOs has been described for systems in four dimensions 
<cit.>, systems with two slow variables and two fast variables <cit.>, and generalized systems in arbitrary finite dimensions <cit.>. The folded saddle node of type I (FSN I) and folded saddle node of type II (FSN II) have been identified as mechanisms for the onset of MMOs in fast-slow systems <cit.>. We leverage these results to determine the parameter regimes where canards and MMOs may be present in the directed two-FN system, which has two slow variables and two fast variables.	
	
Canard-induced MMOs have been studied analytically in numerous systems including chemical reactions <cit.>, the Hodgkin-Huxley neuronal model <cit.>, cortical grid cells <cit.>, and a self-coupled FN model neuron <cit.>. In a two-FN system, the onset of firing, as coupling strength is increased, can be characterized by the appearance of canard solutions, and by MMOs, as the coupling is increased further. The existence of canard solutions in a two-FN system was first proven using nonstandard analysis in the case of model neurons with identical parameters <cit.>. Necessary conditions were found in terms of a model parameter that controls the slope of the linear nullcline of the system. Conditions for different stability types of folded singularities were found in terms of the same model parameter in a slightly modified, but still symmetric, model <cit.>. Here, we fix the corresponding parameter within the range where canard solutions may be present and find conditions for existence of canard-induced MMOs in terms of two parameters that break symmetry: external input and coupling strength.

A condition for the onset of MMOs in a two-FN system was shown as an application of a method developed to study MMOs in systems with two fast variables and two slow variables <cit.>. There are no symmetry requirements and the main result is a necessary condition for MMO onset in terms of a parameter corresponding to the input to one of the neurons. In the spirit of this work, we prove explicit necessary conditions for existence of canard solutions and MMOs in the directed two-FN system in terms of both the external input and the coupling strength. First, we take the singular limit of the system, and obtain necessary conditions on the bifurcation parameters for existence of transcritical bifurcations. The transcritical bifurcations in the singularly perturbed system delineate regions in parameter space where MMOs exist in the original system.

We show, further, that the original system admits Hopf bifurcations within a distance of order ϵ around the point in the parameter space where the singularly perturbed system admits transcritical bifurcations. 
This we use to derive novel bounds for phase-locking in representative networks of model neurons.
Phase locking is a generalization of synchronization where the phases of oscillating models remain separated by a constant offset, while amplitudes and waveforms may vary <cit.>. A common phenomenon in nature, phase locking has been studied in cardiac rhythms <cit.>, in the firing patterns of squid axons <cit.>, in two coupled phase oscillators <cit.>, in local field potential measurements of neurons in the human brain <cit.>, and in the brain as a mechanism for coordination between groups of neurons <cit.>. 


Finally, we consider the more general problem of n FN model neurons linked by unidirectional gap junction diffusive coupling in a directed tree, with heterogeneous coupling strengths and heterogeneous external inputs. As in the directed two-FN system, this can represent a class of large networks that contain cluster synchronized groups of model neurons and satisfy conditions on graph structure <cit.> and connectivity<cit.> so they can be reduced to a quotient network <cit.>.  An analogous problem with homogeneous coupling strength has been analyzed in detail in the strong coupling limit where the dynamics are reduced using singular perturbation theory<cit.>.   Here, we leverage an analysis of the singular perturbation of the directed two-FN system to provide necessary conditions for the existence of MMOs and sufficient conditions for phase locking in the original n-FN system.

Our contributions towards understanding the dynamics of networked nonlinear model neurons are as follows. First, we explain how the bifurcation structure of the directed two-FN system relates to the bifurcation structure of the reduced, singularly perturbed system that is used to study canard solutions. This is critical because the reduced system can be used to explain features of the original system and the original system can be used to understand the reduced system. Second, we provide necessary conditions for canards and MMOs in the directed two-FN system in terms of two model parameters; this is an extension of the conditions found in terms of one parameter in the literature. Third, we provide a sufficient condition for phase locking given heterogeneous external inputs in the directed two-FN system. We generalize these conditions to directed trees of FN model neurons.

The paper is organized as follows. In Section <ref>, we review the standard analysis of a single FN model neuron
 and give a biophysical rationale for bounds on model parameters used throughout the paper. In Section <ref>, we define the directed two-FN system and find conditions for Hopf bifurcations. In Section <ref>, we compute the singular perturbation of the directed two-FN system. We prove necessary conditions for the transcritical bifurcations in the singularly perturbed system and canards and MMOs in the original system in Section <ref>. In Section <ref>, we generalize the results to directed trees of FN model neurons. We provide a numerical example to illustrate our results.










§ SINGLE FITZHUGH-NAGUMO MODEL NEURON




The FN model is a two-dimensional reduction of the four-dimensional Hodgkin-Huxley (HH) model for the neuronal action potential. By letting y and z represent the membrane potential and a slow gating variable, respectively, its dynamics are given by

    dydt   = ψ(y) -z + I,
    
    dzdt   = ϵ (y-b z),

where ψ(y) is a cubic polynomial. For our purposes, we use ψ(y) = y-y^3/3-a. In this model, I corresponds to an external input, 0<ϵ≪1 is a positive timescale separation constant, and a and b are positive constants. The FN model is far simpler to analyze than the full HH model due to the lower dimension.

Despite the lower dimension, the FN model captures a range of physiologically meaningful regimes and behaviors <cit.>. The bifurcation structure of a single FN model is complex, with a Hopf bifurcation and canard explosion, or transition from small, Hopf-induced limit cycles to large relaxation oscillations<cit.>, marking both the transition from quiescent to firing and the transition from firing to saturated. As I is increased, the single stable equilibrium point transitions to relaxation oscillations, or stable limit cycles. As I is further increased, the large stable limit cycles persist and small unstable limit cycles appear. The unstable limit cycles disappear in a subcritical Hopf bifurcation at I_0 and the equilibrium point becomes unstable. Figure <ref> depicts the bifurcation diagram of the FN model when I is varied. 



The FN model is a suitable choice for network analysis because it is both dynamically rich and analytically tractable.  A necessary condition<cit.> for the FN model to exhibit distinct quiescent, firing, and saturated regimes is the existence of a unique equilibrium point for all values of the bifurcation parameter I. In this paper we assume the following:

Parameters a, b, and ϵ are such that the FN model (<ref>) has a unique equilibrium point for all values of I≥ 0. This results in conditions 0<a<1 and 0<b<1. 

The condition on a corresponds to a simple voltage offset requirement and the condition on b corresponds to a requirement that the slope of the linear nullcline of (<ref>) must be greater than that of the cubic nullcline of (<ref>).

The following proposition from Ref. <cit.>  
describes the stability of the unique equilibrium point of (<ref>) given Assumption <ref> and conditions on I for Hopf bifurcations.

 (Ref. <cit.>)
Let Assumption <ref> hold. Then, there exists I_0<I_1 such that the equilibrium point is stable for I<I_0 and, as I increases, it will undergo a transition to an unstable equilibrium point through a Hopf bifurcation at I_0. As I is increased further it will undergo a transition from unstable to stable through a second Hopf bifurcation at I_1. 


 
We review the approach to analyzing stability of the limit cycles arising from the Hopf bifurcations following  the methods of Chapter 3 of Ref. <cit.> and the application to the FN model in Ref. <cit.>.  We  generalize the approach to networks of FN model neurons in later sections.

The dynamics at a Hopf bifurcation at the origin of a two-dimensional system can be written as

    ([ dx_1dt; dx_2dt ]) = ([  0 -ω;  ω  0 ])([ x_1; x_2 ])+([ F(x_1, x_2);  G(x_1,x_2) ]),

such that F and G satisfy F(0,0) = G(0,0) = 0 and D_𝐱 F(0,0) = D_𝐱G(0,0)=0 where 𝐱 = [x_1,   x_2].



 Consider the system 𝐱̇ = 𝐟(𝐱,γ).
The coefficient of the cubic term of the Taylor expansion of 𝐟, (<ref>) is expressed as

    α   = 1/16(F_x_1 x_1 x_1+F_x_1 x_2 x_2+G_x_1 x_1 x_2+G_x_2 x_2 x_2)|_(0,0)
       +1/16ω(F_x_1 x_2(F_x_1 x_1+F_x_2 x_2)
        -G_x_1 x_2(G_x_1 x_1+G_x_2 x_2)
       -F_x_1 x_1G_x_1 x_1+F_x_2 x_2G_x_2 x_2)|_(0,0),
 
where F_x_1 x_2 denotes ∂^2 F∂ x_1 ∂ x_2, and so on. 






The system 𝐱̇ = 𝐟(𝐱,γ), admits a Hopf bifurcation for the parameter value δ_0 at an equilibrium point 𝐩 if





  * D_𝐱𝐟(𝐩,δ_0) has a pair of pure imaginary eigenvalues and no other eigenvalues with zero real parts.


  * ∂∂δ(λ(δ))|_δ = δ_0≠ 0, where (λ) denotes the real part of the eigenvalue λ.


  * The cubic coefficient of the Taylor expansion of 𝐟, denoted by α and defined in Definition <ref>, is nonzero.


Furthermore, if α<0, the Hopf bifurcation is supercritical, while, if α>0, the Hopf bifurcation is subcritical.





The cubic coefficient is also called the first Lyapunov coefficient. It can be used to determine the location of a solution that separates small unstable limit cycles from larger stable limit cycles and so delineates the canard explosion <cit.>. The value of the bifurcation parameter I at the canard explosion
 is close to the value of I at the Hopf bifurcation, 
 as illustrated in Figure <ref>. 

For the FN model, the cubic coefficient is given by

    α   =18(2b-b^2 ϵ-11-b^2ϵ).

In this paper we choose parameters that ensure Assumption <ref> holds and the bifurcations are subcritical Hopf (α>0); these yield biologically realistic dynamics <cit.>. We fix a=0.875, b=0.8, and ϵ = 0.08, and we consider I≥ 0 as a bifurcation parameter.











§ DIRECTED TWO-FN MODEL NEURON SYSTEM




The directed two-FN model neuron system is shown in Figure <ref>. The first model neuron is denoted A and it receives external input I. The second model neuron is denoted B and it receives no external input.  The coupling is unidirectional from A to B, with coupling strength γ. A and B have the same intrinsic dynamics, i.e., the same values of a, b, and ϵ as defined above. We let I and γ be bifurcation parameters. The directed two-FN system corresponds to the reduced cluster synchronized setting where there are many more model neurons in one cluster (the first) than the other (the second).


The equations for the directed two-FN system are


    dy_Adt   = ζ_A(y_A,z_A, y_B, z_B)= y_A-y_A^3/3-a -z_A+I,
    dz_Adt   = ϵξ_A(y_A,z_A, y_B, z_B)= ϵ (y_A-b z_A),
    dy_Bdt   =ζ_B (y_A,z_A, y_B, z_B)= y_B-y_B^3/3 -a -z_B+γ (y_A-y_B),
    dz_Bdt   = ϵξ_B(y_A,z_A, y_B, z_B)= ϵ (y_B-b z_B).


Here, y_A (y_B) is the membrane potential of A (B) and z_A (z_B) represents a slow gating variable in A (B).

The bifurcation structure of directed and undirected two-FN systems have been studied extensively from a numerical perspective <cit.>. The bifurcation structure of the undirected system has been studied through analytical methods that leverage symmetry-based arguments or assume symmetric or near-symmetric FN models <cit.>. In contrast we examine the system with asymmetry in both the external input and the coupling.

We begin by classifying the behavior of the two model neurons A and B in the I-γ parameter space, as shown in Figure <ref>. Regions (1)-(7) are described as follows:






  * For I less than the value where A undergoes the lower Hopf bifurcation, I<I_0A, both A and B will be quiescent at a stable equilibrium point. This corresponds to region (1) of Figure <ref>. See Section <ref> for details.



  * For I greater than the value where A undergoes the upper Hopf bifurcation,  I>I_1A, A becomes saturated and, as γ varies, there are three distinct behaviors  for B. See Section <ref> for  details. For I>I_1A, there exist I_0B(I) and I_1B(I) such that





  * For γ<I_0B(I), B is quiescent. This corresponds to region (2) in Figure <ref>.

  * For I_0B(I)<γ<I_1B(I), B is firing. This corresponds to region (3) in Figure <ref>.

  * For γ>I_1B(I), B is saturated. This corresponds to region (4) in Figure <ref>.




  * For I_0A<I<I_1A, A is firing and there are three distinct behaviors for B:





  * When γ>1-bϵ, B is phase locked with  A. This corresponds to region (5) of Figure <ref>. See Section <ref> for details.

  * When γ<1-bϵ and I is below a curve, denoted by I_*(γ), MMOs or small canard oscillations may be present. This corresponds to region (6) of Figure <ref>. See Section <ref> for the derivation and Section <ref> for details.

  * When γ<1-bϵ and I is above I_*(γ), B is firing. This corresponds to region (7) of Figure <ref>. See Section <ref> for the derivation and Section <ref> for details.








In Section <ref> we desingularize the system (<ref>) by performing a singular perturbation and changing timescales. We analyze the resulting second-order dynamics and draw conclusions about canards and MMOs for the original dynamics (<ref>) in Section <ref>.



§ FAST-SLOW PHENOMENA IN THE DIRECTED TWO-FN SYSTEM


In this section, we assume A is firing and study the onset of firing in B as γ increases. We begin by providing definitions of canards and MMOs, which are observed numerically at the transition from quiescent to firing in B. For a general fast-slow system expressed as

    d𝐲dt   = f(𝐲,𝐳),
    
    d𝐳dt   = ϵg(𝐲,𝐳),

𝐲∈R^m are fast variables, 𝐳∈R^n are slow variables, and 0<ϵ≪ 1 is the timescale separation. The singular limit corresponds to ϵ = 0.


	Given system (<ref>) with ϵ=0, 
	
    C= { (𝐲,𝐳)∈R^m ×R^n : f(𝐲,𝐳) = 0}

	 is called the critical manifold.




Given system (<ref>), denote the curve of points in C that are not normally hyperbolic as L = {(𝐲,𝐳)∈ C : D_𝐲f(𝐲,𝐳) = 0, D_𝐲𝐲f(𝐲,𝐳) ≠ 0}. L is called the fold curve. The singularities on the fold curve are called folded singularities.


 
A solution of (<ref>) is called a canard if it stays within 𝒪(ϵ) of a repelling branch of the critical manifold for a time that is 𝒪(1) on the slow timescale, τ_1=tϵ.



Periodic solutions of (<ref>) with peaks of substantially different amplitudes are called MMOs. Canard solutions often comprise the small oscillations present in MMOs. 


To find canards in system (<ref>), we first examine the dynamics of (<ref>) in the singular limit. A singular canard<cit.> is a solution that crosses from the attracting to repelling branch of the critical manifold through a folded singularity. In the singular limit, a singular canard predicts<cit.> a family of canard solutions in (<ref>). An example of a canard solution and MMOs found in the directed two-FN system is shown in Figure <ref>.



The remainder of this section details the calculations involved in the transformation of (<ref>) into a two-dimensional desingularized system, bifurcation and phase plane results for the desingularized system, and computation of the types of bifurcations. We implement methods that have been developed to determine the presence of canards and MMOs in general systems with m fast variables and n slow variables <cit.> and methods that have been applied to similar problems in two-FN systems <cit.>. In Section <ref>, we prove that solutions are canards during the onset of firing for sufficiently large coupling strength.



 §.§ The critical manifold and fold curve of the directed two-FN system


In the slow timescale τ_1=tϵ system (<ref>) becomes

    ϵ dy_Adτ_1   = y_A-y_A^3/3-a -z_A +I,
    
    ϵ dy_Bdτ_1   = y_B-y_B^3/3 -a- z_B+γ (y_A-y_B),
    
    dz_Adτ_1   = y_A-b z_A,
    
    dz_Bdτ_1   = y_B-b z_B.

The critical manifold, C, can be found by setting ϵ = 0 in the fast two equations:

    C    = {y_A, y_B, z_A, z_B
    |ζ_A(y_A, y_B, z_A, z_B) = 0, ζ_B(y_A, y_B, z_A, z_B) = 0
    . },
	
where ζ_A and ζ_B are defined in (<ref>).

We denote by S_0 a compact submanifold of the smooth manifold C. S_0 is “normally hyperbolic" if the eigenvalues of the Jacobian evaluated at each point p ∈ S_0 are nonzero. We can select a compact submanifold of C that does not contain the singular points and so will be normally hyperbolic  <cit.>. 

By Fenichel's Theorem, the slow dynamics of the two-FN system (<ref>) will lie 𝒪(ϵ) away from S_0, on a normally hyperbolic slow invariant manifold S_ϵ with the same stability properties as S_0. <cit.>

The fold curve is the set of singular points of C as defined in Definition <ref>, or the points where the determinant of the Jacobian with respect to y_A and y_B is zero. The fold curve for (<ref>) satisfies


y_A-y_A^3/3-a -z_A +I    = 0,


y_B-y_B^3/3 -a- z_B+γ(y_A-y_B)    = 0,


(1-y_A^2)(1-y_B^2-γ)    = 0.





 §.§ The singular limit of the directed two-FN system

In the singular limit of (<ref>) when ϵ = 0, the reduced two-FN system is 


    0    = y_A-y_A^3/3-a -z_A +I, 
    
    0    = y_B-y_B^3/3 -a- z_B+γ (y_A-y_B), 
    d z_Adτ_1   = y_A-b z_A, 
    d z_Bdτ_1   = y_B-b z_B.



Following Chapter 19 in Ref. <cit.>, we differentiate (<ref>) with respect to τ_1 to get

    (D_𝐲ζ) ·d 𝐲dτ_1
    +
    (D_𝐳ζ) ·d 𝐳dτ_1
    =
    0,

where ζ = ([ ζ_A; ζ_B ]), 𝐲 = ([ y_A; y_B ]) and 𝐳 = ([ z_A; z_B ]). Using (<ref>)-(<ref>), we can rewrite (<ref>) as

    (D_𝐲ζ)d 𝐲dτ_1   = 
    - adj(D_𝐲ζ) D_𝐳ζ· (𝐲-b𝐳),

where

    D_𝐲ζ = ([   1-y_A^2         0;         γ 1-y_B^2-γ ]),  
    D_𝐳ζ = - 𝐈_2×2,

 and 𝐈_2×2 is the identity matrix. The determinant of D_𝐲ζ is zero at y_A* = ± 1 and y_B*= ±√(1-γ).
 
 The two-dimensional reduced system (<ref>) can be examined, as follows, by rescaling time to obtain the desingularized system described in Ref. <cit.>. 




 §.§ Desingularization of the directed two-FN system




Rescaling time in (<ref>) by τ_1 = -(D_𝐲ζ)τ_2 yields the desingularized system:

    d y_Adτ_2   = ρ_1(y_A, y_B) 
    = (1-y_B^2-γ)(y_A -bz_A ),
    
    d y_Bdτ_2   = ρ_2(y_A,y_B) 
    = -γ(y_A -bz_A )+(1-y_A^2)(y_B-bz_B),

where z_A and z_B are defined as

    z_A   = y_A-y_A ^3/3 -a+I,
    
    
    z_B    = y_B-y_B^3/3-a+γ (y_A -y_B).


The singularities of (<ref>) are the points where adj (D_𝐲ζ)(D_𝐳ζ) · (𝐲-b𝐳)=0. There are two types of singularities <cit.> that are present in (<ref>). The first is denoted an ordinary singularity, which arises when (D_𝐲ζ)≠ 0.
The second type of singularity is the folded singularity, which arises when (D_𝐲ζ)=0.

The ordinary singularities of (<ref>) correspond to the equilibrium of (<ref>), which is the equilibrium of the original system (<ref>) in the 𝐲 plane. In Figure <ref>, the signs of the eigenvalues associated with the ordinary singularity are depicted as a function of I and γ. These show stability changes in (<ref>) that correspond to Hopf bifurcations in (<ref>). The curves in I-γ space that delineate the different signs of the real parts of the eigenvalues are slightly different from the Hopf bifurcation curves because predictions from singular perturbation analysis are accurate up to 𝒪(ϵ). In regions where A and B are both quiescent or saturated (dark gray), both eigenvalues have negative real parts and the ordinary singularity is a stable equilibrium. In regions where both A and B are firing (light blue), the real parts of both eigenvalues are positive and the ordinary singularity is an unstable equilibrium. In regions where either A or B is firing (light gray), the ordinary singularity is a saddle.






The folded singularities in (<ref>)  correspond to singularities on the fold curve in (<ref>), and they satisfy 

(D_𝐲 ζ) = 0,
  and   adj (D_𝐲 ζ)(D_𝐳 ζ) ·(𝐲-b𝐳) = 0,

equivalently



    (1-y_B*^2-γ)(1-y_A*^2)    = 0,
    
    
    (1-y_B*^2-γ)(y_A* -bz_A* )    = 0,
    
      
    (1-y_A*^2)(y_B*-bz_B*) -γ(y_A* -bz_A* )    = 0.


(<ref>) is satisfied when y_A* = ± 1 or  y_B* = ±√(1-γ).

When y_A* = ± 1, solving (<ref>)-(<ref>) results in y_A* = bz_A*, or γ = 0 and y_B* = ± 1. Combining y_A* = bz_A* with (<ref>) results in I = ±1/b∓2/3+a, which is independent of γ and so we do not consider this case further.

When y_B* = ±√(1-γ), we use (<ref>) to solve (<ref>) for y_A*, which is equivalent to solving the cubic equation

    β_3 y_A*^3+β_2 y_A*^2 +β_1 y_A*+β_0 
    = 0,

where

    β_0 
       = 
    bγ(I-a)+y_B*+b(-y_B*+y_B*^3/3+a+γ y_B*),
    
    β_1
       =
    - bγ,
    
    β_2
       =
    -y_B*+b(y_B*-y_B*^3/3-a-γ y_B*),
    
    β_3
       =
    2b γ/3.


The solutions of (<ref>) for y_A* as a function of γ and I are:

    y_A*,k   = -13β_3(β_2+C_k+β_2^2-3β_1 β_3C_k),

where, for  k = 1, 2, 3,

    C_k    = (√(-3)-1/2)^k-1(σ-√(-27β_3^2 Δ)2)^1/3,
    Δ   = 18β_3β_2β_1β_0-4β_0β_2^3+β_2^2β_1^2-4β_3β_1^3-27β_3^2β_0^2,
    σ   =2β_2^2 -9β_3 β_2 β_1+27 β_3^2 β_0.


If Δ>0, there are three real solutions and, if Δ<0, there is one real solution.

The Jacobian of (<ref>) for the folded singularities with y_B* = ±√(1-γ) has the form:

    D_𝐲ρ (y_A*, y_B*) = ([           0   -2y_B*ξ_A; -γ-2y_A*ξ_B    1-y_A*^2;             ]),

where ρ = (ρ_1, ρ_2)^⊤, ξ_A = y_A*-bz_A*, and ξ_B = y_B* -b z_B*.
The trace and determinant of the Jacobian are used to classify each singularity. They are given by

    (D_𝐲ρ (y_A*, y_B*))   = 1 - y_A*^2,
    
    (D_𝐲ρ (y_A*, y_B*))   = -2y_B*ξ_A(γ+2y_A*ξ_B).

When (D_𝐲ρ (y_A*, y_B*)) >0, the eigenvalues have the same sign, so the singularity is a folded node. The stability of the node can be determined by looking at the sign of the trace. When (D_𝐲ρ (y_A*, y_B*))<0, the eigenvalues have opposite signs, and the singularity is a folded saddle.

In Figure <ref>, the signs of the eigenvalues associated with the folded singularities are depicted as a function of I and γ. The white regions in this figure correspond to values of I and γ where the desingularized system (<ref>) has a single equilibrium point (i.e. Δ < 0).







§ DYNAMICS BY REGION



In this section, we apply the analytical results for the desingularized system (<ref>) to draw conclusions about the original system (<ref>). In so doing we provide details of the computations used to produce Figure <ref> and the characterization of each of the seven regions, as described in Section <ref>. We prove the stability of limit cycles of the Hopf bifurcations in model neuron B. 
 We prove necessary conditions for MMOs and sufficient conditions for phase locking in terms of I and γ. 

In the following proposition, we compute the value of the unique equilibrium point (<ref>), when Assumption <ref> holds, and its stability as a function of I and γ.

 
  Consider the directed two-FN system (<ref>) and let Assumption <ref> hold. For any fixed I and γ, there exists a unique equilibrium point denoted by 𝐩_*=(y_A*(I), z_A*(I), y_B*(I,γ), z_B*(I,γ)). Then,






  * 𝐩_* is nonhyperbolic if I and γ satisfy 

    σ_1(I,γ)    = 1-bϵ-y_A*^2 = 0,
    
    or    σ_2(I,γ)    = 1-bϵ-γ-y_B*^2 = 0.



  * 𝐩_* is hyperbolic if I and γ do not satisfy (<ref>). If σ_1(I,γ)<0 and σ_2(I,γ)<0, 𝐩_* is attracting. If σ_1(I,γ)>0 and σ_2(I,γ)>0, 𝐩_* is repelling. If σ_1(I,γ)σ_2(I,γ)<0,  𝐩_* is a saddle.







 Solving for the equilibrium point of (<ref>), we first compute
 y_A* as a function of I as
 
    y_A*   = 
     (3(I-a)/2 + √((3(I-a))^2/4+b̃^3))^1/3 
    +(3(I-a)/2 - √((3(I-a))^2/4+b̃^3))^1/3,

where b̃ = 1/b-1. Similarly, by leveraging (<ref>), we can write y_B* as the following function of I and γ,

    y_B*   = (3(γ y_A*-a)/2 + √(9(γ y_A*-a)^2/4+(b̃+γ)^3))^1/3  
    +(3(γ y_A*-a)/2 - √(9(γ y_A*-a)^2/4+(b̃+γ)^3))^1/3.


Then z_A*= 1/by_A* and z_B* = 1/by_B*.

We compute the linearization of (<ref>) around 𝐩_*. We let ν_A = (y_A, z_A), ν_B = (y_B, z_B), and ξ = (ζ_A, ϵξ_A,ζ_B,ϵξ_B). The Jacobian of (<ref>) evaluated at p_* is

D_(ν_A,ν_B)ξ(𝐩_*) = ([ 1 - y_A*^2         -1          0          0;          ϵ       -b ϵ          0          0;          γ          0 1-y_B*^2-γ         -1;          0          0          ϵ       -b ϵ ]).



The linearization is block triangular, so the eigenvalues of the Jacobian are union of the eigenvalues of diagonal blocks. This means that local stability can be determined through linearization of each FN model separately. The eigenvalues for the first and second blocks are

    λ_1,2   = 1/2(1 - y_A*^2 -b ϵ)
    ±1/2√((1 -y_A*^2 -b ϵ)^2-4 ϵ (1-b+y_A*^2b)),
    
    λ_3,4   = 1/2(1 - y_B*^2 - γ -b ϵ)
    ±1/2√((1 -y_B*^2 -γ -b ϵ)^2-4 ϵ (1-b+y_B*^2b +γ)).


The sign of the real part of the eigenvalues will be determined by the sign of the first term. The first term of λ_1,2 is zero when σ_1(I,γ) = λ_1 + λ_2 = 1 - bϵ - y_A*^2=0. The first term of λ_3,4 is zero when σ_2(I,γ) = λ_3 + λ_4 = 1-bϵ-γ - y_B*^2=0.  Thus, 𝐩_* is nonhyperbolic when σ_1 = 0 or σ_2=0.

The stability of 𝐩_* when σ_1 σ_2 ≠ 0 is derived from the signs of the eigenvalues of D_(ν_A,ν_B)ξ.






The one-dimensional manifolds of nonhyperbolic equilibrium points in (I,γ) space, {(I,γ):σ_1(I,γ)=0  or σ_2(I,γ)=0 }, correspond to the points where A and B undergo Hopf bifurcations.










 §.§ Quiescence: Region (1)



Given the two-FN system (<ref>) and Assumption <ref>, if I<I_0A then A converges to a single stable equilibrium point, (y_A*,z_A*), which is quiescent. The value y_A*, independent of γ, is too low to induce firing in B, i.e., B is quiescent. The phase plane and nullclines for an (I,γ) pair in this regime are shown in Figure <ref>. 

To fully understand the behavior of B, we examine the desingularized system (<ref>), which has seven singularities for the parameter values in region (1). These include one stable ordinary singularity, one stable folded singularity, two unstable folded singularities, and three saddle folded singularities (only five of the seven singularities are visible in Figure <ref>). Due to the presence of the unstable folded singularities in (<ref>), robust families of canards that compose  small oscillations of B could arise in the original system (<ref>) in region (1).










 §.§ Hopf bifurcations in B: Regions (2), (3), and (4)




For the two-FN system (<ref>), if I>I_1A then A is saturated. We prove conditions for when B will be quiescent, firing, or saturated and provide illustrative examples of the desingularized system nullclines and phase plane for representative parameter values.


Consider the two-FN system (<ref>) and Assumption <ref>.	Let I>I_1A and γ<1-bϵ. There exist two curves of Hopf bifurcations defined by

    I_1B,0B(γ)    =  b̃y_A*±+y_A*±^33+a,

where y_A*± is

    y_A*± = ±1/γ(1/3 (1-ϵ b -γ)^3/2+(b̃+γ)√(1-ϵ b -γ)+a).

B transitions from quiescent to firing through a Hopf bifurcation at I=I_0B(γ) and from firing to saturated through another Hopf bifurcation at I=I_1B(γ).

Moreover, there exists γ_* such that, for I<I_1B(γ_*), the following holds. If γ<γ*, the Hopf bifurcation at I_0B(γ) is supercritical and, if γ>γ*, the Hopf bifurcation at I_1B(γ) is subcritical.





The Hopf bifurcations in B occur at nonhyperbolic equilibrium points, which are y_B*± = ±√(1-γ-bϵ) by Proposition <ref>. 

Substituting y_B*±=±√(1-γ-bϵ) and (<ref>) into the equilibrium solution for (<ref>) gives the critical values

    y_A*± = ±1/γ(
    	1/3 (1-ϵ b -γ)^3/2
    	+(b̃+γ)√(1-ϵ b -γ)+a),

Substituting y_A*± and (<ref>) into (<ref>) gives the values 

    I_1B,0B = b̃ y_A*± +y_A*±^33+a.


For a fixed γ, we check the conditions of Proposition <ref> for the bifurcation parameter I. First, we transform (y_B*,z_B*) to the origin (0,0), by introducing y_0 = y_B - y_B* and z_0 = z_B - z_B*. With this change of variable, the dynamics of B (<ref>-<ref>) can be expressed as

    dy_0dt   = 
    (1-γ-y_B*^2)y_0-y_0^3/3-y_0^2 y_B* -z_0,
    
    dz_0dt   = ϵ (y_0-b z_0),

and the Jacobian of (<ref>) evaluated at the origin is 

J_B(0,0) = ([ 1-y_B*^2-γ         -1;          ϵ       -b ϵ ]).

Now we apply Proposition <ref>.


Condition 1 of Proposition <ref>:
This condition holds because (J_B(0,0)) = 0 at the bifurcation values I = I_0B and I=I_1B.

Condition 2 of Proposition <ref>:
The second condition holds as well because 

    ∂∂ I(λ_3,4(I))|_I = I_0B,1B≠ 0,

where λ_3,4 are given in the proof of Proposition <ref>. 

Condition 2 of Proposition <ref>:
The cubic coefficient α of the Taylor expansion of (<ref>)–(<ref>) (Definition <ref>), which determines whether the Hopf bifurcation is subcritical or supercritical <cit.>, is 

    α = 1/8(2b-2bγ-b^2ϵ-11-b^2ϵ).

At γ=γ_*, α = 0 and B undergoes a “generalized Hopf", or Bautin, bifurcation<cit.>. For γ>γ_*, α<0 and the limit cycles resulting from the Hopf bifurcations are stable (supercritical). Otherwise, the limit cycles are unstable and the bifurcations are subcritical, as for A.










For the two-FN system (<ref>), given Assumption <ref>, if A is saturated, then B transitions from quiescent to firing to saturated as a function of I and γ.




The phase plane and nullclines for (<ref>) in the regions where A is saturated and B is quiescent, firing, and saturated are plotted in Figure <ref>, Figure <ref>, and Figure <ref>, respectively. All three use a fixed value of I given by I=1.7. As γ is varied, changes to the shape, position, and existence of the nullclines and qualitative differences in the phase portraits can be seen.

In region (2), there are seven singularities, of which only five are visible in Figure <ref>. The ordinary singularity is stable, and there are three stable folded nodes and three folded saddles. This is a region where A is saturated and B is quiescent, so the two-FN system has a unique stable equilibrium point.








In region (3) there are five singularities, of which only four are visible in Figure <ref>. For low γ, there are small subregions near the boundary where there are seven singularities. The transition from seven to five singularities occurs through a saddle-node bifurcation between a folded saddle point and a stable folded node. The ordinary singularity is a saddle.

In region (4)  there can be one, three, five, or seven equilibrium points. Parameter choices above γ = 1, which corresponds to γ = 1-bϵ in the two-FN system, result in a single singularity as shown in Figure <ref>. Since γ>1, the folded singularities corresponding to y_B* = ±√(1-γ) no longer exist, and the ordinary singularity is a stable node.







 §.§ Phase-locking: Region (5)



Now, consider the two-FN system (<ref>) with I_0A<I<I_1A and γ>1-bϵ. For all I, γ in this regime, the linearization of (<ref>) around the equilibrium point p_* has two eigenvalues λ_1,2 with positive real part, and two eigenvalues λ_3,4 with negative real part. A is in the firing regime since (λ_1,2)>0. Furthermore, for all I, γ under consideration, (λ_3,4)<0, so B will follow the limit cycle from A and the oscillators will be phase locked. 






For the two-FN system (<ref>), given Assumption <ref>, if A is firing and γ>1-bϵ, B is firing and A and B are phased locked.



The phase plane and nullclines for the two-dimensional, desingularized system (<ref>) with I_0A<I<I_1A and γ>1 (region (5)) are plotted in Figure <ref>. This is the same region discussed above, where γ>1-bϵ becomes γ>1 in the singular limit. In this regime, the singularities at y_B* = ±√(1-γ) no longer exist because γ>1. There is a single ordinary singularity, which is a saddle. Note that there are still folded singularities corresponding to y_A* = ± 1, but that they are the points, I_0A and I_1A, where the stability changes in the ordinary singularity.










 §.§ Mixed mode oscillations: Regions (6) and (7)




Again consider the two-FN system (<ref>) when I_0A<I<I_1A so that A is firing. When γ<1-bϵ, a range of dynamics is observed in simulation. For γ<<1, the influence of A is small, so B exhibits only small oscillations that stay close to (y_B*, z_B*). As γ is increased, the oscillations of A can be large enough to yield a mixed mode oscillation (MMO). As γ is increased still further and the influence of A becomes increasingly strong, B approaches 
firing at the same frequency as A.

To better understand these transitions, we study the bifurcation diagrams and phase planes of the two-dimensional desingularized system (<ref>).  We prove necessary conditions for the existence of MMOs as a function of I and γ.

In region (6) there are seven singularities, of which only six are visible in Figure <ref>. In this region, the ordinary singularity of the desingularized system (<ref>) is a saddle, whereas the folded singularities include three saddles, one unstable node, and two stable nodes. Co-existence of the unstable folded singularity and the ordinary saddle singularity allows existence of canards and MMOs in the two-FN system.





In region (7) there are seven singularities, of which only six are visible in Figure <ref>. In this region, the ordinary singularity of the desingularized system (<ref>) is an unstable node, whereas the folded singularities include four saddles, and two stable nodes. As a consequence, canards and MMOs do not exist in the two-FN system, and both A and B are firing and phase-locked.




We next compute the boundary between regions (6) and (7), shown by I_* in Figure <ref>. The boundary is defined by points at which there is a transcritical bifurcation between the ordinary singularity and a folded singularity, called FSN type II bifurcation, where the ordinary singularity transitions from a saddle to unstable node and the folded singularity transitions from an unstable node to a saddle. This transcritical bifurcation is a known location for the onset of MMOs,  <cit.> 
so computing I_* gives necessary conditions for the existence of MMOs. 



	
A system 𝐱̇ = 𝐟(𝐱,μ), admits a transcritical bifurcation at (𝐱_0,μ_0) if






  * D_𝐱𝐟(𝐱_0,μ_0) has a simple 0 eigenvalue with right eigenvector v and left eigenvector w.


  * w D^2_𝐱μ𝐟(𝐱_0,μ_0)v≠ 0.


  * w D^2_𝐱𝐟(𝐱_0,μ_0)(v^⊤,v^⊤)^⊤≠ 0.







Consider the desingularized system (<ref>) with fixed γ < 1 and let 

    I_*(γ)   = 1/3b^3γ^3(√(1-γ) + 2b√(1-γ)/3-ba)^3+a.
      
Then system (<ref>) admits a transcritical bifurcation at (p, I_*(γ)), where p = (y_A* , y_B*)^⊤ is an ordinary singularity of (<ref>), i.e., p solves (<ref>) and (<ref>).





To show the transcritical bifurcation, we apply Proposition <ref> to (<ref>). The Jacobian, D_𝐲ρ,  of (<ref>) is 

([  (1-y_B^2-γ)(1-b+by_A^2)          -2y_B(y_A-bz_A);       -γ- 2y_A(y_B-bz_A) (1-y_A^2)(1-b+by_A^2+bγ) ]).


Condition 1 of Proposition <ref>: 
Evaluating the Jacobian of (<ref>) at (y_A,y_B)^⊤=p and I=I_* gives

D_𝐲 ρ (p,I_*) 
= ([        0        0;       -γ 1-y_A*^2 ]),

which has a zero eigenvalue with a left eigenvector w = (1,0) and a right eigenvector v = (1, γ/1-y_A*^2)^⊤.

Condition 2 of Proposition <ref>:
Taking the derivative of D_𝐲ρ with respect to I and evaluating at (y_A,y_B)^⊤=p and I=I_* gives 

D^2_𝐲I ρ(p,I_*)
= ([      0 2by_B*;      0      0 ]).

Then, multiplying D^2_𝐲γρ(p,I_*) from left by w and from right by v, we have

w (D^2_𝐲γρ(p,I_*) ) v 
= 
2bγy_B*/1-y_A*^2,

which is always nonzero. 

Condition 3 of Proposition <ref>:
Evaluating D^2_𝐲ρ at p=(y_A,y_B)^⊤ and I=I_* gives

D^2_𝐲 ρ(p,I_*) 
= 
([          0 2y_B*(s-1)          0     -2y_A*; 2y_B*(s-1)          0     -2y_A*     2y_B*s ]),

where s=b(1-y_A*^2). Then, multiplying D^2_𝐲ρ from left by w and from right by (v^⊤,v^⊤)^⊤, we have	

    w( D^2_𝐲ρ(p,I_*)) 
    ([ v; v ]) 
    = 
    -2γ1-y_A*^2(y_B*(1-b(1-y_A*^2))+y_A*),

which is also nonzero as shown in Figure <ref>. 






 

A necessary condition for canard-induced MMOs is the existence of a folded node with a return mechanism, since the small canard solutions that form the small oscillations are only found in this context.<cit.> The folded node has a corresponding family of canard solutions because there are many trajectories that cross from the stable to the unstable branch of the critical manifold through the folded node singularity. Furthermore, the return mechanism is required for MMOs because, after each relaxation oscillation or canard trajectory, the dynamics must return near the singularity in order for the MMO to persist. In this case, MMOs are only possible for I<I_*(γ), since that is where the folded singularity is an unstable folded node and the ordinary singularity is a saddle node.

    

Consider the two-FN system (<ref>). For I_0A<I<I_*(γ), this system exhibits MMOs and if I_*(γ)<I<I_1A, it exhibits phase locking.




     
To highlight the location of the transcritical bifurcation (FSN II) and compare to features in Figure <ref>, we show it as the blue star in the bifurcation diagram of (<ref>) in Figure <ref>a.
Here, by treating I as the bifurcation parameter, and maintaining γ fixed at value 0.4, we can observe that the ordinary singularity transitions from stable (thick green) to saddle (thick gray) at I=I_0A (red star). Simultaneously, an unstable folded node (thin dashed red) becomes a folded saddle node (thin gray) and a folded saddle node becomes a stable folded node (thin green). For I>I_0A, MMOs are possible due to the combination of the unstable folded node and unstable ordinary singularity.

Then, at I=I_*(γ) (blue star), derived in Proposition <ref>, the ordinary singularity and a folded singularity swap stability properties in a transcritical bifurcation which can be classified as an FSN II bifurcation. Also, this FSN II bifurcation in the desingularized system (<ref>) corresponds to the generalized Hopf bifurcation in the directed two-FN system (<ref>). For I>I_*, (<ref>) exhibits phase locking and MMOs are no longer possible.

At I=I_1A (magenta star), the ordinary singularity returns to a saddle and two folded saddle nodes become two stable folded nodes. For I>I_1A, phase locking in (<ref>) is no longer possible. We also note that when I is just above I_1A, one folded saddle node merges with a stable folded node through a saddle node bifurcation.


A special case of the transcritical (FSN II) bifurcations occurs when γ and I satisfy I_*(γ)=I_1A. In this case, there is a codimension two bifurcation where the eigenvalues of the linearization of (<ref>) about the ordinary singularity and the eigenvalues of the linearization of (<ref>) about the folded singularity, y_B* = -√(1-γ), are equal to zero. The codimension two bifurcation is illustrated by the orange star in the bifurcation diagram of (<ref>) in Figure <ref>b, for γ = 0.22. If we fix I_0A < I < I_1A and decrease the value of γ below the codimension two value, then MMOs are always possible.  







§ DIRECTED TREE OF FN MODEL NEURONS




In this section, we consider an extension of the previous results to a directed chain of coupled FN models. We leverage the connection between the desingularized system and the directed two-FN system to find sufficient conditions for phase locking. 

Consider a system of k FN model neurons with dynamics

𝐱̇ = 𝐟(𝐱, 𝐈,γ),
 
where 𝐱∈R^2k, 𝐈∈R^k-1, and γ∈R^k-1. All FN models receive an external input except for the last in the chain. Then, by allowing heterogeneity in the external inputs and coupling strengths, the linearization around the equilibrium point can be expressed as

D_𝐱𝐟
= 
([   J_1 0_2×2 0_2×2 0_2×2     ⋯ 0_2×2;   Γ_1   J_2 0_2×2 0_2×2     ⋯ 0_2×2; 0_2×2   Γ_2   J_3 0_2×2     ⋱     ⋮;     ⋮     ⋱     ⋱     ⋱     ⋱ 0_2×2; 0_2×2 0_2×2     ⋱ Γ_k-2 J_k-1 0_2×2; 0_2×2 0_2×2     ⋯ 0_2×2 Γ_k-1   J_k ]),

where the first diagonal block is given by

J_1 
= 
([ 1-y_1^2      -1;       ϵ     -bϵ ]),

and the subsequent diagonal blocks are given by

J_i 
= 
([ 1-y_i^2-γ_i-1            -1;             ϵ           -bϵ ]),
  i ∈{2,…,k}.

The blocks on the lower diagonal are

Γ_i 
= 
γ_i 
([ 0 1; 0 0 ]),
   i ∈{1,…,k-1}.


Due to the lower block triangular structure of the linearization, local stability of the equilibrium can be determined by studying the eigenvalues of the diagonal blocks. Similar to the analysis at the beginning of Section <ref>, we begin by solving for the equilibrium point. The equilibrium of the first model neuron is given by

    y_1*   
    = (3(I_1-a)/2 + √((3(I_1-a))^2/4+b̃^3))^1/3
    + (3(I_1-a)/2 - √((3(I_1-a))^2/4+b̃^3))^1/3.

The equilibrium of the i-th model neuron is given by

    y_i*   = (3Ĩ/2 + √((3Ĩ)^2/4+(b̃+γ_i-1)^3))^1/3
    + (3Ĩ/2 - √((3Ĩ)^2/4+(b̃+γ_i-1)^3))^1/3,

where Ĩ = γ_i-1 y_i-1*+I_i-a, i∈{2,…,k}. The eigenvalues of the individual diagonal blocks are

    λ_1,2 
    = 
    1/2(1 - y_1*^2 -b ϵ)
    ±1/2√((y_1*^2 +b ϵ-1)^2-4 ϵ (1-b+y_1*^2b)),
    
       λ_2i-1,2i 
    = 
    1/2(1 - y_i*^2 - γ_i-1 -b ϵ)
    ±1/2√((y_i*^2 +γ_i-1 +b ϵ-1)^2-4 ϵ (1-b+y_i*^2b +γ_i-1)),

where i =2,…, k. The Hopf bifurcations in the i-th model neuron occur at

    I_H±   
    = ±1/3(1-γ_i-1-bϵ)^3/2±√(1-γ_i-1-bϵ)(b̃+γ_i-1)-Ĩ.


As a directed tree can be decomposed into a collection of directed chains, these results generalize to directed trees as well. In Figure <ref>, we illustrate with the directed chain that starts with the light orange FN model and is directed to the right to the cyan FN model. The first FN model (light orange) receives an input I=1.2, which ensures that it is in the firing regime. The coupling strength to the second FN model (dark orange) with input I=0.4 ensure that the second FN model is in region (6) where MMOs are possible. However, in this case no MMOs are exhibited. The coupling strength to the third FN model (dark cyan) with zero input ensure that it too is in region (6). In this case, MMOs induced by canards are exhibited. The active signal has frequency half that of the first and second FN models. As a result, the input to the fourth FN model (cyan) is an MMO; this case was not covered in our two-FN system analysis. The fourth FN model responds to incoming canards with almost no activity and incoming spikes with a small canard. The frequency of the small canards in the fourth FN model is the same as the frequency of the active signal of the third FN model.














§ DISCUSSION



In this work, we study a system of two FN model neurons in a setting where the first FN model has a constant external input I, the second FN model has no input, and there is a unidirectional coupling with strength γ from the first FN model to the second. We study and rigorously characterize all of the different regions of dynamic behavior for the two-FN system in I-γ space. We prove new necessary conditions in terms of both I and γ for the existence of canards and MMOs. We leverage this result to find a similarly new sufficient condition for phase locking and extend to systems of FN models in directed tree networks. We illustrate for a directed chain of four FN models, where canards, MMOs, and frequency halving is observed as predicted.

Further investigation of the two-FN system is needed to determine the threshold between MMOs and canard trajectories without MMOs, which have been observed in simulation. This threshold has been studied numerically, as well as the chaotic behavior at the boundaries between types of MMOs, e.g., in <cit.>. An analytical understanding of these phenomena would add significantly to the literature on canards and MMOs.

Future directions include consideration of more diverse graph structures that include loops within the graph and a more detailed analysis of the MMOs in these systems. General results have been found for finite dimensional fast-slow systems, which could be applied in this context <cit.>. Incorporating heterogeneous model parameters is another area of future investigation. Changing ϵ changes the frequency of oscillation and the timescale of the FN model, so a network of FN models with differing values of ϵ would be a compelling system for exploring canard phenomena in three or more distinct timescales.












§ ACKNOWLEDGMENTS



This work was jointly supported by the National Science Foundation under NSF-CRCNS grant DMS-1430077 and the Office of Naval Research under ONR grant N00014-14-1-0635. This material is also based upon work supported by the National Science Foundation Graduate Research  Fellowship under grant DGE-1656466. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.








10

hodgkin1952quantitative
A. L. Hodgkin and A. F. Huxley.
A quantitative description of membrane current and its application to
  conduction and excitation in nerve.
The Journal of Physiology, 117(4):500–544, 1952.

fitzhugh_mathematical_1955
R. FitzHugh.
Mathematical models of threshold phenomena in the nerve membrane.
The Bulletin of Mathematical Biophysics, 17(4):257–278,
  1955.

fitzhugh1961impulses
R. FitzHugh.
Impulses and physiological states in theoretical models of nerve
  membrane.
Biophysical Journal, 1(6):445–466, 1961.

nagumo1962active
J. Nagumo, S. Arimoto, and S. Yoshizawa.
An active pulse transmission line simulating nerve axon.
Proceedings of the IRE, 50(10):2061–2070, 1962.

izhikevich2007dynamical
E. M. Izhikevich.
Dynamical Systems in Neuroscience: The Geometry of Excitability
  and Bursting.
MIT press, 2007.

russo_global_2010
G. Russo and J. E. Slotine.
Global convergence of quorum-sensing networks.
Physical Review E, 82(4):041919, 2010.

schaub2016graph
M. T. Schaub, N. O'Clery, Y. N. Billeh, J. Delvenne, R. Lambiotte, and
  M. Barahona.
Graph partitions and cluster synchronization in networks of
  oscillators.
Chaos: An Interdisciplinary Journal of Nonlinear Science,
  26(9):094821, 2016.

sorrentino2016complete
F. Sorrentino, L. M. Pecora, A. M. Hagerstrom, T. E. Murphy, and R. Roy.
Complete characterization of the stability of cluster synchronization
  in complex dynamical networks.
Science Advances, 2(4):e1501737, 2016.

davison_sync_2016
E. N. Davison, B. Dey, and N. Ehrich Leonard.
Synchronization bound for networks of nonlinear oscillators.
Proceedings of the 54th Annual Allerton Conference on
  Communication, Control and Computing, 2016.

aminzare2018cluster
Z. Aminzare, B. Dey, E. N. Davison, and N. Ehrich Leonard.
Cluster synchronization of diffusively-coupled nonlinear systems: A
  contraction based approach.
Journal of Nonlinear Science, 2018.

hoff_numerical_2014
A. Hoff, J. V. Santos, C. Manchein, and H. A. Albuquerque.
Numerical bifurcation analysis of two coupled FitzHugh-Nagumo
  oscillators.
The European Physical Journal B, 87(7):151, July 2014.

campbell2001multistability
S. Campbell and M. Waite.
Multistability in coupled FitzHugh-Nagumo oscillators.
Nonlinear Analysis, 47:1093–1104, 2001.

benoit1990canards
É. Benoît.
Canards et enlacements.
Publications Mathématiques de l'Institut des Hautes
  Études Scientifiques, 72(1):63–91, 1990.

szmolyan2001canards
P. Szmolyan and M. Wechselberger.
Canards in ℝ^3.
Journal of Differential Equations, 177(2):419–453, 2001.

kuehn2015multiple
C. Kuehn.
Multiple Time Scale Dynamics, volume 191.
Springer, 2015.

wechselberger_existence_2005
M. Wechselberger.
Existence and bifurcation of canards in ℝ^3 in the case
  of a folded node.
SIAM Journal on Applied Dynamical Systems, 4(1):101–139,
  January 2005.

milik1998geometry
A. Milik, P. Szmolyan, H. Löffelmann, and E. Gröller.
Geometry of mixed-mode oscillations in the 3-D autocatalator.
International Journal of Bifurcation and Chaos, 8(03):505–519,
  1998.

benoit1983systemes
É. Benoît.
Systemes lents-rapides dans ℝ^3 et leurs canards.
Astérisque, 109(110):159–191, 1983.

tchizawa_two_2013
K. Tchizawa.
On the two methods for finding 4-dimensional duck solutions.
Applied Mathematics, 05(01):16, December 2013.

wechselberger2012propos
M. Wechselberger.
A propos de canards (apropos canards).
Transactions of the American Mathematical Society,
  364(6):3289–3309, 2012.

desroches_geometry_2008
M. Desroches, B. Krauskopf, and H. Osinga.
The geometry of slow manifolds near a folded node.
SIAM Journal on Applied Dynamical Systems, 7(4):1131–1162,
  January 2008.

krupa2010local
M. Krupa and M. Wechselberger.
Local analysis near a folded saddle-node singularity.
Journal of Differential Equations, 248(12):2841–2888, 2010.

Curtu_Rubin_Canard-MMO
R. Curtu and J. Rubin.
Interaction of canard and singular hopf mechanisms in a neural model.
SIAM Journal on Applied Dynamical Systems, 10(4):1443–1479,
  2011.

desroches2012mixed
M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H. M. Osinga, and
  M. Wechselberger.
Mixed-mode oscillations with multiple time scales.
SIAM Review, 54(2):211–288, 2012.

petrov1992mixed
V. Petrov, S. K. Scott, and K. Showalter.
Mixed-mode oscillations in chemical systems.
The Journal of Chemical Physics, 97(9):6191–6198, 1992.

moehlis2002canards
J. Moehlis.
Canards in a surface oxidation reaction.
Journal of Nonlinear Science, 12(4), 2002.

rubin2007giant
J. Rubin and M. Wechselberger.
Giant squid-hidden canard: The 3D geometry of the
  Hodgkin–Huxley model.
Biological Cybernetics, 97(1):5–32, 2007.

rotstein2008canard
H. G. Rotstein, M. Wechselberger, and N. Kopell.
Canard induced mixed-mode oscillations in a medial entorhinal cortex
  layer II stellate cell model.
SIAM Journal on Applied Dynamical Systems, 7(4):1582–1611,
  2008.

desrochesmixedmode2008
M. Desroches, B. Krauskopf, and H. M. Osinga.
Mixed-mode oscillations and slow manifolds in the self-coupled
  FitzHugh-Nagumo system.
Chaos: An Interdisciplinary Journal of Nonlinear Science,
  18(1):015107, March 2008.

tchizawa2002winding
K. Tchizawa and S. Campbell.
On winding duck solutions in ℝ^4.
Proceedings of Neural, Parallel, and Scientific Computations,
  2:315–318, 2002.

ginoux2015canards
J. Ginoux and J. Llibre.
Canards existence in FitzHugh-Nagumo and Hodgkin-Huxley neuronal
  models.
Mathematical Problems in Engineering, 2015, 2015.

krupa2014weakly
M. Krupa, B. Ambrosio, and M. A. Aziz-Alaoui.
Weakly coupled two-slow–two-fast systems, folded singularities and
  mixed mode oscillations.
Nonlinearity, 27(7):1555, 2014.

rosenblum1996phase
M. G. Rosenblum, A. S. Pikovsky, and J. Kurths.
Phase synchronization of chaotic oscillators.
Physical review letters, 76(11):1804, 1996.

van1928lxxii
V. van der Pol and J. van der Mark.
LXXII. The heartbeat considered as a relaxation oscillation, and an
  electrical model of the heart.
The London, Edinburgh, and Dublin Philosophical Magazine and
  Journal of Science, 6(38):763–775, 1928.

glass1982fine
L. Glass and R. Perez.
Fine structure of phase locking.
Physical Review Letters, 48(26):1772, 1982.

guevara1981phase
M. R. Guevara, L. Glass, and A. Shrier.
Phase locking, period-doubling bifurcations, and irregular dynamics
  in periodically stimulated cardiac cells.
Science, 214(4527):1350–1353, 1981.

matsumoto1987chaos
G. Matsumoto, K. Aihara, Y. Hanyu, N. Takahashi, S. Yoshizawa, and J. Nagumo.
Chaos and phase locking in normal squid axons.
Physics Letters A, 123(4):162–166, 1987.

kopell2002mechanisms
N. Kopell and G. B. Ermentrout.
Mechanisms of phase-locking and frequency control in pairs of coupled
  neural oscillators.
Handbook of Dynamical Systems, 2:3–54, 2002.

rutishauser2010human
U. Rutishauser, I. B. Ross, A. N. Mamelak, and E. M. Schuman.
Human memory strength is predicted by theta-frequency phase-locking
  of single neurons.
Nature, 464(7290):903, 2010.

wang2010neurophysiological
X. Wang.
Neurophysiological and computational principles of cortical rhythms
  in cognition.
Physiological Reviews, 90(3):1195–1268, 2010.

varela2001brainweb
F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie.
The brainweb: Phase synchronization and large-scale integration.
Nature Reviews Neuroscience, 2(4):229, 2001.

medvedev2001synchronization
G. S. Medvedev and N. Kopell.
Synchronization and transient dynamics in the chains of electrically
  coupled FitzHugh–Nagumo oscillators.
SIAM Journal on Applied Mathematics, 61(5):1762–1801, 2001.

rocsoreanu2012fitzhugh
C. Rocsoreanu, A. Georgescu, and N. Giurgiteanu.
The FitzHugh-Nagumo Model: Bifurcation and Dynamics,
  volume 10.
Springer Science & Business Media, 2012.

guckenheimer2013nonlinear
J. Guckenheimer and P. Holmes.
Nonlinear Oscillations, Dynamical Systems, and Bifurcations of
  Vector Fields, volume 42.
Springer Science & Business Media, 2013.

dhooge2003matcont
A. Dhooge, W. Govaerts, and Y. A. Kuznetsov.
MATCONT: A MATLAB package for numerical bifurcation analysis of
  ODEs.
ACM Transactions on Mathematical Software (TOMS),
  29(2):141–164, 2003.

Rinzel1981
J. Rinzel.
Models in Neurobiology, pages 345–367.
Springer New York, Boston, MA, 1981.

holmes_notes
P. Eckhoff and P. Holmes.
A short course in Mathematical Neuroscience, 2015.

hassard1978bifurcation
B. Hassard and Y. H. Wan.
Bifurcation formulae derived from center manifold theory.
Journal of Mathematical Analysis and Applications,
  63(1):297–312, 1978.

baer1986singular
S. M. Baer and T. Erneux.
Singular Hopf bifurcation to relaxation oscillations.
SIAM Journal on Applied Mathematics, 46(5):721–739, 1986.

rinzel1983hopf
J. Rinzel and J. P. Keaner.
Hopf bifurcation to repetitive activity in nerve.
SIAM Journal on Applied Mathematics, 43(4):907–922, 1983.

freire2011stern
J. G. Freire and J. A.C. Gallas.
Stern–Brocot trees in cascades of mixed-mode oscillations and
  canards in the extended Bonhoeffer–van der Pol and the
  FitzHugh–Nagumo models of excitable systems.
Physics Letters A, 375(7):1097–1103, 2011.

wechselberger_canard_2013
M. Wechselberger, J. Mitry, and J. Rinzel.
Canard theory and excitability.
In P. E. Kloeden and C. Pötzsche, editors, Nonautonomous
  Dynamical Systems in the Life Sciences, number 2102 in Lecture
  Notes in Mathematics. Springer International Publishing, 2013.

takens1973unfoldings
F. Takens.
Unfoldings of certain singularities of vectorfields: Generalized
  Hopf bifurcations.
Journal of Differential Equations, 14(3):476–493, 1973.

golubitsky1981classification
M. Golubitsky and W. F. Langford.
Classification and unfoldings of degenerate Hopf bifurcations.
Journal of Differential Equations, 41(3):375–415, 1981.

brons2006mixed
M. Brøns, M. Krupa, and M. Wechselberger.
Mixed mode oscillations due to the generalized canard phenomenon.
Fields Institute Communications, 49:39–63, 2006.





