## Abstract

Accurate timing is critical for a wide range of cognitive processes and behaviors. In addition, complex environments frequently necessitate the simultaneous timing of multiple intervals, and behavioral performance in concurrent timing can constrain formal models of timing behavior and provide important insights into the corresponding neural mechanisms. However, the accuracy of such concurrent timing has not been rigorously examined. We developed a novel behavioral paradigm in which rhesus monkeys were incentivized to time two independent intervals. The onset asynchrony of two overlapping intervals varied randomly, thereby discouraging the animals from adopting any habitual responses. We found that only the first response for each interval was strongly indicative of the internal timing of that interval, consistent with previous findings and a two-stage model. In addition, the temporal precision of the first response was comparable in the single-interval and concurrent-interval conditions, although the first saccade to the second interval tended to occur sooner than in the single-interval condition. Finally, behavioral responses during concurrent timing could be well accounted for by a race between two independent stochastic processes resembling those in the single-interval condition. The fact that monkeys can simultaneously monitor and respond to multiple temporal intervals indicates that the neural mechanisms for interval timing must be sufficiently flexible for concurrent timing.

- interval timing
- concurrent timing
- fixed interval

## NEW & NOTEWORTHY

How accurately animals can estimate the duration of elapsed time for multiple intervals simultaneously has not been investigated rigorously. In this study, we show that monkeys can time two concurrent intervals accurately with weakly dilated subjective time for the second overlapping interval.

interacting with a dynamic environment requires an accurate and flexible sense of timing. Accordingly, humans and other animals estimate, remember, discriminate, and reproduce intervals in the range of milliseconds to min (Lewis and Miall 2009), using their sense of timing for directing attention (Ghose and Maunsell 2002), anticipating and planning upcoming movements (Janssen and Shadlen 2005; Spencer et al. 2003), making decisions about future outcomes (Gibbon et al. 1988; Kim et al. 2008), and foraging optimally (Clayton et al. 2003; Kacelnik and Brunner 2002). Research on timing behavior has generated several neurobiologically plausible models of interval timing (Gibbon 1977; Gibbon et al. 1984; Jazayeri and Shadlen 2010; Karmarkar and Buonomano 2007; Matell and Meck 2004; Simen et al. 2011; Sohn and Lee 2013; Treisman 1963). These models capture the scalar property of timing, namely, the observation that the variability in timing estimates increases proportionately with the duration being timed (Gibbon 1977). Neurophysiological studies have found evidence of timing-related neural activity in multiple cortical and subcortical structures (Janssen and Shadlen 2005; Kim et al. 2013; Knudsen et al. 2014; Leon and Shadlen 2003; MacDonald et al. 2011; Maimon and Assad 2006; Matell et al. 2003; Merchant et al. 2011, 2013a, 2013b; Mita et al. 2009; Niki and Watanabe 1979; Shuler and Bear 2006). However, previous tasks have almost exclusively required subjects to estimate only one interval at a time, despite the observation that many human behaviors, such as timing different finger movements when playing an instrument or preparing and cooking multiple culinary dishes, may require concurrent timing. As a result, a fundamental question about interval timing, namely, whether there is one centralized clock or multiple distributed timers, remains unresolved (Buhusi and Meck 2005; Lewis and Miall 2006; Merchant et al. 2008; van Rijn and Taatgen 2008). Requiring the concurrent timing of multiple intervals permits these two hypotheses to be better differentiated and inferences to be made about the nature of the underlying temporal representations (Leak and Gibbon 1995; Meck and Church 1984).

To study timing of concurrent intervals, we developed a novel behavioral paradigm, the concurrent timing task, in which monkeys were required to time two different fixed temporal intervals. During this task, the interval between the two concurrent fixed intervals, or their onset asynchrony, varied randomly and widely across trials, discouraging the animals from adopting alternative strategies, such as timing intervals sequentially rather than concurrently, that could not be distinguished in previous studies (Leak and Gibbon 1995; Matell et al. 2003; Meck and Church 1984; Pang et al. 2001).

We found monkeys could time their behavioral responses to two different temporal intervals simultaneously even when they were presented concurrently, and we quantitatively show that their behaviors were consistent with a race between two independent stochastic processes, each related to a different interval. A small but systematic bias observed during concurrent timing could be accounted for by subjective time dilation of the second temporal interval, resembling the pattern observed in other timing tasks with competing attentional demands (Allman et al. 2014; Brown 1985; Brown and West 1990) or varied stimulus salience (Tse et al. 2004). These results suggest similar processes underlie single-interval and concurrent-interval timing and constrain neurobiological models of timing.

## MATERIAL AND METHODS

### Animal Preparation and Data Collection

Two male rhesus monkeys (*P* and *Q*; body weight = 11–13 kg) were used. A set of titanium head posts were attached to the skull during an aseptic surgery and used to immobilize the animals' heads during behavioral training and testing. Fluid intake by the animals was regulated to motivate the animals during behavioral training and testing. All procedures used in this study were approved by the Institutional Animal Care and Use Committee at Yale University and conformed to the *Guide for the Care and Use of Laboratory Animals*. A 16-in. cathode ray tube monitor was positioned 57 cm from the eyes of animals seated in a primate chair and used for the presentation of all visual stimuli in the experiment. An infrared eye tracking system was used to monitor the animal's eye position at 225 Hz (ET49; Thomas Recording, Giessen, Germany).

### Concurrent Timing Task

We designed a concurrent timing task that incentivized accurate time estimation for one or two fixed intervals on each trial. At the beginning of each testing session, a small white square (2.3° × 2.3°) with a white unfilled square around it (3.45° × 3.45°, line thickness = 0.38°) appeared at the center of the monitor and remained throughout the entire session (Fig. 1). Continuous fixation of this square (fixation window radius = 3.8° and 3.1° for *monkeys P* and *Q*, respectively) yielded a juice reward (0.1 ml for *monkey P*, 0.15 ml for *monkey Q*) after an unpredictable interval. This inter-reward interval was chosen according to an exponential distribution with a 1-s refractory period and an overall mean of 2 s. Reward from this central target introduced an opportunity cost for any eye movements away from the center. The flat hazard rate for the exponential distribution produced a constant probability of reward at any moment and therefore provided a continuous incentive for the animal to maintain its central fixation. Receipt of reward or eye movement away from the central target lasting longer than 1 s led to a new inter-reward interval to be drawn from the same distribution.

