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–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–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–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 periodically-forced, 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 non-smooth 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 parameter-dependent 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 T R-periodic 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.

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 τ u ̇ A ( t ) = − u A ( t ) + H ( a u B ( t ) − b s B ( t − D ) + i A ( t ) ) , τ u ̇ B ( t ) = − u B ( t ) + H ( a u A ( t ) − b s A ( t − D ) + i B ( t ) ) , s ̇ A ( t ) = H ( u A ( t ) ) ( 1 − s A ( t ) ) / τ − s A ( t ) / τ i , s ̇ B ( t ) = H ( u B ( t ) ) ( 1 − s B ( t ) ) / τ − s B ( t ) / τ i . (1)

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: i A ( t ) = c ∑ j = 0 ∞ χ I A j ( t ) + d ∑ j = 0 ∞ χ I B j ( t ) i B ( t ) = d ∑ j = 0 ∞ χ I A j ( t ) + c ∑ j = 0 ∞ χ I B j ( t ) , (2) 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 T R representing the time between tone onsets. The intervals when A and B tones are on (active tone intervals) are given by I A j = [ α j A , β j A ] and I B j = [ α j B , β j B ] , respectively, and have boundaries given by: α A j = 2 ( j − 1 ) T R , β A j = 2 ( j − 1 ) T R + T D , α B j = ( 2 ( j − 1 ) + 1 ) T R , β B j = ( 2 ( j − 1 ) + 1 ) T R + T D , for j = 1, …, ∞. Let us name the set of active tone intervals R as Φ = { R ⊂ R : R = I j A or R = I j B , ∀ j ∈ N }

In this work we assume that TD < D and TD + D < T R, 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.1. The dynamics of the synaptic variable s _A (s _B ) is dictated by that unit A (B). Indeed s _A (s _B ) turns ON the fast time scale when A (B) is ON and slowly decays to zero when A (B) is OFF.

Remark 2.2 (Model symmetry). Let us consider the map swapping the A and B indexes in the variables of system (1) κ : v = ( u A , u B , s A , s B , i A , i B ) ↦ ( u B , u A , s B , s A , i B , i A ) and the application of the time shift T R. Since i _A (t + T R) = i _B (t) and i _B (t + T R) = 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)).

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].

Theorem 1. Consider an arbitrary interval 2TR time interval I containing intervals I A j = [ α j A , β j A ] and I B j = [ α j B , β j B ] and the subdivision of I in the intervals shown in Figure 1C . Any state of system (1) satisfy the following properties:

1. No unit can turn ON between successive active tone intervals (blue intervals)

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).

Definition 3.1 (MAIN and CONNECT states). A state (solution) of system (1) is:

• MAIN if ∀R ∈Φ, if ∃t* ∈ R turning ON time for A or B, then t* = min(R)

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):

• LONG if ∃ t ∈ R − I when both units are ON

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:

1. A (B) turns ON at time t* ∈ R ⇔ A (B) is ON ∀t ∈ (t*, β]

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 A 1 and I B 1 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.

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 A 1 and I 2 = I B 1 and add a period of 2TR-multiple of silence after I ₁ and/or I ₂ (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 S = { S B , S D } , I = { I B , I D } , 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: L k − = e − ( k T R − T D − D ) / τ i , L k + = e − ( k T R − D ) / τ i , R k − = e − ( k T R − 2 D ) / τ i , R k + = e − ( k T R + T D − 2 D ) / τ i . (5)

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 2kT R periods of silence after interval I ₁.

Figure 2 shows one period of the segregated cycle skipping MAIN state SD ₁. We notice that Definition 4.1 for k = 1 introduces two active tone intervals of silence ( I A 2 and I B 2 , 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 A j and I B i , for j = 1, …, k + 1. By definition, this matrix is obtained from the matrix form W = [V ₁, V ₂] 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 ₁ and V ₂. Therefore the corresponding cycle skipping state have period T = (2k + 2)TR. For example:

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 ₀ = SD).

Theorem 3. There are no cycle skipping segregated states of the 2 T R -periodic states IB , ID , AP , AS or ASD .

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.

V = 111 w p q x y z 111 ,

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: I A 1 I B 1 … I A k + 1 I B k + 1 111 0 ̲ … 0 ̲ x y z x y z 0 ̲ … 0 ̲ 111 (6) 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 β k + 1 B . Therefore the delayed synaptic variable s _B (t − D) starts to decay monotonically from time β k + 1 B + 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 ) = − b L 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 − b L 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 2kT R periods of silence after intervals I ₁ and I ₂.

