Continuous and Discrete Volterra-Laguerre Models With Delay for Modeling of Smooth Pursuit Eye Movements

The mathematical modeling of the human smooth pursuit system from eye-tracking data is considered. Recently developed algorithms for the estimation of Volterra-Laguerre (VL) models with explicit time delay are applied in continuous and discrete time formulations to experimental data collected from Parkinsonian patients in different medication states and healthy controls. The discrete VL model with an explicit time delay and the method for its estimation are first introduced in this paper. The estimated parameters of a second-order VL model are shown to capture the ocular dynamics both in health and disease. The possibility of including the estimated time delay, along with the VL kernel parameters, into the set of the model parameters is explored. The results obtained in continuous VL modeling are compared with those in discrete time to discern the effects due to the sampling enforced by the eye tracker used for data acquisition.


I. INTRODUCTION
V OLTERRA models (VMs) offer a black-box framework for modeling smooth nonlinear dynamical systems [1]. Typically, the kernels of a VM are expanded in an orthonormal basis to obtain a parsimonious and linear in parameters characterization. When the kernels are sufficiently smooth and vanish at infinity, the Laguerre functions constitute a suitable basis, a choice resulting in Volterra-Laguerre (VL) models that are available in both continuous and discrete formulations [1]- [4].
Since VL models are linear in the coefficients of the kernel series, linear estimation techniques can be used for system identification. Complete VL models contain many unknown parameters and demand a high degree of system excitation both in amplitude and frequency to guarantee model identifiability. Therefore, sparse estimation methods such as LASSO [5], [6] and SPICE [7] are applied to identify and estimate only the essential kernel components. To further promote model The authors are with the Information Technology, Uppsala University, Sweden (e-mail: alexander.medvedev@it.uu.se).
Digital Object Identifier 10.1109/TBME.2022.3185669 parsimonicity, an approach combining sparse estimation with estimated functional relations between the VL coefficients has been proposed and applied to the modeling of the human smooth pursuit system (SPS) [8].
The SPS enables tracking moving targets with gaze and involves multiple parts, including the eyes and the extraocular muscles, the optic nerves, and multiple regions of the brain [9]. A dysfunction in any of these parts will thus impair the system as a whole. Alcohol intoxication [10], schizophrenia [11], [12] and Parkinson's disease (PD) [13]- [18] are known to reduce the SPS performance. Thus, mathematical modeling of the SPS may be instrumental in diagnosis of these conditions or symptoms quantification.
Eye tracking is the technology to measure eye movements and has been practiced in various forms since at least the beginning of the 20th century [19]. However, modern video-based eye tracking has allowed inexpensive non-intrusive methods. Here, camera recordings of the eyes are used, and the gaze direction can be inferred from the reflections of infrared illumination in the cornea and lens. Applications of this type of eye tracking can be found in e.g. video games [20], user experience evaluation [21], and marketing research [22].
When a test subject, whose eye movements are measured by an eye tracker, follows presented and designed in advance visual stimuli, the experimental setup fits into the framework of system identification. Then, the problem of mathematical modeling from recorded input-output data consists of model structure selection and model parameter estimation. Since the neural mechanism of the SPS control is not completely known [9], a black-box model is required. Further, to capture the SPS dynamics both in health and disease, the model has to be nonlinear. VM [23] can be employed with minimal assumptions regarding the dynamical properties of the modeled system. With enforced sparsity, quadratic VMs were successfully applied to the problem of SPS modeling [24].
The SPS possesses an intrinsic delay due to the finite speed of neural communication. The digital signal processing of the eye tracker contributes more lag to the gaze direction measurement. As the kernels of a VM are allowed to be infinite dimensional, a measurement delay is readily (albeit implicitly) accounted for by an estimation algorithm. Yet, when the value of the delay is of interest, it has to be explicitly parameterized in the model, which case is seldom addressed in the literature. For continuous VL models with time delay, the problem of estimating the kernels is solved in [25]. A corresponding discrete-time formulation is introduced for the first time in the present paper hinging on recently obtained Laguerre-domain representation of discrete linear time-invariant systems with delay [26]. Laguerre functions constitute a special case of Meixner functions that have been applied to the estimation of discrete VM with time delay [27]. This basis generalization is suitable when Volterra kernels exhibit a slow initial onset. Yet, the problem of delay estimation is not considered in [27].
The main contribution of this article is in a comparison of VL model types with respect to the problem of modeling eyetracking data in order to distinguish between PD patients and healthy controls. Discrete and continuous minimal VL models with and without explicit time delay estimation are treated. All the considered model structures show significant differences in the linear terms of the estimated VL models. The estimated time delay value is generally higher for the PD patients and exhibits more variance than that for the controls. In fact, the present study is the first one to attempt estimating the delay in the SP mechanism itself and not in its initiation. The distributions of the quadratic model term estimates also exhibit dissimilarities, but to a lesser extent. A methodological contribution is the formulation of the discrete VL model with explicit time delay that parallels the continuous-time case.
The rest of the paper is organized as follows. First Laguerre functions and the VL models with and without explicit time delay are introduced. Then some background on the human SPS and the eye tracking experiments is provided to explain the data used in the study. Furthermore, discrete and continuous VL models of the eye tracking data, with and without the time delay as a model parameter, are estimated and analyzed.
where p c is the continuous Laguerre parameter. The functions k (s), k = {0, ∞} constitute a complete orthonormal basis in H 2 with respect to the inner product The k-th Laguerre coefficient of W (s) ∈ H 2 is evaluated as a projection of W (s) onto k (s) where the set {w k } k=0,1,... is referred to as the Laguerre spectrum of W (s).
where p d ∈ (0, 1) is the discrete Laguerre parameter. The corresponding time-domain functions The discrete Laguerre spectrum of f ∈ 2 [0, ∞) is then the discrete sequence {f j } j=0,1,... obtained by calculating In what follows, the same notation for the products in (2) and (4) is utilized, and the difference is implied by the considered modeling framework. Notice that shifting the Laguerre functions in time does not change their orthonormality. To apply timeshifted Laguerre functions to the estimation of kernels with time delay requires the time-delay value be known.

