A generative learning model for saccade adaptation

Plasticity in the oculomotor system ensures that saccadic eye movements reliably meet their visual goals—to bring regions of interest into foveal, high-acuity vision. Here, we present a comprehensive description of sensorimotor learning in saccades. We induced continuous adaptation of saccade amplitudes using a double-step paradigm, in which participants saccade to a peripheral target stimulus, which then undergoes a surreptitious, intra-saccadic shift (ISS) as the eyes are in flight. In our experiments, the ISS followed a systematic variation, increasing or decreasing from one saccade to the next as a sinusoidal function of the trial number. Over a large range of frequencies, we confirm that adaptation gain shows (1) a periodic response, reflecting the frequency of the ISS with a delay of a number of trials, and (2) a simultaneous drift towards lower saccade gains. We then show that state-space-based linear time-invariant systems (LTIS) represent suitable generative models for this evolution of saccade gain over time. This state-equation algorithm computes the prediction of an internal (or hidden state-) variable by learning from recent feedback errors, and it can be compared to experimentally observed adaptation gain. The algorithm also includes a forgetting rate that quantifies per-trial leaks in the adaptation gain, as well as a systematic, non-error-based bias. Finally, we study how the parameters of the generative models depend on features of the ISS. Driven by a sinusoidal disturbance, the state-equation admits an exact analytical solution that expresses the parameters of the phenomenological description as functions of those of the generative model. Together with statistical model selection criteria, we use these correspondences to characterize and refine the structure of compatible state-equation models. We discuss the relation of these findings to established results and suggest that they may guide further design of experimental research across domains of sensorimotor adaptation.


38
The accuracy of saccadic eye movement is maintained through mechanisms of saccade 39 adaptation, which adjust the amplitude [1][2][3] or direction [4][5][6] of subsequent movements in 40 response to targeting errors. As online visual feedback cannot be used to correct the ongoing 41 movement, saccadic eye movements need to be preprogrammed and adaptation must largely rely 42 on past experience and active predictions [7,8] rather than closed-loop sensory information.

51
We recently presented a version of this paradigm in which the ISS (the disturbance responsible for 52 inducing adaptation) follows a sinusoidal variation as a function of trial number ([11,12]; see also 53 [4,13,14]). We reported that gain changes were well described by a parametric functional form 54 consisting of two additive components. One component was a periodic response reflecting the 55 frequency of the ISS that was adequately fitted with a lagged but otherwise undistorted sinusoid.

56
The second component constituted a drift of the baseline toward lower saccade gain (larger 57 hypometria) that was appropriately accounted for using an exponential dependence.

