A database of ten healthy male volunteers (aged 54.0 ± 13.5 years, weight 75.6 ± 10.3 kg, and height 169.9 ± 7.0 cm), non-smoking subjects, normotensive, normal lung function, and with no history of lower limb or cardiopulmonary disorders under a CPET in an electromagnetically-braked cyclorgometer (CardiO2; MedGraphics Corp., St. Paul, MN), was used in this study. The experimental protocol was carried out by a trained medical staff of the Pulmonary Function Laboratory of the Hospital Clínico de Barcelona for evaluation of exercise tolerance. This study, which was performed following the Helsinki declaration regarding the investigation with human subjects, had been previously approved by the Committee on Investigations Involving Human Subjects at the Hospital Clinic, University of Barcelona, Barcelona, Spain. Informed consent was obtained from all individual participants included in the study.

After a warming-up session consisting of five min of stretching and three min of unloaded pedaling on the cycloergometer, the exercise workload was increased by 5 or 10 W/min every minute until the subject stopped due to symptoms (i.e., dyspnea and/or leg fatigue) or no longer maintained the constant pedal rate required. The following signals were registered by the cycloergometer every 15 s: exercise workload (W), minute ventilation (V˙E), tidal volume (V_T), inspiration time (T_I), expiration time (T_E), O₂ consumptions (V˙O2), CO₂ productions (V˙CO2), expired fraction of O₂ (P_etO2) and CO₂ (P_etCO2), and heart rate (HR). Taking into account that no subject had gas exchange impartments, arterial blood pressures of CO₂ and O₂ were adjusted considering average normal values at rest, 39.156 mmHg and 104.37 mmHg respectively (Batzel et al., 2007), as follows:

Figure 1 describes registered data in median and interquartile distance of the database analyzed. It can see that, excepting T_I, P_aCO2 and P_aO2 whose values remained almost constant, all variables increased with the exercise workload. This latter was related to increments of both V˙CO2 and V˙O2 that, in turn, had a linear relationship with V˙E.

RS2 comprises the cardiorespiratory model described in Cheng et al. (2010) and Cheng and Khoo (2012) and it is referred by its authors as “PNEUMA.” This model is the result of the integration of key published models of the respiratory and cardiovascular system. It has been designed to simulate the cardiorespiratory control dynamic during wakefulness and sleep, so that provides realistic predictions of the physiological responses under a wide variety of conditions such as the day-to-day sleep-wake cycle, Cheyne-Stokes respiration in chronic heart failure, obstructive sleep apnea (OSA) and hypoxia-induced periodic breathing. It can be used to investigate several types of interventions: isocapnic and hypercapnic and/or hypoxemic gas administration, the Valsalva and Mueller maneuvers, and the application of continuous positive airway pressure (CPAP). The most recent version incorporates a sub-model of glucose-insulin-fatty acid regulation to simulate also the metabolic control of glucose-insulin dynamics and its interaction with the autonomic control in obese individuals (Cheng and Khoo, 2012). RS2 is available at the USC Biomedical Simulation http://bmsr.usc.edu/software/pneuma/.

RS2 has been developed using a hierarchical structure in such a way that the degree of complexity associated with each level of organization is adapted appropriately to the investigation of physiological processes at each level. This feature allows that the whole model can be presented in a compact and efficient way. RS2 is mainly composed of five principal interconnected compartments: the respiratory system, the cardiovascular system, the central control system, the sleep mechanism, and the metabolic control system. The last one allows simulating the metabolic control over the glucose-insulin dynamics and its interaction with the autonomic control in obese individuals (Cheng and Khoo, 2012). A schematic diagram of the model is given in the Supplementary Material Section (see Figure S2).

In this model, the respiratory subsystem allows the simulation of both gas exchange system and ventilatory mechanics. Unlike RS1, RS2 provides a more detailed description of chemical and physical processes generated during respiration. It includes the progressive decline of inspiratory gas pressure (P_ICO2and P_IO2), due to different sectors of anatomic dead space, using a first-order dynamic approach to calculate the arterial partial pressures of CO₂ and O₂ in the areas closer to the alveoli and five small serial compartments (Khoo, 1990). It also considers the phenomena of convection and dissociation of respiratory gases during cardiovascular mixing by using a second-order dynamic system that relates the arterial pressures to the alveolar pressures (Spencer et al., 1979). In general, RS2 presents more elaborate expressions to describe in greater detail the different processes that comprise respiration, as well as its control and interaction with the cardiovascular system (Cheng et al., 2010). Only the variables directly related to ventilation control are considered here for brevity, as shown below.

