Order-indeterminant event-based maps for learning a beat

The process by which humans synchronize to a musical beat is believed to occur through error-correction where an individual's estimates of the period and phase of the beat time are iteratively adjusted to align with an external stimuli. Mathematically, error-correction can be described using a two-dimensional map where convergence to a fixed point corresponds to synchronizing to the beat. In this paper, we show how a neural system, called a beat generator, learns to adapt its oscillatory behaviour through error-correction to synchronize to an external periodic signal. We construct a two-dimensional event-based map which iteratively adjusts an internal parameter of the beat generator to speed up or slow down its oscillatory behaviour to bring it into synchrony with the periodic stimulus. The map is novel in that the order of events defining the map are not a priori known. Instead, the type of error-correction adjustment made at each iterate of the map is determined by a sequence of expected events. The map possesses a rich repertoire of dynamics, including periodic solutions and chaotic orbits.

The process by which humans synchronize to a musical beat is believed to occur through error-correction where an individual's estimates of the period and phase of the beat time are iteratively adjusted to align with an external stimuli. Mathematically, error-correction can be described using a two-dimensional map where convergence to a fixed point corresponds to synchronizing to the beat. In this paper, we show how a neural system, called a beat generator, learns to adapt its oscillatory behaviour through error-correction to synchronize to an external periodic signal. We construct a two-dimensional event-based map which iteratively adjusts an internal parameter of the beat generator to speed up or slow down its oscillatory behaviour to bring it into synchrony with the periodic stimulus. The map is novel in that the order of events defining the map are not a priori known. Instead, the type of error-correction adjustment made at each iterate of the map is determined by a sequence of expected events. The map possesses a rich repertoire of dynamics, including periodic solutions and chaotic orbits.
Music is an important part of human society. As humans, we have an innate ability to quickly recognise rhythmicity and reproduce it by moving our body to the music. To bring our movements into sync with the music, we make a series of adjustments to speed up or slow down in a process known as error-correction. The time between successive movements, such as claps, serves as an estimate of the beat period, and is compared to the actual time between musical beats. The exact timing of each clap compared to the beat time determines whether or not our clapping phase is aligned with the beat. Using these two comparisons we make a judgement on whether we need to speed up or slow down our clapping. In this paper, we explore how such an error-correction scheme could be implemented by a neuronal network. Our work focuses on the derivation and analysis of a two-dimensional map to describe how an error-correction scheme can bring a neural system into sync with a periodic external drive. Importantly, the map can describe situations in which the order of events changes. This is critical because we often aperiodically alternate between being too early and too late when learning a beat, i.e. on one cycle we may clap just before the beat, while on the next we clap after it. We call this property order-indeterminacy.

