A sensory integration account for time perception

The connection between stimulus perception and time perception remains unknown. The present study combines human and rat psychophysics with sensory cortical neuronal firing to construct a computational model for the percept of elapsed time embedded within sense of touch. When subjects judged the duration of a vibration applied to the fingertip (human) or whiskers (rat), increasing stimulus intensity led to increasing perceived duration. Symmetrically, increasing vibration duration led to increasing perceived intensity. We modeled real spike trains recorded from vibrissal somatosensory cortex as input to dual leaky integrators–an intensity integrator with short time constant and a duration integrator with long time constant–generating neurometric functions that replicated the actual psychophysical functions of rats. Returning to human psychophysics, we then confirmed specific predictions of the dual leaky integrator model. This study offers a framework, based on sensory coding and subsequent accumulation of sensory drive, to account for how a feeling of the passage of time accompanies the tactile sensory experience.

Introduction Every sensory experience is embedded in time, and is accompanied by the perception of the passage of time. The coupling of the perception of the content of a sensory event and the time occupied by that event raises a number of questions: Do these percepts interact with each other? Do they emerge within separate neuronal populations? Which neuronal mechanisms underlie the generation of two distinct percepts? By comparing and contrasting human and rat psychophysics, the present study addresses these questions and thus aims to uncover generalized mechanisms through which the brain represents the passage of time.
Two principal frameworks have been proposed to explain the neuronal bases of the feeling of time in time scales of up to about 1 s: one framework posits a central clock, not connected with any specific sensory modality [1] while a second framework posits that the cortical circuit associated with each modality intrinsically encodes the passage of time for events within that modality [2,3]. There are also mixed models that argue for the existence of a core timing structure that integrates cortical activity in a context-dependent way [4,5].
To determine to what extent the coding of time is embedded within the coding of the stimulus itself, here we examine the relationship between the perceived features of a sensory event and the perceived duration of that same sensory event. As a stimulus feature, we focus on the intensity of tactile vibrations (quantified by their mean speed). The psychophysical experiments reveal a systematic interaction between perceived intensity and perceived duration, both in humans and in rats, leading us to propose that an early-stage sensory representation provides common sensory input to two downstream integrators that generate two corresponding percepts. To test this model, we use neuronal activity recorded from somatosensory cortex of behaving rats to generate neurometric curves for perceived intensity and perceived duration, and derive a close match to observed psychophysical curves. From these findings, we propose a framework, general to humans and to rats, for the construction of the feeling of the intensity of a stimulus and the time occupied by that same stimulus.

Results
We carried out four experiments to investigate the effect of two tactile vibration featuresintensity and duration-on the two percepts directly connected to those features-perceived intensity and perceived duration. Experiments 1 and 2 involve both human subjects, to whom stimuli were delivered to the left index fingertip, and rats, to whom vibrations were delivered to the whiskers on the right side of the snout (Fig 1A). Psychophysical experiments point to a candidate mechanism for the generation of the percepts-two accumulators of sensory drive, operating in parallel with percept-specific parameters of integration-and we probe the feasibility of the posited mechanism in Experiment 3 by generating intensity and duration neurometric curves from recorded rat somatosensory cortex firing. Experiment 4 tests additional predictions of the model in human subjects.
(equivalent to the standard deviation of the Gaussian multiplied by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ð2=pÞ p ). We consider perceived intensity to be the subjective experience related to objective intensity. Each stimulus was also defined by its duration in ms, T. Experiment 1 employed a delayed comparison task, also known as the two interval forced choice task (Fig 1C). On each trial, subjects received two vibrations (Stimulus 1, Stimulus 2), separated by a fixed delay (500 ms for human subjects, 2 s for rats). The experiment was comprised of two distinct tasks: for duration delayed comparison, the subject had to judge which of the two stimuli was longer according to the relative T values (T1 > T2 or T2 > T1). For intensity delayed comparison, the subject had to judge which of the two stimuli was of greater intensity (I1 > I2 or I2 > I1). On trials when the parameter of interest was equal, the correct (rewarded) answer was assigned randomly. Each of 10 human subjects carried out both tasks, on different days, while individual rats were trained on a single task: 7 were intensity rats and 7 duration rats.
To constrain subjects to rely on working memory, we used a set of stimulus pairs referred to as the stimulus generalization matrix (SGM; Fig 2) in which any value of I1 could be followed by a larger or smaller I2 and any value of T1 could be followed by a larger or smaller T2 [6,8]. Since neither stimulus alone provided the information necessary for a correct choice, both stimuli had to be attended to and utilized to solve the task.
In each trial, the two stimuli could differ in I, in T, or both. To quantify stimulus differences, we designated two indices. The normalized I difference (NID), defined as (I2 -I1) / (I2 + I1), ranged from -0.3 to 0.3 for both humans and rats, while the normalized time difference (NTD), (T2 -T1) / (T2 + T1), ranged from -0.3 to 0.3 for humans and from -0.35 to 0.35 for rats. The stimulus set was constructed to present every possible combination of NTD and NID values during the session (Fig 1D). Subjects received the same stimuli whether the task was to The left side shows two traces of sampled speed over time, and the right side shows the folded half Gaussian distribution to which each sample corresponds. The distribution's expected value is shown for each trace. C) Delayed comparison trial structure. Each trial consisted of the presentation of two noisy stimuli, with specified durations and intensities, separated by an interstimulus delay. The response was deemed correct according to the task rule: to compare the relative durations (blue-shaded rule) or relative intensities of Stimulus 1 and 2 (red-shaded rule). D) Representation of all possible stimulus intensities and durations presented to the subjects in the delayed comparison task. Each square in the matrix is color coded according to the NTD and NID of the two vibrations presented. Selected NID/NTD combinations from the top left and bottom right of the matrix are illustrated.
https://doi.org/10.1371/journal.pcbi.1008668.g001 PLOS COMPUTATIONAL BIOLOGY judge intensity or duration (see Fig 2 for the set of intensity and duration values). Thus, any resulting difference in performance of the tasks could not be attributed to differences in tactile input.
The upper plot of Fig 3A shows results from the duration delayed comparison sessions, with human data given as solid lines and rat data as dashed lines. The likelihood of judging T2 > T1 is plotted as a function of NTD, and the resulting steep curves report the capacity of human and rat subjects to extract the elapsed time. The similarity between the curves obtained from rats and human subjects demonstrates that the two species are approximately equally as proficient at discriminating stimulus durations. Matrices used for the intensity delayed comparison task. Lower row: matrices used for the duration delayed comparison task, for human subjects. Each trial's pair of task relevant feature values (I in the intensity task, T in the duration task) was drawn uniform randomly from the set of pairs scattered in the leftmost plots. Each trial's pair of task irrelevant feature values (T in the intensity task, I in the duration task) was drawn uniform randomly from the set of pairs scattered in the rightmost plots. B) Same as A, for rat subjects.