58
Here, we investigate whether a generative algorithm that models saccade gain modifications on a 59 trial-by-trial basis by learning from errors made on previous trials can account for this response. To 60 this end, we implemented and fit a series of state-space models in which a modified delta-rule 61 algorithm updates a hidden or latent variable (for which the experimentally observed adaptation 62 gain is a proxy) by weighting the last experienced visual error, in addition to other error-based and 63 non-error based learning components [8,[15][16][17][18][19][20][21][22][23] 163 Saccade amplitude adaptation is usually described in terms of the changes in saccade gain 172 The adaptation gain represents the residual of the saccade gain with respect to perfect landing.

173
When a saccade lands exactly on the first target step (a perfectly accurate saccade), the saccade 174 gain will be one while the adaptation gain will be zero. Therefore, the adaptation gain uses perfect

185
Assessment of the evidence in favor of a model. In implementing the phenomenological parameter 186 estimation, we adopted a Gaussian likelihood for the data given the model. This likelihood can be 187 maximized with respect to the parameters at a fixed but unknown width. Instead we adopted the 188 following procedure [11,13]. Using Bayes theorem, priors for the parameters to be estimated, and 189 assuming a constant prior probability for the data, we can obtain a joint probability amplitude for all 190 parameters that can be marginalized to extract individual probability amplitudes for each 191 parameter. In this process, the width of the Gaussian likelihood is a nuisance parameter that we 192 integrate out using a non-informative prior [13,28,29]. Once such integration is conducted, the 193 volume of the resulting probability density (given the data) provides an estimate of the odds that 194 the model would provide a reasonable description of the data. Here we provide a full model Page 10 of 48 consisting of six parameters (sinusoidal entraining of the oculomotor response riding over a 196 baseline drift) that we want to compare to a partial model (the drift of the baseline alone) and to a 197 minimal model consisting of the mean of the adapting block with variance equal to the variance of 198 the recorded data over that block. To establish which situation is more likely across different 199 number of parameters, we take the log of the ratio of the odds across the models. The resulting 200 magnitude is the evidence that the data are in favor of a particular model and is measured in 201 decibels (db). When this magnitude is positive, the odds favor the model in the numerator, with 202 evidence higher than 3 db indicating that this model is significantly favored to the one in the 203 denominator. We use this metric to assess the quality of our parameter estimation.

204
Statistics. Throughout the manuscript we report results as mean ± SD for individual data and 205 mean ± SEM when we discuss group data. In the phenomenological fittings, to determine average 206 parameters from the parameter estimation other than the frequency, we computed the mean and

210
Alternative estimators (e.g., the modes of the posterior distributions, with and without weighting) 211 gave qualitatively similar results.

212
Modeling of the sensorimotor learning process: the modified delta-rule state equation.

213
To investigate generative models, we adopt the following rationale. In each trial, the oculomotor 214 system must generate a motor command to produce the impending saccade. This needs to be 215 calibrated against the actual physical size of the required movement [15,20,22,24,30,31].

216
If the saccade fails to land on target, the motor command needs to be recalibrated based on 217 preexisting calibrations, and we will hypothesize that those changes take place in an obligatory 218 manner (cf. [19] ) through additive, error-based modifications attempting to ameliorate post-219 saccadic mismatches between the eyes' landing position and target location.

225
Because saccades are extremely rapid movements that do not admit reprogramming in mid-flight,

226
it is assumed that all gain changes take place in between saccades. In our models, therefore, the 227 error-based correction terms weight errors that were experienced in previous saccades. As a 228 consequence, in the estimation of the forthcoming event, the post-saccadic stimulus gain is not 229 compared against the adaptation gain measured for that trial but against the previous estimate of 230 the gain. To justify these assumptions, it is usually assumed that the motor system sends an 231 efference copy of the motor command to the sensory areas, which enables prediction of the 232 sensory consequences of the movement and therefore avails comparison to experienced post-233 saccadic feedback [7,19,20,31,[33][34][35].

234
We will assume that the values of saccade and adaptation gains observed and extracted from the 235 recorded data (i.e., supplemented with an initial condition that sets the initial value

‫ܣ‬ ൌ 1
). In the case of the initial value ‫ܩ‬ , we obtained an estimate by taking the average 291 of the first five valued of the gain. We proceeded in this way because the initial value of the state of 292 the system is unknown and, while the first recorded value of the gain could be considered a proxy 293 for such initial state, execution and motor errors could yield a value of the gain significantly 294 different than the actual initial state of the system; we averaged over 5 trials to alleviate this 295 problem. In models where the initial value of the gain was left free to become a fitting parameter, 296 this average over the first five saccades was used as an initial value for the fitting routine for that

328
The adaptation gain of the oculomotor response to a sinusoidal disturbance is best described by a 329 phenomenological function consisting of a decaying exponential added to a lagged but otherwise We use here the same denominations used in our previous report [11], except for changing the including the next-to-last feedback-error term is given by: In models without next-to-last feedback term ‫ܦ‬ should be set to zero; in models that do not have ‫ܣ‬ 356 as a fitting parameter, its value should be set to 1 in Equation 9.

357
The periodic component of the response to a sinusoidal disturbance in models where the next-to-358 last feedback is included can be written as:

360
Equation 10 Equations 9 to 12 clarify the effect of the presence of the next-to-last error learning rate ‫ܦ‬ .

371
Baseline drift parameters 372 Following a sinusoidal disturbance, the baseline of the error gain will approach an asymptote at 373 large trial number that can be written as a function of parameters of Equation 5 as (see also S1 374 Appendix): ߣ for the decay of the baseline, has units of 1/trials and it is defined by: provides the weights of the impulse response that generates the integral solution by

402
Analysis of the data at the phenomenological level

403
To obtain a general idea of patterns present in the data, we first collapsed the data for each 404 stimulus frequency and adaptation type across participants (group data). We fit these data using a 405 piecewise continuous function given by the addition of a monotonic (exponential) decay of the 406 baseline spanning both pre-adaptation and adapting trials-and a periodic entraining of the 407 oculomotor response to the sinusoidal stimulus that begins at the onset of the adaptation block.

408
This choice was supported by the fact that we had confirmed using statistical model selection

437
Some features are readily apparent from these plots. First, the frequency of the ISS is reliably 438 estimated (cf . Fig 3b,d). Second, the amplitude and the lag of the periodic components of the 439 adaptation gain decay with increasing frequency of the stimulus (Fig 2a,b). The amplitudes of the . This may be related to the fact that at such low frequency the stimulus 444 resembles more the behavior of a ramp that then turns rather than a truly periodic disturbance.

445
The parameters that affect the observed drift in the baseline (i.e., asymptote and timescale, Fig   446   2c,d) remain rather independent of the experimental condition. This feature is more apparent in the 447 ORIG dataset, but it still seems to hold in the FREQ dataset. An exception arises at the lower Page 20 of 48 frequency (1 cpb) tested in the FREQ dataset. However, the case of frequency one is rather 449 special and should possibly be considered as transitional between periodic and non-periodic

486
Because the best fitted model may differ between individuals, we first computed the AIC weights     that experiment was rather narrow, ranging from 3 to 6 cpb. This segregation in the learning rates 531 between Two-way and Global adaptation is also clearly present in the best models fitted to dataset 532 FREQ.

533
A feature observed in all cases is that in models that learn only from the last experienced error, the 534 (single) learning rate ‫)ܭ(‬ shows a mild increase with the frequency. This changes substantially if 535 learning from the next-to-last feedback is included. In all of these models, the following features are 536 observed. First, the magnitude of ‫ܭ‬ , the learning rate of the last-feedback error-correction term 537 increases by about an order of magnitude with respect to the models that do not have next-to-last 538 error-correction. Second, the magnitude of the next-to-last error learning rate ‫)ܦ(‬ is similar to that Page 23 of 48 of the last error ‫)ܭ(‬ but with opposite sign. This seems to suggest that the next-to-last error is 540 weighted negatively (or actively attempted to be forgotten) in the algorithm. Third, the discrepancy 541 in magnitude between ‫ܭ‬ and ‫ܦ‬ is consistently larger for Two-way than for Global adaptation 542 (compare the separation between corresponding blue and red lines in Fig 7, a and b).

551
The values of the parameters fitted with the best four models are shown in S1 Table (     the generative models (Equations 9 through 14). We attempt matching these predictions to the 589 values fitted using the phenomenological parameter estimation implemented before (see Fig 2). It

598
We start with Equation 13 that provides a relationship between the expected asymptote of the 599 adaptation gain at large trial number and the generative model parameters.

‫ܤ‬
A first significant observation about this expression is that in order to observe a drift in the baseline

656
Fig S1 in S1 Appendix). This ordering relation between the timescales of models with and without 657 ‫ܦ‬ was unknown before conducting the fits. Thus, data collected using a sinusoidal adaptation 658 paradigm suggests that including a second error-correction term yields a significant decrease in 659 the timescale with respect to models featuring a single error-correction term. Therefore, the 660 integration window (i.e., the inverse of the timescale) of models with double error-correction grow 661 significantly larger compared to those that lack the second error sampling.

696
The oculomotor response, quantified by the adaptation gain, followed the disturbance variation with 697 comparable frequency, an amplitude ranging between 10 and 30% of that of the stimulus (i.e., 2.5 698 to 7.5% of the saccade amplitude), and lagging the stimulus by a few trials. In addition, it 699 developed a systematic drift of the baseline towards larger hypometria that reached asymptotes of 700 around 40% of the disturbance amplitude (i.e., 10% of the saccade amplitude) and was largely 701 comparable across conditions. The phenomenological description in Equation 8-composed of a 702 periodic response and an exponential decay-captured this behavior well and we estimated all six 703 parameters pertaining to that description.

704
The present study explored whether the phenomenology described by Equation 8

725
We analyzed 16 models that differed in the specific parameters that were fitted and then used

726
Akaike's information criteria to attempt model selection. Since we were primarily modeling intrinsic conditions. Finally, we considered the plausibility and study the effects of a second learning rate ‫ܦ‬ 739 that tracks errors other than the most recent (here, the next-to-last feedback error).

740
To further discuss the effect of the generative parameters, we split the 16 models into four

747
We recall that in models where

757
The fits of the phenomenological model (Fig 1; Equation 8) suggest that asymptotic behavior of

763
Data from Two-way adaptation is depicted with blue tones in bars increasing towards the right.

764
Global adaptation is shown with bars spanning to the left in red tones. The average weight for each

814
Page 33 of 48 reducing the utility of maintaining a saccade gain close to one. We note that this systematic 815 decrease in saccade gain may in general-albeit to different degrees-pervade the study of 816 saccadic adaptation (but see [7,54]

819
On the other hand, from the point of view of the internal model of the movement that the brain may 820 implement [33][34][35], this bias parameter ݉ may hint to a discrepancy between the experimental

831
The models that best explained the data featured a double error sampling, learning not only from 832 the feedback experienced after the last saccade but also from the movement that occurred in a trial 833 before that. Hence, the best models used a feedback reaching further back in time through the

865
Using this double error sampling,  remained relatively independent of, or exhibited a tendency to grow with, the frequency of the 916 stimulus (Fig 6a). Learning rates for Two-way adaptation roughly ranged between 0.01 and 0.035 Page 37 of 48 fraction of the error across the frequencies tested. The same parameter in Global adaptation was 918 smaller and remained within the range 0.005 to 0.015 (cf . Fig 6a and S1 Table). These in that type of adaptation was always a significant factor while ISS frequency 921 never was (S2 Table). when the frequency increases because its size from one trial to the next changes faster. However,

930
this suggestion seems at odds with the intuition that a more consistent stimulus should drive more 931 efficient adaptation [68,69]. In particular, it has been reported that a smooth gradual variation 932 results in more efficient adaptation [3,70]. If this were the case and reflected onto the model 933 parameters, the learning rate should be higher for smaller frequencies.

934
However, the dependence of the learning rate(s) on the frequency described above changed rather 935 dramatically when double error sampling was included (cf. Fig 7, columns a and b) Table).

953
In contradistinction, the retention rate ‫ܣ‬ (Figs 6b and 7c) and the bias parameter ; see also corresponding entries in S1 Table). ANOVAs run over these parameters shifts. They found that a linear dynamical system (LDS) with a single error term and trial-by-trial 990 state update for variability implemented with an estimation-maximization algorithm successfully 991 described mean reach point and the temporal correlation between reaching errors and visual shifts.

992
They further argued that the learning taking place under a random stimulus generalizes to a 993 situation of constant shifts in a block paradigm and, therefore, that adaptation dynamics does not 994 rely on the sequence (or correlation) of feedback shifts but can be generally described with the 995 LDS model. In contrast to random or block constant ISS, our paradigm featured a disturbance that 996 was fully self-correlated since it followed a sine variation with the trial number. Therefore, it may 997 prove advantageous for the oculomotor system to extract correlations embedded in the