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