After a variable waiting time drawn randomly from a uniform distribution (3–6 s), a disk target (radius = 1.2°) appeared to the left or right of the central target along the horizontal meridian (eccentricity = 9.2°). Spatial position of the target varied pseudorandomly trial by trial. The color of this peripheral target indicated the duration of a fixed interval (FI) that had to elapse before reward could be collected from that target. In the present study, two different FI durations, 8 and 16 s, were used, indicated by red and yellow targets for *monkey P* and green and blue targets for *monkey Q*, respectively. To discourage the animals from making reflexive gaze shifts and to ensure that choice saccades to peripheral targets were distinct and punctate events, we also presented a white unfilled square (3.45° × 3.45°, line thickness = 0.38°) around a peripheral target. Saccade to a peripheral target without a white square around it was not registered as a legitimate choice and thus never delivered reward. A white square was illuminated around a peripheral target when the following conditions were met. First, white squares were absent during the first 0.5-s period after the onset of the peripheral target. Second, white squares were extinguished after receipt of reward from the central fixation target, unrewarded saccade to a peripheral target, or gaze anywhere outside of the central or peripheral target fixation window that lasted >1 s. Following any of these events, fixation of the central target for 0.5 s caused the white square to reappear around any peripheral targets. Fixation of the peripheral target that lasted longer than 0.5 s in the presence of the white square was registered as a choice (peripheral fixation window radius = 3.8° and 3.1° for *monkeys P* and *Q*, respectively). This procedure discouraged animals from making very-short-latency saccades to a peripheral target and then back to the fixation target. If the peripheral target was chosen before the end of the corresponding FI, no reward was delivered and the target remained visible and available for choice again after the requisite fixation of the central target. If the target was chosen after its FI, the target and its surrounding white square were extinguished after delivery of a large juice reward (0.35 ml for *monkey P*, 0.3 ml for *monkey Q*), 3.5 or 2 times larger than the reward from the central fixation reward. This larger reward for accurately timed saccades to peripheral targets than for continued fixation of the central target delivered a maximal reward rate to the animal only with accurate timing of saccades. The required 0.5-s fixation on a peripheral target after a choice saccade and subsequent 0.5-s fixation of the central target before the white square returned imposed an effective minimum interval between successive saccades of 1 s. Since juice reward was delivered in a series of multiple pulses (*n* = 7 and 6 for *monkeys P* and *Q*, respectively) with a fixed inter-reward interval (0.16 and 0.13 s for *monkeys P* and *Q*, respectively), a rewarded peripheral saccade introduced an additional delay of 1.1 s for *monkey P* and 0.8 s for *monkey Q*.

In single-interval trials, only a single peripheral target appeared, and the trial ended when the target was extinguished by an appropriately timed saccade. For concurrent-interval trials, the interval between the onset times of 8-s FI and 16-s FI was defined as the interval onset asynchrony (IOA). The IOA ranged from −7 s (8-s target presented 7 s before 16-s target) to +15 s (8-s target presented 15 s after 16-s target) in 1-s increments. Depending on the IOA, the minimum overlap between FI on each trial ranged from 1 s (i.e., IOA = −7 s with 8-s FI first, or 15 s with 16-s FI first) to 8 s (i.e., 0 s ≤ IOA ≤ 8 s with 16-s FI present for the entire duration of the 8-s FI). Because each peripheral target was extinguished by the saccade produced after the corresponding FI expired, the duration of the period in which the two targets remained visible together was almost always longer than the minimum overlap in the two FI. The wide range of IOA ensured that any potential nontiming heuristic strategies, such as always responding first to the first target to appear or to the shorter FI, were suboptimal. Single-interval and concurrent-interval trials were randomly interleaved within each block, which included 20 single-interval trials (5 of each FI with their positions counterbalanced) and 46 concurrent-interval trials (23 different IOAs with the positions of the 2 FI targets counterbalanced).

### Behavioral Training

Each animal was trained for the concurrent-interval timing task in several stages. Progression from one stage of training to the next was decided with the use of qualitative criteria, as described below. During the first four sessions (*stage 1*), the central fixation target was presented alone, and the animal was rewarded stochastically for maintaining its fixation. Reward volume started at 0.2 ml per fixation in this period and was reduced to 0.15 ml per fixation by the end of this stage. In *stage 2*, single-interval trials with only one target for the 8-s FI were presented. Training continued for 15 sessions for *monkey P* and 8 sessions for *monkey Q*, after which few responses were made early in the FI and response rates increased as time elapsed. In *stage 3*, reward was omitted in 20% of the trials, a procedure often referred to as “peak interval” trials. A peak interval trial included the identical 8-s FI stimulus, but it remained illuminated for 32 s regardless of the animal's responses and then was extinguished without reward. The purpose of peak trials was to determine whether the rate of saccades peaked around the time of reward availability at 8 s and then declined, which would be consistent with a learned association between the cue and reward availability after a particular interval. Peak trials were randomly interleaved with single-interval trials (*monkey Q*: 12 sessions; *monkey P*: 13 sessions). Although not reported in the current study but consistent with previous findings (Rakitin et al. 1998; Roberts 1981), at the end of *stage 3*, the response rate during the last 8 s of the peak trial was lower than the peak response rate, indicating that the animals began estimating the fixed interval, rather than simply increasing the frequency of response monotonically over time. In the next two stages, the animals were trained for single-interval trials with 16-s FI only, first without (*stage 4*) and then with (*stage 5*) peak trials, similar to *stages 2* and *3* for the 8-s FI. The duration of a peak trial was 64 s for the 16-s FI. The number of sessions for *stages 4* and *5* were 8 and 4 for *monkey P* and 4 and 17 for *monkey Q*, respectively.

