## Abstract

During coordinated eye–hand movements, saccade reaction times (SRTs) and reach reaction times (RRTs) are correlated in humans and monkeys. Reaction times (RTs) measure the degree of movement preparation and can correlate with movement speed and accuracy. However, RTs can also reflect effector nonspecific influences, such as motivation and arousal. We use a combination of behavioral psychophysics and computational modeling to identify plausible mechanisms for correlations in SRTs and RRTs. To disambiguate nonspecific mechanisms from mechanisms specific to movement coordination, we introduce a dual-task paradigm in which a reach and a saccade are cued with a stimulus onset asynchrony (SOA). We then develop several variants of integrate-to-threshold models of RT, which postulate that responses are initiated when the neural activity encoding effector-specific movement preparation reaches a threshold. The integrator models formalize hypotheses about RT correlations and make predictions for how each RT should vary with SOA. To test these hypotheses, we trained three monkeys to perform the eye–hand SOA task and analyzed their SRTs and RRTs. In all three subjects, RT correlations decreased with increasing SOA duration. Additionally, mean SRT decreased with decreasing SOA, revealing facilitation of saccades with simultaneous reaches, as predicted by the model. These results are not consistent with the predictions of the models with common modulation or common input but are compatible with the predictions of a model with mutual excitation between two effector-specific integrators. We propose that RT correlations are not simply attributable to motivation and arousal and are a signature of coordination.

## Introduction

Coordinating eye and arm movements is central to our natural behavior. Eye–hand coordination depends on a combination of retinal and extra-retinal signals necessary for accurate movement (Crawford et al., 2004). Saccade reaction times (SRTs) and reach reaction times (RRTs) are highly variable and are correlated from trial to trial with a variable degree of correlation (Herman and Maulucci, 1981; Lünenburger et al., 2000; Boucher et al., 2007b) (for an example, see Fig. 1*A*). Because saccade durations are stereotyped (Bahill et al., 1975), one possible advantage of correlated reaction times (RTs) is to allow the eye to acquire a target at a consistent time with respect to the arm movement. This temporal consistency could enhance the performance of visual processes guiding manual accuracy (Neggers and Bekkering, 2002). Therefore, RT correlations may constitute a signature of coordinated eye and arm movements.

RT correlations, however, are not necessarily a measure of behavioral coordination. They occur in behavioral paradigms such as dual-task paradigms designed to tax cognitive processing and perceptual paradigms with time-consuming sensory processing (Pashler, 1994; Tombu and Jolicoeur, 2002). RT correlations may also result from global changes in responsiveness caused by changes in arousal or motivation (Boucher et al., 2007b). Global changes modulate motor preparation and could lead to RT covariations independently of sensorimotor coordination. Correlations could also be attributable to bottlenecks in processing if longer processing of the first process delays the second process (Pashler, 1984; Navon and Miller, 2002).

Dual-task paradigms can be applied to looking and reaching to test the link between RT correlations and coordination. We can cue movements to be made together and coordinated or to be made separately at different, unpredictable times by introducing a random delay, or stimulus onset asynchrony (SOA), between the cue onsets. Varying the SOA will manipulate coordination and test the relationship between RT correlations and sensorimotor coordination.

Computational models allow us to explicitly distinguish between different hypotheses for the possible mechanism underlying RT correlations. Many modeling approaches of RTs are based on the integrator paradigm, in which information is accumulated until a behavioral response is initiated (Laming, 1968; Luce, 1986). The diffusion model, a type of stochastic integrator model in which activity grows continuously in time (Ratcliff, 1978), can account for the chronometry of RT tasks with single (Smith, 1995) and simultaneous cues (Diederich, 1995), and in two-choice RT tasks (Usher and McClelland, 2001; Eckhoff et al., 2008). The model is also consistent with the neuronal activity during detection (Hanes and Schall, 1996), perceptual discrimination (Roitman and Shadlen, 2002; Gold and Shadlen, 2007), and other sensorimotor tasks (Cisek, 2007) and can be implemented in biophysically plausible network models (Wong and Wang, 2006), whose properties can be related to diffusion model parameters (Roxin and Ledberg, 2008).

Here, we combine computational models with behavioral experiments. We extend the diffusion model to dual-task paradigms using two interacting nonlinear leaky integrators that encode the preparation of saccades and reaches. Model predictions are tested against behavioral data from three monkeys performing an SOA task that involves looking and reaching.

## Materials and Methods

### Training and surgical procedures

Monkeys made visually guided saccades and reaches for juice rewards. Reaches were made to a touch-sensitive screen (ELO Touch Systems) mounted in front of either a liquid crystal display or a custom light-emitting diode (LED) board on which targets were presented. Eye position was monitored with an optical video eye tracker (ISCAN). Behavior was controlled using custom Labview (National Instruments) code running on a real-time PXI platform. After several weeks of training, a head-restraint prosthesis was implanted to maintain stable head position during recordings. All procedures were done under isoflurane anesthesia and sterile conditions. At least 3 weeks after surgery, the monkeys began training in eye movement and reaching tasks. All surgical and animal care procedures were done in accordance with the National Institute of Health guidelines and were approved by the New York University Animal Care and Use Committee.

### Behavioral tasks

Behavioral data were collected from three adult male rhesus macaques (monkey H, monkey J, and monkey S). In the SOA task, each trial started with the illumination of a red and green central target (Fig. 1*B*). Monkeys were required to fixate and touch the screen within 4° of the stimulus, after which a yellow peripheral location cue appeared 20° to the left or right of center. When the location cue appeared, the red center fixation point was simultaneously extinguished, cuing a saccade to the peripheral cue. At a randomized interval after the saccade go cue, the monkey received an auditory cue to make a reach to the same peripheral target. The speaker for the auditory cue was located behind the touch screen and in the same location for all trials and thus did not contain information about the target location on each trial. Each trial was initiated when a monkey touched a pair of proximity sensors placed near the body with his hands. The non-reaching hand was required to maintain touch on a proximity sensor throughout the trial. A brief auditory tone preceding reward delivery served as a secondary reinforcer on all correct trials.

The SOA task requires animals to dissociate the timing of the saccade from the reach. However, when each monkey was initially presented with the saccade and then reach cues, they tended to perform a reach and saccade together and did not dissociate the movements from each other. In pilot experiments, we determined that we could train the monkeys to follow the cues more accurately by using an uneven distribution of the duration of the interval between go cues, SOAs, and by interleaving different task conditions. We obtained the best results when SOA was 0 in 70–90% of the trials each day; in the remaining 10–30% of trials, the SOA was a random number drawn from a uniform distribution with range 0–620 ms (Fig. 1*C*). To discourage monkeys from always preparing reaches as soon as a location cue appeared, we interleaved SOA task trials with delayed saccade and touch task trials in which the location cue was red and signaled a saccade alone. In those trials, subjects were rewarded after a saccade to the target while continuing to touch the center green stimulus. Delayed saccade and touch trials constituted 10–20% of total trials each day and accounted for 11% of successful trials for monkey H, 19% of successful trials for monkey J, and 18% of successful trials for monkey S.

The full set of trials, therefore, included a set of trials in which the reach was made quickly together with the saccade (SOA of 0 ms), a set of trials in which the reach was never made (delayed saccade and touch), and a set of trials in which the reach movement was made at varying delays with respect to the saccade (SOA between 0 and 620 ms).