PLOS COMPUTATIONAL BIOLOGY
By the same layout, the lower plot of Fig 3A show results from the intensity delayed comparison sessions. When the likelihood of judging I2 > I1 is plotted as a function of NID, the resulting steep curves report the capacity of human (solid) and rat (dashed) subjects to extract stimulus intensity. Again, the two species show similar capacities in discriminating stimulus intensities. Fig 3B shows the overall performance achieved by humans (left) and rats (right) in the two delayed comparison tasks. The left bar of each plot depicts the percentage of correct trials obtained when the subjects' performance is analyzed by the intensity rule, while the right bar depicts the correct percentage when the subjects' performance is analyzed by the duration rule. The upper panels show that the two species had similar performance (75-80% correct) in duration delayed comparison sessions when their choices were measured according to the duration rule. However, if their choices were measured according to the relative intensities of the stimuli, performance within the same stimulus set would remain above chance (about 55% correct). The small but significant bias according to the non-relevant stimulus feature means that in both species the higher-intensity stimulus, on average, tends to be judged as longer in duration. That is, stronger feels longer.
Similarly, the lower panels show that the two species had similar performance (75-80% correct) in intensity delayed comparison sessions when their choices were measured according to the intensity rule. However, if their choices were measured according to the relative durations of the stimuli, performance within the same stimulus would remain above chance (about 55% correct). Again, the small but significant bias according to the non-relevant stimulus feature means that in both species the greater-duration stimulus, on average, tends to be judged as higher in intensity. That is, longer feels stronger.
To assess the effect of vibration intensity on perceived duration in more detail, the upper panel of Fig 4A shows results from the duration delayed comparison sessions, with separate psychometric curves plotted for each value of NID. In both species, there was a pronounced Interacting stimulus features in delayed comparison. A) Psychometric curves for 10 humans (left) and 7 rats (right) in the duration (top) and intensity (bottom) delayed comparison tasks. B) Upper plot: Bias caused by the nonrelevant stimulus feature, intensity, in duration comparison. Lower plot: Bias caused by the non-relevant stimulus feature, duration, in intensity comparison. In all plots, dots represent single subjects, bars represent mean of biases across subjects, while error bars represent the standard error of the mean across all subjects. The median value of each bias was significantly different from zero (humans: p = 0.002 for both intensity and duration bias, rats: p = 0.0156 for intensity bias, p = 0.032 for duration bias, Wilcoxon signed rank test). https://doi.org/10.1371/journal.pcbi.1008668.g004

PLOS COMPUTATIONAL BIOLOGY
shift of the duration psychometric curves as NID grows from negative to positive (dark red to yellow), signifying that a greater value of I2 relative to I1 increased the likelihood of the subject judging T2 > T1. The lower panel reveals the effect of vibration duration on perceived intensity. The substantial shift of the psychometric curves as NTD grows from negative to positive (dark blue to cyan) signifies that, in both species, a greater value of T2 relative to T1 increased the likelihood of the subject judging I2 > I1. We quantified the bias in perception as the slope of the linear correlation relating the change in the non-relevant feature to the change in the NTD or NID value at which the subject judged Stimulus 1 and Stimulus 2 as equivalent (the point of subjective equality, PSE). In humans and rats, both in the duration task ( Fig 4B, upper panel) and in the intensity task ( Fig 4B, lower panel), the median value of bias was significantly above zero (humans: p = 0.002 for both intensity and duration bias, rats: p = 0.0156 for intensity bias, p = 0.032 for bias, Wilcoxon signed rank test, S1 Fig). Median values reveal a greater influence of the irrelevant feature in humans than in rats.
Two observations point to the shifts of psychometric curves as a perceptual rather than a decisional phenomenon. First, variations in the non-relevant feature affected choices in the same way whether applied to Stimulus 1 or Stimulus 2 (S2 Fig) even though Stimulus 1 is dissociated from any decisional process; the choice can be generated only after presentation of Stimulus 2. Second, the biases were better explained as a horizontal psychometric curve shift than a vertical shift (S3 Fig). In short, the main finding of Figs 2 and 3 is that humans and rats readily extracted the perceptual feature required by the task, be it duration or intensity, but were biased by the non-relevant feature (intensity and duration, respectively).

