## Abstract

Extensive research has described two key features of interval timing. The bias property is associated with accuracy and implies that time is overestimated for short intervals and underestimated for long intervals. The scalar property is linked to precision and states that the variability of interval estimates increases as a function of interval duration. The neural mechanisms behind these properties are not well understood. Here we implemented a recurrent neural network that mimics a cortical ensemble and includes cells that show paired-pulse facilitation and slow inhibitory synaptic currents. The network produces interval selective responses and reproduces both bias and scalar properties when a Bayesian decoder reads its activity. Notably, the interval-selectivity, timing accuracy, and precision of the network showed complex changes as a function of the decay time constants of the modeled synaptic properties and the level of background activity of the cells. These findings suggest that physiological values of the time constants for paired-pulse facilitation and GABAb, as well as the internal state of the network, determine the bias and scalar properties of interval timing.

**SIGNIFICANCE STATEMENT** Timing is a fundamental element of complex behavior, including music and language. Temporal processing in a wide variety of contexts shows two primary features: time estimates exhibit a shift toward the mean (the bias property) and are more variable for longer intervals (the scalar property). We implemented a recurrent neural network that includes long-lasting synaptic currents, which cannot only produce interval-selective responses but also follow the bias and scalar properties. Interestingly, only physiological values of the time constants for paired-pulse facilitation and GABAb, as well as intermediate background activity within the network can reproduce the two key features of interval timing.

## Introduction

Two critical independent variables are used to characterize behavior in a wide range of timing tasks: constant error and temporal variability. Constant error is defined as estimated duration minus target duration; hence, it is the parameter that corresponds to the accuracy of the estimate (how close is the subjective interval to the actual interval; Woodrow, 1934). In contrast, temporal variability corresponds to the precision of the subjective estimate, namely, how reproducible are the estimates for a specific duration (Gibbon et al., 1997). In computational terms, bias is the constant error and SD is the temporal variability (Jazayeri and Shadlen, 2010). A myriad of psychophysical studies across timing tasks and species have shown a systematic shift in constant error across durations, with an overestimation and underestimation of intervals for shorter and longer intervals in a global set of durations, and with an indifference interval that corresponds to the intermediate duration associated with constant error equal to zero. We called this the bias property (Woodrow, 1934; Jones and McAuley, 2005; Gu et al., 2011). In addition, temporal variability (i.e., SD) increases linearly as a function of duration. This phenomenon has been called the scalar property of interval timing and is a form of Weber's law (Gibbon et al., 1997; Merchant et al., 2008b; Mendez et al., 2011; García-Garibay et al., 2016). These two parameters are critical for describing temporal processing, yet few modeling studies have addressed the question of how neural networks process temporal information in give rise to bias and scalar properties.

Previous modeling work has shown that interval timing can arise from changes in the network state. One of the critical aspects of these state-dependent models is the fact that the simulated neurons develop interval-selectivity (Laje and Buonomano, 2013). This selectivity is due to the delicate balance between excitation and inhibition inputs defined by paired-pulse synaptic excitatory plasticity and slow synaptic currents (i.e., paired-pulse facilitation and GABAb, respectively; Buonomano and Merzenich, 1998; Buonomano, 2000; Karmarkar and Buonomano, 2007). The notion of interval-selectivity is similar to interval tuning, where cells present an orderly change in activity reaching a peak in the preferred interval. Recently, it has been shown that cells in medial premotor areas and putamen are tuned to the duration of produced intervals (Merchant et al., 2013b, 2015; Bartolo et al., 2014; Crowe et al., 2014). The task performance of monkeys in these neurophysiological studies followed the bias and scalar properties of interval timing (Zarco et al., 2009), suggesting a relation between interval tuning and the two hallmark properties of temporal processing. Now, both paired-pulse facilitation (Varela et al., 1997; Markram et al., 1998) and GABAb (Benardo, 1994; Wang et al., 2010) are time-dependent properties that extend from hundreds of milliseconds to a second and can be modified by experience (Schulz et al., 1994; Misgeld et al., 1995). Yet, little is known about the role of the time constants of these neuronal synaptic properties in interval-selectivity and the two hallmarks of timing. Here, we provide a formal framework on how constant error and temporal variability depend on interval-selectivity, background activity, and synaptic properties of neurons in recurrent neural networks. A Bayesian decoding method was used to optimally read the neural network activity and estimate the constant error and the temporal variability across durations. We found that physiological values of the time constants for paired-pulse facilitation and GABAb can determine the bias and scalar properties.

## Materials and Methods

##### Network model.

We simulated a network of 800 excitatory and 200 inhibitory integrate-and-fire neurons that were sparsely (*p* = 0.05; Buonomano, 2000) and randomly connected (Fig. 1). Each neuron (Amit and Brunel, 1997a; Brunel and van Rossum, 2007) was characterized by a membrane potential *V* that obeyed the equation:
where τ is the time constant of the neural membrane, which was τ = 10 and τ = 5 for excitatory and inhibitory neurons, respectively. When *V* reached the threshold value of *V _{T}* = 20 mV, an action potential was triggered, which was followed by a membrane potential of

*V*= 0 mV during a refractory period of

_{r}*t*= 1 ms. Both

_{r}*I*,

_{Fac}*I*were the input induced currents that provided information about the interval duration and sequence order to each neuron (Fig. 1

_{GABAb}*B*).

*I*was an excitatory current that included paired-pulse facilitation, whereas

_{Fac}*I*was a slow inhibitory current. The internal dynamics of the recurrent network were driven by

_{GABAb}*I*and

_{AMPA}*I*, corresponding to fast excitation and inhibition, respectively (Amit and Brunel, 1997a; Fig. 1

_{GABAa}*C*). Hence, the facilitation and GABAb plasticity occurred at input and not at recurrent synapses. Finally,