To further discourage early or late movements, each monkey was trained to wait until the go signal for each movement before promptly making that movement. To discourage late movements, monkeys had to respond to the appropriate go signal with a saccade within 350 ms and a reach within 600 ms (700 ms for monkey S). To keep animals from trying to predict the go signal and moving early, a trial was aborted if the SRT was shorter than 120 ms or if the RRT was shorter than 150 ms. These constraints on maximum and minimum RT encouraged monkeys to prepare each movement immediately after the relevant cue.

To control for the influence of an auditory go cue for the reach compared with a visual go cue for the reach, we also ran a separate set of trials with a modified SOA task in one animal (monkey H). In the modified, visual–visual SOA task, we delivered a visual go cue for the saccade and a visual go cue for the reach. All other details of the visual–visual SOA task matched the SOA task described above.

### Data collection and basic analysis

Eye position and touch position on the screen were sampled at 1 kHz. Each signal was time stamped and streamed to disk along with data about each trial from the Labview behavioral control program. The time of the visual cue presentation was recorded as the time at which a photosensor detected a simultaneous stimulus change on the monitor or LED board. Reach start was recorded as the time at which the touch was no longer detected within a 4° window around the central touch point. Eye position was filtered and instantaneous velocity calculated. Saccade start was defined as the time at which velocity first exceeded 50°/s at the beginning of the movement after the saccade go cue that resulted in peak velocity >200°/s.

### Behavioral analysis

Reach and saccade RTs were binned according to the SOA using 13 bins of 50 ms from 0 to 650 ms. The trial counts per bin were 19,567, 268, 899, 754, 940, 770, 577, 542, 954, 1158, 941, 720, and 478 for monkey H, 1155, 406, 298, 262, 307, 274, 301, 283, 330, 396, 352, 236, and 81 for monkey J, and 16,242, 333, 312, 355, 335, 354, 362, 372, 274, 310, 317, 261, and 98) for monkey S.

Overlap in movement preparation for each trial was computed as the elapsed time between the reach go cue and the saccade start (Fig. 1*C*). Thus, overlap depends on both the SOA and the monkey's SRT on that trial. If the reach go and saccade go cue are given simultaneously in a trial, the overlap is the SRT. Because the reach go cue always occurred with or after the saccade go cue, the largest overlap value is the longest SRT when SOA was zero. The overlap is negative when the reach go is given after the saccade begins. Because SRT varies from trial to trial, a given SOA duration is associated with a distribution of overlap durations. We tested SOA intervals uniformly distributed over more than two times the range of SRT so that, when we grouped trials by overlap, the distribution of SRTs for a given overlap condition was not distorted.

To graph the effect of overlap duration on the RT correlation, trials were binned by overlap time using 9 bins of 50 ms from −250 to 200 ms. The resulting trial counts per bin were 931, 679, 558, 542, 1401, 567, 691, 1561, and 7691 for monkey H, 397, 343, 259, 395, 211, 285, 289, 546, and 3104 for monkey J, and 298, 317, 336, 494, 242, 291, 405, 1064, and 9374 for monkey S. Bins with <10 trials of data were not subject to additional analysis.

The correlation between SRT and RRT for each group of trials was computed using Pearson's correlation coefficient. The 95% confidence intervals were computed as tanh(arctanh(*R*) ± 1.96 × SE), where *R* is the correlation coefficient, and SE = 1/*N* being the number of trials in the bin (Zar, 1996).

### Computational models

#### Leaky nonlinear integrators for saccade and reaches

In the computational models, preparation and initiation for saccades and reaches is encoded in the activity of two neuronal populations, *s* and *r* (Fig. 2*A*). A particular movement is initiated when the mean firing rate of the corresponding population reaches a threshold (Fig. 2*B*). The mean firing rate of each population obeys a differential equation of the type
where τ is the recruitment time constant of the neuronal population, of the order of 100 ms, ϕ(·) is an input–output transfer function, and *I _{i}* is the net input received by the integrator, which consists of a sum of external and recurrent inputs induced by the coupling within and between populations (Wilson and Cowan, 1972; Amit and Tsodyks, 1992; Pinto et al., 1996). The first term on the right side of Equation 1 accounts for the decay of network activity in the absence of inputs. The function ϕ(

*I*) is the steady firing rate with which the population responds when driven by a constant input

*I*and is in general a nonlinear and monotonically increasing function of the input. A parsimonious and mathematically convenient choice for ϕ(

*I*) is a threshold-linear function: where

*g*> 0 is a gain modulation factor, and θ is a threshold related to the rheobase current of neurons. When inputs are below threshold, the transfer function is 0 and the population activity

*r*(

_{i}*t*) decays to zero with time constant τ. When they are above threshold,

*r*(

_{i}*t*) grows at a rate dictated by the magnitude of the input

*I*, also with characteristic time constant τ. In its simplest form, the input received by each population is a sum of recurrent inputs, proportional to the mean firing activity of the population and population-specific external inputs: where α >0 is the strength of the recurrent coupling, and

_{i}*E*(

_{i}*t*) is the movement-specific external input, which provides the available moment-by-moment information relevant for initiation of movement

*i*. This information is brought about by bottom-up sensory signals as well as top-down urgency signals and may vary over the course of the trial. We lump together these two contributions and model the net external input as a single time-dependent signal

*e*(

_{i}*t*) perturbed by additive noise. The noise reflects instantaneous fluctuations in the form of neuronal and synaptic noise, both along the sensory pathway as extrinsic from the incoming signal, as well as intrinsic noise attributable to the spontaneous background activity from other brain areas. The external input has therefore the form where ς is the noise intensity, and η

*(*

_{i}*t*) is a source of Gaussian noise with white spectrum, that is, the values of the noise are uncorrelated at different times. Except for the model with shared fluctuations, described below, the noise sources η

*(*

_{s}*t*) and η

*(*

_{r}*t*) are uncorrelated with one another. We choose a very simple form for the time-dependent movement-specific signal

*e*(

_{i}*t*), which switches from 0 to 1 at the time onset of the go cue for movement

*i*. This sudden increase in the input at the go cue is sufficient to drive the unit from the subthreshold regimen,

*I*< θ, where the unit does not integrate inputs, to the suprathreshold regimen,

_{i}*I*> θ, where it does. The dynamical equation of the integrator is obtained by plugging Equations 2–4 into the differential Equation 1. In the suprathreshold regimen, the resulting equation has the form with

_{i}*k*≡ (1 −

*g*α)/τ,

*e′*≡