Experiment 2: Interaction of vibration intensity and duration in a direct estimation task
Delayed comparison (Experiment 1) involves a number of steps: Stimulus 1 encoding, storage in working memory, Stimulus 2 encoding, and comparison of current Stimulus 2 to the memory of Stimulus 1. Interaction between intensity and duration could occur during Stimulus 1 storage, recall, or during the Stimulus 2 versus 1 comparison. In Experiment 2, perception in human subjects was measured through direct intensity estimation and direct duration estimation ( Fig 5A). If an intensity/duration confound were to persist in direct estimation, it would strengthen the argument that mixing emerges within the initial percept, before the percept is committed to or retrieved from working memory.
On a given trial, the subject received a single vibration, defined (as before) by I and T. A slider image appeared on the monitor 500 ms after the end of the vibration. By choosing the mouse-click position along the slider, the subject reported the perceived intensity of the vibration or else the perceived duration of the vibration. Scale normalization procedures are detailed in Methods. The slider did not display any landmarks, numbers, or ticks. To help subjects create two separate subjective scales, the orientation of the slider was specific to the task (e.g. horizontal for the intensity session and vertical for duration session), with task/orientation association randomized among subjects. The test stimulus set was comprised of 10 durations (linearly spaced from 80 to 800 ms in 80 ms steps) and 10 intensity values (linearly spaced from 9.6 mm/s to 67.2 mm/s in 6.4 mm/s steps). All 100 possible combinations of intensity and duration were presented in each session ( Fig 5B). Fig 5C, left panel, shows the results of direct duration estimation, averaged across subjects. Perceived duration increased not only with T, as expected, but also with I (designated by colors from dark red to yellow), confirming the main result from Experiment 1. From the same data set, Fig 5C, middle panel, plots the mean reported duration (pooling all presented durations) as a function of I (in log-scale), highlighting the intensity-induced bias. The rightmost panel of Fig 5C shows that the intensity bias, calculated as the linear coefficient of the fit between mean perceived duration and I in log scale, is significantly different from zero (Wilcoxon signed rank test, p = 0.002). Fig 5D, left panel, shows the results of direct intensity estimation, averaged across subjects. Perceived intensity increased not only with I, as expected, but also with T (from dark blue to cyan). Thus, longer stimuli were perceived as stronger in direct estimation, as found earlier in delayed comparison. Again, the mean reported intensity increases with the non-relevant feature T (Fig 5D, middle panel); the duration bias ( Fig 5D, right panel), calculated as the linear coefficient of the fit between mean perceived intensity and T in log scale, is significantly different from zero (Wilcoxon signed rank test, p = 0.002).
The interaction between the intensity and duration implies the existence of perceptual "invariants" in humans, and presumably also in rats. Considering the set of slider curves in the left plot of Fig 5C, a horizontal line at perceived duration 6 would intersect the yellow-to-dark red curves at a sequence of actual durations ranging from about 400 to 600 ms. Once adjusted for intensity, all such stimuli would be perceived as having the same duration. Notice that the spacing in time would decrease as intensity increases, highlighting that the iso-duration percept curves would not be linear: the effect of intensity on perceived duration saturates at the upper end of the intensity scale. Similarly, the slider curves in the left plot of Fig 5D imply intensity invariants. In rats, the existence of "iso-intensity percept curves" have been previously posited in rats and humans based on the influence of vibration duration on perceived intensity [7,9].

Experiment 3: A sensory integration account for time and intensity perception
Next, we asked which computations the brain might use to construct the two percepts from a common stream of sensory input. Neurons in rat somatosensory cortex precisely encode instantaneous speed [7], and analogous coding mechanisms exist in primates [10]. Because the vibration was stochastic, no instantaneous value could provide an intensity estimate for the PLOS COMPUTATIONAL BIOLOGY entire vibration [11]. A subject could, in theory, achieve optimal performance in the intensity task by linearly integrating (summating) the output of I-coding neurons over the entire vibration and then normalizing the integrated value by elapsed time. The denominator of this normalizing operation-elapsed time-could itself be the basis of the estimate of stimulus duration. But this computation would not explain the observed perceptual confound between the relevant and the irrelevant stimulus features, inasmuch as the time counter in the normalization is conceived of as an independent term.
As an alternative, we posit that the brain constructs the percept of both stimulus duration and stimulus intensity by nonlinear accumulation of the sensory signal over time. Nonlinear relations between stimulus features I, T and percepts are hinted at by the psychophysical functions of Experiment 2, (left panels of Fig 5C and 5D; also see S4 Fig).
A renowned model of accumulation in perceptual decision-making is the leaky integrator, in which some form of input is summated across time, while the accumulator simultaneously diminishes by some proportion of its accumulated quantity [12]. Leaky integration of sensory input can be formulated as: where γ is the percept, f(I t ,t) is the external drive to the integrator, λ is the leak rate and C/λ = τ specifies a time constant of integration. We now link the model to brain activity, hypothesizing that the leaky integrator's source of drive may be sensory cortex. To test the feasibility of this hypothesis, we ask whether neuronal firing from rat vibrissal somatosensory cortex (vS1) can be inserted into Eq (1) in place of the term f(I t ,t), to generate neurometric functions consistent with the observed rat psychometric functions (see Methods, S5 and S6 Figs). If so, then the parameters of integration that optimize the similarity between neurometric and psychometric functions would be informative about the underlying brain mechanisms. Fig  Here, a small rise in firing rate at vibration onset is visible, but subsequent firing rate was not significantly correlated with I.
Could inputs like those shown in Fig 6B be accumulated to generate the neuronal substrates of the intensity and duration percepts? We implemented two leaky integrator models as an attempt to account for duration and intensity perception. The two integrators differ at three levels: the leak time constant τ, the proportion of I-coding versus non I-coding neurons that are integrated, and the neuronal noise, which quantifies the variability in the firing pattern within the neuronal population that serves as the input to the integrator. To replicate the duration delayed comparison performance of one example rat, the input to the duration leaky integrator consists of 34% I-coding neurons (66% non I-coding neurons), accumulated with a time constant τ of 666 ms and a noise level of +/-3.1 standard deviations (see Methods for the details on the quantification of noise). Fig 6C shows how the accumulated quantity, γ duration , of the duration leaky integrator grows over time for different vibration I values. The accumulated quantity builds up almost linearly, due to the long time constant, and increases modestly with stimulus I, due to the drive from I-coding neurons. The left and right panels of Fig 6D show the strong similarity between the psychometric curves and the simulated neurometric curves, respectively, obtained for this example rat (see S5 and S6 Figs and Methods for the model fitting procedure).
To replicate the intensity delayed comparison performance of one example rat, the input to the intensity leaky integrator consists of 90% I-coding neurons (10% non I-coding neurons), accumulated with a time constant τ of 90 ms and a noise level of 1.6 standard deviations. As seen in Fig 6E, the intensity leaky integrator, γ intensity , saturates earlier than the duration leaky integrator, due to the short time constant. Moreover, the integrator grows more strongly with I than does the duration leaky integrator, a consequence of the predominant input from I-coding neurons. The left and right panels of Fig 6F show strong similarity between the psychometric curves and the simulated neurometric curves, respectively, obtained for this example rat. Psychometric/neurometric comparisons for all rats are given in S5 and S6 Figs. Fig 6G shows the optimal values of the 3 variables of the leaky integrator model obtained for each individual duration rat (in blue) and intensity rat (in red). The optimal parameters for the two tasks were clearly separated in two clusters in all 3 dimensions. To replicate the behavioral results of the duration task, the sensory drive must be integrated with a long time constant (average τ, 592 ms; range 282-794 ms) compared to that of the intensity task (average τ, 83 ms; range 74-90 ms). If the duration percept, γ duration , were to integrate the sensory drive with the time constant suitable for the intensity integrator, it would vary too greatly in relation to I, giving an unrealistically large intensity bias (S7 Fig, left plot). Symmetrically, if the intensity percept, γ intensity , were to integrate sensory drive with the time constant suitable for the duration integrator, it would vary too greatly in relation to T, giving an unrealistically large duration bias (S7 Fig, right plot).

