Cascades of Periodic Solutions in a Neural Circuit With Delays and Slow-Fast Dynamics

We analyse periodic solutions in a system of four delayed differential equations forced by periodic inputs representing two competing neural populations connected with fast mutual excitation and slow delayed inhibition. The combination of mechanisms generates a rich dynamical structure that we are able to characterize using slow-fast dissection and a binary classification of states. We previously proved the existence conditions of all possible states 1:1 locked to the inputs and applied this analysis to the tracking of the rhythms perceived when listening to alternating sequences of low and high tones. Here we extend this analysis using analytical and computational tools by proving the existence a set of n:1 periodically locked states and their location in parameter space. Firstly we examine cycle skipping states and find that they accumulate in an infinite cascade of period-incrementing bifurcations with increasing periods for decreasing values of the local input strength. Secondly we analyse periodic solutions that alternate between 1:1 locked states that repeat after an integer multiple of the input period (swapping states). We show that such states accumulate in similar bifurcation cascades with decreasing values of the lateral input strength. We report a parameter-dependent scaling constant for the ratio of widths of successive regions in the cascades, which generalises across cycle skipping and swapping states. The periodic states reported here - emergent behaviours in the model - can be linked to known phenomena in auditory perception that are beyond the original scope of the model’s design.


INTRODUCTION
Differentiating between sound sources that overlap or are interleaved in time is a fundamental part of auditory perception (the so-called cocktail party problem) [1,2]. The neural computations underpinning this segregation of sound sources has been the subject of dynamical systems models in various frameworks that focus on feature separation (e.g., differences in pitch) between interleaved, periodic sequences of sounds [3][4][5]. However, these studies failed to address neural computations where temporal mechanisms interact with inputs on the timescale of the input period. A recent study by the authors [6] addressed these computations in a framework inspired by the structure of the auditory processing pathway and by common features of neural computation in cortical brain areas (excitatory-inhibitory interactions, time-scale separation, transmission delays). Whilst our previous study focused on 1:1 locked states, here we report on more exotic dynamics that organise in cascades of bifurcations leading to solutions with increased period.
The encoding of sensory information in cortex involves subpopulations of tens of thousands to millions of neurons that are suitably represented by coarse-grained variables representing e.g., the mean firing rate of the population. The Wilson-Cowan equations [7] for localised populations of neurons describes the firing rates of neural populations, and they are widely used in small networks with excitatory and inhibitory coupling, intrinsic synaptic dynamics, neural adaptation and a nonlinear gain function [8][9][10].
Timescale separation is a common feature of models at the single neuron level [11,12], and in populations of neurons [13]. Slow-fast analysis including singular perturbation theory has been instrumental in revealing the dynamical mechanisms behind spiking and bursting [12,14] and in explaining complex dynamics in population models of neural competition [10,15]. Extensions of these methods have been applied to systems with delays [16], in non-smooth settings [17] and in networks capable of instrinsically generating patterns of rhythmic behaviours (so-called central pattern generators, CPGs), such as locomotion, breathing, sleep [18].
Delayed inhibition modelled with slow variables in systems of ODEs representing CPG circuits with spiking units is important for the generation of patterns of activity reminiscent of those considered in this paper [19,20]. Delays in small neural circuits modelled using DDE equations can lead to many interesting phenomena including inhibition-induced Hopf oscillations, oscillator death, multistability and switching between oscillatory solutions [21][22][23]. Two other key features of our study are that 1) units are not instrinsically oscillating and 2) periodic forcing of the units drives oscillations. Periodically forced systems with timescale separation have been explored in models of perceptual competition [24,25], but not in the presence of delays. Periodic solutions in autonomous delay differential equations with Heaviside and monotonic gain functions have been studied analytically in a recent paper from the authors [26].
Periodic orbits can undergo cascades of bifurcations under parameter variation. This phenomena leads to an increase in period and/or complexity of the emerging orbits. The mostly widely studied of these is the period-doubling cascade, a typical route to chaos that can be found in simple dynamical systems such as the smooth, discrete logistic map [27], or the periodicallyforced, continuous-time duffing oscillator [28]. Continuous-time dynamical systems with a separation of timescales can produce mixed-mode oscillations with a sequence of bifurcations leading to additional low-amplitude excursions via a canard-induced mechanisms [29]. Analysis of so-called spike adding bifurcations have been instrumental in understanding bursting phenomena in excitable cells [30]: additional large-amplitude excursions can appear in a smooth bifurcation sequence. Here we show analytically that a cascade of periodic solutions with discrete steps in period that exist in non-overlapping, contiguous parameter regions. As the cascade progresses (here decreasing a parameter) these regions become narrower and the period increases in discrete jumps. These features are akin to a period-incrementing bifurcation sequence as reported in nonsmooth maps [31]. We note that 1) overall these solutions exist in large parameter regions that overlap with 1:1 locked states previously reported, 2) in a non-smooth setting we cannot specify the bifurcation occurring at boundaries between regions and 3) nevertheless the solution branches persist when we relax the non-smooth assumptions as confirmed numerically.
We consider a periodically-driven competition network of two localised Wilson-Cowan units with lumped excitation and inhibition, generalised to include dynamics via inhibitory synaptic variables. The units A and B are driven by stereotyped input signals representative of neural responses in early auditory areas. Whilst mutual excitation between the units is fast, inhibition activates slowly and transmission between the units is delayed. The units' activation function is taken as a Heaviside function restricting the possible state of each unit to be active (ON), inactive (OFF) or in a fast transition between these two states. This combination of modelling assumptions (fixed delays, slow-fast timescale separation, heaviside activation function) allows for all possible model states to be conveniently represented in a binary matricial form with entries specifying the state of the system in time intervals relating to the input timecourse, delay and activation timescales. This approach allows for exact parameterdependent existence conditions to be derived, as for all 1:1 locked states in [6] and for more exotic states as reported here.
The paper is organised as follows. In Section 2 we propose a discontinuous system of delay differential equations to encode auditory perception. In Section 3 we report some results from our previous study of this system [6] that are used here to analyse new network states of interest. Specifically, we classify the possible dynamics of states in the intervals when tones are active and propose a binary matrix representation for each state in the system. In Section 4 we define cycle skipping states with period a multiple of the input period. We also derive analytically their existence conditions and show that they accumulate in cascades of period-incrementing bifurcations for decreasing values of the local input strength. A similar analysis is carried out in Section 5 to study switching states: a class of states in which the units' dynamics alternate between TRperiodic states (thus forming solutions with larger period). Finally, in Section 6 we construct a smooth version of the model consisting of ordinary differential equations and use numerical simulations and continuation to show that our findings are not contingent on the presence of non-smooth components.