_{i}(t)*g*([

*e*

_{i}(t)−θ)/τ > 0, and ς′ ≡

*g*ς. Equation 5 describes a stochastic integration of the signal

*e′*in which the rate of accumulation depends linearly on the accumulation variable [Ornstein-Uhlenbeck process (Gardiner, 1985)]. The value of the coefficient

_{i}(t)*k*reflects the net balance between the positive feedback generated by the gain and the recurrent excitation on the one hand, and the negative feedback attributable to the leakage on the other. When

*k*> 0, negative feedback dominates and the activity of the accumulator grows until eventually approaching a steady state determined by the signal strength

*e′*. Conversely, when

_{i}(t)*k*< 0, positive feedback is dominant and the firing activity blows up exponentially unless some mechanism for threshold detection comes into play. In either case, the characteristic time constant of the dynamics is 1/

*k*. The particular case

*k*= 0 corresponds to an infinitely long exponential decay (or growth) and so to a lossless integration of the external input

*I*=

_{i}*e′*+ς′

_{i}(t)_{i}η

*. The integrated variable*

_{i}(t)*r*(

_{i}*t*) follows in such case a Brownian trajectory with a constant drift toward the threshold that is proportional to the signal strength

*e′*, like the classic diffusion model.

_{i}(t)The parameters are chosen so that the integrator remains silent in the absence of stimulation, i.e., when *e _{s}* =

*e*= 0. At the onset of the go cue signal for movement

_{r}*i*, at

*t*=

*T*

_{onset}

*, the deterministic component of the input*

_{i}*e*switches from 0 to 1, and, as a result, the activity of

_{i}*i*th unit starts increasing. This rise in activity persists until the activity of the unit reaches a threshold value

*H*at some time

*t*=

*T*

_{cross i}. After movement initiation, the sensory signal for a particular movement is switched off when the corresponding unit has reached threshold (Fig. 2

*B*), that is, The total RT is the sum of a constant residual time

*T*

_{0}, which represents the net effect of delays attributable to transduction and transmission of signals, and the time of integration to threshold: RT =

*T*

_{0}+

*T*

_{cross i}−

*T*

_{onset i}. The duration of the integration process,

*T*

_{cross i}−

*T*

_{onset i}, is a random variable by virtue of the noise present in the sensory signal

*E*(

_{i}*t*).

#### Models of interaction between integrators

##### Gain modulation model.

Effects of slow, trial-by-trial changes in excitability, such as might occur as a result of arousal or motivation, were modeled with a gain *g* that was randomly drawn on each trial from a normal distribution with mean 1 and variance ς_{g}^{2} ≪ 1 (Fig. 3, left). The value of the gain was identical in both integrators and was fixed throughout the trial. Note that changes in gain are not independent from changes in other parameters because modifying the gain is equivalent to changing appropriately the self-excitation α, the SD ς, and the threshold of the input–output function, θ, as can be seen in the definitions following Equation 5.

##### Common noise model.

The condition that fluctuations in the inputs be uncorrelated among units is relaxed in the common noise model (Fig. 3, middle). In that model, the inputs received by the saccade and the reach units at any given time have correlated variabilities of magnitude *c*ς^{2}, where *c* is the correlation coefficient of the two inputs and satisfies 0 ≤ *c* ≤ 1. The covariance matrix of the inputs is then nondiagonal and can be written as
where angle brackets denote average over trials, δ(*t* − *t*′) is the impulse function, and ρ* _{ij}* = 1 if

*i*=

*j*, whereas ρ

*=*

_{ij}*c*if

*i*≠

*j*. The impulse function on the right side of this equation expresses the fact that fluctuations at different times are uncorrelated. Fluctuations with this covariance matrix are generated from a sum of a shared noise source,η

*(*

_{c}*t*), and an unit-specific noise source, ξ

*(*

_{i}*t*), both with zero mean and unit variance: The scaling in

*c*of the coefficients ensures that the variance of the net fluctuations remains 1 as we vary the correlation coefficient.

##### Shared signal model.

Forward excitation is implemented by adding a fraction *f* of the signal to the counterpart unit (Fig. 3, right). The net signals received by both units are then
where *e _{i}*(

*t*), with

*i*=

*r*,

*s*, is given by Equation 6 and 0 <

*f*< 1.

##### Mutual excitation model.

The integrators can be coupled by adding a cross-term to the inputs:
The parameters β* _{r}* and β

*are positive, not necessarily equal, and represent the excitatory coupling strength between integrators. In the following, we refer to the case in which β*

_{s}*and β*

_{r}*are equal as the symmetric mutual excitation model, whereas we call the case in which β*

_{s}*= 0 the asymmetric excitation model.*

_{s}The models with gain modulation, common noise, shared signal, and mutual excitation (symmetric and asymmetric) have the same number of free parameters. Given that the complexity of the models is the same, we used the goodness of fit of these different architectures to select the model that accounts best for the data.

### Simulations

The dynamical Equations 1 were integrated with a time step *dt* = 0.5 ms using an extension of a Heun algorithm to include stochastic terms (Honeycutt, 1992) and drawing pseudorandom numbers from a Mersenne–Twister generator. First-passage times to threshold were estimated with a first-order interpolation between time steps. RTs were modeled as a sum of the first-passage time to threshold and a fixed residual time *T*_{0}, which was a free parameter of the model. Although we are using eight different parameters in the single integrator model [namely, τ, α, *g*, θ, *H*, ς, and *T*_{0}, and the scale of input switch *e _{i}*(

*t*), Eq. 6], the first-passage times of the integrator depend only on the combination of parameters (Smith, 1995): which represents the distance between the effective drive and the threshold scaled by the intensity of the fluctuations. Parameter values are therefore undetermined with respect to first-passage times and, consequently, to RTs. We set

*H*= 1, ς = 0.1, θ = 0.5,

*e*= {0,1} throughout and took the remaining parameters as free. Also, the gain

_{i}*g*was 1 except in the model with gain modulation.

To estimate the mean and the Pearson's correlation coefficient of the simulated RTs for each particular parameter set, we ran a sample of 10,000 trials using different realizations of the noise. We used such high number of simulated trials to guarantee a good estimate of the Pearson's correlation coefficient, with an SD smaller than 0.05. The confidence intervals of the Pearson's correlation coefficient were estimated from a bootstrap of 1000 replicates.

## Results

Saccade and reach RTs in a simple eye–hand coordinated movement are correlated. Figure 1*A* presents SRTs plotted against RRTs during an example behavioral session when monkey J was required to make both a saccade and a reach from an initial central touch and fixation point to an eccentric target. Each data point represents SRT and RRT for a single trial. The RTs for these coordinated movements are variable (SRT = 219 ± 37 ms, RRT = 294 ± 41 ms, mean ± SD) and correlated (*R* = 0.55, *p* < 0.05).

To understand the significance of RT correlations, we studied a dual reaction time task in which subjects had to respond with a saccade and a reach to their corresponding go cues. Our first aim was to study how saccade and reach RTs are affected by the degree of coordination between movements. A simple way to manipulate the degree of coordination is by cuing the two movements with a temporal separation between cues, or SOA, that is varied on a trial-by-trial basis (Fig. 1*B*) (see Materials and Methods). The idea is that movements that are cued separately and unpredictably in time will be less coordinated than those that are cued simultaneously. We used the dependence of the statistical properties of RTs on SOA to constrain the set of possible mechanisms underlying eye–hand coordination. We focused in particular on the dependence of the means and correlations between saccades and reach RTs. We first hypothesize about mechanisms that could give rise to correlations in the RTs of saccades and reaches and explore these mechanisms using a family of diffusion models of RT. The basic assumption of these models is that RTs are dependent on the time taken for an activity-dependent random walk to reach some prescribed threshold. Although these models are purely phenomenological, recent experimental work suggests that such an accumulation process may be encoded in the activity of neural populations in cortical areas (Hanes and Schall, 1996; Shadlen and Newsome, 2001; Roitman and Shadlen, 2002). Inspired by this link to the neurobiology, we explore several plausible mechanisms that could give rise to correlated RTs in a dual reaction time task. Using a framework of two integrate-to-threshold units, we have studied how RTs are affected when simultaneous gain modulations, common noise, shared signal inputs, and mutual coupling affect the integration process (see below, Model predictions). Each model makes testable predictions for how the mean RT of saccades and reaches and the correlations between the two should vary as movements are increasingly separated in time. We then analyze behavioral data collected from three monkeys performing the SOA task and compare the results with the predictions of the models (see below, Behavioral results). The models with simultaneous modulation and shared input fail to account for the behavioral data. Only the model with mutually excitatory units captures the main trends in mean RTs and correlations. Based on the behavioral results, we then develop interacting integrate-to-threshold models and compare the model predictions with behavior (see below, Interacting accumulator models). Finally, we mathematically fit the parameters of an interacting integrator model to the behavioral data (see below, Model fits to behavioral data).

### Model predictions

We studied models of two neural integrators that determine the RT of saccades and reaches. Each neural integrator is associated with one particular movement and responds selectively to the associated go cue signal by undergoing a stochastic rise of activity (Fig. 2). The behavioral response is assumed to initiate when the buildup of activity reaches a threshold. Immediately after the integrator hits threshold, the selective input responsible for the buildup is switched off, causing the integrator to decay to zero. In this framework, there exist several mechanisms that may potentially give rise to correlations between saccade and reach RTs. We first examine the predictions of each model for mean saccade and reach RT and their correlations as SOA varies.

#### Independent integrators

If saccades and reaches are generated by two parallel sensorimotor pathways and two independent neuromodulation mechanisms, an integrator model with two independent integrators would describe their behavior (Fig. 2*A*). This hypothesis gives the simplest dual-integrator model. The lack of interaction between integrators results in SRTs and RRTs that are uncorrelated and do not depend on SOA (results not shown). The presence of RT correlations when movements are instructed together rules out this possibility.

#### Gain modulation

The neuronal populations involved in the preparation of saccades and reaches may experience a common neuromodulatory signal that affects their sensitivity to inputs (McCormick, 1989), akin to the effect of motivation or arousal. It has been shown that arousal effects are mediated by monoaminergic and cholinergic ascending pathways from the brainstem and the basal forebrain (Sarter et al., 2005; Herrero et al., 2008), which generally change over hundreds of milliseconds. Given these relatively long durations, neuromodulation is likely to alter response times on a trial-by-trial basis, but it is probably too slow to induce appreciable changes within the course of a single trial. Despite their relative slowness, these simultaneous and global changes in the circuitry encoding movement preparation could result in correlations between RTs.

A simple model of common neuromodulation can be implemented with two integrators whose sensitivity to inputs is varied identically in every trial (Fig. 3*A*, left). Sensitivity to inputs is controlled by the gain of the integrator, which ultimately reflects the excitability of the neural population encoding movement preparation. The gain, although variable from trial to trial, is identical for both integrators and constant throughout a trial to account for the global and slow character of neuromodulatory signals. With such a model, neither mean RTs nor correlations will depend on SOA (Fig. 3*B*,*C*, left) because there is no interaction between integrators and no temporal dependence of the passive parameters. Common neuromodulation can, however, give rise to RT correlations. These arise because both responses will be fast when the integrator gains are high and slow when gains are low. Consistent with this intuition, the magnitude of RT correlations depends on the variance of intertrial gain modulations (Fig. 3*C*, left).

If we relax the assumption of static gains by letting the gain change with some characteristic time constant throughout the trial, the RT properties depend only on SOA at the time constant of the gain changes (data not shown). This result reinforces the intuition that models with arousal-like gain modulation predict that changes in RT correlations do not occur quickly (i.e., within a reaction time).

The gain modulation model shows that effector-specific mechanisms of coordination are not necessary to generate RT correlations. Nonspecific arousal-like mechanisms can generate RT correlations, but the correlations do not depend on SOA.

##### Fast interactions.

RT correlations may arise from common inputs driving the areas preparing and initiating reaches and saccades. We will distinguish between common inputs in the form of shared signals and of common noise (Gawne and Richmond, 1993; Averbeck et al., 2006). By signal, we mean the component of the input that encodes information specific to the preparation of movements. In the model, the signal is the deterministic component of the input, which drives the unit to threshold. Noise is lumped as the stochastic component of the input. Although one might expect correlated inputs to be a mixture of shared signals and common noise, we consider each factor separately to isolate their role in the dependence of RT correlations on SOA.

#### Common noise model

The areas responsible for the initiation of saccades and reaches may have shared fluctuating inputs. This common noise can arise from the activity of brain regions that innervate the areas guiding both movements but do not carry information about the signals specific for the task, such as sensory information about the go cues. This form of nonspecific activity can be modeled as correlated noise in the inputs to both integrators generating the movements.

To examine the effect of common noise fluctuations in isolation from signal changes, we varied the correlation coefficient, *c*, between the external input driving both units (Fig. 3*A*, middle) (see Materials and Methods). The coefficient *c* can be interpreted as the fraction of fluctuating inputs that are shared among units when the mean and the variance of the inputs are kept constant. Because the deterministic component of the input is unaffected by the amount of correlated noise, mean RTs are not modulated by *c* (Fig. 3*B*, middle). In addition, mean RTs do not depend on the SOA because the signals of both integrators are effectively independent of one another. In other words, fluctuations in the input to one unit do not convey any information about the signal of the other unit.

Common noise introduces RT correlations (Fig. 3*C*, middle). These correlations arise because both units are driven to threshold by the same fluctuations when they are activated close in time, as occurs at short SOA. Also, the correlation coefficient between SRTs and RRTs is larger when the incoming noise among units is more correlated. Correlations in RT decay monotonically as SOA increases and approach zero when SOA is longer than the mean SRT. This trend results from the shortening of the time window during which both units share the same fluctuating inputs.

The common noise model illustrates that effector-specific mechanisms of coordination are not necessary to generate RT correlations that decrease with SOA. Nonspecific mechanisms can generate RT correlations that decay with SOA but do so in a manner that influences mean SRTs and RRTs equally and independently of SOA.

#### Shared signal model

The preparation stages of saccades and reaches may interact through some type of feedforward excitatory shared signal. For example, a common sensory stage in the processing pathway would feed into areas responsible for preparing the saccade and reach. This would influence subsequent RTs.

To model shared signals, we introduce a deterministic signal received by a unit that includes not only the signal specific to it but also some fraction of the signal driving the other unit (Fig. 3*A*, right) (see Materials and Methods). Shared signal interaction provides an implementation for the so-called “energy integration” models, according to which the “energy” of different sensory modalities is assumed to sum somewhere in the nervous system (Todd, 1912; Hershenson, 1962; Diederich, 1995).

In the shared signal model, mean RTs for reaches and saccades are the same, as in the common noise model. Forward excitation speeds up the integration to threshold if there is some overlap in the stimulation periods of both units, resulting in shorter RTs. Therefore, mean RTs depend on SOA, unlike the common noise model. RTs are on average faster, and therefore movements are facilitated, at short SOAs when the overlap between saccade and reach preparation is longer (Fig. 3*B*, right). Facilitation disappears for SOA longer than the mean RT of the independent integrator.

Interestingly, unlike the common noise model, shared signal excitation gives rise to negative RT correlations when reach and saccade preparation overlap in time (Fig. 3*C*, right). The values of the RT correlations are relatively constant for short SOAs, unlike RT correlations caused by common noise, which progressively decay as SOA increases. The reason for anticorrelation between the RTs is that faster saccadic responses lead to shorter overlaps between preparation stages. This leads to shorter windows of facilitation that result in slower reach responses and hence anticorrelated RTs.

Comparing model predictions across the three integrator models demonstrates that each model makes specific predictions about the mean RTs and the RT correlations and their dependence on SOA. Next, we analyze behavioral RTs during an SOA reach and saccade task and compare the results with the model predictions.

### Behavioral results

The SOA task aimed to manipulate coordination between the saccade and reach by separating movements by a random time interval, the SOA duration (Fig. 1). The reasoning behind the SOA task is that saccade and reach movements made separately at large, unpredictable SOA intervals should be less coordinated than movements made together at short or zero SOA intervals. Therefore, if RT correlations reflect coordination, they should, at a minimum, fall as SOA increases. However, as the integrator models above demonstrate, RT correlations could arise from other mechanisms that we would not necessarily consider to be specific to coordination, such as common gain modulation, common noise, or shared signals. We analyze RT correlations in the behavior of three monkeys performing the SOA task to test the predictions of each model and determine whether the pattern of RT correlations with SOA and overlap reflect any of the potential nonspecific mechanisms modeled above or if another mechanism is necessary to explain the pattern of results.

#### SOA

Mean RTs for saccades and reaches vary with SOA (Fig. 4*A*,*B*). We consider SRT first. For all monkeys, SRT is generally fastest when SOA is short (Fig. 4*A*). Monkey H and monkey J show clear evidence of facilitation at the shortest SOA. They also show stable SRT for SOA greater than the mean SRT. Monkey S shows facilitation with a faster mean SRT for SOAs shorter than 250 ms except for a single interval at SOA = 75 ms. At longer SOA, the mean SRT for all three monkeys is generally constant. Faster SRT for shorter SOA was not attributable to the use of an auditory cue for the reach. SRTs for a modified visual-visual SOA task (see Materials and Methods) were the same as SRTs in the main visual–auditory SOA task (monkey H: visual–visual SRT = 195 ms; visual–auditory SRT = 194 ms; *p* = 0.53; rank sum test). The facilitation of SRT at short SOA is consistent with the shared signal model but not with either the gain modulation or the common noise models.

Next, we consider mean RRT. For all three monkeys and all SOA intervals, RRTs are longer than SRTs (rank sum test, *p* ≪ 0.05). Mean RRTs also vary with SOA significantly differently than do mean SRTs. On average, saccades are systematically faster at short SOA. In contrast, RRT varies non-monotonically with SOA: reaches are faster at SOA shorter than 100 ms, slower at intermediate SOAs between 125 and 225 ms, and faster again at long SOA (Fig. 4*B*). In monkey H and monkey J, there was an exception to this trend at SOA = 25 ms, which elicited the slowest RRT. Here, the effect of overlap that we expect to see on RRTs could be masked by other factors. RTs have been shown to be faster when more time is given to prepare before a go cue. Because the location of the target is given when the saccade is cued, longer SOAs give the reach a longer preparation time. The short SOA leaves the monkey little additional time to prepare, leading to longer RTs. Monkey S is then an anomaly, but note that the mean RRT for monkey S is much longer. Monkey S may take longer to prepare to reach, making the short preparation time much less important.

None of the nonspecific/input models predict the pattern of RT variations with SOA that we observe.

Next, we consider RT correlations. RT correlations are significantly greater than zero at some short SOA (<100 ms) in all three monkeys (Fig. 4*C*). Importantly, RT correlations decrease with increasing SOA, becoming statistically insignificant after several hundred milliseconds.

Precisely how RT correlations vary with SOA differs between monkeys. For monkey H, correlations start high (*R* = 0.40, *p* ≪ 0.001), drop with SOA, and stay significantly positive for SOA shorter than 350 ms (Fig. 4*C*, left). Correlations start high (*R* = 0.49, *p* ≪ 0.001) for monkey J, fall with SOA, but become significantly negative (*R* = −0.13, *p* < 0.001) for SOA around 300 ms, decaying to zero for longer SOA (Fig. 4*C*, middle). For monkey S, RT correlations are highest when the SOA is shortest (*R* = 0.27, *p* ≪ 0.001) but quickly drop to zero before staying slightly but significantly positive for longer SOA durations around 300 ms and returning to zero (Fig. 4*C*, right).

The variations in behavioral RT with SOA significantly differ from the predictions of all three integrator models presented above. The pattern of RT correlations is not consistent with the gain modulation model or shared signal model because RT correlations vary with SOA and are positive at short SOA. The general decrease in RT correlations with SOA is consistent with the common noise model, but the more complex dependence of RT correlations with SOA, such as the negative RT correlation around SOA = 300 ms seen in monkey J, is at odds with the common noise model. The variation of mean SRT and RRT with SOA is also inconsistent with the common noise model.

Because nonspecific/input mechanisms do not predict the RTs observed in the SOA task, effector-specific interactions between the integrators may provide more accurate models. To gain insight into appropriate model architectures, we further examined the behavioral data. We observed that the time elapsed between the auditory reach cue and the time of the saccade, or overlap, may play an important role in understanding the behavioral RTs. For example, the behavioral data showed that the mean RRT increased at intermediate SOAs when the auditory reach cue was delivered around the time of the saccade. In addition, the distribution of SRTs for trials with overlap of 125 ms, when the reach cue was delivered just before the saccade was initiated, is bimodal (not shown here). Hence, we examined SRT and RRT and their correlations as a function of overlap. In the model, the overlap corresponds to the period during which both integrators are activated (Fig. 2*B*).

#### Overlap

SRT is bimodal for 125 ms overlap in all three monkeys. Therefore, we divided trials into fast and slow subgroups by choosing a cutoff SRT that was at the minimum between the two groups of RTs. This cutoff was 180 ms for monkeys J and H and 190 ms for monkey S. The fast SRTs represented 40% of trials for monkey J, 61% of trials for monkey H, and 68% of trials for monkey S. Figure 5 presents the SRT, RRT, and their correlations as a function of overlap. For each measure, data for 125 ms overlap is presented for the fast (solid) and slow (dashed) SRT groups. In all three monkeys, SRT was stable for overlap durations shorter than 125 ms and was fastest for 175 ms overlap (Fig. 5*A*). A bimodal distribution of SRTs is consistent with interference between the reach cue and saccade initiation on a fraction of trials. Note that, when we separate fast and slow SRTs for an overlap of 125 ms, all three subjects are faster than average for the fast SRTs. For the slower SRTs, monkeys H and J have mean SRTs that are no different than the mean SRTs for shorter overlaps. Monkey H, however, has a slower-than-average SRT. Adding the reach soon after the saccade go is given may disrupt saccade planning in monkey S, leading to the anomalous mean SRT at 75 ms. The faster SRTs in all three monkeys may be attributable to an influence of the reach go that speeds initiation of the saccade on some trials.

Consistent with the dependence of the mean RRT on SOA (Fig. 4*B*), the mean RRT reaches a maximum around overlap 0 ms (Fig. 5*B*). That is, the RRT is on average slower when the reach cue is delivered at the time of saccade initiation. A first plausible explanation for this slowing is that there is bottleneck in the processing of sensory information when the saccade is initiated. The bottleneck would only affect the initiation of reaches, because we do not observe any consistent increase of mean SRT at that time. The explanation is, however, incomplete because bottleneck effects predict consistent increases in the mean RRT for positive overlaps, not only for 0 overlaps, and cannot explain the non-monotonic dependence of the mean RRT on SOA. An alternative explanation is that this non-monotonic dependence results from a combination of two independent factors: facilitation at short SOAs (large, positive overlaps) and reduction of uncertainty at long SOAs, which results from choosing a uniform distribution for SOA durations. We discuss these two effects in detail in the next section.

The correlation of SRTs and RRTs varies with overlap consistently across the three monkeys. RT correlations are different from zero for positive overlaps and vanish for negative overlaps (Fig. 5*C*). There are, however, some differences across subjects. In monkey H and monkey S, RT correlations remain positive for all positive overlap values, but in monkey J, correlations become negative for overlaps shorter than 125 ms. The crossover of RT correlations for positive overlaps is notable because it is not predicted by any of the models described above.

#### Interacting integrator models

Interactions between integrators associated with each movement may be required to explain these behavioral results because the overlap in time during which both integrators are active affects RTs and their correlations. In particular, we find that RT correlations are only significantly different from zero when overlap is positive. We consider two models with direct coupling between integrators (Fig. 6) (see Materials and Methods). A similar mechanism, based on inhibitory coupling and a stochastic race model, has been proposed as a model for movement inhibition in a countermanding task (Boucher et al., 2007a). In the symmetric mutual excitation model, reach and saccade integrators directly excite each other. In the asymmetric excitation model, the reach integrator directly excites the saccade integrator but the saccade integrator has no influence on the reach integrator. In the following, we examine the predictions of the interacting integrator models for SRT, RRT, and their correlations.

##### Symmetric mutual excitation.

The net input to each unit includes a positive coupling term proportional to the activity of the partner unit. The coupling term induces excitatory interactions between the units, leading to a positive feedback loop that causes the activity to grow without bound unless there is saturation or a threshold detection mechanism (for an example of the latter, see Lo and Wang, 2006). In contrast to the common noise and shared signal models, coupling induced by mutual excitation depends on the instantaneous activity level of the units. This results in more complex SOA dependencies (Fig. 6, left).

Mutual excitation increases the inputs to both integrators. The input facilitates both saccades and reaches and shortens the mean SRT and RRT (Fig. 6*B*,*C*, left). As in the shared signal model, the degree of facilitation depends on how much overlap there is between the stimulation periods of both units. Importantly, however, with symmetric mutual excitation the mean RRT and SRT depend on SOA differently. Mean SRTs decrease as SOA becomes shorter. Thus, saccades, which are instructed first, occur faster when the reach is activated soon after the saccade. In contrast, mean RRT does not decrease monotonically with decreasing SOA. Instead, mean RRT reaches its minimum value for SOA at around the mean SRT, when the saccade unit reaches its peak of activity and drives the reach unit most strongly. Mean RRT is therefore shortest when the SOA is approximately equal to the mean SRT.

Symmetric mutual excitation leads to RT correlations that are positive for short SOA, negative at intermediate SOA, and approach zero at long SOA (Fig. 6*D*, left). At short SOA, RT correlation is positive and increases with coupling strength. RT correlation is positive because the integrators associated with each movement effectively share the same fluctuating inputs through mutual excitation. The effect is similar to that seen in the common noise model, although the nature of the fluctuations differs between the two models. In the common noise model, fluctuations shared with the other integrator enter simply as an additive term. In the symmetric coupling model, shared fluctuations are low-pass filtered and enter indirectly through the firing activity of the partner unit. This difference explains why correlations stay significantly different from zero for longer SOA periods in the symmetric mutual excitation model compared with the common noise and shared signal models, which both reflect correlated inputs. Another difference between the symmetric mutual excitation model and both correlated input models shows up at longer SOA. As SOA increases, RT correlations decrease and eventually become negative for SOA longer than the mean SRT. The reason for this change in sign is analogous to that given for the shared signal model. Fluctuations leading to shorter-than-average SRT lead to shorter-than-average windows of activity for the saccade integrator and, consequently, to less facilitation for the reach unit. Faster-than-average saccadic responses are therefore paired with slower-than-average reach responses, leading to anticorrelated RTs. RT correlations vanish for sufficiently long SOA because the reach and saccade units become effectively decoupled when they are activated far apart in time.

When plotted as a function overlap, RT correlations are zero for negative overlap and show a negative-then-positive pattern for positive overlap (Fig. 6*E*, left). Thus, mutual excitation critically depends on the relative timing between the reach cue and SRT. If the reach cue is delivered after the saccade, there is no correlation between the SRT and the RRT.

##### Asymmetric excitation.

To better understand the role of mutual excitation in the pattern of RT predictions, we examined an asymmetric excitation model in which reach activity affected saccade activity but there was no reciprocal connection (Fig. 6*A*, right). The parameter β is the coupling strength from the integrator for reaches to the integrator for saccades. Because the reach unit projects to the saccade unit but does not receive any connection from it, saccade activity has no influence on reach activity. Therefore, only SRTs are shorter as a result of simultaneous coactivation when SOA is sufficiently short (Fig. 6*B*, right); mean RRTs are not modulated by β or the SOA (Fig. 6*C*, right). The same effect can be restated in terms of the overlap. The effect of facilitation is stronger with more overlap between preparatory stages (data not shown). As in all other models, facilitation disappears when the SOA is sufficiently long. Moreover, RT correlations increase as SOA gets shorter (Fig. 6*D*, right) and overlap increases (Fig. 6*E*, right).

In summary, when the reach unit projects to the saccade unit but does not receive any connection from it, only saccades will benefit from the simultaneous coactivation occurring at short onset asynchronies, and therefore only SRTs will show facilitation. Analogously, saccade facilitation will be stronger when the overlap in preparation is longer. Facilitation decreases monotonically with the SOA and vanishes for SOA longer than the mean SRT. When SOA is longer than the mean SRT, the extra input provided by the reach unit comes too late to help the saccade unit reach threshold earlier. Moreover, with asymmetric excitation between units, correlations are positive at very short SOA, decrease monotonically as SOA get longer, and approach zero when SOA is on the order of the mean SRT.

### Model fits to behavioral data

Of all the models considered, only the model with asymmetric mutual excitation captures the dependencies of the behavioral measures on SOA and overlap. More specifically, it is the only model predicting the facilitation of saccades at short SOA and the decrease of positive RT correlations with SOA. Given that all the models have the same number of parameters and hence similar complexity, we favor the model with asymmetric mutual excitation as the one best describing the data. We fit the empirical data to the parameters of this model, considering a more general case of the mutual excitation model in which the coupling strengths between integrators take arbitrary values (Fig. 7). This allowed us to study the continuum superset of models comprised between the symmetric and the asymmetric excitation models. Note that, in doing so, we are adding another parameter in the model and we are therefore increasing its complexity. It is then legitimate to ask whether other models of similar complexity could fit the data equally well. A combination of the shared noise and the shared signal models may also give rise to the observed dependences on SOA with both facilitation at short SOAs and decrease of correlations for increasing SOA. However, we have found that it is not possible to have both facilitation and a monotonic decrease of RT correlations with SOA in a model combining shared noise and shared signal. Instead, we have seen that the model with mutual excitation provides a better accounting for the data (data not shown). A more quantitative comparison between models, based on the calculation of the so-called Bayesian factors (the ratio of the probabilities of the data given each model), lies beyond the scope of this work.

The experimental data used for the fits comprised the data in Figure 4, *A* and *C*, which characterize the dependence of mean SRT and correlations on the SOA. We estimated model parameters using an optimization algorithm that minimized the difference between the experimental data and the predictions from the model, weighing each difference according to the precision of each experimental data point (for details, see Materials and Methods). The values of the best-fitting parameters are summarized in Table 1.

Consistent with the predictions of the mutual excitation models, the best fit coupling parameters between the integrators are unequal. In all three monkeys, reach-to-saccade coupling β* _{r}* is at least three times as high as saccade-to-reach coupling β

*. The asymmetry in the coupling strengths obtained from optimizing the data was expected given the lack of significant negative correlations for SOA longer than the mean SRT. Recall that, for models with symmetric coupling, correlations are positive at short SOA and negative at longer SOA, with comparable magnitudes.*

_{s}In general, the mutual excitation model accounts for the facilitation of saccades observed at short SOA. According to the model, mean SRT should start decreasing from its baseline value as the SOA decreases below the mean SRT. This is because, on average, the reach can speed up the accumulation process for the saccade only when overlap is positive. The magnitude of the facilitation should be accentuated for shorter SOA (longer overlap) for the same reason. These two predictions are consistent with the observations for all three subjects, although there are some differences in the magnitude of facilitation as well as in the SOA at which facilitation becomes noticeable. For instance, the mean SRT for monkey S does not show a clear monotonic dependence on the SOA, in contrast to the other two monkeys. This explains why the best-fit value for the recruitment time constant τ of the integrator is considerably shorter than for monkeys J and H and, therefore, why the model predicts a mean SRT and an RT correlation that shows no clear dependence on SOA except at very short SOA (<150 ms).

Our modeling efforts have been focused on how well the models account for the dependence on SOA of the mean SRT and the RT correlations. We have not considered the dependence of the mean RRT on SOA because RRT are likely to be affected by expectancy effects, which are not accounted for by our simple models. In the dual RT task studied here, in which SOAs are random and uniformly distributed, the saccade go signal provides information about when the reach go signal is going to occur. This is because the probability of appearance of the signal, given that it has not yet occurred, increases as time progresses. Subjects are therefore more certain about when the second signal is going to appear as time elapses since the appearance of the first signal. It has been reported that such increases in expectancy consistently reduce RTs (Klemmer, 1957; Nickerson, 1965), and that, in general, the mean RT is a monotonic decreasing function of the conditional probability of appearance of the signal, or hazard function (Stilitz, 1972). It is thus reasonable to expect that the mean RRT at long SOAs is reduced with respect to the mean RRT that would result in the absence of expectancy effects.

When combined with the predictions of our model, expectancy effects can help explain the observed non-monotonic dependence of RRT on SOA. According to our model, RRT are mildly facilitated at short SOAs. Conversely, RRTs are shortened at long SOA because of the increase of certainty brought about by the passage of time. Combined together, these two factors leave a window of SOA values within which reaches are not speeded up, leading to a non-monotonic dependence of mean RRT on SOA.

The reduction of RRT induced by increases in expectancy could be simply implemented in our model by an extra input to the reach integrator. This input would be a monotonically increasing function of the hazard function of the SOA, which for a uniform distribution is approximately constant for short and mid-range SOA and grows increasingly faster at long SOA, until it blows up when the SOA reaches its upper value. Because the exact dependence of the input on the hazard function is unknown, introducing such input in the model would require adding more free parameters to the description, which would increase the complexity of the model.

## Discussion

To better understand the mechanisms behind eye–hand coordination, we designed an SOA task in which reach and saccade go cues were separated in time by a variable time interval. We then studied RT correlations and mean RTs as SOA varied. RT correlations were high when movements were cued at the same time. When movements were separated by a few hundred milliseconds, those correlations disappeared. SRTs were on average faster than RRTs when reach and saccade go cues were simultaneous, suggesting facilitation when saccades are made with a reach.

We compared the behavioral data with predictions from a family of network models implementing integrate-to-threshold dynamics. In all models, each movement RT is defined by the time needed for neuronal activity to reach a prescribed level, consistent with the finding that RTs are strongly correlated with the firing activity of areas coding for motor planning and preparation (Riehle and Requin, 1993; Hanes and Schall, 1996; Roitman and Shadlen, 2002), as well as with the phenomenological diffusion models of RT (Ratcliff, 1978). We have considered a family of models with two neuronal populations that encode saccade and reach preparation. The working hypothesis is that facilitation and correlation of RTs reflects interaction between neural populations involved in preparing movements. We have considered several plausible mechanisms of interaction, including common neuromodulation, shared inputs, and direct coupling, and have proposed simplified models for each of them, to disentangle the effect of each type of interaction on the dependence of the statistics of RTs on the SOA.

Of all the models we analyzed, only the model with two mutually coupled integrators can account for the behavioral data. The model is sufficiently specific to capture the decrease of RT correlations with SOA and the increase of mean SRT with SOA but also sufficiently flexible to account for the differences across subjects. The undershoot in correlations, which is observed in monkey J for SOA longer than the mean SRT, provides an interesting test of the model architecture. According to the integrator model, the undershoot arises because, for SOAs slightly longer than the mean SRT, saccades that are slower than average lead to delayed peaks of saccade-related activity that speed up the reach. As a result, SRTs and RRTs are anticorrelated. Importantly, this analysis depends on trials for which SOA is greater than the mean SRT across trials. When RT correlations are recalculated in terms of overlap, defined as the time elapsed between the reach cue and the saccadic response, RT correlations are negligible at negative overlaps, i.e., when the reach cue is delivered after the SRT. RT correlations vanish at negative overlaps because the time course of activity of the reach unit, which ultimately determines RRT, is independent of saccade activity that determined the SRT. Zero RT correlations for negative overlaps are a strong prediction of the integrator model. This prediction is borne out by the behavioral data in all three monkeys.

The models proposed here are idealized descriptions of the dynamics of neural populations and their potential interactions. We chose a minimal description that incorporates the basic elements of neural processing without attempting to link model parameters to the underlying neurophysiology. More biologically realistic models, such as a network of spiking neurons with explicit synaptic and neuronal dynamics, will not necessarily change the behavioral predictions if the main assumption, that behavioral responses are triggered when neuronal population activity reaches a threshold, remains the same. Predictions also rely on neuronal activity smoothly decaying after reaching threshold. Smooth decay is consistent with the observation that the neuronal activity correlated with RTs decays gradually after reaching its peak (Hanes and Schall, 1996) and is implemented in our model by turning off the external inputs, which makes the after-threshold activity decay. Although this particular choice shapes the dependence of mean RRTs and RT correlations on SOA, we do not expect these dependences to change dramatically with other reset mechanisms.

### Reaction times and eye–hand coordination

SRTs and RRTs are often correlated when the movements are made together, although the degree and underlying cause of correlations in various eye–hand tasks is unclear. Previous studies have reported RT correlations ranging from very high (*R* = 0.8) (Herman and Maulucci, 1981) to moderate (Prablanc et al., 1979a; Gielen et al., 1984; Fischer and Rogal, 1986) to relatively low (*R* < 0.5) values (Biguer et al., 1982). Bekkering et al. (1994) found that eye movements were faster when made separately from a hand movement to the same target, suggesting interference instead of facilitation. However, the study used very large stimuli to emphasize speed instead of accuracy and reported that, unlike previous eye–hand movement studies (Prablanc et al., 1979a,b; Herman and Maulucci, 1981; Gielen et al., 1984; Fischer and Rogal, 1986; Abrams et al., 1990), the eye did not tend to move before the hand and RT correlations were small. A saccade and pointing study also found saccades slowed when combined with pointing, although manual responses were faster when made with saccades (Mather and Fisk, 1985). Hence, many factors can influence relative eye–hand RTs.

Correlations can also result from common sensory processing. In a sensory, perceptual paradigm, if the sensory cues used to initiate movement are difficult to detect or discriminate, variations in the duration of sensory processing will result in significant trial-to-trial variations in RT. When sensory processing drives two movements, the trial-to-trial variations in sensory processing will influence RT for each movement, giving rise to stronger RT correlations. The shared signal model we present might describe shared inputs attributable to sensory processing.

Shared input signals may not only be important during sensory, perceptual paradigms but may also describe subjects that are uncertain about task contingencies, for example, during early stages of behavioral training or exposure to a task. Task uncertainty could lead to a processing stage common to both effectors whose variability drives RT correlations. During initial training for this task, we observed large RT correlations (*R* > 0.8) that may be attributed to lack of practice. Once subjects became familiar with the task, RTs became faster and RT correlations fell, stabilizing at the values presented here (0.4–0.6).

Our results speak to the question of whether reaches are more strongly influenced by coordinated saccades or vice versa, that is, whether we look to where we reach or whether we reach to where we look (Carey, 2000; Horstmann and Hoffmann, 2005). We demonstrate the best fit when the reach-specific integrator drives the saccade-specific integrator. Asymmetric coupling suggests that eye movements can be tied to hand movements. Fisk and Goodale (1985) studied RTs of eye and arm movements to visual targets presented ipsilaterally or contralaterally to the starting stimulus. Interestingly, SRTs were slower when made with a contralateral reach, although there was no difference in RTs when saccades were made without a reach.

Coordination is flexible, and eye–hand interactions may be highly task specific. Here, the instruction to reach after the saccade may have led to asymmetry in the underlying behavioral mechanism. If so, instructing subjects to reach and then saccade could yield the reverse pattern of results, effectively switching SRT and RRT in the results above and implying asymmetric coupling of the saccade to the reach. This would demonstrate dynamic functional coupling, but additional experiments are needed to directly test this alternative.

### Dual-task paradigms

Dual-task paradigms have a long history in the study of psychological processes (Pashler, 1998). Sometimes, these tasks interfere with each other as a result of bottlenecks in sensorimotor processing or central processing. According to the single-channel hypothesis (Welford, 1952; Broadbent, 1958), there is a delay in the processing of a second signal while the first signal is being processed. The delay is taken to be a sign of limited capacity central processing. Evidence for capacity limitations comes from results that show the RT of the second of two closely spaced signals is slower and inversely proportional to the time between signals (Smith, 1967). Slowing of the second response is called the psychological refractory period (Welford, 1952). The increased RRT at short SOA that we observe shares similarities to the traditional psychological refractory period. Consequently, slowing around the SRT could be attributable to bottleneck effects in a channel that processes incoming sensory information. However, unlike classic psychological refractory effects, the longest RRT occurs around the saccade and is not inversely proportional to the SOA.

Facilitation of the SRT for short SOA is a prominent effect in our behavioral data. However, facilitation is not necessarily explained by a capacity limit. Dual-task paradigms also depend on task difficulty, because concurrent performance of easy tasks is speeded compared with individual performance, whereas concurrent difficult tasks are slowed (Keele, 1968). Coordinated eye–hand movements comprise a relatively easy task and could lead to speeded RTs through mutual excitation. A more complete model of the RT behavior in the SOA task could integrate the single-channel models for sensory input processing with other mechanisms that capture facilitation and correlation in RTs.

## Appendix: parameter fitting

The data points in Figure 4, *A* and *C*, were fitted with the mutual excitation model described by Equations 5 and 9, using the *p* = 5 free independent parameters of the model: τ, α, β* _{r}*, β

*, and*

_{s}*T*

_{0}, which we denote collectively as θ. To simplify the optimization method, we binned the raw dataset into

*N*= 6 bins that summarized the experimental data. The parameter values shown in Table 1 were estimated by minimizing the weighted sum of squared errors where

*S̅R̅T̅*

_{exp}

^{i}and

*S̅R̅T̅*

_{model}

^{i}are, respectively, the experimental and simulated mean SRTs within bin

*i*, whereas

*R*

_{exp}

^{i}and

*R*

_{model}

^{i}refer, respectively, to the experimental and simulated correlation coefficients of RTs in bin

*i*. The symbols ς̂

_{S̅R̅T̅expi}

^{2}and ς̂

_{Rexpi}

^{2}denote the half-lengths of the 95% confidence intervals of the experimental estimates of

*S̅R̅T̅*

_{exp}

^{i}and

*R*

_{exp}

^{i}, respectively. The function S(θ) quantifies the differences between the predictions of the model and experimental data by summing the squared error between measures, weighing more heavily data points that are measured with higher precision. The estimates for the mean RTs are more precise than the estimates for the correlation coefficient. From Equation 10, this implies that errors for the mean SRT contribute much more to the cost function than the errors for the correlation coefficient. To compensate for the overrepresentation of mean SRT errors to the cost function, we have introduced the coefficient

*b*> 1, whose value is estimated by requiring the two terms in Equation 10 to have the same order of magnitude. Because the range of observed SRTs and correlations is ∼50 ms and 0.5, respectively, and the average error for the estimated mean SRT and correlations is ∼2 ms and 0.1, the ratio between relative errors of the mean SRT and of RT correlations is This ratio could be used as a value for the compensation factor

*b*. Introducing a compensation factor comes at the price of biasing the fits toward less precise experimental measures. We need therefore a compromise between giving correlation errors a decent share in the cost function (

*b*≃ 5) and weighing data points proportionally to their precision (

*b*= 0). In our optimization algorithm, we set

*b*= 1.5.

There are no simple closed-form expressions for *S̅R̅T̅*_{exp}^{i}(θ) or *R*_{exp}^{i}(θ) and hence for the cost function *S*(θ). For this reason, we minimized *S*(θ) over the parameters using a direct search algorithm in parameter space. We used the downhill simplex algorithm by Nelder and Mead (1965), which does not require the evaluation of derivatives of the cost function and which is thus appropriate for optimizing stochastic models (Bogacz and Cohen, 2004). In short, the downhill simplex algorithm evaluates the cost function at *p* + 1 different vertices of a simplex in the *p*-dimensional parameter space and replaces the vertex with highest cost by a new vertex with lower cost. This new vertex is found by shrinking, expanding, or reflecting the high-cost vertex of the original simplex. Each evaluation requires computing the values of *S̅R̅T̅*_{model}^{i}(θ) and *R*_{model}^{i}(θ), which were estimated from a sample of 10^{4} realizations of the dual integration process with parameter θ and with SOA durations drawn from a uniform distribution within the *i*th bin interval. The algorithm stops when the parameters do not change from one step to the next for more than a predefined tolerance value. On average, the algorithm required ∼150 iterations before stopping.

Direct search methods, such as the downhill simplex algorithm, can find an area with a local minimum relatively fast but cannot determine accurately the precise location of the minimum. This lack of convergence is worsened by the fact that the cost function includes sample estimates, which are inherently stochastic. More importantly, these algorithms may converge to a local, rather than a global, minimum. To minimize these problems, we ran the downhill simplex algorithm 100 times with different noise realizations. The *p* + 1 lowest minima resulting from the 100 simulated trials were then selected as the vertices of a new simplex, which was used as the initial simplex in a new algorithm execution. For this last execution, we increased the sample size to 10^{5} realizations per parameter setting to get more accurate estimates. The optimization method is computationally expensive: each evaluation of the cost function at one particular vertex of the simplex lasts a few seconds, which implies a total process time of several days per subject. These times were reduced to a few tens of hours using the New York University cluster facilities.

## Footnotes

This work was supported, in part, by National Science Foundation CAREER Award BCS-0955701, a Fellowship in Brain Circuitry from the Patterson Trust (H.L.D.), National Institutes of Health Training Grant T32 MH19524 (H.L.D.), the Leonard Bergstein Award from the Swartz Foundation (D.M.), a Career Award in the Biomedical Sciences from the Burroughs Wellcome Fund (B.P.), a Watson Program Investigator Award from New York State Foundation for Science, Technology, and Innovation (B.P.), a McKnight Scholar Award (B.P.), and a Sloan Research Fellowship (B.P.).

- Correspondence should be addressed to Dr. Bijan Pesaran, 4 Washington Place, Room 809, New York, NY 10003. bijan{at}nyu.edu