In *stage 6*, we presented single-interval trials for both FI within the same session, first in blocks of 10–40 trials and then randomly interleaved. This stage persisted for 15 sessions for *monkey P* and for 1 session of interleaved trials for *monkey Q*. During this stage, the fixation reward volume was reduced to 0.1 ml for *monkey P*, whereas that for *monkey Q* remained at 0.15 ml per fixation. These reward volumes were used for the remaining sessions reported in this article. In *stage 7*, randomly interleaved peak trials were then introduced to ensure the animals reliably distinguished the FI even when they were not presented in blocks. At the end of *stage 7*, which lasted 5 sessions for *monkey P* and 14 sessions for *monkey Q*, response rates declined toward the end of peak trials for both FI in both animals (data not shown). In *stage 8*, concurrent-interval trials were introduced for a subset of IOAs (−8, −4, 0, 4, 8, 12, and 16) for 5 sessions for *monkey P* and 13 sessions for *monkey Q*. We next introduced a peak interval version of concurrent-interval trials (*stage 9*), in which the FI that would normally elapse later was a peak interval target for that trial. Reward was omitted only for one of the two targets in the peak trials, to prevent the omission of reward from one target from indicating the omission of reward from the other target, since this would have rendered analyses of response rates in peak trials less meaningful. For trials with IOA = 8 s, in which reward became available for both FI simultaneously, reward was omitted for one of the FI selected randomly. This stage lasted 9 sessions for *monkey P* and 41 sessions for *monkey Q*. Finally, in *stage 9*, the animals were trained for a full range of IOAs from −7 to 15 in 1-s increments, for 11 sessions for *monkey P* and 44 sessions *for monkey Q*, before the collection of the data analyzed in this study began (34 and 75 sessions for *monkeys P* and *Q*, respectively). The number of training sessions for *monkey Q* was larger, especially in the last two stages, because we explored alternative subsets of IOAs during this period.

### Data Analysis

The eye position data were smoothed using a five-point median filter followed by a Gaussian filter (σ = 16 ms). Saccade onsets were defined as the first time points in which instantaneous eye position velocity exceeded 30°/s. Peak saccade velocity for each saccade was defined as the maximum velocity reached at any time between when the velocity first rose above and subsequently slowed to below this. Saccade time was defined as the interval between the onset of a FI target and the onset of each choice saccade directed to it. In this report, the term “saccade time” is used, rather than saccade latency or reaction time, because the animals were allowed to make multiple saccades to the same target. Frequency histograms of saccade times were constructed by combining the data from all sessions to obtain average response rates as a function of time following the onset of a particular FI. In addition, the first saccade time in each trial and subsequent inter-saccade intervals (ISI) were analyzed separately. Because of the requirement for the fixations of the peripheral and central targets, as well as the time for reward delivery, there was a lower limit in the interval between successive choice saccades to peripheral targets. This task-imposed minimum interval (1.0∼2.1 s) was subtracted from the ISI for all the analyses. The Kolmogorov-Smirnov test was used to determine the statistical significance of the difference between cumulative frequency histograms of saccade latencies for the 8-s and 16-s FI, normalized by the FI duration, to test whether latencies conformed to the scalar property. For descriptive statistics used to compare saccades directed to the targets corresponding to the 8-s and 16-s FI, we calculated the median values from individual sessions and report their mean and standard error computed for the entire set of sessions, unless otherwise stated. For statistics related to the ISI for subsequent saccades in single-interval trials, we excluded sessions for which the animal had fewer than three trials with at least two saccades for each FI's target, to reduce the variability of the average ISI in individual sessions. For concurrent-interval trials, the trials were included in the analysis of first saccade times if there were at least two trials for a particular IOA in a given session, and in the analysis of ISI if there were at least two trials with multiple saccades to each FI.

### Probabilistic Models of Saccade Timing

#### Basic assumptions.

To characterize and compare the timing of saccades generated during single-interval and concurrent-interval trials, we treated each saccade as a stochastic process (Luce 1986). For each of the models described below, the likelihood L_{i} of observing the sequence of saccades with their corresponding saccade times in the *i*-th trial was given by the joint probability distribution of all *n* saccades in that trial,

where *N*_{i} is the number of saccades in the *i*-th trial, *T*_{k} refers to the time of the *k*-th saccade in the trial, *X*_{k} = 0 and 1 when the corresponding saccade was directed to the target for the 8-s or 16-s FI, respectively, and **θ** denotes the parameters of the model. For all the models tested in the present study, we assumed that the probability of observing the *k*-th saccade to a given target in a given trial was determined entirely by the elapsed time from the onset of that FI and the time, but not target, of the previous saccade. This allows the likelihood function given above to be written as the following:

where *P*(*T*_{k}, *X*_{k}|*T*_{k−1}; **θ**) denotes the conditional probability of observing the *k*-th saccade directed to the target *X*_{k} at time *T*_{k} given that the previous saccade was observed at *T*_{k−1}, and *P*(*T*_{1}, *X*_{1}; **θ**) is the probability of observing the first saccade to the target *X*_{1} at time *T*_{1}. The following sections describe how the probabilities *P*(*T*_{1}, *X*_{1}; **θ**) and *P*(*T*_{k}, *X*_{k}|*T*_{k−1}; **θ**), are specified for each of the models we tested, in both single-interval and concurrent-interval conditions. To simplify the description, we will omit **θ** in the notations below for these probabilities. The goodness of fit of a particular model was quantified by summing the natural logarithm of each saccade probability across trials. Namely, the log likelihood for each model across all trials, **LL**, is given by the following equation:

#### Single-interval condition: one-stage model.

In one-stage models, the probability of a saccade for the target *X* at any given time is based on a single underlying probability distribution function for all saccades, PDF_{All}, parameterized by a set of parameters denoted as **θ**_{X,All}. Thus

In addition, the conditional probability of the subsequent saccade at time *T*_{k} was given by the same distribution after it was normalized by the probability of observing another saccade after *T*_{k−1}, since there was always another saccade if the previous saccade occurred before the corresponding FI expired. Namely,

