Contribution of Cardiorespiratory Coupling to the Irregular Dynamics of the Human Cardiovascular System

Abstract

Irregularity is an important aspect of the cardiovascular system dynamics. Numerical indices of irregularity, such as the largest Lyapunov exponent and the correlation dimension estimated from interbeat interval time series, are early markers of cardiovascular diseases. However, there is no consensus on the origin of irregularity in the cardiovascular system. A common hypothesis suggests the importance of nonlinear bidirectional coupling between the cardiovascular system and the respiratory system for irregularity. Experimental investigations of this theory are severely limited by the capabilities of modern medical equipment and the nonstationarity of real biological systems. Therefore, we studied this problem using a mathematical model of the coupled cardiovascular system and respiratory system. We estimated and compared the numerical indices of complexity for a model simulating the cardiovascular dynamics in healthy subjects and a model with blocked regulation of the respiratory frequency and amplitude, which disturbs the coupling between the studied systems.

1. Introduction

The origin of the irregularity or complexity of the heart period has been actively discussed for more than 30 years [1,2]. The answer to the question about the origin of irregular dynamics in the cardiovascular system (CVS) is fundamentally important for understanding the CVS behavior [3,4,5]. Moreover, the measures of complexity correlate with the physiological condition of the CVS [6,7,8,9,10] and are important for medical diagnostics of several cardiovascular diseases, including arterial hypertension [11], coronary heart disease and mortality from other causes [12]. Several authors have shown that the largest Lyapunov exponent and fractal dimension are perspective measures for medical diagnostics and therapy of the CVS [13,14,15]. However, researchers still debate the source of complexity in the CVS.

The discussion is held around three possible sources of complexity in the CVS: deterministic chaotic dynamics of the elements of the CVS, noise of a “central” origin and the irregularity of respiration and cardiorespiratory coupling.

One of the possible sources of chaotic dynamics in the CVS is the deterministic chaotic activity of the autonomic control of circulation [7,16,17,18]. Ernst [19] hypothesized that the irregularity of the heart period is caused by the influence of the higher nervous centers. Earlier, we conducted a study on a mathematical model of the CVS [20] and showed that the autonomic control of the circulation exhibits deterministic chaotic dynamics even without the influence of the higher nervous centers.

Another possible source of the complexity is the heart itself. Experimental studies have proved that disturbance of myocardial conductivity could lead to the generation of chaotic dynamics. However, this behavior is caused by pathological changes in the heart and is not typical for patients without cardiovascular diseases [21,22,23,24].

Noise also affects the dynamics of the CVS. It was shown that the power of noise in the sequence of intervals between the R peaks of the ECG signal (RR intervals) decreases in deep sleep. Since deep sleep is associated with a weaker influence of the central nervous system on the circulation [2,25], the central nervous system is hypothesized to be the source of this noise.

Cardiorespiratory coupling is another important factor in the dynamics of the CVS [26,27,28,29,30]. Several authors suggested [3] that the irregularity of respiration and the nonlinear nature of cardiovascular coupling are potential sources of the complexity in the heart period.

In the literature, there are no quantitative estimations for the contribution of each aforementioned factor, and many researchers are skeptical about the validity of complexity measure estimation from experimental signals of biological origin, which are typically nonstationary and noisy. These properties of biological signals complicate the interpretation of the results and the experimental search for the source of complexity in the heart period [6,31,32,33].

Under such conditions, it seems promising to supplement experimental studies with investigations of mathematical models of the CVS. Mathematical models based on physical and physiological principles can be used to simulate experiments which are impossible in vivo. When investigating the source of complexity in the heart rate, a model could be used to simulate a series of active experiments with successive blockades of possible sources of complexity, and after each test, the resulting changes in the complexity indices can be measured.

Detailed models of particular elements and aspects of the CVS have been proposed, which model the autonomic control of mean arterial pressure [17,18], impulse propagation in the myocardium tissue [34], gas exchange in the respiratory system [35] and cardiorespiratory coupling [36]. There are also complex large-scale models that combine the autonomic and humoral control of the circulation, respiratory system and central nervous system [37,38]. Drawbacks of large-scale models are the problems with parameterization and the simplified simulation of particular elements of the CVS. These shortcomings limit the ability of such models to simulate several important effects, such as coupling between the low-frequency (~0.1 Hz) oscillations in RR intervals and arterial pressure [39].

The model we developed for this study is intermediate in its complexity. The model provides a detailed simulation of the autonomic control of the circulation, following the models of [40,41], and also takes into account the nonlinear, self-exciting nature of the autonomic control [17,18]. The results of Ben-Tal et al. [42] were the basis for our model of the respiratory system. The model takes into account the influence of CO₂ and O₂ concentrations in the arterial blood on the amplitude and frequency of respiration. Cardiorespiratory coupling was implemented following the Magosso and Ursino model [38], in which the activity of the autonomic control is a function of the activity of not only the baroreceptors but also the carotid and central chemoreceptors.

Earlier, we conducted a model study [20] which provided further arguments towards the hypothesis that the complexity of the heart period originates from chaotic dynamics of the autonomic control and “central” noises. We also found no significant effects of respiration on the complexity of the heart period. However, this earlier model simulated respiration as an external sine wave generator with a randomized frequency, which was independent of the dynamics of the CVS.

In the present model study, a more detailed simulation of the respiratory system allowed us to analyze the contribution of nonlinear bidirectional coupling between the respiratory system and the CVS to the complexity of the heart period. We investigated this problem by blocking the regulation of the respiratory frequency and amplitude, thus disturbing the coupling, and then measuring the changes in the complexity measures.

2. Materials and Methods

2.1. Mathematical Model

Figure 1 illustrates the structure of our model. Section 2.1.1, Section 2.1.2 and Section 2.1.3 of this article provide a detailed description of the equations, and

Table 1. Parameters of the model.

