=Paper= {{Paper |id=Vol-3909/Paper_36.pdf |storemode=property |title=Analysis of The Accuracy of Simulating the Human Eye Movement System Based on Volterra Models |pdfUrl=https://ceur-ws.org/Vol-3909/Paper_36.pdf |volume=Vol-3909 |authors=Vitaliy Pavlenko,Denys Lukashuk |dblpUrl=https://dblp.org/rec/conf/iti2/PavlenkoL24 }} ==Analysis of The Accuracy of Simulating the Human Eye Movement System Based on Volterra Models== https://ceur-ws.org/Vol-3909/Paper_36.pdf
                                Analysis of the accuracy of simulating the human
                                eye movement system based on Volterra models
                                Vitaliy Pavlenko1,*, , Denys Lukashuk1

                                Odesa Polytechnic National University, Shevchenko Avenue 1, 65044 Odesa, Ukraine
                                1




                                                Abstract
                                                Integral nonlinear models are used to simulate the human eye movement system (EMS) while accounting
                                                for its nonlinear dynamics and inertial properties. Multidimensional transient characteristics (MTCs) of
                                                the EMS were identified based on experimental input-output data obtained from eye-tracking responses
                                                to visual test stimuli. These transient characteristics include first-order and diagonal cross-sections up to
                                                the second and third orders of MTCs. The study aimed to evaluate the accuracy of EMS simulation models
                                                by analyzing the calculation errors of transient characteristics using nonlinear dynamic identification
                                                methods based on Volterra integro-power series (IPS) and integro-power polynomial (IPP). Computational
                                                methods, including the least squares method (LSM), approximation, and compensation, were used to
                                                derive the models. Models developed using the LSM and approximation methods produced consistent
                                                transient characteristics when the same test signals were applied, highlighting the convergence of the
                                                Volterra series within the identified region. The findings showed that increasing the number of test
                                                signals enhanced the accuracy of the EMS models. Quadratic models were identified as the most reliable,
                                                providing a balance between precision and computational efficiency. Cubic models closely matched EMS
                                                responses but exhibited instability in their transient characteristics, making them less practical for EMS
                                                application. The compensation method, while computationally less demanding, proved unsuitable for
                                                tasks requiring high accuracy due to significant errors in the resulting models. Quadratic IPP models
                                                developed with LSM based on three response datasets are recommended for future studies, as they
                                                provide a stable and precise framework for modeling EMS dynamics and exploring psychophysiological
                                                state assessment.

                                                Keywords
                                                eye movement system, simulation, integro-power series and polynomials, multidimensional transient
                                                characteristics, eye-tracking, accuracy of simulation, a neurophysiological condition.1



                                1. Introduction
                                Eye-tracking technology [1] is widely applied in the assessment of neurophysiological conditions
                                [2] [5], cognitive research, and memory studies [6], as well as in monitoring student behavior and
                                learning processes [7]. This technology provides valuable insights into both conscious and
                                subconscious human actions. Understanding eye movements is crucial for expanding research in
                                various professional fields, ultimately enhancing the efficiency of work activities.
                                   Despite the advancements in eye-tracking systems, there is a growing need for new
                                mathematical models to accurately simulate the human eye movement system (EMS). Additionally,
                                specialized equipment is necessary to support experimental research involving these systems. Eye-
                                tracking technology relies on sophisticated devices, known as eye trackers, to precisely determine
                                the coordinates of eye movements.



                                Information Technology and Implementation (IT&I-2024), November 20-21, 2024, Kyiv, Ukraine
                                 Corresponding author.
                                 These authors contributed equally.
                                   pavlenko_vitalij@ukr.net (V. Pavlenko); user111228322@gmail.com (D. Lukashuk)
                                   0000-0002-5655-4171(V. Pavlenko); 0009-0003-3149-3741 (D. Lukashuk)
                                         © 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).




                                                                                                                                                                                   455
CEUR
                  ceur-ws.org
Workshop      ISSN 1613-0073
Proceedings
   In medical applications and experimental psychology, the ability to simulate the EMS is vital for
effective control, monitoring, and diagnostics. Developing a comprehensive mathematical model of
the EMS, which takes into account individual human variability, is essential for creating advanced
treatment methodologies [2]. This includes a wide range of personalized applications, such as
medical and sports simulators, human-machine interface testing [8], secure data access [9], [10],
and more.
   This paper presents the findings on the accuracy of simulating the dynamic characteristics of
the human EMS. These characteristics were derived from experimental "input-output" data,
capturing responses to visual stimuli using cutting-edge eye-tracking technology (a simulation
task). The study focuses on the nonlinear and inertial properties of the EMS, which are crucial for
developing accurate and reliable models.

2. Problem statement
To simulate the human eye movement system, integral nonlinear models [11] [13] are employed,
which consider both the nonlinear and inertial properties of the system being studied. The EMS is
simulated by determining multidimensional transient characteristics (MTCs) based on "input-

enabling accurate recording of eye responses to visual stimuli. The construction of the model
involves an approximation simulation method using Volterra integro-power series (IPS) [14], [15]
and a least squares method (LSM) [14] to create the model based on Volterra integro-power
polynomials (IPP). The simulation methods for nonlinear dynamic systems (NDS) based on IPS and
IPP vary in their computational approaches, providing distinct methodologies for NDS
simulation [15].
   The aim of this research is to evaluate the accuracy of EMS simulation by examining the errors
in calculating MTCs when using nonlinear dynamic simulation methods based on IPS and IPP
models. This study also addresses the development of algorithmic and software tools for extracting
EMS dynamic characteristics from eye-tracking data and assessing the accuracy of different
simulation methods.

3. Theoretical background
In this study, the approximation method [15] and compensation method [14] are employed to
develop models using Integro-Power Series, while the least squares method (LSM) [13] is utilized
for constructing models based on Integro-Power Polynomials.
     Approximation Identification Method. The approximation identification method for NDS (method
of linear combinations of responses) in the time domain is grounded in isolating the n-th partial
component (PC) of the NDS response by constructing linear combinations of responses to test
signals with different amplitudes. This approach is an adaptation of methods originally based on
the Volterra series. It is proved in [14] that:
     Assertion 1. Let test signals a1 (t), a2x(t) aNx(t) be sequentially applied to the input of the NDS,
where N is the degree; a1, a2 aN are different real numbers, non-zero, satisfying the condition
|aj| 1 for j=1,2,...,N; x(t) is an arbitrary function. Then, the linear combination of the system's
responses to these inputs equals the n-th PC of the response to the input signal x(t) with an
accuracy up to the discarded terms of the Integro-Power Series of order N+1 and higher:
                                N

                                c y[a x(t )] = y [ x(t )] + ,
                                j =1
                                       j   j         n                                               (1)

where
   y n [ x(t )] = y n (t );


                                                                                                     456
                                     t      t                                     n
    y[a j x(t )] =  a nj    wn (t − τ1 ,..., t − τ n ) x( τ i )dτ i ;
                        n =1          0      0                                    i =1
           N        
    =  с j  y n [ x(t )],
          j =1    n = N +1

if j are real coefficients such that
                                                                      AN c = b ,                                                          (2)
where
                                                   a1          a2                    aN          c1       b1 
                                                   2                                    2                   
                                                    a           a 22                  aN            c          b
                                             AN =  1                                         , c =  2 , b =  2 ,
                                                                                                       
                                                   N                                                        
                                                  a1          a 2N                    N
                                                                                       a N        c N    b N 
here bl = 1 when l = n; bl = 0 when l n, l            N}.
    The system (2) always has a solution, and it is unique since its determinant differs from the
