Supplementary information on the computational model of optimal transsaccadic integration
The model receives three input variables, representing neural signals in the brain: r is the change in the target’s retinal location, before and after the saccade; v is an independent estimate of that change in retinal location derived from visual velocity detectors; c is the estimated change in eye position, conveyed by efference copy or proprioception. These signals are influenced by two other variables s and j, representing saccades and jumps, i.e. the amplitudes and directions of sudden movements of objects in the visual field (an illustration of the five variables is given in Fig. Suppl1). The model deduces j*, the optimal estimate of j, from c, r and v. This j* is not necessarily the most likely value of j given c, r and v; it is the estimate that will yield the smallest misjudgement of j on average. In other words, the model computes the j* that minimizes the integral
dj (j – j*)2 p(j|crv)(1)
where the range of integration covers all possible values of j, and where p(j|crv) is the conditional probability density function of j given c, r and v.
Derivation of the optimal estimate
As the integral (1) is differentiable with respect to j*, its minimum will occur at a point where its partial derivative with respect to j* equals 0:
dj (j – j*) p(j|crv) = 0.(2)
Rearranging terms, we find that
j* = dj j p(j|crv) / dj p(j|crv).(3)
In practice, we approximate this formula numerically, replacing the integrals by discrete sums and the probability density p(j|crv) by the discrete probability P(j|crv):
j* = j jP(j|crv) / j P(j|crv).(4)
In effect, we pretend that our variables s, j, c, r and v are discrete, with a finite number of possible values. If we consider a lot of discrete values spaced closed together, the sums in (4) closely approximate the integrals in (3). (Passing between probabilities and probability densities also requires some simple scaling, detailed in the computer code provided below.)
The unknown quantity in (4) is P(j|crv). What we want, then, is a formula for computing the conditional probability
P(j|crv) from P(s), P(j), P(c|s), P(r|js), and P(v|js)(5)
The last five variables are simpler quantities which the visuomotor system could plausibly already know; they might be built into the brain genetically or learned by experience.
Our derivation of P(j|crv) relies on three simple laws of probability. One is Bayes’ rule, that for any events a and b, P(ab) = P(a)P(b|a), i.e. the probability of a and b both occurring is the probability of a times the conditional probability of b given a. (Some Bayesians consider all probabilities to be conditional; instead of P(a), for instance, they write P(a|h) where h is some background hypothesis, but in the present derivation this is merely a notational variant.) The second simple law, a special case of Bayes', is the product rule, which says that if events a and b are independent, then P(ab) = P(a)P(b). The third law, which we will call additivity, is that if b1, b2 ... bn are mutually exclusive events whose summed probability P(b1) + ... + P(bn) equals 1, then for any event a, P(a) = P(ab1) + ... + P(abn) = P(a|b1)P(b1) + ... + P(a|bn)P(bn).
By Bayes' rule,
P(j|crv) = P(jcrv) / P(crv)(6)
Computing first the denominator of this quotient, we have, by additivity and Bayes' rule,
P(crv) = the sum, over all possible values of j and s, of P(crv|js)P(j)P(s).(7)
And by the product rule,
P(crv|js) = P(c|js)P(r|js)P(v|js)(8)
This equation is valid because c, r and v, which are carried on physically distinct channels in the brain, are independent signals when j and s are given. (Obviously they are not independent when either j or s is not given, because, for instance, both r and v depend on s, so learning v would help you guess r. That is why, in our derivation, we could not simply write P(crv) = P(c)P(r)P(v).) As c does not depend on j, we can simplify (8) to
P(crv|js) = P(c|s)P(r|js)P(v|js)(9)
And with (7) and (9) we can compute the denominator in (6), given just the five known quantities in (5).
It remains only to compute the numerator in (6). By Bayes' rule and additivity,
P(jcrv) = s P(crv|js)P(j)P(s)(10)
We now have all we need to compute the optimal estimate, j*, from our five known quantities: Eqs. (9) and (10) yield P(jcrv); (7) and (9) yield P(crv); (6) says that the ratio of these two quantities is P(j|crv); and from P(j|crv), (4) yields j*. Collected into one equation, the resulting discretized formula for j* in terms of the five known probabilities in (5) is
j* = j,s jP(j)P(s)P(c|s)P(r|js)P(v|js) / j,s P(j)P(s)P(c|s)P(r|js)P(v|js).(11)
The continuous version of this equation appears in the Methods section of the paper.
Probabilities distributions
The probability density functions p(s), p(c|s), p(r|js) and p(v|js) were approximated using gaussians of the form
G(x, μ, σ) = exp(–((x – μ)/σ)2 / 2) / ((2π)1/2σ)(12)
The probability density of saccades is given by
p(s) = G(s, μs, σs).(13)
Here, saccade scatter σs is an affine (i.e. first-order, linear-plus-bias) function of the distance to the target:
σs = as t + bs.(14)
where t is expressed in degrees. For the scatter orthogonal to the saccade we set as and bs to 0.01 and 0.05˚ – values derived from van Opstal and van Gisbergen29. To simulate scatter parallel to the saccade we multiplied the orthogonal scatter values by the eye-position scatter ratios 2.5, 3, 3.5, 4 and 4.5 (cf. Fig. 2a). The variable μs – the mean saccade amplitude on the task at hand – equals the distance to the target minus the average undershoot, u, which we found, based on our own data, to be well approximated by a function of the parallel saccade scatter σs:
u = au σs + bu.(15)
where au and bu are 0.23 and 0.41˚.
For the signal c, which codes changes in eye position, the probability density is
p(c|s) = G(s, c, σc).(16)
As we found a close relation between the brain’s uncertainty about its signals and saccade scatter (as we describe in the final part of the Methods section of the paper), we approximated σc with σs; in reality these variables may not be so tightly related, but at least this assumption does not introduce extra flexibility into the model, but on the contrary removes a free parameter.
The variable c represents the perceptual system’s estimate of the actual saccade. Immediately after the saccade, this signal is delayed and therefore hypometric, i.e., the eye movement is underestimated22. The reason for this delay is not relevant here (it may arise because the brain estimates eye position by combining efference copy with later proprioceptive signals), but its effect is to make c smaller than s by a factor which we call del. Del varies with saccade size, between about 0.75 and 0.9.
The signal r represents changes in the retinal locations of the target image, so it depends on both s and j. Its probability density is
p(r|js) = G(j – s, r, σr)(17)
where j and s represent the components of the saccade and the target jump in the same dimension as the retinal change. For instance, the mean horizontal change in retinal-image location when the target jumps 2º right and the eye saccades 20º right is 2 – 20 = –18º right, or 18º left. When the saccade is 20º rightward and the jump is 2º up, then the jump has a zero component parallel to the saccade, and so the mean horizontal change in retinal-image location is 0 – 20 = –20º right, or 20º left; and by similar reasoning the mean vertical change in retinal location is 2 – 0 = 2º up. As for σr, it should increase with the retinal eccentricity of the target image because retinal localization declines with retinal eccentricity9,10. But in our experiments, subjects estimated changes in retinal location by comparing images at different times, so it is hard to extract a quantitative value for σr from published data. We chose a simple function,
σr = ar ecc + br(18)
where ecc is the postsaccadic or postjump eccentricity of the target image. Replacing (18) by moderately non-affine functions made little difference to the model's predictions. In our experiments, initial eccentricity was constant, mostly at 15º, so its influence can be expressed in the constant term br. ar and br were largely free parameters, constrained only by reasonableness; in the simulations in Figures 1 – 3 of the paper they are set to 0.02 and 0.15º. Hence, br is smaller than ar times 15º, which means that the initial eccentricity in br contributes less to σr than postsaccadic eccentricity. This is realistic because the model (like our subjects) performed many trials with the same presaccadic target eccentricity, so after the first few trials, that eccentricity could be estimated with confidence.
Here arises an important conceptual point. All these probability densities describe physical processes in sensors or motor systems, but we are also proposing that the same densities are accurately represented in the central neural processors responsible for transsaccadic integration. In other words, the part of the brain that carries out the computation in Eq. (11) 'knows' the correct probability density functions of its sensors and motor systems – this is part of what it means to have optimal integration. But in the case of p(r|js), the central processors have no direct access to the variable ecc in Eq. (18). Instead, ecc must be estimated from the initial target distance and from the available information about saccade amplitude and retinal changes. In our simulations, the resulting estimates of ecc and σr were close to the real values and yielded almost identical predictions.
The probability density of the velocity signal is
p(v|js) = G(v, μv, σv).(19)
We are interested in v as a source of information about retinal displacement. Of course the brain has to compute that displacement from velocity and time, but in our optimized model we assume that the computation is performed correctly, so we simply set μv = r; i.e. we assume that the brain deduces retinal displacement accurately, on average, from its velocity signals, when there is no saccade. During saccades, motion perception is markedly reduced or even turned off13-16, so we set μv = 0 (though the optimal integrator knows that v is uninformative during saccades, so it ignores that signal entirely, and therefore the value of μv actually plays no role). For σv we relied on data from Tynan and Sekuler12 who showed that the threshold for motion perception varies only slightly around 0.1º/s within about 15º of the fovea. For our gaussian, this threshold corresponds to a σv of about 0.005º. But the precise value turns out not to matter greatly; the important point is that σv is much smaller than σc and σr.
The probability density p(j) represents the central neural processor's assumptions about the likelihood of target jumps in the presence or absence of saccades, based on innate knowledge and the visuomotor experience of a lifetime. Little guidance regarding p(j) is available from the literature. We reasoned that the function should slope up steeply to a peak at zero, because objects that are otherwise stationary seldom jump abruptly. We chose the Laplace distribution, which applies to rare events whose probabilities vary across a range of different situations. We set
p(j) = exp(–|j| / w) / 2w.(20)
Here w is a free, width parameter which we set to 0.065º. To simulate neural variability in the integrator itself, we added gaussian noise of standard deviation 0.004 radians to j*. Other values for any of these parameters yield the same qualitative predictions but change the quantitative details, such as the slope in Figure 2a.
Code
Below is the computer code for a simple simulation that yields estimates of target jumps, as in Fig. 3e of the paper. The code is written in Delphi Object Pascal because for our computations that language yielded programs at least twice as fast as any of the usual alternatives, and speed is critical for our larger simulations where the model is put through lengthy perceptual experiments, as in Fig. 1 of the paper. As Delphi closely resembles Borland-Inprise C++, the code should be comprehensible to most programmers. For this sample program we chose a simple perceptual task – judging jumps parallel to the saccade – so that the essential code dealing with the transsaccadic integrator would not be obscured amidst other instructions. With this program we calculated the sensitivity curves in Fig. 3g-j.
procedure TForm1.jstar1Click(Sender: TObject);
//Optimal transsaccadic integration: estimating target jumps in one dimension
var
a_c, a_u, a_r, b_c, b_r, b_u, c, cSD, del, denom, dj, dr, ds, ecc: double;
j, jStar, jump, w, numer, p, pc_s, pj, pjscrv, pr_js, ps, pv_js, r, rSD: double;
s, sac, sSD, stim, targ, u, v, vSD, xScale, yScale: double;
ji, jn, ri, rn, si, sn, ti, tn, xZero, yZero: integer;
colour: array[0..5] of TColor;
begin
tn := 2; //# of saccade targets = 2*tn + 1
jn := 100; sn := jn; //larger jn and sn give slower, more exact integration
rn := 20; //# of nodes in each plotted curve = 2*rn + 1
dr := 10*pi/(180*rn); dj := 0.5/jn; ds := dj; //increments
//Set up graphics
colour[0] := clWhite; colour[1] := clRed; colour[2] := clYellow;
colour[3] := clLime; colour[4] := clBlue; colour[5] := clPurple;
Canvas.Brush.Color := clBlack;
Canvas.Font.Color := clWhite;
Canvas.FillRect(Rect(0, 0, ClientWidth, ClientHeight));
Canvas.Pen.Color := clGray;
xZero := ClientWidth div 2; xScale := 0.5/(rn*dr)*ClientWidth;
yZero := ClientHeight div 2; yScale := -0.5/(rn*dr)*ClientHeight;
Canvas.Pen.Color := clWhite;
Canvas.MoveTo(round(xZero - xScale), yZero);
Canvas.LineTo(round(xZero + xScale), yZero);
Canvas.MoveTo(round(xZero - xScale), round(yZero - yScale));
Canvas.LineTo(round(xZero + xScale), round(yZero + yScale));
//Set parameters
a_c := 0.01; b_c := 0.05*pi/180; //in radians
a_u := 0.23; b_u := 0.41*pi/180; //in radians
a_r := 0.02; b_r := 0.15*pi/180; //in radians
vSD := 0.005*pi/180; //in radians
w := 0.065*pi/180; //in radians
del := 0.75;
//Test transsaccadic integration
for ti := -tn to tn do begin //different saccade targets
//For this simulation, we consider typical saccades with typical eff copy, so we set c = del*sac,
//because c does equal sac on average - except for some delay in the signal.
//The neural estimator, below, that computes j* knows c, but
//can't assume that c = del*sac, so it must consider many possible s's.
targ := ti*7.5/180*pi; //saccade target in radians
sSD := a_c * abs(targ) + b_c; //motor variability
if ti = 0 then begin //stimulus jumps between saccades
stim := 15/180*pi; //during fixation stimulus distance set to 15 deg
u := 1; //‘undershoot’ is total when no saccade is made
end
else begin //stimulus jumps during saccade
stim := targ; //the saccade is aimed at stimulus
sSD := 2.5*sSD; //saccade scatter is larger parallel to the saccade
u := (a_u*sSD + b_u) / abs(targ); //proportional undershoot rises linearly with saccade scatter
end;
cSD := sSD; //variability in efference copy
sac := targ * (1 - u); //actual saccade undershoots the target
c := sac*del; //immediately after the saccade the efference copy is hypometric
Canvas.Pen.Color := colour[ti + 3];
for ri := -rn to rn do begin //different target jumps
jump := ri*dr; //different jump sizes are tested
r := jump - sac; //retinal-image displacement in radians
if ti = 0 then v := r else v := 0; //velocity signal is accurate when there is no saccade
ecc := abs(stim+r); //postsaccadic stimulus eccentricity as estimated by the system
//the true ecc = abs(u*stim + jump) is unknown, as jump is unknown to the system
rSD := ecc*a_r + b_r; //retinal acuity
//Compute jStar
numer := 0; denom := 0;
for ji := -jn to jn do begin //possible jumps entertained by estimator
j := ji*dj; //possible jump in radians
pjscrv := 0;
for si := -sn to sn do begin //possible saccades entertained by estimator
s := si*ds; //possible saccade in radians
//In the following gaussians, we ignore the constant scalings by 1/(SD * sq rt 2pi)
//because they cancel out in j*
p := s/sSD; ps := exp(-0.5*p*p)/sSD; //probability of this saccade
p := (c - s - sac)/cSD; pc_s := exp(-0.5*p*p)/cSD; //cond'l prob of c given s
p := (r - j + sac + s)/rSD; pr_js := exp(-0.5*p*p)/rSD; //cond'l prob of r given j & s
if si=0 then p := (v - j + sac + s)/vSD //between saccades, avg v = r
else p := v/vSD; //during saccades, avg v = 0
pv_js := exp(-0.5*p*p)/vSD; //cond'l prob of v given j & s
pj := exp(-abs(j)/w)/(2*w); //prior probability of jump
pjscrv := pjscrv + pj*ps*pc_s*pr_js*pv_js;
end; //si
numer := numer + j*pjscrv;
denom := denom + pjscrv;
end; //ji
jStar := numer/denom;
//Plot
if ri=-rn then Canvas.MoveTo(round(xZero + xScale*jump), round(yZero + yScale*jStar))
else Canvas.LineTo(round(xZero + xScale*jump), round(yZero + yScale*jStar));
end; //ri
end; //ti
end;