Parameter	Value	Parameter	Value	Parameter	Value
${\tilde{g}}_{nap}$	133.33	$p_o^0$	40	$n_v$	1.5
${\tilde{g}}_L$	98	$L$	1.712 × 108	$k_1^b$	0.55
${\tilde{b}}_L$	0.6667	$K_T$	1 × 104	$k_2^b$	6.5 × 10−2
${\tilde{E}}_L$	0.212	$K_R$	3.6 × 106	$k_3^b$	100
${\tilde{\mathsf{\theta }}}_{hp}$	0.313	$f_{om}$	0.21	$k_4^b$	40
${\tilde{\mathsf{\sigma }}}_{hp}$	0.04	$f_{cm}$	0	$k_5^b$	0.47
${\tilde{\mathsf{\tau }}}_{hp}$	10	$k_1^c$	0.585	$p_0$	50
${\tilde{\mathsf{\theta }}}_{mp}$	0.367	$k_2^c$	61	$v_{\mathsf{\alpha }s}^0$	1.2
${\tilde{\mathsf{\sigma }}}_{mp}$	−3.3 × 10−2	$k_3^c$	20.7	$v_{\mathsf{\beta }s}^0$	0.94
$T{r}_1$	0.35	$k_4^c$	0.442	$k_{\mathsf{\alpha }s}^b$	2.5
$T{r}_2$	0.35	${\mathsf{\theta }}_{vc}$	5	$k_{\mathsf{\alpha }s}^c$	1
$T{r}_3$	0.5	$k_1^o$	−0.5	$k_{\mathsf{\alpha }s}^r$	0.1
$T{r}_4$	1 × 10−4	$k_2^o$	−50	$k_{\mathsf{\beta }s}^b$	2.5
ε	1 × 10−5	$k_3^o$	20	$k_{\mathsf{\beta }s}^c$	1.5
$I_L$	0.99999	$k_4^o$	0.496	$k_{\mathsf{\beta }s}^r$	0.1
$k_1^x$	2	$k_5^o$	−5 × 10−5	$v_p^0$	0
$k_2^x$	1	$k_6^o$	5.3 × 10−2	$k_p^b$	1.3
$P_m$	760	$k_o^g$	60	$k_p^c$	0.6
$P_{L0}$	4.5	$k_c^g$	20	$k_p^r$	0.3
$k_p$	2.5	$k_c^K$	25	$k_p^{\xi }$	5 × 10−2
$E$	2.5	$k_o^K$	7.5	$a$	0.97
$R_a$	1	$g_0$	16	${\mathsf{\tau }}_c$	1.8
$D_c$	7.08 × 10−3	$K_0$	0.75	${\mathsf{\tau }}_v$	2.0
$D_o$	3.5 × 10−4	$k_1^{I_k}$	0.1	$k_c^S$	1.2
$P_w$	47	$k_2^{I_k}$	0.1	$k_v^S$	1.2
${\mathsf{\sigma }}_c$	3.3 × 10−5	$T_0$, s	0.78	${\mathsf{\theta }}_c$	2.65
$\mathsf{\sigma }$	4 × 10−2	$T_{sys}$, s	0.125	${\mathsf{\theta }}_v$	2.65
$V_c$	7 × 10−2	$k_p^B$	0.5	$k_{\mathsf{\phi }}^c$	2.1
$C_u$	25.45	$S_0$	25	$k_{\mathsf{\phi }}^p$	9.5
$T_h$	2 × 10−3	$\widehat{S}$	70	${\widehat{c}}_c$	2
$r_2$	0.12	$k_S^c$	40	${\widehat{v}}_p$	2.5
$l_2$	1.64 × 105	$k_S^t$	10	$n_s$	2
$H$	10−7.4	${\mathsf{\tau }}_v^{(0)}$	2.2	$n_p$	2
$\mathsf{\delta }$	101.9	${\overline{\mathsf{\tau }}}_v$	1.5	${\mathsf{\theta }}_p$	0.5
$p_c^0$	46	${\widehat{c}}_v$	10

All parameters listed in Table 1 are explained in Section 2.1.

lists the parameters of the model.

We propose a model which consists of two main parts, where one part models the cardiovascular dynamics, and the other part models the respiratory system. The equations modeling the cardiovascular dynamics are based on the models proposed in [20,41,42]. We used two equations to model the heart. The “integrate and fire” generator set the heart frequency, which was modulated by the sympathetic and parasympathetic loops of the autonomic control of the circulation. The second equation describes the heart contractility, i.e., how much blood it pumps during one cardio cycle. The heart contractility is also a function of the autonomic control. The model does not simulate the blood hydrodynamics; it only simulates blood pressure in the aorta using two equations. The first equation describes the rapid growth of the blood pressure during the first 0.125 s of the cardio cycle, which is caused by the blood being pumped into the aorta. The second equation describes the slow decrease in blood pressure, caused by the blood flowing into the peripheral arterial vessels.

The focus of the model was on the simulation of the autonomic control of the circulation, which includes the sympathetic control of the heart period and arterial pressure, and the parasympathetic control of the heart period. The sympathetic control is activated when the blood pressure drops below normal values. Activation of the sympathetic control increases the frequency of heart contractions and heart contractility, and contraction of arterial vessels. The combination of these factors causes the blood pressure to rise. Activation of the parasympathetic control has the opposite effect and is caused by the arterial pressure exciding the normal levels. The arterial pressure is measured by the baroreceptors, located in the aorta.

The sympathetic control is simulated using delayed negative feedback loops. The general structure of the autonomic control was adopted from [40,41]. However, we introduced a more detailed nonlinear transfer function for the baroreceptors. We also introduced separated control loops of the α-sympathetic and β-sympathetic regulation of the arterial pressure and heart period.

The model of the respiratory system is based on the model proposed in [42]. It consists of the central pattern generator, which dictates the frequency and amplitude of respiration; gas exchange between the alveoli and arterial blood; and the loops of the autonomic control of respiration. The autonomic control of respiration functions similarly to the autonomic control of circulation and is activated by the increase in the CO₂ concentration and the decrease in the O₂ concentration in the blood. Changes in the gas concentrations are sensed by the central and carotid chemoreceptors. We modified the models for the chemoreceptors by introducing more detailed transfer functions [43].

We introduced the cardiorespiratory coupling in the same way as it was done by Magosso and Ursino [38], where the activation of the central chemoreceptors excites both the sympathetic and parasympathetic control of the circulation.

The proposed model is a non-autonomous system of thirteen differential equations, several of which contain delayed arguments, and one equation contains a stochastic term. The equations were solved numerically, using the Euler approach with an integration step $\Delta t$ = 0.0001 s, as it was recommended in [42]; when switching from the integration step of 0.1 ms to the integration step of 0.2 ms, the average values of the considered indices change by less than 1%.

2.1.1. Central Pattern Generator of the Respiratory Rhythm