Vandermonde determinant only by the factor a1, a2 aN. Thus, for any real numbers aj, which are
non-zero and pairwise distinct, it is possible to find numbers cj such that the linear combination (1)
of the NDS responses equals the n-th term of the Integro-Power Series with an accuracy up to the
discarded terms of the series. By satisfying the conditions for forming the system of linear
algebraic equations (2), we obtain the relation (1).
    When test signals in the form of step functions (Heaviside functions        t)) with amplitudes a1,
a2 aN are applied to the input of the system being identified, we obtain estimates of the diagonal
cross-sections of the NDS multidimensional transient characteristics:
                                      N
   hˆn (t ,..., t ) = yˆ n (t ) =  c (jn) y(a j θ(t )) = c1( n) y(t | a1 ) + c2( n) y(t | a2 ) + ... + c N( n) y(t | a N ), n = 1, N ,   (3)
                                      j =1

where y(t|aj) = y(aj t)) are the NDS responses to the test signal with amplitude aj .
   Identification of NDS using the Least Squares Method. The method of NDS identification based on
the Volterra polynomial model in the time domain relies on approximating the NDS response y(t)
to an arbitrary deterministic signal x(t) in the form of an integro-power polynomial of N-th order
(N the order of the approximation model):
                                             N             N t                t                             n
                             y N (t ) =  yˆ n (t ) =  ...  wn (t − τ1 ,..., t − τ n ) x( τ i )dτ i .                                 (4)
                                                                    n times
                                             n =1          n =1 0             0                         i =1

   Valid assertion [14].
   Assertion 2. Let test signals a1x(t), a2x(t) aLx(t) be sequentially applied to the input of the