I. INTRODUCTION
Humans can easily recognize rhythmicity within speech and music which spans the range of 0.5 to 10 Hz 1 . The a) Electronic mail: aine.byrne@ucd.ie ability to discern and track a periodic structure in music is called beat perception. While a piece of music may be quite complex in terms of its rhythmic structure, experimental studies have shown that neurons in various parts of the brain exhibit oscillations in their voltage profiles that match the beat frequency 2,3 . A primary way to assess beat perception is through finger tapping experiments where participants tap at what they perceive to be the beat or along to an isochronous (evenly spaced in time) metronome 4,5 . Intertap intervals are compared to interbeat intervals to measure the participants ability to match the period. The exact timing of taps at each cycle is compared to beat times to measure the phase of tapping. These two measures, period and phase, have led to a set of error-correction models 6,7 . These iterative algorithms attempt to describe how humans use error measurements at each cycle to correct their tap times in subsequent cycles. In essence, the schemes defines a two-dimensional, event-based map.
Event-based maps refer to a class of dynamical systems where the set of dependent variables are updated at the occurrence of a particular event, such as a trajectory crossing a Poincaré section. Event-based maps arise in other contexts besides error-correction. In neural systems, maps based on spike timing are often derived. For example, using a phase response curve in the presence of a weak-coupling assumption between neurons is used to assess the existence and stability of phase-locked solutions 8,9 . In cardiac systems, maps have been used to study the repolarization of the left and right ventricles 10 . In the field of robotic movement, event-based maps are used to implement specific control strategies for stabilization of walking 11 .
In this work, we derive a map that originates from our recently developed mathematical model for how humans learn to generate a beat 12 . We proposed a biophysical framework whereby a neuronal oscillator, called a beat generator (BG), learns the period and phase of an exter-nal, isochronous tone sequence (S). We derived an errorcorrection scheme that adjusts a biophysical parameter on a cycle-by-cycle basis allowing the BG to synchronize its voltage spikes to the tone times. Updates to this parameter occur at BG spikes and stimulus tones, as such they are characterized as events, and form the basis of the event-based map. Our previous results 12 are largely numerical, but the simulation results give rise to set of mathematical questions that we shall address here.
In the ideal case, the learning process can be described in a 1-to-1 manner in which every BG spike is followed by an S tone and vice versa. This leads to a situation of monotone convergence to the synchronized solution in which the order of events is preserved. However, when tapping to a beat, our tap times typically jitter in a neighborhood of the actual stimulus tone onset times. For example, if at one cycle, the tap is too late, we might compensate by tapping at a faster rate. This may then result in the second tap occurring before the second tone. These situations can be described as order-indeterminant because at the moment of any one event, say a tap, it is not a priori known whether the next event will be a tone or a tap. The event itself depends on the adjustments made on the fly by the listener in an attempt to synchronize to the tone sequence. Modeling this situation mathematically is a challenge as one can no-longer make a 1-to-1 (order preserving) assumption about events. Instead, we derive a novel two-dimensional order-indeterminant event-based (OIEB) map that can predict whether and how the BG will synchronize with the stimulus, even when the order of events is not prescribed from cycleto-cycle. The main difficulty in deriving this map is that different adjustments are made at stimulus tones and BG spikes. Thus when the order it unknown, the sequence of updates is unknown. One of the important contributions of this paper is to show how to systematically overcome this obstacle to make the appropriate number and type of updates per cycle.
An interesting property of the OIEB map is that it is piecewise-smooth, with a discontinuity arising due to the update rule for phase learning. As such, our study falls within the category of low-dimensional discontinuous maps, which are seen in a variety of contexts. From a theoretical point of view, some of the earliest studies of this type involved piecewise linear one-dimensional maps 13 . Discontinuous maps arise more generally in switching systems, where the vector field changes discontinuously across a lower dimensional manifold 14 . Biological applications have provided an abundance of examples. In particular, systems defined on periodic phase spaces such as a circle or torus lead to discontinuous maps. Applications of such maps include cardiac dynamics 15 , circadian rhythms 16 , and sleep rhythms 17 , to name just a few.
The organization of this paper is based on deriving the OIEB map from a series of intermediate steps to build maps of increasing complexity. The simplest such map is one-dimensional and accounts just for period learning. We will show the conditions under which this onedimensional map has a stable fixed point which corresponds to a learned period and then describe how the loss of stability leads to higher periodicity orbits. The next level of complexity involves phase learning, leading to the two-dimensional order preserving map. Here, we conduct a linear stability analysis for the synchronized (learned) solution. We demonstrate that while there are regions of parameter space in which the synchronized solution is stable, only part of the parameter space corresponds to order preserving solutions. This leads us to derive the OIEB map where information from the stability analysis in the order preserving case will be used to explain the observed dynamics of the OIEB map. This map possesses a rich set of dynamics that include not just the existence of fixed points, but also various periodic orbits, as well as chaotic solutions.

II. MODEL FOR BEAT GENERATION
We briefly introduce the important aspects of our beat generation model 12 . We first define a periodic stimulus tone sequence with an initial tone time, t 0 and interstimulus onset interval of T stim , such that tones occur at t = t 0 + kT stim , where k = 0, 1, 2 . . .. We then define a beat generator (BG) as a neuronal oscillator with easily identifiable spike times whose interspike interval can be controlled by parameters. There are a variety of biophysical models based on the Hodgkin-Huxley formalism that we could choose. Here, we utilize perhaps the simplest, the leaky integrate-and-fire (LIF ) model, to allow us to focus on the mathematical results. The LIF is given by with the reset condition v → 0, when v crosses the firing threshold, set (without loss of generality) at 1. At the moment of reset, the BG is said to exhibit a spike. The variable v represents voltage, the parameter I provides an external drive, and τ is a time constant. For I < 1, the system has a stable fixed point at v = I. The discontinuous reset condition leads to a bifurcation at I = 1, and the system exhibits oscillatory behaviour for I > 1. The period of these oscillations is given by which is a one-to-one invertible function on the domain I > 1, with inverse Thus each value of the period results from a unique value of the drive and vice versa. In the map, the value of I is updated at each BG spike and at each stimulus tone time, until the period and spike times of the BG exactly match those of the stimulus. Hence, I is a variable of the two-dimensional map. To update I we must keep track of the time between consecutive BG spikes (period) and the time between the BG spike and stimulus tone (phase). The phase of the BG φ is the second variable of the map. Phase is updated at every tone time. Thus, the two-dimensional map will exhibit a type of asynchronous updating where one of the variables, I, is updated more frequently, than the other, φ.