*N*(

*t*) is white noise with a SD σ

*and zero mean.*

_{N}##### Input stimuli.

Input *I _{Fac}*,

*I*currents were activated by four pulses separated by a duration

_{GABAb}*d*, generating an isochronous sequence with three serial-order (

_{s}*So*) elements.

*d*covered a range from 100 to 1500 ms, on steps of 50 ms (29 values; Figs. 1

_{s}*A*, 3

*B*). This input mimics the sensory metronome used to drive the motor response of subjects during a metronome synchronization task (Merchant et al., 2011). Because the recurrent network properties were time-varying and depended on stimulus history (i.e., depended on

*d*, the neuronal responses to the first pulse were not included in the analysis. Consequently, the responses to the stimulus across the serial-orders

_{s}*So*1,

*So*2, and

*So*3 corresponded to the activity linked to the second, third, and fourth input stimuli, respectively. Last,

*r*(

_{i}*d*) was the response of neuron

_{s}*i*to the input stimulus

*d*, and had a value of 1 if the neuron generated at least one action potential within 20 ms after each input pulse, or 0 otherwise (Fig. 1

_{s}*D*). Overall, this is a neural network that processed temporal sensory information and produced interval selective responses.

##### Input currents.

Four input currents (*I _{s}*) changed the membrane potential of each cell in the neural network. The temporal dynamics of each current consisted of two coupled linear differential equations (Destexhe et al., 1994):
where τ

*and τ*

_{r,S}*are the rise and decay time constants, and*

_{d,S}*R*is an auxiliary variable. w

_{S}*is the synaptic weight that determines the efficacy of synaptic transmission coming from neuron*

_{i,S}*i*. t

*is the time neuron*

_{k,S}*i*emitted spike number

*k*. We used two input driving currents corresponding to the paired-pulse facilitation and

*GABAb*, and two recurrent network currents:

*AMPA*and

*GABAa*. The time constants were as follows: τ

*= τ*

_{r,AMPA}*= 0.2*

_{r,GABAa}*ms*, τ

*= τ*

_{d,AMPA}*= 0.7*

_{d,GABAa}*ms*, τ

*= 21*

_{r,GABAb}*ms*.

Additionally, the values of *w _{i,AMPA}* and

*w*maintained the recurrent network response stable. These internal recurrent currents (AMPA and GABAa) had a smaller effect on interval processing, because the fast currents decayed before duration

_{i,GABAa}*d*. However, the balance of AMPA and GABAa was essential for maintaining the stability of the neural network (Brunel, 2000; Vogels et al., 2011). Thus, the GABAa weights (

_{s}*w*) were fourfold the AMPA weights (

_{i,GABAa}*w*) because we used a 4:1 ratio between the number of excitatory and inhibitory neurons. In contrast, the driving

_{i,AMPA}*w*and

_{i,Fac}*w*were randomly distributed (Fig. 2

_{i,GABAb}*A*,

*C*) and were responsible for the selectivity to interval durations and serial-order within the network activity.

##### Paired-pulse facilitation.

The input current *I _{Fac}* was an AMPA current where the amplitude increased upon the repeated presentation of stimuli. Paired-pulse facilitation changed over time the synaptic efficacy and reflects the history of presynaptic inputs. Facilitation in actual neurons last hundreds of milliseconds and is caused by calcium influx inside the axon terminal after a presynaptic action potential. Thus,

*I*had a progressive amplification for repetitive presynaptic activity, and it was modeled as an increase in transmitter release probability

_{Fac}*p*(

*t*) (Tsodyks et al., 1998). Between presynaptic action potentials, this probability followed the equation: where τ

*was the time constant and*

_{Fac}*P*

_{0}= 0.17 was the stable release probability. In this model,

*p*(

*t*) changed its value immediately after the appearance of an input action potential following the rule

*p*→

*p*+

*f*(1 −

_{P}*p*) (

*f*= 0.62). Therefore, the synaptic weight

_{P}*w*was multiplied by factor (Varela et al., 1997; Dayan and Abbott, 2001). Notably, τ

_{i,Fac}*and τ*

_{d,GABAb}*were the key independent parameters of our analysis.*

_{Fac}##### Neural encoding.

Principal component analysis (PCA) was used to reduce the high dimensional activity of the recurrent network activity (Jolliffe, 2002; Figure 3*A*). Thus, instead of using 800 excitatory time-varying responses, we used the two principal components as the network responses * r_{PCA}* = (

*r*

_{PCA}_{1},

*r*

_{PCA}_{2}). The two principal components explained >40% of the variability across simulations. The PCA was performed on the covariance matrix (800 × 800) of the excitatory activity of the network using the 2900 observations that corresponded to 100 simulations for each of the 29

*d*. A bivariate normal distribution was fitted to

_{s}*to calculate*

**r**_{PCA}*p*(

*|*

**r**_{PCA}*d*) which is the probability of evoking response

_{s}*given the stimulus*

**r**_{PCA}*d*(Fig. 3

_{s}*B*, ellipses). Figure 3

*C*shows the generalized SD [

*gSD*(

*d*)] for each

_{s}*d*. The

_{s}*gSD*(

*d*) is the square root of the determinant of the covariance matrix and is a measure that is directly associated with the area of the ellipse of

_{s}*p*(

*|*

**r**_{PCA}*d*) of Figure 3

_{s}*B*. The network activity showed an increase in

*gSD*(

*d*) as a function of interval, until

_{s}*d*reached ∼1000 ms (Fig. 3

_{s}*C*). After this duration,

*gSD*(

*d*) was asymptotic. We also measured the percentage of overlap between the

_{s}*p*(

*|*

**r**_{PCA}*d*) ellipses using the adjacent matrix ℳ

_{s}*(*

_{overlap}*d*,

_{s,test}*d*) for every ellipse of

_{s}*d*against

_{s,test}*d*. The diagonal of ℳ

_{s}*(*

_{overlap}*d*,

_{s,test}*d*) corresponds to a perfect overlap (100%), because

_{s}*d*and

_{s,test}*d*are the same (Fig. 3

_{s}*D*). Interestingly, ℳ

*(*

_{overlap}*d*,

_{s,test}*d*) showed an increase in width as a function of

_{s}*d*(Fig. 3

_{s,test}*D*, dotted line). Based on previous neurophysiological studies we defined

*d*as 450, 550, 650, 750, or 850 ms (Merchant et al., 2011, 2013b, 2014).

_{s,test}##### Network decoding.

Here, we assumed a normal prior distribution *p*(*d _{s}*) with a mean of 650 ms and SD σ

*= 172 ms (Fig. 4*

_{prior}*B*). This prior created a clear bias property on the constant error that was difficult to obtain from a uniform prior. Indeed, the mean at 650 ms produced an indifference value close to 650 ms. Bayesian inference uses the posterior probability

*p*(

*d*|

_{s}*) for optimal decoding. This can be estimated using likelihood probability*

**r**_{PCA}*p*(

*|*

**r**_{PCA}*d*) (Fig. 3

_{s}*B*), prior distribution

*p*(

*d*), and Bayes' theorem: In this study, we calculated the Bayes least square decoded estimator (

_{s}*d*), a scalar that was dependent on a particular