where CDF_{All} denotes the cumulative distribution function for PDF_{All}.

Various forms of density functions were tested for PDF_{All}, including lognormal, gamma, inverse gamma, and Weibull distributions. Each of these functions is specified by two parameters, with the saccade times for each FI fit separately, resulting in four free parameters for each model. The definition of each probability density function is given below:

For the gamma and inverse gamma probability density function, is the gamma function and is defined as follows:

#### Single-interval condition: two-stage model.

As shown in results, we observed distinct distributions for the first and subsequent saccade times. To investigate the difference between them quantitatively, we developed a two-stage model, in which saccade times were modeled by combining two different distinct stochastic processes. The time of the first saccade (FS) was determined by its own probability distribution function, PDF_{FS}, defined as a function of the elapsed time from the onset of a particular FI, whereas subsequent saccade times were determined by a mixture of PDF_{FS} and PDF_{ISI}, another probability distribution defined as a function of the ISI. In other words, the second process is reset after each saccade. Thus the probability of the first saccade at time *t* in a given trial was given by the following equation:

where **θ**_{x,FS} denotes the parameters of PDF_{FS} for the interval of that target *x*. Next, the probability of saccade times for subsequent saccades within a trial was given by the following:

where **θ**_{x,ISI} and ω_{x} denote the parameters and weight of PDF_{ISI}, respectively. Thus ω_{x} = 0 reduces this model to the one-stage model, whereas ω_{x} = 1 implies that the process of saccade generation undergoes a deterministic transition after the first saccade to a stage where all subsequent saccades within a trial follow the same distribution of ISI. The same four functional forms used in the one-stage model were tested for PDF_{FS} and PDF_{ISI}, and the mixture probability ω_{x} was fit separately for each FI, which resulted in 10 free parameters for each model.

#### Assumptions specific to the concurrent-interval condition.

For concurrent-interval trials, we assumed that each saccade observed in a given trial resulted from a race between two independent stochastic processes associated with the two alternative targets (Logan and Cowan 1984; Logan et al. 1984; Smith and Ratcliff 2004). The process that finished earlier generated a saccade to its target, after which both processes began a new race for control of the next saccade. In the following, the finishing times of each of these latent stochastic processes are indicated by *P*_{X}(*T* = *t*), denoting the probability that a saccade to the target for a given FI, *X* (equal to 0 and 1 for the 8-s and 16-s FI, respectively), would occur at time *t*, unless it is preceded by an earlier saccade to the target of the alternative FI, ∼*X* (i.e., ∼*X* = 0 and 1 if *X* = 1 and 0, respectively). It should be noted that for the single-interval condition, these latent stochastic processes entirely determine the probability for the observed saccade sequence, namely, *P*(*T*_{1}, *X*_{1}) = *P*_{X1}(*T*_{1}) and *P(T*_{k}, *X*_{k}|*T*_{k−1}) = *P*_{Xk}(*T*_{k}|*T*_{k−1}). By contrast, for the concurrent-interval condition, these two probabilities, *P*(*T*_{1}, *X*_{1}) and *P*(*T*_{k}, *X*_{k}|*T*_{k−1}), can be computed by a product of the two probabilities characterizing these underlying stochastic processes. The first is the probability that a saccade would be directed to *X*_{k} at time *T*_{k} without the race and corresponds to *P*_{X1}(*T*_{1}) and *P*_{Xk}(*T*_{k}|*T*_{k−1}) for the first and subsequent saccades in a given trial, respectively. The second is the probability that no saccade would be generated to the other target and is denoted as *P*_{∼X1}(*t* > *T*_{1}) and *P*_{∼Xk}(*t* > *T*_{k}|*T*_{k−1}) for the first and subsequent saccades, respectively. With these two additional probabilities, *P*(*T*_{1}, *X*_{1}) and *P*(*T*_{k}, *X*_{k}|*T*_{k−1}) are now given by the following:

and

#### Concurrent-interval condition: one-stage model.

Analogous to the one-stage model for the single-interval condition, the latent process for the saccades directed to the target of a given FI in the one-stage model for the concurrent-interval condition is derived from a single probability density function. However, the probability for each saccade is computed differently, depending on whether it is the first saccade in a given trial or not. First, the latent process for the first saccade and subsequent saccades in each trial is given by the following two equations:

and

where *T*_{k−1} was the time of immediately preceding saccade to either target.

Second, the probability that no other saccade would occur to the other interval before the first saccade in the trial is given by the following:

The probability that no saccade would occur to the other interval during the ISI preceding a subsequent saccade at *T*_{k} is given by the following:

As in the single-interval condition, we tested each of the four distributions (lognormal, gamma, inverse gamma, and Weibull). We fit one set of parameters, **θ**_{X,All}, for the PDF_{All} for each FI across all IOAs, resulting in four free parameters for each model.

#### Concurrent-interval condition: objective two-stage model.

Similar to the one-stage model described above, it is also assumed in two-stage models for the concurrent-interval condition that each saccade is determined by a race between two stochastic processes. In contrast to the one-stage model, however, each of these stochastic processes in the two-stage model is characterized by two different probability density functions for the first saccade and subsequent ISI, as in the two-stage model for the single-interval condition. When the first saccade to a given target was also the first saccade in a concurrent-interval trial, its latent process was determined by the following:

If the first saccade to a given target was preceded by another saccade to the other target, its latent process was determined by the following:

For all other saccades, the latent process was determined by the following:

where the mixture parameter ω_{Xk} controls how much the time of the *k*-th saccade was determined by the probability distribution for the ISI.

Next, the probability of no saccade to the other interval before the first saccade in a concurrent-interval trial does not depend on the probability distribution for the ISI and therefore was determined as follows:

Similarly, even for subsequent saccades within a trial, if no saccade had yet occurred for the other interval, the probability of no saccade to the other target during the previous ISI was given by the following, since this also does not depend on the probability distribution for the ISI:

By contrast, if at least one saccade had already occurred for the other FI, then this probability now depended on the probability of the ISI and was given by the following:

where CDF_{ISI} denotes the cumulative density function of the ISI for the target *X*. As with the single-interval two-stage model, we tested each of the four distributions for both PDF_{FS} and PDF_{ISI}, with separate parameters **θ**_{X,FS} and **θ**_{X,ISI}, and mixture parameters, ω_{X}, for each FI, resulting in 10 free parameters for each model.

#### Concurrent-interval condition: subjective two-stage model.

To test how a distortion in the internal representation of time could influence responses in concurrent-interval trials, we tested a variant of the two-stage model described above in which the probability of saccade is determined by a function of subjective time, *u*(*t*), rather than the physical time, *t*, itself. This transformation from physical time to subjective time was applied only to the saccades directed to the second interval in the concurrent-interval condition for two reasons. First, as shown in results, changes in the IOA during the concurrent-interval condition mainly influenced the distribution of saccade times for the second interval to appear. In addition, a large portion of the first interval in the concurrent-interval condition was physically indistinguishable from the single-interval condition.

We included both additive and multiplicative distortions in the subjective time. First, subjective time could be lengthened or shortened by a fixed amount proportional to the absolute value of the IOA in a given trial. Second, subjective time could be dilated or contracted increasingly as time elapsed, compared with physical time. These two effects were combined into the following equation:

where λ and φ are parameters governing the degree of multiplicative time dilation/contraction and additive translation, respectively, and Δ_{T} = 0 for the first interval to appear and |IOA| for the second interval. The value of λ was constrained to be greater than and for the 8-s and 16-s FI, respectively, so that *u*(*t*) increased monotonically with *t* for all values of Δ_{T} used in the task.

The two latent processes responsible for saccades directed to each FI in a concurrent-interval trial with subjective time were determined similarly as in the objective two-stage model. Thus the latent process for the first saccade in a given trial was given by the following:

where the denominator and the coefficient *c* = 1 + λΔ_{T} ensured that *P*_{X1}(*T* > 0) = 1.

If the first saccade to a given target was preceded by another saccade to the other target, its latent process was given by the following:

For all other saccades,

The probability of no saccade to the other interval before the first saccade was determined similarly to the objective two-stage model but using the subjective time, *u*(*t*), as given by the following:

For subsequent saccades within a trial, the following equation was used to determine the probability of no saccade to the other interval during the previous ISI, if no saccade had yet occurred for the other interval:

Finally, if at least one saccade had already occurred for the other FI, then this probability was given by the following:

We tested eight different models that incorporated subjective time. Each had the same 10 free parameters as in the objective two-stage concurrent-interval model, along with 1–4 additional parameters testing whether the subjective time parameters λ and φ were necessary and distinct for the 8-s and 16-s FI. Each model tested a different combination of λ and φ terms: fixed at zero, shared between FI, and separate for different FI.

### Model Selection

Maximum-likelihood estimation was implemented in MATLAB (The MathWorks) using the fminsearch function and custom-written code and was used to fit the parameters for each of the models described above. Model selection used the Bayesian information criterion (BIC) = −2**LL** + *h·*ln(*N*), where **LL** is the log likelihood of the model, *h* is the number of free parameters, ln is the natural logarithm, and *N* is the total number of saccades across all sessions. For visualization of each best-fitting model in concurrent-interval trials, we used the best-fitting parameters to simulate 1,000 trials for each IOA. The predicted response rate, median first saccade time, and median ISI values were calculated from these simulated trials.

## RESULTS

### Effects of Single Fixed Intervals on Saccade Rates

Behavioral data analyzed in this study (*monkey P*, 2,714 trials in 34 sessions; *monkey Q*, 7,260 trials in 75 sessions) were obtained after all training stages (see methods) were complete. The timing task used in this study (Fig. 1) included single-interval trials in which only one peripheral target, corresponding to either an 8-s (429 and 1,147 trials for *monkeys P* and *Q*, respectively) or a 16-s FI (421 and 1,152 trials for *monkeys P* and *Q*, respectively), was presented. The results from these single-interval trials allowed us to characterize the baseline timing behavior in the task and showed that both animals responded differentially to the targets for the 8-s and 16-s intervals (Fig. 2). Averaged over all trials, response rates to the 8-s target rose sooner and more sharply than to the 16-s target (Fig. 2, *A–D*). Accordingly, mean response rates in a 0.5-s window beginning at 8 s from interval onset were significantly higher for the target for the 8-s FI than for the target for the 16-s FI. We calculated the mean response rate within each session and report the mean ± SE across sessions. For *monkey P*, the mean response rate in this window was 0.176 ± 0.041 and 0.07 ± 0.02 saccades/s for the 8-s and 16-s target, respectively (paired *t*-test, *P* < 0.05). For *monkey Q*, the mean response rate in this window was 0.174 ± 0.019 and 0.044 ± 0.011 saccades/s for the 8-s and 16-s target, respectively (paired *t*-test, *P* < 10^{−6}). The time courses of saccade rates approximately conformed to the scalar property, in that they were more similar when plotted as a function of the elapsed time normalized by the corresponding FI (Fig. 2, *E* and *F*). However, the difference was still significant in one animal (Kolmogorov-Smirnov test, *monkey P*: *P* = 0.93; *monkey Q*: *P* < 0.005).

We found that the difference in the saccade time distributions for the two FI largely arose from the difference in the onset times of the first saccades, and that the ISI for the subsequent saccades within a trial were more similar for the two FI. For both animals, first saccades tended to occur earlier to the 8-s FI target than to the 16-s FI target (Fig. 3, *A* and *B*). For *monkey P*, the median first saccade times for the 8-s and 16-s FI were 9.38 ± 0.397 and 11.84 ± 0.548 s, respectively (paired *t*-test, *P* < 0.005). For *monkey Q*, the median first saccade times for the 8-s and 16-s FI were 10.65 ± 0.233 and 14.53 ± 0.281 s, respectively (paired *t*-test, *P* < 10^{−23}). Variability in first saccade times was also smaller for the 8-s FI. For *monkey P*, the standard deviation of first saccade times for the 8-s and 16-s FI was 6.12 ± 0.489 and 7 ± 0.429 s, respectively (paired *t*-test, *P* < 0.05). For *monkey Q*, the standard deviation of first saccade times for the 8-s and 16-s FI was 4.12 ± 0.173 and 4.5 ± 0.209 s, respectively (paired *t*-test, *P* < 0.05). The coefficient of variation (CV) for first saccade times for each FI, namely, the ratio of the standard deviation to the mean, was not significantly different for *monkey P*, with mean CV for the 8-s and 16-s FI of 0.58 ± 0.037 and 0.57 ± 0.034, respectively (paired *t*-test, *P* = 0.81). For *monkey Q*, the mean CV for the 8-s and 16-s FI was 0.39 ± 0.015 and 0.3 ± 0.011, respectively, and this difference was significant (paired *t*-test, *P* < 10^{−5}).