III. ORDER PRESERVING 2-DIMENSIONAL MAP
We begin by deriving the order preserving map under the assumption of 1-to-1 firing. Namely, every spike of the BG is followed by a stimulus tone and vice versa. The goal is to show the process by which the BG learns the exact spike times of the stimulus as well as the interonset interval between tones. To do so, we introduce two learning rules, a period correction rule and a phase correction rule. Both rules update the value of I, which we treat as a variables of the map. Period correction updates are done at every spike of the BG, while phase correction updates are performed at every stimulus tone time. In Section IV, we will extend these results to derive the OIEB map.

A. Period correction leads to a one-dimensional map
Consider a periodic stimulus oscillating with a period T stim with initial spike at t 0 = 0. For the BG to match the stimulus period, it must learn the value of I associated with that period. We define the following onedimensional, period correction map, which iteratively increases (decreases) I to decrease (increase) the oscillatory period of the BG, until it is the same as the stimulus period T stim . The parameter δ T is the strength of the period correction rule. Updates are associated with the spike times of the BG. Namely, once the update to I j+1 is made, the next interspike interval is uniquely determined by this new value. Using (2) we rewrite (4) as and The map has a unique fixed point at I * = 1/(1 − exp(−T stim /τ )), whose stability is determined using f ′ (I * ) = 1−δ T τ /(I * (I * −1)). The fixed point is stable for |f ′ (I * )| < 1, and, as such, δ T < 2I * (I * −1)/τ . Given that T is a monotonically decreasing function of I, it follows that to converge to a stable fixed point at larger periods, smaller δ T is required. The graph of f (I) has a vertical asymptote at I = 1, a local min at I = (1+ √ 1 + 4τ δ T )/2. It is concave up and in the limit as I → ∞ has a slant asymptote at I − δ T T . The slope of f (I) is less than 1 for finite I and tends to 1 as I → ∞. Thus the graph intersects the diagonal at exactly one point, corresponding to the unique fixed point. As T stim increases, the map shifts down in the I-f (I) plane ( Fig. 1(a)), and the unique fixed point becomes unstable when the slope at the intersection decreases below -1. For a fixed period, increasing δ T decreases the slope of the curve at the fixed point ( Fig. 1(b)). Combining these two findings, one notes that there are privileged parameter pairs that lead to fast convergence. For example, by choosing a δ T = 0.005, the slope at the fixed point associated with the T stim = 500 ms map is zero. For this same value of δ T , the slope of the T stim = 250 ms is positive, but less than 1. Thus a tempo change from T stim = 250 to 500 ms will converge much faster than a change from 500 to 250 ms ( Fig. 1(c)).
The fixed point remains stable so long as f ′ (I * ) > −1. The δ T value for which stability is lost can be calculated by solving f ′ (I * ) = −1, for δ T (Fig. 2). Convergence to the fixed point is fastest when f ′ (I * ) = 0. Hence, solving f ′ (I * ) = 0 for δ T , gives the optimal δ T value for a given period (dashed line). Notice that these values are highly dependent on the stimulus period, with smaller periods allowing for significantly larger δ T values consistent with what is shown in Fig. 1. The fixed point loses stability when f (I * ) < −1 and in that case it is possible to obtain periodic and bounded oscillatory solutions of the map provided that the local minimum of f (I) is greater than 1 (blue region). This necessary condition on the minimum of f (I) guarantees that iterates of the map fall in the domain of the map, (1, ∞). Periodic solutions are obtained in the usual manner, namely by searching for fixed points of suitable composition of the map with itself, e.g. period-2 solutions arise as solutions of f 2 (I) = f (f (I)) = I. For parameter values that lie in the blue region we were able to find period-2, 4 etc points. Moreover we found values of parameters at which the solution does not converge for any initial value, indicative of chaotic behavior. To illustrate the existence of chaos for this map, consider the period doubling route to chaos 18 . We calculated the ratios of the differences of parameters of successive period doubling bifurcation values (period-1 to period-2, period-2 to period-4 , etc). These ratios turn out to be approximately F 3 = 4.328, F 4 = 4.619 and F 5 = 4.655. Thus the period doubling behavior suggests that the value F n converges to the value F ≈ 4.669, the so called Feigenbaum ratio 19 , giving strong indication that the one-dimensional period correction map does exhibit chaos. The findings presented here, will be instructive when we consider the dynamics of the OIEB map in later sections.