_{e}*. In fact,*

**r**_{PCA}*d*was computed as the mean of the posterior probability

_{e}*p*(

*d*|

_{s}*).*

**r**_{PCA}Because *d _{e}* was a deterministic function of

*[*

**r**_{PCA}*d*=

_{e}*f*(

*)], then*

**r**_{PCA}*p*(

*d*|

_{e}*) = δ(*

**r**_{PCA}*d*−

_{e}*f*(

*)), which is a Dirac function. Instead of using the Dirac distribution, which is difficult to calculate numerically, we used a normal distribution with mean*

**r**_{PCA}*f*(

*) and small SD (5 ms), as follows: Finally, we calculated the conditional probability of evoking the duration*

**r**_{PCA}*d*given that tested duration

_{e}*d*was presented (Fig. 4

_{s,test}*A*). This probability depends on population response

*; nevertheless, we removed the dependence of the response by marginalization: Then, we substituted Equation 6 on Equation 7 and obtained the relation: The decoded mean*

**r**_{PCA}*d̂*and SD σ

_{e}*of the decoder density*

_{e}*p*(

*d*|

_{e}*d*) were used as the key dependent parameters in the present study. Indeed, the classical psychometric parameters of constant error and temporal variability correspond to the bias

_{s,test}*CE*(

*d*) =

_{s,test}*d̂*−

_{e}*d*and SD

_{s,test}*TV*(

*d*) = σ

_{s,test}*, respectively (Fig. 4*

_{e}*C*).

##### Mutual information.

We computed the mutual information (MI) between the network decoding and the stimuli. It is a measure of the statistical dependency between the behavioral variable, in this case *d _{s,test}*, and a neural response parameter

*d*(i.e., the output of the decoder). First, we calculated the joint probability

_{e}*p*(

*d*,

_{e}*d*) and

_{s,test}*p*(

*d*) using the decoder probability

_{e}*p*(

*d*|

_{e}*d*). Then, we calculated MI as follows:

_{s,test}## Results

### Network model

We simulated a neural network of 800 excitatory and 200 inhibitory integrate-and-fire neurons that were sparsely and randomly connected (Fig. 1; see Materials and Methods). The neural response of each neuron was driven by two induced input currents: paired-pulse AMPA facilitation and slow GABAb inhibition. Both input currents were time-dependent properties defined by their decay time constants (τ* _{d,GABAb}* and τ

*), which had postsynaptic effects on network neurons lasting hundreds of milliseconds (Fig. 1*

_{Fac}*B*). Within the recurrent network, the cells communicated through AMPA with no paired-pulse plasticity and fast GABAa currents (Fig. 1

*C*). The induced input currents were modulated by four input pulses separated by a duration

*d*(100 to 1500 ms). This produced an input stimulus sequence with three serial-order elements (

_{s}*So*1,

*So*2,

*So*3; Fig. 1

*A*). The neural selectivity to different

*d*was determined by the history of input stimuli in the sequence and balance of weights of the paired-pulse facilitation and GABAb (Fig. 2). Hence, this is an unsupervised recurrent network whose time selective properties depend on the stimulus history and the equilibrium between induced slow excitatory and inhibitory currents. In contrast, the responses to the first stimulus were independent of

_{s}*d*because they were not shaped by slow synaptic currents and paired-pulse synaptic plasticity. Consequently, we focused on cell responses triggered by the last three stimuli in the input sequence. The neural response

_{s}*r*(

_{i}*d*) of each neuron

_{s}*i*was 1 when at least one action potential was fired within the first 20 ms after each stimulus, and 0 otherwise (Fig. 1

*D*).

Neurons in the network showed differences in their average response across the spectrum of *d _{s}*, because the balance between synaptic weights that undergo paired-pulse facilitation and GABAb generates temporal selectivity (Buonomano, 2000). In fact, the response of excitatory neurons to temporal stimuli changed as a function of these synaptic weights. The color gradients shown in Figure 2

*A*correspond to the probability of a neuronal response during

*So*2 as a function of

*w*(abscissa) and

_{i,Fac}*w*(ordinate), using the decay time constants of τ

_{i,GABAb}*= 650 ms and τ*

_{Fac}*= 650 ms. It is evident that the response map is shifted to the right and has a larger slope for a*