NDS; a1, a2 aL are different real numbers satisfying the condition 0<aj1 for j=1,2,...,L; x(t) is
an arbitrary deterministic signal, then
                               N                    N           t             t                         n                 N
        y N (a j x(t )) =  yˆ n (a j x(t )) =  a nj  ...  wn (t − τ 1 ,..., t − τ n ) x(τ i )d τ i =  a nj yˆ n (t )
        ~
                               n =1                 n =1        0
                                                                    n times
                                                                              0                        i =1               n =1            (5)
     for j, j = 1, L.
   The partial components in the approximation model yˆ n (t ) are found using the least squares
method. This allows obtaining estimates for them, where the sum of squares of deviations of the
NDS responses being identified, y[a j x(t )] from the model responses ~y N [a j x(t )] , is minimal, thus
ensuring the minimum mean square criterion
                                                                                                                  2
                                (                                
                                                                        )                       
                        L                                    L                  N
                 J N =  y (a j x(t )) − ~y N (a j x(t )) =   y (t | a j ) −  a nj yˆ n (t )  → min .
                                                         2
                                                                                                                                          (6)
                       j =1                                 j =1              n =1             
                                                                                                                                          457
   Minimizing criterion (6) boils down to solving a system of normal equations Gauss, which in
vector-matrix form can be expressed as:

                                                                AAŷ = A y ,                                                             (7)
where
                              a1 a12  a1N            y (t | a1 )          yˆ1 (t ) 
                                    2       N         y (t | a )          ˆ         
                               a a2  a2                                        y (t )
                         A= 2                   , y=            2 
                                                                       , ŷ =  2  .
                                                                       
                                                                                    
                                                                               yˆ N (t )
                                     2       N
                             a L a L  a L          y (t | a L )
   If the identified system is supplied with test signals in the form of step functions with
amplitudes a1, a2 aL, we obtain estimates of the transient characteristics hˆ1( N ) (t ) and the diagonal
cross-sections of the transient characteristics of the eye movement system hˆ ( N ) (t , t ) , hˆ ( N ) (t , t , t )               2   3

hˆN( N ) (t ,..., t ) [14].
     Responses of the investigated EMS models are generally calculated based on expressions:
                              ~
                              y j (t | a j ) = a j yˆ 1 (t ) + a 2j yˆ 2 (t ) + ... + a Nj yˆ N (t ), j = 1, L ,                           (8)
or
                  ~
                  y (t | a j ) = a j hˆ1( N ) (t ) + a 2j hˆ2( N ) (t , t ) + ... + a Nj hˆ N( N ) (t ,..., t ), j = 1, L .                (9)
    Compensation Identification Method. The formalism of the method for determining the
intersections of n-th order transient characteristics of nonlinear dynamic systems is based on the
following assertion [14].
    Assertion 3. Let the test inputs be the sum of n step signals x k (t ) = a k θ(t − τ k ) (k=1,2,...,n), shifted
in time by 1, ..., n. Then, for an NDS with a single input and a single output, the estimate of the
intersection of the n-th order transient characteristic is
                                                                                                   n
                                                                         −1
                                                                                                n+  δk
                                                        n                         1
                        hˆn (t − τ 1 ,..., t − τ n ) =  n!  a k                  (−1)         k =1
                                                                                                          y(t | δ1 ,..., δ n ) .           (10)
                                                        k =1                δ1 ,...,δ n = 0

where y(t | δ1 ,..., δ n ) is the NDS response at time t when subjected to a multi-step input signal with
amplitudes ak, obtained as a result of processing experimental data based on (11). If δ k = 1 , the test
input contains a step signal shifted by k; otherwise, if δ k = 0 , it does not contain it.
   In certain cases, we have:
for n=1
                                                         y(t |a1 )
                                              hˆ1 (t ) =           ;                                                                       (11)
                                                           a1
for n=2
                                                 1
                                hˆ2 (t , t ) = 2 y(t |a2 ) − 2 y(t | a1 ) ;                                                              (12)
                                               2a1
or
                                                            1
                                         hˆ2 (t , t ) =         y(t |a3 ) − y(t | a1) − y(t |a2 ) ;                                      (13)
                                                          2a1a2
for n=3
                                                   1
                               hˆ3 (t , t , t ) =       y(t |a3 ) − 3 y(t | a2 ) + 3 y(t | a1 ) .                                        (14)
                                                  6a13
where a1, a2 = 2a1, a3 = 3a1 are the amplitudes of test signals.
                                                                                                                                           458
   The responses of the second and third-order models are calculated accordingly using the
expressions:
                                       ~
                                       y (t | a j ) = a j hˆ1 (t ) + a 2j hˆ2 (t , t ) .                     (15)
                       ~
                       y (t | a j ) = a j hˆ1 (t ) + a2j hˆ2 (t , t ) + a3j hˆ3 (t , t , t ) , j = 1, L .   (16)
4. Research results
The EMS's responses to the test step signals, defined as x(t) = aj t) with amplitudes aj (j=1, 2, 3):
a1=1/3, a2=2/3, a3=1 were analyzed. These responses formed the basis in the construction of
Volterra models [13]. Horizontal visual stimuli displayed at varying distances from the starting
position on a monitor were utilized as test signals, effectively simulating the application of step
signals with different amplitudes to the EMS. The responses of the EMS were recorded using eye-
tracking technology, integrating both hardware and software components. In the simulation
process, when applying the approximation method, models based on IPS are identified as
M1.N/x: (N order of approximation, x number of test signals; a1 aL amplitudes of
the test step signals). For models based on IPP, the least squares method (LSM) is used, resulting in
models designated as M2.N/x:. Additionally, when employing the compensation method
of simulation, models are determined as M3.N/x:.
    For EMS simulation in work [15], experimental "input-output" data were gathered using three
test step signals with amplitudes a1=1/3, a2=2/3 and a3=1. The Tobii Pro TX300 eye tracker was
employed to collect these experimental data (Fig. 1), from which transient characteristics were
determined for models M1.N, M2.N, and M3.N for N=1 (linear model), N=2 (quadratic model), and
N=3 (cubic model). The transient processes of EMS responses to visual stimuli with varying
amplitudes are depicted in Fig. 2.
    The software tools were developed using the Python programming environment.




Figure 1:EMS responses to visual stimuli of Figure 2:Transient processes of EMS responses
different amplitudes                        to visual stimuli of different amplitudes

    To evaluate the accuracy of the developed models for varying amplitudes of the test signals a1,
a2 and a3, the metric applied is the normalized root mean square error (NRMSE):
                                                                                 1/ 2
                                        M
                                          (                 ~               )
                                                                          2 
                                         y (t m | a j ) − y (t m | a j ) 
                               ε a =  m =0                                  , j=1, 2, 3;           (17)
                              j
                                                M
                                                                            
                                                y (t m | a j )
                                                      2
                                                                            
                                               m =0                        
where y(t m | a j ) and ~y (t m | a j ) are the responses of the EMS and the model of the EMS to the test
signal in the form of a step function with amplitude aj, measured/ computed at the time instant tm
(tm is the observation time of the EMS responses); j=1,2,3.
                                                                                                     459
   For N=1, transient characteristics hˆ1 (t | a j ) (j=1, 2, 3) were obtained based on the responses
y (t | a1 ) or y (t | a 2 ) or y (t | a3 ) , as shown in Fig. 3. Fig. 4 shows transient characteristics graphs of
the M2.1/2 models, which were calculated based on two responses: y (t | a1 ) and y (t | a 2 ) , or
 y (t | a1 ) and y (t | a 3 ) , or y(t | a2 ) and y(t | a3 ) and the transient characteristic graph of the model M2.1/3.
For the M1.1 identification method does not allow to calculate the transient characteristics based
on two or three responses [14].