B. Phase correction leads to a 2-dimensional map
The period rule brings the stimulus and BG periods into alignment, but pays no regard to the phase of the BG events relative to the stimulus events. As such, an additional rule is required, a phase correction rule. This rule also targets the drive I in order to speed up/slow down the BG oscillations. It works in concert with the period correction rule, sequentially making adjustments to I. We define the phase of the BG φ as the time since the last BG event divided by the stimulus period. If φ < 0.5 when a stimulus event occurs, the BG is interpreted to be firing before the tone and needs to be slowed down, decreasing both I and φ. If φ > 0.5, the opposite interpretation is taken and the BG should speed up, increasing I and φ. The maximal correction occurs close to 0.5. The phase update rule occurs at every stimulus event and will update the I value by an amount sgn(x) is the sign function and δ φ > 0 is the strength of the update rule.
The two-dimensional map requires that the I value is updated at both BG spikes (period correction) and S tones (phase correction). The phase φ is only updated once per cycle, at a BG spike. A cycle is defined as an oscillatory cycle of the BG, the nth cycle begins at the nth BG spike and ends at the n+ 1th spike. A complete cycle of the 2-dimensional map takes (I n , φ n ) to (I n+1 , φ n+1 ) through a sequence of intermediate steps. At the S tone, phase correction is applied The BG voltage (1) evolves with this new value I temp until its next spike. At the BG spike, I is updated again using the period rule. To apply the period rule, we must calculate the BG cycle period T n . This cycle period is comprised of two parts, the time between the nth BG spike and the nth tone, and the time between the nth tone and the n + 1th BG spike, where v(I, φ) = I 1 − e −Tstimφ/τ .
and the second term in (8) is obtained by integrating (1) from the initial value v(I n , φ n ) to 1. Then the period correction is given by where . At the (n + 1)th BG spike, we also update the phase for the next cycle, The system can then be written as a two-dimensional map, with updates to I and φ on every iteration defined by the paired event of an S tone and a BG spike, The functions F 1 and F 2 are given by A fixed point of the map satisfies the algebraic conditions I = F 1 (I, φ) and φ = F 2 (I, φ). For T stim = 500, we graph the surfaces F 1 (Iφ) and F 2 (I, φ) in separate three dimensional spaces (Fig. 3(a) and (b)), and examine their intersection with the two-dimensional planes z = I and z = φ, respectively. The ensuing curves of intersection are then projected onto the φ − I domain (Fig. 3(c)) and their intersection yields the fixed points of the map. As expected, the fixed points are located at (I * (T stim ), 0) and (I * (T stim ), 1). Though the fixed points at φ = 0 and φ = 1 correspond to the same synchronous solution, we will show later that the stability properties of each differ. Note that φ = 0.5 constitutes a line of discontinuities of the map, due to the sgn function that appears in the phase update rule. This discontinuity is clearly seen in the projections onto the φ − I plane (Fig. 3C). The blue curve separates regions of φ − I space where the surface z = F 2 (I, φ) lies above the diagonal plane z = φ or below it. In the projection the region above the blue curve represents where the surface lies above the diagonal plane. Similarly, the region above the red curve is where the surface z = F 1 (I, φ) lies below the diagonal plane z = I.
It is useful to recast the φ−I phase plane by placing the origin at I = I * and at the dual values φ = 0 and 1 (Fig.  4A), as both phase values correspond to the synchronized solution. Thus the entire vertical axis corresponds to both values φ = 0 and 1 and, as such, the BG spike and stimulus tone occurring at the same time. The horizontal axis is ordered to show how the phase changes, keeping in mind the dual role of the origin. In the left half-plane, phase decreases (left to right) toward the origin. in the right half-plane, phase increases (right to left) moving toward the origin. This recast phase plane allows us to easily identify how I and φ should be updated in each region, and in which direction the iterates should move. Namely, in the upper right first quadrant (Q1), the phase indicates that the BG spiked after the stimulus tone, and is thus late. As a result the phase rule will try to speed up the BG. However, the current value of I n is too large, thus the BG is too fast and will have an interspike interval less than T stim . Thus the period rule will tend to slow down the BG. For example, the iterates of the map with initial value (φ 0 , I 0 ) = (0.75, 2.62) (green) systematically decrease the I n value to slow the BG down, while simultaneously increasing the phase towards the value 1 (i.e. convergence to the origin). Note that the iterates quickly converge to the red nullcline along which I n+1 = I n . In a vicinity of this nullcline, the strengths of the phase and period rules are relatively balanced allowing the convergence towards synchrony to be monotone in phase. The initial condition (φ 0 , I 0 ) = (0.25, 2.47) (purple) lying in Q3 (lower left) corresponds to the BG being initially too slow, but also too early. Thus the period rule needs to increase the I n value, to speed up the BG, while the phase rule seeks to slow it down to decrease the phase. Again, the two-step updates eventually come into balance as can be seen the iterates move back and forth between Q2 and Q3 for a number of iterates, before converging towards the origin from Q2 (upper left) along the red nullcline.
The time courses for the two examples clearly display how the system convergences towards the synchonized solution in an order preserving manner (Fig. 4(a) and (c)). Initial condition lying in Q3 shows a systematic phase delay in its time course (Fig. 4(b)); the BG spike time (purple) moves towards the stimulus tone time (black) that follow it. Note that due to the way phase is defined, a delay implies that the phase decreases. If the initial condition lies in Q1, then it is in the basin of attraction of the φ = 1 fixed point, and the BG spike times systematically phase advance (Fig. 4(c)); the BG spike time (green) moves towards the stimulus tone time (black) that precedes it.