_{d,GABAb}*d*of 850 (Fig. 2

_{s}*A*, right, purple) than 450 ms (Fig. 2

*A*, left, magenta). Consequently, the dark gray dots at the bottom of these response maps correspond to a neuron whose combination of weights generate a tuning function with a short-preferred interval, whereas the light gray dots at the top correspond to another neuron with a long-preferred interval (Fig. 2

*B*). The neural network included neurons covering a wide range of

*w*and

_{i,Fac}*w*, and hence, showing an extensive range of preferred durations (Fig. 3

_{i,GABAb}*A*,

*B*). These simulations suggest that artificial neural circuits were capable of encoding temporal information as recorded in the medial premotor cortex (Merchant et al., 2013b). In addition, the interval-selectivity also changed as a function of the time constants of the synaptic currents. At τ

*= 350 ms and τ*

_{Fac}*= 350 ms the response probability weight maps showed a shift to the right and an increase in slope (Fig. 2*

_{d,GABAb}*C*). This shift is accompanied by a change in the preferred interval of the neurons in the network. For example, the tuning functions of the two cells in Figure 2

*B*were shifted to the left when the time constants of the long-lasting currents are smaller (Fig. 2

*D*) and in general, the width of the tuning functions of interval selective neurons had the tendency to increase as a function of preferred interval (Fig. 3

*A*). However, this measure was highly variable across the neurons in the network (data not shown).

The wide range of preferred intervals in the network was a basic feature of the model. Nevertheless, our main goal was to compute an estimate of temporal variability and constant error from a decoder density function, which in turn was calculated from the population activity of the network (Fig. 3). Thus, the synaptic underpinnings of the bias and scalar properties were studied using the decoded estimates of temporal variability and constant error across durations.

Bayesian inference was used to compute the decoder density function. First, we characterized the network encoding of duration *d _{s}* using the response vector

**= (**

*r**r*

_{1},

*r*

_{2}, …,

*r*

_{N}), where

*N*is the number of excitatory neurons. It has been shown that different behavioral parameters can be encoded in the activity of cell populations (Ma et al., 2006; Naselaris et al., 2006) using a low dimensional representation (Stopfer et al., 2003; Churchland et al., 2012). Here, we found that some of neurons in the network were interval selective; and therefore, we used PCA on the covariance matrix of the excitatory activity of the recurrent network

**across all durations (29 durations × 100 simulations; see Materials and Methods) to extract the time-related activity from the simulated population responses. In fact, the average response variability explained by the first two principal components**

*r**= (*

**r**_{PCA}

**r**_{PCA}_{1},

**r**_{PCA}_{2}) was ∼40% (Fig. 3

*B*). The encoding relationship between

*d*and

_{s}*was quantified by the conditional probability*

**r**_{PCA}*p*(

*|*

**r**_{PCA}*d*) that is the likelihood probability of observing the response

_{s}*given that duration was*

**r**_{PCA}*d*. In fact,

_{s}*p*(

*|*

**r**_{PCA}*d*) was computed from the bivariate normal distribution of

_{s}*for each*

**r**_{PCA}*d*(Figs. 3

_{s}*B*, ellipses, 4

*A*). The variability of neural network responses was determined using the generalized SD [

*gSD*(

*d*)], which is the square root of the determinant of the covariance matrix of the

_{s}*and, therefore, is a measure that is directly associated with the*

**r**_{PCA}*p*(

*|*

**r**_{PCA}*d*) and area of the ellipses of Figure 3

_{s}*B*. The network activity showed an increase in

*gSD*(

*d*) as a function of interval, until

_{s}*d*reached ∼1000 ms. After this duration, the

_{s}*gSD*(

*d*) was asymptotic (Fig. 3

_{s}*C*). These results support the notion of different neural timing mechanisms for sub and suprasecond scale (Grondin, 2012, 2014; Méndez et al., 2017), and suggest an important role of paired-pulse facilitation and slow inhibition in the former time range. In addition, we computed the percentage of overlap between the

*p*(

*|*

**r**_{PCA}*d*) ellipses using the adjacent matrix ℳ

_{s}*(*

_{overlap}*d*,

_{s,test}*d*) for all the ellipses of

_{s}*d*against

_{s,test}*d*. The diagonal corresponds to a perfect overlap because

_{s}*d*and

_{s,test}*d*are the same. Notably, the width of the overlap in ℳ

_{s}*(*

_{overlap}*d*,

_{s,test}*d*) increased as a function of

_{s}*d*(Fig. 3

_{s,test}*D*, dotted line), with a larger spread for a

*d*

_{s}_{,test}= 850 ms against surrounding

*d*(Fig. 3

_{s}*D*, green curve) than for the

*d*= 450 ms (Fig. 3

_{s,test}*D*, blue curve).

Then, we performed Bayesian decoding from the population response * r_{PCA}* to estimate the duration

*d*that evoked this response. Classical Bayesian inference depends on posterior probability

_{e}*p*(

*d*|

_{s}*) (Eq. 5), which is the probability that the interval*

**r**_{PCA}*d*was presented given that the population response

_{s}*was recorded. This posterior probability was computed as the product of*

**r**_{PCA}*p*(

*r*|

_{PCA}*d*) (Figs. 3

_{s}*B*, 4

*A*) by

*p*(

*d*), which was the prior probability of (

_{s}*d*), divided by

_{s}*p*(

*), that, in turn, was the probability of the population responses across all*

**r**_{PCA}*d*. The prior probability reflects the knowledge about the occurrence of previous durations, and several studies have successfully used a Gaussian prior distribution to model motor timing across different tasks in humans (Acerbi et al., 2012; Shi et al., 2013). Thus, we used as prior probability a Gaussian distribution with a mean of 650 ms and a SD of σ