PLOS COMPUTATIONAL BIOLOGY
Moreover, the two integrators differ in the nature of their sensory drive. Duration delayed comparison neurometric curves replicate actual psychometric curves only when the sensory drive includes a low proportion of I-coding vS1 neurons with high neuronal noise, whereas intensity neurometric curves replicate actual psychometric curves only when the sensory drive includes a high proportion of I-coding vS1 neurons with low neuronal noise (Fig 6G). The duration leaky integrator's tolerance for non I-coding neurons and for noise implies that it is "open" to inputs unrelated to the vibration sensory code. This is consistent with the fact that time perception is a supramodal process; in the perceptual experience of an event, all sensory channels are connected with one unique feeling of time [13]. Furthermore, multimodal (audio-visual) stimuli are perceived as longer in duration than unimodal stimuli, suggesting the convergence of separate streams [14]. One possible interpretation of our data is that the duration leaky integrator normally accumulates neuronal activity from sensory areas besides vS1, reflected in the integrator's requirement to receive non I-coding activity with high noise within its sensory drive. On this basis, we predicted that the percept of duration (but not intensity) could be affected by input carried within a sensory modality other than that whose duration must be judged. Psychophysical experiments in human subjects support the prediction (see S8 Fig and S1 Text), revealing that an acoustic stimulus delivered together with tactile vibration is accumulated by the duration integrator but not the intensity integrator.

Experiment 4: Dual leaky integrators
Having seen that the leaky integrator framework offers a plausible account for both percepts, when acting on sensory cortical spike trains, we next asked whether a single integrator is at work, one which operates with parameters tuned on each trial according to the current task (Model 1: Fig 7A). Alternatively, do two integrators, characterized by different parameter values, operate in parallel, with the subject retrieving the percept from the task-specified integrator (Model 2: Fig 7A)? To select between the two models, we designed Experiment 4, in which human subjects performed direct stimulus estimation, however the instruction cue indicating which stimulus feature to report was given to the subject either before or after stimulus delivery ( Fig 7B). According to the single integrator Model 1 of Fig 7A, performance would be good on trials with cue delivery prior to vibration onset, inasmuch as the integration Test of single versus dual integrators. A) Alternative hypotheses for the leaky integration process underlying the construction of both intensity and duration perception. Model 1 is represented by a single integrator that receives tactile drive but switches between task-specific values for the parameter τ. Model 2 is represented by dual integrators that receive the same tactile drive. Each integrator has task-specific values for the parameter τ. B) Schematic representation of cue-before versus cue-after direct estimation experiment. On half the trials, the cue providing trial instruction (symbolized by red or blue box) was provided before the vibration (above dashed line), and on the remaining half, the cue was presented after the vibration (below dashed line). C) Comparison of median perceived duration (upper row) and median perceived intensity (lower row) when the cue was presented before (left column) versus after (right column) the vibration, for 8 human subjects. Time of cue did not affect estimation.
https://doi.org/10.1371/journal.pcbi.1008668.g007 PLOS COMPUTATIONAL BIOLOGY parameters could be correctly "pretuned" for the relevant feature. Performance would be lower with cue delivery after the conclusion of the stimulus inasmuch as the integration parameters could not be optimally pretuned prior to stimulus presentation. According to the dual integrator Model 2 of Fig 7A, performance would be nearly the same on trials with cue delivery after vibration offset versus cue delivery prior to vibration onset. This is because the two integrators operate in parallel, each with optimized parameters. The results, illustrated in Fig 7C,

Rat tactile perceptual capacities
This study extends the range of high-level perceptual functions for which rats can serve as models. In Experiment 1, rats were able to perform a tactile intensity delayed comparison task with capacities similar to those of human subjects, replicating recent studies [6,7]. In Experiment 1, rats were also able to perform a tactile duration delayed comparison task, again with capacities strikingly similar to those of human subjects (Fig 3). Extracting stimulus duration and committing it to memory for future reference now enters for the first time the rodent perceptual repertoire.

Intensity/duration confound reflects intertwined mechanisms
Having established rats as an appropriate target for investigation, the study turned to the problem of understanding how, from a single tactile stimulus, rats and humans can extract multiple percepts-the vibration's intensity together with the time it occupies. In both tasks, judgments were biased by the non-relevant stimulus feature-duration in the intensity delayed comparison task and I in the duration delayed comparison task (Fig 4). As neither rats nor humans were able to generate independent representations of intensity and duration, a systematic perceptual interaction that transcends the species and the physical configuration of the receptors (fingertip and whiskers) must be at play. Experiment 2 replicated the perceptual interaction using a direct estimation task in human subjects, indicating that the interaction takes place regardless of whether the percept is manipulated in working memory (Fig 5).
Our results derive from the case in which the duration that must be measured is the elapsed time of the sensory stimulus itself. In conditions where the time to be judged lies between two discrete events-a start and stop signal, for example-the integration mechanisms remain to be worked out and the involvement of the dorsal striatum has been highlighted [15,16].
In earlier psychophysical work, an increase in perceived intensity with increased stimulus duration was reported in touch [7,17], audition [18][19][20] and vision [18]. The duration effect on perceived intensity was thought to arise from a temporal integration process: the sensory system summates input over time linearly, following Bloch's Law [21], or else nonlinearly [7,20,22]. Independently, an increase in perceived duration with non-temporal stimulus features has been reported in touch [23,24], audition [25] and vision [24,26]. However, the effect of stimulus features other than duration itself on perceived duration has been interpreted in a separate framework, which posits that an internal, central clock keeps track of time [1]. In this framework, the increase of perceived duration with increasing vibration I would be an attentional phenomenon, where a stronger stimulus leads to an increase in arousal, resulting in augmented speed of the central clock [27,28].
By uncovering both perceptual confounds on the basis of a single set of stimuli (Figs 1D and 2), our experiments made it possible to configure perceived intensity and perceived duration within a unified framework. Opposing the view of two independent confounds, our experiments suggest that the influence of stimulus duration on perceived intensity and the influence of vibration I on perceived duration constitute inextricable phenomena.

