Abstract
Recent singlecell studies in monkeys (Romo et al., 2004) show that the activity of neurons in the ventral premotor cortex covaries with the animal's decisions in a perceptual comparison task regarding the frequency of vibrotactile events. The firing rate response of these neurons was dependent only on the frequency differences between the two applied vibrations, the sign of that difference being the determining factor for correct task performance. We present a biophysically realistic neurodynamical model that can account for the most relevant characteristics of this decisionmakingrelated neural activity. One of the nontrivial predictions of this model is that Weber's law will underlie the perceptual discrimination behavior. We confirmed this prediction in behavioral tests of vibrotactile discrimination in humans and propose a computational explanation of perceptual discrimination that accounts naturally for the emergence of Weber's law. We conclude that the neurodynamical mechanisms and computational principles underlying the decisionmaking processes in this perceptual discrimination task are consistent with a fluctuationdriven scenario in a multistable regime.
 somatosensory
 vibrotactile discrimination
 decision making
 Weber's law
 multistability
 stochastic neurodynamics
 probabilistic behavior
Introduction
A recurrent challenge in neuroscience is to understand how the brain represents the external world through the stimulation impinging on our sensory epithelia. A crucial attempt in this direction was the conceptualization introduced by E. H. Weber (1795–1878), who described the relationship between the physical magnitude of stimuli and the perceived (subjective) intensity of the sensation. Given the difficulties inherent in estimating the absolute strength of sensations, Weber operationalized the perception of stimulus change in intensity by mental comparison of two or more sensations. Perceptual comparison allows one to discern the minimal physical difference in stimulus intensity that produces a detectable change in sensation (i.e., the difference threshold). Weber's law states that the ratio between the difference threshold and the background stimulus intensity is a constant. This law has been extremely successful in describing human reactions to a wide range of sensory stimuli (Oberlin, 1936; Gaydos, 1958; Consweet and Pinsker, 1965; Fechner, 1966; Indow and Stevens, 1966; Panek and Stevens, 1966; Leshowitz et al., 1968) (see also Laming, 1986), but its neural underpinnings are still far from clear. Here, we address the potential neurophysiological implementation of Weber's law in relation to decision making.
Recent neurophysiological studies provide a good description of the neural correlates associated with the discrimination of vibrotactile stimuli within the flutter range (5–40 Hz). Vibrotactile stimulation has been often used to study the correlates of the sensation of roughness in texture perception (Johnson and Hsiao, 1992) as well as implemented in highly controlled laboratory settings to study neural coding and decision making. In particular, Romo and colleagues (Romo and Salinas, 2001, 2003; Hernandez et al., 2002; Romo et al., 2002, 2003, 2004) have studied the neural mechanisms underlying perceptual comparison by measuring singleneuron responses in monkeys trained to compare two mechanical vibrations applied sequentially to the fingertip. As in any psychophysical task, this discrimination implied a certain amount of training (both in humans and monkeys), but it has been observed that after a sharp initial improvement attributable to cognitive processing, performance usually attains an asymptotic level that most likely reflects true discrimination abilities based on perceptual representations (Hernandez et al., 1997). Using this task, Romo et al. (2004) found neurons in the ventral premotor cortex (VPC) of the monkey the firing rate of which was dependent only on the difference between the two applied frequencies, the sign of that difference being the determining factor for correct task performance. These neurons reflect the implementation of the perceptual comparison process and, to a large extent, seemed to underlie the cognitive process of decision making.
Deco and Rolls (2006) designed a neurodynamical model to account for the neurophysiological data of Romo et al. (2004) that can account for the most relevant characteristics of the decisionrelated neural activity. The behavior of the model predicts that Weber's law will apply to perceptual discrimination. This is not a trivial prediction, given that contrary to most other demonstrations, here Weber's law is predicted in a domain different from stimulus intensity, namely vibrotactile frequency. The results of previous research (in humans and monkeys) are somehow mixed with respect to whether vibrotactile frequency discrimination does or does not abide by Weber's law (Goff, 1967; Mountcastle et al., 1972; Gescheider et al., 1990; Mahns et al., 2006). This is especially relevant given that such attempts have not always considered the important physiological differences between the sensory pathways encoding vibrotactile stimuli in the flutter range and those encoding higher frequency vibrations. Here we performed behavioral tests of vibrotactile discrimination in humans and propose a computational explanation of perceptual discrimination that accounts naturally for the emergence of Weber's law. The human behavioral data complements the current neurophysiological evidence, helping to constrain the underlying neurodynamics and computational principles involved in perceptual discrimination.
Materials and Methods
Historical background
Johnson (1980a,b) carefully analyzed and formulated the relationship between neurophysiological measurements (reflecting the neuronal representation of sensory stimulation along the somatosensory path) and the psychophysically observed probabilistic subject's behavior. In particular, he proposed a rigorous theoretical framework including as explicit variables the experimental setup, the neural representations of the stimuli, variance of the neural representations, and functional dependence on the stimuli that they represent. In this theoretical framework, he distinguished between: (1) neuronal processes encoding peripheral stimuli information on which the discrimination is based, and (2) neuronal processes underlying the decisionmaking discrimination according the experimental setup that he considered. In this way, Johnson (1980a,b) was able to relate Weber's law with the underlying neuronal mechanisms.
Nevertheless, Johnson (1980a,b) assumed linear mechanisms and assumed that the neural representation of the stimuli results from a stimulusevoked activity and noise. In this study, we complement and extend the theory of Johnson (1980a,b) by offering a concrete realistic biophysical neuronal network implementing decision making in a discrimination task. We address how the dynamics of the neuronal network can transform explicitly the underlying neuronal activity and noise in the specific subject's behavior. Even more, we will show that the dynamical working point of the neuronal network is fundamental and intrinsically related with the probabilistic behavior of the subjects.
In other words, we assume, like Johnson (1980a,b), that a hallmark of neurophysiological measurements is the high degree of variability in the neuronal activity both within and between trials. Rather than attributing this variability to poor sampling, we suggest that the stochastic fluctuations in the neuronal dynamics may actually have a functional role. Indeed, the fact that perception and behavior during certain types of tasks can be well described by probabilistic models suggests a link between the stochasticity at the cellular and behavioral level. The main aim of the current study is to elucidate the mechanisms underlying this link by constructing computational models that account for measurements both at the cellular and behavioral levels.
Theoretical framework: neurodynamical model
Theoretical models can help understand the underlying neurodynamical and computational mechanisms responsible for the evaluation of perceptual evidence that lead to decision making (Brody et al., 2003; Machens et al., 2005; Deco and Rolls, 2006). Such models are constrained by, and consistent with, extant neurophysiological data (Romo and Salinas, 2001, 2003; Romo et al., 2004). They typically involve two populations of excitatory neurons engaged in competitive interactions mediated by inhibition. Sensory input may bias the competition in favor of one of the populations, potentially resulting in a gradually developing decision in which neurons belonging to the winning population exhibit increased activity, whereas activity in the other population is inhibited.
We modeled the neural, synaptic, and cortical dynamics underlying the computation of perceptual discrimination using the theoretical framework of attractor networks (Brunel and Wang, 2001) based on the principle of biased competition/cooperation [see previous applications (Rolls and Deco, 2002; Deco and Rolls, 2005a,b)]. This model was designed to account for the activity of neurons implementing the comparison process found in the VPC, as evidenced by singlecell recordings and the corresponding behavioral measures in monkeys. The response of these neurons during the comparison period of the task depends on the frequency difference between the two applied vibrotactile stimuli, the sign of that difference being the determining factor for correct task performance. The computational model of Deco and Rolls (2006) implements a neuronal network that can reproduce the decisionmaking response selectivity of the mentioned VPC neurons. Competition and cooperation mechanisms are implemented in an attractor network consisting of two recurrently connected populations of excitatory neurons, mutually connected with a common inhibitory population.
The model enables a proper description of the transients (nonstationary) and probabilistic nature of behavior (performance) by the explicit use of spiking and synaptic dynamics of onecompartment integrateandfire (IF) neuron models (Tuckwell, 1988) at the microscopic level. This allowed us to use realistic biophysical time constants, latencies, and conductances to model the synaptic current, enabling the thorough study of the time scales and firing rates involved in the evolution of the neural activity. The IF neuronal cells modeled here have three types of receptors mediating the inflow of synaptic currents: the EPSCs are mediated by AMPA (fast) and NMDAglutamate (slow) receptors, whereas external EPSCs imposed onto the network are driven by AMPA receptors only. IPSCs to both excitatory and inhibitory neurons are mediated by GABA receptors [for details of the mathematical formulation, see Brunel and Wang (2001) and Deco and Rolls (2005a, 2006)].
The attractor network implementing the comparison mechanism (Fig. 1) is composed of interacting neurons organized into a discrete set of populations (i.e., groups of excitatory or inhibitory neurons sharing the same inputs and connectivities). The network contains 800N_{E} (excitatory) pyramidal cells and 200N_{I} inhibitory interneurons, in accordance with the proportion of 80% pyramidal cells versus 20% interneurons typically observed in physiological studies (Abeles, 1991; Rolls and Deco, 2002). In this minimal model, specific populations encode the categorical result of the comparison between the two sequentially applied vibrotactile stimuli, f1 and f2 (i.e., f1 > f2 or f1 < f2). Each specific population of excitatory cells contains rN_{E} neurons (here we use r=0.1). There is also a nonspecific population (labeled “Nonspecific”), which groups all other excitatory neurons in the modeled cortical area, not specifically involved in the present task, and one inhibitory population (labeled “Inhibitory”) grouping the local inhibitory neurons in the modeled brain area. The latter population regulates the overall activity and implements competition in the network by spreading a global inhibition signal.
The conductance values for the synapses between pairs of neurons is specified by connection weights, which can deviate from their default value, 1. The structure and function of the network is achieved by differentially adapting these synaptic strengths within and between populations of neurons. The labeling of the weights is defined in Figure 1. We assume that the connections are already formed (e.g., by earlieroccurring selforganization processes such as Hebbian learning). We assume that the two possible outcomes of the decision, f1 > f2 and f1 < f2, corresponding to the two response categories, are already encoded, in the sense that the monkey or human is already trained to push one or the other button, but not both (to obtain some reward). As a consequence of this, neurons within a specific excitatory population are mutually coupled with a strong weight w_{+}. Furthermore, the populations encoding these two decisions are likely to have anticorrelated activity in this behavioral context, resulting in weaker than average connections between them. Consequently, we choose a weaker value w_{−}=1 − r(w_{+} − 1)/(1 − r), so that the overall recurrent excitatory synaptic drive in the spontaneous state remains constant as w_{+} is varied (Brunel and Wang, 2001). Neurons in the inhibitory population are mutually connected with an intermediate weight w=1. They are also connected with all excitatory neurons in the same layer with the same intermediate weight, which for excitatorytoinhibitory connections is w=1 and for inhibitorytoexcitatory connections is denoted by a weight w_{I}. Neurons in a specific excitatory population are connected to neurons in the nonselective population in the same layer with a feedforward synaptic weight w=1 and a feedback synaptic connection of weight w_{−}. Each individual population is driven by two different kinds of input. First, all neurons in the model network receive spontaneous background activity from outside the module through N_{ext}=800 external excitatory connections. Each connection carries a Poisson spike train at a spontaneous rate of 3 Hz, which is a typical value observed in the cerebral cortex. This results in a background external input with a rate of 2.4 kHz for each neuron. Second, the neurons in the two specific populations additionally receive external inputs encoding the stimulusspecific information. These inputs are assumed to originate from the second somatosensory cortex (S2) and from the prefrontal cortex (PFC), encoding the frequency of both stimuli f1 (stored) and f2 (present) to be compared during the comparison period (i.e., when the second stimulus is applied). As described in neurophysiological studies, there are two main types of S2 and PFC neurons, namely neurons whose firing rate exhibits a positive monotonic relationship with stimulus frequency and others in which this relationship is negative. Based on the experimental results (Romo et al., 2004), we model the firing rate of positive monotonic neurons as f_{+}^{x}=5 + 2.3f_{x} Hz, and the firing rate of negative monotonic neurons as f_{−}^{x}=25 − 0.6f_{x} Hz, where f_{x} is the frequency of the vibrotactile stimulation in Hz (i.e., f_{x} is equal to f1 or f2). When stimulating, the rate of the Poisson train to the neurons of both specific populations f1 > f2 and f1 < f2 is increased by an extra value of λ_{1}=f_{+}^{1} + f_{−}^{2} and λ_{2}=f_{−}^{1} + f_{+}^{2}, respectively, coding the two vibrotactile stimuli to be compared.
The stationary states of a network of integrateandfire neurons can be exhaustively studied using a reduced consistent mean field, to simplify the integrateandfire equations by replacing, after the diffusion approximation (Tuckwell, 1988), the sums of the synaptic components by the average component and a fluctuation term. The stationary dynamics of each population can be described by the population transfer function, which provides the average population rate as a function of the average input current. The set of stationary, selfreproducing rates ν_{i} for the different populations i in the network can be found by solving a set of coupled selfconsistency equations. This enables the selection of the parameter region that reveals the desired emergent behavior in a bifurcation diagram.^{a} In the present case, the essential requirement is that, for the stationary conditions, different possible attractors are stable. The attractors of interest for our task correspond to the activation (high spiking rates) or nonactivation (low spiking rates) of the neurons in the specific populations f1 > f2 and f1 < f2. The activation of the specific population f1 > f2 (f1 < f2) and the simultaneous lack of activation of the complementary population f1 < f2 (f1 > f2) correspond to an encoding “single state” associated with a motor response reporting the categorical decision f1 > f2 (f1 < f2). The lack of activation of both specific populations (“spontaneous state”) would correspond to an encoding state that cannot lead to a behavioral decision (i.e., there is no answer, or a motor response is generated completely randomly). The same happens if both specific populations are activated to the same degree (“pair state”). Because responses in animals are probabilistic in nature, the operating working point of the network should be such that both possible categorical decisions (i.e., both possible single states) are bistable. In addition, we will show that the model predicts a behavior consistent with Weber's law if, and only if, the spontaneous state is also a stable state (i.e., when the dynamical working point of the network is in a regime of multistability). In this way, the empirical confirmation of Weber's law informs about the dynamical working point of the network. We use meanfield techniques for analyzing the nonstationary asymptotic states via a reduced model that allows us to select easily the regime of multistability. The mathematical formulation of the integrateandfire neurons and synaptic currents as well as the corresponding consistent meanfield is summarized in the supplemental material (available at www.jneurosci.org) (see also Deco and Rolls, 2006).
Experimental human behavior
Participants.
We tested a group of eight human participants (range, 18–35 years of age) who volunteered for the study. All were right handed and reported normal tactile sensitivity.
Apparatus and materials.
Participants sat ∼60 cm away from a computer cathode ray tube monitor in a dimly lit, soundattenuated testing room. They were instructed to direct their gaze to a fixation cross at the center of the screen. Stimulation was delivered using a bone conduction vibrotactile stimulator (OticonA bone conduction vibrators; 3.8 cm^{2} vibrating surface; Oticon, Hamilton, UK). The protocols were programmed using Expe6 software (Pallier et al., 1997) and run on a Pentium (Intel, Santa Clara, CA) PC computer. The vibration waves were generated with the PC sound card and fed to the vibrators through an amplifier (PioneerA307; frequency response, 5 Hz to 100 kHz). The stimulator was attached to the distal pad of the left annular finger using a velcroed cloth loop, so the vibrator location and pressure was held constant throughout the experimental session. The subjects wore headphones (HD435; Sennheisser, Wedemark, Germany) delivering white noise at constant intensity throughout the experiment, sufficient to mask any sound produced by the vibrotactile stimulators.
Procedure.
The frequency discrimination experiment consisted of a twointerval, forcedchoice task in which two 500 ms vibrotactile stimuli (the base and the comparison frequencies, order randomized, 500 ms interstimulus interval) were presented to the participant for classification. At the beginning of each trial, the computer screen showed a “−” sign at the center, which was temporarily replaced by an “*” sign during the presentation of each vibration. After the second vibration, a “?” sign appeared in the center of the screen, until a response was made. The task was to press the key (1 or 2) corresponding to the stimulus of higher frequency (first or second, respectively). Five different base frequencies were tested (20, 30, 40, 60, and 80 Hz), although only the frequencies <50 Hz (flutter) were used to test the model. This was done in an attempt to engage the same type of somatosensory receptors and pathways (Meissner receptors) throughout all conditions (Werner and Mountcastle, 1965; Mountcastle et al., 1967; Freeman and Johnson, 1982). Each particular base frequency was confronted with eight different comparison frequencies (±2, ±4, ±6, and ±8 Hz) in a separate block of 112 trials (eight combinations of base and comparison frequency presented 14 times each). No feedback was provided to participants during the experimental phase. Before this vibration discrimination test, the intensities of the different vibrotactile stimuli were carefully adjusted individually for each participant in the following way.
Given that perceived intensity of vibration varies across frequencies, certain adjustments were necessary to ensure that participants based their judgments on frequency. First, we determined the simple detection threshold for each base frequency using interleaved staircases based on the parameter estimation by sequential testing (PEST) algorithm to adjust the gain of the sound card for 50% detection performance. After this procedure, we raised the output level of the amplifier by 15% and ran a second phase to find the point of subjective equality in intensity across vibrotactile frequencies. In particular, we gradually adjusted the output of the sound card [sound pressure level (SPL)] individually for each comparison frequency in order for it to feel as intense as its corresponding base frequency (we used interleaved PEST staircases). Then, for each participant and base frequency, we found the best fitting line (linear regression) to describe the relationship between frequency (for all comparison frequencies) and SPL and selected the appropriate intensity values to be used in the discrimination experiment for each comparison frequency.
Once the particular intensity values for each frequency and subject had been adjusted, and before the frequency discrimination task, participants received extensive training with visual feedback presented on the computer screen. In an initial training block, we used large frequency differences (>30 Hz) to acquaint participants with the task. In a second training block, participants were introduced to smaller frequency differences (20 ± 2 Hz). The third training block was divided into subblocks corresponding to each base frequency, and the tests always included the base frequency paired with the ±8 Hz comparison frequencies. Each of these subblocks was repeated until a performance of 90% was achieved. Then, each subblock was repeated once or twice more again, just before the experimental block corresponding to its base frequency. It took between 5 and 6 h to test each individual, and the experiment was divided into two sessions. All the experimental blocks took place within the second session.
Results
Model results
Figure 2A shows the probability of classifying the frequency of the comparison stimulus f2 as higher than the frequency of the base stimulus f1, as a function of the frequency to be compared. Each data point corresponds to average across 200 trials of the full spiking simulations [for the plotting of typical evolutions of the spiking network of VPC neurons during the comparison period and additional details, see Deco and Rolls (2006)]. The lines were calculated by fitting the data points using a Weibull function. A correct classification occurs when during the 500 ms comparison period, the network evolves to a “singlestate” attractor that shows a high level of spiking activity (>10 Hz) for the population f1 > f2 and simultaneously a low level of spiking activity for the population f1 < f2 (at the level of the spontaneous activity). One can observe from the different panels corresponding to different base frequencies of f1 that to reach the 85% performance threshold (top horizontal dashed line), the difference between f1 (base frequency) and f2 (comparison frequency) must become larger as f1 increases. Figure 2B plots the justnoticeable difference (JND) corresponding to the difference limens calculated as onehalf the difference between the frequency identified as higher than the standard on 85% of the trials and on 15% of the trials. The JND increases linearly as a function of the base frequency. This linear increase corresponds to Weber's law for the present vibrotactile discrimination task.
Behavioral results
We assessed the proportion of comparisonhigher responses (f2 > f1) as a function of base frequency for each participant and then adjusted a Weibull function to the observed data (for individual results, see supplemental material, available at www.jneurosci.org). From there, we calculated the average of individual JND thresholds (difference limen calculated as onehalf the difference between the frequency identified as higher than the standard on 85% and on 15% of the trials). The relationship between JND and base frequency is plotted in Figure 3.^{b} The plot reveals a positive correlation between stimulus magnitude (base frequency) and JND, thereby illustrating Weber's law in vibrotactile flutter frequency discrimination for humans. In particular, the JND for base frequency 20 Hz (2.4 Hz) was smaller than the JND for 30 Hz (4.1 Hz; p=0.035) and than the JND for 40 Hz (6.5 Hz; p=0.024), whereas the JND for 30 Hz was smaller than the JND for 40 Hz (p=0.025). The average Weber's fraction (k) values were 0.122, 0.138, and 0.162 for the base frequencies of 20, 30, and 40 Hz, respectively. None of the paired comparisons (uncorrected) between these k values were significant (all, p > .25). Therefore, the estimated value of (k) for vibrotactile (flutter) discrimination is 0.140.
Computational principles underlying Weber's law
To extract the neurodynamical mechanisms underlying a Weber's law behavior in decision making, we analyzed the computational capabilities of the multistable network model by means of a reduced model. If circuits exhibiting multistability are comprised of large numbers of spiking neurons, the fluctuations needed to drive the transitions arise naturally through noisy input and/or disorder in the collective behavior of the network. However, such dynamics can also be qualitatively captured in a system of nonlinear coupled differential equations that describe the evolution of the average firing rate of each population (meanfield reduction). In this case, a fluctuation term must be added to drive the transitions. Such a minimal probabilistic decisionmaking network consists of two distinct populations of neurons whose activity encodes the two alternative choices (Fig. 4). Neurons within a specific population interact via strong recurrent excitation with weight w_{+}. Neurons in one population are mutually coupled to all other neurons in the other population in an inhibitory manner with a weight w_{I}. Throughout the study, the excitatory weight is set to w_{+}=1, and the inhibitory weight is set to w_{I}=1.
The temporal dynamics of the firing rates of the neuronal populations can be qualitatively captured via a system of firstorder differential equations of the WilsonCowan type (Renart et al., 2003; La Camera et al., 2004). For the decisionmaking architecture shown in Figure 4, the firing rate equations are given by the following: where ν_{i} denotes the firing rate of population i, w_{ij} denotes the synaptic strength between populations i and j (i.e., w_{11}=w_{22}=w_{+}, and w_{12}=w_{21}=−w_{i}), and τ=1 ms. The external, sensory input to the population i is denoted by λ_{i}. Throughout the study, λ_{1}=Aλ + B, whereas input to the second population is set to λ_{2}=A(λ + Δλ) + B (here we use A=1/60 and B=8/3, so that the multistable regime in the model falls in the flutter range). We will refer to the case of Δλ=0 Hz as the unbiased case. The nonlinear transfer response function φ(·) is sigmoidal. In this work, we use with a=15, and b=0.25. Fluctuations are modeled via an additive Gaussian noise term denoted by ξ_{i}. Here 〈ξ_{i}(t)〉=0 and 〈ξ_{i}(t)ξ_{i}(t′)〉=β^{2}δ_{ij}δ(t − t′), where the brackets 〈…〉 denote the average over stochastic random variables. Here we used the value of β=0.5. This noise term represents finitesize effects that arise because of the finite number N of neurons in the populations.^{c}
In the absence of noise, the fixed points of Equation 1 can be determined by setting the time derivative equal to zero and solving for ν_{i}. The noisefree case corresponds to the limit N→∞, which leads to the socalled classical meanfield approximation (Tuckwell, 1988; Amit and Brunel, 1997; Brunel and Wang, 2001), and a standard bifurcation analysis of the fixed points can be performed. Figure 5 plots the bifurcation diagram showing the fixed points of Equation 1 as a function of the input parameter λ for the unbiased case. For decision making in the multistable regime, the relevant fixedpoint solutions are the spontaneous state and the two states representing a decision, henceforth referred to as decision states. In particular, in the unbiased case, as one increases λ, the spontaneous state in Equation 1 loses stability at a critical value λ_{c1}, and the system undergoes a bifurcation to the decision state solutions. There exists, then, a region of multistability between the spontaneous state and the decision states defined by the interval (λ_{c1}, λ_{c2}), where λ_{c2} is the value of λ at which the stable decision state annihilates with the unstable fixed point in a saddlenode bifurcation.
The different dynamical regimes that can be observed in the bifurcation diagram are associated with different ways of computing decision making (Fig. 6 represents schematically the different computational principles associated with the different dynamical regimes). The xaxis represents the neuronal activity of one of the populations (ν_{i}), and the landscape represents an energy landscape regulating the evolution of the system. Three different dynamical regimes are shown: stable spontaneous state, multistable, and bistable. For values of λ < λ_{c1} (Fig. 6, left), only the spontaneous state is stable, whereas the decision states do not appear (for the unbiased case). For increasing Δλ (biased case), one decision state (corresponding to the choice where the increased value of λ + Δλ is applied) emerges, attracting the dynamics toward this decision states. For values between λ_{c1} and λ_{c2} (Fig. 6, middle), there is a region of multistability between the spontaneous state and the decision states. In this interval, the fluctuations are responsible for driving the system from the initial stable spontaneous state to one of the two decision states corresponding to the two possible response choices. Thus, in this scenario, fluctuations play a crucial role in the computation of decision making. For values of λ > λ_{c2} (Fig. 6, right), we find a region of bistability in which the initial spontaneous state is unstable, and only the two decision states are stable. In this regime, the spontaneous state destabilizes, so that the dynamics of the network rapidly evolves toward one of the two decision states, resembling therefore a pure diffusion model integrating the relative evidence for one choice over another.
These three different dynamical regimes involving different computational principles can be distinguished at the behavioral level. Figure 7 characterizes the responses of the theoretical model associated with a decisionmaking task for the different dynamical regimes. The figure shows the critical discrimination Δλ value corresponding to an 85% correct performance level (“difference threshold”) as a function of the base frequency λ. For the first region (“spontaneous state stable,” λ < λ_{c1}), the decision states do not exist for the unbiased case. Basically, the network is not able to reach a decision in the first regime. Accordingly, Δλ has to be large enough to bring the network into the second regime. The larger the base frequency λ, the smaller the necessary value of Δλ. As a result, JNDs decrease with increasing base frequency. Thus, in this region, the difference threshold shows an experimental inconsistent linear but negative correlation with the base frequency λ. For the second region (“multistable,” λ_{c1} < λ < λ_{c2}), corresponding to a fluctuationdriven multistable scenario, a perfect linear correlation between the difference threshold and the base frequency is observed. This behavior corresponds to Weber's law, and as we demonstrate experimentally in this paper, is consistent with the observed behavior. For the third region (“bistable,” λ > λ_{c2}) corresponding to a pure diffusion process, a deviation from Weber's law is observed. The difference threshold starts to show a nonlinear dependence with the base frequency (for initial behavioral evidence, see supplemental material, available at www.jneurosci.org).
It is important to stress that the effect of the noise is particularly relevant in the multistable regime, because the fluctuations are the driving force that allow the system to escape the decision barriers around the stable spontaneous state. In the multistable scenario, the choices are associated with stable attractors, and the starting condition is also given by a stable spontaneous state. To make a decision, the system has to escape the stable spontaneous state toward one of the choice attractors [this is related to the socalled “Kramers escape problem” (Kramers 1940)]. In contrast, in the bistable regime (the socalled “ballistic” regime; see Discussion), the noise is of course relevant as the basis of the diffusion process, but it is not the main driving force. This is because in the bistable scenario the spontaneous state is not a stable state, and therefore with or without noise, the system will necessarily evolve to one or the other decision attractors just because of the neurodynamical flow. We investigated the effect of noise on both regimes numerically and found that for weaker noise, the Weber's fraction increases in the multistable regime, whereas the nonlinearity (namely the saturation observed for large λ) is accentuated. Additional effects of the noise on the bifurcation and an indepth discussion about analytical ways of studying it for the cases of neuronal networks implementing decision making have been described by Deco and Marti (2007a,b).
Discussion
The main finding to emerge from the theoretical analysis of the model is that Weber's law behavior for the decision making involved in the discrimination of two vibrotactile stimuli is consistent with an underlying fluctuationdriven scenario in a multistable regime. One nontrivial prediction of the model (see Deco and Rolls, 2006) is a Weber's law behavior for vibrotactile frequency discrimination. An interesting aspect of this prediction is that, if the connectivity parameters of the network are tuned using meanfield techniques (so that the network has two possible stable stationary final attractors respectively related to the two possible decisions), then the firing rate of the neurons in the winning attractor reflects the difference in the frequencies (Δf) being compared, but not the absolute frequencies. Therefore, Weber's law for frequency comparison is not explicitly encoded by the firing rate of these attractors. An analysis of the nonstationary evolution of the dynamics of the network model shows that Weber's law is implemented in the probability of transition from the initial spontaneous firing state to one of the two possible attractor states. In this way, statistical fluctuations caused by finite size noise play a crucial role in the decisionmaking process. The model establishes that a Weber's law for frequency discrimination could be implemented not by the firing rate of a given population of neurons (which reflects just Δf), but by the probability that that particular population will be activated, which is influenced by both Δf and the absolute value of the frequency. Thus, the model shows that the neural correlates underlying perceptual discrimination are consistent with a decisionmakinglike scenario of fluctuationdriven computation that causes a probabilistic transition between multistable states.
The finding of Weber's law in terms of vibrotactile frequency discrimination is not trivial because, in contrast to most previous implementations of the Weber's law, frequency does not involve a perceptual correlate of the physical intensity (or strength) of a magnitude, but rather of the rate of stimulation. Therefore, the prediction for Weber's law when performing vibrotactile discrimination is not directly obvious. Moreover, the finding that the Weber's law applies in the domain of vibrotactile frequency can be very informative as to the neural correlates of perceptual comparison. Indeed, previous attempts to establish Weber's law for vibrotactile frequency have provided mixed results. For instance, Gescheider et al. (1990) reported that intensity difference thresholds on the finger in humans did not vary as a function of frequency (25 vs 250 Hz). Similarly, Tommerdahl et al. (2005) reported that the Weber fraction (the ratio between difference threshold and absolute stimulus frequency), which should remain relatively constant across frequencies if Weber's law were to hold, varied considerably from 25 to 200 Hz base frequencies. In contrast, however, Mahns et al. (2006) (see also Rothenberg et al., 1977) found that the Weber fraction remained relatively constant across different base frequencies (they investigated glabrous as well as hairy skin, with frequencies from 20 to 200 Hz) [but see Goff (1967) for a slightly different result]. Nevertheless, one potential challenge when comparing threshold across such different base frequencies could be that the somatosensory channels being excited by the different stimuli are different, thereby confounding the psychophysical results. Mountcastle et al. (1972) reported the Weber fractions of humans and monkeys performing a frequency discrimination task with base frequencies ranging from 30 to 200 Hz, but that included several base frequencies in the range of flutter. Their human Weber fractions were relatively variable across base frequencies even within the flutter range, but they followed Weber's law in that difference limens did not have a particular trend toward increasing or decreasing with base frequency. However, the study by Mountcastle et al. (1972) might have been subject to a methodological problem whereby the observers (monkey or human) could have solved the task without necessarily performing a mental comparison (Hernandez et al., 1997). In a more recent psychophysical study with monkeys, Hernandez et al. (1997) addressed this problem, but the pattern of discrimination thresholds was not perfectly clear as to how much the Weber's law held. Here we found a clear positive relationship between difference threshold and base frequency within the flutter range in humans. Note, however, that the slightly worse discrimination values reported in the present study, as compared with other previous studies, are the result of using a threshold value of 85% to calculate our just noticeable differences, rather than the 75% threshold values used in some earlier studies. In fact, the Weber fraction across the frequencies tested remains statistically stable.
The attractor model described here computes decision making via competitivebased bias in a bistable neurodynamical system. We show experimental evidence supporting the prediction of this model, namely that a decisionmaking task follows Weber's law. In general, the dynamics relevant for decision making in these networks depend on the stability of the spontaneous activity state (i.e., the state in which no decision has yet been made). If, once the second stimulus is presented, the spontaneous state destabilizes, then the dynamics rapidly evolves toward one of the two decision states (Wong and Wang, 2006). This scenario is consistent with recent theoretical studies of behavioral data using the socalled diffusion models (Smith and Ratcliff, 2004). In these models, it is assumed that information that drives the decision process is accumulated continuously over time until a decision boundary is reached. Given the success of diffusion models in explaining behavioral data, it seems likely that some decisionmaking processes in the nervous system indeed rely on a similar accumulation of evidence. Alternative phenomenological models have been developed in which the effective dynamics is equivalent to an OrnsteinUhlenbeck process with fixed boundaries (Usher and McClelland, 2001). Such connectionist models differ from the classical diffusion model in that the drift of the decision variable is proportional to the value of the variable itself (i.e., it can be leaky or repelling; ballistic model) (Brown and Heathcote, 2005). The diffusion thus occurs not on a flat landscape but on a curved one. Furthermore, through fine tuning of higherdimensional systems to encode lowdimensional choices, one can also construct an approximate line attractor (Machens et al., 2005). In this case, the dynamics is close to pure diffusion and is thus a true integration of the relative evidence for one choice over another. In fact, the nonlinear dependence of JND as a function of λ in this bistable regime is perhaps not surprising, given that it results from an underlying nonlinear diffusion process.
On the other hand, we show here that both neurophysiological and behavioral evidences suggest that decision making results from an alternative fluctuationdriven multistable scenario. This scenario occurs when the spontaneous state does not lose stability but is instead bistable with the decision states, hence leading to multistability between three possible fixed points. Locally, then, the dynamics is similar to the “leaky” or attracting OrnsteinUhlenbeck process. Such multistability only occurs if the recurrent excitation within each neuronal group or population is strong enough. In this case, only a sufficiently strong perturbation would drive the system from the stable spontaneous state to one of the two decision states. This differs from the earlier ballistic scenario, in which the system will evolve toward one of the two choices even in the absence of fluctuations. Thus, in the multistable regime, fluctuations, perhaps noisedriven, are essential for decision making. Such a multistable region has also been evoked in a model of decision making as a means of holding the decision state in working memory once it has been made (Wong and Wang, 2006). In the model of Wong and Wang (2006), the input reflecting the evidence for the two alternatives destabilizes the spontaneous state, thus leading to a decision. Removal of the input brings the system into the multistable regime in which hysteresis ensures that the elevated level of activity in the neuronal group encoding the decision will be maintained.
The computational analysis of the network model showed that a Weber's law behavior for decision making is only consistent with a fluctuationdriven scenario in a multistable regime. Hence, the experimental behavioral evidence showing a Weber's law behavior for the decision making involved in the discrimination of two vibrotactile stimuli in humans suggests that the neurodynamical mechanisms and computational principle underlying this process is consistent with a fluctuationdriven scenario in a multistable regime. More specifically, the computational analysis performed here allows us to reveal that the way in which a decision is computed depends strongly on the dynamical regime in which the system is functioning. Neurophysiological data are compatible with both types of scenario being analyzed (i.e., Kramers multistable scenario and ballistic diffusion scenario). On the other hand, the psychophysical results obtained here with humans seem to suggest that Weber's law behavior is more consistent with an underlying multistable Kramers scenario. Nevertheless, we do not exclude the bistable scenario. In fact, we believe that, depending on the range of the input used (λ) and learning (which can shift the dynamical working regime, by modifying w_{+} for example), the working regime and related underlying scenario could change. This would explain a breaking of Weber's law (i.e., a change in the slope of the linear relation) for larger λ (as evidenced experimentally), and at the same time offers a concrete prediction: the reaction time distribution should be different in both scenarios. Thus, integrating behavioral and neurophysiological data, we are able to further increase our understanding of the computational principles underlying decision making.
Footnotes

