Heart Rate Variability Control Using a Biofeedback and Wearable System

Heart rate variability is an important physiological parameter in medicine. This parameter is used as an indicator of physiological and psychological well-being and even of certain pathologies. Research on biofeedback integrates the fields of biological application (physiological behavior), system modeling, and automated control. This study proposes a new method for modeling and controlling heart rate variability as heart rate acceleration, a model expressed in the frequency domain. The model is obtained from excitation and response signals from heart rate variability, which through the instrumental variables method and the minimization of a cost function delivers a transfer function that represents the physiological phenomenon. This study also proposes the design of an adaptive controller using the reference model. The controller controls heart rate variability based on the light actuators designed here, generating a conditioned reflex that allows individuals to self-regulate their state through biofeedback, synchronizing this action to homeostasis. Modeling is conducted in a target population of middle-aged men who work as firefighters and forest firefighters. This study validates the proposed model, as well as the design of the controllers and actuators, through a simple experiment based on indoor cycling. This experiment has different segments, namely leaving inertia, non-controlled segment, and actively controlled segment.


Introduction
Heart rate variability (HRV) is a parameter widely used as an identifier of physiological and mental well-being in people [1]. This cardiac parameter is used in systems in the area of physical well-being, as well as in the monitoring of vital signs. The monitoring of physiological variables has been expanded with the use of mobile devices, allowing indices of physical well-being, stress, etc. to be established [2].
The use of HRV is widespread and is used in areas of medicine for the detection and recognition of stress [3,4]; cardiovascular health [5]; insomnia [6], and autoimmune diseases, such as ulcerative colitis [7].
This parameter is obtained from electrocardiograms (ECGs) measuring the duration of the R-R intervals [8]. Having the heart rate, the HRV measurement needs a time window (both in the time domain and in the frequency domain) for its estimation; this presents an inherent delay in its processing time [9]. This is a problem, as the responses in the feedback system must be tolerant of these [10] delays.
On the other hand, HRV control represented by an adequate response time implies a direct benefit in the reestablishment of the balance of the autonomic nervous system, which is associated with involuntary actions of the body, including pain and vasovagal syncope [11][12][13][14]. Currently, the search for self-control of physiological functions in the human body has led to the application of feedback on these variables in the individual. This is shown in the work that seeks to train the responses of individuals to changes in HRV [15]. This concept has also evolved by integrating automatic control [2,16].
However, the raised solutions do not yet solve a fundamental problem, the response time of the HRV. The principal efforts associated with time series processing have optimized the response to obtain this parameter [17,18]. Solving the self-control of HRV using biofeedback over continuous time [2] contributes to physical and psychological well-being and safety in individuals such as firefighters and forest firefighters, who are constantly subjected to environments with extreme and aggressive conditions. Having a consequence of fatigue, dizziness, weakness, fainting, and so on [11,19,20].