MODEL DESCRIPTION
We consider a periodically forced neural network of two Wilson-Cowan populations coupled by fast direct excitation and slow delayed inhibition ( Figure 1A). Inhibition from both units act on a slower time scale and its dynamics is described by synaptic variables. The forcing terms mimic stereotyped synaptic inputs in the auditory cortex [32] to alternating sequences of A and B tones and has been used to study the auditory streaming paradigm (see [6] for more details). The combination between the timescale separation, forcing and delays give rise to a rich repertoire of dynamical states.
The model is described by the following system of four DDEs representing two periodically forced neural populations (A and B populations) coupled by fast direct excitation and slow delayed inhibition The variables u A and u B represent the activity of the A and B populations and have timescale τ, while s A and s B represent the activity of the inhibitory synaptic variables, have timescale τ i and are delayed by fixed D. We will assume τ i to be large and τ to be small. This condition poses system (1) in a slow-fast regime which will enable us to apply analytical tools to study model states. We consider a Heaviside gain function with threshold θ ∈ (0, 1): H(x) 1 if x ≥θ and 0 otherwise, which is a common choice in Wilson Cowan models [10,33]. Parameters a ≥ 0 and b ≥ 0 respectively represent the strengths of the excitatory and inhibitory coupling.
The forcing, square-waved periodic inputs (shown in Figure 1B) i A and i B , are defined by: where c ≥ 0 and d ≥ 0 represent the local and lateral input strengths, respectively; χ I is the indicator function χ I (t) 1 for t ∈ I and 0 otherwise. Throughout this work we will assume that local inputs are stronger than lateral ones, i.e., c ≥ d. We define the parameter TD representing the duration of each tone's presentation and TR representing the time between tone onsets. The intervals when A and B tones are on (active tone intervals) are given by , respectively, and have boundaries given by: for j 1, . . ., ∞. Let us name the set of active tone intervals R as In this work we assume that TD < D and TD + D < TR, which are guaranteed for moderate delays (motivated in [6]). We further assume that a − b < θ and c ≥ θ to avoid trivial saturating or resting dynamics in the system (as discussed in [6]). Throughout this paper system 1 was numerically integrated using Matlab's built-in delay differential equation solver dde23 with fixed time step dt 10 -5 .
Remark 2. and the application of the time shift TR. Since i A (t + TR) i B (t) and i B (t + TR) i A (t), the model is symmetrical under the composition of these two maps (Z 2 symmetry). Thus, periodic solutions of the model are of two types: symmetrical cycles (κ-invariant) and asymmetrical cycles. Any asymmetrical cycle v(t) coexist with its κ-conjugate cycle κ(v (t + TR)).

BACKGROUND
In our previous paper [6] we split system (1) into slow and fast subsystems to characterise their quasi-equlibria. This analysis enabled us to classify all the possible 2TR-periodic states and determine their existence conditions in the parameter space. We now provide a summary of the main ideas, which are important for the new analysis presented here. The following theorem summarises three key properties characterizing the dynamics of each network state during any 2TR interval ( Figure 1C), reported from Theorem 1, Lemma 1 and Lemma 2 in [6].  Figure 1C. Any state of system (1) satisfy the following properties: As shown below, this theorem enabled us to split the dynamics of all network states in four different classes: the disjoint classes of MAIN and CONNECT states and the disjoint classes of SHORT and LONG states ( Figure 1C). Point 1. in Theorem 1 guarantees that the turning ON times for the units can only occur during an active tone interval. Since the A and B active tone intervals (green and purple intervals in Figure 1C) are included in intervals where the delayed synaptic variables are decaying the total inputs to the units in these active tone intervals is monotonically decaying, too. Thus if a unit turns ON at any time in an active tone interval it must remain ON at least until the offset of the interval. Given an active tone interval R ∈Φ we can therefore classify the set of states where the units may turn ON at onset of R (MAIN states) and the set of states where the units turn ON with some non-infinitesimal delay after the onset of R (CONNECT states). The following definition extends this classification for an arbitrary number of active tone intervals (see Figure 1D for examples of MAIN/ CONNECT states in R). • MAIN if ∀R ∈Φ, if ∃t* ∈ R turning ON time for A or B, then t* min(R) • CONNECT if ∃R ∈Φ and ∃t* ∈ R, t* > min(R) turning ON time for A or B, with t * K 0 for τ → 0 Moreover, due to point 2. in Theorem 1 if both units are ON at the offset of an active tone interval R they may continue to be ON for some time after such offset time. However, due to point 3. in the same theorem, there exists a subsequent interval where both units turn OFF (i.e., before the onset of the next active tone interval). This leads to two possible classes of states: SHORT states for which both units are/turn OFF between the offset of R and the onset of its successive active tone interval, and LONG states for which the units remain ON for some time after the offset of R. We now extend this classification for an arbitrary number of active tone intervals (see Figure 1D for examples of SHORT/LONG states in R).
Definition 3.2 (LONG and SHORT states). We define any state of system (1): • SHORT if both units are OFF ∀t ∈ R − I Moreover, in our earlier paper [6] we showed that the dynamics a unit's jump up points during an active tone interval determine the subsequent dynamics in the same interval, as stated in the following theorem, which corresponds to Lemma 5 in [6].
Theorem 2. Given an active tone interval R [α, β] ∈ Φ we have: We used this theorem to define the matrix form (a binary matrix) of each state, which describes the units' dynamics in one active tone interval R. We extended this representation for two active tone intervals I 1 A and I 1 B to describe the dynamics of 2TR-periodic states and use it to determine the existence conditions for these states. Here we extend this tool to periodic solutions with higher periods using an intuitive definition of the matrix form. This tool helps us to visualise the dynamics of the states of interest and to determine their existence conditions.