The ventilatory controller in RS2 incorporates the contribution of the central (D_c) and peripheral (D_p) chemoreceptors. In this model, the central chemoreceptors only respond to variations in the brain partial pressure of CO₂ (P_bCO2) while the peripheral chemoreceptors are influenced by the arterial partial pressures of CO₂ (P_aCO2) and the oxygen saturation in arterial blood (SAO₂) and their multiplicative interaction (Khoo, 1990). During wakefulness, the total ventilatory demand (D_T) is defined by the sum of the central and peripheral chemoreceptors responses as follows:

where, I_c, I_pCO2, and I_pO2 represent the central and peripheral chemoreceptors activation threshold and, they are equal to 45, 38, and 102.4, respectively. In this case, P_bCO2 is controlled by the metabolic rate (MR_bCO2) and the brain blood flow (Q_B) and it is defined as a function of P_aCO2 (Read and Leigh, 1967).

To adjust the breathing pattern, RS2 computes the respiratory frequency in function of ventilatory demand by using the following expressions (Duffin et al., 2000):

where, F_b is the basal frequency and T_D and T_p are thresholds of D_T that determine the behavior of respiratory rate. For the last two options of Equation (9) the respiratory rate varies linearly respect to D_T with a slope established by the scaling factors S_1F and S_2F, which indirectly determine the adopted ventilatory pattern by the subject (frequency and depth) depending on the level of ventilation.

Once f_R is determined, the neural control, derived from the respiratory centers, establishes the muscular activity integrating the total ventilatory demand and modulating it, in turn, by an auto-rhythmic and square signal. Such a signal determines each breathing cycle (T_TOT) by using a relation 1.5:4 to define T_I and T_E, as follows:

In the case of assisted mechanical ventilation, internal neural activity is decreased. Depending on the type of ventilatory assistance, the respiratory period would be determined by the ventilator, the subject or the interaction between both.

The cardiovascular subsystem allows simulating the heart nature pulse and blood flow through the pulmonary and systemic circulations. Unlike RS1, this subsystem includes several processes like atria-ventricular mechanics, circulatory hemodynamics, SA node, change of total peripheral resistance and baroreflex (see Figure S2). Through these mechanisms, the system calculates the arterial blood pressure, ABP, heart period, HP, cardiac output, CO, and blood flow to lung for gas exchange depending on inputs from the autonomic control system, the respiratory system, and the sleep control system (Cheng et al., 2010).

Additionally, considering Equations (7, 8) of RS1, metabolic dynamics of CO₂ and O₂ were incorporated in RS2 to simulate exercise stimuli.

A third model called RS3 has been proposed in this study to get a completed and detailed model with a more appropriate dynamic response to exercise stimuli. This model is based on RS2 to take advantages of its associated subsystems, but mainly two key features from RS1 are replaced: (a) the estimation of ventilatory demand as a function of P_bCO2, P_aCO2, P_aO2, and MRV and (b) the adjustment of breathing pattern by using the optimization criteria set by Otis et al. (1950).

Schematic diagram of RS3 is shown in Figure 2. The main changes, highlighted in orange color, are associated with the replacement of “ventilatory drive” and “respiratory rhythm” blocks of the original model RS2 by a clock pulse generator, similar to RS1 (see Figures S1, S2). The latter allows determining the onset and the end of each breath as follows: (a) the “mean value detector” block determines average values of P_bCO2, P_aCO2, and P_aO2; (b) the “ventilation controller” block uses previous mean values and MRV to calculate V˙A by using Equation (3) and, finally, (c) the “frequency optimizer” block computes f_R using the Otis' Equation (Equation 9).

Moreover, three significant changes were added to the gas exchange subsystem, which is also highlighted in orange color in Figure 2, according to:

During inspiration,

During expiration,

where V_LCO2 and V_LO2 denote the volume storage of CO₂ and O₂ in the lungs. Likewise, Ṗ_ACO2 and Ṗ_AO2 represent the alveolar partial pressures of CO₂ and O₂; C_aCO2, C_aO2, C_vCO2, and C_vO2 symbolize the arterial and venous concentrations in body tissues of these gases, respectively; C_BCO2, C_BO2 denote the brain concentrations and Q_T and Q_B the blood rate in body and brain tissues.

where V_BCO2 and V_BO2 denoted volume storage of CO₂ and O₂ in the brain.