_{s}*= 173 ms (Fig. 4*

_{prior}*B*; see the next section for the validation of this σ

*). Given the Bayesian least-squares,*

_{prior}*d*was a scalar that was dependent on a particular

_{e}*and corresponded to the mean of the posterior probability*

**r**_{PCA}*p*(

*d*|

_{s}*) (Fig. 4).*

**r**_{PCA}From the scalar decoded value *d _{e}*, we calculated the decoder density function

*p*(

*d*|

_{e}*d*), which measures the probability of observing an estimated interval

_{s,test}*d*given that the stimulus

_{e}*d*was presented (see Materials and Methods; Eq. 8; Fig. 4

_{s,test}*C*). We computed the mean

*d̂*and SD σ

_{e}*of the decoded density function*

_{e}*p*(

*d*|

_{e}*d*). Importantly, the decoded constant error was calculated as

_{s,test}*CE*(

*d*) =

_{s,test}*d̂*−

_{e}*d*and the temporal variability as

_{s,test}*TV*(

*d*) = σ

_{s,test}*(Fig. 4*

_{e}*C*).

*d _{e}* was a decoded measure of the

**linked to a specific**

*r*_{PCA}*d*; however, some responses

_{s,test}**were evoked by adjacent stimuli**

*r*_{PCA}*d*, corresponding to the overlap of

_{s}*p*(

*|*

**r**_{PCA}*d*) between this specific

_{s}*d*and the tested

_{s,test}*d*. Thus, the width of decoder density [

_{s}*TV*(

*d*)] depended on the ℳ

_{s,test}*(*

_{overlap}*d*,

_{s,test}*d*), because a larger ℳ

_{s}*(*

_{overlap}*d*,

_{s,test}*d*) produced a larger uncertainty that a particular

_{s}*r*was linked to a

_{PCA}*d*. For example, ℳ

_{s,test}*(*

_{overlap}*d*= 850,

_{s,test}*d*= 700) (Fig. 4

_{s}*A*, bottom) was larger than ℳ

*(*

_{overlap}*d*= 450,

_{s,test}*d*= 600) (Fig. 4

_{s}*A*, top). Therefore, the decoder density of

*p*(

*d*|

_{e}*d*= 850) (Fig. 4

_{s,test}*C*, green density) was broader than

*p*(

*d*|

_{e}*d*= 450) (Fig. 4

_{s,test}*C*, blue density).

We systematically measured the effects of the time constants of the slow inhibitory current (τ* _{d,GABAb}*) and paired-pulse facilitation (τ

*) on the bias and scalar properties resulting from the decoded behavior of the constant error and temporal variability across*

_{Fac}*d*, respectively.

_{s}### Hallmark timing properties in the network

Figure 5*A* shows the decoder density function *p*(*d _{e}*|

*d*) for the different

_{s,test}*d*for the serial order

_{s,test}*So*

_{2}. Three key features can be extracted from these data. First, the temporal estimator

*d̂*was tuned to

_{e}*d*because of the balance of weights of the paired-pulsed facilitation and GABAb (Fig. 2), which in turn produced systematic changes in network responses for each stimulus

_{s,test}*d*(Fig. 3

_{s}*B*). Second,

*d̂*was overestimated for the short durations and underestimated for longer intervals (Fig. 5

_{e}*A*, arrows vs lines). Indeed, as seen in Figure 5

*B*, there was a monotonic decrease in constant error (negative

*slope*) as a function of duration, and the indifference interval was observed at 650 ms (with intercept

_{CE}*M*650

*equal to 0 at 650 ms). Therefore, the network decoding followed the bias property observed in psychophysical studies. Last, the temporal variability followed the scalar property of interval timing, namely, there was an increase in variation as a function of*

_{CE}*d*(Fig. 5

_{s,test}*C*), with a positive

*slope*and a positive intercept

_{TV}*M*650

*at 650 ms.*

_{TV}The neural network showed also a response modulation as a function of the serial order of input stimuli. The mean squared error (MSE) of the decoded serial-order reached a minimum when the likelihood *p*(* r_{PCA}*|

*d*) and decoder

_{s}*p*(

*d*|

_{e}*) probabilities corresponded to the same serial order (e.g., the MSE of*

**r**_{PCA}*So*1 was minimum when

*p*(

*|*

**r**_{PCA}*d*) and

_{s}*p*(

*d*|

_{e}*) came from*

**r**_{PCA}*So*1 [right, red] and increased for

*So*2 [right, green], and

*So*3 [right, blue]; Fig. 5

*D*). Nevertheless, the bias and scalar properties were similar across

*So*'

*s*. Hence, in the remaining of the paper we focused our analysis on the data from

*So*2.

We determined the σ* _{prior}* of the Gaussian prior distribution that produced the largest mutual information of duration and that optimally generated the scalar and bias properties. Figure 6 shows that a σ

*= 172 ms (green dotted line) was associated with the largest mutual information (Fig. 6*

_{prior}*D*), and produced a decoded constant error and temporal variability across

*d*that robustly followed the bias and scalar properties (Fig. 6

_{s,test}*B*,

*C*, respectively). Indeed, from the linear fits of Figure 6

*B*on the constant error as a function of

*d*, we extracted the slope (

_{s,test}*slope*) and the intercept at 650 ms (

_{CE}*M*650

*) across σ*

_{CE}*to build Figure 6*

_{priors}*E*, which defines four quadrants. This Figure demonstrates that an increase in σ

*produced both an increase of*

_{prior}*slope*and a drift of

_{CE}*M*650

*from negative to positive values, acquiring a value of zero (which is the indifference interval) close to σ*

_{CE}*= 172 ms (green point). With this σ*

_{prior}*the magnitude of the constant error and temporal variability as a function*

_{prior}*d*is similar to what we found empirically in previous studies (Zarco et al., 2009; Donnet et al., 2014). Note that a small σ