Intuitive Construction of the Matrix Form
The matrix form for a given state in an active tone interval R [α, β] ∈Φ is the 2 by 3 binary matrix where the first (second) row follow the A (B) unit's dynamics in R. Due to Theorem 2 if a unit turns ON at some time t* ∈ R if must be ON in (t*, β]. Therefore we have only four possible cases to consider for the first row of V. If unit A turns ON at time α then it is ON in all the interval R and we have x A y A z A 1. If unit A is OFF at time α but turns ON after a small delay α + δ due to the excitation from unit B (δ is an infinitesimal delay on the fast time scale, δ ∼ τ) then we have x A 0, y A 1 and z A 1. If unit A turns ON at some intermediate time t* ∈ (α + δ, β) then we have x A 0, y A 0 and z A 1. Otherwise, if unit A is OFF in R all entries in the first row are zero. Analogous arguments lead to the construction of the second row via the dynamics of unit B.
Since the units' dynamics dictate the dynamics of the synaptic variables (Remark 2.1) the matrix form uniquely describes the complete four dimensional dynamics of system (1) in all active tone intervals. However, the matrix form does not provide information about the dynamics between each pair of successive active tone intervals (thus establishing if the state is LONG or SHORT, see Remark 3.1 below).
The matrix form extension to 2TR-periodic states during two active tone intervals I 1 [0, TD] and I 2 [TR, TR + TD] is the 2 × 6 binary matrices of the form where V 1 and V 2 are the matrix forms in the intervals I 1 and I 2 , respectively in Eq. (3). The matrix form of these states is shown in Table 1. One can easily check that the matrix form of each state provides a visualization of the associated dynamics in the active tone intervals I 1 and I 2 shown in Figure 1E. In our previous work [6] we defined the matrix form Eq. (4) rigorously using analytical tools. Here we used an intuitive approach to facilitate the reading. The matrix form definition enabled us to derive the existence conditions for all viable 2TR-periodic states (and to rule out which states are impossible) and to visualise their existence regions in the space of parameters. A similar technique was used to find all CONNECT states and to define their existence conditions. Table 2 shows the matrix form of these states. We omit time histories because these can be visualised from the matrix forms of each state. Figure 1F shows the regions of existence of 2TR-periodic SHORT MAIN and CONNECT states in the (c, η) parameter plane, where η is a newly introduced parameter for scaling the lateral input as d ηc and all other parameters are fixed.
Remark 3.1. The matrix forms shown in Tables 1, 2 do not provide information on the units' dynamics outside the active tone intervals. Therefore, to establish if a state is SHORT or LONG we need to impose additional conditions at the offset times of the active tone intervals. As shown in [6] if both units are ON at the end of an active tone interval, they may continue to be ON after the offset of the active tone interval and turn OFF before the onset of the next active tone interval (LONG condition). The 2T R-periodic states in which at least one unit is OFF at the end of both intervals I 1 and I 2 are those states with matrix form given by S, AP, INT and ZcAP in Table 1 and Table 2 (these states cannot therefore be LONG). The remaining states can be LONG. The analysis of all the combinations of SHORT/LONG and MAIN/CONNECT states is given in [6]. 1 | Matrix form of all 2T R-periodic MAIN states (* asymmetrical states with coexisting conjugate). States were named following their proposed link with percepts in the auditory streaming paradigm ( [6]). Specifically, names starting with S, I and AS represent segregation, integration and bistability, respectively. The ending letter B indicates states in which both units turn ON at the onset of the same active tone interval. The ending letter D indicates states in which a unit turns ON an infinitesimal delay after the other unit for at least one active tone interval. In this section we define and study the conditions for the existence of cycle skipping states with 2TR-multiple period and show that they accumulate in cascades for decreasing values of the local input strength. For cycle skipping states we start from 2TR-periodic MAIN states with active tone intervals I 1 I 1 A and I 2 I 1 B and add a period of 2TR-multiple of silence after I 1 and/or I 2 (silence refers to active tone intervals that do not generate a response). Let us define the sets of 2TR-periodic segregated states S and integrated states I as where SB, SD, IB and ID are the states having matrix form given by their corresponding name ( Table 1). As we shall soon prove, cycle skipping solutions of all other states do not exist.
The existence conditions for cycle skipping states of order k depend on the following quantities: We note that for these quantities L + j ≤ L − j for any j, and that L +/− j ≤ L +/− m for each pair of indexes j, m such that j ≥ m. We now proceed by defining kth-order cycle skipping segregated and integrated states and their matrix form. For simplicity in the notation of the next sections we introduce the vector 0 (0, 0, 0). Definition 4.1. A cycle skipping segregated state s k of order k ∈ N is a periodic solution of system (1) obtained from a segregated state s ∈ S by adding 2kTR periods of silence after interval I 1 . Figure 2 shows one period of the segregated cycle skipping MAIN state SD 1 . We notice that Definition 4.1 for k 1 introduces two active tone intervals of silence (I 2 A and I 2 B , compared with SD in Figure 1E). The matrix form of each cycle skipping segregated state s k consists of 2 rows and n 2k + 2 three dimensional row vector elements per column that describe the dynamics of the two units during the active tone intervals I j A and I i B , for j 1, . . ., k + 1. By definition, this matrix is obtained from the matrix form W [V 1 , V 2 ] of state s defined in Eq. (4) by adding 2k matrices with 2 rows and 3 columns composed of vectors 0 between the matrices V 1 and V 2 . Therefore the corresponding cycle skipping state have period T (2k + 2)TR. For example: 2 | Matrix form of 2TR-periodic CONNECT states (* asymmetrical states). The existence regions of each state are located between the existence regions of specific MAIN states ( Figure 1F). CONNECT states were named to indicate the branches of such MAIN states interspersed by the letter c. For example, the existence region of state APcAS is located between states AP and AS. Names containing the letter Z indicate that they connect branches of zero solutions (which cannot appear in the system for our parameter constraints).  where the A and B active tone intervals are indicated using the matrix form. We note that for k 0 the cycle skipping states correspond to the original segregated state (i.e., SD 0 SD). Theorem 3. There are no cycle skipping segregated states of the 2T R-periodic states IB, ID, AP, AS or ASD.