1) Continuous Time: A general continuous VM is written as
and is given by the kernels {k 1 (·), k 2 (·), . . .}. Defining the continuous Volterra functional of degree n as model (6) can be written in the compact form of Now consider an input signal u with a delay τ , i.e.
where H(t) is the Heaviside step function. Following [25], the system may then be written as Thus, the delay can as well be attributed to the kernels.

2) Discrete Time:
The discrete model matching (6) is Similarly, the discrete Volterra functionals are defined as and the discrete model is then written as C. The Volterra-Laguerre Model

1) Continuous Time:
Assuming kernel separability and expressing the n-th kernel k n in the basis functions yields where the coefficients γ n (j 1 , . . . , j n ) are referred to as the VL coefficients. Now consider the convolution integral between the input signal u(t) and the j-th Laguerre function The function ψ j (t) can be seen as the output signal of a filter whose transfer function is j (s) and that is driven by the input u(t), i.e. the Laguerre filter output. Using (15), combined with (8), the Volterra functionals in the Laguerre basis become The full VL model corresponding to (6) is then (17) 2) Discrete Time: As in the continuous case of (14), the n-th separable discrete kernel in terms of Laguerre functions is Denoting the convolution between the input u(k) and the jth discrete Laguerre function, i.e. the discrete Laguerre filter output, as the n-th discrete Volterra functional can be written as The full discrete VL model corresponding to (10) is then

D. The VL Model With Time Delay
A general VL model may readily incorporate time delay in the kernels. The continuous kernels are square-integrable functions but not necessarily represent finite-dimensional operators. In discrete time, the kernels have to be square-summable and remain finite dimensional despite a time delay. 1) Continuous Time: Consider a second-order VL model where ψ (τ ) From [25], the same output is obtained from where γ (τ ) n denote the coefficients of a VL model with implicit delay. These coefficients are obtained by (24) and γ (τ ) Find the VL coefficientsγ using ordinary least squares, or using some sparse estimation method (e.g. [7], [30]). Notice that, given τ i , the VL model is linear in the unknown parameters. 6: Compute the output loss ε (τ i ) = y −ŷ(τ i ,γ) 2 7: end for is the associated Laguerre polynomial [28] defined by 2) Discrete Time: Similarly to the continuous case, a VL model with explicit delay as well as one with implicit delay can be formulated. For coefficients in the linear part, the relation between the explicit and implicit case is where The properties of the polynomials Ł (τ ) m (·) are studied in detail in [29]. Assuming kernel separability leads to

E. Estimation of VL Models With Time Delay
Estimation was performed on the individual data set and on a population level. Estimation of VL parameters for a single data set is performed as described in Algorithm 1.
The estimation procedure on the population level is described in [8] and can be summarized as in Algorithm 2.

III. THE SMOOTH PURSUIT SYSTEM
Smooth pursuit (SP) is the type of eye movement arising when the gaze tracks a moving object. The SP movement is initiated and driven by the continuous movement of the tracked target, the visual stimulus. In SP, the angular velocity of the eye is matched with that of the stimulus [31]. The maximum angular velocity that the human SP system (SPS) is able to handle is around 80 − 100 • /s [32]. Above this, the saccadic system takes over the gaze control, to catch up with the moving target in a ballistic manner. Previous studies on delay in SPS have focused mainly on the movement onset, i.e. the delay between the start of target movement and the start of SP. Simpler stimuli such as velocity steps and ramps are usually employed to facilitate delay estimation. The onset delay reported in these studies are typically in the range of 0.1-0.2 s [33]- [35].