C. Assessing stability of fixed points through linearization
As demonstrated in Fig. 4 4. System can phase delay or phase advance during convergence to the synchronous solution. (a) Recast phase plane with dual origin (φ = 0 and 1) corresponding to the synchronized solution. Note the non-standard ordering of the phase along the horizontal axis. Iterates for two different initial conditions show that the phase is systematically decreased to φ = 0 when the BG is too early (purple) and increased toward φ = 1 when the BG is too late (green). (b) BG spike times (purple) and stimulus tone times (black) for the sequence of iterates that converges to φ = 0 via phase delay. (c) BG spike times (green) and stimulus tone times (black) corresponding to the iterates that converges to φ = 1 by phase advance. Parameter values given in Appendix A. φ = 1 fixed point. The convergence properties depend on the learning rule parameters δ T and δ φ as well as T stim . To better understand how the dynamics depend on these parameters, we assess the stability of the fixed points of the two-dimensional map by computing the eigenvalues of the Jacobian matrix, where (I * , φ * ) denotes a fixed point and g(I) = dT (I)/dI. The eigenvalues of J are given by A fixed point of the map is stable if both eigenvalues at the linearization lie inside the unit circle. We computed the stability boundaries for a fixed stimulus period T stim in the δ T −δ φ parameter plane (Fig. 5) by solving |λ| = 1 (17). This condition is met when an eigenvalue is real and equals ±1 or is complex with magnitude equal to 1. Although the fixed points at φ * = 0 and φ * = 1 correspond to the synchronous solution, they have different stability characteristics. The solid green lines correspond to the fixed point φ * = 1 for T stim = 500. The line emanating from the origin is where the eigenvalues are complex with magnitude equal to one; the other when one of the eigenvalues equals -1. The solid purple line correspond to φ * = 0 for T stim = 500 when an eigenvalue equals -1. The shaded green (purple) region correspond to parameter values for which the φ = 1 (φ = 0) fixed point is stable. It is straightforward to show that dg/dI > 0 which implies that the stability boundary curves for λ = −1 shift to the right and for |λ| = 1 have smaller slope as T stim is decreased. In turn, this implies that there is a larger set of parameters for which the φ = 0 and 1 fixed points are both stable when T stim is decreased (dashed lines in Fig. 5).
Contrary to convention, stability of the fixed points is not sufficient to determine whether iterates of the map will actually converge to them. The reason is that the 2-dimensional map is built under the order preserving assumption, namely that every BG spike is followed by a stimulus tone. When one of the fixed points, say φ = 1, is a stable spiral point, the iterates of the map spiral in towards φ = 1, but when they cross the vertical axis there is a switch to convergence towards the φ = 0 fixed point. Such a switch corresponds to two consecutive BG spikes and hence, a violation of the 1-to-1 assumption of the map. In the next section we address how to handle these types of situations.