ZcS
Proof. We assume that these states are SHORT (the proof for LONG states is analogous). All of these states have a general matrix form given by.
for some binary values x, y, z, w, p and q. Application of Definition 4.1 would give the following matrix form for the corresponding kth-order cycle skipping state: Where 0 . . . 0 is a sequence of 2k repetition of the vector 0 . We highlighted in colors some key entries of the matrix form. Since the green entry is 1 the B unit is ON and instantaneously turns OFF at time β B k+1 . Therefore the delayed synaptic variable s B (t − D) starts to decay monotonically from time β B k+1 +D. Due to the periodicity of the solution and the red entry being 1 the total input to the A unit turns ON at time 0 and its total input is j A (0) −bL − 1 +c≥θ. Moreover the A unit turns OFF at time TD. Therefore the first 0 in the blue entry shown in this matrix gives the condition −bL + 1 +c<θ, which is absurd since L + 1 ≤L − 1 . This proves that Def 1 cannot be applied to states AP, AS, ASD, ID and IB, which conludes the PROOF.
We now introduce a new class of cycle skipping states obtained from integrated states.
Definition 4.2. A cycle skipping integrated state i k of order k ∈ N is the periodic solution of system (1) obtained from a integrated state i ∈ I by adding 2kTR periods of silence after intervals I 1 and I 2 . Figure 2 shows an example dynamics of the integrated cycle skipping MAIN state ID 1 . In this case the Definition 4.2 for k 1 introduces four active tone intervals of silence (I 1 ). Definition 4.2 can be applied to segregated states s ∈ S (not just to integrated states). The resulting cycle skipping state s k is equal to the segregated cycle skipping state S 2k defined by 4.1. Therefore, Definition 4.1 include segregated cycle skipping states and we can omit an equivalent of Definition 4.2 for segregated states.
The matrix form of any cycle skipping integrated state i k has 2 rows and n 4k + 2 three dimensional vector elements per column that describe the dynamics of the two units during the active tone intervals I i A and I i B , for i 1, . . ., 2k + 1. Similar to cycle skipping segregated states, this matrix is obtained from V [V 1 , V 2 ] by adding k matrices formed by 2 row vectors 0 after matrices V 1 and V 2 . For example: where the A and B active tone intervals are indicated using the matrix form. Therefore corresponding cycle skipping states have period T (4k + 2)TR. We note that for k 0 the cycle skipping states correspond to the original integrated state (i.e., ID 0 ID). Proof. In Theorem 3 we showed that cycle skipping segregated states AP, AS and ASD cannot be defined by Definition 4.1. Thus, by contradiction suppose that there exist cycle skipping integrated states of AP, AS and ASD defined by Definition 4.2. The matrix forms of each of these cycles skipping states extends the matrix form of AP, AS and ASD shown in Table 1 and can written in the general form: . Therefore, for the periodicity of these states the delayed synaptic variable s B (t − D) evaluated at time α 1 A 0 (red entry) is equal to L − 2k+1 and evaluated at time β k+1 A 2kTR+TD (blue entry) is equal to L + 4k+1 . Since the red entry is 1 we have that c−bL − 2k+1 ≥θ and since the last 0 in the blue entry is 0 we have that c−bL + 4k+1 <θ, which is absurd since L − 2k+1 ≥L 4k+1 . Lastly, we prove that Definitions 4.1 and 4.2 cannot be applied to states S and INT. Suppose that one these definitions can be applied, which we will show to be absurd. The extended matrix form of the cycle skipping states resulting from Definitions 4.1 and 4.2 have respectively the general forms the blue entry in the matrix forms above. This is absurd since the values of these entries are different, which concludes the PROOF.
We have therefore proven that the only possible cycle skipping MAIN states are obtained from segregated states s ∈ S or integrated states i ∈ I by adding 2TR-multiple silent periods. It is possible to derive the existence conditions of cycle skipping states in dependence on the model parameters. Table 3 shows the matrix form and the existence conditions of all cycle skipping segregated and integrated MAIN states. We now prove how these existence conditions are obtained in one example: ID k (shown in Figure 2 right for k 1). The proofs for the remaining cases are analogous.
Theorem 5. The existence conditions for state ID k are Proof. For this state the A and B unit turn and remain ON during the interval I 1 A and I k+1 B on the fast time scale, as shown by its matrix form in Table 3 and are OFF during all the remaining active tone intervals. The delayed synaptic variables s A (t − D) and s B (t − D) therefore are ON in the intervals I 1 [D, TD+D] and I 2 [(2k+1)TR +D, (2k+1)TR+ TD+D] (marked by the green shaded area in Figure 2). At the end of these intervals s A (t − D) and s B (t − D) monotonically decay following the slow subsystem s′ − s/τ i . One can therefore calculate their values at any point. The total input strength to the A unit at time β k+1 A is therefore equal to c−bL + 2k <θ, where the inequality holds because the first row of the a matrix form below I k+1 A is zero. Due to the monotonic decay of the delayed synaptic variables all the entries of the matrix below I is equal to a−bL + 2k+1 . Thus the condition a−bL + 2k+1 <θ guarantees that point (1, 1) is not an equilibrium for the fast subsystem at these times (see [6] for further details), which completes the PROOF.
Remark 4.1 (Analysis of SHORT CONNECT and LONG states). Cycle skipping SHORT CONNECT states can be defined using Definitions 4.1 and 4.2 by adding periods of silence multiple of 2T R after I 1 and/or I 2 to the states in Table 2. More precisely, Definition 4.1 can be applied to states ZcS, ScSD and Definition 4.2 can be applied to states ZcI, APcI. It can be proven that these definition cannot be applied to all remaining CONNECT states (ZcAP, ScAS, SDcAS, ZcAS and APcAS). This proof is analogous to the one of Theorem 4. In Supplementary Appendix A1 we derived the existence conditions of all possible cycle skipping SHORT CONNECT states using a similar approach to that of Theorem 5.