This work was supported by the European Union Grant EC005024 (Specific Targeted Research Project“Decisions in Motion”) and by the Spanish Research Projects TIN200404363C0301/02. G.D. and S.S.F. were supported by Institució Catalana de Recerca i Estudis Avançats. We acknowledge the referees and Reviewing Editor for their helpful comments to previous versions of this manuscript.

↵a Note that the dynamics of the whole system depend on two parameters regulating the level of competition/cooperation, namely w_{i} and w_{+}. We characterized the network's different modes of operation corresponding to different parameter regimes by exploring these two connecting weights. In other words, these connecting weights are fixed by selecting the dynamical region in the bifurcation diagram that is consistent with the data [for more details, see Deco and Rolls (2006)].

↵b Figures 2 and 3 should be compared qualitatively. Both figures (i.e., model and experimental results) show a behavior consistent with Weber's law. It was not the aim of this study to perform data fitting. We used realistic biophysical parameters, and we have adjusted the scaling of the input in an arbitrary way such that the multistable regime in the model is in the flutter range.

↵c We note that there are two sources of noise in such spiking networks: the randomly arriving external Poissonian spike trains and the fluctuations caused by the finite size of the network. Here, we concentrate on finitesize effects attributable to the fact that the populations are described by a finite number N of neurons. In the meanfield framework (Mattia and Del Giudice, 2002, 2004), “incoherent” fluctuations caused by quenched randomness in the neurons' connectivity and/or by external input are already taken into account in the variance, and “coherent” fluctuations give rise to new phenomena. In fact, the number of spikes emitted in a time interval dt by the network is a Poisson variable with mean and variance Nν(t)dt. The estimate of ν(t) is then a stochastic process ν_{N}(t), well described in the limit of large Nν by ν_{N}(t) ≃ ν(t) + [ν(t)/N]^{1/2}γ(t), where γ(t) is Gaussian white noise with zero mean and unit variance, and ν(t) is the probability of emitting a spike per unit time in the infinite network. Such finiteN fluctuations, which affect the global activity ν_{N}, are coherently felt by all neurons in the network. This approach leads to the additive Gaussian noise corrections adopted in Equation 1.
 Correspondence should be addressed to Gustavo Deco, Department of Technology, Computational Neuroscience, Institució Catalana de Recerca i Estudis Avançats, Universitat Pompeu Fabra, Passeig de Circumvallació 8, 08003 Barcelona, Spain. gustavo.deco{at}upf.edu