IV. ORDER-INDETERMINANT 2-DIMENSIONAL MAP
As suggested in the previous section, convergence to the synchronous solution is not always order preserving. For example, if the BG spikes just before the stimulus tone, the phase correction may decrease I too much and slow the BG down such that it does not fire again until after the next stimulus tone (Fig. 6(a); cycle 1 and 6). Similarly, if the period rule over-corrects when speeding up the BG, the BG may fire again before the next stimulus tone (Fig. 6(a); cycle 4).
As with the order preserving map, in the orderindeterminant case we say that a cycle starts and ends at a BG spike. However, the number of stimulus events per cycle can vary. For example, if there are two consecutive BG spikes, then there would be no stimulus spikes in this cycle, and if there are two consecutive stimulus tones, then both tones are said to be in the same cycle. The solution shown in Fig. 6(c) and (d) is a period-3 solutions, where the first cycle (green) contains two stimulus events, the second (blue) contains 1 stimulus event and the third (pink) contains no stimulus events.
At the beginning of each cycle, we compare the expected period of the BG (the period should no further updates happen) to the amount of time to the next stimulus spike. If the expected period is greater than this amount of time, the BG spike is followed by a stimulus tone, if not, it is followed by another BG spike. This leads us to define a function H S which will be zero if the expected period is less than the time to the next stimulus tone, and 1 otherwise, where Θ(x) is the Heaviside function equal to 0 for x < 0 and to 1 for x ≥ 0 and T (I) is given by (2). If there is a second stimulus tone in a cycle, we do not make another phase update, as we do not have sufficient information about the phase at this time. The BG has not spiked yet, so although we know that it is too late, we do not know by how much. Once the BG spikes again the period rule will act to speed up the BG. As such, there is at most one phase update per cycle.
Using (18) we write the order-indeterminant version of (14) as For the phase variable φ, we note that if there are two or more stimulus spikes in a particular cycle, (15) updates the phase to a negative number. We include a modulo 1 function to convert this to the proportion of time until the next stimulus tone rather than from the previous tone, With the OIEB map now defined, we return to describe the observed dynamics of Fig. 6. During the transient to synchronization, the order of BG spikes and stimulus tones may alternate (Fig. 6(a)). The BG is initially too slow and early, at the first stimulus spike (t = 0ms) phase correction acts to slow down the BG further. At the next stimulus event (t = 500ms), no updates occur as the phase has not been updated. Shortly afterwards, the BG spikes and I is increased to speed up the BG and the phase is updated to φ ≈ 0.8. This marks the end of the first cycle (yellow). The second cycle contains 1 stimulus tone, and as φ ≈ 0.8, phase correction increases I, speeding up the BG. On the 4th cycle (red), the BG is too fast and the BG spikes again before the next stimulus event. Hence, there is no phase correction in this cycle. The order switches again in the 6th cycle (orange), when the BG is too slow and there are two consecutive stimulus tones. As with the first cycle, phase correction occurs at the first stimulus tone in the cycle, but not the second. The iterates of the map can be viewed on the φ − I-plane ( Fig. 6(b)). In general, when an order switch occurs, the iterates cross from the left plane to the right plane for consecutive stimulus tones or from right to left for consecutive BG spikes. Thus the number of transitions between half-planes is equal to the number of order switching events. The iterates I n and φ n correspond to the value of I and φ at the start of the nth cycle, i.e. the value at the nth BG spike after period correction has been applied. The 8 cycles shown in the time course ( Fig.  6(a)) are colored accordingly in the phase plane and the remaining iterates before convergence to the synchronous solution are shown in grey.
The map also exhibits limit cycle behaviour (Fig. 6(c)-(d)). A period-3 solutions is shown, where the first cycle of the periodic orbit (green) contains 2 stimulus tones, the second (blue) contains 1 stimulus tone and the third doesn't contain any stimulus tones (pink). At the first stimulus tone, the BG is too early and slightly too slow (I 1 < I * , φ 1 < 0.5). The phase correction rule decreases I and slows down the BG. As a result, the BG doesn't spike again until after the next stimulus tone, and there are two stimulus tones in this cycle. The cycle ends at the next BG spike, after period correction is applied and I is increased such that I 2 > I 1 . The BG is now too fast, but also late (φ > 0.9). In the second cycle, phase correction increases the value of I and period correction decreases it. We note that the phase correction is stronger as I 3 > I 2 . With a large phase φ ≈ 1 and I 3 > I * , the BG spikes again before the next stimulus tone. Hence, phase correction is not applied and the period rule acts to slow down the BG. The period-3 cycle then repeats.