Analysis of the Remaining Classes of Cycle Skipping States
The analysis of LONG MAIN and CONNECT states follows analogously to the SHORT case. One caveat of this analysis is the need to extend the matrix form for short states during an active tone interval R (3) by adding one binary entry at the end of both rows in this matrix (see [6]). This additional entries are either both 0 or both 1. If they are both 1 both units are active after the offset of the interval R (LONG condition), if they are both 0 the units turn OFF at the end of the interval R. The possible LONG cycle skipping states and their existence conditions are shown in Supplementary Appendix A2 (for LONG MAIN cycle skipping states) and Supplementary Appendix A3 (for LONG CONNECT cycle skipping states). We note that these conditions depend on both quantities L ± k and R ± k defined in Eq. (5). Using the analysis above we can group all cycle skipping states of order k in two classes. The class S k of all (2k + 2)TR-periodic segregated cycle skipping states having the same number of turning ON times as state s k defined by (4.1): S k {SB k , SD k , ZcS k , ScSD k , SDL k , SL k , ScSDL k , ZcSL k }, and the class I k of all (4k + 2)TR-periodic integrated cycle skipping states having the same number of turning ON times as state i k defined by (4.2):    Figure 3A shows an example (for other parameters fixed) of the regions of existence for grouped states S k and I k when varying c and η.

Organization in the Parameter Space
As c → θ cycle skipping solutions S k and I k accumulate in infinite cascades with period tending to infinity and are separated by exclusive transition boundaries ( Figure 3B, there is no bistability between these states). The existence regions for each class S k and I k form a repeating pattern as k → ∞ with c decreasing. We note that at k 0 these regions correspond to the regions of existence of 2TR-periodic MAIN and CONNECT states (compare with Figure 1F; same parameter values). Figures 3C,D show the zoom of regions in panel A highlighting the features of subdivisions of existence conditions for states in S k and I k . Panel C shows an example of a pattern formed by SHORT states in S 1 and I 1 . As k increases (c decreases) this pattern continues for states in groups S k and I k and with decreasing (increasing) width (height). The lower borders separating ZcI k (ZcS k ) and APcID k (ScSD k ) from other non cycle skipping states shifts downwards and reaches the axis η 0. At the intersection with η 0 existence regions for cycle skipping states change to match a pattern formed by LONG states. Figure 3D shows an example of this pattern for k 7, a pattern that repeats as c → θ.

Parameter Region Relationships
We now show that the width ratio of the region of existence of successive cycle skipping states in group S k converges to a constant in the limit as k → ∞. A similar proof holds for states in I k and is therefore omitted. The existence region of each state S k is bounded horizontally by vertical lines ( Figure 3A). The right vertical boundary c c R k separates states SD k and ZcI k , while the left vertical boundary c c L k separates states ZcS k and ID k+1 ( Figure 3C). From the existence conditions of these states we can derive these lines analytically and obtain the expressions c R k bL 2k−1 +θ and c L k bL 2k +θ. The width of group S k is therefore given by <S k > c R k −c L k . It is straightforward to obtain the limit of the width between successive parameter intervals:

SWITCHING STATES
The remaining two classes of states considered in our study consist of SHORT states with a period that is multiple of TR.
Following the approach outlined in the previous section we define these states using their matrix form and use the matrix form to derive their existence conditions. To simplify the notation we introduce vector 0 (0, 0, 0) and 1 (1, 1, 1). We also define two positive integers values m, k ∈ N such that k ≥ 1, m ≥ 2 even. Definition 5.1. Let m be even. The Q-switching state MQ m,k is the periodic SHORT MAIN state with period (m + 2k)TR described by the matrix form: We note that this state is characterised by silence in m repetitions of the half period TR (blue boxes) and k repetitions of the B unit's turning ON at every other tone (red boxes). A time history example of one Q-switching state Q 2,3 is shown in Figure 4. In the next theorem we used an similar analysis to the one in the previous section to prove the existence conditions for MQ m,k states.
Hence the A unit is also OFF during all of these intervals. The condition d − bL − m+2k ≥θ guarantees that the B unit turns ON at the onset the interval I 1 A , and the condition a + c − bL − 1 ≥θ guarantees that the A unit turns ON with some small delay from the onset of the interval I 1 A . Lastly, the condition a − bL + 1 <θ guarantees that the state is SHORT by making sure that the units turn OFF at the offset of the interval I 1 A ((1, 1) is not an equilibrium for the fast subsystem, see [6] for more details), which completes the PROOF. In Supplementary Appendix A4, 5 we provide a definition and analyses of the existence conditions for SHORT CONNECT Qswitching states of order m and k, which have the same period as MQ m,k . We classify these states as: 1) simple CONNECT Q-switching states (cQ m,k , dQ m,k and gQ m,k ) for which one unit turns ON during only one interval at a delayed time from the interval's onset, and 2) complex CONNECT Q-switching states (cdQ m,k , cgQ m,k , dgQ m,k and cdgQ m,k ) with delayed turning ON of the units across multiple intervals. These classes complete CONNECT Q-switching states. For a fixed k and m, CONNECT Q-switching states preserve the number and the order of the units' turning ON times during each active tone interval from MQ m,k , except for delaying at least one unit's turning ON time during at least one active tone interval.
Next we proceed to define the final class of states considered in this paper: the W-switching states.
Definition 5.2. Let m ≥ 0 be odd. A W-switching state MW m,k is the periodic SHORT MAIN state with period 2 (m + 2k)TR described by the matrix form: where p(k) k + (m + 1)/2. Figure 4 shows the dynamics of state W 1,2 . We note that this state demonstrates switching between unit A being active at every other tone whilst unit B is inactive, and unit B being active at every other tone whilst unit A is active. With respect to auditory streaming this state may be interpreted as a segregated percept with periodic switching between the A tone in the foreground and the B tone in the foreground (see Discussion). We define SHORT CONNECT W-switching states of order m and k as those states that preserve the number and the order of the units' turning ON times during each active tone interval of state MW m,k , and that delay a unit's turning ON time during one or more active tone intervals. The matrix form and existence conditions for all W-switching states are similar to those of Q-switching states and are presented in Supplementary Appendix A6. Also in the case of CONNECT Q-switching states we differentiate between simple CONNECT W-switching states (cW m,k , dW m,k and gW m,k ) and complex CONNECT W-switching states (cdW m,k , cgW m,k , dgW m,k and cdgW m,k ) based on the number of active tone intervals for which the units' turning ON times are delayed.
From the analysis of Q-and W-switching states we may group all switching states of order m and k based on having the same number and order of the units' turning ON times during each active tone intervals. Specifically, we have the class of SHORT Q-switching states formed by all Q-switching states of order m and k: Q m,k {MQ m,k , cQ m,k , dQ m,k , gQ m,k , cdQ m,k , cgQ m,k , dgQ m,k , cdgQ m,k }, and the class of SHORT W-switching states formed by all Wswitching states of order m and k: W m,k {MW m,k , cW m,k , dW m,k , gW m,k , cdW m,k , cgW m,k , dgW m,k , cdgW m,k }. Figure 5A shows example of the dependence of the regions of existence of grouped states Q m,k and W m,k when varying parameters c and η (same parameters as Figures 1F, 3). For decreasing values of η switching groups Q m,k and W m,k accumulate in cascades with increasing period so that k → ∞ ( Figure 5B) until they reach the curve η θ/c, which is the upper boundary of state S shown in Figure 1F. By definition d ηc. Therefore, this accumulation occurs for decreasing values of the lateral input strength d for d → θ. We note that there is no bistability between these states (i.e., coexistence between any pair of states in Q m,k ). Figures 5C,D show details from panel A. These panels show example patterns of the existence regions for states in groups Q 2,4 and W 1,3 . These patterns repeat at different sizes as d → θ. In this limit the height of the intervals of existence of groups Q m,k and W m,k decreases to zero.
Remark 5.1. We note that the some Q-and W-switching do not exist for the selected parameter values in Figure 5 (for example states MQ m,k ). However, we numerically checked that these states exist in other parameter regions (see Figure 4 for state MQ 2,3 ).
We now show that the ratio of the height between successive W-switching states converges to a constant. An analogous proof holds for Q-switching states and is therefore omitted. The transition boundary between groups W m,k and W m,k−1 shown in Figure 5A is marked by the upper existence boundary of states MW m,k , dW m,k and dgW m,k expressed as a function of c as η k (c) (bL + m+2(k−1) +θ)/c, where L + j e (jTR−TD −D)/τ i . Therefore the height of the region of existence for the group W m,k is given by η k+1 (c) − η k (c). It is straightforward to find the limit of the height between successive tone intervals: This constant is independent from c and equal to the width ratio of successive cycle skipping states Eq. (8).

COMPUTATIONAL ANALYSIS IN A SMOOTH AND CONTINUOUS SYSTEM
In this section we extend the analysis of switching states W m,k using a continuous version of model (1), described by the following system of six ordinary differential equations (ODEs): where the gain function S is defined as a sigmoid S σ (z) 1/(1 + e −σ(x−θ) ), which approximates a Heaviside function for large σ.
The interpretations of parameters a, b, θ, τ, τ i and of variables u A , u B , s A , s B are the same as the ones described in Section 2 for model (1). The new parameters α x and β x represent the rate of activation and inactivation of the variables x A and x B . The parameters σ 180 and λ 20 are the selected slopes of the sigmoid. System (12) differs from the discontinuous delayed model studied in the previous sections because the delay in the synaptic variable s A (s B ) is introduced by a second variable x A (x B ) instead of the fixed amount D, following [19] as motivated by indirect synapses. For example, if u A turns ON then x A will activate when u A crosses the threshold θ. In this case the synaptic variable s A activates when x A crosses the threshold θ. The delays produced by the indirect synapses depend on the rates α x and β x . Preliminary simulations revealed that the variation of these two parameters allows for only small delays to be produced whilst guaranteeing the convergence to 1 when variables s A and s B turn ON (delays measured to the peak of s A and s B , in the range 0-20 ms). Since our previous analysis assumes that D > TD we consider values of α x and β x that give a delay of approximately 15 ms and fixed TD 10 ms.
The input functions i A and i B are redefined as: where x sin (δt), y cos (δt), z sin (δ · TD)y − cos (δ · TD)x and δ π · PR. We recall that parameter η is a scaling of the input strength c so that the lateral inputs have amplitude d ηc. The new input function i A and i B are similar to the discontinuous inputs shown in Figure 1B but with smoothed, square waveform.
We simulate model (12) at fixed parameters and find a set of states W 1,k analogous to the ones shown in Figures 4, 5 for model (1). We numerically continue each branch of periodic orbit W 1,k for k 2, 3, 4, 5 varying the parameter η ( Figure 6). We find large regions of stability that accumulate with k increasing as η decreases (in agreement with Figure 5B). These regions of stability are non overlapping (no bistability) and are separated by unidentified bifurcation points (the continuation algorithm fails to converge at these points). We suspect that these bifurcation correspond to period-incrementing bifurcations.
Overall, this analysis confirms the existence and of Wswitching states and their organisation in a cascade for decreasing values of η using a smooth and continuous version of the original model (1). The numerical continuation of branches of W-switching states reveals discontinuous transition points that suggest bifurcation phenomena that cannot be detected using continuation software. A similar analysis of system (12) could be applied to study cycle skipping and Q-switching states.