The model for the central pattern generator of the respiratory rhythm consists of two main parts. The first one—the pre-Bötzinger complex (pre-Bötc)—is a population of neurons in the medulla that activates and generates the spikes during inspiration. The second part—the rostral ventral respiratory group (rVRG)—integrates the spikes from the pre-Bötc complex and transmits them to the spinal cord, by which the electric signal reaches the neurons of the diaphragm. Then, the diaphragm contracts and initiates inspiration. During expiration, the central pattern generator of the respiratory rhythm is inactive, and elastic forces relax the diaphragm [43]. The model did not simulate the activity of the Bötzinger complex and also did not account for the pons and retrapezoid nucleus effects on the respiratory center, since these factors only play a major role during forced breathing [43], which is beyond the scope of the model. For the same reason, the model did not simulate the intercostal muscles. When describing the equations, we used terminology similar to that used in [42].

We adopted a model for pre-Bötc in such a way that it ensures an independent control of the respiratory frequency and depth. The variable $B(t)$ in Equation (1) denotes the mean instantaneous firing frequency of the pre-Bötc neuron population. $B(t)$ = 0 corresponds to the absence of any neuronal activity, and $B(t)$ = 1 corresponds to the synchronous firing of each neuron from the pre-Bötc population at the highest frequency. $B(t)$ is described by the following equation:

$$\dot{B}(t)=\mathsf{\alpha }(t)(1-B(t))-\mathsf{\beta }B(t)+\mathsf{\gamma },$$

(1)

where $\mathsf{\alpha }(t)(1-B(t))$ is the activation rate, $\mathsf{\beta }B(t)$ is the inactivation rate, $\mathsf{\gamma }$ is the external drive, which can inhibit or excite the population of neurons, $\mathsf{\alpha }(t)={\tilde{g}}_{nap}{\overline{m}}_p(t)h_p(t)$, $\mathsf{\beta }={\tilde{g}}_t(t)+{\tilde{g}}_L$ and $\mathsf{\gamma }=\tilde{b}{\tilde{g}}_t\left(t\right)+{\tilde{E}}_L{\tilde{g}}_L$ are the functions of the excitatory and inhibitory drives, ${\tilde{g}}_{nap}$,${\tilde{g}}_L$, ${\tilde{b}}_L$ and ${\tilde{E}}_L$ are parameters that link $B(t)$ to the average voltage, produced by the pre-Bötc neuron population [42], ${\tilde{g}}_t(t)$ is the variable regulating the respiratory frequency (a higher ${\tilde{g}}_t(t)$ leads to more frequent inspirations and a somewhat lower amplitude of respiration), $h_p(t)$ is the variable representing the inactivation gating of the persistent sodium current and ${\overline{m}}_p(t)$ is the variable representing the activation gating of the persistent sodium current.

$${\dot{h}}_p(t)={\mathsf{\alpha }}_{hp}(t)\left(1-h_p(t)\right)-{\mathsf{\beta }}_{hp}(t)h_p(t),$$

(2)

where

$${\mathsf{\alpha }}_{hp}(t)=\frac{\exp \left(-\overline{B}(t)\right)}{2{\tilde{\mathsf{\tau }}}_{hp}},$$

(3)

where $\overline{B}(t)=\left(B(t)-{\tilde{\mathsf{\theta }}}_{hp}\right)/2{\tilde{\mathsf{\sigma }}}_{hp}$,

$${\mathsf{\beta }}_{hp}(t)=\frac{\exp \left(\overline{B}(t)\right)}{2{\tilde{\mathsf{\tau }}}_{hp}},$$

(4)

where ${\tilde{\mathsf{\theta }}}_{hp}$, ${\tilde{\mathsf{\sigma }}}_{hp}$ and ${\tilde{\mathsf{\tau }}}_{hp}$ are parameters.

$${\overline{m}}_p(t)=\frac{1}{1+\exp \left(\tilde{B}(t)\right)},$$

(5)

where $\tilde{B}(t)=\left(B(t)-{\tilde{\mathsf{\theta }}}_{mp}\right)/2{\tilde{\mathsf{\sigma }}}_{mp}$, ${\tilde{\mathsf{\theta }}}_{mp}$ and ${\tilde{\mathsf{\sigma }}}_{mp}$ are parameters.

The output signal of rVRG is modeled using the following expressions:

$$R_p(t)=\{\begin{array}{l}K(t),\mathrm{if}B(t)>T{r}_1\mathrm{and}{R}_p(t-\mathsf{\Delta }t)<T{r}_4,\\ I_p\left(t\right),\mathrm{if}B(t)>T{r}_1,\\ 0,\mathrm{if}B(t)\le T{r}_2,\\ R_p(t-\mathsf{\Delta }t),\mathrm{if}{R}_p(t-\mathsf{\Delta }t)>T{r}_1\mathrm{and}\left|\frac{B(t)-B(t-\mathsf{\Delta }t)}{\mathsf{\Delta }t}\right|<\epsilon ,\end{array}$$

(6)

where

$$I_p(t)=I_L\left(\left(B\left(t-\mathsf{\Delta }t\right)+B\left(t\right)\right)\frac{\mathsf{\Delta }t}{2}+R_p\left(t-\mathsf{\Delta }t\right)\right)$$

(7)

is the integral, which is calculated using the trapezoidal rule of integration, $\mathsf{\epsilon }$ is a small positive value and $K(t)$ is the variable affecting the amplitude of respiration [42]; the equation for $K(t)$ is shown in Section 2.1.3.

2.1.2. Mathematical Model of the Lungs

The diaphragm is modeled using a linear equation for a spring, which is compressed by the external force $R_p(t)$ representing the output of the rVRG complex:

$${\dot{x}}_m(t)=-k_1^x{x}_m(t)+k_2^x{R}_p(t).$$

(8)

The pleural pressure $P_L(t)$ is described by the following expression:

$$P_L(t)=P_m-P_{L0}-k_p{x}_m(t),$$

(9)

where $P_m$ is the mouth pressure, $P_{L0}$ is the residual air pressure in the lungs after expiration and $k_p$ is the transfer coefficient. Gas pressure in the alveoli $P_A(t)$ is described by the following expression:

$${\dot{P}}_A(t)=\frac{P_{m}E}{P_A(t)}{Q}_A(t)+{\dot{P}}_L(t),$$

(10)

where $E$ is the lung elasticity, and $Q_A(t)$ is the total gas flow into the alveoli:

$$Q_A(t)=q(t)+D_c\left(p_c(t)-p_{ac}(t)\right)+D_o\left(p_o(t)-p_{ao}(t)\right),$$

(11)

where $q(t)=\left(P_m-P_A(t)\right)/R_a$ is the gas flow in the mouth ($q(t)>0$ corresponds to inspiration, and $q(t)<0$ corresponds to expiration), $R_a$ is the resistance of the airways, $D_c$ and $D_o$ are the diffusion capacities, $p_c(t)$ and $p_o(t)$ are the partial pressures of CO₂ and O₂ in the blood, respectively, and $p_{ac}(t)$ and $p_{ac}(t)$ are the partial pressures of the same gases in the alveoli:

$$p_{ac}(t)=f_c\left(P_A\left(t\right)-P_w\right),$$

(12)

$$p_{oc}(t)=f_o\left(P_A\left(t\right)-P_w\right),$$

(13)

where $P_w$ is the pressure of water vapor at 37 °C. The partial pressures $p_c(t)$ and $p_o(t)$ were set to values $p_c^0$ and $p_o^0$, respectively, (see Table 1) after each heart contraction, and in between the contractions, they were modeled using the following equations:

$${\dot{p}}_c(t)=\frac{D_c}{{\mathsf{\sigma }}_c{V}_c{C}_u}\left(p_{ac}\left(t\right)-p_c\left(t\right)\right)+\delta \left(\frac{l_{2}z\left(t\right)H}{{\mathsf{\sigma }}_c}-r_2{p}_c\left(t\right)\right),$$

(14)

$${\dot{p}}_o(t)=\frac{D_o}{\mathsf{\sigma }{V}_c{C}_u}\left(1+4\frac{T_h}{\mathsf{\sigma }}\frac{dF}{dp}|{}_{p_0}\right)\left(p_{ao}\left(t\right)-p_o\left(t\right)\right),$$

(15)

where ${\mathsf{\sigma }}_c$ and $\mathsf{\sigma }$ are the solubilities of CO₂ and O₂, respectively, in blood plasma, $V_c$ is the total volume of the lung capillaries, $C_u$ is the molar volume of the gas at 37 °C and 760 mm Hg and $T_h$ is the concentration of hemoglobin. The blood concentration of bicarbonate $z(t)$ was set to the fixed value $z(t)=z_0=\left(46{\mathsf{\sigma }}_c{r}_2\right)/\left(l_{2}H\right)$ after each heart contraction, and in between the contractions, it was modeled using the following equation:

$$\dot{z}(t)=\delta \left(r_2{\mathsf{\sigma }}_c{p}_c\left(t\right)-l_{2}z(t)H\right),$$

(16)

where $r_2$ and $l_2$ are the speed of dehydration and hydration reactions, respectively, $H$ is the concentration of $H^{+}$ and $\mathsf{\delta }$ is the speed of the chemical reaction of bicarbonate (HCO₃⁻), H₂O and CO₂ [42].

$$F(p)=\frac{L{K}_T\mathsf{\sigma }p{\left(1+K_T\mathsf{\sigma }p\right)}^3+K_R\mathsf{\sigma }p{\left(1+K_R\mathsf{\sigma }p\right)}^3}{L{\left(1+K_T\mathsf{\sigma }p\right)}^4+{\left(1+K_R\mathsf{\sigma }p\right)}^4},$$

(17)

where $L$, $K_T$ and $K_R$ are parameters.

The dynamics of the relative concentration of CO₂ and O₂ are described as follows:

$$\begin{array}{l}{\dot{f}}_o=[D_o\left(p_o\left(t\right)-p_{ao}\left(t\right)\right)+q\left(t\right)\left(f_{oi}\left(t\right)-f_o\left(t\right)\right)-\\ -f_o\left(t\right)\left(D_o\left(p_o\left(t\right)-p_{ao}\left(t\right)\right)+D_c\left(p_c\left(t\right)-p_{ac}\left(t\right)\right)\right)]/V_A\left(t\right),\end{array}$$

(18)

$$\begin{array}{l}{\dot{f}}_c=[D_c\left(p_c\left(t\right)-p_{ac}\left(t\right)\right)+q\left(t\right)\left(f_{ci}\left(t\right)-f_c\left(t\right)\right)-\\ -f_c\left(t\right)\left(D_o\left(p_o\left(t\right)-p_{ao}\left(t\right)\right)+D_c\left(p_c\left(t\right)-p_{ac}\left(t\right)\right)\right)]/V_A\left(t\right),\end{array}$$

(19)

where $f_{oi}$ and $f_{ci}$ are the relative concentrations of CO₂ and O₂, respectively, in the inspired air:

$$f_{oi}\left(t\right)=\{\begin{array}{l}{f}_{om},q\left(t\right)>0\\ f_o\left(t\right),q\left(t\right)\le 0\end{array},f_{ci}\left(t\right)=\{\begin{array}{l}{f}_{cm},q\left(t\right)>0\\ f_c\left(t\right),q\left(t\right)\le 0\end{array},$$

(20)

where $f_{om}$ and $f_{cm}$ are the relative concentrations of CO₂ and O₂, respectively, in the mouth. The lung volume $V_A$ was calculated as follows:

$$V_A\left(t\right)=\left(P_A\left(t\right)-P_L\left(t\right)\right)/E.$$

(21)

2.1.3. Control of Respiration

When modeling the control of respiration, we preserved the general structure of the control system but replaced the linear models of the chemoreceptors, which were used in [42,44], with the nonlinear sigmoidal transfer functions [43]. Following [37,38], we refer to the chemosensitive area of the respiratory center as the central chemoreceptors. Since these chemoreceptors are mainly stimulated by the changes in the CO₂ concentration in the arterial blood, we refer to them as the central CO₂ chemoreceptors. The activity of the central CO₂ receptors in the medulla is modeled as follows:

$$v_c(t)=k_1^c\mathrm{th}\left(\left(p_{ce}\left(t-{\mathsf{\theta }}_{vc}\right)-k_2^c\right)/k_3^c\right)+k_4^c,$$

(22)

where $k_1^c$, $k_2^c$, $k_3^c$ and $k_4^c$ are parameters, and ${\mathsf{\theta }}_{vc}$ is the time delay due to the finite speed with which blood flows from the aorta to the brain.

The activity of carotid O₂ receptors is modeled as follows:

$$\{\begin{array}{l}{v}_o(t)=k_1^o\mathrm{th}\left(\left(p_{co}\left(t\right)-k_2^o\right)/k_3^o\right)+k_4^o,p_{ce}\left(t\right)\le 100,\\ v_o(t)=k_5^o{p}_{oe}\left(t\right)+k_4^o,p_{ce}\left(t\right)>100.\end{array})$$

(23)

Then, following the results of Ben-Tal et al. [40], the activity of the receptors controls the frequency and depth of respiration:

$${\tilde{g}}_t(t)=g_0+k_c^g{v}_c(t)+k_o^g{v}_o(t),$$

(24)

$$K(t)=K_0+k_c^K{v}_c(t)+k_o^K{v}_o(t)+I_k(t),$$

(25)

where $I_k(t)$ is the integral sum of chemoreceptor activity:

$$I_k(t)=k_1^{I_k}\mathsf{\Delta }t\left(\frac{v_c(t)+v_c(t-\mathsf{\Delta }t)}{2}\right)+k_2^{I_k}\mathsf{\Delta }t\left(\frac{v_o(t)+v_o(t-\mathsf{\Delta }t)}{2}\right).$$

(26)

In this study, we used the proposed model to simulate the cardiovascular dynamics in healthy subjects and subjects with blocked control of the respiratory frequency and amplitude and, therefore, disturbed cardiorespiratory coupling. When simulating the blockade, we set the parameters $k_c^g$, $k_o^g$, $k_c^K$, $k_o^K$, $k_1^{I_k}$ and $k_2^{I_k}$ to zero.

2.1.4. Model of the Cardiovascular System

The model of the cardiovascular system is based on the model proposed in our earlier study [20]. Here, we closely follow this description, except for the aspects of the model which were changed to integrate with the model of the respiratory system.

We used an “integrate and fire” model to simulate the beating heart:

$$\dot{\mathsf{\phi }}\left(t\right)=\frac{1}{T_0}{f}_s\left(t\right)f_p\left(t\right),$$

(27)

where $\mathsf{\phi }\left(t\right)$ is the phase of the sinus node, and $T_0$ is the heart period, which is modulated by the sympathetic and parasympathetic factors $f_s\left(t\right)$ and $f_s\left(t\right)$, respectively. Without the autonomic control, the sinus node generates periodic saw-like impulses with a constant period of $T_0$ seconds.

The arterial pressure is modeled using the following two expressions. Expression (28) describes the arterial pressure during the first T_sys seconds (0.125 s) after the initiation of the heart cycle, and during that time, the arterial pressure rapidly grows:

$$p_{sys}\left(t\right)=D_{i-1}+S\left(t\right)\frac{\left(t-T_{i-1}\right)}{T_{sys}}\exp \left(\frac{\left(t-T_{i-1}\right)}{T_{sys}}\right)+k_p^B{V}_A\left(t\right),$$

(28)

where $D_{i-1}$ is the systolic pressure at the end of the previous cardiac cycle, $T_{i-1}$ is the moment in time when the previous cardiac has ended, T_sys is a fixed duration of this systolic phase and S(t) is the cardiac contractility [40] defined as

$$S\left(t\right)=S^{\prime }\left(t\right)+\left(\widehat{S}-S^{\prime }\left(t\right)\right)\frac{{S^{\prime }}^{n_c}\left(t\right)}{{\widehat{S}}^{n_c}+{S^{\prime }}^{n_c}\left(t\right)},$$

(29)

where $S^{\prime }\left(t\right)=S_0+k_S^c{c}_c\left(t\right)+k_S^t{L}_{i-1}$. Here, $S_0$ is the heart contractility without the autonomic control, $\widehat{S}$, $k_S^c$ and $k_S^t$ are parameters, $c_c\left(t\right)$ is the concentration of noradrenaline in the blood that circulates in the heart and $L_{i-1}$ is the duration of the previous cardiac cycle.

After T_sys seconds, the systolic phase of the cardiac cycle ends and the arterial pressure begins to slowly decrease following the Windkessel model [45]:

$${\dot{p}}_{dia}\left(t\right)=-\frac{p_{dia}\left(t\right)}{{\mathsf{\tau }}_v\left(t\right)},$$

(30)

where ${\mathsf{\tau }}_v\left(t\right)$ is the function of the vascular concentration of noradrenaline ($c_v\left(t\right)$):

$${\mathsf{\tau }}_v(t)={\mathsf{\tau }}_v^{(0)}-{\overline{\mathsf{\tau }}}_v\left[c_v(t)+\left({\widehat{c}}_v-c_v(t)\right)\frac{c_v{(t)}^{n_v}}{{\widehat{c}}_v{}^{n_v}{c}_v{(t)}^{n_v}}\right],$$

(31)

where ${\mathsf{\tau }}_v^{(0)}$, ${\overline{\mathsf{\tau }}}_v$, ${\widehat{c}}_v$ and $n_v$ are parameters.

The absolute value and rate of change in arterial pressure [40] are sensed by the baroreceptor located in the carotid sinus nodes. To incorporate the nonlinear properties of the baroreceptors, reported in [46], we also introduced the following sigmoidal nonlinearity. To preserve the stationarity of the model, we also excluded the Brownian noise introduced into this equation by Kotani et al.:

$$\{\begin{array}{l}{v}_b\left(t\right)=k_1^b\mathrm{th}\left(p\left(t\right)+k_2^b\dot{p}\left(t\right)-k_3^b\right)/k_4^b+k_5^b,p(t)\ge p_0\\ v_b\left(t\right)=0,p(t)<p_0=50\end{array},$$

(32)

where $k_1^b$, $k_2^b$, $k_3^b$, $k_4^b$ and $k_5^b$ are parameters of the nonlinearity, fitted to the data in [46], and $p_0$ is the lowest pressure baroreceptors react to [47].

The higher nervous centers process the inputs of the baroreceptors and adjust the activity of the heart period autonomic control and the vessel tone autonomic control. Seidel et al. [40] introduced an expression to calculate the efferent sympathetic nervous activity $v_s\left(t\right)$ as an averaged firing rate of this population of neurons, similar to how the neurons of the pre-Bötc complex were modeled in [42]. We mostly followed this approach but separated the α-sympathetic activity $v_{\mathsf{\alpha }s}\left(t\right)$, mostly affecting the vessel tone, and the β-sympathetic activity $v_{\mathsf{\beta }s}\left(t\right)$, mostly affecting the heart period [37,38,46]:

$$v_{\mathsf{\alpha }s}\left(t\right)=v_{\mathsf{\alpha }s}^0-k_{\mathsf{\alpha }s}^b{v}_b(t)+k_{\mathsf{\alpha }s}^c{v}_c(t)-k_{\mathsf{\alpha }s}^r{V}_A\left(t\right),$$

(33)

$$v_{\mathsf{\beta }s}\left(t\right)=v_{\mathsf{\beta }s}^0-k_{\mathsf{\beta }s}^b{v}_b(t)+k_{\mathsf{\beta }s}^c{v}_c(t)+k_{\mathsf{\beta }s}^r{V}_A\left(t\right),$$

(34)

where $v_{\mathsf{\alpha }s}^0$ and $v_{\mathsf{\beta }s}^0$ describe the activity of the sympathetic control under resting conditions, and $k_{\mathsf{\alpha }s}^b$, $k_{\mathsf{\alpha }s}^c$, $k_{\mathsf{\alpha }s}^r$, $k_{\mathsf{\beta }s}^b$, $k_{\mathsf{\beta }s}^c$ and $k_{\mathsf{\beta }s}^r$ are the coefficients that reflect the influence of baro- and chemoreceptor activity. Following the models from [37,38], we simulated the excitatory influence of the central CO₂ chemoreceptors on both the sympathetic and parasympathetic branches of the autonomic control of the circulation. To simplify the model, we did not simulate the effect of carotid O₂ receptors on the autonomic control. According to the experimental data, this influence is much weaker than the influence from the central CO₂ chemoreceptors when the blood concentration of O₂ is normal or above the norm. Activation of the parasympathetic autonomic control is described as [42]

$$\{\begin{array}{l}{v}_p\left(t\right)=v_p{}^{\prime }\left(t\right),v_p{}^{\prime }\left(t\right)>0\\ v_p\left(t\right)=0,v_p{}^{\prime }\left(t\right)\le 0\end{array},$$

(35)

where $v_p{}^{\prime }\left(t\right)=v_p^0+k_p^b{v}_b(t)+k_p^c{v}_c(t)-k_p^r{V}_A(t)+k_p^{\xi }\xi (t)$, $v_p^0$ is the parasympathetic activity under resting conditions, $k_p^b$, $k_p^c$ and $k_p^r$ are the coefficients that reflect the influence of the higher nervous centers on the parasympathetic control of the heart period and $\xi$ is red noise of a “central” origin [2,25], which is related to the influence of the higher nervous centers on the parasypathetic branch of the autonomic control. To preserve the simplicity of the model, we introduced the noise only to the parasympathetic branch of the autonomic control, but the attribution of the noise of a central origin to the parasympathetic loop only is an incomplete interpretation [48]. To generate the noise, we used the first-order AR model: ${\xi }_t=a{\xi }_{t-1}+{\xi }_n,$ where a = 0.97, and ${\xi }_n$ is the normal noise with zero mean and dispersion of unity. Changes in the sympathetic activation affect the cardiac and vascular concentrations of noradrenaline:

$${\dot{c}}_c\left(t\right)=-\frac{c_c\left(t\right)}{{\mathsf{\tau }}_c}+k_c^S{v}_{\mathsf{\beta }s}\left(t-{\mathsf{\theta }}_c\right),$$

(36)

$${\dot{c}}_v\left(t\right)=-\frac{c_v\left(t\right)}{{\mathsf{\tau }}_v}+k_v^S{v}_{\mathsf{\alpha }s}\left(t-{\mathsf{\theta }}_v\right),$$

(37)

where ${\mathsf{\tau }}_c$ and ${\mathsf{\tau }}_v$ are the time constants, and ${\mathsf{\theta }}_c$ and ${\mathsf{\theta }}_v$ are the time delays caused by the finite speed of the neural transition and the finite speed of noradrenaline turnover. Changes in the concentration of noradrenaline and the activity of the parasympathetic control affect the heart period:

$$f_s\left(t\right)=1+k_{\mathsf{\phi }}^c\left(c_c\left(t\right)+\left({\widehat{c}}_c-c_c\left(t\right)\right)\frac{c_c^{n_s}\left(t\right)}{{\widehat{c}}_c^{n_s}+c_c^{n_s}\left(t\right)}\right),$$

(38)

$$f_p\left(t\right)=1-k_{\mathsf{\phi }}^p\left(v_p\left(t-{\mathsf{\theta }}_p\right)+\left({\widehat{v}}_p-v_p\left(t-{\mathsf{\theta }}_p\right)\right)\frac{v_p^{n_p}\left(t-{\mathsf{\theta }}_p\right)}{{\widehat{v}}_p^{n_p}+v_p^{n_p}\left(t-{\mathsf{\theta }}_p\right)}\right)F\left(\mathsf{\phi }\left(t\right)\right),$$

(39)

where ${\mathsf{\theta }}_p$ is the delay time. No separate equation is introduced to model the changes in the concentration of acetylcholine (parasympathetic transmitter) because the rate of its turnover is faster than the dynamics on which the model is focused. The phase effectiveness curve represents the changes in the sinus node sensitivity to the parasympathetic control throughout the cardiac cycle:

$$F\left(\mathsf{\phi }\right)={\mathsf{\phi }}^{1.3}\left(\mathsf{\phi }-0.45\right)\frac{{\left(1-\mathsf{\phi }\right)}^3}{0.008+{\left(1-\mathsf{\phi }\right)}^3}.$$

(40)

To block the regulation of the respiratory frequency and amplitude, the parameters $k_o^g$, $k_c^g$, $k_c^K$, $k_o^K$, $k_1^{I_k}$ and $k_2^{I_k}$ were set to zero, while the other parameters were left unchanged.

2.2. Complexity Indices

To measure the complexity of RR intervals in the model, we used the correlation dimension [49] and the largest Lyapunov exponent [50]. In our earlier study [20], we already used these measures to estimate the complexity of the model and experimental RR intervals. Here, we closely follow the description of the measures provided in [20].

To reconstruct the phase space, we used the method of delays as recommended in [50]. The correlation dimension d was estimated from the correlation integral C(l) [49]:

$$C\left(l\right)=\underset{n\to \infty }{\lim }\frac{n\left(l\right)}{N^2},$$

(41)

where n(l) is the number of points in the reconstructed attractor, for which the Euclidean distance to the nearest neighbor is smaller than l. The values of l were varied from 0.1 to 0.3 of the standard deviation of the RR intervals. N = 4000 is the number of points used for calculation. For dynamical systems, C(l) ~ l^d and d can be calculated as

$$d=\frac{\ln \left(C\left(l\right)\right)}{\ln \left(l\right)}.$$

(42)

To estimate the largest Lyapunov exponent, we used the Rosenstein algorithm [50] since it can be applied to short time series. In the first step of the Rosenstein algorithm, we found the nearest neighbor of each point on the attractor, and then neighbors close in time were excluded [50]. For the dynamical system, the distance between nearest neighbors increases as follows:

$$\ln \left(L\right)\approx \ln \left(L_0\right)+{\lambda }_{0}t,$$

(43)

where $L_0$ is the initial distance, λ₀ is the largest Lyapunov exponent and t is the time. Then, λ₀ can be calculated as

$${\lambda }_0=\frac{\langle \ln \left(L\right)\rangle }{t}.$$

(44)

Parameters of the methods were estimated in our previous study [20]. They are listed in

Table 2. Parameters of the methods for calculating the correlation dimension and the largest Lyapunov exponent.

Correlation Dimension		The Largest Lyapunov Exponent
D	13	D	13
τ (s)	0.04	τ (s)	1
l	0.1–0.3	t (s)	0.6
Window length (s)	1000	Window length (s)	1000

, where D is the embedding dimension, τ is the delay time used for the reconstruction of the embedding space and l is the maximal Euclidean distance between the nearest neighbors in relation to the standard deviation of the sequence of RR intervals.

Since the model is limited to the simulation of processes with frequencies of 0.04–0.4 Hz, as this is the time scale of the autonomic control of the circulation, the RR intervals were filtered using a 0.04–0.4 Hz bandpass filter.

3. Results

In this study, we used the proposed model to simulate the cardiovascular dynamics in healthy subjects and subjects with blocked regulation of the respiratory amplitude and frequency and, therefore, disturbed cardiorespiratory coupling.

Figure 2, Figure 3 and Figure 4 show the time series and power spectra of the model RR intervals, arterial pressure signals and lung volume signals for a healthy subject and a subject with blocked regulation of the respiratory amplitude and frequency.

Table 3. Spectral and statistical indices and complexity measures estimated from the model data.

Indices	Healthy Subject	Control of Respiratory Amplitude and Frequency Is Absent	Significance Level
Minute volume of ventilation, l	5.09 ± 0.002	12.7 ± 0.0001	p < 0.05
Blood concentration of CO2, mmHg	40.8 ± 0.0005	34.7 ± 0.0006	p < 0.05
Blood concentration of O2, mmHg	100 ± 0.01	126 ± 0.003	p < 0.05
Depth of respiration, l	0.17 ± 0.0001	0.40 ± 0.0001	p < 0.05
Heart frequency, Hz	1.05 ± 0.0003	1.07 ± 0.0006	p < 0.05
SAP, mmHg	138 ± 0.02	136 ± 0.01	p < 0.05
DAP, mmHg	94 ± 0.02	94 ± 0.007	p < 0.05
LF, ms2	954 ± 33	118 ± 4	p < 0.05
HF, ms2	725 ± 20	630 ± 6	p < 0.05
λ0	0.07 ± 0.002	0.07 ± 0.003	p = 0.10
d	4.70 ± 0.27	4.89 ± 0.22	p = 0.27

SAP is the systolic arterial pressure in the aorta, DAP is the diastolic arterial pressure in the aorta, ${\lambda }_0$ is the largest Lyapunov exponent, d is the correlation dimension, LF is the power in the 0.04–0.15 Hz range of RR intervals and HF is the power in the 0.15–0.4 Hz range of RR intervals. The values of the indices are listed as the mean ± standard error of the mean. We used the Mann–Whitney test to compare the data.

lists the spectral indices, statistical indices and complexity measures estimated from these signals.

The power spectra and statistical and spectral indices agreed well with the experimental data in the 0.04–0.4 Hz range, which is related to the autonomic control of the circulation [46,51,52]. The model did not simulate the very low frequency rhythms (0.003–0.04 Hz) related to the humoral control of the circulation, which Kotani [40] modeled by adding the Brownian noise to Expression (32). We excluded the noise because it led to the nonstationarity of the model at the time scale of hundreds of seconds. Additionally, the simulation of the humoral control is beyond the scope of the model.

The power spectra in Figure 2b and Figure 3b show that the blockade of the regulation of the respiratory frequency and amplitude led to more narrow respiration-related peaks at approximately 0.3 Hz. A similar change can be seen in the power spectrum of the lung volume time series (Figure 4b). The respiration became very regular because, after the blockade of the regulation of the respiratory amplitude and frequency, the RR intervals were no longer affecting the frequency or amplitude of the respiratory signal.

Another qualitative change caused by the blockade was the decrease in the power in the LF frequency range of RR intervals (Figure 2b), which is also evident from the decrease in the LF index (Table 3). This effect was caused by an increase in the depth and frequency of respiration (Figure 4a), which led to a higher than normal concentration of O₂ in the arterial blood (125.80 ± 0.003 mm Hg versus 100 mm Hg) and a lower than normal concentration of CO₂ (34.65 ± 0.0006 mm Hg versus 40 mm Hg) (Table 3). This inhibited the activity of the central CO₂ chemoreceptors and, therefore, the activity of the sympathetic and parasympathetic control of the heart period and arterial pressure, because the activity of these systems is directly proportional to the activity of both the carotid baroreceptors and the central CO₂ chemoreceptors [38] (Equations (33) and (34)). From a physiological standpoint, we interpret these results as follows. The blockade of the regulation of the respiratory amplitude and frequency inhibited the cardiovascular system’s ability to adapt to the physiological changes. As a result, the frequency and amplitude of the respiration exceeded the required values, which decreased the concentration of CO₂ in the arterial blood below the normal values and increased the concentration of O₂. These changes led to the deactivation of the chemoreceptors, which, in turn, decreased their excitatory influence on both the sympathetic and parasympathetic control of the circulation.

A similar effect inhibited the activity of the parasympathetic control of the heart period (Equation (35)), and, since noise of a “central” origin affects the autonomic control through the parasympathetic control loop [41], the influence of noise on the circulation was also inhibited. This change is evident from the power spectra of the model signals (Figure 2b, Figure 3b and Figure 4b), as the baseline power is lower in the case of the blocked regulation of the respiratory frequency and amplitude.

The blockade of the regulation of the respiratory frequency and amplitude also slowed down the heart and decreased the heart rate variability and blood pressure in the aorta. However, seemingly paradoxically, the estimations of the correlation dimension and the largest Lyapunov exponent did not change (Table 3). The results shown by our previous model [20] may suggest an explanation for this intriguing effect. In [20], we showed that the blockade of the autonomic control of the circulation increased the λ₀ and d to higher than normal values. At the same time, when the dispersion of the noises was set to zero, both λ₀ and d became lower than normal values. In the active test with the blockade of the regulation of the respiratory frequency and amplitude, both the activity of the autonomic control and the power of the “central” noise decreased, compensating for each other’s influence on the complexity measures.

4. Discussion

Investigation of mathematical models of biological systems has several advantages over the analysis of experimental data. The models are stationary and can simulate active tests which are dangerous for the health of a real subject. However, it also has limitations. Every model of a complex biological system is a compromise between the complexity of the model and its ability to reproduce the dynamics of the system. Simpler models cannot simulate and explain certain features of the dynamics of the studied systems, but overcomplicated models cannot be properly calibrated and produce unreliable results.

The proposed model focuses on the nonlinear properties of the autonomic control of the circulation, because these nonlinear properties are necessary to explain several important features of cardiovascular coupling [20,53,54]. For the same reason, we adapted the nonlinear control loops for the respiratory frequency and amplitude.

At the same time, our earlier results show that detailed modeling of the heart and the blood hydrodynamics is not essential to simulate the important features of cardiovascular coupling [20,53,54]. Based on these results, we kept simplified models for these elements.

The main limitation of the model of the respiratory system is its inability to simulate deep and fast breathing, which requires modeling of additional sections of the respiratory center and additional respiratory muscles. This simplification limits the model’s ability to simulate breathing under stressful conditions, during exercises or in diseases such as sleep apnea. However, simulation of these effects is not required to investigate cardiorespiratory coupling.

Despite the adopted simplifications, the proposed model still has a high number of free parameters. This complicates the process of their calibration. However, as shown in Section 3, the model still quantitatively reproduces the average characteristics of a healthy subject.

Earlier, we used a mathematical model to investigate the contributions of the deterministic chaotic dynamics of the autonomic control, noises of a central origin and the irregularity of respiration to the complexity of the heart period [20]. The model we investigated in [20] was proposed from the “first principles”, and we consider it an adequate representation of the cardiovascular system because it previously simulated several important effects, including a passive orthostatic test [53], arterial hypertension [54], blockade of the autonomic control [54] and cardiorespiratory synchronization [54]. However, this model has one important limitation: it is a simplified model of respiration. We modeled the respiration using a sinusoidal signal, and the model supported two types of respiration: with a constant frequency, or with a randomized, normally distributed frequency of respiration. The model did not account for the regulation of respiration, and because of that, the cardiorespiratory coupling was only one-directional.

In the model with a simplified simulation of respiration, we observed several interesting effects: the complexity of the heart period persisted even after the dispersion of the “central” noises was set to zero; blockade of the autonomic control led to higher than normal values of complexity measures; the values of complexity measures became lower than normal when the dispersion of the “central” noises was set to zero. Finally, we did not observe any changes in the complexity measures in the presence of switching from respiration with a randomized frequency to respiration with a constant frequency.

In the present study, we aimed to investigate the contribution of nonlinear bidirectional coupling to the complexity of the heart period and used a more detailed model of the respiratory system from [42]. The modified model accounted for the control of respiration, which is based on the heart period and concentrations of CO₂ and O₂ in the arterial blood. However, even in this modified model, we did not observe any changes in the complexity measures after the blockade of the regulation of the respiratory frequency and amplitude. At the same time, we detected changes in spectral and statistical characteristics of the model signals. We think that the results we obtained in this study and our previous study contribute to the discussion about the advantages of complexity measures over the classical statistical and spectral indices. We present a summary of this debate in the following two paragraphs, and the role of our results is discussed after that.

Both classical indices and complexity measures have their advantages and disadvantages. It is mostly accepted that the LF and HF indices are related to the modulation of the sympathetic and parasympathetic branches of the autonomic control, and some active tests, such as the passive orthostatic test [52], support this hypothesis. However, spectral and statistical indices do not linearly correlate with other measurements of autonomic activity, such as sympathetic nervous outflow [55] or norepinephrine spillover [56]. This fact suggests that the relationships between the RR intervals and the activity of the autonomic control of the circulation are nonlinear, and that nonlinear measures of complexity applied to the analysis of RR intervals can be a better marker of overall cardiovascular health and its adaptability to external stress and internal challenges [57,58].

However, several researchers advise being cautious when interpreting the complexity indices estimated from experimental data of a biological origin, because at least some elements of the underlying system are not deterministic, but stochastic. Some authors suggest that stochastic elements of the cardiovascular system cause inconsistency in complexity measures [31]. For example, the experimental studies [58,59] showed no clear correlation between the blockade of the sympathetic control and complexity measures. Our results suggest that the inconsistency in the complexity measures estimated from experimental data can be caused by different regulatory responses taking place simultaneously, and that under such conditions, the complexity measures should be used in conjunction with the statistical and spectral indices that provide different important information about the dynamics of the autonomic control of the cardiovascular system. Parallel investigation of mathematical models can also help with the interpretation of the results.

5. Conclusions

We investigated the contribution of bidirectional cardiorespiratory coupling to the complexity of the heart period, using a mathematical model of the cardiovascular and respiratory systems to simulate an active experiment, which involved the blockade of the regulation of the respiratory frequency and amplitude and disturbance of the bidirectional cardiorespiratory coupling. Well-known spectral indices, statistical indices, the largest Lyapunov exponent and correlation dimensions were estimated to measure the changes that occurred in the dynamics of the model and in the complexity of the heart period.

Blockade of the regulation of the respiratory frequency and amplitude slowed down the heart, decreased the heart rate variability and aortic arterial pressure, inhibited both the sympathetic and parasympathetic control of the circulation and inhibited the effects of “central” noise on the circulation. However, we detected no changes in the correlation dimension or the largest Lyapunov exponent.

Our earlier model study showed that the values of complexity measures increase after the blockade of the autonomic control of the circulation. On the contrary, the values of complexity measures decrease when the power of central noises decreases. These effects took place simultaneously after the blockage of the autonomic control of respiration and compensate for each other’s influence on the complexity of the heart period.

The obtained results indicate that complexity measures should be used together with classical statistical and spectral indices when analyzing the dynamics of the cardiovascular and respiratory systems under different physiological conditions, and that mathematical modeling helps to interpret the results.