Figure 3: Transient characteristics of the EMS Figure 4: Transient characteristics of the EMS
models M1.1/1 and M2.1/1, built using test models M2.1/2, built using test signals with
signals with amplitudes 1; 2; 3                amplitudes: 1 and 2; 1 and 3; 2 and 3; and
                                               M2.1/3

   Table 1 shows the values of the percentage NRMSE estimates of the responses built using the
identification methods of the EMS models M1.1/1 and M2.1/1; in Table 2 of the models M2.1/2
and M2.1/3.

Table 1
Percentage Normalized Root Mean Square Error of the EMS models M1.1/1 and M2.1/1, %
                                       Amplitudes of test signals                    Mean            Maximum
          Models
                                    a1             a2              a3                value            value
         M1.1/1:a1                  0             18.2            20.3                19.2             20.3
         M1.1/1:a2                 17.2            0               5.5                11.4             17.2
         M1.1/1:a3                 18.5           5.3               0                 11.9             18.5

     For N=2, based on the two responses y(t|a1) and y(t|a2), or y(t|a1) and y(t|a3), or y(t|a2) and y(t|a3),
the corresponding multidimensional transient characteristics                      1(t|aj,ak) and 2(t,t|aj,ak),
j, k = 1,2,3; j k were obtained. In the models M1.2/2 and M2.2/2, identical MTCs
  1(t|aj,ak) and 2(t,t|aj,ak) were obtained for the same experimental data.


Table 2
Percentage Normalized Root Mean Square Error of the EMS models M2.1/2 and M2.1/3, %
                                      Amplitudes of test signals                    Mean             Maximum
         Models
                                   a1             a2             a3                 value             value
      M2.1/2:a1, a2               13.8           3.6             7.1                 8.2               13.8
      M2.1/2:a1, a3               16.7           4.9              2                  7.9               16.7
      M2.1/2:a2, a3                18            3.7             1.7                 7.8                18
     M2.1/3:a1, a2, a3            16.7           3.5             2.5                 7.6               16.7
                                                                                                                   460
    Fig. 5 displays the primary-order transient characteristics graphs, while Fig. 6 shows the
diagonal cross-sections of the second-order MTCs for the EMS models M1.2/2 and M2.2/2. These
intersections were computed based on the responses: y(t|a1) and y(t|a2), or y(t|a1) and y(t|a3), or
y(t|a2) and y(t|a3).




Figure 5:Transient characteristics of the first Figure 6:Diagonal cross-sections of the second-
order of the EMS, models M1.2/2 and M2.2/2      order transient characteristics of the EMS, models
                                                M1.2/2 and M2.2/2

   The transient characteristics of the EMS model M3.2/2 are illustrated in Fig. 7. The
corresponding responses for the same model are given in Fig.8.
   The graphs of the multidimensional transient characteristics 1(t) and 2(t,t), which were
determined based on the three responses y(t|a1), y(t|a2), y(t|a3) for the EMS model M2.2/3, are shown
in Fig. 9. The corresponding graphs comparing the M2.2/3 responses with the EMS responses to
identical test signals are presented in Fig. 10. Analogous results were obtained for the M3.2/3 and
are shown in Fig. 11 and Fig. 12.
   For N=3, the graph of the primary-order transient characteristic, along with the graphs of the
diagonal cross-sections of the second and third-order multidimensional transient characteristics for
the EMS model M3.3, are displayed in Fig.13. The responses of the EMS model M3.3 are shown in
Fig.14. Analogous results for the EMS models M1.3 and M2.3 are shown in Fig. 15 and Fig.16.




Figure 7:Multidimensional transient                Figure 8:Responses of the EMS and the model
characteristics of the EMS, model M3.2/2           M3.2/2




                                                                                                 461
Figure 9:Multidimensional transient        Figure 10:Responses of the EMS and the model
characteristics of the EMS, model M2.2/3   M2.2/3




Figure 11:Multidimensional transient       Figure 12:Responses of the EMS and the model
characteristics of the EMS, model M3.2/3   M3.2/3




Figure 13:Multidimensional transient       Figure 14:Responses of the EMS and the model
characteristics of the EMS model M3.3      M3.3




Figure 15:Multidimensional transient       Figure 16:Responses of the EMS and models

                                                                                   462
characteristics of the EMS, models M1.3 M1.3 and M2.3
and M2.3

   Table 3 provides the values of the percentage normalized root mean square error of the