DISCUSSION
In this paper we analysed a system of four periodically forced delay differential equation encoding rhythm perception. The system represents two neural populations mutually connected by fast excitation and slow delayed inhibition receiving square-wave sound inputs consisting of two alternating pure tones. We previously used singular perturbation techniques to classify all viable states 1:1 locked to the inputs and to define their existence conditions using a binary matrix representation [6]. Here we extended this analysis to study two classes of n:1 locked states called cycle skipping and switching states. Cycle skipping states were derived from 1:1 locked states by adding intervals of silence (i.e., without response) equal to a multiple of the input period. Switching states are defined by dynamic alternation between 1:1 locked states. We generalised the matrix form for these states to visualise their dynamics and derive their existence conditions. We then analysed the existence regions in the space of parameters (c, η), where c is the strength of the local inputs and η the scaling strength of the lateral input. We found that cycle skipping and switching states accumulate in cascades with diverging periods for decreasing values of c and η respectively, and that this accumulation derives from an infinite period-incrementing cascade of bifurcations. Lastly, we confirmed the analytical predictions using numerical integration and continuation in a smooth ODE version of the model.

Period Incrementing Cascades
Sequences of bifurcations leading to an increased period are a common features of dynamical systems, a period-doubling cascade being the most commonly observed. Other periodincreasing cascades have been most widely studied in nonsmooth maps: the period-adding scenario and periodincrementing scenario (see [35] for a recent review). In the period-incrementing scenario the period increments from a base value p 0 by a fixed amount Δp in integer n steps: in each successive increment of the cascade the period is p n p 0 + nΔp. In the period-adding scenario, the sequence of periods observed as a parameter is varied is more complicated, with a Cantor-set like structure [35]. Either scenario can be found in 1D piecewiselinear maps [31,36] (and their regularised, smooth counterparts [37]). Successive intervals typically exist in wedge-shaped regions converging to a point in the parameter plane (a so-called big bang bifurcation) and an infinite cascade of period-adding or periodincrementing transitions is found by following a radial path in parameter space around these organising centers in the parameter plane [31,36]. In some cases, bands of chaotic dynamics (like periodic windows found in the period-doubling route to chaos) can be found between parameter intervals with periodic states [31]. Period-incrementing and period-adding cascades are found in a broad range of applications in settings that include 1D or nD maps (non-smooth or smooth), and include continuous time dynamical systems with delays or with non-smooth characteristics [35][36][37]. Our discussion here will focus on a comparison of our findings with canonical examples [31,36] and with models from mathematical neuroscience (below). In our model we identified several cascades characteristic of period-incrementing sequences. In each successive parameter region the period increases in steps associated with the period of the model inputs. In the case of skipping states the sequence progresses with increasing period for decreasing values of the local input strength c. The boundaries of these regions are parallel in the (c, η)-plane (and thus do not converge at a big bang bifurcation). Interestingly there are two interleaved periodincrementing sequences associated with integrated-like states (both units respond together) and segregated-like states (only one unit responds). We further found that the ratio of width in c of successive intervals to be a parameter dependent constant e −2TR/τ i (rather than a universal constant as for the Feigenbaum constant in period-doubling sequences). The structure in the (c, η)-plane was similar for switching states but with the sequence occurring as η is decreased. Interestingly, the ratio of height in η of successive existence intervals of switching states is also e −2TR/τi . Noting that the existence regions for periodic states in our model are disjoint and contiguous, we did not find windows with chaotic dynamics. While the transition boundaries between successive skipping states are vertical lines in the space of parameters (c, η), the boundaries of switching states are curved. We note that by plotting the existence regions in the space of parameters (c, d) instead of (c, η) we would obtain horizontal transition boundaries for switching states and vertical transition boundaries for skipping states.

Cascades in Neural Models
Models of single neurons typically feature a time-scale separation with rapid activation and inactivation dynamics that generate spikes and slower recovery dynamics leading to a refractory period between spikes [12]. Methods from singularperturbation theory have been central to the development of our understanding of more complex spiking phenomena such as bursting [12,13]. Successive spikes appear on the bursting plateau as a parameter is varied through a complex series of homoclinic bifurcations [30,38,39]. However, the period does not change much [40] in these spike-adding sequences (therefore not classified as period-adding or period-incrementing). However, there are several examples of single-neuron models with nonsmooth spike resetting dynamics where such cascades have been reported. In [41], a piece-wise linear approximation of the Fitzhugh-Nagumo's nonlinearity was taken. The periodically forced system exhibits rich dynamics that includes a periodadding sequence and chaotic dynamics. An integrate and fire model (with spike reset) also features a period-adding sequence and chaotic dynamics [42]. More recently a similar model was reported to have a period-incrementing structure [43] as revealed through a map on the model's adaptation variable. In the present study we did not find evidence of chaotic dynamics. Indeed our analysis precludes any 'gaps' in the contiguous parameter intervals with cycle skipping and switching states where such dynamics are found in other models [42,43].
In canonical models of neural competition at the population level, timescale separation leads to relaxation oscillations in models of perceptual rivalry [8,44] (i.e., perceptual bistability, as for binocular rivalry or the necker cube). Mixed-mode oscillations can emerge in these models through a canard mechanism [15,45,46]. A cascade of non-local bifurcations leads to disjoint, contiguous parameter regions with branches that have increased period [46], though the period is not fixed on each segment in this smooth model. This class of models can Frontiers in Applied Mathematics and Statistics | www.frontiersin.org October 2021 | Volume 7 | Article 716288 exhibit a cascade sharing many features of a periodincrementing cascade: when periodic forcing is introduced a sequence of non-local bifurcations lead to a family of solutions, similar to the W-switching states reported here, with stepincrements in the period matching the period of the forcing [46]. Periodic inputs to competition networks have been considered in several other models [5,24,25,47,48], but these typically only report 1:1 locked states or oscillatory states modulated by forcing but not locked with the forcing. A model featuring a Heaviside activation function and stepfunction or sinusoidal inputs exhibited cycle-skipping states [25], much like the cycle skipping states found in the present study. Our study reveals richer dynamics than these earlier studies and demonstrates a link between the regions with the switching and cycle skipping states through a single, parameterdependent ratio for existence intervals.