A. Dynamics of the OIEB map
The values (I * (T stim ), 0) and (I * (T stim ), 1) are also fixed points of the OIEB map. Due to the discontinuity induced by the Heaviside function, it is not possible to linearize about the fixed points. Thus we cannot directly obtain stability information. However, the OIEB map reduces to the order preserving map of Section III whenever a sequence of BG spikes and stimulus tones occur consecutively. Hence, we can estimate the stability of the fixed points of the OIEB map using the order preserving map.
When the fixed point is a stable node, we see monotone convergence to φ = 0 or φ = 1, but when the fixed point is a spiral the iterates jump from converging towards φ = 0 to converging towards φ = 1 when there are two consecutive stimulus spikes and from converging to φ = 1 to converging to φ = 0 when there are consecutive BG spikes. To quantify this we calculate where the discriminant of the eigenvalue given in (17) is 0 for T stim = 500 ms (Fig. 7). The solid lines correspond to the stability boundaries, as in Fig. 5, while the dashed curves represent the real-complex boundaries. For φ = 0 (purple) and φ = 1 (green), the eigenvalues are real to the right of the dashed curves. The curves separate the δ T -δ φ plane into 9 regions, labelled I-IX (Table I).
Region I is the only part of parameter space where the eigenvalues of both fixed points are real and of magnitude less than one. As a result, it is the only region where convergence to a fixed point is consistently order-preserving, as seen in Fig. 4. In region II, initial condition in the basin of attraction of φ = 0 will exhibit order-preserving convergence, but those in the basin of attraction of φ = 1 will not. As φ n increases through 1, it is reset to close to 0 (two consecutive BG spikes) and then converges monotonically to φ = 0 from there ( Fig. 8 (a)). Both fixed points are stable spirals in region III, and iterates initially spiral in towards one of the two fixed points, but when two consecutive BG spikes or stimulus tones occur the iterates cross into the other plane and spiral towards the other fixed point as seen in Fig 6(a)-(b). These crossings continue to occur until the system ultimately converges to the synchronized solution.
Moving into region IV and V, the fixed point at φ = 1 becomes unstable. In region IV, for parameter values suitably close to the stability boundary for φ = 1, the system behaviour is similar to that seen in region III, the iterates switch back and forth between the two fixed points but ultimately converge to the synchronized solution. If δ φ is increased further, we observe periodic orbits ( Fig. 8(b)). This is a period-5 orbit with 2 order switches. The iterates orbit φ = 0 in a counter-clockwise direction (starting with iterate in Q2), until phase correction acts to slow down an already too slow BG, resulting in two consecutive stimulus event (iterate moves from left-half to right-half plane). The iterates switch to orbiting φ = 1, also in a counter-clockwise direction. Also in this region, periodic orbits with unequal numbers of BG spikes and stimulus tones exist (Fig. 8(c)). This solution corresponds to a period-4 orbit with 4 BG spikes, but only 3 stimulus tones. In general, a periodic orbit can have the number of BG spikes differ from the number of stimulus tones by at most one. In region V, the system exhibits behaviour akin to that of region II close to the φ = 1 stability boundary and limit cycle behaviour as we move away from the boundary.
Both fixed points are unstable in region VI. Close to the stability boundaries limit cycle behaviour exists. However, as δ φ is increased, chaotic solutions arise (Fig.  8(d)). The iterates follow similar trajectories, but never return to the same point twice. In this region, further increasing δ φ leads to divergence, where I is decreased below 1 and the BG stops oscillating (Fig. 8(e)).
In regions VII and VIII, the φ = 1 fixed point is stable, but the φ = 0 fixed point is not. For parameter choices close to the φ = 0 stability boundary the system con-verges to φ = 1 with a number of order switches in the transient. As we choose parameters away from the stability boundary, limit cycles emerge. In particular, we see large period limit cycles (Fig. 8(f)), with periods on the order of 100 or larger. Further increasing δ T or δ φ in these regions once again leads to chaotic attractors and then divergence when the minimum of F 1 (I, φ) falls below 1. Both fixed points are unstable in region IX. As in region VI, we observe limit cycles close to the φ = 1 stability boundary. As δ T is increased the limit cycles become chaotic attractors. Indeed, this follows naturally from our results on chaotic dynamics in the one-dimensional map of Section III A. Recall that when the fixed point of that map lost stability, there exists a set of δ T values over which period doubling bifurcations occur leading to chaos. For fixed T stim that set lies on the δ φ = 0 axis in Region IX of Fig. 7. Thus for δ φ sufficiently small it is reasonable to intuit that such period doubling routes to chaos persist. Finally, when δ T is increased further the system diverges as minimum of F 1 (I, φ) falls below 1.
The above analysis demonstrates that the linearization at the fixed points of the order preserving map are useful in predicting and explaining the dynamics of the OIEB map. In particular, the characterization of the stability characteristics shown in Table I allowed us to explain why certain iterates spiralled around the dual origin and either ended up converging as Fig. 6(a) or remained bounded away from the origin as in Fig. 8(b). While we have not exhaustively explored parameter space, a more detailed parsing of parameter space would reveal additional boundaries that separate different types of behaviour within the same regions, e.g. a boundary curve separating the behaviour in Region V from chaotic 6(d) to divergent 6(e).

