A neurocomputational model for optimal temporal processing

Abstract

Humans can estimate the duration of intervals of time, and psychophysical experiments show that these estimations are subject to timing errors. According to standard theories of timing, these errors increase linearly with the interval to be estimated (Weber’s law), and both at longer and shorter intervals, deviations from linearity are reported. This is not easily reconciled with the accumulation of neuronal noise, which would only lead to an increase with the square root of the interval. Here, we offer a neuronal model which explains the form of the error function as a result of a constrained optimization process. The model consists of a number of synfire chains with different transmission times, which project onto a set of readout neurons. We show that an increase in the transmission time corresponds to a superlinear increase of the timing errors. Under the assumption of a fixed chain length, the experimentally observed error function emerges from optimal selection of chains for each given interval. Furthermore, we show how this optimal selection could be implemented by competitive spike-timing dependent plasticity in the connections from the chains to the readout network, and discuss implications of our model on selective temporal learning and possible neural architectures of interval timing.

1 Introduction

Our world changes in time, and our brain faces the challenge to cope with these changes. Sequences of stimuli often convey information in their order and timing, e.g. in speech or music. Our sense of causality requires knowledge about the natural temporal order in which events happen. Our brain can use this knowledge together with information about the typical duration of events to predict the evolution of sequences of events. Also on the level of behavior, timing is obviously crucial, as a given action can be right or wrong based only on the time of its execution.

Recognizing this importance, many researchers have posed the question of how time is represented in the brain. For some time-related stimulus features such as the speed of motion in the visual field or the frequency of a tone, such representations could be found in distinct brain areas such as the middle temporal region of visual cortex for speed and the inner hair cells for tone frequency. For the duration of a stimulus or the interstimulus interval, however, no single neural correlate has been identified (Buonomano and Karmarkar 2002). Instead, lesion and imaging studies have revealed the possible participation of a whole network of structures such as the cerebellum, the basal ganglia, the thalamus and various cortical regions such as prefrontal cortex (PFC) and the supplementary motor area (Buonomano and Karmarkar 2002; Lewis and Miall 2003; Ivry and Spencer 2004). This spectrum of brain structures is accompanied by a variety of possible timing mechanisms. Basically, any process in the brain could be used to represent the time that has elapsed while the process unfolds. Models of temporal processing have exploited neural structures that range from single neurons (Tieu et al. 1999; Grossberg and Schmajuk 1989), neural oscillators read out by coincidence detectors (Matell and Meck 2004) and short-term synaptic plasticity (Buonomano 2000) to reverberating loops within the cerebellum (Medina et al. 2000; Yamazaki and Tanaka 2006), slowly climbing activity in PFC neurons during working memory tasks (Durstewitz 2006) and stochastic decay of memory traces (Kitano et al. 2003).

These neurobiological models focus on the ability of a given neuronal circuit to represent temporal information. However, most of them are not sufficient to explain behavioral performance in timing tasks, as they do not discuss errors in the representation. Temporal precision and the change of timing performance under various conditions are one major subject of psychological experiments in both animal and man (Grondin 2001; Gibbon et al. 1997). Results from such experiments can constrain neuronal models regarding their predictions on timing errors. A typical class of experiments is given by the task of interval reproduction (Allan and Kristofferson 1974). The participants are presented with an interval of duration T, e.g. a continuous tone or a series of flashes of light. Afterwards, they are asked to reproduce this interval for example by pressing a button as long as they think the interval was. This experiment results in a set of reproduced intervals, usually clustered around some mean m with a standard deviation σ. These two measures are used to characterize the participant’s ability to reproduce the interval T. To explain the results from such experiments, information-processing models are used. They are composed from a set of functional processing stages which interact with each other and provide some understanding of the mechanisms behind the experimental results. The most popular of these models is the pacemaker-accumulator system (Creelman 1962): An oscillating or random process (the pacemaker) generates pulses with a fixed frequency, while another system (the accumulator) counts the number of these pulses. The number of pulses accumulated during an interval is used to estimate the duration of the interval. This theory has been formalized in several ways (Gibbon 1977; Killeen and Fetterman 1988) and can explain a wide range of phenomena (Grondin 2001; Gibbon et al. 1997). Another widely used concept is called interval timers (Ivry 1996) or labeled lines (Buonomano and Karmarkar 2002). Models within this framework assume a number of units which operate with different time constants, such that each of these units is tuned to a specific interval of time.

While information-processing models are designed to explain timing errors in behavioral experiments, they are only rarely connected to some neural substrate and thus can not identify the source of these errors. To connect the experimental results with a proposed timing mechanism in the brain, it is useful to assume that the brain performs an estimation of the time elapsed since the press of a button, such that the participant releases the button when the memorized duration of the target interval is reached. The estimation will be subject to timing errors σ and possibly also some bias m − T if the participant responds systematically too early or too late. Taking this view, a neuronal model of temporal processing is supported by psychophysical experiments if its estimation statistics are compatible with the response statistics found in the experiment.

One prominent finding of timing experiments is that timing errors increase monotonically with the duration of the interval to be processed. This increase is characterized by the Weber fraction, defined by σ/ m. The shape of this fraction as a function of the interval T is then to be reproduced by a neuronal model. According to “scalar expectancy theory” (SET, Gibbon 1977), temporal perception obeys Weber’s law, which means that the Weber fraction does not depend on T. Indeed, many studies were conducted that seem to confirm this “scalar property” (Gibbon et al. 1997). However, a model that tries to explain a constant Weber fraction faces an inherent problem: If the timing errors arise from noise affecting the timing units without correlations, or with a finite correlation length, these timing errors will increase as $\sqrt{T}$, as the variances σ ² add up linearly over time. Thus, the observed linear increase of the timing error needs another error source. SET, and also some other information-processing theories solve this problem by assuming that the observed scaling of the errors is built in one of the processing stages such as the counter (Killeen and Fetterman 1988) or a memory stage (Gibbon 1977). Some of the neuronal models also use ad hoc assumptions about scalar variability, e.g. in the distribution of synaptic weights (Matell and Meck 2004) or in the rise times of firing rates (Grossberg and Schmajuk 1989), to generate the scalar property. Others assume special properties of the noise, such as a low-pass filtered frequency band (Tieu et al. 1999) or independence of the stimulus duration (Staddon and Higa 1999). There are also a few attempts to explain Weber’s law by inherent properties of the model (Okamoto and Fukai 2001). However, none of these studies discusses the relation of the steeply increasing errors with the principle of optimality: If it is possible to represent time with an error of the order of the square root of the interval duration, why should a brain under evolutionary pressure use some mechanism that is worse than that? Moreover, there is also evidence that the scalar property is not universally valid in temporal processing. The Weber fraction increases both at short and long intervals, with a minimum in between (Bizo et al. 2006; Drake and Botte 1993; Getty 1976). While the increase at short intervals can be explained by an additional, time-independent error (“generalized Weber’s law”), the increase at longer intervals was mostly ignored in information-processing models (but see Staddon and Higa 1999), because it is not compatible with the predictions from SET. To date, there is no neurocomputational model that accounts for the entire Weber fraction with its decreasing, constant, and also increasing parts.