A. Eye Tracking Data
The eye tracking data used in this paper were collected at Uppsala University hospital by Clinical Trial Consultants AB (CTC) in 2015. The study was a part of the project "MuSyQ: Multimodal motor symptoms quantification platform for individualized Parkinson's disease treatment" funded by the Swedish Innovation Office (Vinnova), and a comprehensive description of the performed experiments is found in [36].
A total number of 144 data sets were collected from a group of seven patients diagnosed with Parkinson's Disease (PD), and 499 data sets were collected from 22 healthy, age-matched control subjects. A video-based desktop eye tracker from SmartEye AB was used for acquiring the data.
In order to evaluate response to levodopa (LD) treatment throughout the day, the PD patients followed a protocol. After an eight-hour long washout period, the patients received 150% of their usual morning dose of LD. Subsequent to drug administration, the patients performed eye-tracking and motor tests in 20 min intervals. The experiment design aimed at capturing the transition in the patient state from Parkinsonism (off), to symptom-free (on), further, probably, to a state of dyskinesia due to the exaggerated dose, and back to the off state, when the therapeutical effect of levodopa ends. To quantify SP performance, a dot moving on the computer screen according a pre-defined optimized pattern (see [15] for the stimuli design algorithm) was presented to the test person. It is developed to possess excitation properties in both frequency and amplitude to reveal possible nonlinear oculomotor dynamics. An example of the horizontal part of the stimulus is shown in Fig. 1.
During the eye-tracking experiment, the gaze position (the screen coordinates at which the gaze is directed) was recorded, resulting in an input-output pair of two-dimensional signals reflecting the horizontal and vertical coordinates of both the stimulus (the moving dot) u, and the response (the gaze position) y. The duration of the stimulus is ca. 25 s and the sampling frequency of the eye tracker is 60 Hz, thus yielding data sets of about 1500 samples.

IV. RESULTS
The discrete and continuous VL model structures were chosen by means of the procedure described in [8]. Laguerre and nonlinearity orders N = 2, L = 2 were used initially, and in the explicit delay case the delay τ was also considered as a parameter. The models were subsequently reduced by sparse estimation and linear approximation of relations between the VL coefficients estimated across population.

A. Discrete Time
The following model structure was estimated in both the implicit and explicit delay case with p d = 0.4 where the factors θ Thus, four coefficients (denoted by γ i (·), i = 0, 1, 2) are to be estimated from a given data set.  The obtained values of the linear coefficient γ 1 (0) for the explicit and implicit delay cases, as well as the estimated delay values in the explicit delay case, are shown in Fig. 2. In the case of explicit delay model, the delay estimate was limited to 10 samples (0.1667 s), at most. The estimates of the quadratic coefficients are shown in Fig. 3. There is a clear distinction between the histograms of the linear coefficients' estimates for patient and control groups, while the distributions of the quadratic coefficient are very similar and apparently independent of the delay parametrization.
The distribution of parameters in the quadratic part seems centered around zero. Since the data include records collected from both healthy controls and Parkinsonian patients, this is not unexpected. In Fig. 4, the uncertainty of the quadratic parameter estimates from one patient data set is shown to demonstrate that while the population distribution is centered around zero, the individual estimates are non-zero when the symptoms are significant.
Examples of the kernel functions, including the estimated delay, are given in Fig. 5. Note that while the second order kernel    is a two-dimensional function h 2 (i 1 , i 2 ), the plot shows the one-dimensional function h 2 (i, i). This is however equivalent to the Volterra functional (H 2 δ)(k), where δ denotes the discrete impulse function. Furthermore, an example of the stimulus and measured gaze position (input and output signals) is provided in Fig. 6, together with the simulated model output. Fig. 7 illustrates the model output error for a given data set and a wide range of p d . The model quality is evidently robust with respect to the choice of the Laguerre parameter.