V. DISCUSSION
In this paper, we have derived and analyzed a twodimensional map that addresses how a beat generator (BG) neuronal network can learn the period and phase of an external periodic signal, such as a metronome. This map represents a type of error-correction algorithm 4-7 in which we define separate learning rules for period and phase that in turn adjust a biophysical parameter, I, of the BG to allow it to synchronize to the stimulus tones. The map originates from our prior work 12 where we introduced an error correction algorithm that relied on the dynamics of neuronal oscillators. Here, we have simplified that presentation to create an event-based map, where iterations of the map correspond to cycles of the BG. Given that the order of BG spikes and stimulus events may alternate, we devised the map independent of the order of events. This led to what we termed an orderindeterminant event-based (OIEB) map. At each event, whether it be a BG spike or stimulus tone, the current values of I and φ are used to determine what the next event is. Period correction is said to occur at BG spikes, while phase correction occurs at stimulus tones. As an iteration of the map is defined as a cycle of the BG, period correction occurs on every iteration, but phase correction may not. Our analysis showed that there exists only a small region of (δ T , δ φ ) (learning rate) parameter space (Region I in Fig. 7) in which the synchronous solution was stable and approachable monotonically. In this region, order preserving convergence signifies that the system is learning a beat by sequentially speeding up (phase advance) or slowing down (phase delay). In other regions (e.g. IV), the BG converges to the synchronous solution but does so in an order-indeterminant way, much in the manner a human participant would when learning a beat. As the strength of the period or phase rule becomes stronger, the BG loses the capability of convergence to the synchronous solution. This loss has two key implications: 1. there is a balance between how much either the phase or period rule can contribute to stable synchrony, and 2. the learning rates cannot be too large, and as such, convergence cannot be too fast. Clearly the ability and rate at which humans learn to keep a beat improves with training. Thus it is plausible that the learning rules proposed here may be subject to some sort of longer term plasticity that is shaped by experience.
In the one-dimensional period correction map as well as the two-dimensional period and phase correction map, we established the existence of both periodic and chaotic solutions. Chaos can arise in strictly continuous systems, such as the logistic map, but also occurs due to discontinuities in the definition of the map, as in switching systems. The discontinuity at φ = 0.5 constitutes a switching manifold. It would be of interest to see if the well-developed bifurcation theory 14 of such systems is applicable to OIEB maps. Chaos also occurs in systems that exhibit border collision bifurcations in which a discontinuity of the map passes through a border on the domain of the map 20,21 . In the order preserving map, the discontinuity of the map is fixed at φ = 0.5 independent of parameters. The OIEB map presents further discontinuities as a result of the Heaviside and modulo functions. These discontinuities appear to occur at interior points on the domain of the map. However, as can be noted from Fig. 8, the orientation of the nullclines changes with parameters. Perhaps such changes in orientation are analogous to border collisions.
An alternate strategy for keeping a beat involves direct periodic forcing to entrain a set of mutually excitatory resonant oscillators 22,23 . These systems typically rely on weak coupling assumptions 24 and can be analyzed by classical methods of phase response curves 25 . An-other recently introduced possibility for learning a beat is to employ a pulse-forced adaptable competition system consisting of units that ramp up at adjustable rates to a fixed threshold to produce specific time intervals 26 . This approach follows the error-correction paradigm in that the ramping rate is subject to a learning rule. Our approach here is distinct from either of these strategies. We propose that the BG is a type of adaptive oscillatory system that is trying to develop an internal representation to an external dynamic source, the metronome. Our approach may be considered as a temporal analogue in the case of timing of repetitive events to Kalman filtering in visual systems that corrects the internal presentation so its representation matches the external world 27 .
The map considered in this paper assumed that the event times were known and could be calculated exactly. However, in the context of humans estimating time, exact time information may not be available. An alternate strategy that we proposed in previous work 12 relies on subdividing time intervals into smaller reference intervals and counting the number of reference intervals between events; a form of discrete, integer-based estimation of timing. From a neuronal network point of view, this could be implemented by counting the number of oscillatory cycles of a fast spiking neuron or neural population between beat generator spikes and stimulus tones. This method of discrete-time estimation could be included in the map with the introduction of an additional variable for each time interval to be estimated. These kinds of systems remain to be studied in future mathematical work.