Second and later saccades occurred with shorter ISI than the average latency of first saccades, and median ISI was similar for both FI (Fig. 3, *C* and *D*). For *monkey P*, the median ISI for the 8-s and 16-s FI was 1.81 ± 0.204 and 2.13 ± 0.319 s, respectively (paired *t*-test, *P* = 0.23). For *monkey Q*, the average median ISI for the 8-s and 16-s FI was 1.25 ± 0.173 s and 1.07 ± 0.117 s, respectively (paired *t*-test, *P* = 0.39). With these results taken together, first saccade times to the 8-s and 16-s FI targets were clearly distinguishable, whereas second and later saccades were emitted with similar ISI for both FI.

### A Two-Stage Model for Single-Interval Trials

To investigate more quantitatively how the observed patterns of saccade times might reflect the internal representation of elapsed time, we tested two main categories of models. In one-stage models, it was assumed that all saccades were generated by a single underlying stochastic process, whereas in two-stage models, it was assumed that first saccades in a given trial to each FI and subsequent ISI might be characterized by two different stochastic processes. The probability density functions that describe these two stochastic processes are denoted as PDF_{FS} and PDF_{ISI}, respectively (see methods). For all functional forms used for these two probability density functions, the two-stage models performed better than the one-stage models (Table 1). The Bayesian information criterion (BIC) for the worst two-stage and best one-stage model was 9,8991.1 and 10,005, respectively, for *monkey P* and 19,902.4 and 20,230.8, respectively, for *monkey Q*. For the two-stage models, both the gamma and Weibull distributions performed similarly well for PDF_{FS}, whereas the inverse gamma function for PDF_{ISI} performed best overall for both animals (Table 1). The mixture parameter ω_{X} determined the contribution of PDF_{ISI} to second and later saccades for each FI. For the best-fitting two-stage model, the value of ω_{X} for *monkey P* was 0.84 and 0.67 for the 8-s and 16-s FI, respectively, and the value for *monkey Q* was 0.7 and 0.76 for the 8-s and 16-s FI, respectively (Table 2). These results indicated that saccades after the first to each FI were determined more strongly by the distribution for the ISI than by the distribution referenced to the onset of the fixed interval for both animals.

### Effects of Concurrent Fixed Intervals on Saccade Rates

The majority of trials (*monkey P*: 1,864; *monkey Q*: 4,961) in this experiment were concurrent-interval trials, in which the 8-s and 16-s FI at least partially overlapped (Fig. 1). For both animals, the distribution of saccade times differed between the 8-s and 16-s FI for each IOA, similar to the difference in saccade times for each target in single-interval trials (Fig. 4). The trial-averaged response rates were qualitatively similar to those for single-interval trials (Fig. 5). As in single-interval trials, the observed difference in saccade times was due to earlier onset of the first saccade to the 8-s than to the 16-s target (Fig. 6; paired *t*-test, *monkey P*: *P* < 10^{−11}; *monkey Q*, *P* < 10^{−50}), whereas subsequent ISI were not significantly different between the two targets (paired *t*-test, *monkey P*: *P* = 0.529; *monkey Q*: *P* = 0.065).

In comparing single-interval with concurrent-interval timing behavior, we found systematic changes in saccade times for a subset of IOAs (Fig. 7). Specifically, the time of the first saccade for the second FI target tended to decrease as the absolute value of the IOA increased (Fig. 7, *A* and *B*). We calculated median first saccade times for each IOA within each session and then examined the correlation between IOA and median first saccade time across sessions. When the 16-s FI appeared first (i.e., IOA > 0), first saccade times to the 8-s FI were significantly and negatively correlated with IOA (Pearson's correlation, *monkey P*: *r* = −0.21, *P* < 10^{−4}; *monkey Q*: *r* = −0.2, *P* < 10^{−9}). Similarly, when the 8-s FI appeared first (i.e., IOA < 0), first saccade times to the 16-s FI were negatively correlated with the absolute value of IOA, although this was only observed for *monkey P* (Pearson's correlation, *monkey P*: *r* = −0.2, *P* < 0.01; *monkey Q*: 16-s FI: *r* = 0.023, *P* = 0.62). By contrast, the ISI for subsequent saccades within a trial were more consistent across IOAs for both the 8-s and 16-s FI targets (Fig. 7, *C* and *D*). Median ISI exhibited no significant systematic differences across IOAs for either FI for either animal (all *P* > 0.1).

### Two-Stage Models for Concurrent-Interval Trials

We explored whether the same modeling framework that accounted well for saccade times for single-interval trials could be extended to account for the behaviors in concurrent-interval trials. Specifically, we assumed each saccade in concurrent-interval trials resulted from a race between two independent processes that were each individually responsible for saccades to 8-s and 16-s targets. Each process was defined as in the corresponding single-interval trials, without introducing any additional processes necessary for timing multiple intervals.

First, we fit the one-stage and objective two-stage models to concurrent-interval trials, testing the same four distributions used to analyze the results from single-interval trials. The one-stage model assumed a single distribution described all saccade times, whereas the objective two-stage model assumed that the timing of saccades after the first to each FI was determined by a mixture of two distributions. Both models fit one set of parameters for all IOAs and therefore assumed no systematic change in the underlying process for each FI as a function of IOA. For all distributions, the two-stage model, with one distribution, PDF_{FS}, determining first saccade times and a second distribution, PDF_{ISI}, that contributed to second and later saccades to each FI, fit better than the one-stage model, in which all saccades were determined by one PDF (Table 3). The gamma distribution and inverse gamma fit best for PDF_{FS} and PDF_{ISI}, respectively, for both animals.