_{s,test}*did not produce a perfect alignment of constant error ∼650 ms, suggesting that not only normal prior distribution but also the neural network population activity determined the bias alignment as a function of*

_{prior}*d*. Similarly, Figure 6

_{s}*F*shows the four quadrants resulting from plotting the

*slope*against the intercept

_{TV}*M*650

*at 650 ms. In this case, the increase in σ*

_{TV}*produced an increase in both the*

_{prior}*slope*of the scalar property and the intercept

_{TV}*M*650

*. It is important to notice that the σ*

_{TV}*in the Bayesian decoding modified the slope of the scalar property (Fig. 6*

_{prior}*C*).

Overall the scalar property depended on both the likelihood probability that was associated with the representation of duration in the activity of the neural network and the prior probability of the decoding. However, two features of the likelihood probability had a substantial impact on the decoded temporal variability *TV*(*d _{s,test}*): the generalized SD

*gSD*(

*d*) (Fig. 3

_{s,test}*C*, Pearson correlation:

*TV*vs

*gSD*,

*r*= 0.23,

*p*< 0.0001) and, more importantly, the percentage of overlap between the

*p*(

*r*|

_{PCA}*d*) ellipses, corresponding to ℳ

_{s}*(*

_{overlap}*d*,

_{s,test}*d*) [Figs. 3

_{s}*D*, 4

*A*,

*C*; Pearson correlation:

*TV*vs ∫

*ℳ*

_{ds}*(*

_{overlap}*d*,

_{s,test}*d*),

_{s}*r*= 0.95,

*p*< 0.0001].

### Synaptic effects

The properties of the decoder density function *p*(*d _{e}*|

*d*) depend on the time-varying input synaptic currents. Hence, we changed the temporal dynamics of

_{s,test}*I*and

_{Fac}*I*so that their decay time constants (τ

_{GABAb}*and τ*

_{d,GABAb}*) varied from 50 to 1000 ms (Benardo, 1994; Markram et al., 1998). We investigated how the decay time constants affected time accuracy and precision for different durations. For this purpose, we characterized how the slopes and intercepts of constant error and temporal variability, with respect to*

_{Fac}*d*, varied as a function of network parameters. Figure 7

_{s,test}*A*shows

*M*650

*and*

_{CE}*slope*values across time constants. The bias property corresponds to negative

_{CE}*slope*(Fig. 7

_{CE}*A*, inset, light green and yellow regions). There is a wide range of values of τ

*and τ*

_{Fac}*where the bias property was satisfied because the mean of the decoding prior distributions is 650 ms. However, there is a sector with low τ*

_{d,GABAb}*where the*

_{Fac}*slope*was positive and that did not follow the bias property (inset, red and blue regions). On the other hand, Figure 7

_{CE}*B*shows the

*M*650

*and*

_{TV}*slope*values across time constants. The blue and red regions located mainly along the diagonal of Figure 7

_{TV}*B*are associated with positive

*slope*; hence, they are the sectors where the scalar property was followed (Fig. 7

_{TV}*B*, inset, light red and blue regions). The gray shadow regions in Figure 7,

*A*and

*B*, show the combinations of decay constants associated with no linear relation between the constant error or the temporal variability as a function of

*d*. Last, we used mutual information as a global measure of information about duration

_{s,test}*d*contained in the decoder. Because we used five different values for

_{s,test}*d*, mutual information close to

_{s,test}*log*

_{2}(5) indicated that

*d*was strongly dependent on stimulus duration, whereas mutual information close to

_{e}*log*

_{2}(1) implied random decoding with no interval-selectivity. Figure 5

*C*shows mutual information as a function of τ

*and τ*

_{d,GABAb}*. It increased monotonically as a function of the time constants, although it reached asymptotic values at large and similar τ*

_{Fac}*and τ*

_{d,GABAb}*. Notably, the largest values of mutual information were observed at low τ*

_{Fac}*and intermediate to high τ*

_{d,GABAb}*, as well as at intermediate to high τ*

_{Fac}*and low τ*

_{d,GABAb}*.*

_{Fac}Together, these findings support the notion that interval-selectivity, as well as the bias and scalar properties showed complex relations when generated by a recurrent neural network with long-lasting input currents. Only specific ranges of the decay time constants produced networks with proper duration tuning that also followed the two hallmark properties of interval timing (Fig. 7*A*,*B*, insets, bright areas; Fig. 7*D*, dark area). For example, the area around τ* _{Fac}* = 650 ms and τ

*= 650 ms (Fig. 7, black circle in each panel) met these criteria and corresponds to physiological values of these decay constants.*

_{d,GABAb}There were combinations of τ* _{Fac}* and τ

*where the bias and/or scalar properties were not met. For example, when τ*

_{d,GABAb}*= 100 ms and τ*

_{Fac}*= 300 ms the decoding followed the scalar but not the bias property (Figs. 7, black triangles, 8*

_{d,GABAb}*C*,

*D*, black triangles), whereas with a τ

*= 600 ms and τ*

_{Fac}*= 150 ms the decoding results followed the bias but not the scalar property (Fig. 7, black squares, 8*

_{d,GABAb}*C*,

*D*, black squares). Finally, there is a region around τ

*= 950 ms and τ*

_{Fac}*= 50 ms where both bias and scalar property were not satisfied (Fig. 7, black diamonds, 8*

_{d,GABAb}*C*,

*D*, black diamonds). Two important notions emerge from these findings. First, the prior distribution centered ∼650 ms induced the bias property in specific ranges of the decay time constants (Fig. 7

*C*). However, there were regions where the constant error did not follow the bias property, suggesting that not only prior distribution but also the activity of the neural network (Fig. 8

*A*) was able to shape the constant error. Second, the increase in width of ℳ