These changes provide to RS3 a more detailed gas exchange plant, and although they do not improve its transient response, they provide additional information that can be useful in future studies connected to, for example, analysis of brain-tissues relationship implicated in the pulmonary gas exchange.

Furthermore, to guarantee that ventilation generated by the mechanical plant matched ventilatory demand, the neural signal N(t) in RS3 was adjusted through the following expression:

where,

RC denotes the muscle constant time (0.060s) and VC the vital capacity (5 L).

Like in RS1 and RS2, metabolism dynamic of CO₂ and O₂ was incorporated in the gas exchange plant with the aim of simulating exercise stimuli, see Equations (7, 8).

Responses and features of the three models were evaluated considering different levels of exercise. In this stimulus, the consumption of O₂ and the production of CO₂ rise significantly increasing the ventilated values of CO₂ and O₂. For this reason, a step input of V˙CO2 and V˙O2 from rest (0.20 and 0.25 L/min), to moderate exercise (1.20 and 1.50 L/min), under conditions of normoxia, was considered to analyze their transient responses. Then, similarly, 11 equidistant step inputs among such intervals were taken into account to evaluate their stationary responses. These values were selected considering those published in Mañanas et al. (2002) and Guyton (2015) for moderated exercise and experimental data obtained in the CPET test (see Figure 1). Additionally, a sensitivity analysis was carried out with RS3 to assess the individual roll of the neurogenic and neuro-humoral mechanisms implemented to simulate exercise.

To assess the prediction capability of each model, experimental values of V˙CO2 and V˙O2 were used to simulate exercise stimulus. Then, the output cardiorespiratory variables predicted by the models were analyzed concerning the ones obtained experimentally. Comparison of each model response regarding experimental data was evaluated quantitatively by the prediction error (PE) calculated from the following variables:

PE was calculated by measuring percent differences between simulated, SIM, and experimental, EXP, variable as follows:

where k denotes the number of samples. Overall prediction error was obtained averaging the PE for all variables.

Regarding HR, due to RS1 does not provide direct information about it, this variable was obtained indirectly from the cardiac output (Q), considering a constant stroke volume (SV = 70 mL) and using the expression HR = Q/SV (Fincham and Tehrani, 1983; Batzel et al., 2007). Although this is a simple approximation, it allows comparing the cardiac response of RS1 with the other models.

Non-parametric tests, Friedman and Wilcoxon-Mann-Whitney (WMW), were used to identify statistical differences between prediction capability with a significance level of ρ = 0.05. The former was used in order to find differences between the model errors, and the latter to identify the model with the best fitting to experimental data. Each simulation was run once due to models are determinists (i.e., their responses do not change if the initial conditions and stimulus step size remain unchanged).

The models RS1, RS2, and RS3 can be interactively tested through a Matlab app, which is available at https://bioart.upc.edu/en/virtual-laboratories/modules upon query. Additionally, median values and interquartile distances of experimental data analyzed in this study as well as the stimulus levels used during simulation models are also provided.

Figure 5 shows the final values of analyzed variables in function of different levels of exercise. Eleven step inputs from 0.2 to 1.2 L/min for V˙CO2, with P_IO2 fixed to its sea level value (150 mmHg), were applied to simulate the models and analyze their steady-state responses under moderated exercise.

Excepting P_aCO2, variables of all models increased followed similar trends: T_I and P_aO2 decreased while f_R, V_T, and V˙Eincreased with higher stimuli. Static performance in RS1 and RS3 for T_I, f_R, V_T, and V˙E was quite similar since both models include the same neural ventilatory control to calculate V˙E, optimize f_R and, therefore, set T_I and V_T.

The slight differences found among the model responses were mainly due to: (a) model basal values no related to exercise, (b) differences among the gas exchange plant, and (c) the control systems implemented in each model. The steady-state values of P_aCO2 and P_aO2 further evidence these differences. RS1 and RS3 were able properly to regulate P_aCO2 and P_aO2, while RS2 converged to values very different from those expected during moderate exercise (P_aCO2 ≈ 40 mmHg and P_aO2 ≈ 104 mmHg) (Guyton, 2015). For RS1 and RS3, the proper regulation of such gases could be because V˙E was adjusted to match the ventilatory demand generated by exercise.

Regarding HR, in RS2 and RS3, this variable raised slightly and linearly with the stimulus level, reaching approximately an increase of 20% at the highest stimulus. On the contrary, in the RS1 model, HR increased considerably (up to 130%), moving away from the expected range for this type of stimulus.