We next tested the subjective two-stage model, which allowed subjective time to deviate from physical time. We examined whether dilation or expansion of subjective time could account for the progressively earlier first saccade times observed for the second interval in concurrent-interval trials as the absolute value of the IOA increased (quantified by λ). We also hypothesized that the estimate of the time of the second interval onset may be influenced by the remembered onset time of the first interval (quantified by φ). Model comparison by BIC showed that the results from both animals were best explained by the models including subjective time that deviated from physical elapsed time (Fig. 8, Table 4). Both animals had positive λ parameters, with a single λ term sufficient for both FI in both monkeys, indicating similar subjective time dilation for both FI (Table 5). *Monkey Q* additionally exhibited subjective time translation that was equivalent for both FI, with a shared φ parameter in the best-fitting model. *Monkey P* showed no such translation, with saccade times best fit by a model with φ = 0. Similarly to the objective two-stage model, the gamma and inverse gamma fit best for PDF_{FS} and PDF_{ISI} in the subjective two-stage model, respectively, for both animals (Table 3). Taking these results together, we found that concurrent-interval timing was characterized by similar underlying processes as single-interval timing, requiring only the addition of subjective time dilation during concurrent timing.

## DISCUSSION

### Concurrent Timing

We developed a novel behavioral task to study whether monkeys can time two intervals concurrently. The animals distinguished between 8-s and 16-s FI, even when the intervals overlapped temporally. This was reflected in earlier initial responses to the target for the shorter FI than for the longer. The timing of the first saccade to each target in a trial was a reliable measure of timing discrimination, whereas subsequent saccades were emitted with short and less flexible ISI. This pattern of responding suggests the animals withheld responses until the subjective estimate of elapsed time crossed some threshold, after which responses were emitted continuously and stereotypically until reward was received, consistent with previous examinations and models of timing behavior in individual trials (Church et al. 1994; Mello et al. 2015; Rakitin et al. 1998). Despite response rates being approximately scalar (Fig. 2, *E* and *F*), the distribution of first saccade times to each FI was neither scalar nor normally distributed (Fig. 3, *A* and *B*) but was positively skewed and significantly different for each FI (Hanson and Killeen 1981; Mello et al., 2015). As measured by the time of the first saccade, temporal estimates for the shorter FI tended to be slightly overestimated, and those for the longer FI tended to be slightly underestimated, consistent with Vierodt's law of timing (Lejeune and Wearden 2009) and the temporal context of the experiment leading to prior expectations about durations (Jazayeri and Shadlen 2010; Jones and McAulery 2005). In concurrent-interval trials, the times of first saccade were progressively earlier for the second target to appear as the delay between first and second interval onset lengthened.

We found that responses were best described by two-stage models, with the gamma distribution determining first saccade times and subsequent ISI determined by a mixture of the gamma and inverse gamma distributions in which the inverse gamma predominated. The gamma distribution has been widely used to model time-dependent processes (Anderson 1974; Killeen and Fetterman 1988; Ratcliff and Murdock 1976), including the waiting time to respond after onset of a fixed interval (Hanson and Killeen 1981). However, too much theoretical significance should not be attached to the identity of the function, because different distributions with appropriately chosen parameters can be indistinguishable given the sample sizes and behavioral measurements used (Van Zandt and Ratcliff 1995). In the present study, multiple functions were fitted and compared to ensure the model features were not an artifact of an inappropriate distribution.

In concurrent-interval trials, modeling saccades as a race between the same processes operating in single-interval trials provided a good fit to the data. Despite differences in training and the reward volumes delivered for the central and peripheral targets for each monkey, the behavior of both animals was nevertheless fit well by remarkably similar models, with a two-stage process generating saccades and deviations of subjective time from objective time that occurred for the second interval to appear in each concurrent timing trial. Rather than evidence of a disruption specific to concurrent timing of a second interval, the subjective time dilation observed in the present study is consistent with the disruption expected from another competing, nontemporal demand by attentional or perceptual factors, such as visual stimulus salience or motion (Tse et al. 2004; Wittmann et al. 2010) or the difficulty of the task being performed (Brown 1985; Brown and West 1990). Behavior during concurrent timing can therefore be modeled without requiring any fundamental change to the processes of interval timing between single and concurrent conditions. These results provide evidence that similar processes underlie timing in single-interval and concurrent-interval conditions, which imposes strong constraints on neurobiological models of timing.

### Relationship to Previous Studies

In previous studies on concurrent timing of multiple intervals, behavioral tasks were based on fixed temporal relationships between the timed intervals (Kirkpatrick and Church 2000; Leak and Gibbon 1995; Matell et al. 2003; Meck and Church 1984; Pang et al. 2001). In many studies, one FI was always embedded within another, longer FI, usually with simultaneous onsets. Without the inclusion of variable onset asynchronies between the intervals, the possibility that animals might have learned habitual patterns of responses could not be excluded. For example, a task presumably requiring concurrent timing, with simultaneous onset of a 10-s FI and a 40-s FI, could possibly be timed with satisfactory precision by sequentially timing a 10-s interval followed by a 30-s interval (Matell et al. 2003). Similarly, varying the interval between onsets of the FI was required to uncover the observed systematic effect of IOA on first saccades. In single-interval timing tasks, the end of each interval usually co-occurs with the expected reward and behavioral response, complicating the dissociation of reward anticipation and movement preparation from time estimation (Lewis and Miall 2006). Because many brain regions in which timing-related activity has been reported also show responsiveness to reward and motor planning (Fiorillo et al. 2008; Janssen and Shadlen 2005; Knudsen et al. 2014; Maimon and Assad 2006; Matell et al. 2003; Merchant et al. 2011; Mita et al. 2009), introducing concurrent timing is useful for allowing responses and reward for one interval while timing must be continued for another interval. Neural activity specifically related to timing of an interval must be robust to such intervening responses and rewards. Our novel task therefore provides a useful test for positing links between timing and neural activity.