Figure 2 shows an example dynamics of the integrated cycle skipping MAIN state ID ₁. In this case the Definition 4.2 for k = 1 introduces four active tone intervals of silence ( I B 1 , I A 2 , I A 3 and I B 3 ).

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 A i and I B i , for i = 1, …, 2k + 1. Similar to cycle skipping segregated states, this matrix is obtained from V = [V ₁, V ₂] by adding k matrices formed by 2 row vectors 0 ̲ after matrices V ₁ and V ₂. 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 ₀ = ID).

Theorem 4. There are no cycle skipping solutions of states AP , AS , ASD , S and INT .

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: I A 1 I B 1 … I A k + 1 I B k + 1 I A k + 2 … I B 2 k + 1 111 0 ̲ … 0 ̲ x y z 0 ̲ … 0 ̲ 000 0 ̲ … 0 ̲ 111 0 ̲ … 0 ̲ (7) where 0 ̲ … 0 ̲ is a sequence of 2k repetition of the vector 0 ̲ and x, y and z are binary values. The red entry in this matrix is 1, meaning that the total input to the A unit at time 0 is j _A (0) ≥ θ. Since the B unit is ON and turns OFF at the end time of interval I B k + 1 , i.e., time β B k + 1 = ( 2 k + 1 ) T R + T D and it remains OFF throughout the following active tone intervals I A k + 2 , … , I B 2 k + 1 . Therefore, for the periodicity of these states the delayed synaptic variable s _B (t − D) evaluated at time α A 1 = 0 (red entry) is equal to L 2 k + 1 − and evaluated at time β A k + 1 = 2 k T R + T D (blue entry) is equal to L 4 k + 1 + . Since the red entry is 1 we have that c − b L 2 k + 1 − ≥ θ and since the last 0 in the blue entry is 0 we have that c − b L 4 k + 1 + < θ , which is absurd since L 2 k + 1 − ≥ L 4 k + 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 I A 1 I B 1 … I A k + 1 I B k + 1 111 0 ̲ … 0 ̲ x y z 000 0 ̲ … 0 ̲ 000 and I A 1 I B 1 … I A k + 1 I B k + 1 I A k + 2 … I B 2 k + 1 111 0 ̲ … 0 ̲ x y z 0 ̲ … 0 ̲ 000 0 ̲ … 0 ̲ 000 0 ̲ … 0 ̲ where x = y = z = 0 for state S and x = y = z = 1 for state INT. In both cases the B unit is OFF at all times. Therefore the synaptic variable s _B is constant and equal to zero and the total inputs to unit A at the start of intervals I A 1 and I A k + 1 are both equal to c. These inputs correspond to the red entry and the first zero of 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 c − b L 2 k + < θ c − b L 2 k + 1 − ≥ θ d − b L 2 k + 1 − < θ a + d − b L 2 k + 1 − ≥ θ a − b L 2 k + 1 + < θ