Figure 6 shows the results obtained in the sensitivity analysis carried out in RS3. V˙E, P_aCO2 and P_aO2 are shown in two different conditions: first, with all mechanisms working and second, without the neurogenic mechanism MRV (exercise neural control component). Results are compared with those obtained with RS2.

Selective elimination of neurogenic mechanism MRV leads to a slower increase of V˙E during exercise with a consequent and important drop of P_aO2 at the onset of exercise (up to 50 mmHg) and an exponential increase in P_aCO2 (up to 12 mmHg from baseline). During the stimulus, the final ventilation value of RS3 differs only by 5 L/min from the value reached when all the mechanisms are considered, case in which P_aCO2 and P_aO2 vary slightly from their values at rest. This is because the absence of a neurogenic mechanism is compensated by central and peripheral chemoreceptors, which are greatly stimulated by the large decompensation of arterial pressures, especially P_aCO2. On the other hand, the action of the humoral mechanism (central chemoreceptors response) of RS2 generates a ventilation inferior to that obtained by RS3, a fact that contributed to a greater increase of P_aCO2 for the same stimulus level.

Figure 7 shows the contribution from each controller component of RS3 regarding its basal value. In this case, it can be seen that, during exercise, the presence of the exercise neural control component (neurogenic mechanism) generates a change in ventilation of up to 16 L/min (67%) while central and peripheral chemoreceptors contribute only with 8 L/min (33%) in the total ventilation. For the latter, this contribution is mainly generated by the changes given in P_aCO2 due to the relationship between V˙A and P_aCO2 (i.e., changes in V˙A by P_aO2 only are significant when this is lower than ≈ 60 mmHg, see Equation 4).

Given the nature of the stimuli used to simulate the models from experimental data (a progressive increase, contrary to step function), only responses in “transient regime” were analyzed. Figure 8 shows simulation results in median and interquartile distance. In this case, both experimental data and model responses of V˙E, f_R, and V_T increased with the increment of exercise stimulus whereas T_I decreased slightly. Regarding P_aCO2 y P_aO2, RS1 and RS3 showed an appropriate regulation. For RS2, these latter variables took values a little apart from experimental data at the highest stimulus levels. Respect to HR, RS2 and RS3 showed a good fitting whereas RS1 presented values far away of experimental data.

Figure 9 shows the median values and interquartile distance of prediction error calculated from Equation (24). Friedman test showed statistically significant differences in all variables except f_R and V_T, where the ρ-values were 0.149 and 0.122 respectively. In general, RS3 was the model with the best adjustment. It reached the lowest overall PE (13.76 vs. 20.92% for RS1 and 16.57% for RS2). Wilcoxon signed-rank test showed statistical differences between RS3 and the other models (ρ < 0.01) whereas RS1 and RS2 did not have any statistically significant difference (ρ = 0.432). Regarding each variable, RS3 presented the best fitting for V˙E, T_I, P_aO2, P_aCO2, and HR with ρ < 0.01 for the first three variables. RS1 and RS2 showed a better adjustment for f_R and V_T, respectively, but without any significant statistical difference respect to RS3 (see Figure 9).

On the other hand, HR in RS1, which was calculated from Q˙T and a stroke volume of 70 mL, presented a very high prediction error. In this case, the absence of a more complex cardiovascular subsystem impeded a suitable simulation of HR response to exercise.

In RS3, the abrupt increase of V˙E at the onset of exercise and, therefore, the initial changes in P_aCO2 and P_aO2, was mainly determined by the “MRV” component incorporated in its controller (see Equation 4). This component played an important role in the determination of anticipatory ventilatory response and allowed simulating, with a straightforward approach, the neurogenic mechanism of ventilation to exercise, i.e., the so-called central command (direct activation of the respiratory control centers by the locomotor stimulus), chemoreceptors in large vessels and mechanoreceptors located in exercising muscles (Dempsey and Smith, 2014), see Figures 6, 7. Although, an only “exercise” component in the respiratory controller could not be enough to simulate each underlie mechanism to this type of control, similar results to those obtained in RS3, especially with regard to the behavior of V˙E, P_aCO2, and P_aO2, have been reported in several sources (Krogh and Lindhard, 1913; Pearce and Milhorn, 1977; Whipp et al., 1982; Mateika and Duffin, 1995; Turner et al., 1997; Whipp and Ward, 1998; Magosso and Ursino, 2005; Wasserman et al., 2011; Parkes, 2013; Guyton, 2015).