Integration of sensory drive and alternative models
We propose leaky integration of one common source of sensory drive (Eq 1) as the mechanism underlying both sets of psychophysical data. In Experiment 3, we recorded vibration-evoked neuronal activity from vibrissal somatosensory cortex (vS1) of awake, behaving rats and inserted that firing as the sensory drive term of the leaky integrator formulation. After setting the parameters of sensory drive integration separately for the two tasks and then simulating choices by applying a decision rule to the leaky integrator's accumulated quantity, the resulting neurometric curves mimicked observed intensity and duration psychophysical curves (Fig 6).
The three parameter settings that optimized neurometric/psychometric overlap are informative about mechanisms. First, the leak time constant, τ, was greater for the duration integrator (282-794 ms) than the intensity integrator (74-90 ms). These values fit with the intuition that the accumulated quantity underlying the perception of the passage of time must grow in a manner approaching linearity, while the accumulated quantity underlying the perception of vibration intensity must quickly reach a steady state. It is interesting that the intensity time constant is similar to that for the accumulation across touches of vibrissal kinematics-evoked activity in the formation of a texture percept [29,30], suggesting some common integrative mechanism or circuit underlying both percepts.
Second, the degree to which stimulus-unrelated neuronal activity is integrated differed for the two integrators. Intensity neurometric curves fully overlapped psychometric curves only when the leaky integrator input was comprised of 80-100% I-coding neurons, suggesting that the neuronal population underlying the final intensity percept gives more weight to synapses transmitting stimulus-specific information and tends to exclude non-responsive neurons. Whereas the tactile features of a multimodal event may be experienced distinctly from the concurrent visual or acoustic features, there is not a distinct sense of elapsed time for each sensory stream [14]. At the level of neuronal processing, this implies the convergence of diverse channels onto a single duration integrator, promiscuous to all inputs. Our study uncovered three correlates of this presumed convergence. First, duration neurometric curves matched psychometric curves only when the proportion of I-coding vS1 neurons in the sensory drive was reduced to 14-68%. Second, neurometric matched psychometric only when the noise level within the sensory drive was boosted from 1.2-1.7 (for the intensity neurometric) to 2.6-3.9 standard deviations. Third, the duration percept but not the intensity percept of tactile vibrations was influenced by the addition of acoustic noise alongside the trial (S8 Fig), again consistent with the expected multimodal convergence upon the duration integrator.
The hypothesis that perceived duration is achieved through leaky integration of stimulusrelated and stimulus-unrelated sensory input is particularly relevant to the debate on models for the perception of time in the scale that ranges from tens of milliseconds up to a few seconds. Our results are in line with other work that assumes temporal integration processes behind the encoding of the passage of time [31,32] and are also in agreement with the idea that sensory-specific areas contribute to time perception [3,33]. They are harder to reconcile with the amodal central clock theory [1]. The "state-dependent networks" model [34] and the "striatal beat frequency" model [35] have not yet been challenged as to the dependence of time perception on non-temporal stimulus features. In the first type of model, time is encoded by the evolving temporal pattern of activity of a recurrent neuronal network, so that almost any network could in principle represent the elapsed time without the need of an explicit representation of duration [34,36]. In the second type of model, time is encoded by a pattern of relative phases of different oscillators, thought to be present in both thalamic and cortical neuronal activity, and is read out by striatal neurons, which act as coincidence detectors. Both models would explain the influence of vibration I on its perceived duration as an increase or decrease of the speed by which the activity of the connected neuronal population evolves in time or else by which the oscillators follow the pattern of their relative phases. Whether intensity-dependent modulation could be reliably implemented in these models is not known.
Beyond the integrative mechanisms acting on the stimulus within the single trial, another factor that can cause nonlinearity in the percept or choice is the recent history of stimuli [37,38].

Integration timescales point to target regions and mechanisms
Temporal leaky integration of sensory information can be performed by a recurrent neuronal circuit [39][40][41][42]. Recurrent neuronal circuits have been widely used to explain decision-related neuronal activity in areas of the brain such as lateral intraparietal cortex (LIP) [43], premotor cortex [44], prefrontal cortex [45], and dorsal striatum [46,47]. The different level of leakiness for intensity versus duration could be achieved by a difference in the strength of recurrent connections in the network, and also on the different levels of background input [39,48,49]. Experiment 4, where human performance did not depend on preparatory knowledge about the percept to be reported, argues that the two percepts are generated through two computations that operate in parallel. We speculate that two different neuronal populations receive a common input, and the final state of these populations can be interrogated separately in order to produce the required judgment (Fig 6A, Model 2). The hierarchical ordering in the timescales of intrinsic fluctuations in the primate cortex, increasing from posterior to anterior [50,51], suggests that duration may be computed in a more anterior region with respect to intensity, perhaps at a downstream stage. Our rat experiments were limited to individuals reporting one percept, either duration or intensity. We speculate, however, that rats also simultaneously experience two distinct percepts, an assumption that would require discrete brain populations to simultaneously integrate sensory drive with different time constants. Primary somatosensory cortex is characterized by a short intrinsic timescale, and does not show temporal integration [7,11,29,50,52]. Intrinsic timescales for intensity integration may be found in primary and secondary vibrissal motor cortex (vM1 and vM2) and intrinsic timescales for duration integration may be found in farther anterior or medial regions. The convergence of tactile and acoustic input for the generation of the duration percept might occur in a population that processes sensory inputs from several modalities, such as premotor cortex [53] or dorsal striatum [54].
The two integrators are not a literal portrayal of a physiological mechanism, but are a characterization of some network property that participates in the conversion of the vibration sensory code to the conscious experiences of intensity and duration. The dual integrators constitute a unified framework inasmuch as both key features of the stimulus-the coding of vibration I and elapsed time occupied by that stimulus-contribute to both integrators. While the framework put forward in this study cannot exclude the feasibility of all other models, it does create a set of predictions that can serve to alert us as to which network properties should be sought in future physiological work. For instance, the successful generation of neurometric curves to replicate both duration and intensity perception suggests vS1 as a common input, a hypothesis that could be directly tested by optogenetic control over vS1. Our model also makes the straightforward prediction that the neuronal population implementing the readout of stimulus duration must be modulated by stimulus I.