Proof. For this state the A and B unit turn and remain ON during the interval I A 1 and I B k + 1 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 , T D + D ] and I 2 ̄ = [ ( 2 k + 1 ) T R + D , ( 2 k + 1 ) T R + T D + 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 β A k + 1 is therefore equal to c − b L 2 k + < θ , where the inequality holds because the first row of the a matrix form below I A k + 1 is zero. Due to the monotonic decay of the delayed synaptic variables all the entries of the matrix below I B 1 , … , I A k + 1 and I A k + 2 , … , I B 2 k + 1 are zero, since the total input is lower than c − b L 2 k + . The value of the entries x A 1 and x B 2 are 1 and indicate dynamics of unit A and B at times 0 and (2k + 1)TR, respectively. Calculating the value of the delayed synaptic variables at these times leads to the total inputs to both units c − b L 2 k + 1 − , and thus to the condition c − b L 2 k + 1 − ≥ θ . The value of the entries x B 1 and x A 2 are 0. Thus a similar reason leads to the condition d − b L 2 k + 1 − < θ . Since the values of y B 1 ( y A 2 ) is 1 one unit B (A) turns ON a small delay after A at time 0 ((2k + 1)TR) due to the excitation with strength a from unit A (B). This occurs because a + d − b L 2 k + 1 − ≥ θ . Lastly, the condition a − b L 2 k + 1 + < θ guarantees that the state is SHORT. Indeed the total input to the units at the end times of the active tone intervals I A 1 and I B k + 1 is equal to a − b L 2 k + 1 + . Thus the condition a − b L 2 k + 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 ₁ and/or I ₂ 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.

The remaining two classes of states considered in our study consist of SHORT states with a period that is multiple of T R. 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 T R (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.

Theorem 6. The existence conditions of state MQ _m,k are given by: c − b L m + < θ c − b L m + 1 − ≥ θ d − b L m + 2 k − 2 + < θ a + d − b L m + 1 + < θ d − b L m + 2 k − ≥ θ a + c − b L 1 − ≥ θ a − b L 1 + < θ . (9)

Proof. Due to the monotonic decay of the delayed synaptic variables in the m between intervals I B 1 and I A m / 2 + 1 (Theorem 1), the condition c − b L m + < θ guarantees that the units are OFF during these intervals (entries of the matrix form are 0 ̲ ). This condition also guarantees two properties. Firstly, that the A unit is OFF at the onset of the interval I A 1 , since its total input is c − b L 1 − < c − b L m + < θ (since m ≥ 2 from the definition of MQ _m,k ). Secondly, that the A unit is OFF during the A tone intervals I A m / 2 + 2 , … , I A m / 2 + k and the B tone intervals I B m / 2 + 2 , … , I B m / 2 + k . Indeed the total input to the A unit at the end of the first set of intervals is c − b L 1 + ≤ c − b L m + < θ , and A unit’s input at the end of the second set of intervals is d − b L 2 + ≤ c − b L m + < θ . The condition c − b L m + 1 − ≥ θ guarantees that the B unit turns and remains ON in the interval I B m / 2 + 1 . Due to the decay of s _A (t − D) this is valid also for intervals I B m / 2 + j , for j = 2, …, k. Overall, this leads the entries in the second row below each of the intervals I B m / 2 + j , for j = 1, …, k to be equal to 1 ̲ . The condition d − b L m + 2 k − 2 + < θ guarantees that the entry in the second row below the interval I A m / 2 + k is equal to 0 ̲ . The decay of the synaptic variable s _A (t − D) implies that all second row entries below the intervals I A j + k are 0 ̲ , for j = m/2 + 2, …, m/2 + k. The condition a + d − b L m + 1 + < θ guarantees that the A unit is OFF during the intervals I B m / 2 + 1 . Moreover, since the total input to the A unit at the end of the intervals I B m / 2 + j , for j = 2, …, k is equal to a + d − b L 2 + < a + d − b L m + 1 + , since m ≥ 2. Hence the A unit is also OFF during all of these intervals. The condition d − b L m + 2 k − ≥ θ guarantees that the B unit turns ON at the onset the interval I A 1 , and the condition a + c − b L 1 − ≥ θ guarantees that the A unit turns ON with some small delay from the onset of the interval I A 1 . Lastly, the condition a − b L 1 + < θ guarantees that the state is SHORT by making sure that the units turn OFF at the offset of the interval I A 1 ((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 Q-switching 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 = { M Q m , k , c Q m , k , d Q m , k , g Q m , k , c d Q m , k , c g Q m , k , d g Q m , k , c d g Q m , k } , and the class of SHORT W-switching states formed by all W-switching states of order m and k: W m , k = { M W m , k , c W m , k , d W m , k , g W m , k , c d W m , k , c g W m , k , d g W m , k , c d g W 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 ) = ( b L m + 2 ( k − 1 ) + + θ ) / c , where L j + = e ( j T R − T D − 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: lim k → ∞ η k + 1 ( c ) − η k ( c ) η k ( c ) − η k − 1 ( c ) = e − 2 T R / τ i .

This constant is independent from c and equal to the width ratio of successive cycle skipping states Eq. (8).

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): τ u ̇ A ( t ) = − u A + S σ ( a u B − b s B + c i A ( t ) ) , τ u ̇ B ( t ) = − u B + S σ ( a u A − b s A + c i B ( t ) ) , x ̇ A ( t ) = α x ( 1 − x A ) S λ ( u A ) / τ − β x x A , x ̇ B ( t ) = α x ( 1 − x B ) S λ ( u B ) / τ − β x x B , s ̇ A ( t ) = S λ ( x A ) ( 1 − s A ) / τ − s A / τ i , s ̇ B ( t ) = S λ ( x B ) ( 1 − s B ) / τ − s B / τ i , (10) 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: i A ( t ) = S σ , 0 ( x ) S σ , 0 ( z ) + η S σ , 0 ( − x ) S σ , 0 ( − z ) i B ( t ) = S σ , 0 ( − x ) S σ , 0 ( − z ) + η S σ , 0 ( x ) S σ , 0 ( z ) , (11) 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. To allow numerical continuation with periodical forcing in AUTO-07p [34] we append the following two ODEs to system (13), for which x = sin (δt) and y = cos (δt) is a solution: x ′ = x + δ y − x ( x 2 + y 2 ) , y ′ = − δ x + y − y ( x 2 + y 2 ) .

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 W-switching 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.

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.

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].