response estimates for the constructed EMS models M1.2/2 and M2.2/2, and the model M2.2/3.
Table 4 presents the percentage NRMSE values for the models M3.2/2 and M3.2/3, and the model
M3.3/3.

Table 3
Percentage Normalized Root Mean Square Error of the EMS models M1.2/2, M2.2/2, M2.2/3, %
                                Amplitudes of test signals            Mean         Maximum
        Models
                             a1             a2             a3         value         value
  M2.2/2:a1, a2              0              0              19           19            19
  M2.2/2:a1, a3              0             9.1             0           9.1           9.1
  M2.2/2:a2,a3              17.3            0              0           17.3          17.3
  M2.2/3:a1, a2, a3         8.2            4.3             1           4.5           8.2

Table 4
Percentage Normalized Root Mean Square Error of the EMS models M3.2/2, M3.2/3, M3.3, %
                                 Amplitudes of test signals           Mean         Maximum
         Models
                              a1             a2              a3       value         value
  M3.2/2: a1, a2             17.2           18.2            37.4       24.3          37.4
  M3.2/3:a1,a2,a3            6.1            11.3            07.8        8.4          11.3
  M3.3/3: a1,a2,a3           9.3            19.9            48.6       25.9          48.6

    For N=3, based on the three responses y(t|a1), y(t|a2), y(t|a3) the transient characteristics
 1(t), 2(t,t), 3(t,t,t) were obtained. For the models M1.3/3 and M2.3/3, identical MTCs were
obtained, and the responses of the models practically coincide with the responses of the EMS for
the same input signals. At the same time, Figs. 13 and 15 demonstrate that the transient
characteristics of the third-order models exhibit instability.
    On Fig. 17, a comparative analysis diagram of errors based on the percentage NRMSE criterion
constructed using identification software tools for EMS models: M2.1/1, M2.1/2, M2.1/3, is
presented. On Fig. 18, the same analysis is provided for models M2.2/2, M3.2/2, M2.2/3, M3.2/3 and
M3.3 (based on average values). The EMS models M1.3 and M2.3 are not shown in the diagram
because they have negligible deviations from the EMS responses.

5. Conclusion
Python-implemented tools for nonlinear dynamic identification were applied to develop
mathematical models of the human eye movement system (EMS) based on Volterra integro-power
series (IPS) and integro-power polynomials (IPP). This study aimed to evaluate the accuracy of the
developed models in the form of first-order transient characteristics and multidimensional
transient characteristics of the second and third orders, derived from eye-tracking data under the
influence of visual test stimuli (step signals of varying amplitudes). To construct models based on
empirical data, computational identification methods were employed, including the compensation
method, approximation method, and least squares method (LSM). The identification errors of EMS
models were assessed using normalized root mean square error (NRMSE), reflecting deviations
between model responses and EMS responses to identical test signals.


                                                                                               463
Figure 17:Comparative analysis of the average Figure 18:Comparative analysis of the average
percentage error values of the EMS models: percentage error values of the EMS models:
M2.1/1, M2.1/2, M2.1/3                        M2.2/2, M3.2/2, M2.2/3, M3.2/3, M3.3

   Models obtained using the approximation and LSM methods on the same test signals exhibited
identical transient characteristics, as these characteristics coincided within the convergence region
of the Volterra series. It was found that the accuracy of models, represented by transient
characteristics, improved with an increasing number of test signals. For the linear model, the
average error decreased from 11.4% with one test signal to 7.6% with three test signals when LSM
was applied. For the quadratic model, an average error of 9.1% was observed with two test signals,
which decreased to 4.5% when three signals were used. For the cubic model, using three test signals
resulted in nearly identical responses between the model and the EMS. However, it was established
that third-order models exhibited instability in their transient characteristics, limiting their
practical applicability.
   The compensation method required fewer computational resources but produced models with
significant errors, rendering them unsuitable for diagnostic studies. The best result for the
quadratic model, achieved using three test signals, yielded an average error of 8.4%.
   This study, for the first time, analyzed the errors of EMS mathematical models in the form of
IPS and IPP derived from eye-tracking data using three test step signals of varying amplitudes. The
analysis employed compensation, approximation identification methods, and the LSM. It was
determined that the most accurate EMS model is the quadratic Volterra polynomial, identified
based on three EMS responses. Therefore, integral quadratic models are recommended for
diagnostic studies of the human psychophysiological state.