Modelling Assumptions
Our model features several mechanisms that are known to lead to complex dynamics (in isolation, or in combination with each other): slow-fast timescale separation, delays, periodic inputs, a non-smooth (heaviside) activation function and non-smooth steps in the input profile. In Section 6 we explored whether the skipping and switching cascades persist when we relax the non-smooth characteristic of the activation function and the inputs. We demonstrated that a sequence of solutions selected from the skipping cascade persists and that, these non-smooth components in the model are not a requirement to find the period-incrementing cascades reported. We further explored the skipping states in a version of the model where the delays were introduced via synaptic coupling [19] in an extended ODE rather than DDE formulation. Tracking branch segments with numerical continuation allowed us to confirm that the skipping states persist in this case. Whilst it may be possible to produce switching states [46] and cycle skipping states [25] without delays, removing the delays in the present model dramatically simplified the dynamics. Indeed excitation and inhibition would instantaneously turn ON when the units turn ON. In the slow-fast limit there are only two possible states: a state where both units instantaneously turn and remain ON for each active tone interval (corresponding to IB or ID described in Section 3) or a state where unit A (B) turns ON for each A (B) active tone interval (corresponding to AP described in Section 3).
In summary the rich dynamics produced by our model result from an interaction between slow-fast mechanisms with delays and periodic forcing. The condition TD < D assumed in our study guarantees that the delayed synaptic variables of system 1 are monotonically decaying in each active tone interval. This simplifies the system under study by reducing the number of possible states and their existence conditions (see also [6]). A more detailed analysis of the case TD ≥ D leads to switching and cycle skipping states similar to those studied here (not shown). Under this condition for integrated and segregated cycle skipping states the units' turn ON for some active tone interval and must turn OFF after the delay D, instead of TD as in the case TD < D (see Figure 2). This leads to additional existence conditions for these states.
Overall, a similar organization of the parameter space with cascades of period adding bifurcations can be proven analogously under condition TD ≥ D. The latter condition may be relevant when modelling auditory perception, as the potential factors for generating delayed inhibition in brain circuits would most likely lead to short or moderate delays that could be less than TD.

Interpretation for Neural Encoding of Auditory Streams
The solutions found in the switching and skipping cascades each have a meaningful interpretation for the model. In the motivating auditory streaming paradigm, two sequences of interleaved tones can be perceived by human listeners in different configurations: integrated into a single stream, segregated with A in the foreground or segregated with B in the foreground [49]. Some experiments consider only bistable perception (integrated vs segregated [5,50]), whilst others consider tri-stability [49] or multi-stability [51], allowing for two different segregated states (A or B in the foreground). In the present model, switching states correspond to periodic states with alternations between an equal number of A-foreground and B-foreground cycles separated by a single integrated cycle during the transition. Tri-stability in perceptual competition has been investigated elsewhere in competition models [52].
In the perception of a regular beat, we can follow along (e.g., clapping, or tapping a foot) in time with different elements in the time-structure. Musical rhythms are perceived to have a pulse or basic beat in the range 0.5-4 Hz that is further subdivided by higher beat frequencies above 4 Hz, nevertheless, we are able to follow along at the pulse frequency. The skipping states can be linked to the encoding of neural rhythms at periods that span multiple cycles of the input (say tracking every 2nd or 3rd A tone). Such cycle skipping behaviour has been studied in m: k Arnold resonance tongues in a model of beat perception [53]. By contrast we find 1: k locked states appearing in a periodincrementing cascade in our model. We further find much larger regions of parameter space (i.e., range of c-values) with skipping states when the presentation rate is large. The model therefore predicts that the tendency to skip cycles is more likely when the inputs are too fast to follow.

CONCLUSION
In this study we address the formation of exotic dynamics via several complex mechanisms (see "Modelling assumptions" above). The architecture of the proposed model is probably not unique in producing such dynamics. We expect that similar network dynamics could be obtained by modification of the connectivity (i.e., by introducing global inputs and removing the mutual excitation). Whilst the presented methods should generalise to larger populations of neurons, this would result in significantly increased in complexity for the number of possible states and their existence criteria.
The states analysed here, and in our previous study [6], do not cover all the possible dynamics the model can produce. In preliminary simulations we detected states of infinite period (quasi-periodic), which remains a topic for future study. A further topic for future analysis is to modify the inputs by considering 4TR-periodic inputs representing the repetition of the triplet ABA-consisting of A and B pure tones and by a silent gap -, which has been widely used in behavioural experiments [54].
Our approach generalises to the study of complex oscillatory dynamics in neural population models featuring delays and different architectures, such as competition models in the visual or tactile domains. A model is currently being investigated to study perceptual competition with vibrotactile stimulation [55].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.