A neural circuit model for human sensorimotor timing

Humans and animals can effortlessly coordinate their movements with external stimuli. This capacity indicates that sensory inputs can rapidly and flexibly reconfigure the ongoing dynamics in the neural circuits that control movements. Here, we develop a circuit-level model that coordinates movement times with expected and unexpected temporal events. The model consists of two interacting modules, a motor planning module that controls movement times and a sensory anticipation module that anticipates external events. Both modules harbor a reservoir of latent dynamics, and their interaction forms a control system whose output is adjusted adaptively to minimize timing errors. We show that the model’s output matches human behavior in a range of tasks including time interval production, periodic production, synchronization/continuation, and Bayesian time interval reproduction. These results demonstrate how recurrent interactions in a simple and modular neural circuit could create the dynamics needed to control timing behavior.

S ensorimotor coordination in humans is remarkably flexible. We can anticipate events based on few observations and use that information to adjust our movements. For example, musicians can use a metronome to adjust the tempo of their movements, and children can rapidly coordinate their movements during a clapping game. However, we still lack an understanding of how networks of neurons generate such coordinated movements.
Recent studies have proposed that the neural basis of sensorimotor coordination may be understood using the language of dynamical systems [1][2][3] . The key intuition is that recurrent neural networks in the motor cortex form dynamical systems 4,5 whose output can be controlled by sensory inputs 2,[6][7][8][9] . This idea has been explored in large-scale distributed recurrent neural network models 10 . However, the complexity of these network models often makes it difficult to understand their behavior from an algorithmic perspective.
Timing provides a prime example of sensorimotor coordination and is crucial in behaviors demanding the generation of delayed motor responses, generation of rhythmic movements with a desired tempo, and synchronization of movements to anticipated events. An early model proposed that the brain controls action timing by adjusting the speed of an internal clock whose ticks are integrated toward a fixed threshold 11 . Consistent with this proposal, experiments in animal models have found that neural activity in anticipation of a delayed response reaches a fixed threshold [12][13][14] at a rate that is inversely proportional to the delay period 7,15,16 . These results suggest that the brain supports flexible timing by controlling the speed at which neural activity approaches a movement initiation threshold.
Recently, it was shown that flexible control of speed can be achieved through nonlinear interactions within a simple model consisting of a pair of units with reciprocal inhibitory connections 7 . In this model, the speed at which the output evolves toward a movement initiation threshold can be adjusted flexibly via a shared input (Fig. 1a). From a dynamical systems perspective, this model can be viewed as an open-loop controller that converts an instruction (i.e., shared input) to the desired dynamics (i.e., speed).
However, the larger utility of this model depends on whether it can be extended to accommodate additional constraints associated with temporal control of movements. In particular, temporal control of movements must accommodate noise in the nervous system 17 , prior expectations 18 , and sensorimotor delays imposed by internal 19 and external 20 temporal contingencies. These factors limit the capacity of open-loop systems to achieve robust temporal control. Efficient control systems often rely on sensory feedback to combat noise, and when sensory feedback is delayed, they rely on a mechanism to predict future sensory and motor states 19,21,22 . Here, we augment the original open-loop system with a sensory anticipation module (SAM) and a feedback mechanism, and show that the resulting neural circuit model can accommodates internal noise, prior expectations, and sensorimotor delays. Finally, we demonstrate that the model can capture key features of human behavior in a number of classic timing tasks (Fig. 1b-e).

Results
We will describe the full model in four steps. We start by introducing a basic circuit module (BCM) that acts as a flexible open-loop controller for producing desired time intervals. We extend the BCM to a motor planning module (MPM) capable of producing isochronous rhythms. We then introduce a SAM that provides the means for anticipating and predicting upcoming temporal events. Finally, we introduce the full model that combines the MPM and SAM to create a system that can dynamically coordinate motor plans and actions with anticipated and unanticipated stimuli in a range of behavioral tasks.
BCM for interval production. The BCM has been described in detail previously 7 . Briefly, the BCM includes three units, u, v, and y, each representing the average activity of a population of neurons (Fig. 1a, top). u and v inhibit one another and receive input, I, that is tonic, or constant over time. y receives excitatory input from u and inhibitory input from v, and leverages the nonlinear dynamics of mutual inhibition between u and v to generate ramp-like activity. Finally, the model initiates a "movement" when y crosses a fixed threshold, y 0 . The rate dynamics of u, v, and y are as follows: W uI and W vI denote the strength with which I drives u and v, respectively. W uv and W vu denote the strength of inhibitory coupling from v to u, and from u to v, respectively. τ is the time constant of each unit. θ(x) is a sigmoidal function that maps the input to an output between 0 and 1 (see "Methods"). Finally, η u , η v , and η y are stochastic synaptic inputs to each unit and are modeled as independent white noise with standard deviation σ n . We assume W uI = W vI = 6 (identical shared excitatory input), W uv = W vu = 6 (symmetric mutual inhibition), W yu = W yv = 1, and τ = 100 ms for all units, a value consistent with previous models. With these parameters, the model functions in a dynamical regime with three fixed points: an unstable point at u = v, and two stable fixed points with either u dominating or v dominating (Fig. 1a, middle). In this regime, u and v evolve, or change over time, toward one of the stable fixed points depending on the initial conditions. These dynamics lead to a ramp-like activity in y whose rate is inversely related to I (Fig. 1a, bottom). In other words, the BCM acts as an open-loop controller that converts an instruction conveyed by I to a dynamic output, y, that evolves toward the threshold, y 0 , at an appropriate speed (Fig. 1a). Therefore, by adjusting I, the BCM can flexibly adjust the movement initiation time.
MPM for periodic production. We next considered the production of multiple movements with a fixed inter-productioninterval (IPI). This is challenging with the BCM which can only produce a single action. One solution is to use multiple concatenated BCMs, but that seems unrealistic as the number of modules would grow with the number of movements. We tackled this problem using a modification of the BCM, which we call the MPM. The MPM is identical to BCM except for an additional reset mechanism after each movement that allows the model to generate an arbitrary number of timed outputs (Fig. 2a). From a neurobiological perspective, the reset signal can be generated by the corollary discharge associated with the movement command 23,24 , or by sensory feedback triggered by the movement. We implemented the reset as a transient 10 ms pulse, I p , activated right after the threshold crossing (Fig. 2a, dashed line). Mathematically, this can be formalized as follows:

ARTICLE
In the MPM, u p , v p , and y p evolve identically to the BCM between consecutive threshold crossings (Fig. 2b). In addition to producing an output at each threshold crossing, however, the model activates I p , which resets u p and v p , and causes a rapid drop in y p . This allows the circuit to restart the dynamics and generate another output. As expected, IPIs increase monotonically with I (Fig. 2c, d; r 2 = 0.84; F(1, 160) = 828.6; p ≪ 0.01) within a suitable range of inputs ( Supplementary Fig. 1).
In humans, IPIs are variable with a standard deviation that increases linearly with the mean 25 . To evaluate IPI variability in the model, we performed simulations in the presence of Gaussian noise (see "Methods"). When noise levels were not too high ( Supplementary Fig. 2), the model exhibited a qualitatively similar behavior to that of humans ( Fig. 2e): IPI variability increased monotonically with IPI (one-tailed F test; I = 0.75 to I = 0.76: p = 0.014, F(80, 70) = 0.60; I = 0.76 to I = 0.77: p < 0.01, F(57, 70) = 2.95; I = 0.77 to I = 0.78: p < 0.01, F(40, 57) = 10.93). However, this relationship was nonlinear in the model, which is unsurprising given the simplicity of our assumptions regarding the nature of noise in the brain and the underlying dynamics.
SAM for predicting future events. In many circumstances, the IPI is not known in advance and has to be adjusted based on external timing cues. Here, we considered a case in which the IPI has to be adjusted by the interval between external temporal events. In general, measuring time between events can be achieved in two ways. One way is to implement a circuit whose output is a ramp with a fixed slope. In such a circuit, the output level provides a moment-by-moment estimate of elapsed time ( Supplementary Fig. 3a). Alternatively, the circuit may function predictively, and adjust the slope of the ramp such that the output reaches a certain expected level at the anticipated time of the next event ( Supplementary Fig. 3b).
A recent study in monkeys found evidence in support of the predictive mechanism 20 : neural dynamics in the frontal cortex were adjusted so that responses reached an expected state at the expected time of the stimulus. Based on this finding, we developed a SAM that measures time predictively (Fig. 3a). The SAM behaves identically to the MPM except that it does not generate an action when its output, y s , reaches y 0 . Instead, the SAM adjusts the input so that y s would reach y 0 exactly at the expected time of the sensory feedback. When I is too high, y s would increase at a slower pace than needed, and its value at the time of the sensory feedback would be below y 0 . Conversely, when I is too low, y s would go beyond y 0 at the time of the feedback. Accordingly, the system must decrease I when y s < y 0 , and increase I when y s > y 0 . The SAM implements these , and drive an output unit, y, with excitatory and inhibitory connections, respectively. Excitatory and inhibitory connections are shown by triangles and circles, respectively. The speed control mechanism can be understood by analyzing system dynamics in the phase plane of u and v (middle; after Wang et al. 7 ). Dots represent u and v over time for a large (pink) or small (blue) input for no noise and with initial conditions set so that that u eventually dominates v. The system's saddle and stable fixed points are indicated by larger circles and squares, respectively. Dashed and solid lines indicate the nullclines of u and v, respectively. Inputs configure the positions of the nullclines for u and v, and therefore, control the speed. Inset shows the relationship between speed and input. Because u excites and v inhibits y, operation of the system in this regime increases y (bottom) until it reaches a threshold for action initiation (dashed line). b-e Classic timing tasks used to study human timing behavior. b Periodic production requires the subject to produce a series of actions over time (vertical lines), with a constant inter-production-interval (IPI). Top and bottom are examples of two different IPIs. c Synchronization requires the subject to time a series of actions (vertical black lines) such that they are simultaneous with a series of sensory inputs (e.g. flashes; vertical magenta lines) with a set inter-stimulus-interval (ISI). d Interval reproduction requires the subject to measure an interval, t s , demarcated by two stimuli (flashes; vertical magenta lines) and to produce an interval, t p , by initiating an action (dashed red line). t p has to matche t s as accurately as possible. e Synchronization/ continuation requires the subject to synchronize actions to a series of inputs with an ISI selected at random from a prior distribution and then, after the stimulus is extinguished, continue to produce actions with an IPI matching the ISI selected on a given trial.
adjustments by updating the value of I dynamically at a rate proportional to the error signal, y s −y 0 : In this update rule, K scales the rate of change of I with the error signal, and s is a gating parameter that is set to 1 when the feedback is on and 0 otherwise. This formulation allows the model to update I only when the feedback is on (s = 1). Fig. 3b illustrates the response of the SAM to three equidistant sensory inputs with a short (400 ms) or long (1000 ms) inter-stimulusinterval (ISI). The first stimulus (S1) triggers a transient input, I s . This input does not alter I since the system has not yet measured the ISI, but resets u s and v s ( Fig. 3a; see "Methods"). Following S1, the output, y s , increases over time. For a short ISI, y s at the time of the second stimulus (S2) is lower than y 0 causing I to decrease (Fig. 3b, left). In contrast, for a long ISI, y s is greater than y 0 causing I to increase (Fig. 3b, right). The transient input, I s , triggered by S2 also resets u s and v s , allowing the dynamics to evolve with the updated speed (Fig. 3b). By iterating this process after each sensory input, the SAM dynamically adjusts its output such that y s will eventually match y 0 at the precise time of each stimulus input. These results are robust with respect to when the first stimulus is presented (750 ms in Fig. 3b  Using simulations of the SAM, we analyzed the effect of K and I 0 on the root-mean-squared-error (RMSE) between the predicted time (the time at which y s crosses y 0 ) and the actual time of the stimuli ( Fig. 3c; see "Methods" for simulation details). When I 0 is large and K is small, the SAM generates large errors. Intermediate values of I 0 and K lead to smaller errors, and a specific combination of (K Ã ; I Ã 0 ) minimizes the RMSE (see Fig. 3c). Exploration of the parameter space revealed three key results. First, I Ã 0 is largely independent of noise level and the number of stimuli (Fig. 3d, left). This is consistent with I 0 serving as an initial estimate before any feedback is integrated. Second, K * decreases with increasing σ n (Fig. 3d, right). This is expected because, when internal noise increases, the error signals generated by the SAM are less reliable and should be appropriately discounted. Third, as N increases, K * decreases, which reduces the weight given to each error and allows the model to integrate across inputs (Fig. 3d, right). Finally, we verified that the SAM updating process is robust with respect to the units' initial conditions (u 0 and v 0 ) (Supplementary Fig. 4a-d), and can be generalized to different dynamical regimes ( Supplementary Fig. 5). Therefore, the SAM provides a plausible mechanism for measuring time intervals predictively.
Sensorimotor updating by combining the SAM and MPM. So far, we found that the MPM can produce different IPIs depending on the level of input, and the SAM can adjust the input level based on the anticipated time of sensory events. Accordingly, we reasoned that the SAM and MPM might together be able to generate timed outputs that are in register with incoming sensory events: the SAM would adjust the input based on sensory events and the MPM would use that input to adjust IPI. We therefore connected the SAM to the MPM by having them share the input, I (Fig. 4a), and measured IPI as we changed the ISI randomly in a blocked fashion (Fig. 4b). The circuit was able to successfully track ISI throughout the run ( Fig. 4c; see "Methods"). After each block transition, the SAM detects the error between IPI and ISI and adjusts I so that the MPM can gradually bring IPI closer to the new ISI (Fig. 4c). This mechanism allows IPIs to increase lawfully with the ISI (r 2 = 0.53; p ≪ 0.01; Fig. 4d).
Full circuit model for synchronization. Coupling between the MPM and SAM allowed the circuit to adaptively match its output frequency (IPI) to the input frequency (ISI). However, it failed to match the input and output in terms of phase, as evidenced by the uniform distribution of phase differences between inputs and outputs ( Fig. 4e; Rayleigh test of uniformity, p = 0.10). The source of this problem is that the circuit has no mechanism to determine whether the MPM output leads or lags the stimulus. Previous models have proposed highly nonlinear systems for phase adjustment [26][27][28][29] . We found that our model provides a simple solution to this problem. Since the output of the SAM (y s ) and MPM (y p ) are synchronized to sensory and motor events, respectively, the difference between them provides a graded signal reflecting their relative phase. Specifically, when y s > y p , the MPM is lagging and should be sped up, and when y s < y p , the MPM is leading and should be slowed down. To implement this solution, we augment the input I with a signal ΔI of the form where α controls the learning rate. This adjustment can be realized by a unit which receives excitatory and inhibitory input form y p and y s , respectively (Fig. 5a, cyan). Fig. 5b demonstrates this correction scheme. Initially, y p slightly lags behind y s . This asynchrony generates a biphasic ΔI whose value is transiently positive (between the stimulus onset and motor output) and then negative until the next stimulus onset. The addition of this error signal reduces the total tonic drive to the MPM and increases the output speed relative to the SAM. After this adjustment is applied over several ISIs/IPIs, the MPM becomes increasingly synchronized with the stimulus. This strategy allows the circuit to gradually reduce asynchrony and narrow the distribution of phase errors (Rayleigh test of uniformity, p ≪ 0.01; Fig. 5c).
We examined the behavior of the model for different values of α and K. Larger values of α make the model more sensitive to phase errors and allow the model to cancel asynchronies more effectively (Fig. 5d, left panel). However, larger values of α cause rapid beat-by-beat adjustments in IPI, which can compromise rhythmicity (Fig. 5d, middle panel). As such, the optimal value of α depends on the relative cost of minimizing phase asynchrony and IPI error. In our simulations, values between 0.1 and 0.15 allow the model to minimize the sum squared error of phase and IPI (Fig. 5d, right panel). The parameter K, on the other hand, has little effect on how α influences phase asynchrony and IPI error (Fig. 5d, different colors).
This mechanism allowed the model to capture two counterintuitive observations in human behavior. First, as shown in Fig. 5c, the MPM had a persistent phase lead relative to the stimulus (−27.14°± 71.45°) in a manner similar to humans performing analogous tasks 30 . Second, the interaction of α with the level of noise caused the model to occasionally skip a beat ( Supplementary Fig. 6), which occurs occasionally in humans performing similar tasks 31 .
Model responses to perturbations are similar to humans. We compared the behavior of the model to that of humans in  Fig. 2 and 3. b Example inter-stimulus-interval (ISI) sequence for a trial of the ISI tracking task. Colors indicate different ISIs. We initiated each trial with a block of twenty consecutive 800 ms ISIs. Each subsequent block of 20 ISIs was randomly selected from a discrete uniform distribution between 600 and 900 ms. Each trial consisted of 5 blocks and 100 total ISIs. c Inter-production-intervals (IPIs; black circles) associated with each ISI (colors) for two example trials of the ISI tracking task. d Density of IPIs (grayscale) as a function of the associated ISI. Dashed line indicates perfect tracking. e Distribution of the phase of motor output. A phase of 0 indicates output that is synchronous with the stimulus input.
synchronization tasks during which the rhythmic input was perturbed in one of three different ways: a step change in ISI (Fig. 6a, top), a phase shift (Fig. 6b, top), and a jitter in the timing of a single event (Fig. 6c, top). When facing a sudden change in ISI (Fig. 6a, top), response times have to be adjusted in both phase and frequency. For example, after an uncued increase in ISI, the first motor response would lead the stimulus. To regain synchrony, one has to delay the response to cancel the phase lead and additionally increase the IPI to match the new ISI. Human subjects concurrently reduce asynchrony and adjust IPI such that after a few samples, actions are in sync with sensory inputs 30,[32][33][34] . A hallmark of this errorcorrecting strategy is that IPIs exhibit a transient overshoot relative to the new ISI (Fig. 6a, middle) [35][36][37] .
To test the behavior of the model in response to a change in ISI, we used the following protocol: we allowed the model to reach steady-state for an ISI of 800 ms, stepped the ISI to 1000 ms, and measured subsequent IPIs produced by the model. We repeated this procedure 1000 times and measured how the average IPI changed over time after the step change. Qualitatively, the model exhibits the transient IPI overshoot relative to ISI that is observed in humans, and gradually adjusts the IPI to match the new ISI (Fig. 6a, bottom). Quantitatively, the degree of overshoot depends on the model's two learning rates, K for frequency (Fig. 6a, bottom), and α for phase ( Supplementary Fig. 7).
Next, we analyzed the behavior of the model in response to a phase shift, in which a single ISI is increased or decreased, causing all subsequent stimuli to occur at a different phase without any change to the subsequent ISIs (Fig. 6b, top). In this case, humans adjust their response times so that the initial mismatch is gradually reduced (Fig. 6b, middle) 38 . This model was able to capture this response pattern, and the speed of recovery could be adjusted by α (Fig. 6b, bottom).
Finally, we considered a perturbation in which a single stimulus is jittered temporally. This perturbation alters two consecutive ISIs in equal and opposite directions without To synchronize the MPM with SAM output, a second circuit pathway was introduced to measure the difference between these two outputs. The difference, weighted by α, augments the input to the MPM with input ΔI. b MPM output, y p , SAM output, y s , and their difference weighted by α = 0.1 (i.e. ΔI = α(y p − y s )). Augmenting I with ΔI adjusts controller output to the MPM such that the time of production tends to match the time of flashes (vertical dashed lines). c Distribution of the phase of production for the circuit with (black) and without (red) augmented input. d Optimization of α. Increasing α decreases asynchrony (left; as defined in the Methods). In contrast, increasing α increases the sum of the squared IPI errors (middle). As a result, a limited range of α values minimize Pythagorean errors, defined as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðIPI À ISIÞ 2 þ ðAsynchronyÞ 2 q (right). changing either the phase or the ISI of the subsequent stimuli (Fig. 6c, top). In response to the perturbation, human subjects exhibit a characteristic change in IPI relative to ISI. For example, when a single stimulus is delayed, subjects detect the error and delay their next response accordingly. However, since the perturbation is transient, subjects have to then undo their corrective response, which is done gradually over the course of the subsequent stimuli 39,40 (Fig. 6c, middle). Again, the model was able to capture this response pattern, and the dynamics of the error correction was moderately dependent on α (Fig. 6c, bottom,  circles). These results demonstrate that the circuit model can capture key qualitative aspects of human behavior in response to various perturbations of stimulus timing. Behavioral studies have found individual differences in how subjects respond to these perturbations [35][36][37][38][39][40] . Accordingly, we verified that our model can capture this behavioral diversity through adjustments of model parameters (K and α) with an accuracy that is comparable to previously proposed algorithmic models 30,34,36 ( Supplementary  Fig. 8).
Circuit model implements Bayesian interval reproduction. In the absence of any prior knowledge about the ISI of the first few beats, there is no way for the circuit (or a human) to choose an informed initial speed, and only sensory feedback can guide motor timing. However, when an observer has some prior expectation about the possible values of ISIs, the optimal strategy, as prescribed by Bayesian integration, is to integrate sensory inputs with that prior knowledge. This integration enables the system to reduce variability due to internal noise 25,41 and to estimate the ISI with greater precision.
This Bayesian integration strategy is a hallmark of human behavior in time interval reproduction experiments 18,42 . In these experiments, subjects are typically provided with a sample interval, t s , drawn from a fixed prior distribution, and are asked to produce a matching interval, t p . The behavior of an optimal Bayesian observer performing this task exhibits two characteristic features. First, t p values are biased toward the mean of the t s distribution. Second, the magnitude of biases becomes smaller when measurements of t s are more reliable.
Recently, we tested humans in a time interval reproduction experiment in which t s was sampled from a fixed discrete uniform distribution ranging between 600 and 1000 ms 31 . The experiment involved two trial types. In one type (Fig. 7a, left), which we refer to as "1-2-Go," subjects were presented with the first two beats of an isochronous rhythm (1)(2) and were asked to produce the third omitted beat (Go). In the second type (Fig. 7a, right), referred to as "1-2-3-Go," subjects were presented with the first three beats (1-2-3) and had to synchronize their response with the fourth omitted beat. Results provided clear evidence that subjects relied on their prior knowledge, since t p values were biased toward the mean of the t s distribution (Fig. 7b, purple). Moreover, the presentation of two ISIs in the 1-2-3-Go compared to one ISI in the 1-2-Go condition allowed subjects to measure t s more accurately and reduce the bias (Supplementary Fig. 9). We simulated the model to test if it can emulate these characteristics. The SAM received input pulses representing the beats of an isochronous rhythm (2 pulses for the 1-2-Go task and 3 for the 1-2-3-Go task) with the ISI sampled from the same distribution that was used in the human experiment. We defined the production interval t p of the model as the interval between the final input pulse and the next time the SAM's output exceeded y 0 .
The behavior of the model in this task depends on three key parameters: the initial input, I 0 , the update constant, K, and the level of noise in the model, σ n . I 0 determines the speed of dynamics before the first measurement of t s . This reflects the model's initial guess about t s based on prior expectations. K determines the strength with which the SAM updates the input based on the difference between anticipated and observed t s . Finally, σ n determines the variability of the model's behavior due to internal noise. Here, we assume that the statistics of the environment are stationary (i.e., the prior distribution is fixed), an assumption consistent with the experimental design 31 . Under these conditions, a single set of parameters is required to optimize the performance of the model. Note that we do not explicitly model the learning process that optimizes model parameters for a specific prior distribution.
Together, the parameters I 0 , K, and σ n allow the model to capture a range of behaviors observed in human subjects (see Supplementary Table 1). Results in Fig. 7b show the behavior of the model that was fit to the data from two human subjects in the 1-2-Go and 1-2-3-Go tasks (for parameter values, see Supplementary Table 1). Evidently, the model captured several key features of Bayesian integration present in the human behavior (Fig. 7b, magenta): (1) Average t p increased monotonically with t s (Fig. 7b, black circles), (2) noise in the model caused t p to vary on a trial-by-trial basis for the same t s (Fig. 7b, black error bars), (3) average t p was biased toward the mean of the prior distribution (Fig. 7b, deviation from the dashed unity line), and (4) the magnitude of bias was smaller in the 1-2-3-Go compared to the 1-2-Go trials (z = 2.7, p < 0.01, one-sided Wilcoxon signed rank test; Supplementary Fig. 9). An alternative implementation of this task using the full circuit architecture yields similar production times ( Supplementary Fig. 10).
To further compare the model's behavior to that of the human subjects, we partitioned the root-mean-square error (RMSE) between t p and t s to two terms, a BIAS term that measures the overall bias, and a VAR term that measures average variability across all values of t s (see "Methods"). The model was able to capture the observed BIAS and VAR across subjects in both 1-2-Go and 1-2-3-Go (Fig. 7c), indicating that it can correctly implement the bias-variance trade-off exhibited by humans during interval reproduction. In this model, the update constant, K, plays a central role in determining the bias-variance trade-off (compare Fig. 7b, top and bottom). For a Bayesian observer, the magnitude of bias is determined by noise in the measurements of t s and by the imposed cost function. Accordingly, the value of K in the model was adjusted to achieve a level of bias-variance trade-off that is inversely related with the inherent noise in the model (σ n ) and the operative cost function (r = − 0.94, p = 0.0002; see Supplementary Table 1). These results provide a potential mechanistic account of linear updating algorithms proposed previously to describe Bayesian timing behaviors 43 ( Supplementary Fig. 11).
Bayesian synchronization/continuation. As our final test, we examined the model in a Bayesian synchronization/continuation task that demands sensory anticipation, motor timing, and Bayesian integration. In this task, subjects are asked to tap synchronously with a metronome (synchronization phase) and continue to tap with the same tempo without the metronome until instructed to stop (continuation phase). In the classic version of this task, subjects have no uncertainty about the ISI of the metronome as it is kept constant across trials. Recently, humans were tested on a variant of this task in which ISI for each trial was   sampled randomly from a prior distribution 44 . While performing this task, humans exhibit distinct patterns of biases in their IPIs. During synchronization, responses were biased toward the mean of the ISI distribution and the magnitude of the bias decreased with the number of synchronization stimuli (Fig. 8a, purple lines). These observations are consistent with the Bayesian integration scheme that we discussed previously in the context of a time interval reproduction task (Fig. 7). During the continuation phase, the bias persisted and its magnitude increased substantially (Fig. 8a, red lines) for all subjects (Fig. 8b). This increase in bias was due to a drop in the sensitivity of IPI to ISI (shallower slope relating IPI to ISI) as well as an overall bias toward shorter IPIs (Fig. 8c, left). We fitted the model parameters (I 0 , K, and α) to individual subject's behavior in this task (see "Methods"; Supplementary Table 2). Model fits generated IPIs that were biased towards the mean ISI. Consistent with subjects' behavior, the magnitude of the bias decreased during synchronization and increased during continuation (Fig. 8a, black line). Moreover, like humans, the increase in bias during continuation was due to a combination of decreased sensitivity to the ISI (decrease in slope of IPI-ISI relationship, z = 2.1, p = 0.02, one-sided Wilcoxon signed rank test) and an overall shift towards shorter IPIs (decrease in IPI of the middle interval, z = 2.1, p = 0.02, one-sided Wilcoxon signed rank test; Fig. 8c, right; Supplementary Fig. 12).
In the model, this pattern of biases is explained by the augmented input, ΔI, to the MPM, which is necessary for the full model to synchronize its outputs with sensory inputs. Following the final stimulus input in the synchronization phase, the SAM no longer resets and instead proceeds toward the terminal fixed point. As a result, its output, y s , becomes fixed at a value greater than y 0 . However, y p continues to oscillate between values smaller than y 0 because it is reset by the motor output. This mismatch makes ΔI negative and decreases the overall input to the MPM. The decreased input, in turn, speeds up the dynamics of the MPM and shortens IPIs during the continuation phase across all ISIs. Further, the impact of this negative ΔI is integrated over time, leading to an asymmetric impact on the speed of MPM dynamics associated with a short ISI compared to a long ISI. Specifically, longer ISIs lead to a greater increase in speed than shorter ISIs, accounting for the decreased sensitivity to the ISI during continuation. It is important to note that these results require that the SAM does not reset. Therefore, an alternate circuit configuration in which the output of the MPM resets the SAM will not be able to capture the increase in bias during continuation. This finding validates our assumption that resetting the SAM relies on sensory inputs and not motor outputs. Further, this feature of the circuit model allowed it to more accurately capture the observed pattern of bias than simple linear updating algorithms that have previously been used to model human synchronization behavior (Supplementary Fig. 13).

Discussion
To coordinate movements with external stimuli, the dynamics of neural activity must be adjusted based on sensory inputs 7,15,45 . However, precise control of dynamics is challenging because neural signals are subject to internal noise, and timing cues are often discrete and delayed 21 . Here, we proposed a simple neural circuit model that can address these challenges in a wide range of timing tasks.
The model is comprised of two modules, an MPM that generates timed actions, and a SAM that predicts the time of upcoming external events. Both modules were engineered to emulate a recently discovered speed control mechanism in the frontal cortex of monkeys during a flexible motor timing task 7 . The MPM uses this speed control mechanism in conjunction with a reset to produce rhythmic outputs. The SAM uses the same speed control mechanism in conjunction with an error-correction scheme that enables it to predict the timing of external events.
Coupling the MPM and SAM creates a closed-loop control system with versatile timing capacities in the presence of noise and delayed sensory feedback. For example, the model was able to capture human behavior in classic beat synchronization tasks 30,[32][33][34]37 . In this case, the function of the SAM was to detect and correct discrepancies between the model's output and external beats so that the MPM could adjust the outputs accordingly. The model was also able to emulate the biases associated with Bayesian time interval reproduction in the presence of noisy measurements when sample intervals were drawn from a fixed prior distribution 18,31 . To do so, the initial input to the model had to be adjusted based on the prior distribution, and the model's updating parameters had to take the magnitude of the noise into account. Finally, a phase-correction mechanism between the SAM and MPM enabled the model to capture several non-trivial features of human behavior during a Bayesian synchronization/continuation task 44 that previous algorithmic models could not. Previous studies have proposed algorithmic models for how humans coordinate their motor outputs with sensory inputs. The basic algorithm is to compare the time of each motor output to the time of the corresponding sensory input, and use the error between the two to adjust the timing of the next movement [31][32][33][34]43 . This algorithm is conceptually similar to how our model functions, and its error-correcting behavior in beat synchronization tasks is comparable to our model 30,43 (Supplementary Fig. 8). However, by virtue of implementing this algorithm using a neural circuit, our model's predictions can be readily compared to neural data. Moreover, our model provides a mechanistic understanding of hypothesized predictive processes that algorithmic models have attributed to human behavior during perceptual timing tasks 43,46,47 .
Previous circuit models of timing largely fall into three classes, each with their own strengths and weaknesses. The first class is based on the accumulation of ticks of a central clock 11,[48][49][50] . Similar to our circuit model, these clock-accumulator models rely on ramping activity that is observed in individual neurons 12,13 . However, these models do not explain how the recurrent circuit interactions lead to such ramping activity. The second class uses large recurrent neural circuits capable of producing rich dynamics 2,7,51-53 . These models can produce activity patterns that are strikingly similar to those observed in local populations of neurons, but it is unclear how they can flexibly integrate sensory and motor feedback. The third class uses a system of coupled oscillatory units [26][27][28]54,55 . These models can produce sophisticated timing behaviors and can integrate sensory information across a range of time scales 26 . However, the activity profile of neurons in the brain regions causally involved in timing is typically not oscillatory 7 . Our model provides an understanding of the link between these model classes. First, it explains the ramping activity in terms of recurrent dynamics due to interactions between neurons. Second, the model's dynamics can be flexibly adjusted based on incoming sensory input and motor feedback. Third, the model generates oscillations at its output. However, unlike abstract models comprised of inherently oscillatory units, our model generates oscillations through circuit-level interactions that are consistent with the speed control mechanism observed in single neurons in the frontal cortex 7 .
Finally, the wiring of our circuit model might provide insight into the individual functional contributions of the cortical and subcortical systems important to action and perceptual timing 56 . One node, the basal ganglia, has been heavily implicated in action timing through lesion studies 57 and physiological evidence 7,45,58 . An important principle of the basal ganglia function is the inhibition of downstream neural activity 59 . Given their proposed function as a competitive selection mechanism 60 , these inhibitory pathways may be the substrate for implementing the mutual inhibitory interactions needed for the temporal control of movements 6 . Another key node, the cerebellum, has also been linked to timing through lesion studies 61 , physiology 62,63 , and modeling 44 . A hallmark of cerebellar function is the detection and correction of sensory errors during sensorimotor control 64,65 , which is the key function of the SAM in our model. Finally, the output of both the basal ganglia and cerebellum are sent transthalamically to regions of the frontal cortex involved in sensory and motor timing 66 . Accordingly, we speculate that the integration of the sensory anticipation and motor production modules in our model may rely on interaction between the basal ganglia, cerebellum, and frontal cortical areas 66 .

Methods
Modeling and analyses were performed in Python 3.5.4 and Matlab R2017a.
BCM for interval production. The fundamental circuit architecture consists of three rate units, u, v, and y, which are governed by the following set of equations: where τ, the time constant of each unit, was set to 100 ms and θ(x), the activation function of each unit, was specified by θðxÞ ¼ 1=½1 þ expðÀxÞ, and I is a tonic input. W uI , W uv , W vu , and W vI specify the weighting of the interactions between units and each was set to 6 for all simulations. Similarly, we set W yu = W yv = 1. Noise in the system was modeled by 0 mean Gaussian white noise inputs, represented by the variables η u , η v , and η y for unit u, v, and y, respectively. These variables were independently sampled at every time point from a Gaussian distribution with mean set to 0 and standard deviation specified by σ n . The units were initialized at u = 0.7, v = 0.2, y = 0.5, and I = I 0 , where I 0 is a free parameter. All simulations were carried out using Euler's method with a step size of 10 ms.
Motor planning module. The MPM consists of four rate units I, u p , v p , and y p , which are simulated based on the following equations: I p specifies a transient input to u p and v p that serves to reset the system when y p > y 0 . We set I p to 50 and y 0 to 0.7 for all simulations. The activity level of I controls the speed of the dynamics of the MPM. The units were initialized at u p = 0.7, v p = 0.2, y p = 0.5. Noise, represented by the variables η u p , η v p , and η y p , was injected into each unit as in the BCM. All other parameters are the same as the BCM.
Sensory anticipation module. The SAM consists of four rate units I, u s , v s , and y s . The system evolves according to the equations where K is a free parameter and η indicates noise input to each unit as in the BCM. The module receives a binary input sI s , where s represents the visual stimulus (s = 1 when the stimulus is on and s = 0 when the stimulus is off) and I s is set to 50. This input serves to reset the values of u s and v s at the time of each stimulus. All other parameters are the same as the BCM. Sequences of stimuli were implemented by setting s = 1 periodically. Each stimulus presentation lasts 10 ms. The units u s , v s , and y s were initialized as in the MPM and I was initialized at I 0 , a free parameter. The dynamics were allowed to run forward for 750 ms before the first stimulus was presented. This initialization period serves to put the circuit closer to steady-state so that the effect of reset signals is consistent across different ISIs. I does not update during the first presentation of the stimulus (by setting K = 0 during the first stimulus presentation).
Optimization of SAM parameters was achieved by simulating the module with 100 different pairs of I 0 and K for each σ n . K was uniformly sampled from values between 1 and 8.0, and I 0 was sampled from values between 0.77 and 0.79. N stimuli were presented to the model, where N = 2, 3,…,10. ISIs were sampled from a discrete uniform distribution with five values between 600 and 1000 ms.
For each set (σ n , I 0 , K, N), the model was run as described above and the RMSE was calculated according to where N s = 500, the total number of iterations across all ISIs. The pair (I 0 , K) that results in the smallest RMSE was recorded. The procedure was repeated 10 times for each σ n to obtain a distribution of the optimal (I 0 , K). In Fig. 3c, I 0 and K were picked from a 30 × 30 grid spanning the the same range as the optimization procedure, and N s = 5000.
Behavioral tasks. We used a suite of behavioral tasks to test circuit model timing performance. In all tasks, we performed numerical simulations of the dynamical equations expressed above using Euler's method and a time step of 10 ms. All behavioral experiments were performed with the approval of the Committee on the Use of Humans as Experimental Subjects at MIT after receiving informed consent.
Periodic production task. To achieve periodic production, we simulated the MPM in isolation at four different levels of input, I, uniformly spaced from 0.75 to 0.78, for 40 s at each level. To measure the performance of the model, we calculated the IPI as where t n was the time of the nth action produced by the model. For all levels of input, we set σ n = 0.01. To compare the timing behavior of the model across levels of input, we calculated the mean and standard deviation of the first 40 IPIs generated by the circuit.
ISI tracking task. Circuit model synchronization was tested in a randomized ISI tracking task. The ISI was defined as where m n was the time of the nth stimulus. On each trial of the task, we initialized the ISI at 800 ms and simulated a sensory input with 20 consecutive 800 ms intervals. After the first block of 20 ISIs, the ISI for the next 20 intervals was selected at random from a discrete uniform distribution between 600 and 900 ms with four possible values. This process was repeated four times, to generate five blocks of ISIs, corresponding to a sequence of 100 total intervals.
Circuit parameters were selected as follows: I 0 was set to 0.771; K was selected to be 2; α was chosen to be 0 for simulations with no augmented input (e.g. Fig. 4) and 0.1 for simulations with augmented input (e.g. Fig. 5). I 0 was set such that the MPM generated an 800 ms IPI on average, matching the initial ISI of every trial. This ensured that differences in circuit output resulted from updates to the value of I, rather than the value of I 0 selected.
The nth IPI was quantified as in Eq. (19) and compared to the nth ISI. The asynchrony between an action and the stimulus presented at time m n was defined as a n ¼ t n À m n ; ð21Þ where t n was the time of the action which was closest in time to m n . The relative phase of the nth action, ϕ n , was calculated as where the phase is in units of degrees. For all simulations, we set σ n = 0.01.
ISI perturbation task. The basis of the mechanisms behind human synchronization behavior are generally probed using perturbations of the ISI after an experimental subject has reached steady-state synchronization performance to a given ISI (see Repp 30 for a review). To compare the circuit model synchronization performance to that of humans, we explored the behavior of the model in response to three common ISI perturbations: (1) a step change in the ISI, where the circuit synchronized to a stimulus with an ISI of 800 ms before the ISI was stepped to 1000 ms for all subsequent ISI. (2) A phase shift in the stimulus timing, where the circuit synchronized to a stimulus with a 500 ms ISI before the timing of stimuli was shifted by increasing a single ISI to 600 ms. All subsequent ISIs were 500 ms. (3) Stimulus jitter, where the circuit synchronized to an ISI of 500 ms before the timing of a single stimulus is perturbed, while subsequent stimuli remained in phase with the stimulus before the perturbation. To accomplish this, we perturbed two successive ISIs; the first was increased to 600 ms and second was decreased to 400 ms. We simulated 1000 trials of each perturbation type and calculated the asynchrony associated with each action as in Eq. (21) and the mean IPI between actions as in Eq. (19). To ensure the circuit model was fully synchronized before a perturbation, we simulated 30 ISIs of the same duration before applying the perturbation. The circuit model was simulated with I 0 = 0.771 and σ n = 0.005. The values of K and α were varied to explore the behavior of the model with different sensitivities to errors in simulation and synchronization.
1-2-Go and 1-2-3-Go tasks. We compared the circuit model interval reproduction behavior to that of humans performing a timing task we refer to as 1-2-Go and 1-2-3-Go. The methods used for testing human interval reproduction are summarized briefly here. Please see our behavioral paper for a full description of the task and methods 31 .
The behavior of nine human subjects was analyzed. Subjects sat in a dark, quiet room at a distance of approximately 50 cm from a display monitor. The display monitor had a refresh rate of 60 Hz, a resolution of 1920 by 1200, and was controlled by a custom software (MWorks; http://mworks-project.org/) on an Apple Macintosh platform.
The interval reproduction task consisted of two randomly interleaved trial types "1-2-Go" and "1-2-3-Go". On 1-2-Go trials, two flashes demarcated a sample interval (t s ). On 1-2-3-Go trials, three flashes demarcated t s twice. For both trial types, subjects had to reproduce t s immediately after the last flash by pressing a button on a standard Apple keyboard. On all trials, subjects had to initiate their response proactively and without any additional cue. Subjects received graded feedback on their accuracy.
Each trial began with the presentation of a 0.5°circular fixation point at the center of a monitor display. The color of the fixation cued the trial type. After a 2 second delay, a warning stimulus and a trial cue were presented. The warning stimulus was a white circle that subtended 1.5°and was presented 10°to the left of the fixation point. The trial cue consisted of 2 or 3 small rectangles 0.6°above the fixation point (subtending 0.2°× 0.4°, 0.5°apart) for the 1-2-Go and 1-2-3-Go trials, respectively. After a random delay with a uniform hazard (250 ms minimum plus an interval drawn from an exponential distribution with mean of 500 ms), flashes demarcating t s were presented. Each flash lasted 6 frames (~100 ms) and was presented as an annulus around the fixation point with an inside and outside diameter of 2.5°and 3°, respectively. t s was sampled from a discrete uniform distribution ranging between 600 and 1000 ms with a 5 samples. To help subjects track the progression of events throughout the trial, after each flash, one rectangle from the trial cue would disappear (starting from the rightmost). The produced interval (t p ) was measured as the interval between the time of the last flash and the time when the subject pressed the response key. We discarded any trial when the subject responded before the second (for 1-2-Go) or third (for 1-2-3-Go) flash and any trial where the response was 1000 ms after the veridical t s . We further used a model based approach to identify "lapse trials," or trials in which t p was not related to t s 31 . We then calculated the mean t p conditioned on t s and trial type.
The SAM was run for 100 simulated trials for each t s (600 ms, 700 ms, 800 ms 900 ms, 1000 ms). To simulate the 1-2-Go task, the SAM was presented with 2 stimuli that were separated by the selected t s . The 1-2-3-Go task was simulated similarly but with 3 stimuli. The production interval t p of the model was defined as the interval between the final input pulse and the next time y s exceeded y 0 .
For each subject, we calculated the mean and standard deviation of t p for each t s in both the 1-2-Go and 1-2-3-Go tasks (10 task conditions in total). We call these values μ 1,subject , STD 1,subject , . . . , μ 10,subject , STD 10,subject (one pair of values for each condition).
For each set of parameters (σ n , I 0 , K), the model was run as described above. We calculated the mean and standard deviation of the model's t p for each t s , and for both the 1-2-Go and 1-2-3-Go tasks. We call these values μ 1,model , STD 1, model , . . . , μ 10,model , STD 10,model .
The fitting was done by alternating between fitting σ n and jointly fitting (I 0 , K). The parameters were uniformly and independently sampled from the following intervals: [1.0, 8.0] for K, [0.77, 0.79] for I 0 and [0.005, 0.4] for σ n . For each step, 100 sets of parameters were sampled and the optimal set was used to update the parameters.
To evaluate the anticipated t s associated with I 0 and σ n , we initialized the model with I 0 , presented a single stimulus to the model 750 ms after initialization, and determined the interval between this stimulus and the next time y s exceeded y 0 .
We defined BIAS and VAR for each subject and model simulation as where N = 5 is the number of target intervals, t p i and σ 2 i are the mean and variance of production intervals that correspond to the target interval t s i .
Synchronization/continuation task. We compared the circuit model to the behavior of humans performing a synchronization/continuation task in which the inter-stimulus-interval (ISI) was selected at random each trial from a set distribution. The methods used for testing human experiments are summarized briefly here. Please see Narain et al. 44 for a full description of the task and methods 44 .
We analyzed the behavior of six human participants. Stimuli were viewed from a distance of approximately 67 cm in a dark, quiet room. The display monitor had a refresh rate of 60 Hz, a resolution of 1920 by 1200, and was controlled by a custom software (MWorks; http://mworks-project.org/) on an Apple Macintosh platform.
Each trial began with the presentation of a red circular fixation stimulus (diameter = 0.75°visual angle) in the center of the screen. After a variable delay (200 ms plus an additional amount of time which was exponentially distributed with mean = 300 ms and a maximum value of 2300 ms), a synchronization stimulus was flashed four times with an ISI chosen from a discrete uniform distribution (five intervals, minimum = 550 ms, maximum = 817 ms). The flashing stimulus was a gray annulus centered around fixation with an inner diameter of 1°and outer diameter of 1.25°. Participants were instructed to tap a response key synchronously with the third and fourth synchronization flashes and continue tapping to reflect the ISI until instructed to stop. The number of continuation taps required was three plus an exponentially distributed number with a mean of nine and a maximum of 22. The first IPI was defined as the interval between the middle of the second flash and the first key press. Subsequent IPIs were defined as the interval between successive key presses.
Biases in IPIs were calculated according to the following procedure. For each interval, ISI i , we found the mean IPI for the kth interval in the sequence ðIPI i;k Þ. The mean squared bias for each k was defined as where N = 5 is the number of distinct ISI values. We simulated 21 trials of the full synchronization model for each ISI. In each trial, the stimuli presented to the model were a series of three flashes separated by the ISI of that trial. After the three flashes, the stimulus was set to s = 0 and the model was run and allowed to produce 17 more productions. The IPI of the circuit was defined as the time difference between successive productions.
We fixed σ n at 0.01 and varied three free parameters K, I 0 , and α. For each subject, we simulated 100 different parameter combinations with K randomly selected from the interval [0.01, 5], I 0 from the interval [0.76, 0.78] and α from the interval [0.01, 0.1]. For each combination of parameters, the model was run as described above and bias was calculated as in Eq. (27). We then found the combination of parameters that minimized the mean squared errors between the BIAS k observed in the subject and the BIAS k of the circuit, across k = 1, …, 20.
Synchronization bias was quantified as the bias for k = 3, and continuation bias was the bias for k = 7.
Linear behavioral algorithm for sensorimotor timing. For comparison purposes, we also developed a linear updating model for synchronization based on previous work 30 . In the absence of any error, the model computes the (n + 1)th production time, t n+1 , by adding the expected ISI, T n , to the nth production time, t n . When there is a discrepancy, the updating rule incorporates two additional error terms. The first error term measures the asynchrony, which is the difference between t n and the time of the nth metronomic input, m n . The second error term measures the difference between T n and the observed ISI, ISI n . In the model, these two error terms are weighted by β Asynch and β ISI , respectively. t n þ 1 ¼ t n þ T n À β Asynch ðt n À m n Þ À β ISI ðT n À ISI n Þ: ð28Þ We simulated this timing model for different values of β Asynch and β ISI to compare this standard algorithm to the behavior of the circuit model.
In simulations of the 1-2-Go and 1-2-3-Go tasks, we added a noise term, η, to t n+1 . η was sampled independently at every time point from a Gaussian distribution with mean 0 and standard deviation σ n . For fitting, β Asynch was set to 0 and parameters were uniformly and independently sampled from the following intervals: [0.0, 1.0] for β ISI , [700, 900] for T 0 , and [30,120] for σ n . The fitting procedure is the same as described for the circuit model.
In simulations of the synchronization/continuation task, the production times t n during the synchronization phase were simulated like the 1-2-3-Go task. For the continuation phase, t n+1 = t n + T N + η, where η was sampled independently at every time point from a Gaussian distribution with mean 0 and standard deviation σ n , and T N is the value of T at the final stimulus presentation. For fitting, σ n was set to 50, and parameters were uniformly and independently sampled from the following intervals: [600, 1200] for T 0 , [0, 2] for β ISI , [0, 1] for β Asynch . The fitting procedure is the same as described for the circuit model.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
No new data were collected. Data for previously published work is available at https:// jazlab.org/resources/.