Ethics statement
Human subjects were tested after giving their written informed consent. Protocols conformed to international norms and were approved by the Ethics Committee of SISSA (protocol number 6948-II/7). 14 male Wistar rats (Harlan, San Pietro al Natisone, Italy) were housed individually or with one cage mate and maintained on a 14/10 light/dark cycle. Daily access to water was restricted to promote motivation in the behavioral task, yet weight gain followed a standard Wistar-specific curve, indicating that the quantity of water obtained during training and testing was comparable to the ad lib quantity. After each session, rats were placed for several hours in a large, multistory enriched environment with other rats. Protocols conformed to international norms and were approved by the Ethics Committee of SISSA and the Italian Health Ministry (license numbers 569/2015-PR and 570/2015-PR).

Human and rat subjects
Thirteen healthy human subjects (age range 22-35 yrs) were tested in the delayed comparison task. Only subjects that reached better than 70% performance in both intensity and duration delayed comparison (10 out of 13) were included in the analysis. The same 10 subjects were then recruited for the direct estimation tasks. All subjects were recruited on-line through the SISSA Sona System. (https://sissa-cns.sona-systems.com/). The number of subjects was 8 in Experiment 3 and 9 in Experiment 4. Individual rats were trained on a single delayed comparison task: 7 were intensity rats and 7 duration rats.

Stimulus generation
Vibrations were generated by stringing together sequential velocity values (v t ) at 10,000 samples/s, taken from a normal distribution. Converting v t to its absolute value, sp t , the distribution takes the form of a folded, half-Gaussian (see Fig 1B). A Butterworth filter with 150 Hz cutoff was then applied to yield low-pass filtered noise. The sp t time series for a given trial was taken randomly from among 50 unique sequences of pseudo-random values. Because stimuli were built by sampling a normal distribution, the statistical properties of an individual vibration did not perfectly replicate those of the distribution from which it was constructed. As a vibration's actual mean speed could not precisely match the distribution from which it was sampled, the assigned value may be considered "nominal" mean speed, referred to as I.

Delayed comparison task for rats
Each trial began when the rat positioned its nose in the nose-poke (equipped with optic sensor) and placed its whiskers on the plate (Fig 1A). After a delay of 500 ms, Stimulus 1 was presented, characterized by intensity, I1, and duration, T1 (Fig 1B). After the inter-stimulus delay of 2 s, Stimulus 2 (with I2 and T2) was presented (Fig 1C). The rat remained in the nose-poke throughout both stimuli and could withdraw only when the "go" cue sounded at the end of the post-stimulus delay of 500 ms. Early withdrawal was considered an aborted trial and went unrewarded; it was not scored as correct or incorrect. After the go cue, the rat selected the left or right spout; reward location depended on the relative values of I1 and I2 for rats doing the intensity delayed comparison task, while it depended on the relative values of T1 and T2 for rats doing the duration delayed comparison task. Incorrect choices went unrewarded. Trials with I1 = I2 or T1 = T2, according to the task, were rewarded randomly.
Rats learned the delayed comparison task by generalizing the comparison rule across the entire stimulus range, referred to as the stimulus generalization matrix (Fig 2), whereby neither stimulus alone provided the information necessary for a correct choice (1,2). Seven rats were trained to discriminate T1 vs T2 and another 7 were trained to discriminate I1 vs I2. The stimulus range used during both duration and intensity delayed comparison task was from 112 to 1000 ms for stimulus duration, and from 25 to 160 mm/s for stimulus I.

Delayed comparison task for human subjects
Experiments 1 and 4 employed a delayed comparison design. Stimulus 1 was characterized by intensity I1 and duration T1. After the interstimulus delay of 1 s, Stimulus 2 (with I2 and T2) was presented. Stimuli delivered to human subjects on the fingertip were the same as those used in rats except that the velocity values were halved. Each subject went through both an intensity and a duration delayed comparison session on different days. Subjects were verbally instructed to report which of the two stimuli was perceived as longer in duration or stronger in intensity, according to the behavioral task, by pressing the left (for Stimulus 1) or right (for Stimulus 2) arrow on the computer keyboard. They received visual feedback (correct/incorrect) on each trial through a monitor. A total of 1,456 trials were presented at each session.
The stimulus range used for intensity delayed comparison session was from 161 to 557 ms for stimulus duration, and from 9.28 to 110.36 mm/s for stimulus intensity. The stimulus range used for duration delayed comparison session was from 87 to 1034 ms for stimulus duration, and from 17.2 to 60 mm/s for stimulus intensity.

Direct estimation task
The same human subjects that went through the duration and intensity delayed comparison task participated in the estimation task. Each subject went through both a duration estimation and an intensity estimation session, held on different days. Each session began with a training phase. In this phase, subjects received 40 stimuli, sampled randomly from the 100 possible stimuli (10 possible I values from 9.6 mm/s to 67.2 mm/s and 10 possible T values from 80 to 800 ms, linearly spaced), to become confident with the task and to sample the stimulus range. In the test phase, a single stimulus was presented, characterized by intensity, I, and duration, T. After a post-stimulus delay of 500 ms, a slider appeared on the screen. The slider did not present any landmark, ticks or numbers. The orientation of the slider was different between the two sessions, either vertical or horizontal. Subjects were instructed to report the perceived intensity or else the perceived duration of the vibration on a subjective scale, in which the extreme left/bottom position indicated a very weak or a very short stimulus, and the extreme right/top position indicated a very strong or a very long stimulus. The report was done by mouse-clicking on the chosen position along the slider. A total of 1,300 trials was presented at each session.
Five durations linearly spaced from 80 to 800 ms, and 5 I values from 9.6 to 67.2 mm/s were used. A visual cue, either a blue or a red square, was presented for 1 second, either before or after the delivery of the vibration. The orientation of the slider was kept horizontal for both intensity and duration estimation trials, so that the orientation could not be used as a cue for the trial type and subjects were forced to attend to the visual cue. A total of 1,300 trials was presented at each session.

Analysis of human and rat delayed comparison data
To characterize the performance of the intensity delayed comparison task, we computed the proportion of trials in which subjects judged Stimulus 2 greater than Stimulus 1 on stimulus pairs characterized by a fixed I1 (I1 = 32 mm/s for human subjects, I1 = 64 mm/s for rats) and different I2 values, separately for each normalized time difference (NTD) value, defined as (T2 -T1) / (T2 + T1). We fit the data with a four-parameter logistic function using the nonlinear least-squares fit in MATLAB (MathWorks, Natick, MA): where NID is normalized intensity difference, (I2-I1) / (I2+I1), δ is the lower asymptote, λ is the upper asymptote, 1/ν is the maximum slope of the curve and μ is the NID at the curve's inflection point.
For the duration delayed comparison task we computed the proportion of trials in which subjects judged Stimulus 2 > Stimulus 1 on stimulus pairs characterized by a fixed T1 (T1 = 300 ms for human subjects, T1 = 300 ms for rats) and different T2 values, separately for each NTD value by fitting: where δ is the lower asymptote, λ is the upper asymptote, 1/ν is the maximum slope of the curve and μ is the NTD at the curve's inflection point.
In order to quantifying the bias of stimulus duration on the percept of stimulus intensity in the intensity delayed comparison task, we then computed a linear correlation between the PSE values fitted for different NTD values, and the actual NTD values. The additive inverse of the regression coefficient, was defined as duration bias. Symmetrically for the duration delayed comparison task, we computed a linear correlation between the PSE values fitted for different NID values, and the actual NID values. The additive inverse of the regression coefficient, was defined as intensity bias.

Delayed comparison task: perceptual versus choice bias
Fig 4A reveals that both species show a pronounced shift in their psychometric curves on both duration and intensity discrimination tasks due to the task irrelevant feature, NID and NTD respectively. Horizontal and vertical shifts of the psychometric curve are frequently attributed to different steps in the cognitive process, a perceptual shift and a choice shift, respectively [55]. To quantify whether the shifts are better explained as purely horizontal or vertical we fit the following two models of decision probabilities: where P (Stim2 > Stim1) is the probability of choosing Stimulus 2 greater than Stimulus 1, Δp and Δb are the normalized differences of the relevant and irrelevant task features, respectively (i.e. NTD and NID for the duration delayed comparison task, and NTD and NID for the intensity delayed comparison task). p L is the probability of a lapse trial, β p is the linear weight that the task relevant feature has on choice, β b is the linear weight that the task irrelevant feature has on choice, and s(x) is a logistic function, ε p is the constant perceptual bias, and ε L is the lapse trial constant choice bias. Both models in Eqs 4 and 5 assume that on each trial there is a probability p L that the choice will not be determined by the task-relevant features. A non-sensory error is called a lapse. On the remaining trials, choice is determined by a generalized linear model (GLM), specifically the logit link function, of the task relevant feature Δp. The two models differ in the role of the task irrelevant feature. In the model of Eq 9, Δb is linearly combined with Δp on the non lapsed trials, whilst in the model of Eq 10, Δb goes through a separate GLM that determines the choice on lapsed trials. This means that the model of Eq 9 assumes that the task-irrelevant feature biases the effective percept yielding only a horizontal shift, whilst the model of Eq 5 assumes that the task-irrelevant feature biases choice on lapsed trials, yielding only a vertical shift. We fit both models to the human and rat data using a Python software package PyMC3 [56]. In both cases, we pooled all the subjects together. We then compared the goodness of fit of both models using the Widely Applicable Information Criteria (WAIC, [57]). WAIC shows that the entire dataset is better explained by a perceptual bias (Eq 4) as compared to a choice bias (Eq 5). In the case of the humans, the normalized WAIC difference (ndWAIC which is equal to WAIC divided by its standard deviation) was equal to 6.62. In the case of the rats, ndWAIC = 0.34 (S3 Fig). A detailed comparison can be found in Table 1.

Analysis of direct estimation task data
Results from the direct estimation task were analyzed taking into account two main factors. First is that the subjects tended to be more variable in their responses at the beginning of the session before setting a consistent subjective scale. In order to keep this variability from affecting our results, we calculated the mean response for each of the possible combinations of I and T and excluded outlier responses, those more than 1.5 SD displaced from the mean. Second was that not all subjects used the whole range of the slider; every participant set their minimum and maximum responses at a different position in the scale. In order to make each subject's subjective scale comparable, we used a min-max normalization algorithm: Normalized x n ¼ 9 x n À minðxÞ maxðxÞ À minðxÞ þ 1 ð6Þ where x n is the non-normalized response on trial n, x is the range of total responses, and Normalized x n is the normalized response. We then multiplied by 9 and added 1, so that the normalized responses range from 1 to 10.
In order to estimate duration bias for intensity estimation trials in Experiment 3, we first computed the average normalized response across all possible I, for each T. We then computed a linear regression between the average normalized responses and stimulus T, and defined the regression slope as duration bias.

PLOS COMPUTATIONAL BIOLOGY
Similarly, in order to estimate intensity bias for intensity estimation trials in Experiment 3, we first computed the average normalized response across all possible T, for each I. We then computed a linear regression between the average normalized responses and stimulus I, and defined the regression slope as intensity bias.

Electrode implantation and data acquisition
Trained rats (n = 2) were anesthetized with 2%-2.5% Isoflurane delivered with oxygen under controlled pressure through a plastic snout mask. They received an implant in the primary vibrissal sensory cortex (vS1), which was accessed by craniotomy, using standard stereotaxic technique (centered 2.8 mm posterior to bregma and 5.8 mm lateral to the midline). The multielectrode array (Tucker Davis Technologies) was inserted via micromanipulator. The extracellular activity of vS1 was sorted into 31 single units and 92 multiunits, as verified through the spike waveform and the refractory period observed in interspike interval histogram using a MATLAB-based software, UltraMegaSort 2000 [58,59]. In total, 123 neurons were recorded in 7 recording sessions. Part of the data set analyzed in the present study has been analyzed to different purposes in previous work [6].

Neuronal population response
All analyses and statistical tests were done with custom codes written in MATLAB. We propose that the population activity of the vS1 that represents vibration I is integrated by downstream areas to produce the percepts of intensity and duration. An observer is imagined to make a perceptual choice based on the neuronal population integrated value, γ. The resulting choice, when plotted against stimulus difference (NID or NTD), is referred to as the neurometric curve. A high degree of similarity between the neurometric curve and the psychometric curve observed in the behavioral task supports the feasibility of the posited neuronal operation as a mechanism underlying behavior. In order to obtain an adequate number of trials per neuron per condition, we focused on trials with fixed and long duration and variable I values. In these sessions, stimuli lasted 200 ms, 600 ms, and in some sessions 1s. Only trials in which the stimulus duration was 600 ms or longer were used in the generation of neurometric curves.
In order to construct the neurometric curves from equivalent (I, T) values to those used in the psychometric analysis, pseudo random trials were created by measuring the responses of individual neurons from t = 0 to T (1 ms bin size) and averaged over equi-I trials. In order to have the same number of trials per combination of I and T, 200 resampled pseudo random trials per (I, T) stimulus class were generated by resampling [60]. The data set was then divided in half with the first half serving as the training set to estimate parameters of interest, and the second half serving to test how accurate the model with selected parameters predicts the observed psychometric curve.

Leaky integration of neuronal firing
To quantify the coding properties of individual neurons, the Spearman correlation strength between stimulus I and the average response during the first 600 ms of the vibration was measured. Neurons with significant correlation (p < 0.05) were considered I-coding neurons (66 out of 123, 54%). The integrated neuronal response, γ, for each (I, T) stimulus class was calculated in 1 ms steps and γ at time T (ms) was taken as input to the observer's choice. Neuronal activity drove the leaky integrator through a modified form of Eq (1): The external drive f(r t , n t ) to the integrator is: where the first sum is the summation of the response of N neurons randomly sampled from the I-coding neuronal population. The second sum is the response of non I-coding neurons, and n t is the neuronal fluctuation, sampled from a Gaussian distribution with 50 different values of variance and mean. The neuronal fluctuation within the sensory drive is generated from the Gaussian distribution with zero mean and standard deviation σ, Nð0; s 2 t Þ. The variance of the noise term at a given time-point is proportional to the average response of input channels pooled over all I-coding vs. non I-coding neurons at the same time-point (R t ).
The proportion of coding vs. noncoding neurons was chosen among 50 different (N, M) combinations such that N NþM spanned 0 to 1 linearly. The leak time constant τ was chosen randomly among 50 values that spanned 50-800 ms linearly.

Neurometric curves
For each combination of (τ, proportion coding neurons, noise level) parameters the neurometric curve is calculated as: where δ and λ are the probability of upper and lower lapse (incorrect decision unrelated to perceptual processing) and are estimated from the psychometric curve of each subject. Integrator outputs for Stimulus 1 and Stimulus 2 (γ 2 and γ 1 ) are compared on a trial by trial basis (pseudorandom trial) for each (I, T) stimulus pair. The neurometric curves should replicate two key features of the psychometric curves. The first is the bias, namely, the shift in the psychometric curve caused by the feature that should be exclude. This is quantified by the slope of the linear correlation relating the change in the non-relevant feature to the change in the NTD or NID value at which the subject judged Stimulus 1 and Stimulus 2 as equivalent-that is, the point of subjective equality (PSE). In the case of the neuronal analysis, we found that the noise in the bias measure (caused by neuronal response variability) could be diminished by computing the area under the neuromeric curve rather than the shift in PSE. The second key feature is the overall performance achieved by an observer who compares two stimuli on the basis of the integrated neuronal firing. We selected the values of leaky integrator parameters (τ, proportion coding neurons, noise level) such that the bias and performance resulting from the output of Eq (11) fell within 5 percent of the actual psychometric values that were meant to be replicated (S5 and S6A Figs). The neurometric curves constructed with the selected values of leaky integrator parameters produce similar performance and bias to those observed in behavior. For each selected value (S5 and S6A Figs, right, filled yellow dots), the neurometric curves were constructed using the test trials. The neurometric values for each (I, T) stimulus pair was then compared to the behavioral response. The parameter set that produces the least error was then used as the optimized integrator values (Figs 6D and 6G and S5 and S6B) Supporting information S1 Fig. Quantification of the perceptual biases. Upper panel: Duration delayed comparison task. Each bar is the standard error, centered on the mean, of the psychometric curve PSE for each NID value across all 10 human subjects (solid) and all 7 rats (dashed), relative to the PSE for the NID = 0 condition. Lower panel: Same analysis for intensity delayed comparison task. The downward slanting distribution of data indicates that, for the duration task, PSE shifted to the left as NID grew while, for the intensity task, PSE shifted to the left as NTD grew. neurons that yield leaky integrator performance replicating the actual performance of the rat (black dots). Lower panel shows the values of τ and percent I-coding neurons yield leaky integrator bias replicating the intensity bias of the same rat (black dots). Middle panel shows the parameters that gave a match in performance (blue dots) and intensity bias (orange dots) in the 3d parameter space. Yellow dots indicate the parameters values that produced a match in both features. Among those parameter values, the ones that minimized the difference between the choice of the rat and the choice of the ideal observer, for all NTD and NID values, were used to generate the neurometric curves (rightmost panel). B) Psychometric curves (upper row) and neurometric curves (middle row) for all individual rats. (TIFF) A) Schematic representation of a trial with non-informative acoustic noise delivered through headphones. B) Bars denote mean performance on the duration and intensity tasks for trials with noise on and off, across 9 human subjects. Each dot represents a single subject's mean performance. Orange error bars are standard error of the mean across subjects. The presence or absence of noise did not affect accuracy (Kruskal-Wallis test, p = 0.72, Bayes Factor = 3.07 for the duration task, p = 0.66, Bayes Factor = 2.02 for the intensity task). C) Effect of acoustic noise on the bias caused by the task-irrelevant feature. For the duration task, the presence of noise reduced the bias normally caused by intensity (one sample, one-tailed Wilcoxon signed rank test, p = 0.0273). For the intensity task, noise did not affect the bias caused by duration (one sample, one-tailed Wilcoxon signed rank test, p = 0.5). Each dot represents a single subject's bias difference, whilst the bar represents the average across subjects. Error bars are standard error of the mean across subjects. (TIFF)

S9 Fig. Observed bias in Experiment 4.
Left plot: Bias caused by the irrelevant feature (I) in the duration estimation task. Each dot corresponds to a single subject, while the bar is the mean across subjects. Right plot: Bias caused by the irrelevant feature (T) in the intensity estimation task. Each dot corresponds to a single subject; the bar is the mean across all 8 subjects.
In both plots, orange error bars are standard error of the mean across subjects. Cue delivery did not affect the bias of the non-relevant feature (Kruskal Wallis test: for duration estimation p = 0.83, Bayes Factor = 2.92; for intensity estimation p = 0.75, Bayes Factor = 2.88). (TIFF) S1 Text. Discrimination task with supplementary auditory noise. (PDF)