### Implications for Pacemaker-Accumulator Models

The finding that monkeys accurately time simultaneous intervals with trial-by-trial changes in the onset time of each FI imposes substantial constraints on neurobiological models of timing. Any proposed neural mechanism for timing must account for the near independence observed between timing of each interval while allowing for subjective time dilation to occur selectively (Tse et al. 2004). One of the most experimentally supported timing models is the pacemaker-accumulator model (Church et al. 2003; Gibbon et al. 1984; Treisman 1963). A pacemaker emits pulses, which are accumulated during interval timing. The accumulated pulse count is compared with stored/remembered pulse counts that were previously associated with a specific duration. When these counts are sufficiently similar, the remembered interval is considered elapsed and behavioral responses begin to be emitted (Church et al. 1994). Physiological studies have found evidence of ramplike neural activity that could reflect such accumulators (Kalenscher et al. 2006; Knudsen et al. 2014; Komura et al. 2001; Matell et al. 2003; Merchant et al. 2011; Mita et al. 2009; Niki and Watanabe 1979), although to our knowledge such ramping activity has not been investigated during concurrent timing.

Several modifications to the pacemaker-accumulator model would render it adept at incorporating concurrent timing of independent intervals. First, for each timed interval, a separate accumulator is necessary. Neurophysiological correlates of this accumulation, such as ramping activity, should therefore be specific to only one of the intervals. Second, time dilation that was interval-specific could result from superadditive accumulation of pulses in the second interval's accumulator. Theoretical work has postulated the presence of an arousal signal that could act on the pacemaker to modulate its pulse rate (Treisman 1963). This mechanism would not support interval-specific time dilation in concurrent-interval trials, but another modulatory signal might affect a later stage that is accumulator specific, such as a retrieval or comparator process. Therefore, pacemaker-accumulator models of timing can be revised to accommodate concurrent timing.

### Implications for Alternative Timing Models

Like the pacemaker accumulator, alternative models of timing have been supported by behavioral and physiological results. However, concurrent timing poses a challenge to current formulations of many of these models. One alternative model of timing in neural networks uses state-dependent network trajectories to encode duration (Buonomano and Merzenich 1995; Karmarkar and Buonomano 2007). At the start of an interval, the network begins a stereotyped trajectory associated with that duration. The elapsed time is read out on the basis of where the population is in the state space. One feature of this model is that temporal information is read out from a population of neurons that individually might not track time with high fidelity (Buonomano and Laje 2010), providing an attractive mechanism by which sensory, motor, and association areas represent elapsed time while nontemporal information continues to be processed.

Another proposed timing mechanism relies on neural oscillators with variable frequencies, which are synchronized at the beginning of an interval and convey timing information to downstream neurons based on the relative phases or activity levels of the oscillators (Allman and Meck 2012; Matell and Meck 2004; Meck et al. 2008; Miall 1989). In this so-called striatal beat-frequency model, cortical oscillators are synchronized by dopamine release after a reward-predicting stimulus and read out by the striatum (Matell and Meck 2004). Although no direct evidence for this striatal beat-frequency model has been provided, pharmacological and lesion studies suggest the dopaminergic system has some role in tracking time in the seconds range (Meck 2005, 2006; Meck et al. 2008).

State-dependent networks and the striatal beat-frequency model rely on state or synchrony reset at the beginning of each timed interval. The onset of the second interval presents a critical problem. For the state-dependent network, at this time the first interval has been encoded by the ongoing population state trajectory, which is now interrupted by an additional interval that must be encoded in the trajectory without disrupting the encoding of the first interval. For the beat-frequency model, accurate time readout is dependent on the phase reset of the oscillators, but resetting them to begin timing the second interval abolishes the representation of the first interval.

In theory, several additions to these models could allow them to perform simultaneous timing. First, concurrent intervals could be represented by independent populations of neurons, predicting orthogonal state representations in state-dependent networks and separate oscillators in the beat-frequency model. This would sidestep the problems related to timing multiple intervals with a single population but would require new populations of neurons to be recruited to time additional intervals. Physiological experiments should then observe distinct populations of neurons encoding concurrent intervals, or a partial remapping of a subset of neurons representing the first interval to begin representing the second interval to appear. State-dependent networks might alternatively represent simultaneous intervals as a compound cue. The state trajectory representing a single FI would be unrelated to the trajectory representing the same FI after the onset of another FI. Similarly, the trajectory after the second FI onset would be distinct for different delays between first and second interval onset. Only by distinguishing between each IOA would such a model be able to accurately time the simultaneous FI. Although learning and decoding such trajectories may not be trivial, this mechanism could encode simultaneous intervals. Similarly, the beat-frequency model could time two intervals by allocating subsets of oscillators to each interval. This predicts reduced timing precision during concurrent intervals and more localized dopamine release and phase resetting for the second timing cue relative to the first timing cue.

Neurophysiological experiments using this and similar tasks will be necessary to uncover neural correlates of elapsed time per se, independent of activity related to other behaviorally relevant variables that typically covary with time. In addition, further experimental and theoretical work will be critical to determine whether there is neurobiological support for modified timing models that incorporate concurrent timing. Indeed, uncovering more generally the neural basis of concurrent timing will be crucial for fully understanding the neural representation of time.

## GRANTS

This work was supported by National Institutes of Health Grants F31MH099706 (to M. R. Kleinman) and R01DA029330 and R01MH108629 (to D. Lee).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

M.R.K. and D.L. conception and design of research; M.R.K. performed experiments; M.R.K. and D.L. analyzed data; M.R.K., H.S., and D.L. interpreted results of experiments; M.R.K. prepared figures; M.R.K. and D.L. drafted manuscript; M.R.K., H.S., and D.L. edited and revised manuscript; M.R.K., H.S., and D.L. approved final version of manuscript.

- Copyright © 2016 the American Physiological Society