Simulating Dynamics of Circulation in the Awake State and Different Stages of Sleep Using Non-autonomous Mathematical Model With Time Delay

We propose a mathematical model of the human cardiovascular system. The model allows one to simulate the main heart rate, its variability under the influence of the autonomic nervous system, breathing process, and oscillations of blood pressure. For the first time, the model takes into account the activity of the cerebral cortex structures that modulate the autonomic control loops of blood circulation in the awake state and in various stages of sleep. The adequacy of the model is demonstrated by comparing its time series with experimental records of healthy subjects in the SIESTA database. The proposed model can become a useful tool for studying the characteristics of the cardiovascular system dynamics during sleep.


INTRODUCTION
The study of the dynamics of the circulatory system during sleep attracts a lot of attention (de Zambotti et al., 2018;Tall and Jelic, 2019;Kontos et al., 2020;Nano et al., 2020). It was shown that scaling properties of the cardiac dynamics are different during sleep and wake periods (Ivanov et al., 1999). Characteristics of rapid eye movement (REM) sleep are correlated with exacerbation of ischemic heart disease, sometimes even leading to myocardial infarction (Schumann et al., 2010). These events are most common in patients with coronary artery disease (Nowlin et al., 1965). Of particular importance is the change in the dynamics of the loops of the autonomous control of blood circulation, due to the influence of the higher nervous centers on them (van Roon et al., 2004;Ivanov, 2006). In experimental studies (Nowlin et al., 1965;King et al., 1973), it was shown that in animals with aortic stenosis, sympathetic activation can lead to decreased myocardial perfusion. The development of pathologies such as apnea is associated with dysfunction of the autonomic control of the cardiovascular system (Aydin et al., 2004;Khoo and Blasi, 2013;Roder et al., 2018).
Despite numerous studies, the functioning of the cardiovascular system (CVS) during sleep is essentially terra incognito. This is due to the complexity of CVS, which includes a large number of interacting non-linear elements, as well as technical difficulties and ethical limitations of experimental studies. Therefore, the development of mathematical models based on "first principles, " i.e., using physical and physiological laws, is of great importance (Wolkenhauer, 2014). Such models allow one to simulate complex CVS dynamics and generate stationary time series of any length (Karavaev et al., 2019). Using models, it is possible to simulate various pathologies (Karavaev et al., 2016), drug effects (Karavaev et al., 2016), physiological tests (Ishbulatov et al., 2020), changes in physiological conditions (van Roon et al., 2004), and the influences of brain activity (van Roon et al., 2004) and respiration (Wessel et al., 2009;Cheng et al., 2010). The development of such models is promising for solving the problems of personalized medicine, when some parameters of the model can be estimated directly from the experimental data of a particular patient. The use of personalized mathematical models expands the possibilities of medical diagnostics and therapy, making it possible to predict the course of diseases and simulate the patient's response to medications.
A number of mathematical models of CVS are known, which take into account the dynamics of the loops of autonomic control of blood circulation (Guyton et al., 1972;Ivanov et al., 1998;Seidel and Herzel, 1998;Ursino, 1998;Ottesen, 2000;Kotani et al., 2005). However, only a few studies are aimed at modeling the dynamics of autonomic control loops and CVS during sleep (van Roon et al., 2004). The most well-known model of CVS during sleep is PNEUMA (Cheng et al., 2010). It is the eleventhorder system, which contains more than 80 algebraic terms and more than 200 parameters, most of which are estimated empirically and have no physical or physiological meaning. The PNEUMA model has shown its adequacy and importance . However, its complexity complicates the interpretation of the results. Moreover, the model is too cumbersome to be personalized and fitted to a particular patient.
A number of studies suggest that adequate simulation of CVS dynamics in different stages of sleep requires taking into account the effects associated with autonomic control. In experimental studies, it was shown that non-REM sleep (especially stage 4) is characterized by a decreased tone of sympathetic activity, a decrease in heart rate and average level of blood pressure, and a decrease of blood pressure variability (Sommers et al., 1993). Other studies reported pronounced respiratory sinus arrhythmia caused by parasympathetic activity (Zemaityte et al., 1984).
Rapid eye movement sleep is characterized by more complex dynamics, for which epochs of sharp increase in sympathetic activity, alternating with intervals of decreased sympathetic activity, are typical (Sommers et al., 1993). As shown by direct measurements from sympathetic nerves, on average, there is an increase in sympathetic activity, leading to an increase in blood pressure variability, but the heart rate corresponds to a waking person at rest. These results were confirmed in active experiments with blocking of sympathetic or parasympathetic autonomic control, which compensated for the corresponding changes in the dynamics of CVS (Zemaityte et al., 1984).
It is known that the transition to sleep is associated with the activity of the higher nervous system. In model studies (van Roon et al., 2004) and (Cheng et al., 2010), potential ways for the influence of higher nerve centers on autonomic control during sleep are presented.
In the present paper, we propose a more compact model describing the dynamics of CVS in different stages of sleep and during wakefulness. In contrast to our earlier models, the new model for the first time takes into account the effect of the cerebral cortex on blood circulation in the awake state and in sleep. The proposed model consists of four differential equations with time delay. The model has 55 parameters, 39 of which have physiological meaning and can be estimated experimentally. Despite the relatively simple structure, the model takes into account the non-linear properties of the autonomic control and simulates with good accuracy the time series of real arterial blood pressure and interbeat intervals of healthy resting subjects. As well as our earlier models, it reproduces the pathological changes that lead to arterial hypertension (Karavaev et al., 2016), reaction to the passive orthostatic test (Ishbulatov et al., 2018), and reaction to the autonomic blockade due to administration of Arfonad (Karavaev et al., 2016). Moreover, the model simulates the chaotic dynamics of heart rate (Karavaev et al., 2019) and phase synchronization between the autonomic control loops (Ishbulatov et al., 2017), which is observed in humans (Bartsch et al., 2012) and is important for diagnostics and understanding of some circulatory diseases (Prokhorov et al., 2003;Kiselev et al., 2016aKiselev et al., ,b, 2020. The paper is devoted to the description of the proposed model and the study of its dynamics in the awaking state, during REM sleep, and during non-REM sleep. We also compared the simulated data with the experimental data from the known SIESTA database (Klösch et al., 2001).

Study Participants
Our study included recording of 20 healthy subjects that fulfilled the following criteria: no regular shift work, usual bedtime before midnight, and no acute depressive or anxious symptoms. Each subject was monitored by wristworn actigraphs (Actiwatch, Cambridge Neurotechnology, England) one-week prior and one-week after the recording session. Quality control ensured that there are no outliers in the final data set, and all signals were recorded with no protocol violations.
For each subject, we analyzed a set of three 20-minute ECG signals, one recorded in the waking state, one during REM sleep, and one during stage 4 of non-REM sleep. Each signal was recorded at the sampling rate of 200 Hz. The high-pass cutoff frequency was set between 1.6 Hz and 16 Hz.

Mathematical Model
To modify the previously proposed model (Karavaev et al., 2019), we used the ideas proposed in Cheng et al. (2010) and van Roon et al. (2004). In the equations of autonomic control of blood circulation, we added eight parameters characterizing the influence of the higher nervous centers. These parameters were set equal to zero for the awake state, but they took non-zero values during REM and non-REM sleep in order to simulate the increase and decrease in autonomic control activity. Figure 1 depicts the scheme of the model, where the places where the parameters are added are shown in red.
We used integrate-and-fire model to simulate the heart rate: where ϕ (t) is the phase of the sinoatrial node, T 0 = 1.5 s is the heart rate, f s (t) and f p (t) are the sympathetic and parasympathetic factors, respectively, which modulate T 0 , and ξ is the 1/f noise (Ivanov et al., 1998), which is added to simulate myocardial and humoral control. Spectral properties of noise were chosen to match experimental signals (Bezerianos et al., 1995;Goldberger, 1996). In the absence of noise or autonomic control, the sinoatrial node generates periodic saw-like impulses with the period T 0 . Within the first T sys seconds (0.125 s) after the initiation of heart cycle, the arterial pressure rapidly grows and is described by the following equation: where D i−1 is the systolic pressure at the end of the previous cardiac cycle, T i−1 is the moment of time when the previous cardiac cycle ended, T sys is the duration of the current systolic phase, t i is the time since the beginning of the current cardiac cycle, B(t) is the respiratory signal, k B p is a non-dimensional amplitude of respiration, and S(t) is the cardiac contractility (Wessel et al., 2009) defined as follows: whereŜ is the resting heart contractility, n c is a fitting parameter, and where S 0 is the contractility of the denervated heart, c c (t) is the concentration of noradrenaline in blood that circulates in the heart, c v (t) is the concentration of noradrenaline in vessels, k c S and k v S are the parameters that characterize sensitivity of contractility to changes in the noradrenaline concentration, k t S is the parameter that characterizes sensitivity of contractility to changes of the heart rate, and L i−1 is the duration of the previous cardiac cycle. The values of all model parameters used for simulation are presented in Table 1.
Respiratory signal B (t) is defined as follows: where T br is the mean duration of respiratory cycle and ζ is a noncorrelated Gaussian zero-mean noise that is used to modulate the rate of respiration after each cycle. After T sys seconds, the systolic phase of cardiac cycle ends and arterial pressure begins to slowly decrease in accordance with the Windkessel model (Frank, 1899): where C is the parameter that reflects elastic properties of aorta and R(t) is the peripheral resistance, which is modulated by the sympathetic activity: where R 0 is the resistance of unstressed vessels and k v R is the parameter that describes sensitivity of peripheral resistance to changes of noradrenaline concentration in vessels walls c v (t).
The absolute value and rate of change of arterial pressure (Warner, 1958) are sensed by two baroreceptor nodes located in carotid sinus and in the major vessels of the lower body. Activities of these baroreceptor nodes [v b (t) and v l b (t), respectively] are defined as follows: where p 0 and p l 0 are the lowest pressure the baroreceptors react to and k 1 , k 2 , k l 1 , and k l 2 are the sensitivity of arterial baroreceptors to arterial blood pressure and the rate of its change.
The higher nervous centers process the inputs of baroreceptors and adjust the activity of heart rate autonomic control v s (t) and vessel tone autonomic control v l s (t). To describe their activity, we used equations with time delay and sigmoidal non-linearities as in Ringwood and Malpas (2000): where C b s , C lb s , C v s , and C lv s are the parameters, which characterize the influence of higher nervous centers, a s , b s , y s , a l s , b l s , and y l s are the parameters of the non-linear transfer function in the feedback loops of the baroreflex control, v 0 s is the resting afferent tone of the heart rate sympathetic control, v l0 s is the resting afferent tone of the vascular tone sympathetic control, and k r s and k lr s are the parameters characterizing the influence of respiration. Activity of parasympathetic autonomic control v p (t) is described as follows: where v 0 s , v l0 s , and v 0 p are the sympathetic and parasympathetic activity under resting conditions, C v p and C k p are the parameters that reflect the influence of higher nervous centers on the parasympathetic control of heart rate, k b p characterizes the influence of baroreflex activity, and k r p characterizes the influence of respiration. Changes in the levels of sympathetic activation affect cardiac concentration of noradrenaline c c (t) and vascular concentration of noradrenaline c v (t): where τ c and τ v are the time constants, θ c and θ v are the time delays caused by finite speed of neural transition and secretion of noradrenaline, and k S c , k S v , and k v are the transfer coefficients. Changes of concentration of noradrenaline and activity of parasympathetic control affect the heart rate through the sympathetic factor f s (t) and parasympathetic factor f p (t): where θ p is the time delay, C s ϕ and C p ϕ are the parameters that reflect the influence of higher nervous centers on the parasympathetic and sympathetic factors of heart rate control, respectively, and k c ϕ , k p ϕ ,ĉ c , n s ,v p , and n p are the dimensionless factors. No separate equation was introduced to model the changes in concentration of acetylcholine (parasympathetic transmitter) because the rate of its secretion is faster than the dynamics on which the model is focused. The phase effectiveness curve represents the changes in the sinoatrial node sensitivity to parasympathetic control throughout the cardiac cycle: where ϕ is the phase of the cardiac cycle.
To simulate the sleep stages and awake state, we set the corresponding parameter values and generated long time series with the fixed parameters. To exclude the transient process, we did not analyze the first 1,000 s of time series. The length of the model time series was equal to the length of the experimental time series from the SIESTA database.

Methods
The model data were compared with the experimental data from the SIESTA database using statistical measures, spectral analysis, and calculation of the largest Lyapunov exponent. The spectral analysis was carried out for the interbeat interval time series, which were extracted from the ECG signals by detecting the duration of time intervals between successive R-peaks. Since the dynamics of the heart in our model is described by a simple integrate-and-fire Eq. 1, in which the moment of sinoatrial node excitation corresponds to the reset of the cardiac cycle phase to zero, the location of R-peaks was defined by detecting the moment of such reset of the phase. We analyzed the spectral power in the ranges 0.05-0.15 Hz and 0.15-0.4 Hz to calculate the low-frequency (LF) and high-frequency (HF) indices, which reflect the activity of the sympathetic and parasympathetic control, respectively (Heart Rate and Variability., 1996).
Complexity is another important characteristic of CVS dynamics (Porta et al., 2017). To estimate the complexity, we calculated the largest Lyapunov exponent (Rosenstein et al., 1993) from the interbeat interval time series filtered in the range 0.05-0.4 Hz. The largest Lyapunov exponent was estimated using the Rosenstein algorithm (Rosenstein et al., 1993), which is well suited for the analysis of short time series. The first step was to find a nearest neighbor for each point of the reconstructed phase space, but the neighbors close in time were excluded from the analysis (Rosenstein et al., 1993). In dynamical systems, Sensitivity of contractility to changes in the heart rate the average rate of divergence of the trajectories of the nearest neighbors obeys the following expression: where L 0 is the initial distance between the nearest neighbors, λ 0 is the largest Lyapunov exponent, and t is the time of calculation. Then, λ 0 is defined as follows: To reconstruct the phase space, we used the method of time delays. The dimension of the phase space was equal to 13. The time series were analyzed in windows having the length 1,000 s. A detailed explanation of the choice of this set of parameters is presented in Karavaev et al. (2019).

RESULTS
The mathematical model was used to simulate the healthy subjects in the awake state, during REM sleep, and stage 4 of non-REM sleep. The parameters of autonomic control were constant throughout all stages. We only changed the values of those parameters that were introduced to take into account the inputs from the higher nervous centers. The values of these parameters for each studied stage are presented in Table 1. Figure 2 shows the model and experimental time series of arterial pressure and interbeat intervals. It can be seen in Figure 2 that arterial pressure variability and interbeat interval variability are most pronounced in the awake state. These results agree well with the data from literature (Zemaityte et al., 1984;Heart Rate and Variability., 1996) and reflect the higher activation of sympathetic control during the wakefulness.
The power spectra of model and experimental interbeat intervals are shown on Figure 3. In the model power spectra, the peaks at the frequencies of about 0.1and 0.29 Hz are significantly more narrow, sharp, and pronounced. This is explained by the fact that our model does not take into account the humoral control of CVS. However, the model qualitatively simulates the difference between the experimental spectra during wakefulness and at different stages of sleep. The power of the 0.1-Hz spectral component, which is associated with the sympathetic activation, takes the highest value in the awake state. The power of the 0.29-Hz component, which is associated with the parasympathetic control, takes the highest value during the stage 4 of non-REM sleep. These results are in good agreement with the known results of other studies (Zemaityte et al., 1984;Sommers et al., 1993).
The spectral and statistical indices calculated from model and experimental data are presented in Table 2.

DISCUSSION
PNEUMA is one of the most famous models that simulate the influence of higher nervous centers on the dynamics of autonomic control during sleep stages and awake state (Cheng et al., 2010). The PNEUMA model takes into account many factors, such as respiration, chemoreceptors, and blood hydrodynamics. As a result, the PNEUMA model is highdimensional and complex, which makes it difficult to interpret the results of model studies and observed effects. The high order of the model makes it almost impossible to reconstruct the model parameters from the data of a particular patient and to individually fit the model. The authors themselves stated that there is no available database that can be used to verify their model, and many parameters do not have physiological meaning (Cheng et al., 2010).
We have proposed a mathematical model that allows one to simulate the heart rate variability and arterial pressure oscillations under the influence of respiration and autonomic control loops. Earlier, it was shown that taking into account   the non-linear properties and self-oscillatory dynamics of these loops it is possible to explain the experimentally observed complex chaotic dynamics of the heart rate and synchronization processes within the CVS (Ringwood and Malpas, 2000;Prokhorov et al., 2003;Karavaev et al., 2009Karavaev et al., , 2016Karavaev et al., , 2019Kiselev et al., 2016aKiselev et al., , 2020. In the present paper, a modified model is proposed that takes into account the influence of the activity of the cerebral cortex on autonomic control in the states of sleep and wakefulness. The respiratory process and some other factors are included into the model in a simplified form. Nevertheless, the proposed model allows us to simulate complex oscillatory dynamics of the interbeat intervals associated with the autonomic control of circulation. For simplicity of the model, we have accepted a number of limitations: the model does not account for humoral control (Ganten and Stock, 1978;Chopra et al., 2011), blood hydrodynamics, and local control mechanisms (Clifford, 2011;Hong et al., 2020), such as the Bowditch effect (Usman et al., 2020). This explains the lower λ 0 in the model compared to the experimental data. The narrower peaks in the model power spectrum compared to the experimental one (Figure 3) can be explained in a similar way. Moreover, the model is stationary, while real CVS is not.
Despite the aforementioned limitations, we consider the model to be an adequate representation of the dynamics of autonomic control during sleep and wakefulness. For experimental interbeat intervals, LF index, which is associated with the activity of the sympathetic autonomic control, typically takes the larger value in the awake state than in REM sleep and the smaller value in non-REM sleep than in REM sleep (Kantelhardt et al., 2002). The proposed model shows similar properties of the LF index. The model HF index, which is associated with the activation of the parasympathetic control, takes the largest value in the 4th stage of non-REM sleep. We obtained similar results for many experimental signals.
Our experimental data did not contain records of arterial pressure, but the known results (Sommers et al., 1993) report that both systolic and diastolic arterial pressure during REM sleep is slightly lower or the same as in the awake state, while during non-REM sleep, the arterial pressure is significantly lower. The results of our model agree well with this experimental observation.
Estimations of the largest Lyapunov exponent λ 0 from the model and experimental data in both REM and non-REM sleep are in a good agreement. However, in the awake state, the estimations of λ 0 differ for the model and experimental data. This disagreement is due to the simplicity of the model compared to the real system. Some more compact models use stochastic components to simulate the elements of blood circulation (Ivanov et al., 1998;Kantelhardt et al., 2003). Such approaches may have advantages in the compactness of equations, but they do not allow one to study the non-linear dynamics of blood circulation and build personalized models.
In further research, we plan to develop a method for reconstructing the parameters of our model from experimental signals of ECG, blood pressure, and respiration. For example, the delay times of model equations can be estimated using different techniques (Prokhorov et al., 2005;Sysoev et al., 2016). However, to propose such a technique, it is necessary to have a good understanding of the capabilities and limitations of the model. In the present study, we took a step in this direction by taking into account the influence of brain activity on CVS during sleep and wakefulness. Such personalized mathematical model can be used in the future for the development of personalized medicine and will provide a new tool for the study of autonomic control.

CONCLUSION
We have proposed a mathematical model of autonomic control in a form of fourth-order system of non-autonomous differential equations with time delay. Despite its compact structure, the model is able to simulate the dynamics of autonomic control during sleep stages and in awake state and related changes in interbeat intervals. The adequacy of the proposed model is demonstrated by comparing its time series with experimental records of healthy subjects in the SIESTA database.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: These are data which belong to medical faculties and are not publicly available. Requests to access these datasets should be directed to TP (thomas.penzel@charite.de), he will provide a link to people at The SIESTA group in Vienna, Austria.

ETHICS STATEMENT
For testing the model, it was compared to the experimental electrocardiogram (ECG) signals from the SIESTA database (Klösch et al., 2001). The studies involving human participants were reviewed and approved by the Ethics Committee of Klinikum der Philipps-Universität Marburg, Germany. The study participants, all above age of 18 years, provided written informed consent to participate in the study.

AUTHOR CONTRIBUTIONS
AK and TP finalized the manuscript for submission. All authors contributed to writing and discussion on the manuscript.