Particularly, the changes in ventilation, blood flow, pulse rate, respiratory exchange and alveolar CO₂ tension, which take place in men during the first minute of light (moderate) or heavy work, presented in Krogh and Lindhard (1913) agree with the transient responses obtained in RS3, see Figure 4. They found in all the cases examined (six subjects, three of them subject trained to sudden exertions) a sudden rise in ventilation when an exercise step was used to stimulate the subjects. This abrupt increase was greater with heavier work. Although differences were found between trained and no-trained subjects, the latter presented a higher increase in V˙E that, although was not as pronounced as in trained subjects, was always present even during moderate exercise. Regarding gases regulation, as in RS3 a considerable fall in alveolar CO₂ tension was detected, which was, likewise, followed by a CO₂ increase. A similar behavior was also reported in Parkes (2013). They found a nearly exponential increment of ventilation during exercise, while P_aO2 decreased regarding its rest value at the beginning of stimulus for, then, growing up exponentially to its steady state. Moreover, the simulated data reported in Magosso and Ursino (2005), which were found according to experimental dynamic responses of human subjects to bicycle exercise published in Pearce and Milhorn (1977), were also consistent with RS3.

It is important to note that respiratory drive and therefore, ventilatory response (V˙A) in RS3 (and RS1) is related to P_bCO2, P_aCO2, P_aO2 and MRV through Equation (3). Thus, there is a linear relationship with P_bCO2, P_aCO2, and MRV and an exponential dependence with P_aO2 (see Equation 4). So, unitary changes in P_bCO2, P_aCO2, and MRV would modify V˙A by a quantity of 0.2332, 0.2025, and 1 of its basal value (V˙Abasal). By contrast, decreases of P_aO2 would exponentially increase V˙A only if P_aO2 < 104.9 (see Figure 7).

On the other hand, in this study, exercise stimulus was used to generate different levels of ventilation and evaluate how the models adjust ventilation and breathing pattern to accomplish metabolic demand and regulate arterial blood gases. RS1 and RS3 used an optimization principle that allows adjusting f_R by minimizing work of breathing. In this sense, the minimization of WOB has been extensively considered as a control criterion to adjust the breathing pattern (Yamashiro and Grodins, 1971; Poon et al., 1992; Serna Higuita et al., 2014; Serna et al., 2016). Moreover, recent formulations seek to describe how sensory information influences the dynamics of respiratory rhythm under the hypothesis that “respiratory rhythms arise from the interplay of central rhythm generation circuits, biomechanics and feedback from peripheral signaling pathways” (Butera et al., 1999a,b; Diekman et al., 2017). Particularly, rhythmogenesis is investigated in Diekman et al. (2017) in a simple model of close-loop control, incorporating biomechanics, oxygen handling, metabolism and chemo-sensation. In such study, the Butera-Rinzel-Smith model (Butera et al., 1999a) of bursting pacemaker neurons in the preBötzinger complex is adopted as their central pattern generator. Although, peripheral processes modeled in Diekman et al. (2017) are not as comprehensive as RS1, RS2, and RS3, this approach to simulate the dynamics of rhythm respiratory could be interesting to provide RS3 a more realistic pattern generator.

Finally, highly regulated neural inputs are critical to maintaining normal cardiovascular function. Although the cardiovascular central command during exercise is typically associated with a perception of effort, there is not a clear understanding of the role of central command in the integration of sensory information that can define more completely the relevance of central command for the neural control of exercise (Williamson, 2010).

Regarding the model goodness of fit, RS3 presented the overall best adjustment to experimental data with the lowest prediction error and an improvement of 17% respect to RS2 (overall PE = 13.51%, see Figure 9). RS3 also presented the lower prediction errors for T_I, V˙E, P_aCO2, P_aCO2, and HR with statistically significant differences for the first three.

While the reduction of the prediction error in RS3 (13.76%) is lower concerning RS2 (16.57%) than to RS1 (20.92%) when using experimental data, in both cases, these reductions were statistically significant (ρ < 0.01). This improvement is more evident if only the respiratory variables T_I, f_R, V_T, and V˙E are considered. In this case, the prediction errors are 23.0, 26.4, and 16.3%, for RS1, RS2, and RS3 respectively, and RS3 presents an improvement of 38% regarding RS2. Such difference is due to prediction errors of P_aCO2 and P_aO2in all models were relatively small by the magnitude of these variables.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.