*(*

_{overlap}*d*,

_{s,test}*d*) is necessary to observe an increase in temporal variability as a function of

_{s}*d*(Fig. 8

_{s}*B*).

### Noise effects

We also determined the effect of noise within the recurrent network on the decoding of duration. The purpose was to determine how the background activity of the network, in combination with the long-lasting synaptic input, modulates the constant error and temporal variability across intervals. We changed the membrane potential *V* of each neuron with white noise of magnitude σ* _{N}* (measured as percentage of threshold

*V*; see Materials and Methods). As expected, larger noise produced a decrease in accuracy, precision, and mutual information (Fig. 9). However, the effects of noise on these parameters had an interaction with the decay time constants of the input currents. Network noise produced an increase in the negative

_{T}*slope*, indicating that the decoding accuracy decreased with noise, although the bias property was held across noise levels (Fig. 9

_{CE}*A*). Figure 9

*B*shows that the temporal variability increased as a function of noise, and that the slope of the scalar property (

*slope*) showed an increase as a function of noise. Notably, at large noise levels the temporal variability remained similar across

_{TV}*d*, and therefore, the scalar property was lost. Finally, mutual information decreased asymptotically as a function of noise (Fig. 9

_{s}*C*).

## Discussion

The present study supports three conclusions. First, recurrent neural networks with time-dependent properties show time-selective responses that can follow the bias and scalar properties of interval timing in the range of hundreds of milliseconds when read by an optimal Bayesian decoder. Second, interval-selectivity, accuracy, and precision showed complex behaviors as a function of the decay time constants of the modeled synaptic properties, suggesting that the constant error and the temporal variability, present in numerous species and timing tasks, depend on specific combinations of time constants for paired-pulse facilitation and GABAb. Finally, the increase in the variability of cell membrane potential (i.e., noise) of the recurrent network produced a generalized decrease in mutual information, accuracy, and precision of temporal processing. At large noise levels the scalar property is lost.

Several neurocomputational models have been implemented to describe the main operations behind temporal processing in the subsecond scale. These models include ramping activity (Durstewitz, 2003; Miller et al., 2003; Méndez et al., 2014), drift diffusion (Simen et al., 2011; Merchant and Averbeck, 2017), synfire chains (Bienenstock, 1995; Hass et al., 2008), neural oscillators (Treisman, 1963; Medina et al., 2000; Matell and Meck, 2004; Merchant and Bartolo, 2018), coincidence detectors (Karmarkar and Buonomano, 2002; Matell and Meck, 2004), and state-dependent networks (Buonomano, 2000; Karmarkar and Buonomano, 2007; Hardy and Buonomano, 2018). Most of them implicitly or explicitly produce time selective responses and can replicate the scalar property of interval timing. However, none of the previous models were designed to reproduce, as a whole, the neural selectivity to durations, as well as the bias and the scalar properties (Hass and Durstewitz, 2014; Addyman et al., 2016). We accomplished this by using recurrent neural networks that simulated time-varying synaptic properties and membrane potential noise variations and then using Bayesian inference to optimally decode the activity of the network.

The scalar property is a hallmark of temporal processing and had an enormous impact on the study of the psychophysical and neural basis of interval timing (Gibbon et al., 1997; Grondin, 2001; Merchant et al., 2008b,c; Merchant and Averbeck, 2017). In addition, bias property, also known as the Vierordt's law, has been described since the end of the 19th century (Vierordt, 1868; Woodrow, 1934) and has had a key role in the study of the neurobiology of temporal processing. Thus, the preferred internal interval (Fraisse, 1963), the temporal memory (McAuley and Jones, 2003; Gu and Meck, 2011; Taatgen and van Rijn, 2011), the distribution of the intervals used in an experiment (i.e., global temporal context; Jones and McAuley, 2005; Jazayeri and Shadlen, 2010), and the timing abilities of the subjects (e,g, musicians vs non-musicians; Cicchini et al., 2012) have an important influence on the indifference interval and the slope of the bias property. Both the scalar and bias properties are commonly seen in the temporal performance across tasks and species, and both can be explained using the classical scalar timing model (Treisman, 1963; Gibbon et al., 1997; Grondin, 2001; Shi et al., 2013). From a Bayesian perspective, these properties emerge from the integration of two independent functions: likelihood and prior distribution. A recent study has shown a Bayesian encoding system in V1, where the activity evoked by complex visual stimuli corresponds to the likelihood response, and the spontaneous activity that changes with the animal visual experience represents the prior distribution (Berkes et al., 2011). Using this probabilistic framework, the temporal performance can be explained properly (Jazayeri and Shadlen, 2010; Shi et al., 2013; Gu et al., 2016), although it has been shown that the bias feature is not optimal when a high accuracy is needed, such as in the case of professional drummers reproducing a beat (Cicchini et al., 2012). A previous study successfully used the Bayesian scheme to model time reproduction in humans, although it assumed sensory and motor likelihood functions that were Gaussian and that explicitly followed the scalar property (Jazayeri and Shadlen, 2010; Shi et al., 2013; Roach et al., 2017). In contrast, the present study obtained the likelihood probabilities from the synaptic properties of cells in recurrent neural networks without further assumptions. Indeed, the neural population activity in the network follows the scalar property when the proper long-lasting synaptic constants are used (Figs. 4, 7). However, it is important to mention that we had to include a normal prior distribution to reproduce the change in timing accuracy as a function of duration, consistent with previous studies on time processing (Acerbi et al., 2012; Shi et al., 2013; Roach et al., 2017). A uniform prior did not generate a robust bias property. Furthermore, the width of the Gaussian prior (σ* _{prior}*) influenced the slope of the scalar feature. Hence, an optimal reader that has access to this sensory timing information needs to include prior knowledge about the task conditions to reproduce the bias property and modulate the slope of the scalar feature (Ivry and Hazeltine, 1995; Merchant et al., 2008c, 2013a).