B. Continuous Time
A continuous VL model with L = 2, N = 2 was used with the Laguerre parameter p c = 15. It was established by sparse estimation that, in the case of the implicit delay model, the non-zero parameters occurring in a majority of the cases were γ 1 (0), γ 1 (1) and γ 2 (0, 0), γ 2 (2, 2). In the explicit delay case the most occurring non-zero parameters were γ 1 (0), γ 1 (2) and γ 2 (0, 0), γ 2 (2, 2). That is, the implicit delay model is and the explicit delay model is Introducing the delay as a parameter changes the linear part of the estimated model significantly in this case. Recall that all kernels of the VL model are delayed equally. The coefficients were estimated as follows. Continuous-time Laguerre filter outputs were computed through the three-term formula which is obtained from (1). These were subsequently sampled, and the VL parameters were estimated according to Algorithm 1. The linear coefficient estimates are shown in Fig. 8, together with the estimated delay in the explicit delay model. As in the discrete case, a constraint was put on the delay estimate. The estimated continuous delay is constrained to be under 0.166 s corresponding to 10 samples in the discrete case. The coefficients of the quadratic part are shown in Fig. 9. In contrast to the discrete case, there was no linear relation is found between the coefficients of the linear kernel. However, dissimilarities in the estimate distributions evaluated through histograms may be observed here as well, as in the discrete time case. While the distribution of the estimated delay for the patients  is shifted to the right compared to the controls, the difference is not as obvious as in the discrete case.
Accuracy of the four considered model types is illustrated in Fig. 10. Although the differences between the models are small, the error in continuous modeling is slightly smaller.

V. DISCUSSION
The differences arising in continuous and discrete VL modeling of the human SPS were investigated, making use of models with both implicit and explicit time delay. The experimental data for model estimation comprise recorded eye movements of PD patients and healthy controls. In order to highlight the distinction between the ocular dynamics in health and disease, parsimonious models describing the data with a minimal set of free parameters are sought for. A minimal structure for the VL models was established by estimation on a population level. Thereafter, coefficients were estimated and analyzed for individual data sets.
Second-order VL models capture the oculodynamics of the SPS well both in the patients and controls. This is in line with previous research on VL modeling of the SPS [8], [24], [25]. Therefore, higher-order nonlinearities appear to be insignificant in the oculodynamics. Even in the PD patients, where nonlinearities are expected to be significant, the estimated quadratic term is small compared to the linear one, as Fig. 3 and Fig. 9 demonstrate. This relationship holds notwithstanding the way the delay is included in the model and the considered time frame. As Fig. 3 and Fig. 9 show, the quadratic parts of the estimated VL models do not change much with the introduction of delay as an explicit parameter. While the individual parameter estimates are different, the distributions are highly similar.
The parameter distributions, including the distribution of the delay parameter in the explicit delay models, suggest that the kernel estimates change both with the introduction of the delay as a parameter, and differ between discrete and continuous time formulations. Therefore, the delay is indeed essential in the VL model structure. Notice that the delay value around 0.1 − 0.2 s that is usually mentioned with respect to the SPS is the one in initiation of SP. This is naturally not the same delay that arises when the test subject is engaged in the process of following the visual stimuli. Generally, the established methods for modeling SP that do not rely on system identification technology are too coarse to capture nonlinearity and delay [14]. This study is the first one that attempts to estimate the delay in smooth pursuit from eye-tracking data. Unfortunately, the delay imposed by the data acquisition and digital processing in the eye tracker is not documented and apparently not constant. This complicates distinguishing between the SP delay specific to the patient and the additional delay introduced by the eye-tracking technology.
In the discrete time setting, similar kernel parametrizations were found in the implicit and explicit delay models, as specified in (29)- (31). However, functional relationship between the Laguerre coefficients of the linear kernel is not the same. In the continuous time case, different kernel parametrizations were found for the implicit and explicit delay models, since different functions are present in the Laguerre spectra of linear kernels. The result is expected since the delay operator is infinite-dimensional in continuous time and significantly changes the dynamics. In discrete time, the introduction of a delay only increases the order of dynamics by adding multiple zero poles.
In the discrete case, the estimated delays are clearly larger in a majority of cases for the patient group compared to the control group, as is seen in Fig. 2. The continuous modeling does not show a difference as clear as the discrete one. This can be explained by the fact that the experimental data are discrete, to start with. Further, the calculations in continuous time are more prone to numerical issues than those in discrete one. However, the histogram in Fig. 8 supports that there is a larger proportion of patients than controls that tend toward higher delay estimate values. It should be kept in mind that the eye-tracking data for the PD patients correspond to different degrees of medication and some data set were registered on patients not showing apparent symptoms [36]. This also explains the obvious overlap in the parameter distributions for the patients and controls.

VI. CONCLUSION
Three main conclusions have been drawn from the results above, which are supported by the results of both the discrete and continuous time modeling. First, the structure of the linear part of the VL model is changed by introducing the delay as a parameter. There is also a clear difference between the patients and controls in the linear parameter estimates. Second, the quadratic part of the VL model does not change much with the introduction of the delay as a parameter. Third, the estimated delay is consistently larger for the patient group than for the control group in both the discrete and continuous time cases. However, the estimated delay values are in general smaller than the commonly reported 0.1 − 0.2 s in studies of SP onset.