In this paper, we offer a model in which the U-shaped form of the Weber function emerges directly as a result of an optimization process. The model consists of a number of timing units with different time constants. These units project onto a set of readout neurons, which shows a unique spike pattern for each interval to be represented. In this framework, synaptic plasticity implements an optimal selection of timing units under limited resources. For the units itself, we demand high temporal precision and reliability to obtain optimal results, and also a sufficiently wide range of time constants. In a noisy brain, arguably the most precise temporal code is provided by synfire chains (Abeles 1991). A synfire chain is a feed-forward network with strongly converging and diverging connectivity. Such a network is able to stably propagate a wave of neural activity from pool to pool at a precision in the range of milliseconds (Herrmann et al. 1995; Diesmann et al. 1999) even in the presence of usual biological inaccuracies. This property makes synfire chains exquisite timing devices: If an activity volley is injected at stimulus onset, the time elapsed after the onset is reliably converted into the position of the volley in the chain. In this way, synfire chains can integrate two of the main concepts of timing: If each of the pools is used, the chain is equivalent to a delay line and can be compared to a pacemaker-accumulator system, in the sense that each pool corresponds to a pulse that is accumulated and the transmission speed from one pool to the next corresponds to the pulse frequency. On the other hand, it is also possible to use several chains as interval timers, if they have different transmission speeds and only the last pool of each chain is used for time estimation. And finally, it is also possible to connect the final pool of a chain with the first one, which results in a neural oscillator with a frequency given by the number of pools and the transmission speed.

For a single synfire chain, timing errors again increase like the square root of the interval length. However, we arrive at the result that the timing errors increase superlinearly with the delay of transmission from one pool to next. As the number of pools in a chain can be assumed to be limited, this constitutes an additional error source for longer intervals, as they can only be represented by chains with a larger delay. Under these conditions, we show that the observed U-shaped form of the Weber fraction arises from optimal selection of the chain with the lowest possible timing error for any given interval to be represented. Finally, we propose a combination of spike-timing dependent plasticity [STDP (Bi and Poo 1998)] and homeostatic plasticity (Turrigiano 1999) for the connections from the chains to the readout neurons implements an optimal and unique selection of chains. This selection is based on the fact that the effective learning rate of STDP depends on the temporal distribution of the input pattern.

The paper is organized as follows: In Section 2, we describe the model structure and provide the equations that are used to simulate the neurons, synapses and plasticity rules. Section 3 defines the notions of transmission delay and timing error and demonstrates how the temporal statistics of synfire chains affect temporal processing. Furthermore, we discuss the effect of variations in the model parameters, especially the rise time of the postsynaptic potentials (PSPs) on the delay and the timing error. In Section 4, we combine these results to a theory of optimal temporal processing and show how the U-shaped form of the Weber fraction emerges. Section 5 discusses the implementation of the optimal and unique selection of chains by synaptic plasticity. Finally, in Section 6 the results are discussed regarding their implications for selective learning, effects of attention, and also possible extensions of the model.

2 Neuron and network model

2.1 Network structure

The model consists of neurons which are described by their membrane potential V _i , and connected by synapses of strength w _ij , where i denotes the postsynaptic and j the presynaptic neuron. The neurons are organized in different networks (see Fig. 1, left): Synfire chains consist of L pools denoted by ${\cal P}_l$ which contain N neurons each. Each neuron in a pool ${\cal P}_l$ has a probability of p _S to be connected to any neuron in the subsequent pool P _l + 1 with strength w _S