Declaration on Generative AI
The authors have not employed any Generative AI tools.

References
[1] M. Khamis, Y. Sugano, L. Sidenmark, Proceedings of the 2024 Symposium on Eye Tracking
    Research and Applications, ETRA 2024, Glasgow, United Kingdom, June 4-7, 2024. ACM 2024,
    pp. 1-525.
[2] J. Opwonya, D. N. T. Doan, S. G. Kim et al., Saccadic Eye Movement in Mild Cognitive
                                                                      -Analysis, Neuropsychol
    Rev 32 (2022), pp. 193 227. doi:10.1007/s11065-021-09495-3.
[3] D. Jansson, O. Rosén, A. Medvedev, Parametric and nonparametric analysis of eye-tracking
    data by anomaly detection, IEEE Transaction control system technology 23 (2015), pp. 1578
    1586. doi:10.1109/TCST.2014.2364958.


                                                                                                 464
[4] V. Bro, A. Medvedev, Continuous and Discrete Volterra-Laguerre Models with Delay for
     Modeling of Smooth Pursuit Eye Movements, IEEE Transactions on Biomedical
     Engineering 70 (2023), pp. 97 104. doi:10.1109/TBME.2022.3185669
[5] L. Lanata, L. Sebastian, F. Di Gruttola, S. Di Modica, E. P. Scilingo, A. Greco, Nonlinear
     Analysis of Eye-Tracking Information for Motor Imagery Assessments, Frontiers in
     Neuroscience 13 (2020) 1431. doi:10.3389/fnins.2019.01431.
[6] B. Keehn, P. Monahan, B. Enneking et al., Eye-Tracking Biomarkers and Autism Diagnosis in
     Primary         Care,       JAMA           Netw         Open         7(5):e2411190       (2024).
     doi:10.1001/jamanetworkopen.2024.11190.
[7] K. Weiss, M. Kolbe, Q. Lohmeyer, M. Meboldt, Measuring teamwork for training in healthcare
     using eye tracking and pose estimation, Front. Psychol. 14:1169940 (2023), pp. 1 12.
     doi:10.3389/fpsyg.2023.1169940.
[8] W. Sun, Y. Wang, B. Hu, Q. Wang, Exploration of Eye Fatigue Detection Features and
     Algorithm Based on Eye-Tracking Signal, Electronics 13, No.10:1798 (2024), pp. 1 19.
     doi:10.3390/electronics13101798.
[9] H. Griffith, D. Lohr, E. Abdulin, O. Komogortsev, GazeBase, a large-scale, multi-stimulus,
     longitudinal eye movement dataset, Scientific Data, Nature, 8(13) (2021), pp. 1 9.
     doi:10.1038/s41597-021-00959-y.
[10] J. Yin, J. Sun, J. Li, K. Liu, An Effective Gaze-Based Authentication Method with the
     Spatiotemporal Feature of Eye Movement, Sensors, 22(8) (2022), pp. 1 18.
     doi:10.3390/s22083002.
[11] S. Solodusha, Y. Kokonova, O. Dudareva, Integral Models in the Form of Volterra Polynomials
     and Continued Fractions in the Problem of Identifying Input Signals, Mathematics, 11(23),
     4724 (2023). doi:10.3390/math11234724.
[12] F. J. Doyle, R. K. Pearson, B. A. Ogunnaike, I, Identification and Control Using Volterra
     Models, Communications and Control Engineering. Springer, London, 2001, 314 p.
[13] V. Pavlenko, M. Milosz, M. Dzienkowski, Identification of the oculo-motor system based on
     the Volterra model using eye tracking technology, in: Proceedings of the 4th Int. Conf. on
                                                                                          Journal of
     Physics: Conference Series, Vol. 1603 (2020), pp. 1 8. doi:10.1088/1742-6596/1603/1/012011.
[14] V. Pavlenko, S. Pavlenko, Deterministic identification methods for nonlinear dynamical

     (2018), pp. 9 29. doi://10.15276/aait.01.2018.1.
[15] V. Pavlenko, T. Shamanina, V. Chori, Nonlinear dynamic model of the oculo-motor system
     human based on the Volterra series, in: J. Awrejcewicz (Ed.), Perspectives in Dynamical

                                                                         442. doi:10.1007/978-3-031-
    56496-3_27.




                                                                                                 465