Interval tuning occurs early in the processing of sensory stimuli, as in the case of the midbrain for auditory stimuli in the range of tens of milliseconds (10–100 ms; Sayegh et al., 2011). In contrast, the primary auditory cortex (A1) shows cells that are tuned to a wider range of durations (10–500 ms) than the inferior colliculus, suggesting that duration selectivity in A1 results from integration along the time domain across the auditory ascending pathway (He et al., 1997). Conversely, the first node of the visual pathway that shows duration tuned cells is the primary visual cortex or V1, with cells tuned to a range of 50 to 400 ms (Duysens et al., 1996). In these neurophysiological studies, the width of the tuning function does not seem to scale as a function of the interval (A1: He et al., 1997; Fig. 5; V1: Duysens et al., 1996; Figure 4), although no specific analysis were performed to determine the scalar property of neurons. Accordingly, the present study predicts that using an Integrate and Fire recurrent neural network it is difficult to see the scalar property from the tuning curves of individual neurons.

However, at the population level the generalized SD *gSD*(*d _{s,test}*) (Fig. 3

*C*) and, more importantly, the percentage of overlap between the ellipses [ℳ

*(*

_{overlap}*d*,

_{s,test}*d*)], had a large correlation with the decoded scalar property. On the other hand, it has been shown recently that single striatal cells follow the scalar property with a spread in the discharge rate that increases as a function of the onset response latency (Mello et al., 2015). The medial premotor areas also show temporal scaling in their response profile, namely, the responses are stretched or compressed in accordance with the produced interval (Crowe et al., 2014; Wang et al., 2018). Altogether these studies support the idea that temporal processing depends on the top-down reading of sensory timing information and the generation of predictive signals in the motor system that obey the scalar property (Merchant and de Lafuente, 2014; Merchant and Yarrow, 2016; Petter and Merchant, 2016).

_{s}With the model reported here, the two key timing properties showed a complex behavior and could be reproduced simultaneously for a limited set of neural noise magnitudes and decay time constants, which correspond to physiological values. The magnitude of the dynamic synaptic facilitation (Varela et al., 1997; Markram et al., 1998), as well as the slow inhibition due to the activation of GABAb receptors (Benardo, 1994; Tamás et al., 2003; Wang et al., 2010) in cortical and subcortical areas show time constants that affect the processing of intervals in the scale of hundreds of milliseconds. In addition, the dynamics of both time-dependent properties can be modulated by behavioral experience (Schulz et al., 1994; Misgeld et al., 1995). As far as we know, this is the first modeling study describing how the input synaptic properties of a recurrent network can modulate the accuracy and precision of time processing. Thus, this paper shows that the ability of cortical networks to intrinsically process temporal information and follow the scalar and bias properties depends on reaching decay time constants ∼650 ms for paired-pulse facilitation and 650 ms for GABAb metabotropic receptors, both of which are physiologically meaningful values.

Large neural noise in the recurrent network, which can be linked to the spontaneous activity of cells in the mammalian brain (Shadlen and Newsome, 1994; Amit and Brunel, 1997b; Vogels et al., 2011), produced an asymptotic decrease in mutual information and an increase of constant error and temporal variability. Crucially, at large noise levels the scalar property was lost because the temporal variability was similar across *d _{s}*. These findings are consistent with the observation that a subgroup of Parkinson's disease (PD) patients showed an increase in variability and

*slope*during a time production task, whereas others lost the scalar property of interval timing (Merchant et al., 2008a; Gu et al., 2016). Hence, the heterogeneity of temporal performance in PD patients can be associated with different magnitudes of damage in the dopaminergic input, which in turn can generate different neuronal noise levels in the corticobasal ganglia circuit, as demonstrated previously (Levy et al., 2000; Timmermann et al., 2003). To end, it is important to notice that we limited the effect of fast recurrent synapses (AMPA and GABAa) on interval processing on our network. Thus, the reported results were associated to synaptic noise and not to the network recurrence. In further studies we will determine the role of the recurrent connections within the neural network on temporal processing, because it is well known that asynchronous and irregular dynamics within a neural network, due to a complex balance between fast recurrent synapses, is important to induce broad tuning activity in response to the transient input stimuli (Vogels et al., 2011; Hennequin et al., 2017).

_{TV}In conclusion, learning a timing task may be associated with the appearance of interval tuning, the reaching of particular levels of spontaneous activity, as well as the adoption of specific time constants on paired-pulse synaptic plasticity and slow synaptic currents in cortical or subcortical recurrent circuits. Although the experiments needed to test these hypotheses are technically challenging, our laboratory has already demonstrated interval tuning in cells of the primate medial premotor areas and the putamen during interval production tasks (Merchant et al., 2013b, 2015; Perez et al., 2013; Crowe et al., 2014).

## Footnotes

This work was supported by Consejo Nacional de Ciencia y Tecnología Grants 236836 and 196, and Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica Grant IN202317, Secretaria de Ciencia, Tecnología e Innovación to H.M. We thank Victor de Lafuente, Warren Meck, Fernando Peña, and Dobromir Dotov for their fruitful comments on earlier versions of the paper, Luis Prado and Raul Paulín for their technical assistance. Oswaldo Pérez is a doctoral student from Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional Autónoma de México (UNAM) and received fellowship 204516 from CONACYT.

The authors declare no competing financial interests.

- Correspondence should be addressed to either Dr. Hugo Merchant or Oswaldo Pérez, Institute of Neurobiology, UNAM, Campus Juriquilla, Boulevard Juriquilla 3001, Querétaro, Qro 76230, México. hugomerchant{at}unam.mx or doswaldo.perezm{at}gmail.com