Related Work
Controlling HRV has an important meaning in the safety and well-being of individuals subjected to work under hostile conditions. Currently, it is a challenge to find a model of HRV as a system, either dynamically in a discrete or continuous way, and it requires control strategies and actuators that allow, through their actions, to train individuals so that they respond to an external stimulus achieving self-control HRV.
The HRV is used in several areas of knowledge as an indicator of metabolic function and stress levels, as shown in [5,21]. This parameter was shown to be correlated with insomnia using polysomnographic studies that correlated changes in HRV with this pathology [6]. Similarly, this parameter is correlated with stress and mental health, as shown in [7,[22][23][24].
From a cardiovascular point of view, HRV is currently used to correlate its behavior with severe and chronic heart disease [4] In the area of health wellness, the use of HRV is associated with physical performance metrics, which it is used for forms of physical training [25], and even behavioral training is associated with the controllability of this parameter [26].
In addition, the HRV has been used to determine the cardiac status of people affected by the COVID-19 virus [27]. The latter shows and supports HRV as a variable to control, has great challenges, and new potential areas of knowledge based on its control.
In feedback systems, the dynamic representation of the system under study is essential to achieve this goal. However, the differential representation of a physiological system is a really complex task, as models will always be subject to nonlinearities [28,29]. In this context, the identification of systems through the observation of an experiment to correlate variables that are not directly related but that are integrated by the principle of causality makes the use of instrumental variables (IV) an interesting method to apply in the estimation of a dynamic model of the HRV based on experimental observations, as a mechanical system, for example [30][31][32][33].
Considering controlling HRV can be considered as an input to a biofeedback system. Therefore, it makes it possible to perform actions that control the behavior of this variable [26]. An example of this is shown in conventional medicine in patients with joint and muscle pain problems [34][35][36], where biofeedback techniques generate self-correction of body positions that benefit study subjects.
This physiological and psychological control effect has been applied in the area of video games [37], also in physical training, training people to improve breathing, and even self-control of HRV in physical activities [15,[38][39][40]. In the area of medicine, applications are being studied in anorectal recovery therapies, integrating the IoMT (Internet of Medical Things) concept [41,42].
Moreover, the use of biofeedback has proven results in the medical area with applications in the control of HRV showing that through external excitations such as harmonic sounds it is possible to improve HRV, and thus the homeostasis or self-regulatory process of the organism [43], in this solution HRV and biofeedback are integrated to generate self-control of mental stress in drivers [3].
Referring to automatic control, some controller applications based on fuzzy logic have been applied in biofeedback, such as [16,[44][45][46]. In addition, other types of applications associated with automatic control are associated with the control of sleep apnea and physical behavior in sports, among other related areas, controllers have been applied in continuous and discrete time [47][48][49].
Regarding control signals in a biofeedback system, light signals are signals that practically have no time delay compared to the HRV response. Considering the psychological effect of light on individuals that works associated with the study of the effects of light intensity and its flickering frequency has a direct effect on the states of alertness of individuals [50,51]. The same goes for effects on individual stress as evidenced in [52,53], and even in the circadian cycle [54].
Note that in this work different areas are involved regarding the control of physiological variables; therefore, we resume the main contributions as follows.
• We propose a model identification method via a transfer function for HRV focused on a specific population of individuals using the method of instrumental variables. • We propose the design of an MRAC controller (Model Reference Adaptive Controller) using the estimated HRV model. • We describe the design of light actuators to establish a conditioned reflex. • We propose an experiment to validate the HRV model obtained through biofeedback.
The objective of our study was to investigate the ability to control the viability of heart rhythm through biofeedback, generating conditioned reflexes based on the training of individuals responding to external stimuli, such as light signals. Our hypothesis is that the use of an HRV-focused biofeedback system allows HRV regulation.

Materials and Methods
This research focused on the search for a dynamic and controllable model based on biofeedback. This is to generate a controlled system that allows individuals to self-monitor their HRV.
This research used, as input data, segments of electrocardiograms of a specific population: firefighters and forest firefighters who are exposed to aggressive and extreme environmental conditions that make the process of homeostasis very dynamic and reflected in the behavior of HRV [19,55]. The characteristics of this studied population are given by a universe of 12 individuals who have the following characteristics: This research was approved by the Institutional Ethics Committee of the University of Santiago de Chile. Each participant in this research provided their written informed consent. This specified the tests that were carried out and also specified that the data obtained during the tests were used in this investigation. On the other hand, the document also specified that participants could withdraw from the research for no specific reason.
From these ECG segments, the representative model of the described population is obtained using instrumental variables for this. With the design of controllers and actuators, the biofeedback system shown in the block diagram in Figure 1 was modeled, designed, and implemented.

System Modeling
The modeling method of the dynamic system is based on the relationship between excitation (input) and response to the same (output). After confirming that these signals can be correlated, a representative model of HRV dynamic behavior in its transfer function form using the instrumental variables (IV) method.
To obtain and develop the input and output signals, 12 ECG samples are collected from individuals with the characteristics mentioned above. ECG segments have a period equal to T = 60 s with a sampling frequency of f = 360 Hz and a normal sinus rhythm. These segments were processed to obtain the HRV behavior of each individual.
HRV can be in relation to time defined as Equation (1) [17].
If we assume from Equation (1) and Figure 2 that the R1 peak travels a ∆x distance in ∆t time, this implies that the peak has a v speed that in an ECG corresponds to HR [56]. According to this analogy, the ratio of the speed change in the peak Ri can be obtained when this "travels" an x distance, that is, the acceleration of the heart rate shown in Equation (3). Then, in the form of HR acceleration, an HRV expression is obtained almost immediately, which does not depend on the time windows (T-period of analysis) but only on changes in two points or velocities (HR), which makes it different from the temporary HRV (Equation (1)).
This action is carried out by the data acquisition block, specifically by the sensor, as seen in Figure 1.
The ECG collected is processed using the algorithm described in [56], obtaining the acceleration HR of each ECG.

Input Signal
The input variable x(t) is a step function based on the HRV of the samples described above. The objective of this signal is to establish the level of HR acceleration, which is representative of HRV from the ECG samples of the population under study. The sample mean µ is taken as a sampling estimator of the data described, which requires a known probability distribution similar to its origin in all samples as its origin to be valid. This condition is tested through the Kolmogorov-Smirnov test [57,58]. Once confirmed that all samples originate from the same type of probability distribution, the step of λ mean is generated from these experimental data, with Equation (4) as the hypothesis. Hypothesis 1.F originates from an exponential distribution with parameters EXP(x, λ).

Hypothesis 2.F Does not originate from an exponential distribution with parameters EXP(x, λ).
Taking this into account, the contrasting probability distribution is F 0 .
Here, the empirical and theoretical distributions X i ∈F and x i ∈ F 0 are compared through Kolmogorov-Smirnov's distance, obtaining the deviation of two empirical cumulative distributions, or a function of a reference theoretical cumulative distribution. This parameter was used because of its sensitivity to differences in location, as in the cumulative distribution function. Kolmogorov-Smirnov's distance D KS is defined by Equation (5).
The function to counteract is composed of the vectorF = a t = [a 1 , a 2 , . . . , a12]. This vector contains all ECG data from the population studied, as shown in Figure 3. In the case of the theoretical or reference cumulative distribution F 0 , it is estimated to be an exponential distribution in the form of x ∼ EXP(λ) for λ = 0.0323 mm/s 2 . Figure 4 represents the superposed accumulative distributionsF for the experimental or empirical distribution and F 0 for the theoretical or reference distribution. To accept the hypothesis H 0 , the parameter D KS is set, which is the largest absolute difference between both distributions (theoretical and experimental) (D + KS , D − KS ). The condition D KS ≤ D α should be met according to Table 1. It is concluded that all the samples of experimentF = a t originate from an exponential distribution of parameters x ∼ EXP(λ) and D = max{D + KS , D − KS } and a confidence interval α = 0.05 [57,58].
Having proven that the sample distributions are equal and known, it can be stated that the mean is representative and, therefore, the value λ = 0.0323 mm/s 2 can be used as the maximum value of the step as an input signal (Equation (6)).
For a better understanding and reference of the results, millimeters (mm) are used as the measurement unit for length, since the normalized distance for the measurement and interpretation of ECG is the normal distance traveled by an R peak during 1 mm with respect to 0.04 ms [59].

Output Signal
The λ of amplitude was obtained, and considering the duration of the different changes in HR acceleration from the collection of ECGs, an average threshold of the transient period is observed (tending to zero to maintain homeostasis [60]). Based on the average time of the ECG samples, in which the permanent regime period tends to zero, a response representative of the HR acceleration change process is established. This is shown in Figure 5, in which the HR acceleration data collected from the ECG superposed. The twelve samples that form the HR acceleration vector a t , are shown in Table 2, in which the average duration of the transient period before an HRV corresponds to ∆t = 4.5 ms since the acceleration a 2 is the most representative sample of the experiment and transient behavior of all samples, the signal y(t) was used as a response to the final model.

Identification Model
Model identification is conducted to convert the model into a transference function. To this end, the instrumental variables method is used as an algorithm to determine the parameters of the sought model, employing the MSE (Mean Squared Error) index as a cost function to minimize error or differences between the responses of the sought model (y(t)) and the known response (y m (t)). This is based on the adjustment IVs conduct in each parameter interaction (poles and zeros) of the sought transfer function, thereby identifying the dynamic model for HR acceleration in the studied population.
The input x(t) and output y(t) variables of the system are known and represent the excitation and response of the theoretical model; however, there is a causality relationship between this input signal x(t) and the output signal y(t). This allows the correlation between x(t) and y(t) based on the transfer function. The IV method uses the regression model given by Equation (7) [31,33,61]. Here,ŷ(k) is a known parameter of the system and the output of the φ(k) system is a vector of m × n dimensions that have known parameters and are related to the input of the system, θ is an m × n vector of unknown parameters (transfer function). In turn, k ∈ N is the discrete sampling time of the variables under study.
With data fromŷ(k) and ϕ(k) we attempt to minimize the cost function given by the least squares (Equation (8)), assessing y ( i) from the known signal, and y m (i) corresponding to the output of the model under construction.
The estimation of the unknown parameter is conducted based onθ, which is replaced by the vector Φ, a set of parameters correlated to the regression of transfer function, i.e., the instrumental variables [61]. This means Φ T = [ϕ(1), ϕ(2), . . . , ϕ(k)]. Thus, the function to estimate the transfer function is defined in Equation (9).
where W(k) = [φ(1),φ(2), . . . ,φ(k)] are the regression parameters given by Equation (10) The iterations executed to obtain the desired transfer function begin with a known transfer function form of the 1st and 2nd order, respectively. Shortening the transfer function parameter search based on the error yielded by (k) = y(t) − y m (k). This process is summarized in Figure 6. When applying the IV method in x(t) and y(t), transfer function H(s) is obtained, this function changes its parameters according to the values of a minimum MSE calculated based on the output of the desired system and the reference output, modifying the poles and zeros of the transfer function [31], which are contained inθ. As previously mentioned, to limit the search, the following two types of transfer functions were estimated.

•
Transfer function type 1: Transfer function type 2: H 2 (s) = s + a s s + s + a In this case, the results obtained for the transfer function are presented in Figure 7, in which the transfer function corresponds to Equation (11).
(11) Figure 7 shows the acceleration response signal of the known y m (t) and the response signal of the identified model y(t), These signals present a minimal error of MSE = 16.93% and a representativity of 83.03% with respect to the known signal y m (t).

Controller Design
To design the controller, the MRAC controller is used through an adaptation mechanism that modifies the adaptation parameters of the PID controller based on a reference model [10]. The architecture in Figure 8 was considered for the development of the controller. The controller corrects the deviation e m (t) associated with the previously obtained response of the model y(t) with respect to the output of the reference model y r (t). Based on the behavior of both outputs, the cost function j is minimized. In this way, the parameters of θ(t) are adjusted using the adaptation mechanism. The cost function is defined by the minimum integrative error of the output error expressed in Equation (12), in which γ i is the gain of the adaptation of the parameters θ i .
Parameter adaptation is based on an MIT (Massachusetts Institute of Technology 1958) rule [62]. This seeks to rapidly reduce the error towards the gradient of the controller θ(t), parameters, i.e., these parameters are adjusted based on the reduction of e m (t) by means of the minimization of the cost function j. This is shown in Equation (13).
The temporary variation in the controller parameters ∆θ i with respect to their last adjustment, including the adaptation gain γ i is expressed in Equation (14).
This gain was selected to reduce the settling time and, at the same time, decrease the transition of the adaptation process of the adaptive system. By including the calculation of the gradient and adaptation parameters into the error of the responses e m (t), Equation (15) is obtained.
As the output is a single variable, in which the gradient of the term y m (t) tends to zero, and as the reference model y r (t) does not depend on θ i (t) the calculation of the tracking error gradient is reduced to the partial derivative of the heart rate acceleration model kG(s) in the closed loop with respect to the parameter of the control law θ i (t), as shown in Equation (16).
Referencing these equations to the control signal u in Figure 8. The control law based on the adaptable parameters θ i and the control signal or set point u c is defined according to Equation (18).
From this, these parameters are shown in terms of error in Equation (19).
If kG(s) is the heart rate acceleration (HRV) model, k o G(s) is the reference model with its corresponding gain. Calculating the sensitivity derivative ∂y r ∂θ i , indicates how the error changes according to the adaptation parameters θ. Its negativity indicates the minimization of the adaptation time that defines the robustness of the system (Equation (19)).
The variation of the adaptation parameters over time can be defined using Equation (20).
Finally, the convergence point k = k 0 is obtained based on Equations (20) and (21).
Given the parameters implemented in the controller, it is necessary to know the homologous reference level of the HRV model in the form of a transfer function (Equation 11).
In the estimation of this model, the transfer function of the reference model H r (s) in Equation (23) is established using a settling time t s ≤ 1s of and a maximum overshoot of M p < 20%, and applying the criterion of 2% of t s [62].
The controller handles 5V signal levels. Therefore, the general system is normalized in all the blocks of the system under study. The PID algorithm was adjusted by the Ziegler & Nichols reactive curve method [10], In this way, the following parameters were estimated in the controller: P = 0.0138, I = 0.241 and D = 0.0012. Tuning this controller, the gain of the adaptation method can be empirically adjusted using γ T (i) = [0.1 3.3 6.6 9.9 13. 3 20]. Figure 9a shows the response of the output of the system y(t), as well as the responses to different γ T (i) values.
(a) (b) Figure 9. System output response and regulation. (a) Behavior of the adaptation mechanism γ i . (b) Response of the HRV system y(t) respect to the reference model y r (t). Figure 9a shows several gains of the adaptation mechanism γ i . One of these is chosen based on the energetic effort made by the controller when modulating the control signal in order to maintain the stability of the system. This is achieved with γ i = 20, obtaining a response that is controlled in terms of speed and overshoot of the system response. The above described is shown in Table 3 for the adaptation gains and their dynamic characteristics [62].  Figure 9b shows the response of the system y(t) with an overshoot superior to that of the reference model. However, this value is 16.77% and the settling time is 0.3604 s; therefore, it is an acceptable response according to the requirements established for the system response.

Actuator Design
The actuator of the system is based on light transduction; this means the conversion of the electric control signal into a light signal. In this case, the effect of light signals is considered within a spectrum range of visible light, as well as its effect on the behavior of individuals [63], being one of the first studies in the field. In this field, the works in [52,64,65] follow this line, demonstrating the psychological and physiological effect of light and its inherent properties, such as intensity and flickering frequency.
The actuator considers the work in [52,64] uses the effects of a short wavelength on the alertness and stress state of individuals.
When using a short wavelength λ o = [380-450] nm, combined with flickering effects in the range of f c = [0-10] Hz; blue light affects the attention and alertness of individuals, which can become a conditioned reflex that allows maintaining attention during the performed activity, as well as regulating a "calmness" state that helps body homeostasis (regulation of heart rate acceleration or HRV to keep speed constant).
The control signal works within a range of u c = [0-5] v, This tension is converted into the f c , frequency, which corresponds to the flickering frequency of the blue light in the actuator.
This voltage-frequency conversion is carried out via a VCO (Voltage Controlled Oscillator), which drives the continuous component of a modulated signal to vary in a linear and periodical way at the output. with a frequency that is directly related to the input of voltage v(t) that is modulated by the VCO. In this way, a response without time delay τ is obtained given the rapid response of the voltage-light intensity conversion to the responses of the dynamic system under study.
The VCO modulates the PWM signal response, which controls the overshoot current i γ of an LED diode of high brightness. After using Kirchhoff's voltage laws to analyze the circuit in Figure 10 an experimental model is obtained that relates the overshoot current i γ to the relative light intensity of the LED light i L as shown in Equation (24).
According to the circuit in Figure 10, the total voltage is v t = v d + v l , as a function of current. This relationship is expressed as v(t) = i t (R d (t) + R l ). The LED diode is not a linear load. Due to its low internal resistance, this parameter is not considered because an internal resistance R in (t) < 0.01R th is assumed to be 10% of the total circuit resistance. Therefore, the modulated current is the total current of the circuit, i.e., i t = i γ . However, the input voltage is modulated in PWM, which implies that the tension is not constant and that its variation is based on the modulation frequency v t = v( f (t)). In this way, a work cycle of %D = t 1 t 1 −t 0 is obtained. The modulated voltage makes the amplitude of its continuous component vary based on the PWM frequency. Direct voltage (mean component) is the input signal to generate the conduction current i γ according to Equation (25).
From Equation (25), and considering the signal amplitude of the constant PWM voltage and constant charge resistance R l constant, it is ensured that the direct current component is expressed according to Equation (26). This causes the desired effect on the brightness and frequency modulation of the LED diode, with %D = [0.001-0.0].
Since the action of the actuator cannot interrupt the vision of individuals, the actuator was attached to lenses or safety glasses with specially designed adapters that cover the "monocular vision" spectrum, which is between 94 • and 110 • from the center of the left and right eyes, respectively. This is achieved using the diffusion angle of the LED in the mentioned adapters, as shown in Figure 11.

Results
The experimental setup was designed to validate the proposed HRV acceleration model, as well as testing and monitoring the functioning of HRV control through the design of the proposed controllers and actuators, thus, confirming the hypotheses of this study.

Experiment Design
The experiment was conducted using indoor cycling. This activity has been selected because it has characteristics of intense cardiovascular effort, which produces a high HRV variance according to the intensity changes that occur during the experiment. A magnetic induction roller is used, which generates a controlled force in the opposite direction to the traction of the bicycle's back wheel. This implies a conversion of the torque generated by the biker's legs to counter the opposite force generated by the roller, producing high physical wear that leads to high cardiovascular activity and, therefore, a high HRV. The specifications for the force generated by the roller are as follows. In this test, the HR of the individuals was monitored based on the duration of the R-R interval of the cardiac cycle [56], bicycle speed, algorithm occurrence time, and HRV over time. The last variable is expressed as the acceleration of HR with respect to the above parameters. The MRAC controller generated control actions in the light actuators to maintain HRV acceleration. The sensors and actuators installed on the individual are shown in the testbed scheme in Figure 12.
The room temperature was 25°C and its relative humidity was 54%. Both variables were constant during the whole experiment.
The experiment was divided into four stages {A, B, C, D}: The complete test had a duration period of T = 630 s. This period was divided into the stages mentioned above, as shown in Table 4. The test was conducted with a different individual from those who participated in obtaining the model in Equation (11), but his physical characteristics belong to the universe of the obtained model, i.e., male, 42 years of age, the weight of 91 Kg and height of 173 cm.

Experimental Results
The results of the experiment are now presented, with a focus on the validation of the model and the study of the operation These results were polished using the LOESS (locally estimated smoothed scatterplot) method, which uses linear regression by least squares to adjust models over their local data subsets, creating a function that describes the event. In this case, a 2% scope of data was used for all variables.
When presenting data that are statistically significant, it should be proven that the data have the same distribution; in this case, a normal distribution is validated with p > 0.05 using the Kolmogorov-Smirnov method. Figure 13 graphically shows the results obtained divided by each experimental stage in terms of time and stage duration.
The variables shown are as follows • y exp Response as acceleration of heart rate mm/s 2 . • y bpm Response of heart rate as speed BPM. • v Response of speed of the bicycle in km/h.
Analyzing the graphs in Figure 13a-c, and using the {A, B, C, D} segments of the experiment as a reference, we have the following: in segment A, there is an increase in all variables due to the inertia output of the bicycle, with a maximum effort that reaches the peak of HR acceleration. A maximum HR of 129 BPM, an HR acceleration of 0.0078 mm/s 2 , and a speed of approximately 20 km/h are reached. In segment B, the individual performs free cycling based on his physical resistance, achieving a maximum HR of 140 BPM, an HR acceleration of 0.009 mm/s 2 , and a maximum speed of 51 km/h. In segment C, with the active controller, it is observed that due to the light modulation, HRV is controlled; therefore, HR has a linear increase but an acceleration tending to zero, controlling the roller speed and consequently the applied effort. In this case, a maximum HR of 101 BPM is reached, as well as an HR acceleration of 0.009 mm/s 2 and a speed peak in the last segment of 38 km/h. Finally, the D segment corresponds to free cycling.  Table 5 shows the descriptive statistics of the entire sample universe obtained in the experiment (including the results of all segments of the study). This aims to show a general summary of the variables under study, as well as the measurements of the data.

Analysis Result
After obtaining the results, two analyzes are conducted that compare the behavior of the segments in the experiment with the proposed model. The first model contrasts experimental and theoretical data, while the second one assesses controller performance during the experiment.
In this experiment, one of the main characteristics that can be seen corresponds to the high deviation that is obtained between the sample mean referring to standard deviation. These percentage differences found are y exp = 18.41% of HR acceleration, heart rate of h exp = 16.48%, and roller speed of v exp = 37.5%. This is mainly due to differences in the behavior of the variables for each of its segments.

Model Analysis
This analysis compares the time-based acceleration responses of the proposed HR with respect to the transient response obtained in the experiment. The experiment defined segment A for the transient study aimed at validating the model. In this segment, the individual emerges from inertia and reaches a maximum effort point. This translates into the effect on the response of HRV as HR acceleration and generated speed.
To establish the response differences between the proposed model y(t) and the response obtained experimentally y exp (t), the RMSE (Root Mean Square Error) index was used; this represents the absolute fit of two functions and has the property of being in the same units as the response variable. Lower values of RMSE fit the data vectors better. This index is a good measure of the precision with which the model predicts the response (evolution in t). In addition, the R 2 fitness index and Pearson and Kendall correlation coefficients are used to compare the responses. The former is used to test the direct relationship of two continuous variables, while the latter is used to demonstrate the statistical dependence of two variables.
To test the model, the experimental response y exp (t) is considered in the time domain. In turn, the proposed model, as a transfer function, is in the frequency domain and therefore needs to transit to the time domain, so both responses belong to the same domain.
As the proposed theoretical model is a transfer function H(s), and X(s) = u(s) is an input signal (step type) in this system, it is feasible to obtain the response Y(s), based on Equations (6) and (11). It is then observed that the theoretical time response is developed according to Equation (27 From this, Equation (28) is obtained.
By applying Laplace's inverse transform to function Y(s), it is assumed that this function can be considered a closed curve function and the poles of this function are singularities of the same. Then, the residue theorem can be applied to L(Y(s)) −1 . Consequently, the time response y(t) was obtained, which is necessary to compare this output with the experimental response y exp (t). The definition is presented in Equation (29).
where Res(Y(s), Z k ) is the residue of the function Y(S) in the pole Z k .
After calculating Equation (30), the response of the system y(t) is obtained (Equation (31)).
Knowing the temporary response of the HR acceleration model, the transient analysis period should be established. This period is defined as the point where the amplitude differences of both responses are null, as shown in Equation (32), where y ss (k) was the difference response.
Applying the condition to Equation (32) the transient period for the validation of the theoretical model y(t) with respect to the experimental data y exp (t) as T = 8.5 s. Figure 14 shows the comparison between the signals y(t) and y exp (t) in relation to the transient period T. The experimental response for the HRV behavior in the A segment of the experiment has been obtained.  Table 6 shows the results of the assessment indices used in the analysis of the model. The RMSE obtained from the comparison of the signals y(t) and y exp (t) was RMSE = 4.35%, indicating the difference between the amplitudes of both functions and the periods of maximum energy consumed in the change of HR acceleration. This is good due to the statistical basis on which the HR acceleration model H(s) was developed in Equation (11), and represented in terms of time by y(t) in Equation (31). The fit of both curves, R 2 = 94.65% shows the coherence of directions and duration periods of the transient in both signals. Regarding the correlation, Pearson and Kendall ρ show a strong and positive correlation between both responses. Being a negative correlation, zero indicates no correlation and 1 pointing to a perfect correlation. In this case, both are higher than 60%, which implies a strong correlation and coherence of the responses, reinforced by a positive covariance. Therefore, the proposed HRV model is validated, as well as the method to obtain the same.

Controller and Actuators Analysis
In this section, the responses of the experiment's segments B and C are compared. In segment B, the individual was active depending on his effort, trying to increase speed. Conversely, in segment C, the individual will try to increase the speed, but with the HRV controller activated.
The analysis consists of comparing the descriptive statistics with the experimental response in segment B with that in segment C to study the behavior of each response by segment, superposing them. The analysis is complemented by the study of the error over time e(t) of the responses. Table 7 shows the descriptive statistics of segments B and C of the experiment, including the calculation of the continuous component of each variable, which will be useful for the calculation of the error over time function e(t).  Figure 14 and Table 7 show that in segment B, the magnitude of the means with respect to HR acceleration is 91.1% larger than the magnitude of acceleration in the controlled segment C. This implies that the acceleration means used in segment B are much higher than in segment C, resulting in high dispersion in HR acceleration and consequently higher power dissipation. Regarding the C segment, minimal dispersion that reaches a deviation from the mean of only 3.54% is observed. Conversely, the dispersion of segment B corresponds to 7.2%, demonstrating the action of the controller and the regulation of HR acceleration through biofeedback.
The behavior of the error over time e(t) was individually calculated for the responses y exp (t) of the B segment and for y exp (t) of the C segment for a period equal to T = 158 s, which is the duration time of the controlled segment C.
Equation (33) is the expression of e(t) where y exp are the independent experimental responses in segments B and C. The acceleration value of HR that represents the population under study, defined in Equation (6), corresponds to u(t) = 0.00323 mm/s 2 , which is a reference value for the study of error behavior in experimental responses (value of reference variable).
The error behavior is shown in Figure 15, in which the responses y exp (t) for the segment B and C, as well as the reference variable u(t) superposed with one another. The difference between the behavior of the system response and the estimated reference variable u(t) = 0.00323 mm/s 2 is highlighted in this figure. The output of the B y exp−B segment exhibits a 4% difference with u(t). Likewise, y exp−C , which is the output of segment C, presents a 15% difference with u(t), which is generated by variations of the mean of the reference variable (continuous) and the difference of the former with the reference variable u(t). In segment B, the output has a sample mean of µ B = 0.079 mm/s 2 and a standard deviation of σ B = 0.133 mm/s 2 , with µ B < σ B . This implies a high variance in the response of this segment compared to its principal component due to the absence of control mechanisms for the acceleration of HR (HRV) in this part of the experiment. However, segment C and its output have a sample mean of µ C = 0.016 mm/s 2 and a standard deviation of σ C = 0.0114 mm/s 2 , with µ C > σ C . This results in low variance with respect to the principal component in the response of this segment, which is a consequence of the action of the controller on the system and reveals a change in the HRV response of the individual.
Another piece of evidence supporting the above corresponds to the difference between the means of each output, which corresponds to 79.7%. In addition, the difference between the standard deviations is practically 100%, i.e., in order of magnitude, the non-controlled segment B does not control the energy consumed, which implies that there is no control besides homeostatic control to maintain HR when experiencing signs of tiredness and fatigue. Instead, segment C, as it has lower variability in its central statistics, shows that the energy applied is much lower than that of segment B to increase HR speed slowly and linearly without sudden changes in HRV expressed as the acceleration of HR.

Conclusions
Through this research, a reliable and novel method has been presented for the estimation and identification of a dynamic model for the representation of HRV and its control, using light actuators and an MRAC-type control strategy.
We believe that this research is very valuable in the area of biofeedback, providing new mechanisms for the self-control of physiological variables such as HRV to avoid and improve the well-being and safety of individuals. Having the ability to use this technique as an anticipatory method of HRV behavior patterns under given conditions makes this research an important contribution to the area of engineering, medicine, and psychology.
From the results obtained, the behavior of the statistically obtained HRV model was compared with the estimated model based on the experimentally obtained data. This dynamic model represented in the form of a transfer function was subjected to reactive type tests such as the step response test. These responses were counteracted with the results applied to an individual not belonging to the initial group used for the system modeling, using a step test, where the individual emerged from inertia until a moment of stability with the maximum effort applied, obtaining a direct correlation R 2 = 89.5%, verifying that the model is representative of a population group.
In reference to the self-control of the HRV, it is established that for a given test, the error in the time reached by using the controller and the signals from the light actuators is only 13.3% compared to the uncontrolled response that reached a dispersion close to 100%, verifying the effectiveness in the use of biofeedback through a continuous time controller with fast, noninvasive actuators such as light actuators.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Conflicts of Interest:
The authors declare no conflict of interest.