|Home | About | Journals | Submit | Contact Us | Français|
Artificial sensations can be produced by direct brain stimulation of sensory areas through implanted microelectrodes, but the perceptual psychophysics of such artificial sensations are not well understood. Based on prior work in cortical stimulation, we hypothesized that perceived intensity of electrical stimulation may be explained by the population response of the neurons affected by the stimulus train. To explore this hypothesis, we modeled perceived intensity of a stimulation pulse train with a leaky neural integrator. We then conducted a series of two-alternative forced choice behavioral experiments in which we systematically tested the ability of rats to discriminate frequency, amplitude, and duration of electrical pulse trains delivered to the whisker barrel somatosensory cortex. We found that the model was able to predict the performance of the animals, supporting the notion that perceived intensity can be largely accounted for by spatiotemporal integration of the action potentials evoked by the stimulus train.
Microstimulation of the sensory regions of the cerebral cortex has been shown to mimic the transmission of normal sensory inputs by producing artificial sensory perceptions (Schafer 1888; Dobelle and Mladejov 1974; Tehovnik 1996; Rousche and Normann 1997; Romo et al. 1998; Talwar et al. 2002; Rousche et al. 2003; Troyk et al. 2003a, b; Tehovnik et al. 2004, 2005). Cortical stimulation paradigms can take advantage of the spatial topography of cortical representations, which underlie visual (Dobelle and Mladejov 1974; Bak et al. 1990; Schmidt et al. 1996; Troyk et al. 2003a, b; Tehovnik et al. 2004, 2005), somatosensory (Talwar et al. 2002; Rousche et al. 2003; Ohara et al. 2004), and auditory sensations (Rousche and Normann 1997; Rauschecker and Shannon 2002; Scheich and Breindl 2002; Rousche et al. 2003; Otto et al. 2005). The artificial sensations induced by electrical stimulation have been shown to elicit responses that are indicative of the mapped topographical location of the corresponding sensory inputs (Romo et al. 2000; Troyk et al. 2003b; Tehovnik et al. 2004, 2005).
The perceptual effects of modulating electrical stimulation delivered to a single location on the cortex have been explored to a lesser extent than the spatial mapping of the stimulation. Romo et al. (1998, 2000) have previously reported that a monkey was able to discriminate between pulse rates of different frequencies when one was delivered with electrical stimulation into the finger region of the somatosensory cortex and the other using a vibro-tactile stimulator attached to the corresponding finger. Stimulation pulse amplitude discrimination experiments by Rousche et al. (2003) conducted in the rat auditory cortex have also demonstrated the ability of a rat to perceive differences in the intensity of auditory-cortex stimulation. More recently, the discrimination of stimulus frequency of two electrical stimulation pulse trains delivered to the somatosensory cortex was explored by London et al. (2008), for the purpose of providing feedback for brain-machine interfaces. The results of their study showed that a monkey was able to discriminate frequencies of stimulus trains at a variety of cortical depths of electrode placement. The effects of varying the frequency, amplitude, and duration of cortical-stimulus trains have also been explored in the context of developing a human visual prosthesis. These subdural (Dobelle and Mladejov 1974; Dobelle et al. 1974) and intracortical (Schmidt et al. 1996) stimulation experiments in the human primary visual cortex suggest that adjusting the amplitude, duration, or frequency of stimulus trains can produce correspondingly brighter or dimmer percepts.
We hypothesized that the perceived intensity of the delivered stimulation pulse train can be largely explained by spatiotemporal integration of the action potentials evoked by the stimulus train, i.e. increased amplitude stimulates more cells, while increased frequency evokes more spikes over the stimulated population of neurons. To test this hypothesis, we modeled perceived intensity of a stimulus train by a leaky neural integrator based on prior studies of the spread of neural excitation due to electrical stimulation (Motz and Rattay 1986; Rattay 1989, 1990, 1999; Tehovnik 1996; Butovas and Schwarz 2003; Rattay et al. 2003). We then conducted two-alternative forced choice behavioral experiments adapted from the olfactory discrimination task of Uchida and Mainen (2003), in which we tested the ability of a rat to discriminate duration, frequency, and amplitude of electrical stimulus pulse trains delivered to the somatosensory cortex. We found that the performance of the animals, given a stimulus pulse train presentation was in agreement with the model predictions.
The model used to predict performance on the two-alternative forced choice discrimination task from electrical stimulation train parameters consists of two blocks (Fig. 1). The first block is used to map a stimulation pulse train to a perceived intensity. The second block maps the perceived intensity to a measure of performance on a psychophysical discrimination task.
The first block of the model computes the perceived intensity of a train of stimulation pulses delivered to the brain at an arbitrary frequency, pulse duration, and pulse intensity. This spike integrator model (SI) computes the perceived stimulus intensity using a leaky temporal integrator that sums evoked neural activity over the volume of stimulated tissue. This integrator can be conceptualized as a “container” into which neural activity is “deposited” following each stimulation pulse. The amount of activity deposited into the container following each pulse is directly proportional to the volume of tissue activated by the pulse. Once deposited, neural activity leaks out from the container at a rate specified by the time constant of integration.
Assuming that injected current radiates spherically away from tip of the stimulating electrode, the volume of the neural tissue that is stimulated by a single biphasic pulse, denoted Vs, is given by
where Va is the volume of a sphere (centered at the electrode tip) whose edges correspond to the furthest neurons that respond to the stimulation pulse, and Vx is the volume of a smaller concentric sphere of tissue that does not respond to the stimulus pulse because of damage or gliosis surrounding the electrode tip. If the radii of these two spheres are given by ra and rx, then Eq. 1 yields
Current spread away from the electrode tip obeys an inverse square law, so that the distance to the furthest responding neuron (ra) is proportional to the square root of the stimulation amplitude (Tehovnik 1996) and can be expressed as
where I is the stimulation amplitude of one phase of the stimulation pulse and C is a constant. Let Ia be the stimulation amplitude necessary to depolarize the region of radius ra and let Imin conceptually be the minimum stimulus amplitude necessary to get past the unstimulatable region (e.g., the electrode encapsulation of radius rx). Because the stimulation pulses are biphasic, we define Ia as the absolute value of the current-pulse train. Quantitatively, we define Imin as equal to the threshold of stimulation during each pulse presentation and equal to 0 between stimulation pulses. We can now rewrite Eq. 2 with respect to the stimulation amplitudes as
We can simplify Eq. 4 by accumulating all constants into one Ks = 4C3π/3, which yields the following expression for the volume of stimulated tissue Vs in terms of stimulation currents:
Equation 5 describes the volume of neural tissue that responds to each stimulation pulse. We assume that the stimulation pulse, which activates the neurons in Vs, will increment the activity level Ps that encodes the perceived intensity of the electrical stimulation. It is important to note that the pool of neurons which encodes Ps is not necessarily identical with the neurons residing within the volume Vs. This is because the neurons in Vs may send efferent projections to other populations of neurons that participate in encoding the perceived intensity of the electrical stimulation and because Vs may contain neurons that do not participate in the encoding of the sensation. Hence, Vs defines the amount of activity that is added to Ps following each stimulus pulse. The activity level Ps can be expressed as
In this equation, Ps represents the sustained activation of a neural population that decays with a time constant τ. We introduce a constant cpv in order to allow a change in dimensional units from the volume Vs to the perceived intensity Ps. For simplicity, we set the conversion factor cpv = 1. The same equation can be expressed in terms of the stimulation current as
In Eqs. 6 and 7, the time constant τ determines how long the perceived intensity of each stimulus pulse will last after each pulse. In monkey somatosensory cortical stimulation experiments, this time constant has been estimated to be 250–500 ms (Luna et al. 2005). The parameter Ks is a scale factor that regulates the growth of Ps in response to each stimulus pulse. For the simulations, we convert the differential Eq. 7 to a difference equation, given by
where i denotes the simulation time step. The sampling time interval dt in the simulations was chosen to match the duration of the biphasic pulse interval used in the rodent experiments, dt = 400 μs.
For simplicity, we set Ks = 1 because all of the simulations involved comparing the relative intensity of one stimulus train versus another. Because Eqs. 7 and 8 are linear, homogeneity property of linearity ensures that scaling the input would produce a proportionally scaled output. Therefore, Ks proportionally scales any perceived intensity.
When Ia < Imin, Ps produces negative values, for which there is no physical interpretation, as the stimulus train cannot be perceived if the delivered stimulus is composed of low amplitude pulses, characterized by Ia, which are below perceptual threshold, indicated by Imin. In this case, we set Ps = 0, as indicated in Eq. 8.
Perceived intensity increases each time the pulse is encountered and the function decays between the pulses (Fig. 2b). The function plateaus for a long duration stimulus, in which the added activation due to pulse presentation is equal to the decay between pulses (Fig. 2a). When we compare the perceived intensity for various pulse trains, we see a more rapid increase in the function for larger amplitudes and higher frequencies of stimulation. The same Ps can be obtained by manipulating amplitude, frequency, or duration of the stimulus train. As an example, we used the function to make a specific prediction that a 20-Hz, 130-μA/phase pulse train will produce the same perceived intensity as a 50-Hz, 70-μA/phase pulse train (Fig. 2c, see arrows).
In the model (Fig. 1), the output of the SI block (computed by Eq. 8) is passed as input to the behavioral decision block. The behavioral decision block maps the perceived intensity of brain stimulation onto a probability of generating a behavioral response to the stimulation in a simple two-choice discrimination task that we used in the rodent experiments. In this discrimination task, rats were trained to perform one response (denoted Rmax) to a high-intensity stimulus train and another response (denoted Rmin) to a low-intensity stimulus train. When a well-trained subject is presented with an arbitrary test stimulus, we denote the response to the test stimulus as Rtest. Over multiple presentations of the test stimulus, the fraction of responses, FHS, for which Rtest = Rmax, will be close to 1 if the test stimulus is perceived as similar to the high-intensity training stimulus (HS), and close to zero if the test stimulus is perceived as similar to the low-intensity stimulus. To model this behavioral choice process, the behavioral function FHS is represented by the following sigmoidal function:
The parameters a, b, c, and d in Eq. 9 allow the behavioral function to fit any dataset that has the sigmoidal characteristics. In the experiments, these parameters were chosen by fitting Eq. 9 to experimental data obtained from the animals during the frequency discrimination task. The response value R, used as an input into the equation, is the final value of perceived intensity R = Ps(tmax) for any given pulse-train stimulus.
In the present study, we tested the model predictions of performance of a rat on two-alternative forced choice experiments. The model parameters were obtained from each animal’s performance on frequency discrimination and threshold determination experiments.
The model predictions were experimentally tested by training rodents to perform a simple, two-choice stimulus-discrimination task, which was adapted from the olfactory discrimination task of Uchida and Mainen (2003). All surgical and behavioral procedures were approved by the UCLA Animal Research Committee in accordance with Federal regulations and guidelines for animal research.
Male Long-Evans Rats weighing 350–450 g were ordered from a commercial breeder and housed singly in a room with a 12-h on/off light–dark cycle. The rats were placed on a limited feeding schedule reducing them to 90% of the ad-lib weight, and they were maintained at this weight for the duration of the experiments. After reaching the target weight, they were surgically implanted with monopolar stimulating electrodes in the right hemisphere whisker barrel somatosensory cortex (S1BF). The electrode array was constructed using the technique described in (Williams et al. 1999), using 50 μm diameter Pt–Ir electrodes (California Fine Wire, Grover Beach, California, USA). The animal was deeply anesthetized using isoflurane, the skull exposed, and a 2 mm × 4 mm window was opened in the skull over the S1BF cortex. The exact location of the first electrode was based on stereotaxic coordinates originating at the bregma. In relation to the bregma, the location was 2.56 mm posterior, 4.2 mm lateral, and 2.5 mm ventral, which resulted in the electrodes being lowered to a depth of approximately 0.5 mm into the cortex. A small burr hole was drilled in the skull for insertion of a stainless steel return electrode. A set of 6–8 stainless-steel securing screws were attached to the skull, and bone cement was then used to seal the wound and securely mount the electrode connector plug on the skull.
For redundancy, the electrode array contained two rows of 7 electrodes each; however, in each animal we only stimulated the first electrode in the array for the duration of the experiments. The electrode was always stimulated in a monopolar fashion with respect to a distant electrode with large surface area.
The behavioral testing cage contained two infrared sensors, two feeding bowls, an ambient light source, and a tether connected to the headstage connector (Fig. 3). The Plexiglas cage measured 24 cm × 32 cm × 27 cm. It contained two feeding hoppers (open plastic boxes attached to the cage floor), each 9 cm × 7 cm × 2 cm. The experimental cage was positioned in an enclosed cabinet to reduce distractions during experiments. Two feeders (Med Associates, Georgia, Vermont, USA) positioned on top of the cabinet were controlled by a computer to deliver a 20-mg food pellet reward into the feeding bowls. Approaches to the feeding hoppers by the animal were recorded by the computer using two infra-red (IR) emitter–detector pairs. The detectors were positioned 3 cm above the cage floor to avoid inadvertent triggering by the tail. The computer communicated with the feeders and the IR detectors via a data acquisition card (Model DT300, Data Translation, Marlboro, Massachusetts, USA). A computer controlled house light and a camera were positioned above the cage, on the ceiling of the cabinet.
To send stimulation pulse trains, the head mounted plug was connected via a tether and a commutator to a neural stimulator and associated hardware (HR90K cochlear implant stimulator, Advanced Bionics Corp, Valencia, California, USA). The cochlear implant was used because it was readily available in our laboratory, provided the current-based stimulation characteristics, and had a ready-made computer interface. Custom software, written in Visual Basic 6 (Microsoft Corp. Redmond, Washington, USA) and Matlab (The MathWorks, Inc, Boston, Massachusetts, USA), was used to generate the stimulation pulse train parameters, control the feeders for behavioral reinforcement, and score the behavioral performance sensed by the IR detectors.
The pulse-timing for the stimulus pulse trains were generated using Matlab for each stimulus presentation. In all except for the last experiment, the biphasic pulses were generated periodically at 1/F time intervals for a given frequency F for the duration of the pulse-train presentation. In the last experiment, where we compared the responses to stimulus trains with randomized versus un-randomized inter-pulse intervals, a Poisson random generator was used to determine the duration of the randomized inter-pulse intervals for the given frequency, and the number of pulses presented were counted to deliver the same number as in the case of un-randomized stimulus presentation. In both cases, we used the time interval of 1 ms for the construction of the pulse trains.
In the experimental tasks, the rats were trained to obtain food from one feeder when they received a low frequency brain stimulation cue (LOW-STIM), and from the other feeder when they received a high-frequency brain stimulation cue (HIGH-STIM). For a given animal, the two cues were identical in amplitude and duration. Table 1 provides the training stimulus parameters used for each animal.
During training sessions, the animal was required to remain in the center of the box (without triggering the left or right IR detector) for a minimum of 0.5 s prior to the start of each discrimination trial. When either detector was triggered before 0.5 s, the 0.5 s timer was restarted. At the onset of the trial, either the LOW-STIM or the HIGH-STIM cue was delivered (the cue was selected randomly with equal probability on each trial).
The stimulus presentation was followed by a 1.5-s wait period during which the detectors were inactive. We noticed early on during the acquisition of the task that the animals were switching back and forth between the feeders during the first second or so from stimulus onset. To reduce the contribution of impulsive choices and increase stimulus control over choice behavior, we implemented a 1.5-s wait time after stimulus onset during which feeder activity was not reinforced. We began recording feeder activity beginning at 1.5 s after stimulus onset. Following this wait period, the detectors became active for a 10-s decision period. During this period, the first detector to be triggered was used to indicate the feeder selection. If the correct feeder was selected, then animal received a pellet reward in the chosen feeder. If the incorrect feeder was selected or the animal did not select either feeder for 10 s, the animal received no reward and the house light was turned off for a 3-s period. The next trial began when the animal had returned to the middle position in the cage so that neither detector was triggered for a minimum of 0.5 s. Each training session lasted between 1 and 2 h per day, depending on the level of cooperation of the test subject.
The pulse amplitude was initially set to 70 μA/phase for each subject, but if performance on the task remained at chance level (50%) for 7 days, then the amplitude was increased by ~100 μA/phase and another week of training was given. Once the training amplitude had been determined, the animal was considered to be fully trained when the success rate for choosing the correct feeder in response to the stimuli exceeded a criterion of 75% correct choices. After animals had been fully trained, discrimination experiments began. During the discrimination experiments, 75% of the trials were identical to the training trials designed to maintain the trained behavior, and the other 25% of the trials were non-rewarded test trials which provided the experimental data.
To determine threshold, the animals were trained to approach the HIGH-STIM feeder using the same training stimulus as used for the frequency, duration, and amplitude tasks (training parameters for the HIGH-STIM training stimuli are provided in Table 1 for each animal). The training method reinforced the association of the HIGH-STIM training stimulus with a reward pellet delivered to the HIGH-STIM feeder. The difference between this training paradigm and the one used for frequency, amplitude, and duration tasks was that in the training for threshold determination, we eliminated the LOW-STIM training stimulus. Instead, the LOW-STIM feeder delivered the pellet at random times unassociated with any stimulus. Thus, the training consisted of either delivering a pellet to the LOW-STIM feeder or equally likely, delivering a HIGH-STIM stimulus. As in the training used for frequency, duration and amplitude tasks, the animal had to remain in the center of the cage for 0.5 s before receiving the training HIGH-STIM training stimulus or receiving no stimulation (equivalent to receiving a LOW-STIM stimulus with 0 μA amplitude). After the 0.5-s interval, we implemented a 1.5-s wait period during which the animal had to decide on the choice of the feeder. We began recording feeder activity beginning at 1.5 s after stimulus onset. Following this wait period, the detectors became active for a 10-s decision period. During this period, the first detector to be triggered was used to indicate the feeder selection. Reinforcement pellet was delivered to the HIGH-STIM feeder if the animal approached the HIGH-STIM feeder after receiving the HIGH-STIM stimulus. If the animal instead approached the LOW-STIM feeder after receiving the HIGH-STIM stimulus, the house light turned off for 3 s. Alternatively, if the HIGH-STIM stimulus was not delivered and the animal approached the HIGH-STIM feeder, the light would be turned off for 3 s.
Using this training method, the animals learned to approach the HIGH-STIM feeder only when the HIGH-STIM training stimulus was delivered. Otherwise they were more likely to approach the LOW-STIM feeder instead. The animals learned to get into the middle of the cage, wait approximately 2 s if they did not receive the HIGH-STIM stimulus, and then approach the LOW-STIM feeder because it was likely to deliver a pellet. However, if, while it was in the middle of the cage, the HIGH-STIM stimulus was delivered, the animal learned that it had to approach the HIGH-STIM feeder in order to get a reward.
After animals achieved success rate of 75% of approaching the HIGH-STIM feeder, when they received the HIGH-STIM stimulus, threshold determination tests began. In these tests, the pellet-reinforced training trials (either a HIGH-STIM stimulus, or equally likely a pellet delivered to the LOW-STIM feeder) were delivered in 75% of the trials to maintain task performance, and the other 25% of the trials were non-rewarded test trials, in which the test stimuli were delivered. The test stimuli consisted of the same frequency and duration as the HIGH-STIM training stimulus, but the amplitude was varied from 0 μA/phase to that of the HIGH-STIM training stimulus in 12 evenly spaced increments. The amplitudes for the test stimuli were selected at random from the pool of the possible amplitudes without replacement and then the pool was reinstated after all amplitudes had been used for stimulation and the amplitudes were again chosen from the original set.
This testing procedure resulted in low probability for the animal to approach the HIGH-STIM feeder when it did not perceive the test stimulus and high probability when it was presented with a test stimulus that was similar in its intensity to the HIGH-STIM stimulus. As the amplitude of the test stimuli increased from 0 μA/phase, the probability that the animal would approach the HIGH-STIM feeder remained low for lower amplitude stimuli and then rose toward 1 for the amplitudes above threshold.
To quantify finding the threshold using this methodology, we first found the stimulus with the lowest amplitude from the discrete set of tested amplitudes which resulted in probability of the animal approaching HIGH-STIM feeder that was greater than two standard deviations above the probability of approaching HIGH-STIM feeder, given a 0 μA/phase test stimulus. We then approximated the stimulation threshold to be just below that amplitude from the set of possible amplitudes. The standard deviations were estimated from the binomial distribution statistics of the responses to the test stimuli.
During these non-rewarded test trials, the animal was presented with test stimuli of different frequencies, while other electrical stimulation parameters were held constant. By interleaving the test trials within many standard training trials, it was possible to maintain the performance accuracy on the task while probing the effect of varying the electrical stimulation parameters on the perceived intensity of the stimulation. The LOW-STIM and HIGH-STIM training trials were used to keep a real-time estimate of the trained performance during the experiment. The test trials were collected and analyzed after the experimental data was collected.
During each test trial, the animal was presented with stimulus of the same amplitude and duration as the HIGH-STIM and the LOW-STIM training stimuli, and characterized by a novel stimulus frequency, selected randomly without replacement from a fixed pool of test frequencies (10, 20, 30, 40, 50, 60, 70, and 80 Hz). We selected these frequencies because they spanned the range of physiological spike rates of cortical neurons. Once the entire set of test stimuli had been used for stimulation, the set was reinstated and the test stimuli again started to be chosen without replacement from the original set. Each daily test session was terminated when the subject had received over 20 iterations of each stimulus in the pool, if the subject’s performance on the trained trials decreased or if the subject had been in the apparatus for over 2 h. The testing resumed the following day. The number of accumulated trials of each test stimulus is provided in column labeled frequency in Table 1 for each animal. Several days were sometimes necessary to accumulate ~20 iterations for stimuli of each frequency for an animal because the total number of test trials for 20 iterations of each frequency is 8 frequencies * 20 trials/frequency = 160 test trials.
The performance of the animal in response to a test stimulus of any given frequency was quantified by observing how often the animal would approach the HIGH-STIM feeder when the test stimulus of that frequency was delivered.
We used the model to predict the performance on a task which tested the ability of the animal to determine the perceived intensity of electrical stimulation trains which differed in pulse train duration. Duration testing was conducted in a way that was similar to the frequency testing experiments, except that instead of varying the frequency of the test stimuli, we kept the frequency and amplitude constant and varied the duration of the test stimulus pulse trains. The durations used for the test stimuli were selected from a pool of discrete values from 0.1 to 1.1 s in 0.1-s intervals. The amplitude and frequency for these test stimuli were kept at or near the HIGH-STIM training stimulus level for each animal. The frequencies and amplitudes were 40 Hz, 70 μA/phase; 70 Hz, 300 μA/phase; 80 Hz, 500 μA/phase; and 70 Hz, 300 μA/phase for animals 1–4, respectively. The number of trials accumulated for each test stimulus is provided for each animal in Table 1 in the duration column.
We tested the ability of the model to predict the performance of the animal on a task in which the amplitude of the test stimulus pulses was varied while the frequency and duration were held constant. This task was conducted in a way that was similar to frequency and duration testing, but instead of varying frequency or duration, we varied the amplitudes of the test stimuli, while keeping all other parameters constant. The duration was kept at 1 s for all test stimuli.
For the amplitude testing experiments, we had to select constant frequency test stimuli with amplitudes that spanned a range, which according to the model prediction would produce the behavior associated with LOW-STIM stimulus for the low amplitudes and would produce the behavior associated with HIGH-STIM stimulus for the higher amplitudes. Because animals 2, 3, and 4 had thresholds which were high compared with animal 1, we used the HIGH-STIM training frequency for the test stimulus trains for these animals. Using this pulse frequency meant that we could find lower amplitudes for the test stimuli that, according to the model prediction, would be generalized as the LOW-STIM stimulus. Alternatively, because animal 1 had low stimulation threshold, we used the LOW-STIM training frequency and higher amplitudes for the test stimuli to obtain the experimental performance associated with the HIGH-STIM stimulus. Because we hypothesized that the model would predict the behavioral performance on this experiment, we used the model to predict the range of amplitudes for the test trials for each animal.
The SI model is a leaky integrator. It predicts that the perceived intensity of brain stimulation depends only upon the accumulated activity within the integrator and nothing else. One implication of this prediction is that for pulse trains with inter pulse intervals that are significantly shorter than the time constant of integration (several hundred milliseconds), the perceived intensity of brain stimulation should depend only upon the mean and not the variance of the inter pulse intervals. We obtained the model prediction for the performance of animal 1 on a frequency testing experiment, in which the test stimuli were composed of 2-s periodic pulse trains with 0 coefficient of variance (CV = 0) and pulse trains with randomized inter-pulse intervals (CV = 1) of the same frequencies, amplitude, and duration. To test the model predictions, we conducted the corresponding experiment, in which we varied the CV for test pulse trains with varied frequency for animal 1. In these experiments, the test stimuli varied in the pulse frequency from 10 to 80 Hz in 10 Hz intervals and in whether or not the inter-pulse intervals were randomized (Poisson distributed with CV = 1) or fixed (constant interpulse interval of 1/frequency with CV = 0) (Koch and Segev 1998). The amplitude of these test stimuli was the same as used for the HIGH-STIM and LOW-STIM training stimuli. The duration used for these test stimuli was 2 s.
The parameters that describe the behavior of the SI perceived intensity model, defined in Eq. 8, consist of the stimulation threshold Imin and the time constant τ. Parameters a, b, c, and d in Eq. 9, describe the sigmoid of the behavioral response function. We experimentally found the threshold of stimulation Imin as discussed earlier. The SI parameter τ and the behavioral model parameters a, b, c, and d were obtained by finding best least-squares function fit to the experimentally determined performance on the frequency testing experiment. The initial values for the unconstrained nonlinear optimization routine implemented in Matlab were: a = 0.1, b = 0.01, c = 20, d = 0.5, and τ = 0.25 s for all animals.
Once the complete model was determined using the threshold determination experiment and fitting the model curve to the frequency testing experiment, the model for each animal was used to predict that animal’s performance on the duration and amplitude testing experiments. These experiments test the perceived intensity of electrical stimuli, which vary in duration while keeping the amplitude and frequency of the stimulus trains constant and test amplitude, while keeping the duration and frequency of test stimulus trains constant.
Additionally, we tested the sensitivity of the model predictions to perturbations of the SI model from the optimally derived time constant τ and from the theoretically derived exponent. For each perturbation to the SI model, we obtained the behavioral function parameters by determining an optimal model fit to the performance on the frequency testing experiment, and then we used the model with the complete set of parameters to predict the animal performance on the duration and amplitude testing experiments.
The experimental data was collected over several days for each animal on all experiments, and because there were a large number of testing stimuli for all of the experiments for each animal, we were able to obtain a limited number of trials for each test stimulus (~20 trials per stimulus) as shown in Table 1. For each test stimulus, we used the fraction of trips to the HIGH-STIM feeder when that stimulus was presented as a measure of performance on all of the experimental tasks. This measure provided us with the ratio of the trips the animal made to the HIGH-STIM feeder to the total number of times the particular stimulus was delivered to the animal. Because all experiments were two-alternative-forced-choice, we used the standard deviation of the binomial distribution statistic in order to obtain a measure of uncertainty in the experimental data. This measure is defined as:
Using this standard uncertainty measure in the context of our experiments, n corresponds to the total number of test trials of the particular stimulus and p corresponds to fraction of trips to the HIGH-STIM feeder. This measure of standard deviation decreases with the increased number of trials for any given test stimulus. The error bars in each plot that shows experimental performance correspond to this measure of experimental uncertainty.
In order to investigate the ability of the SI model to predict perceived intensity of an electrical stimulus pulse train delivered to the somatosensory cortex, we first determined the parameters for each animal’s model using the frequency testing experiment and the threshold determination experiment and then used the model with the complete parameters to predict the performance on the duration testing experiment and the amplitude testing experiment.
The training for the animals for the frequency, duration, and amplitude tasks as well as for the threshold determination task resulted in similar training success for both HIGH-STIM and LOW-STIM stimulus presentations for all animals. The individual training performance for each animal is shown in Table 2. During the frequency, duration, and amplitude tasks averaged across all subjects, when the HIGH-STIM stimulus was presented, the animals approached the HIGH-STIM feeder in 82 ± 10% of the training trials and the LOW-STIM feeder in 19 ± 10% of the training trials. When the LOW-STIM stimulus was presented, the animals approached the LOW-STIM feeder in 79 ± 11% of the training trials and the HIGH-STIM feeder in 21 ± 11% of the training trials. During the training for the threshold determination experiment, the animals learned to respond to the HIGH-STIM stimulus presentation by approaching the HIGH-STIM feeder in 75 ± 15% of the training trials and the LOW-STIM feeder in 15 ± 15% of the training trials. In the training trials in which the animals received no stimulation, they approached the LOW-STIM feeder in 70 ± 5% of the training trials and the HIGH-STIM feeder in 30 ± 5% of the training trials.
We obtained the values for Imin using the threshold determination experiment. The values for the thresholds for each animal were determined experimentally to be 10, 225, 248, and 125 μA/phase for animals 1–4, respectively. The standard deviations on the testing trials for the animals on the threshold determination experiments were 9.2 ± 3.7%.
To investigate the ability of the animal to perceive frequency of electrical stimulation train delivered to the somatosensory cortex and to obtain the parameters for the perceptual model (Fig. 1), we conducted a frequency testing experiment. Figure 4 plots the fraction of trips that the animal made to the HIGH-STIM feeder when stimulated with the test stimulus of each possible frequency. Because the fraction of trips to the HIGH-STIM feeder was a single ratio of the number of trips to that feeder to the total number of stimulus presentations of a given frequency, the error bars were calculated using standard deviation of the binomial distribution of the performance for each frequency of the test stimuli. The number of test trials of each frequency is shown in Table 1.
Based on the frequency testing experiment results shown in Fig. 4, the animals consistently generalized progressively higher frequency stimulus trains as being more similar to the HIGH-STIM stimulus associated with the HIGH-STIM feeder. The parameters for the HIGH-STIM and LOW-STIM training pulse trains are shown above the plot for each animal. The amplitude and duration of the test pulse trains are inset in the upper left corner of each plot and the frequencies of the test stimuli are indicated along the x-axis of each plot. For example, for animal 1, there were 8 test stimulus trains presented to the animal: 70 μA/phase, 1 s, 10 Hz; 70 μA/phase, 1 s, 20 Hz; 70 μA/phase, 1 s, 30 Hz; …; 70 μA/phase, 1 s, 80 Hz. The only value that changes for the presented test stimuli in each plot is the frequency of the pulse train that is plotted along the x-axis.
We used the results from the frequency testing experiment shown in Fig. 4 to obtain the parameters needed for the complete neural stimulation perceptual model (Fig. 1). The model is specific to each animal and contains two free parameters (Imin and τ) in the Spike Integrator part of the model (Eq. 8), and four free parameters (a, b, c, and d) in the behavioral decision part of the model (Eq. 9). Imin is the threshold current pulse amplitude obtained from the threshold determination experiment. We used unconstrained nonlinear optimization to obtain values for τ, a, b, c, and d that resulted in FHS values which fit in the least-squares sense each animal’s performance on the frequency task. The optimization procedure used the same initial conditions to fit the performance on the frequency testing experiment for each animal. The model fit of the animal’s performance on the frequency task is plotted in dashed lines on the plots in Fig. 4. The parameters which resulted from the optimization are shown in the lower right of each plot. Because the output of the behavioral model corresponds to a particular SI model output, we indicated the output of the SI model (R = Ps(tend)) along the x-axis below the corresponding frequency values. Animal 1 had lower R values because it had lower stimulus thresholds than the other animals. Since Ps increases proportionally to the 3/2 power of the stimulus current (Eq. 7), the higher the stimulus current used to evoke a response, the higher the values of Ps, given the similar values for the time constant τ. Importantly, the optimal values of the time constant in the duration discrimination experiment were on the order of hundreds of milliseconds, in agreement with previous studies that have measured the time constant of integration using model-weighted spike trains for the accumulated stimulus intensity of pulse trains delivered to the somatosensory cortex in monkeys (Luna et al. 2005).
Once the complete model parameters were determined for each animal, we used the model for each animal to predict their performance on two other experiments. In the duration testing experiment, we varied the duration of test stimulus trains, while the frequency and amplitude remained constant. In the amplitude testing experiment, we varied the amplitude of the pulses of the test stimulus trains, while the frequency and duration remained constant. For these two experiments, the HIGH-STIM and LOW-STIM reinforced training stimuli remained the same as for the frequency testing experiment. Solid black lines in Fig. 5A, B show the model predictions for the individual animal’s performance on the duration and amplitude tasks, respectively. The model parameters for each animal used in the model predictions are shown above the column of the two plots for each animal. The SI model outputs are shown below the x-axis values corresponding to each test stimulus. The test stimuli were determined by two constant parameters indicated in the upper left of each plot and the value of the third parameter which varied as indicated along the x-axis. In Fig. 5A, the durations of test stimuli are indicated along the x-axis, while the amplitude and frequency were kept constant for each test stimulus as indicted in the upper left of each plot. In Fig. 5B, the amplitudes of test stimuli are varied as indicated along the x-axis, while the duration and frequency were kept constant for each test stimulus as indicted in the upper left of each plot.
The performance of the animals on each of the two experiments is plotted in Fig. 5A, B on the same plots as the model predictions. The model predictions agreement with the experimental data is assessed using mean squared error between the animal’s performance and the model prediction of the performance given the same stimuli.
We can gain the perspective on the relationship between the perceived intensity of the test stimuli and the perceived intensity of the HIGH-STIM training stimulus by illustrating the experimental comparison of the perception of two stimuli, the perceived intensities of which were compared in Fig. 2c and the “Materials and methods” section. While every point on the solid line in each panel of Fig. 5A, B represents the model prediction of the animal selecting HIGH-STIM feeder given a test stimulus, the vertical and horizontal dashed lines in Fig. 5B plot for animal 1 correspond to just one such prediction exemplified in Fig. 2C and in the “Materials and methods” section. That specific theoretical example indicated that a 50-Hz, 70-μA/phase, 1-s stimulus train will have the same perceived intensity (same SI model output) as a 20-Hz, 130-μA/phase, 1-s train. The HIGH-STIM stimulus used for animal 1 is the same as the first stimulus in the theoretical example, 50 Hz, 70 μA/phase, 1 s. As we discussed in the “Materials and methods” section, the animal’s behavior was maintained at approximately 75% of the trips to the HIGH-STIM feeder, given the HIGH-STIM training stimulus. In the particular case of animal 1, this accuracy was measured at 67% during the amplitude testing experiment. The dashed horizontal line in Fig. 5B for animal 1 indicates this performance when the HIGH-STIM training stimulus is delivered. If the animal perceives the 20-Hz, 130-μA/phase, 1-s test stimulus train as being the same intensity as the HIGH-STIM training stimulus (which was 50 Hz, 70 μA/phase, 1 s) as predicted by the example then we would expect the same fraction of trips to the HIGH-STIM feeder whenever this test stimulus is presented as when the HIGH-STIM training stimulus is presented. For the amplitude testing experiment for animal 1, we chose 20-Hz, 1-s test stimuli, with varied amplitude from 70 to 170 μA/phase. We can see the fraction of trips the animal made to the HIGH-STIM feeder for each one of these test stimuli as the amplitude is varied along the x-axis in Fig. 5B for animal 1. The one example test stimulus which we focus on here is the 130 μA/phase test stimulus, since it corresponds to the discussed prediction that the perceived intensity of that stimulus will be the same as that of the HIGH-STIM stimulus. The vertical dashed line in Fig. 5B plot for animal 1 corresponds to the 20-Hz, 130-μA/phase, 1-s test stimulus. The experimental results indicate that the animal made approximately 65% of the trips to the HIGH-STIM feeder when this test stimulus was presented. The nearly identical fraction of trips to the HIGH-STIM feeder (67 and 65%) when either the 50-Hz, 70-μA/phase, 1-s training stimulus or the 20-Hz, 130-μA/phase, 1-s test stimulus was delivered respectively suggests the high degree of similarity of perceived intensity of the two stimuli. We would like to point out that the 65% of trips to the HIGH-STIM feeder should be compared not to 50% chance performance. If chance level is defined as the number of times the animal will approach HIGH-STIM feeder when no stimulation is delivered, it will be similar to the fraction of times the animal will approach the HIGH-STIM feeder when LOW-STIM stimulus is delivered (22% of the training trials).
The particular theoretical example presented in Fig. 2c is straightforward in its ability to relate to the experimental results because this example compares two stimuli which should result in identical perceived intensity, and therefore identical behavioral response on the amplitude testing experiment to the HIGH-STIM training stimulus. All 20-Hz, 1-s test stimuli characterized by the other amplitudes along the x-axis do not result in an identical perceived intensity as that of the HIGH-STIM training stimulus and therefore the behavioral responses to these test stimuli are not expected to be the same as the behavioral response to the HIGH-STIM training stimulus. However, given our hypothesis for the model, one would expect the animals to behave the same (make the same fraction of trips to the HIGH-STIM feeder) given the same perceived intensity independent of whether the stimuli that resulted in these perceived intensities differed in duration, amplitude, or frequency. This is the reason we included the perceived intensity values along the x-axis of each of the plots in Figs. 4 and and5.5. For each animal, the model should predict the same fraction of trips to the HIGH-STIM feeder for the same perceived intensity as indicated for any test stimulus in Figs. 4 and and55.
In order to establish the relevance of the derived form of the model on its ability to predict animal behavior, we tested the sensitivity of model prediction to changing the time constant τ from the optimally obtained value of 0.48 s for animal 1 and then changing the derived 3/2 exponent of the SI model.
We compared the predictive ability of three SI models, which varied only by values of τ. We used 1.1, 0.1 s, and the optimally derived value of 0.48 s for τ. For each model, we found the best fit behavioral function such that we obtained least-squares fit to the frequency testing experiment (Fig. 6A, left plot). The model fitting optimized the functions well enough to overlay them on top of each other on the plot. With the three sets of behavioral function parameters corresponding to each SI model, we then used the three complete models to predict behavior on the duration and amplitude testing experiments (Fig. 6A, right plots). We observed that the changes in the time constant caused substantial changes to the ability to predict behavior on the duration testing experiment and to a lesser extent on the amplitude testing experiment.
As shown in the SI model derivation (Eqs. 7 and 8) in the “Materials and methods” section, the model uses 3/2 exponent that was obtained based on the hypothesis that the number of neurons that depolarize in response to a current pulse are proportional to a sphere surrounding the electrode (power of 3) and that diameter of that sphere is proportional to the square root of the pulse amplitude. To examine if the derived 3/2 exponent in the SI model is important in obtaining a more accurate model prediction for the perceived intensity of a stimulation train, we conducted the model predictions for the animal 1 amplitude testing experiment performance in which we changed the exponent used in the SI model, while keeping the model parameters τ and threshold Imin constant. For each of the three models, varied only by their exponents (1, 3/2, and 2), we first found the behavior model parameters (a, b, c, and d) using the best-fit approximation of the frequency testing experimental data, and then we used the complete model for each of the three cases to obtain predictions for the amplitude and duration testing performance. Figure 6B on the left plots the best-fit model along with the frequency testing experimental data. The best-fit models in each of the three cases overlay each other on the frequency testing plot, indicating that the behavioral model could capture the performance on the frequency testing experiment equally well using the three possible SI models. We then obtained the predictions using each of the three models now complete with the behavioral model parameters, for the animal performance on the duration and the amplitude testing experiments. As shown in the Fig. 6B on the right duration and amplitude testing experiments, the models varied in their ability to accurately predict the performance on the amplitude task, with the model that uses the 3/2 exponent being the most accurate of the three.
By changing the time constant and then the exponent in the SI model, we observed that the perturbation to the SI model which predicts the perceived intensity of a stimulus train delivered to the cortex resulted in a decreased ability to predict the experimental behavior. This observation suggests that the derived form of the perceived intensity SI model is at least partially responsible for obtaining more accurate predictions of the animal performance on perceived intensity experiments which we conducted.
We conducted a further experiment, in which we varied the coefficient of variance (CV) for the frequency discrimination test pulse trains for animal 1. In agreement with the model predictions, the results indicate that the performance of the subject did not depend on the coefficient of variation of the pulse train (Fig. 7). These results agree with prior findings (Romo et al. 1998) on the ability of the monkey to discriminate frequencies of cortical stimulation. Additionally, the results are consistent with the prior findings which suggest that firing rate rather than periodicity account for the ability to differentiate frequency of stimulation (Romo et al. 1998; Hernandez et al. 2000; Salinas et al. 2000; Luna et al. 2005).
All other parameters kept equal, perceived intensity of a stimulus train increases with higher frequency, higher amplitude, or longer duration of the stimulus train. We presented a model which suggests that the perceived intensity of electrical stimulation of the somatosensory cortex can be predicted quantitatively by a leaky integrator circuit that sums individual stimulation pulses over time and incorporates a first-order forgetting factor. The model parameters were obtained on the frequency discrimination and threshold determination task, and tested on amplitude and duration discrimination tasks. Consistent with the model predictions, the animals generalized higher amplitude and longer duration test stimuli as being similar to the high-frequency training stimuli and low amplitude and short duration test stimuli as being similar to the low-frequency training stimuli.
We used the highest possible frequency (corresponding to HIGH-STIM stimulus) for test stimuli in determining the threshold of stimulation, Imin. According to the model and experimental results, perceived intensity is a trade-off between frequency and amplitude of the stimulus. Therefore, the lowest possible perceived intensity can be obtained using either high frequency and low amplitude stimulation or high amplitude and low frequency stimulation. Because our goal in determining the value of the Imin parameter for the model was to find the lowest amplitude which in theory would affect the fewest number of neurons in the stimulated region, we used the highest possible frequency, so that this minimal amplitude could be detected.
The modeling results in which we perturbed the theoretically derived form of the perceived intensity SI model, produced less accurate predictions of task performance compared to the predictions obtained using the derived form of the SI model. This suggests that the derivation of the SI model is at least partially responsible for the ability of the model to predict the behavior.
Even though the biophysical model described here can be used to make predictions of perceived intensity of the stimulus, the experiments make no direct claim as to the stimulation mechanism on which the model is based. However, previous work in this field lends support to the assumptions used for the model creation. The model is based on the assumption of homogeneous activation of the neural tissue surrounding the electrode and on the assumption that the activated cells will fire in phase with the pulse presentations.
Electrophysiological work by Butovas and Schwartz (2003) further provides experimental support to both the assumption of homogeneous activation as well as the assumption that the activated neurons respond in phase with the presented stimulation pulses. In this work, the authors describe the effects of electrical stimulation pulses delivered to an electrode positioned in the rat somatosensory cortex. The recorded extracellular action potentials were recorded from the cells surrounding this electrode at a variety of cortical depths and radial distances from the stimulating electrode. The orientation and depth of the stimulating electrode had no effect on the evoked responses. This investigation describes a uniform behavior from 80 to 100% of the cells within the stimulated region (up to 1.3 mm from the electrode tip depending on the stimulus amplitude). If a single pulse was presented, the neurons responded after 2–3 ms with one or two spikes, followed by 100–150 ms of inactivity. When the pulses were delivered in form of a pulse train at pulse rates greater than 10 Hz (in particular, they show data for 20 and 40 Hz), the cell responses locked to the presented stimulus and depolarized in response to each pulse with inactivity between the pulses. At the edges of the stimulated region, the number of the responding cells reduced rapidly with increased distance from the electrode. These results suggest that the activation of the cortical tissue in response to presentation of current pulse trains is homogeneous within the stimulated region and that the activated neurons respond in phase with the presented stimulation pulses against the background of inactivity.
While the neural excitation for individual cortical neurons varies by an order of magnitude (100–4000 μA/mm2) (Tehovnik 1996), Butovas and Schwartz make an argument that the response of the neurons that they observed in the vicinity of the electrode is a result of the preferential stimulation of the fibers of passage near the electrode tip (Rattay 1999; McIntyre and Grill 2000). The activation spreads out from the electrode tip by the large number of randomly oriented fibers near the electrode of both the inhibitory and excitatory neurons. The activation of the inhibitory neurons could largely be responsible for reducing the activity between the pulse presentations and phase-locking the neural response into the presented stimulus pulses.
The radius of this activated region has been shown in numerous studies to be proportional to the square root of the stimulation current (Tehovnik 1996) and the behavior of the cells outside this uniform activation region would be dominated by the differences in excitation threshold. In his review article, Tehovnik found the uniform measure of the number of neurons stimulated within this activation region:
These results suggest that the activation of the cortical tissue in response to presentation of current pulse trains is uniform within the stimulated region and that the activated neurons respond in phase with the presented stimulation pulses against the background of inactivity.
According to the model and the experimental results, high amplitude and low-frequency stimulus train is interpreted behaviorally the same as a lower amplitude, but higher frequency stimulus train. High amplitude stimulates a spatially wide area of the cortex, while low amplitude stimulates a spatially narrow area of the cortex. Frequency of the pulse train conveys the “stimulus intensity” to the neurons that get stimulated. Based on this interpretation, our results suggest that a spatially wide stimulus (analogous to high amplitude) of weak intensity (analogous to low frequency) will be interpreted the same as a spatially narrow stimulus (analogous to low amplitude) of stronger intensity (analogous to high frequency). This description is consistent with human psychophysics studies (Lederman 1974), as well as the results of Pruett et al. (2001), who discovered that a monkey could not discriminate between wide groove-width gratings at low force with narrow groove-width gratings at high force on a psychophysics experiment. In this context, it is feasible that the perception of an electrical stimulus train is consistent with the physiological mechanism which interprets natural somatosensory sensations.
The quantitative model for the perceived intensity of somatosensory-cortex stimulation may prove useful for the development of technologies which would mimic normal sensory processing by electrical stimulation. Such technologies hold great promise for the development of brain–machine interfaces for assisting patients with brain or spinal cord injuries, as well as amputees. For example, artificial stimulation of somatosensory cortex (similar to that used in our present study) might be useful for endowing limb prostheses with artificial somatosensory receptors, so that patients can use their prosthetics not only for movement but also to “touch” objects. The model could provide an intensity function that may be used to calibrate stimulation parameters against the gain of an encoded signal that is delivered to the brain.
In addition to benefiting the development of technologies for artificial sensation, the model and experimental procedures we have introduced here may provide useful tools for further investigation of the neural mechanisms that underlie normal sensory processing. For example, our findings as well as previous findings (Romo et al. 1998; Luna et al. 2005) demonstrate that the brain integrates somatosensory stimulation with a time constant on the order of several hundred milliseconds, but the neural substrate for this slow integration process is presently unknown. The two-choice discrimination task that we used to test the model predictions in this study may provide a useful experimental tool for further investigation of this question. Figure 6A in particular suggests that in order to investigate the physiological nature of the time constant, one should consider specifically the comparison of stimuli that vary in duration. Since it has been hypothesized that NMDA receptors, which have a slow time constant of decay, may participate in slow sensory integration (Wong and Wang 2006), it may be possible to manipulate the value of the integration time constant in the discrimination tasks by delivering NMDA agonists or antagonists to regions off the brain that participate in the sensory-integration process. Such future experiments may facilitate the development of a biologically plausible neural network model of the slow sensory integration process, which would be constrained by the quantitative model we have presented here.
We thank C. Della Santina, and N. Davidovics for helpful comments on this manuscript and S. Badelt for technical assistance and discussion. We thank M. Faltys and Advanced Bionics for providing the use of the cochlear stimulation hardware and software. Dr. Fridman’s current address is at Johns Hopkins University, Ross Bldg. Rm 830, 720 Rutland Ave., Baltimore, MD 21205. This work was supported by UCLA Core Faculty Grant 2005–2007 and NIH NRSA Fellowship F32 DC009917.
Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.