$$p(w_{ij}) = \left\{ \begin{array}{rl} p_S & \qquad \mathrm{for} \; w_{ij}=w_S \quad \\ 1-p_S & \qquad \mathrm{for} \; w_{ij}=0 \end{array} \right. \forall \; i \in {\cal P}_{l+1}, j \in {\cal P}_l.$$

(1)

If all neurons in pool ${\cal P}_l$ fire nearly synchronously with a small temporal jitter, this induces on average N w _S inputs in each neuron in the subsequent pool ${\cal P}_{l+1}$. Thus, the firing times from the preceding pool are averaged and the jitter is reduced in the firing times of pool ${\cal P}_{l+1}$. As each neuron in pool ${\cal P}_{l+1}$ in turn projects on average to N w _S neurons, the activity in pool ${\cal P}_{l+2}$ will be even more synchronized. If all neurons in the chain are disturbed by synaptic noise, the temporal jitter in the spike times will not decrease to zero, but converge to a near-synchronized fixed point where the effect of the connectivity and the noise are balanced (Herrmann et al. 1995; Diesmann et al. 1999, see Fig. 1, right).

Fig. 1

Left: Illustration of the model structure. A readout network ${\cal M}$ receives convergent connections to from different synfire chains such as ${\cal C}_1$ and ${\cal C}_2$. By the competition between the respective weights, w ₁ and w ₂, the network determines which chain optimally responds at a time interval represented by the output unit in ${\cal M}$. Right: Raster plot showing the spikes in the readout network ${\cal M}$ and selected pools from the chains ${\cal C}_1$ and ${\cal C}_2$. Each dot corresponds to a spike. In ${\cal C}_1$, activity propagates faster and with smaller jitter σ _P compared to ${\cal C}_2$

Apart from the synfire chains, there is a readout network $\cal M$ consisting of M neurons with no connections among each other ($w_{ij} = 0 \quad \forall i,j \in \cal M$), but which connections from the synfire chains. A pool P _l is connected to a readout neuron i ∈ M by the rule

$$p(w_{ij}) = \left\{ \begin{array}{rl} p_M & \qquad \mathrm{for} \; w_{ij}=w_S \quad \\ 1-p_M & \qquad \mathrm{for} \; w_{ij}=0 \end{array} \right. \forall \; j \in {\cal P}_l.$$

(2)

The set of all neurons in a given synfire chain is denoted by ${\cal C}_\alpha$, as all the parameters are identical across chains except for α, which is defined below. The values of all parameters regarding network connectivity are listed in the left column of Table 1.

Table 1 Values of all model parameters that are used unless otherwise stated

2.2 Neuron model and synapses

The neurons are modeled as leaky integrate-and-fire units embedded in a stochastic background network. While the membrane potential V _i stays below a threshold V _thr, the time evolution of V _i is given by

$$\tau \frac{dV_i}{dt} = - (V_i-V_\mathrm{rest}) + \sum_{k : w_{ik} \ne 0} \frac{w_{ik}}{\alpha} \sum_{j=1}^{S_k} \mathrm{PSP}\left(t-t_{kj}^\mathrm{sp}\right) + I_\mathrm{noise} + I_\mathrm{ext}.$$

(3)

Without any input (I _noise = I _ext = 0, S _k  = 0 ∀ k), V _i relaxes exponentially to the resting potential V _rest with time constant τ. I _noise represents synaptic input from the stochastic background network, which is collapsed into an excitatory and an inhibitory part

$$I_\mathrm{noise} = \epsilon_+ {\cal N}(\lambda_+, \; \lambda_+) - \epsilon_- {\cal N}(\lambda_-, \; \lambda_+),$$

(4)

where ε ₊, ε − > 0 and ${\cal N}(m, \, \sigma^2)$ denotes a random variable with a Gaussian distribution with mean m and standard deviation σ. In this form, the Gaussians approximate two Poisson processes with rates λ ₊ and λ −, respectively. Using the standard parameter values (see Table 1, middle column), ignoring the threshold V _thr and without further input, the membrane potential converges to a mean of $\langle V \rangle = -46.4$ mV and a standard deviation of σ _V  = 1.4 mV.

Whenever the membrane potential V _i crosses the threshold V _thr from below, the neuron i fires a spike. V _i is then set to the reset potential V _reset and the current time t is included in the set of firing times of neuron i by increasing the number of spikes S _i by one (S _i →S _i  + 1) and setting $t^\mathrm{sp}_{i\,S_i} = t$. A spike in neuron j influences the membrane potential V _i of all neurons $i: w_{ij} \ne 0$ it is connected to presynaptically. The time evolution of the induced PSP in neuron i is described by an α function

$$\mathrm{PSP}(t) = \frac{t}{\alpha} \exp \left( -\frac{t}{\alpha} \right),$$

(5)

where α is the rise time of the PSP. The synaptic weights w _ij in Eq. (3) are normalized by α to ensure that the total impact of a single spike on the postsynaptic membrane potential does not change with α. As mentioned before, different synfire chains denoted by ${\cal C}_\alpha$ will differ only in this parameter. No additional synaptic delays are incorporated, so α is the only parameter that determines the time course of the PSP. Introducing a distribution of delays does not qualitatively change the results.

2.3 Synaptic plasticity

The connections from the chains to the readout neurons $\{w_{ij}: i \in {\cal M}\}$ are subject to two forms of synaptic plasticity: STDP (Bi and Poo 1998) and homeostatic plasticity (Turrigiano 1999). STDP is applied after each spike $t^\mathrm{sp}_{ik}$ in neuron i. For all spikes $t^\mathrm{sp}_{jl}$ in neurons that are connected to neuron i $\big\{t^\mathrm{sp}_{jl}: w_{ij} \ne 0 \lor w_{ji} \ne 0 \big\}$, the time $t = t^\mathrm{sp}_{jk} - t^\mathrm{sp}_{il}$ is calculated and the respective weights are updated by adding

$$\Delta w = \left\{ \begin{array}{rl} A_p \exp(-t/\tau_p) \;\; & \mathrm{if} \;\; t > 0 \\ -A_d \exp(t/\tau_d) \;\; & \mathrm{if} \;\; t < 0 \end{array} \right.$$

(6)

Thus, the synapse between two neurons is strengthened if the presynaptic spike occurs earlier than the postsynaptic spike (t > 0), but if this order is reversed (t < 0), the synapse is weakened. This introduces the notion of causality into the learning rule.

Synaptic weights under control of the STDP learning rule tend to diverge either to zero or infinity, depending on the network’s activity (van Rossum et al. 2000). One way to prevent this is to complement STDP by a homeostatic learning rule, which adjusts the synaptic weights such that the network achieves a certain mean firing rate a _goal. A simple implementation of this mechanism is given by a proportional-integral feedback controller (van Rossum et al. 2000)

$$\frac{dw}{dt} = g_P w (a_{\mathrm{goal}} - a) + g_I w \int_0^t{dt (a_{\mathrm{goal}} - a)}.$$

(7)

This equation is applied to update w after each trial based on the actual mean firing rate a of all readout neurons. This seems more appropriate then using it every time step, as homeostatic plasticity is believed to act slowly (Turrigiano 1999). Parameter value for both forms of plasticity are listed in Table 1, right column.

3 Temporal statistics of synfire chains

3.1 Quasi-spatial representation of time

An interval can be represented by a synfire chain if the stimulus onset triggers a volley of activity in the first pool ${\cal P}_1$ of the chain. This volley can be characterized by the number of spikes a(1) in the first pool and the temporal jitter σ _P (1) of these spikes around the mean t(1), which is defined to be zero. We assume that the volley follows a Gaussian distribution ${\cal N}(t(i), \, \sigma_P(i)^2)$, which turns out to be a good approximation for the simulated data. For a large range of initial conditions (a(1), σ _P (1)), activity will then stably propagate through the pools, approaching a fixed point (a, σ _P ) in both parameters (Diesmann et al. 1999). The time at which the activity wave has reached pool ${\cal P}_i$ can be estimated by the mean spike times of the neurons in that pool

$$t(i) = \frac{1}{a(i)} \, \sum_{n=1}^N \sum_{j=1}^{S_n} t_{nj}^{\mathrm{sp}}, \qquad \mathrm{with} \; a(i) = \big| \big\{ t^\mathrm{sp}_{nj}: \, n \in {\cal P}_i \big\} \big|.$$

(8)

a(i) is the number of spikes in pool ${\cal P}i$, which may also be larger than N if there are neurons that spike more than once. For practical calculation in simulations, only those spikes enter this equation which give rise to a firing rate above a threshold value defined from the background firing rate. The rate is calculated by counting the number of spikes in a time bin of 0.1 ms.

From t(i), we define the transmission delay Δt(i) from pool ${\cal P}_i$ to pool ${\cal P}_{i+1}$ as

$$\Delta t(i) = t(i+1) - t(i).$$

(9)

An interval can only be reliably estimated if Δt(i) has a robust relation to i. For a synfire chain, one would expect that Δt(i) converges to a fixed point value Δt, as a(i) and σ _P (i) do. In Fig. 2 (left panel), Δt(i) is shown for initial conditions (a(1), σ _P (1)) = (100, 4 ms), which is close to the border of the basin of attraction of the fixed point (a, σ _P ). As expected, the mean of Δt(i) converges to a fixed point after a short transient, establishing an approximately linear transformation of time into the position of the volley in the chain. This representation does not need to be literally spatial, as the pools are defined only by the topology of their connections and not their spatial position.

Fig. 2

Left: Transmission delay Δt(i) as a function of the pool number i. After a transient, the mean of Δt(i) converges to a constant Δt, while some finite jitter σ _Δt remains. Right: Timing error σ _T as a function of the stimulus duration T for PSP rise times α = 0.5 ms (lower curve) and α = 1.5 ms (upper curve). The dots represent the simulation data and the line is a plot of Eq. (14) with σ _Δt fitted to the data. The coefficient σ _Δt is fitted to 0.035067 (±28) for α = 0.5 and 0.13093 (±15) for α = 1.5

The exact form of the transient depends on the initial conditions. For instance, starting close to the border of the basin of attraction that surrounds the fixed point results in an overshoot of Δt(i) in the first pools (Fig. 2, left). In the following, we suppress the effect of the transients by choosing initial conditions that are close to the fixed point.

3.2 Timing errors

While the mean of Δt(i) converges to a constant Δt, fluctuations in the actual realizations of Δt(i) remain. They can be considered Gaussian random variables

$$\Delta t(i) = {\cal N}\big(\Delta t, \sigma_{\Delta t}^2\big)$$

(10)

with mean and standard deviation independent of i. The temporal jitter σ _Δt can be derived directly from the parameters a and σ _P for the steady state. Adding two random variables with a Gaussian distribution results in a Gaussian variable where mean and variance are the sum of those from the two original variables. ΔT is the sum of two random variables t(i + 1) and − t(i) which are themselves the sum of a variables $t_{nj}^{\mathrm{sp}}$ each, divided by a. The spike times, in turn, have a Gaussian distribution with standard deviation σ _P (cf. Section 3.1). As a and σ _P are constant in the steady state, both t(i + 1) and − t(i) have a standard deviation of $\sigma_P / \sqrt{a}$, and σ _Δt , so that we find

$$\sigma_{\Delta t}^2 = \frac{2}{a} \; \sigma_P^2 = \mathrm{const}.$$

(11)

The estimate of an interval T can now be written as

$$T(i) = \sum_{l=1}^{i}{\Delta t(l)} = \Delta t \cdot i + {\cal N}\big(0, \, \sigma_T^2\big)$$

(12)

The second equality results from Eq. (10). σ _T is the timing error which appears in psychophysical experiments. Again, the variances of Δt(i) add up and the timing error results in

$$\sigma_T(i)^2 = \sum_{l=1}^{i}{\sigma_{\Delta t}^{2}(l)} = \sigma_{\Delta t}^2 \cdot i,$$

(13)

where the second equality holds at the fixed point value of σ _Δt [Eq. (11)]. Using Eq. (12), one sees that the timing error increases with $\sqrt{T}$

$$\sigma_T(T) = \sigma_{\Delta t} \cdot \sqrt{T}.$$

(14)

Note that the timing error σ _T and the spike time jitter σ _P refer to quite different concepts. σ _P denotes the standard deviation of the spike times within a pool and also within the same trial. This jitter must not increase with time, or else activity can not be stably propagated. On the other hand, σ _T is the accumulated error of the estimating the transmission delay Δt across trials. Increase of σ _T means that the temporal information encoded in the synfire chain is gradually lost at longer times, although it may still stably propagate activity in each individual trial.

In Fig. 2 (right panel), we have plotted the estimate of σ _T for two different synfire chains ${\cal C}_{0.5}$ and ${\cal C}_{1.5}$, based on 1000 simulations each. One sees that both curves can be fitted with a $\sqrt{T}$-law, with 94.7% of the total variance explained for ${\cal C}_{0.5}$, and 95.2% for ${\cal C}_{1.5}$.

Linear representation of time and a $\sqrt{T}$ law for timing errors are also key properties of a pacemaker-accumulator system (Gibbon 1977), so in this sense a synfire chain is equivalent to this class of models.

3.3 Parameter variations

A synfire chain can be seen as a timing device with a time constant of L Δt. In this section, we examine how we can increase Δt to obtain different time scales without changing the pool number L.

3.3.1 Effect of varying α

The position of the fixed point in a and σ _P , and thus, ΔT and σ _ΔT depends strongly on the parameters of the synfire chain. The easiest way to change ΔT is to vary the rise time α of the PSP. If the time until each spike fully affects the postsynaptic membrane potential increases, the transmission time of the entire volley will also be delayed. Moreover, this parameter has already been studied in its influence on the fixed point in a and σ _P (Diesmann 2002).

We assess the effect of α by running simulations with 20 trials each, raising α by 0.05 ms after each set of trials. The dependency of Δt, a and σ _P on α turns out to be well fitted by a polynomial of second order

$$f(x) = A x^2 + B x + C.$$

(15)

As seen from the coefficients in Table 2, Δt and σ _P increase mostly linear with α. As A ≪ B, the quadratic term only becomes relevant as a small correction at higher α. a, on the other hand, decreases quadratically over the whole range of α, but only moderately in total (A,B ≪ C). For σ _Δt , a fit to Eq. (15) only works well for small values of α. Furthermore, when we check whether Eq. (11) holds for the simulated data, it turns out that it does only for α ≤ 2 (Fig. 3, left panel). As we are interested in the relationship between σ _Δt and Δt, we fit a curve to the data points in these two dimensions. As a boundary condition, we demand σΔt(0) = 0, as it makes no sense to assume a timing error if activity travels through the chain infinitely fast. With this constraint, the data turns out the be fitted very well with a third-order polynomial

$$\sigma_{\Delta t}(\Delta t) = A \Delta t^3 + B \Delta t^2 + C \Delta t. $$

(16)

The coefficients are listed in the final row of Table 2. Once again, the dependency is mostly linear, as A,B ≪ C, with a nonlinear correction is of the order of Δt ³. The fit explains 82.4% of the variance in the data (Fig. 3, right panel).

Table 2 Coefficients of Eq. (15) fitted to the data sets of a, Δt and σ _P (first three rows), and coefficients of Eq. (16) fitted to the data set of σ _Δt (final row)

Fig. 3

Left: The variance in the transmission delay σ _Δt as a function of the PSP rise time α. Right: σ _Δt as a function of mean transmission delay Δt. The dots are data from the simulation. In the left panel, the line is a plot of Eq. (11) with the values of a and σ _P from the simulation. In the right panel, the line is a fit to Eq. (16) using the parameters in Table 2, final row

What is the reason for the nonlinearities in Δt, σ _P and σ _Δt that become relevant at higher values of α? Note that α has an upper limit at 2.7 ms. At this point, a bifurcation occurs, i.e. the fixed point becomes unstable or even collides with the saddle point at the border of the basin of attraction, making them both vanish (see Diesmann 2002 for a discussion of these scenarios). We propose that the nonlinearities occur as α approaches the bifurcation point. This also explains why the effect of α is stronger for σ _Δt compared to the other parameters: Close to a bifurcation point, the transients that lead to the fixed point become longer and more variable. That does not affect the jitter of individual pools, nor the mean transmission delay, but it makes the time that the volley spends in the first few pools highly unpredictable, which increases the total timing error σ _T , and thus, indirectly also σ _Δt .

3.3.2 Effect of other parameters

While the focus of this study is on the influence of α, here we briefly discuss how other model parameters affect Δt and σ _Δt . Although σ _Δt increases with Δt as α is increased, it is conceivable that these two measures are anticorrelated as another parameter is changed. We checked whether this is possible, with all those parameters that directly affect the dynamics of the network. Individual parameters to be changed are, apart from α, the synaptic weights w _S , the connection probability p _S , the number of neurons in a pool N and the membrane time constant τ. While changing p _S and N, w _S is normalized to 0.345/p _S and 100/N, respectively. Without this normalization, p _S and N would have similar effects as w _S , as the total number of presynaptic synapses of a neuron is changed.

Furthermore, the statistics of the membrane potential is important for the dynamics, described by $\langle V \rangle$ and σ _V . $\langle V \rangle$ enters the dynamics only by its distance from the firing threshold $V_\mathrm{thr} - \langle V \rangle$. This distance is most easily changed by modifying V _thr. σ _V , on the other hand, can only be modified by jointly changing λ ₊ and λ − such that $\langle V \rangle$ stays constant. This is guaranteed if there is a certain linear relationship between the two rates (Diesmann 2002).

We increased and decreased each parameter individually until the chain either breaks down (synchronizing effect too weak) or activity volleys form spontaneously without external stimulation (synchronizing effect too strong). Then, we calculated Δt and σ _Δt at parameter values slightly before one of the two events occur. Figure 4 shows the results for all cases where Δt increased. In all these cases, σ _Δt increases as well. We also included α in the analysis for comparison. From Fig. 3, as well as from Eq. (16), the coefficients in Table 2 and the upper limit of α at 2.9 ms, one can see that changing α increases Δt up to 5.44 ms and σ _Δt up to 0.44 ms. Thus, changing α makes it possible to increase Δt much more than any of the other parameters, and, as seen from Fig. 4, also at the lowest relative error.

Fig. 4

Effect of the variation of the model parameters on the chain transmission characterized by Δt and σ _Δt . The common starting point of each line corresponds to the combination of standard parameters listed in Table 1. Only one single parameter is varied at a time. Each line is drawn up to the point where the synfire chain becomes unstable. The only exception is the line for α, which actually extends up to Δt = 5.44 ms and σ _Δt  = 0.44 ms

Note that this analysis does not exclude the possibility of increasing Δt without an increase of σ _Δt . For instance, one could increase α and also increase the number of neurons in a pool N to compensate for the increase in σ _Δt . However, such a compensation would always result in a decrease of Δt as well, limiting its possible range. Furthermore, Fig. 4 illustrates that the effect of any other individual parameter then α is rather limited. We conclude that one can only attenuate the increase of σ _Δt , but not abolish it completely over the full range of Δt.

Nevertheless, different combinations of parameters may still extent our results. For instance, it was shown in Wennekers and Palm (1996) that Δt could be increases by a factor of 2 to 5 by changing the membrane time constant τ and the external input I _ext (and thus, $V_\mathrm{thr} - \langle V \rangle$), compared to an increase of merely 25% which we report here. What we have shown is that α is the most efficient parameter in changing Δt in the sense that it induced the largest dynamic range for Δt with the least relative increase in σ _Δt .

4 Optimal temporal processing

If there are different synfire chains ${\cal C}_\alpha$ with different transmission delays Δt, a given interval T could be potentially represented by each of these chains. This redundancy opens the possibility to optimize the representation of time, in the sense of minimizing the timing error σ _T for each interval T. The principle of optimization has already been used before in a more abstract pacemaker-accumulator system (Killeen and Weiss 1987).

An interval T is represented by a pool i in the chain ${\cal C}_\alpha$ with transmission delay Δt. Only one of these parameters can be freely varied, the other one is fixed by the implicit equation T = T(i, Δt). σ _T depends on both parameters, but much stronger on Δt ($\sigma_T = O(\Delta t^3)$, Eq. (16) and Eq. (13)) than on i ($\sigma_T = O(\sqrt{i})$, Eq. (13)). Thus, minimization of σ _T requires using the lowest possible transmission delay min(Δt) and represent all intervals by the different pools i of the fastest chain ${\cal C}_{\mathrm{min}(\alpha)}$.

This optimization results in a timing error that increases as $\sigma_T = O(\sqrt{T})$. This explains only part of the experimental results, namely the decreasing Weber fraction for short intervals. At longer times, optimization must be constrained such that a linear or faster-than-linear increase of σ _T results. We propose that this constraint is given by a limited chain length. If a chain has a maximum of L pools, the range of intervals that can be represented by the fastest chain ${\cal C}_{\mathrm{min}(\alpha)}$ has an upper limit of T = min(Δt) ·L. For longer intervals, a chain with a higher Δt and thus, a higher σ _T must be used. As the timing errors are dominated by the third-order dependency on Δt, the constraint optimization problem can be formulated as

$$\Delta t^*(T) = \{\min(\Delta t) \; | \; T \le \Delta t \cdot L\}.$$

(17)

where Δt ^*(T) is the optimal choice of Δt for representing the interval T. In the simplest case, assuming a smooth distribution of Δt, every interval can be exactly represented by an individual synfire chain, and Eq. (17) simplifies to

$$\Delta t^*(T) = \frac{T}{L}.$$

(18)

Taken together, Eqs. (14), (16), and (18) yield the optimal form of the timing error

$$\sigma_T^*(T)\! =\! \left\{\! \begin{array}{l} \displaystyle \sigma_{\Delta t}(\min(\Delta t)) \cdot \sqrt{\!\frac{T}{\min(\Delta t)}}\;\; \mathrm{for} \; T\! \le \min(\Delta t) \cdot L \\ \displaystyle \frac{A}{L^{5/2}} T^3 + \frac{B}{L^{3/2}} T^2 + \frac{C}{\sqrt{L}} T \;\; \mathrm{otherwise} \\ \end{array} \right.$$

(19)

with the coefficients A, B, and C from Eq. (16) and Table 2 (final row). For T > min(Δt) ·L, the timing error is entirely dominated by σ _Δt . To minimize σ _Δt , always the final pool L of each chain is used [Eq. (18)], so the $\sqrt{T}$ dependency degenerates into $\sqrt{L}$, which is constant for all T and enters the coefficients of the third-order polynomial. $\sigma_T^*$ depends mostly linear on T, with small nonlinear corrections of the order of T ³, just as σ _Δt .

In Fig. 5 (left), we have plotted σ _T (T) for a number of different ${\cal C}_\alpha$ from simulations together with the optimal $\sigma_T^*(T)$ from Eq. (19), which is close to the lower envelope of the simulation data. Both from Fig. 5 and Eq. (19), one sees the three regimes that are found in the experiments: $\sigma_T \propto \sqrt{T}$ for T ≤ min(Δt) ·L, σ _T increasing faster than linear for T ≫ min(Δt) ·L, and an intermediate regime where the linear term in Eq. (19) dominates and σ _T / T is approximately constant. Consequently, the Weber fraction σ _T / T calculated from Eq. (19) follows the experimentally observed U-shaped form (Fig. 5, right).

Fig. 5

Left: Timing error, e.g. variance of the total runtime of a chain σ _T as a function of T for various values of α. The solid curves depict simulation data and the dotted line represents the optimal timing error $\sigma_T^*(T)$ from Eq. (19). It is close to the lower envelope of the simulation data. Right: Weber fraction σ _T / T as a function of T calculated from the lower envelope in Fig. 5 (left). The U-shaped form of the Weber fraction that is known from the psychophysical experiments is reproduced

5 Optimization by competitive STDP learning

In the preceding section, we have shown how the observed timing errors can be explained by an optimal selection of synfire chains. In this section, it remains to be shown how this selection is neuronally implemented and a unique representation of time is formed from the different chains. Here we show how both issues can be resolved if the chains project to a set of readout neurons and the synaptic weights of these projections are subject to both STDP and homeostatic plasticity.

To see how the representation of an interval T can be learned in this framework, consider the following experiment (cf. Fig. 8, top panel). At time t ₀, a stimulus S ₀ (called initiation stimulus) is given that activates the first pool of all chains and makes the volleys travel along the chains with their respective speeds. At a second time t ₁ = t ₀ + T, another stimulus S ₁ (training stimulus) activates the readout neurons. Around this time, inputs from the chains also arrive. If this experiment is repeated, the connections to those pools of the chains will grow that were active slightly before the stimulus S ₁. Eventually, a certain set of readout neurons will fire at time t ₁ even in the absence of stimulus S ₁ if the stimulus S ₀ is given at time t ₀.

Using this paradigm, we first show in Section 5.1 that the effective learning rate of STDP is higher for inputs from a fast synfire chain with a peaked spike time distribution compared to the input from a slower one, meaning that the mean synaptic weight of the connection is higher after a given number of trials. Second, we consider a set of readout neurons that are connected to two different chains, one with a low α (${\cal C}_{\alpha_1}$) and one with a higher one (${\cal C}_{\alpha_2}$). If the input from the two chains arrives at the same time, we can test the optimal selection of chains (Section 5.2): The timing errors σ _T in the readout neuron should be comparable to the smaller errors in chain ${\cal C}_{\alpha_1}$. If the inputs arrive at different times, we can test for the unique representation (Section 5.3): Even if a training stimulus is given at both arrival times, the readout neurons should learn to fire only once. Both of these properties are confirmed and can be extended to a scenario with several input chains.

5.1 Effective learning rate depending on timing errors

Before we explain the mechanism that brings about optimality and uniqueness of the representation, we elaborate on how the effective learning rate of STDP depends on the temporal distribution of the presynaptic spikes. For this purpose, we consider only connections from a single synfire chain and neglect homeostatic plasticity. Furthermore, we assume for simplicity that the postsynaptic spikes in the readout neurons (i.e. the training stimulus S ₁) occur at a fixed time t ₁ without any jitter. Assuming for example a Gaussian distribution of S ₁ around t ₁ does not change the following argument, as it only increases the variance of the relative time between the pre- and postsynaptic spikes.

To illustrate the effects of different temporal spreads in the presynaptic spike times, consider two Gaussian firing time distributions with the same mean spike time t ^*, but different variances (Fig. 6). The temporal asymmetry in STDP, and the exponential decrease of its efficiency with increasing |t| (cf. Eq. (6) and Fig. 6, top panel) are the reasons why the width of the input distribution influences the effective learning rate: The trial average of this rate is given by the convolution of the weight increase at a given value of t and the distribution of t. From Fig. 6, one can see that the area under a peaked distribution largely overlaps with the area under the positive branch of the learning curve up to its half width at half maximum, marked by the two dotted lines. The mean learning rate is therefore larger for a more strongly peaked input distribution than for a broader one, as the latter overlaps more with areas of less positive increase, and even with the negative branch of the STDP curve.

Fig. 6

Illustration of the STDP learning curve (top panel), aligned to spike time distributions from two different chains, one with a faster transmission (middle panel) and one with a slower transmission (bottom panel). t denotes the relative time between the presynpatic and postsynaptic spike (cf. Eq. 6). If the means are aligned, the different standard deviations cause a relative advantage for the narrower distribution. This is because the peaked distribution has a larger overlap within the area in between the dashed lines than the broader one

In Sections 3 and 4, we have shown that the temporal spread of the firing times in synfire chains increases with α. Thus, by the preceding argument, the mean weight from chain neurons to a readout neuron after a given number of learning trials should decrease with increasing α. To test this assumption, we fully connected all neurons in the pool number 50 of a chain to all the readout neurons and performed the simulated learning experiment 50 times for each different value of α. S ₁ was presented 3 ms after the mean firing time of the chain neurons in pool 50. Figure 7 shows that the synaptic weights after the learning trials indeed decrease for increasing α. This result is independent of the pool number and of the α used in the connections to the readout units.

Fig. 7

Mean synaptic weights $\langle w_M \rangle$ of the connections from a pool (shown is pool number 50) to a readout neuron after 50 presentations of the corresponding interval. The graph is obtained by varying the PSP rise time α used in the chain connections. All weights were initialized by the values 0.3

5.2 Optimal selection of synfire chains

If a readout neuron is connected to two synfire chains ${\cal C}_{\alpha_1}$ and ${\cal C}_{\alpha_2}$ (Fig. 1), its firing pattern may be shaped by the input from both of them. The combination of STDP and homeostatic plasticity introduces synaptic competition among the chains: The weights are increased by STDP with a different effective rate (cf. Section 5.1), but the rate of compensation by homeostasis is the same for both chains [cf. Eq. (7)]. Thus, the faster chain ${\cal C}_{\alpha_1}$ with the higher STDP rate will win the competition and dominate the firing pattern of the readout neuron.

First, we study the effect of this competition for the case that both inputs arrive at the same time. This means that the neurons which are connected to the readout neuron have the same mean firing time in both chains, marked by the dotted line in Fig. 8 (top panel). S ₁ is given slightly after this mean. If the input from chain ${\cal C}_{\alpha_1}$ dominates the firing pattern of the readout neuron, its timing error σ _T should be lower compared to the case where it only receives input from the slower chain ${\cal C}_{\alpha_2}$. This corresponds to a selection of the optimal chain ${\cal C}_{\alpha_1}$.

Fig. 8

Top: Illustration of the activation pattern in the readout neurons $\cal M$ that are driven by the synfire chain ${\cal C}_{\alpha_1}$ and ${\cal C}_{\alpha_2}$, cf. Fig. 1. At time t ₀, the initiation stimulus S ₀ activates the first pool of each chains and causes volleys along all the chains with their respective speeds. At time t ₁ = t ₀ + T, another stimulus S ₁ that indicates the termination of the temporal interval activates the readout neurons. Around this time, inputs from the chains also arrive. The learning rule (see Fig. 6) increases the connections to those pools that were active slightly before the stimulus S ₁. Middle: Mean synaptic weight $\langle w_M \rangle$ of all connections from ${\cal C}_{0.5}$ (black curve) and ${\cal C}_{1.5}$ (gray curve) to the readout neurons $\cal M$ as a function of the number of trials. During the first 200 trials, the weights are modified by STDP and homeostatic plasticity, while the final 100 trials are used to calculate $\sigma_T^M$. Ultimately, the connections to the faster chain ${\cal C}_{0.5}$ are much stronger than those to the slower chain ${\cal C}_{1.5}$. Bottom: Timing error $\sigma_T^M$ in the readout neurons $\cal M$ as a function of α ₂ (cf. Fig. 6). The grey curve shows $\sigma_T^M$ for α ₁ = 0.5 and the black curve for α ₁ = α ₂. For the case of α ₁ = 0.5, the timing error is dominated by the input from $C_{\alpha_1}$ and much lower compared to α ₁ = α ₂ for α ₂ > 0.5

To test if this selection takes place, we connected the readout neurons to one pool of chain ${\cal C}_{\alpha_1}$ and to another pool of ${\cal C}_{\alpha_2}$, chosen such that the mean firing times of the two pools are at 68 ms, which is the largest interval in the chain ${\cal C}_{0.3}$. α ₁ was fixed to 0.5 ms, while α ₂ was varied from 0.3 to 2.5 ms. If there is no pool which is activated exactly at 68 ms, we chose the one that is closest to this time and shift the starting time of the respective chain. This shift is always less than a millisecond. Without such a shift, the deviations of the means from 68 ms would be another source of timing errors. For every α ₂, we performed 10 sets of trials with 300 trials each. The first 200 trials were used to modify the synaptic weights to the readout neurons by STDP and homeostatic plasticity, while in the last 100 trials, the timing error $\sigma_T^M$ in the readout neurons was calculated without any further learning. $\sigma_T^M$ is defined in the same way as σ _T for the synfire chains (cf. Section 3), just using the neurons in the readout network instead of those in a pool. The target rate is set to a _goal = 2 spikes per neuron and trial, one spike from S ₁, and another from the chain input. The chain input usually arrives before stimulus S ₁, which is set to be strong enough to make the neurons fire even shortly after a spike induced by the chain input.

Figure 8 (middle panel) shows the synaptic competition between the two chains. Initially, STDP dominates the learning dynamics and increases the two types of connections according to the speed of the input chain (cf. Section 5.1). This produces an overshoot over the firing rate over the target rate a _goal, and homeostatic plasticity is finally strong enough to bring the weights down again. Different from STDP, the homeostatic learning rule is blind to the different speeds of the input chains, but reduces the weights only according to their current strength [cf. Eq. (7)]. Thus, the difference between the connections from $C_{\alpha_1}$ and ${\cal C}_{\alpha_2}$ remains as the mean firing rate approaches a _goal, resulting in partial suppression of the input from the slower chain ${\cal C}_{\alpha_2}$.

In Fig. 8 (bottom panel), we show the timing error $\sigma_T^M$ resulting from this synaptic competition as a function of α ₂. For comparison, we conducted another simulation where α ₁ = α ₂, instead of being fixed to 0.5 ms. In this case, the timing error is fully determined by the input from ${\cal C}_{\alpha_2}$. One sees that the error is much lower in the case of synaptic competition for all α > 0.5. For α = 0.3, on the other hand, the errors in both cases are comparable, as ${\cal C}_{\alpha_2}$ is the faster chain now and thus dominates $\sigma_T^M$.

5.3 Unique representation in readout neurons

We now study the second case that is possible with the connectivity depicted in Fig. 1, ${\cal C}_{\alpha_1}$ having an different mean firing time from $C_{\alpha_2}$. In this case, the question is how to prevent the readout neuron to respond to two different time intervals. If the system is exposed to two training stimuli S ₁ and S ₂ slightly after each of the chain firing means, both intervals could be trained to the same neuron (see Fig. 9, top panel). Thus, one could not tell which of the two times has elapsed upon firing of this neuron.

Fig. 9

Top: Illustration of the activation pattern in the readout neurons $\cal M$. Same as in Fig. 8 (top panel), except for the different mean firing times of ${\cal C}_{\alpha_1}$ and ${\cal C}_{\alpha_2}$ and two training stimuli occurring at t ₁ and t ₂, which are slightly after the temporal means of the two chains. Bottom: Mean firing time of the readout neurons T _M as a function of the mean firing time $T_{\alpha_2}$ in ${\cal C}_{\alpha_2}$. The firing time of ${\cal C}_{\alpha_1}$ is kept constant to $T_{\alpha_1} = 140.2$ ms. For $T_{\alpha_2} < T_{\alpha_1}$, the readout neurons fire at about the same time as the neurons in chain ${\cal C}_{\alpha_2}$

Such a double training of the same neuron is also prevented by synaptic competition. If the same neuron fires two times responding to a single training stimulus, the firing rate is higher than the target rate a _goal. This leads to a compensation that weakens the connections to both chains, but finally leads to a suppression of the input from the slower one, as its connection was weaker in the first place (cf. Section 5.1). The readout neurons only respond to the interval represented by the faster chain and the unique representation is restored.

As a test of this reasoning, we use the same simulation as in Section 5.2, with the only difference that the connecting pool of ${\cal C}_{\alpha_2}$ is chosen to match a certain mean firing time different from the one of ${\cal C}_{\alpha_1}$. This interval between the two mean firing times is now varied instead of the α, which are kept constant to α ₁ = 1.5 ms and α ₂ = 0.5 ms.

Figure 9 (bottom panel) shows the mean firing time of the readout neurons T _M as a function of the varied mean firing time $T_{\alpha_2}$ in $C_{\alpha_2}$. For a unique representation, T _M should either be identical to $T_{\alpha_2}$ (grey dotted diagonal) or to the fixed firing time of ${\cal C}_{\alpha_1}$ (black horizontal line). The figure shows that such a unique representation is maintained for all $T_{\alpha_2} < T_{\alpha_1}$. For $T_{\alpha_2}$ after the mean firing time of $C_{\alpha_1}$, the situation is less clearly described. While T _M faithfully follows $T_{\alpha_2}$ for 20 ms after $T_{\alpha_1}$, it jumps between $T_{\alpha_1}$ and values between $T_{\alpha_1}$ and $T_{\alpha_2}$ for higher $T_{\alpha_2}$, indicating a significant number of spikes at both of these times. This can be explained if one considers that the timing error increases with $\sqrt{T}$ [Eq. (14)]. If $T_{\alpha_1}$ stays constant while $T_{\alpha_2}$ increases, the difference in the timing errors decreases, and thus, the synaptic weights from both chains will be more similar. We conclude that a unique representation is possible for all $T_{\alpha_2} < T_{\alpha_1}$. For α ₂ < α ₁, however, this will be the most common situation, as the range of time intervals in ${\cal C}_{\alpha_2}$ is smaller than ${\cal C}_{\alpha_1}$ (cf. Fig. 5, left) and thus, $T_{\alpha_2} > T_{\alpha_1}$ rarely occurs.

6 Discussion

At first glance, the $\sqrt{T}$ dependency of the timing errors in a synfire chain seems to be incompatible with the experimental results of a constant or even increasing Weber fraction, which is a problem shared by many other models of timing. However, we identified a mechanism that makes the additional error plausible, namely the superlinear increase of the timing error with the transmission delay. Thus, we do not need to postulate any ad hoc assumptions about the scalar property, but could explain both the linear and faster-than-linear error increase from a constraint optimization process. Moreover, we found a neuronal implementation of this optimization by synaptic plasticity that also solved the problem of combining output from the various synfire chains to a unique representation of time.

A central assumption of our work is the limitation of the number of pools in a synfire chain. One possible reason why such a limit should apply is provided by a capacity argument. Synfire chains have been proposed to model the function of the cortical column (Herrmann et al. 1995; Bienenstock 1995), a structure containing 10⁴ to 10⁵ neurons. A pool size of the order of 10², comparable to the size of a minicolumn, has been shown to be necessary for stable propagation of the chain (Herrmann et al. 1995; Bienenstock 1995; Diesmann et al. 1999). Thus, the number of pools in the chain is constrained to the order of 10² to 10³. Of course, each of the neurons in a column could participate in multiple chains, but the capacity of network for synfire chains has been found to be limited (Herrmann et al. 1995; Bienenstock 1995), and it has been proposed that this capacity only allows for representation of events of durations up to 1s (Herrmann et al. 1995). However, all these studies assume a transmission delay of about 1 ms, which is true only for the fastest chains in our framework. Possible delays up to 6 ms do not seem to enable computations much above the range of one second, due to the increase in timing errors. In order to compensate these errors, an increase of the width of the chain were necessary which in turn reduces capacity. The one second range has also been found in physiological experiments with precise spiking patterns (Ikegaya et al. 2004), although the results of this study are disputed.

Another argument for a constrained pool number relates to the formation of synfire-like structures with a distribution of transmission delays. It has been shown that such structures might emerge from STDP learning in recurrent networks (Izhikevich 2006). In this study, the number of neurons in each of the “polychronous groups” was less than 20 in the mean, in a network of 1000 neurons. Much larger networks tended to become unstable. Although it seems to be possible to stabilize such groups by external guidance (Izhikevich 2006; Buonomano 2005), problems of unstable connections are likely to put a further constraint on the length of a chain.

Apart from the limited chain length and the general connectivity of several synfire chains projecting onto readout neurons, many of the assumptions used in this model can be relaxed. First of all, the synfire chains are allowed to contain a certain amount of recurrent connections, which introduce an additional source of error, but do not destabilize the propagation of activity (Herrmann et al. 1995). Second, it is not necessary to prewire the connections to the readout neurons in the way we have used here. Rather, this connectivity will arise spontaneously from an initially random wiring because of the synaptic competition. Consider each readout neuron being initially connected randomly to a certain fraction of the pools in a single chain. STDP will then only enhance the connections to those pools that are active slightly before the stimulus. But at the same time, homeostasis weakens all the connections, including those which were not enhanced. As a result, connections which does not fit into the scheme we have proposed for learning temporal representations end up to be very weak, and mights also be removed by means of synaptic turnover. Finally, the model is also robust to changes in the properties of the noise. Introducing a finite correlation length into the noise only adds a constant to the timing error and does not change the form of the error function. And even if the $\sqrt{T}$ law in the timing errors of the individual chains changed due to some properties of the noise, one would expect that this affected all of the synfire chains alike. So the selection of optimal chains would still work in this case and the U-shaped form of the Weber fraction would be preserved.

The combination of synfire chains and the readout network with plastic connections opens the possibility to explain some further phenomena of temporal processing. For instance, it has been shown that the subjective length of an interval depends on attention: If a timing task has to be performed in parallel with a second, non-temporal task, the duration of the interval is systematically underestimated (Grondin 2001). This can be explained by our model if attention is modeled by the level of activation in the synfire neurons. The mean membrane potential $\langle V \rangle$ is increased, and thus, the difference $\langle V \rangle - V_\mathrm{thr}$ is decreased. This decreases Δt (cf. Section 3.3.2 and Wennekers and Palm 1996). Conversely, decreased attention due to a parallel task decreases $\langle V \rangle$ and slows down the chain, resulting in an underestimation of intervals.

Moreover, temporal representations are subject to selective learning: If a participant is trained with stimuli of a certain duration, discrimination of that duration is improved after training, but this effect does not generalize to different intervals (Buonomano and Karmarkar 2002). This is also readily explained in the framework of our model: Training of a specific duration strengthens the connections of the responsible readout neurons with the pools that are active at this time, and suppresses the random connections to other pools by means of synaptic competitions. The learning experiment described in Section 5 can also be related to the paradigm of classical conditioning, where the initiating stimulus S ₀ corresponds to the conditioned stimulus (e.g., the ring of a bell or a flash of light) which can be learned to predict the unconditioned stimulus (e.g., food or an airpuff), corresponding to the training stimulus S ₁ in our case. This may also solve the problem that learning seems to occur on time scales that are much longer than those of the STDP learning rule (Shors and Matzel 1997). Note that there is no need to assume that S ₀ only activated the synfire chains and S ₁ only the readout network. If there is no such distinction, S ₁ would both mark the end of a first interval and the beginning of another, starting off a new volley of synfire activity. In this way, the apparent “reset” of the timing system could be explained (Buonomano and Karmarkar 2002).

Based on our results and earlier descriptions of neuronal structures and connections that might be relevant for temporal processing (Buhusi and Meck 2005), we sketch a hypothetical architecture of our model in the brain: Synfire chains are present in all areas of the neocortex, performing computational tasks like pattern storage (Herrmann et al. 1995; Bienenstock 1995) or compositional binding (Hayon et al. 2005). They have different transmission delays that might have been shaped during their formation by the time scale of the task they perform. As a by-product of their usual computation, the chains encode the temporal information of a real or imagined event. These distributed time representations are then projected onto a central readout network that is located in the striatum (Buhusi and Meck 2005). Distortions in the level of dopamine, as induced by certain drugs or Parkinson’s disease will strongly affect the function of the readout neurons and thus, also the timing performance (Buhusi and Meck 2005; Rammsayer 1999). The connections from the chains to the readout neurons are initially randomly distributed and are shaped by synaptic plasticity to implement an optimal, unique representation of time. Nevertheless, input from suboptimal chains will not be entirely suppressed, so the random connectivity remains an additional source of errors that can be further reduced by training.

Note that within this framework, it is improbable that there is a separate chain for each conceivable time interval, as we have assumed in Section 4. More likely, there will we a finite set of chains that represents an entire range of durations by using more than just their final pool. Of course, this violates the optimality condition Eq. (18) to some extent, introducing another error source. More specifically, the timing errors will not increase as smoothly as Eq. (19) implies, but there will be jumps in the error whenever a certain chain has reached its final pool and longer intervals must resort to the next chain. Interestingly, such jumps have indeed been observed in psychophysical experiments after excessive training (Kristofferson 1980). It seems that those jumps are normally masked by noise that is reduced by training. One possible source of this noise might be the random connectivity to the readout neurons, which can be refined by plasticity (see above). Furthermore, it is conceivable that the transmission delay of the chains can be fine-tuned by slightly changing the activity level (Wennekers and Palm 1996). This might explain short-term adaptation effects which occur at the presentation of sequences (Blaschke et al., submitted for publication). Mechanisms that are not contained in the current form of the model include memory and decision.

A quantitative view on the Weber fraction calculated from our simulation data (cf. Fig. 5, left) reveals that its value of 0.5 to 4% of the represented interval is too low compared to the psychological experiments, which report values between 2 and 20% (Gibbon et al. 1997; Drake and Botte 1993; Getty 1976). This is due to a relatively low level of synaptic noise (σ _V  = 1.4 mV in our study compared to e.g. 2.85 mV in Diesmann et al. 1999). We conducted tests of whether this level can be increased while maintaining stable propagation of chains. Preliminary results show that this is possible by compensation of the increased noise with increasing both synaptic weights w _S and the firing threshold V _thr. Using these measures, the Weber fraction is increased to values between 3 and 9%. A full exploration of the synfire parameter space is beyond the scope of this research, but it seems that at least the lower range of the Weber fractions experimentally observed can be obtained within the biologically realistic range of parameters. Some additional error sources have been mentioned in this section.

Finally, we note that our framework is not necessarily limited to synfire chains. Any timing system with a limited dynamic range will show a similar effect, given that this range can be extended at the cost of a superlinear increase in the timing error. The optimization scheme and readout network will be the same in this case. It seems worthwhile to check this properties for neurocomputational models of timing such as state-dependent networks (Buonomano 2000, 2005), ramping activity (Durstewitz 2006) or the striatal beat model (Matell and Meck 2004). The convergence of evidence from psychology and neuroscience is likely to decide which classes of models are able to explain how our brain tells time.