A computational model that predicts behavioral sensitivity to intracortical microstimulation

Objective. Intracortical microstimulation (ICMS) is a powerful tool to investigate the neural mechanisms of perception and can be used to restore sensation for patients who have lost it. While sensitivity to ICMS has previously been characterized, no systematic framework has been developed to summarize the detectability of individual ICMS pulse trains or the discriminability of pairs of pulse trains. Approach. We develop a simple simulation that describes the responses of a population of neurons to a train of electrical pulses delivered through a microelectrode. We then perform an ideal observer analysis on the simulated population responses to predict the behavioral performance of non-human primates in ICMS detection and discrimination tasks. Main results. Our computational model can predict behavioral performance across a wide range of stimulation conditions with high accuracy (R2 = 0.97) and generalizes to novel ICMS pulse trains that were not used to fit its parameters. Furthermore, the model provides a theoretical basis for the finding that amplitude discrimination based on ICMS violates Weber’s law. Significance. The model can be used to characterize the sensitivity to ICMS across the range of perceptible and safe stimulation regimes. As such, it will be a useful tool for both neuroscience and neuroprosthetics.


Introduction
Intracortical microstimulation (ICMS) has proven to be a useful tool to explore the neural mechanisms of behavior [1] and has yielded important insights into visual [2][3][4], auditory [5,6], and tactile [7][8][9][10] processing. As ICMS elicits meaningful and repeatable percepts, it may also be used to restore sensation for patients who have lost it [8,9,[11][12][13][14][15][16][17]. Sensitivity to ICMS has been studied by measuring the detectability and discriminability of different regimes of electrical stimulation applied to different cortical regions [1,8,12,[18][19][20][21][22][23][24]. These previous studies provide important insights on how changing different parameters affects the resulting perceptual experience. However, quantitative predictions about the detectability or discriminability cannot be made based on these previous studies, given how sparsely the parameter space is sampled. To fill this gap, we sought to develop a model that predicts the sensitivity to ICMS whose parameters span the range that has been deemed safe [25]. To this end, we developed a simple simulation of ICMS-evoked activation in a neuronal population and then used ideal observer analysis to predict behavioral sensitivity from the simulated population response. The model was developed based on behavioral data we had previously collected that characterize the detectability and discriminability of ICMS applied to the primary somatosensory cortex (S1) of macaques and spanning a wide range of stimulation parameters [23]. We find that the simple model can accurately predict sensitivity to ICMS across all conditions tested.

Behavioral methods
The methods are described in detail in a previous paper [23] and are only summarized here.
Animals and implants. Procedures were approved by the University of Chicago Institutional Animal Care and Use Committee. Each of two male Rhesus macaques (Macaca mulatta, 6 years of age, around 10 kg in weight) was implanted with a Utah electrode array (UEA; Blackrock Microsystems) in the hand representation of area 1 of the right hemisphere. After about 2 years of data collection, the array of one of the animals failed so another array was implanted in area 1 of its left hemisphere. Each UEA consists of 96 electrodes with 1.5 mm long shanks, spaced 400 μm apart, and spanning a 4×4 mm area. The electrode tips are coated with a sputtered Iridium Oxide film using a standard process [26,27]. The electrode shaft is insulated with parylene-C along its length, with the exception of the tip, which has a targeted exposure length of 50 μm. Electrode impedances were measured to be between 10 and 80 kΩ prior to implantation. We mapped the receptive field of each electrode by identifying which areas of skin evoked multiunit activity (monitored through speakers). Stimulation pulses were delivered using a 96-channel constant-current neurostimulator (CereStim R96, Blackrock Microsystems).
Behavioral tasks. Two animals were trained to perform two variants of a two-alternative forced choice task: detection and discrimination. Animals were seated at the experimental table facing a monitor, which signaled the trial progression. Each trial comprised two successive stimulus intervals, each 1 s and separated by a 1 s inter-stimulus interval. The animals responded by making saccadic eye movements to a left or right target presented on the visual monitor and correct responses were rewarded with water or juice. The animals were first trained with mechanical stimuli, consisting of mechanical indentation of the skin with varying depth (from 0 to 2000 μm), until performance stabilized. Their task was to report which stimulus interval contained the only mechanical indentation (detection) or the mechanical indentation of greater depth (discrimination). The animals then performed the same detection and discrimination tasks based on ICMS of S1 (see below for a detailed description of each behavioral task). Unless otherwise specified, ICMS consisted of a train of symmetric biphasic pulses with cathodal phase leading and an interphase interval of 53 μs. Performance was very similar across animals, as previously documented [8,23], so we combined results from the two animals for the purposes of the present modeling effort.
Detection. On each trial, ICMS was delivered in one of two consecutive stimulus intervals and the animal's task was to indicate which stimulus interval contained the stimulus. In different experimental blocks, we investigated the effects of pulse width (from 50 to 400 μs), pulse train frequency (from 50 to 1000 Hz), and pulse train duration (from 2 to 1000 ms) on the animals' ability to detect stimuli varying in amplitude. We also studied the interaction between frequency and duration by varying these two parameters while keeping the pulse amplitude fixed at 40 μA (chosen because it was near threshold) as well as the interaction between pulse width and frequency using the range of amplitudes. Each manipulation was repeated with multiple electrodes to gauge the generality of our results.
Amplitude discrimination. On each trial, two ICMS pulse trains were presented sequentially and the animals' task was to indicate which of the two was stronger. One of the two stimuli was a standard stimulus, with amplitude of 70 μA (approximately half way between threshold and maximum amplitude), and the other was a comparison stimulus whose amplitude varied from trial to trial. We varied ICMS frequency (from 50 to 500 Hz) or pulse train duration (from 50 to 1000 ms). In those conditions, the frequency or duration was the same for both intervals but varied randomly from trial to trial. predicting behavioral responses based on these simulated neuronal responses (figure 1). To simulate the responses of an individual neuron, we compute its probability of firing in response to each ICMS pulse [28][29][30], taking into consideration both the properties of the stimulation pulse and the neuron's recent spiking history (the latter to incorporate absolute and relative refractoriness). The firing probabilities are then summed over all ICMS pulses to determine that neuron's response. To simulate the population response, the current waveform is described as decaying in amplitude as it spreads away from the stimulation site. The behavioral response is determined by comparing the simulated neuronal response to the ICMS pulse train delivered in each stimulus interval using ideal observer analysis. The model comprises 11 free parameters, which are optimized to fit the behavioral data [23]. The model thus provides a succinct description of a complex behavioral data set using a simple description of the process that culminates in the measured behavior.
Neuronal simulation. Fluctuations of the membrane potentials at the nodes of Ranvier of myelinated fibers have been shown to be approximately normally distributed [31]. Accordingly, the threshold current, I th , which bridges the gap between resting membrane potential and threshold (assuming a constant threshold, see [28]), is also distributed normally, so the firing probability of a neuron can be modeled as Bernoulli random variable, p, driven by current I c as follows: and the probability that a spike will be evoked by a given ICMS pulse is given by [30]: where I c is the ICMS current, G is a gain parameter, and I th50 is the mean threshold current, that is, the current that causes the neuron to fire with a probability of 50%. Because the standard deviation of the membrane potential increases approximately linearly with its mean, the standard deviation of the threshold current is also a linear function of its mean σ=sI th50 [28,31,32], where s is the 'relative spread.' After each spike, the threshold current I th50 jumps up then decays exponentially along a time course that tracks absolute and relative refractoriness [28]: where t abs represents an absolute refractory period (fixed at 1 ms) during which the firing probability equals zero (see equations (2), (3)); γ and τ ref reflect the spike-triggered conductance change and its decay time constant, respectively. The effect of pulse width PW on threshold I 0 is described according to the following empirical relationship [30,33]: where I 0,PW=∞ is the threshold current for long pulses, equivalent to the rheobase current, and C is the chronaxie, which is the time required for ICMS that is double the strength of the rheobase current to stimulate a neuron and is necessary to incorporate the effects of pulse width on sensitivity.
Using equations (2)-(4), then, we computed the firing probability of each simulated neuron in response to a train of ICMS pulses. After each spike, the membrane potential returned to its resting state, the threshold increased [34,35] then decayed back to its resting level after the absolute refractory period (see equation (3)), resulting in a decrease in firing probability during the refractory period (equation (2)). After computing the firing probability of the first pulse using equation (2) (with I th50 =I 0 ), the model recursively calculates the firing probability for the rest of pulses.
The number of spikes evoked in each neuron by a pulse train is given by the sum of firing probabilities across pulses. We then calculate the number of spikes generated by neurons across the volume of simulated tissue, over which the ICMS current decays at a rate of 1/r 2 , where r is distance from the electrode tip. We used 21 neurons from nearest (r=1) to farthest (r=3) with an increment Δr=0.1, having established that the results were relatively insensitive to the number (or density) of simulated neurons. The firing probability and corresponding spike counts decrease as distance increases and converges to zero (figure 1).
Because neurons in the stimulated volume fired independently (that is, were not synaptically connected) in this simplified model, spikes were spatially and temporally integrated to calculate the effective number of spikes, R p evoked in the neuronal population (equation (5)). To incorporate the fact that neuronal signals are not integrated indefinitely, but rather over tens or hundreds of milliseconds [36,37], in making a perceptual decision, we used an exponentially decaying time window, w n , whose time course was a parameter.
where p(n, r) denotes the firing probability of a neuron located at distance r by the nth pulse, N is the number of ICMS pulses, and w n is an exponentially decaying function with time constant τ i . The term 4πr 2 computes from the neuronal responses evoked along radius r the responses evoked over a volume consisting of a series of 21 concentric spherical shells of radius r.
Given that the spike count follows the binomial distribution, the variance of the population response is given by: In summary, there are 11 free parameters: I 0,PW=∞ (rheobase current), C (Chronaxie), γ (conductance change), τ ref (time constant of refractoriness), s (proportionality constant to relate variability of the input current to its amplitude), τ i (time constant of the integration time window) (equation (5)), and 5 gain parameters (3 for detection and 2 for discrimination (G in equation (2) and table 1). The gain parameter was fit to each protocol separately to account for observed differences in sensitivity across conditions, likely caused by sensory adaptation and learning (see below). Note, however, that the gain parameter changes the overall sensitivity to ICMS but does not affect the relative sensitivity to different ICMS pulse trains. Parameters were optimized to minimize the root mean squared error (RMSE) between actual and predicted performance (see below) for the detection and discrimination tasks. For this, we used a combination of Matlab function fmincon with multiple initial points and simulated annealing [38].
Ideal observer analysis. From equations (5) and (6), the population response evoked by stimulus 1 on any given trial is given by: where n p1 is the random deviation (i.e., noise) of the actual response from the mean response. The ideal observer analysis assumes that stimulus 1 (the comparison stimulus) will be perceived as more intense than stimulus 2 (the standard stimulus) if r p1 >r p2 , the probability of which is given by: where P stands for the cumulative distribution function of the random variable indicated in the subscript. Since the population response is the sum of independent responses of individual neurons, it is normally distributed, according to the central limit theorem. The response variance was computed from equation (6) where Φ is the cumulative Gaussian function with zero mean and unit variance. In signal detection theory, the argument of Φ is the discriminability, d p ¢ [39]. On detection trials, the responses on a blank interval, R p2 , was calculated from equation (2) with I c = 0, yielding a value near zero..

Results and discussion
Behavioral data As described previously in more detail [23], detection performance improved as ICMS amplitude and duration increased, regardless of what other parameters were manipulated (figures 2(A)-(C)). However, the psychometric functions and resulting detection thresholds were also   The difference in performance across these paired conditions cannot be accounted for based on the small difference in frequency and reflects differences in sensitivity across these different protocols. These differences may be due to different states of sensory adaptation or different levels of training (the two were confounded in these data sets).
dependent on the other parameters [40]. In brief, thresholds decreased with increases in pulse width and pulse frequency, both effects tending to level off (∼200 μs for pulse width and ∼250 Hz for frequency). As might be expected, discrimination performance improved as the difference in amplitude between the two stimuli increased (figures 2(D), (E)). Unlike its detection counterpart, however, discrimination performance was relatively insensitive to ICMS frequency; discrimination performance improved with increased stimulus duration up to about 200 ms and leveled off, a phenomenon that is also observed with the discrimination of mechanical stimuli applied to the skin [41]. Finally, performance varied across experiments for the same set of electrodes and identical stimulation regimes (figure 3). For example, sensitivity was higher for experimental blocks where pulse width (detection) or stimulus duration (discrimination) was manipulated than for blocks where frequency was manipulated. The ICMS parameters are almost the same except for a small difference in frequency (250 Hz versus 300 Hz), which is unlikely to underlie the observed difference. To account for this variability, the model included a gain parameter that was allowed to vary across protocols (equation (2), table 1(B)), with higher gain denoting higher sensitivity.

Neuronal simulation
The goal was to develop a model that predicted performance on the ICMS detection and discrimination tasks. The model comprised two components (figure 1): the first was a simulation of the neuronal activity evoked by each ICMS pulse train; the second was an ideal observer analysis that simulated the perceptual decision of which of two intervals was selected on a given detection or discrimination trial. In brief, each simulated neuron was subjected to the ICMS waveform, with amplitude decaying as an inverse function of the square of the neuron's distance from the electrode tip [42]. Thus, simulated neurons that were farther away from the electrode tip generated fewer spikes ( figure 4(A)). As expected, larger ICMS currents induced more spikes from individual neurons and from the neuronal population. Specifically, the simulated response increased with both pulse width ( figure 4(B)) and pulse frequency ( figure 4(C)). Of course, the population spike count increased with ICMS durations, but the effect of stimulus duration on detection or discrimination performance was limited by the integration time window [36,37]. The variance of the population response was well approximated using a power function ( figure 4(D)). The discriminability was calculated from simulated neuronal responses and used to predict detection and discrimination performance using ideal observer analysis. Circles and triangles indicate detection trials and discrimination performance, respectively. Different colors denote different experimental conditions (see figures 3 and 5). The black line shows the predicted probability (see equation (9)).

Predicting behavioral responses from simulated neuronal responses
As described above, we simulated the population response to each ICMS pulse train then computed the behavioral response using ideal observer analysis (equation (9), figure 4(E)). The parameters of the model were fit so that the resulting predictions matched the large behavioral data set, which consisted of 159 different ICMS combinations, 87 from detection and 72 from discrimination. The parameters of the neuronal model (table 1(A)) were identical across protocols except for the gain parameter (table 1(B)), which was allowed to vary across protocols to account for the varying sensitivity to the same (or similar) ICMS pulse trains ( figure 3). Once model parameters were optimized to fit the behavioral data, we constructed psychometric functions for the pulse width and pulse frequency manipulations in the detection experiments by simulating the population response for ICMS trains whose amplitude spanned the range from 10 to 90 μA in 1 μA steps and computing the resulting detectability ( figures 5(A), (B)).
For the duration/frequency manipulation ( figure 5(C)), the number of ICMS pulses was increased in increments of 1. For the pulse frequency and duration manipulations in the discrimination experiments, we simulated the responses of comparison stimuli ranging in amplitude from 20 μA to 100 μA and compared these to each of the fixed standards (figures 5(D), (E)).
This simple model could predict both detection and discrimination performance with high accuracy (R 2 =0.97, RMSE=0.052, figure 5(F)). That is, the dependence of the simulated population response on stimulation parameters mirrored the dependence of the behavioral sensitivity to these same parameters. In other words, the model can accurately predict the degree to which a given ICMS pulse is detectable or the degree to which a pair of ICMS pulses are discriminable across all conditions tested. The model is able to predict behavior across a wide range of conditions even though it does not take into consideration neuronal morphology [43,44], the biophysics of the stimulation and of neuronal activation [45], and the active spread of neuronal activity through synaptic connections [43,44,46,47].

Model validation
We validated the model by testing it on behavioral data sets that were not used to fit its parameters [23]. One of the main differences between model-fitting and validation datasets was that the polarity of the pulses was reversed in the validation set (from cathodal to anodal phase leading). One validation dataset consisted of behavioral performance on a detection task, with varying pulse width, pulse train frequency, or pulse train duration. The other consisted of performance on a discrimination task with two different reference amplitudes (at 30 and 100 μA). These datasets consisted of 100 different conditions, 90 from the detection experiments and 10 from the discrimination experiments. All the model parameters were fixed, except for the gain parameter, which was fit to individual behavioral sets to compensate for differences in overall sensitivity, especially necessary given the well-documented phenomenon that anodal phase leading pulses lead to lower sensitivity [21,48,49]. The model accounted for 91% in the variance of the validation set (figure 5(G)) (R 2 =0.91, 0.89, 0.79 for the three detection protocols, respectively, and 0.99 for the discrimination protocol). Thus, the model generalizes to a data set that was not used to obtain its parameters (except for the G, which was refit).
We also wished to assess whether a simpler possible model, one in which detection and discrimination predictions are based solely on charge delivery, could predict behavioral performance. To this end, we computed the amount of charge delivered during each stimulus and fit a logistic function to detection performance as a function of charge delivered and to discrimination performance as a function of difference in charge delivered. This model comprised two free parameters for each experiment condition for a total of 10 parameters. Model performance accurately predicted the effects of phase duration on detection and those of pulse train duration on discrimination but fared very poorly in all other conditions (see table S1). Thus, the present model outperforms a marginally simpler model, suggesting that its (relatively low) level of complexity is warranted.

Comparison of simulated ICMS-evoked responses to measured mechanically evoked responses
Having established that the computational model can accurately predict behavioral performance, we compared the simulated responses to ICMS with measured neuronal responses to mechanical stimulation. To this end, we used neurophysiological recordings from S1 obtained while the animals performed detection and discrimination tasks based on mechanical stimulation of the skin that were analogous to those described above for ICMS. From these data, we calculated the effective spike count using the same integration time window used for the simulated responses. For the measured responses to mechanical stimulation, we found that the rate-intensity was a negatively accelerating function of stimulus amplitude and could be well approximated using a power function with exponent less than 1 (α≈0.3) (figure 6(A)); furthermore, variance increased slightly supralinearly with firing rate (β≈1.3) ( figure 6(B)). In contrast, simulated neuronal responses to ICMS exhibited the opposite pattern: The rate-intensity function was accelerating (α≈2.6) (figure 6(A)) and the variance-rate function decelerating (β≈0.9) (figure 6(B)). The shallow variance rate function reflects the highly synchronized population responses to ICMS [50]. According to this model, which accurately captures the behavioral performance of the animal, the neuronal activation evoked by ICMS has different properties than does its counterpart evoked by mechanical stimulation of the skin. These differences in rate-intensity and variance-rate functions across modes of activation might account for the fact that mechanical discrimination follows Weber's law whereas electrical stimulation does not [23].

Model parameters
The optimized parameters for the neuronal model are comparable to equivalent quantities measured in cortical neurons (table 1(A)). In visual cortex, measured rheobase currents, I 0, PW=∞ , span a wide range, from 1 to 500 μA [51] as do chronaxies, C, from 0.03 to 31 ms, depending on the cell type [48,52]. These two parameters are required to capture the dependence of sensitivity on pulse width (figures 2(A), 5(A)). The fitted relative spread, s (0.25), is comparable to that measured in auditory neurons of cats (which ranges from 0.05 to 0.4) [29,53]. The two parameters modulating relative refractoriness, γ and τ ref (2.32 and 112 ms respectively), are in the same order of magnitude as those used in previous models (γ=0.97 in [28][29][30] and τ ref =100 ms in [54]). These two parameters were optimized to capture the saturation of sensitivity at high frequencies shown in figures 2(B)-(D), 5(B)-(D)) Lastly, the time constant of the integrating window, τ i (40 ms), is comparable to that derived in previous studies (5-8 ms in [37], 12-22 ms in [36]) and this parameter is required to capture the effect of the pulse train duration (figures 2(C), (E), 5(C), (E)).
We also tested how sensitive the predictions of the model were to the value of each of its parameters (figure S1). We found that model performance was highly sensitive to the two time constant parameters, τ ref and τ i , but much less sensitive to parameters such as rheobase current, I 0,PW=∞ and spiketriggered conductance change, γ.

Comparison with a previous model
In a previous model [24], the performance of rodents on an ICMS discrimination task was predicted based on a leaky integrator of delivered charge. A logistic function (with 4 free parameters) was fit to the data relating discrimination performance to integrated charge, so variability in the underlying percepts was not explicitly represented in the model. The model could predict the animals' ability to distinguish ICMS pulse trains varying in amplitude, frequency, or duration. Importantly, however, a different logistic function was fit to each condition. Unlike the present model, then, a single set of parameters could not account for discrimination across all conditions. Thus, while the study concluded that changes in ICMS amplitude, frequency, and duration all have consequences on perceived magnitude, the relative contributions of these to discriminability could not be inferred from the model. The present model fills that gap by accounting for the effects of all stimulation parameters with a single set of parameters (in addition to the variable gain). It is important to note that the gain parameter would be unnecessary if all stimulation conditions had been interleaved within each experimental block.

Model limitations
First, simulated neurons are excited by current that is passively diffusing away from the electrode tip, which is at odds with the observation that ICMS excites primarily axons rather than cell bodies [55]. Second, neurons in the model are not interconnected even though circuit-level dynamics may play a role in shaping the detectability and discriminability of ICMS pulse trains. Third, the gain parameter needs to be adjusted to the animal's sensitivity. In practice, this parameter could be estimated based on a few threshold measurements and the resulting model could then be used to predict a much larger data set. Note also that, while the gain parameter determines the overall sensitivity to ICMS, it does not affect the relative sensitivity to different ICMS trains. Fourth, the model does not take into account the neuronal adaptation-the progressive desensitization to ICMS over the course of an experimental block-caused by prolonged stimulation. Indeed, adaptation may have in part underlied the requirement for different gain parameters for different conditions, each with different overall levels of stimulation (the other cause for different gains being learning).

Conclusions
Despite its simplicity, the proposed model accurately captures detection and discrimination behavior across the behaviorally relevant range of stimulation conditions. The model can thus be an invaluable tool to develop ICMS regimes with known psychometric properties without time-consuming preliminary psychophysical experiments. Furthermore, results from the neuronal simulation suggest a reason why ICMS amplitude discrimination violates Weber's law.