An Improved Dynamic Model for the Respiratory Response to Exercise.

Respiratory system modeling has been extensively studied in steady-state conditions to simulate sleep disorders, to predict its behavior under ventilatory diseases or stimuli and to simulate its interaction with mechanical ventilation. Nevertheless, the studies focused on the instantaneous response are limited, which restricts its application in clinical practice. The aim of this study is double: firstly, to analyze both dynamic and static responses of two known respiratory models under exercise stimuli by using an incremental exercise stimulus sequence (to analyze the model responses when step inputs are applied) and experimental data (to assess prediction capability of each model). Secondly, to propose changes in the models' structures to improve their transient and stationary responses. The versatility of the resulting model vs. the other two is shown according to the ability to simulate ventilatory stimuli, like exercise, with a proper regulation of the arterial blood gases, suitable constant times and a better adjustment to experimental data. The proposed model adjusts the breathing pattern every respiratory cycle using an optimization criterion based on minimization of work of breathing through regulation of respiratory frequency.

Introduction

The respiratory system is a complex and feedback system which is responsible for supplying sufficient oxygen (O2) for metabolism and eliminating carbon dioxide (CO2) produced by metabolic reactions in order to keep the homeostasis of arterial blood gases and pH in any situation and, particularly, during exercise (Duffin, 2013). For achieving this goal, the respiratory control system regulates pulmonary ventilation so that at equilibrium gas exchange in the lungs matches metabolism: O2 provision and CO2 elimination in the lungs equal O2 consumption and CO2 production in the tissues (Bell, 2006; Duffin, 2013; Guyton, 2015).

The regulation of O2 is achieved keeping a partial pressure of O2 (P), in the arterial blood, that saturates arterial hemoglobin and provides a sufficient gradient to supply the tissue metabolic. Because of the hemoglobin can be saturated within a wide range of P, oxygen is mainly regulated when its saturation falls to hypoxia limits (P < 70 mmHg) (Roussos and Koutsoukou, 2003; Duffin, 2013). On the other hand, due to CO2 diffuses more easily than O2 and quickly reacts with H2O, generating hydrogen ion concentrations ([H+]), the regulation of CO2 is more difficult to achieve although it can carry out by controlling the partial pressure of CO2 (P) and, therefore, [H+] (Duffin, 2013).

Pulmonary ventilation () is normally generated through controlled contraction and relaxation of respiratory muscles during inspiration and expiration, respectively. The pressures generated by them combined with the airway flow resistance, and the lung elastance determine the airflow and lung volume, correspondingly (Bianchi et al., 1995). For this, the respiratory control system, which is responsible for automatic control of breathing, (a) receives and integrates several afferent inputs from both central and peripheral chemoreceptors and pulmonary sensors, to determine ventilatory demand and efficiently adjust tidal volume (V) and respiratory frequency (f); and (b) provides signals to the phrenic and intercostal moto-neurons, which drive the diaphragm and intercostal muscles (Duffin, 1994; Bianchi et al., 1995). Due to this, the respiratory control system is frequently seen as a pattern central generator that seeks to maintain healthy levels of O2, CO2, and pH in body and brain tissues (Feldman et al., 2013; Richter and Smith, 2014). This control can also be affected by behavioral inputs such as speech, voluntary control of respiratory muscles and wakefulness state (Duffin, 2013).

The central chemoreceptors, which are located in the medulla, sense increases of [H+] in their local environment and, therefore, increases in the brain partial pressure of CO2 (P) (Nattie and Li, 2009). On the other hand, peripheral chemoreceptors, which are located in the carotid bodies, sense changes in the arterial partial pressures of CO2 and O2 (P and P, respectively) (Marshall, 1994). In this way, both the central and the peripheral chemoreceptors allow matching pulmonary ventilation to tissue metabolism via a chemical arc that includes the respiratory control system, the chemoreceptors and respiratory muscles (Duffin, 2013). Recently, changes in P have been related to increases of chemoreceptor sensitivity to [H+] at extreme hypoxia (Blain et al., 2010; Duffin, 2013).

During mild and moderate exercise, metabolic rate and pulmonary exchange increase as a result of muscular activity, such increases quickly produce a higher ventilation which seeks to prevent hypercapnia (due to increase in CO2 production) and hypoxia (due to increase in O2 consumption) (Turner et al., 1997; Whipp and Ward, 1998; Haouzi, 2006). Because of increases nearly immediately at the onset of exercise, levels of CO2 and O2 remain practically unchanged from their rest values so, for this reason, a respiratory drive is not compatible with sensing an error signal transported in the blood (Duffin, 1994).

In recent decades, several clinical and experimental investigations have been carried out to determine the control mechanisms responsible for adjusting ventilation during exercise. A challenging aspect of this ventilatory stimulus is that increases while the brain and arterial partial pressures of CO2 and O2 remain almost unchanged. A generally accepted theory is the so-called neurohumoral theory (Turner, 1991; Mateika and Duffin, 1995; Turner et al., 1997; Whipp and Ward, 1998), which is mainly based on the respiratory system response to a step workload. In this theory, respiratory dynamic response from rest or light to moderate exercise is characterized by three phases: the first one, determined by a sudden increase of , the second one, by a gradual and exponential increase of and the third one, by its value in the steady state. The “abrupt” increase of is usually attributed to neurogenic mechanisms, since this increase is considered too fast to be explained by humoral agents, such as central and peripheral chemoreceptors, due to delayed transport. Until now, such mechanisms are not yet well-understood, because they seem to involve feed-forward control systems or learned processes that have not been clearly figured out (Bell, 2006; Williamson, 2010).

Some studies establish that the behavior of ventilatory response during moderate exercise is related to the frequency of limb movement, and the force carried out by exercising muscles (Duffin, 1994). On the other hand, other studies establish that such response is mainly based on factors related to gas exchange more than factors related to the motor activity (Haouzi, 2006). A broad review on the mechanisms currently implicated in the control of breathing at the onset of exercise from a perspective of an integrated system can be found in Bell (2006) and Duffin (2014).

Many empirical and functional models have been proposed in the literature to describe numerous aspects of the respiratory system (Fincham and Tehrani, 1983; Butera et al., 1999a,b; Cheng et al., 2010; Williamson, 2010; Tsai and Lee, 2011; Cheng and Khoo, 2012; Serna Higuita et al., 2014; Serna et al., 2016; Diekman et al., 2017). Due to the primary goal of this system is to regulate the CO2 and O2 in the brain and body tissues, ventilatory stimuli like exercise, hypoxia and hypercapnia are frequently used to evaluate the performance of such models. Exercise has been one of the most used ventilatory stimuli for validating this kind of models and their control mechanisms (Magosso and Ursino, 2005; Hermand et al., 2016).

Particularly, Hermand et al. (2016) have presented a mathematical model that allows analyzing the mechanisms responsible for the instability of the respiratory control system under simultaneous metabolic (exercise), and environmental (hypoxia) stresses. In this case, the model analysis is mainly focused on variations of and f taking into account several settings to simulate the central and peripheral chemoreceptor responses. On the other hand, Maggoso and Ursino have presented a respiratory model that allows obtaining the transient and steady-state cardiorespiratory response to exercise with a good performance, but it does not include a comprehensive cardiovascular model like that published in Cheng et al. (2010) and the variables related to breathing pattern (i.e., inspiratory time, respiratory frequency, and tidal volume) are not evaluated. Our group has previously analyzed and developed several models and tools in the framework of this research (Mañanas et al., 2003; Hernandez et al., 2008; Serna et al., 2010). However, although the steady-state response of such models has been thoroughly evaluated (Mañanas et al., 2003, 2004), the transient behavior had not been studied enough. For these reasons, there is a need to provide computational models that allow, with a physiological meaning, simulating a comprehensive dynamic response of the respiratory system under ventilatory stimuli like exercise.

The aim of this study is two-fold. Firstly, to propose an improved model of the respiratory system that allows simulating its dynamic response to ventilatory stimuli like exercise. Secondly, to compare its transient and static responses with those obtained from two known respiratory models (Fincham and Tehrani, 1983; Cheng et al., 2010; Cheng and Khoo, 2012) and from which the proposed model is based on. The comparison was performed by simulation and by using experimental data from healthy subjects who carried out the cardiopulmonary exercise testing (CPET). The former allowed to analyze differences among model responses and, the latter, to assess the prediction capability of each one. Due to complex structures and mechanism that may be involved during exercise, this study was focused on the dynamic response analysis of respiratory system under moderate exercise (below the lactate threshold). All models were implemented in SIMULINK/MATLAB®.

Materials and methods

Experimental data

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 (), tidal volume (V), inspiration time (T), expiration time (T), O2 consumptions (), CO2 productions , expired fraction of O2 (P) and CO2 (P), and heart rate (HR). Taking into account that no subject had gas exchange impartments, arterial blood pressures of CO2 and O2 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, P and P whose values remained almost constant, all variables increased with the exercise workload. This latter was related to increments of both and that, in turn, had a linear relationship with .

Figure 1

Experimental data distribution in function of time. (Left) Experimental exercise stimuli used to simulate the models under study: ventilated values of CO2 () and O2 (). Rest and moderate exercise levels are also marked in both variables. (Right) Experimental data of the cardiorespiratory variables analyzed: inspiratory time (T), respiratory frequency (f), tidal volume (V), minute ventilation (), arterial partial pressure of CO2 (P) and O2 (P), and heart rate (HR). Gray dots show original experimental data, bar plots and red lines display the median values obtained in 50-s consecutive time intervals and error bars show the interquartile ranges calculated in each interval.

Respiratory system models

RS1 model

RS1 is the respiratory model described in Fincham and Tehrani (1983) which has been extensively studied by many researchers for teaching and research purposes (Mañanas et al., 2003; Tehrani et al., 2004; Batzel et al., 2007; Hernandez et al., 2008). Additionally, it has been used to predict the effects of ventilatory settings on blood gases of mechanically ventilated patients and to adjust the ventilation parameters to optimize such a treatment both adult and infant (Tehrani and Abbasi, 2012). A schematic diagram of this model is provided in the Supplementary Material Section (see Figure S1).

RS1 includes a self-adaptive and discrete controller (Priban and Fincham, 1965), which self-adjusts the ventilation and the breathing pattern, every respiratory cycle, from signals captured by different afferent pathways, as it happens physiologically. It incorporates complex peripheral processes like the gas exchange between lungs and body tissues and brain, transport delays due to blood circulation and blood gas dissociation.

To adjust , RS1 controller continuously receives information sent by the central (P) and peripheral chemoreceptors (delayed by transport through the circulatory system, P′ and P′), computes the mean value of such pressures in each respiratory cycle and sends these signals to the ventilatory controller which calculates the alveolar ventilation, cycle to cycle, through the following expression:

where,

In this case, is defined by the sum of and dead space ventilation (), where V is, in turn, a function of :

The first term of Equation (3) represents the response generated by the central chemoreceptors, where P is defined as a function depending on the brain venous and cerebrospinal fluid partial pressures of CO2 (PvbCO and PCSFCO, respectively) and, therefore, on P (see Figure S1). The second and third terms define the response of peripheral chemoreceptors to the arterial partial pressure of CO2 and O2, respectively; the fourth term, the neural control related to metabolism or exercise and, the last one β, is a controller constant equals 17.4. is the basal alveolar ventilation for a healthy adult, and it is ~0.0673 L/s (Fincham and Tehrani, 1983).

MRV, the neural impulse derived from metabolism, represents in RS1 the neurogenic mechanism of the respiratory control system during exercise, and it is determined by the following expression:

where MRR is the metabolic ratio defined by the current and basal of metabolic rates of the brain (MR) and tissues (MR), as follows:

and an exercise level increment can be simulated by a change in the input RTT (the magnitude of the step defining the final value of MR) which is given by:

In this way, Equations (6–8) allow determining the metabolic dynamics of body and brain tissues as well as the subject's neural response during exercise. These equations represent first order dynamic systems whose response rates are determined by the time constants τ4 = 50 and τ3 = 30, respectively (Fincham and Tehrani, 1983).

There are essential conditions on which Equation (3) is invalid. Under conditions of acute hypoxia associated with low levels of the arterial partial pressure of CO2, this equation can produce a negative value of , clearly inadmissible. In this case, apnea occurs and immediately equals zero.

Once is determined, which in turn is equal to the product of V and f, breathing pattern is adjusted by the regulation of f through an optimization criterion based on minimizing the work of breathing (Otis et al., 1950), as is shown in the following expression:

where, R and E are the resistance and the elastance of the respiratory system, respectively. For that, the RS1 controller generates a neural signal at the beginning of each cycle such that:

In RS1, the blood flows control, this is, the “brain blood flow controller” and “the cardiac output controller” were modeled through algebraic relationships that allow calculating brain blood flow () and total cardiac output () depending on partial pressures of CO2 and O2 and metabolic rate ratio (MRR).

RS2 model

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 (Pand P), due to different sectors of anatomic dead space, using a first-order dynamic approach to calculate the arterial partial pressures of CO2 and O2 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) and peripheral (D) chemoreceptors. In this model, the central chemoreceptors only respond to variations in the brain partial pressure of CO2 (P) while the peripheral chemoreceptors are influenced by the arterial partial pressures of CO2 (PaCO) and the oxygen saturation in arterial blood (SAO2) and their multiplicative interaction (Khoo, 1990). During wakefulness, the total ventilatory demand (D) is defined by the sum of the central and peripheral chemoreceptors responses as follows:

where, I, I, and I represent the central and peripheral chemoreceptors activation threshold and, they are equal to 45, 38, and 102.4, respectively. In this case, P is controlled by the metabolic rate (MR) and the brain blood flow (Q) and it is defined as a function of P (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 is the basal frequency and T and T are thresholds of D that determine the behavior of respiratory rate. For the last two options of Equation (9) the respiratory rate varies linearly respect to D with a slope established by the scaling factors S1 and S2, which indirectly determine the adopted ventilatory pattern by the subject (frequency and depth) depending on the level of ventilation.

Once f 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) by using a relation 1.5:4 to define T and T, 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 CO2 and O2 were incorporated in RS2 to simulate exercise stimuli.

RS3 model

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, P, P, 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, P, and P; (b) the “ventilation controller” block uses previous mean values and MRV to calculate by using Equation (3) and, finally, (c) the “frequency optimizer” block computes f using the Otis' Equation (Equation 9).

Figure 2

Schematic diagram of the proposed respiratory system model RS3. A model based on RS2 (Cheng et al., 2010; Cheng and Khoo, 2012) that integrates key features of RS1 (Fincham and Tehrani, 1983). The changes introduced in RS2 have been highlighted in orange. In the model, mean values of P′, P′ and P (P) and a metabolically related neural drive component to the ventilation (MRV) are used by the central neural control system to compute the alveolar ventilation () and the total ventilatory drive (). Then, is used to adjust the respiratory frequency (f) each breathing cycle following the optimization principle set in Otis et al. (1950). A detailed description of model variables is provided in the Supplementary Material Section (Table S1).

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

Production of CO2 and consumption of O2 into brain tissues were included into the mass balance equations for the rate of change of the lung CO2 and O2 volumes, because they constitute ~20% of the metabolic rate at the basal level.

During inspiration,

During expiration,

where V and V denote the volume storage of CO2 and O2 in the lungs. Likewise, Ṗ and Ṗ represent the alveolar partial pressures of CO2 and O2; C, C, C, and C symbolize the arterial and venous concentrations in body tissues of these gases, respectively; C, C denote the brain concentrations and Q and Q the blood rate in body and brain tissues.

The exchange in the brain considered as follows:

where V and V denoted volume storage of CO2 and O2 in the brain.

Calculation of V as a function of (see Equation 5).

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 CO2 and O2 was incorporated in the gas exchange plant with the aim of simulating exercise stimuli, see Equations (7, 8).

Simulation

Responses and features of the three models were evaluated considering different levels of exercise. In this stimulus, the consumption of O2 and the production of CO2 rise significantly increasing the ventilated values of CO2 and O2. For this reason, a step input of and 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.

Prediction error

To assess the prediction capability of each model, experimental values of and 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:

, T, f, and V, which provide information about ventilatory strategy or breathing pattern adopted for each model (controller) to adjust ventilation,

HR, which provides information about cardiac activity, and

P and P, which allow assessing the regulation of CO2 and O2 respectively.

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.

Statistical analysis

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).

Availability

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.

Results

Simulations

Transient response

Figure 3 shows a time series breath-to-breath of airflow signal for each model at the onset of exercise when a step input of 1.25 L/min and 1.50 L/min was applied to simulate each model (0.2 and 0.25 L/min values were considered as rest levels, respectively). In all models, it can be seen an increase in ventilation, which is higher and faster in RS3.

Figure 3

Time series breath-to-breath of airflow signal at the onset of exercise when a step input was used to simulate the respiratory system models RS1, RS2, and RS3 (see text).

Figure 4 shows the transient response obtained by the models. An additional interval, after of the stimulus, was also considered to simulate the recovery phase. Values of each variable are shown regarding their basal or rest values. Table 1 shows the constant times obtained in each of them.

Figure 4

Transient responses of analyzed variables—ventilation (), inspiratory time (T), respiratory frequency (f), tidal volume (V), arterial CO2 pressure (P), arterial O2 pressure (P), and heart rate (HR)—when a step exercise input ( 1.20 L/min and 1.50 L/min) was used to simulate RS1, RS2, and RS3. Variations of each variable are showing regarding their respective values in rest.

Table 1

Rise time (τ) of the dynamic response of RS1, RS2, and RS3 obtained when a step exercise input of 1.20 L/min and 1.50 L/min was used to simulate them.

Variable	τ (mm:ss)
	RS1	RS2	RS3
TI	00:38	03:00	00:45
fR	00:52	04:15	01:15
VT	01:24	01:55	02:45
V˙_E	01:22	03:15	02:40
PaCO2	00:57	04:20	03:30
PaO2	00:44	00:55	01:20
HR	00:21	04:55	03:25

Regarding variables related to breathing pattern, T, f, and V, for RS1 and RS2, an exponential behavior was observed both exercise and recovery, being RS1 faster than RS2. For RS3, such responses were characterized by two main temporal phases during exercise: an initial phase determined by an instantaneous increase of almost 100% of and a second phase defined by a gradual increase of to its steady state-value. During recovery, and V presented a similar behavior and opposite to that obtained during exercise while f and, therefore, T were given by the dynamics defined in Equation (9) (Otis et al., 1950).

Regarding the variables related to gas exchange and for RS1, P presented initially a slight overshoot, which was followed by a decrease and subsequent exponential evolution toward a value close to its basal level. By contrast, P exhibited an initial slight drop, which was followed by a positive overshoot of about 10% and an exponential evolution that, as P, converged to a value close to that obtained during rest. This behavior was also reported by the same authors in Fincham and Tehrani (1983). For RS2, P increased exponentially until overcomes its basal value by 40% (≈56 mmHg), while P rapidly declined to 50% (≈50 mmHg), reaching values close to hypoxia (P < 70 mmHg) (Roussos and Koutsoukou, 2003), for then evolves exponentially at a steady-state value of 80 mmHg.

For RS3, P was characterized by an initial decrease, of about 15% from its basal value, followed by an exponential growth toward a steady-state value 10% higher than its basal value. P was defined by, first, an overshoot that, after ~1 min, evolved almost exponentially toward its stationary value with a negative overshoot that was not as critical as that found in RS2. During recovery, both P and P presented similar and opposite behaviors in all models. In addition, evolution of P was slower than P, especially for RS2 and RS3 (see Table 1), possibly due to the larger storage capacity available for CO2 (V = 15 L) compared to O2(V = 6 L). A similar behavior has been reported in Mateika and Duffin (1995).

Finally, for all models, HR exponentially evolved with different velocities toward its final value. Particularly for RS1, the response time of HR was lower than that obtained by the other models, and its increase from baseline was quite higher (≈170 bpm) than that expected for moderate exercise (85–110 bpm), see Figure 1.

Steady state response

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 , with P fixed to its sea level value (150 mmHg), were applied to simulate the models and analyze their steady-state responses under moderated exercise.

Figure 5

Steady state values obtained when step inputs of exercise from 0.25 to 2.50 L/min of were used to simulate the respiratory system models under study.

Excepting P, variables of all models increased followed similar trends: T and P decreased while f, V, and increased with higher stimuli. Static performance in RS1 and RS3 for T, f, V, and was quite similar since both models include the same neural ventilatory control to calculate , optimize f and, therefore, set T and V.

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 and P further evidence these differences. RS1 and RS3 were able properly to regulate P and P, while RS2 converged to values very different from those expected during moderate exercise (P ≈ 40 mmHg and P ≈ 104 mmHg) (Guyton, 2015). For RS1 and RS3, the proper regulation of such gases could be because 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.

Analysis of RS3 control mechanisms

Figure 6 shows the results obtained in the sensitivity analysis carried out in RS3. , P and P 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.

Figure 6

Response of RS3 to a step exercise input under two conditions: firstly, with all mechanism of respiratory control operative and, secondly, without the neurogenic mechanism MRV (see Equation 3). Results are compared with those obtained in RS2 for the same stimulus. Variations of each variable are showing regarding their respective rest values.

Selective elimination of neurogenic mechanism MRV leads to a slower increase of during exercise with a consequent and important drop of P at the onset of exercise (up to 50 mmHg) and an exponential increase in P (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 and P 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. 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 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 due to the relationship between and P (i.e., changes in by P only are significant when this is lower than ≈ 60 mmHg, see Equation 4).

Figure 7

Contribution of each controller component of RS3 on ventilatory response (see Equation 3). Results are shown regarding their respective rest values.

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 , f, and V increased with the increment of exercise stimulus whereas T decreased slightly. Regarding P y P, 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 8

Transient responses of models analyzed in this study. Median values and interquartile distance of experimental and simulated data are presented for each variable and model. Gray shadows show the IQR of experimental data and error bars the IQR of simulated data obtained by each model.

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 and V, 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 , T, P, P, and HR with ρ < 0.01 for the first three variables. RS1 and RS2 showed a better adjustment for f and V, respectively, but without any significant statistical difference respect to RS3 (see Figure 9).

Figure 9

Prediction errors (PE) obtained by each one of analyzed models. Bars graph shows median values and error bars indicate interquartile distance (range of values from 25 to 75% quartile). Wilcoxon signed-rank test was used to find statistically significant differences between the obtained PE. (*) for ρ < 0.05 and (**) for ρ < 0.01.

On the other hand, HR in RS1, which was calculated from 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.

Discussion

Dynamic responses of three respiratory system models under exercise stimuli have been analyzed in this study. The performances of these models were evaluated by using two settings: first, through simulation to assess the realism of the dynamics of their responses and second, using experimental data to estimate the prediction capability of each model.

Dynamic response

Transient response of each model was evaluated by using a moderate exercise step input. In this case, RS1 presented a relatively fast response at the onset of exercise (see Table 1), with an exponential response for and an appropriate regulation for P and P at the end of the stimulus. RS2 presented an exponential growth in that was not enough to regulate the arterial gases. Due to this, P decreased dramatically, during the transient phase reaching hypoxemic values (P ≈ 50 mmHg) and P increased by 40% its basal value reaching hypercapnic values (P ≈ 56 mmHg). This fact invalidates its simulation for this kind of stimulus. Finally, RS3 presented at the onset of exercise an sudden increase of which was followed by an intermediate phase characterized by an exponential growth, and a final phase defined by its static value. The “abrupt” increase of caused P to decrease slightly and that P obtained an overshoot at the beginning of the stimulus. Then, in the intermediate phase, P decreased while P raised as a result of the increase in consumption of O2 and production of CO2, respectively, and its dissociation with (Whipp and Ward, 1991) and, at the final phase, P and P took values near their basal values (see Figure 4).

Regarding variables involved in breathing pattern (T, f, and V), it was found that, except RS3 at the onset of the stimulus, all variables evolved similarly and exponentially toward their steady-state values. Differences found among them were mainly given by the ventilatory controller implemented in each model. On the other hand, behavior of HR in RS2 and RS3 was very similar, and its values were within the range defined for moderate exercise, contrary to RS1, where HR took values far away to those expected (see Figures 1, 4).

Static response

The steady-state model responses were evaluated considering a sequence of step inputs from rest to moderate exercise. In all models, linearly increased with the stimulus level (see Figure 5). Changes between consecutive stimulus were quite similar. Respect to f, this variable also linearly augmented in RS2 according to Equation (14), while in RS1 and RS3 it was determined by Equation (9), which guarantees the minimum respiratory work of breathing (Otis et al., 1950). These adjustments affected the behavior of T and V due to their relationship with and f ().

On the other hand, RS1 and RS3 presented an appropriate regulation for P y P (the final values were closer to normal ones at rest, P ≈ 40 mmHg and P ≈ 104 mmHg) (Guyton, 2015). Particularly in RS3, P and P shown slight variation regarding their basal values such has been reported from several studies (Krogh and Lindhard, 1913; Pearce and Milhorn, 1977; Magosso and Ursino, 2005). RS2 did not show a proper regulation of the arterial gases due to P increased toward hypercapnic values (P ≈ 54 mmHg) while P decreased toward hypoxemic values (P ≈ 80 mmHg), lower than expected ones (Guyton, 2015).

Control mechanisms

In RS3, the abrupt increase of at the onset of exercise and, therefore, the initial changes in P and P, 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 , P, and P, 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 CO2 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 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 CO2 tension was detected, which was, likewise, followed by a CO2 increase. A similar behavior was also reported in Parkes (2013). They found a nearly exponential increment of ventilation during exercise, while P 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 () in RS3 (and RS1) is related to P, P, P and MRV through Equation (3). Thus, there is a linear relationship with P, P, and MRV and an exponential dependence with P (see Equation 4). So, unitary changes in P, P, and MRV would modify by a quantity of 0.2332, 0.2025, and 1 of its basal value (). By contrast, decreases of PaO2 would exponentially increase only if P < 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 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).

Prediction capability

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, , P, P, 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, f, V, and 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 and Pin all models were relatively small by the magnitude of these variables.

Conclusion

Three respiratory system models have been analyzed in this paper. Two of them published Fincham and Tehrani (1983), Cheng et al. (2010), and Cheng and Khoo (2012), named RS1 and RS2 respectively, and the other one is a model proposed in this study, called RS3 and based on the integration of key features of the first two.

The first analyzed model, RS1 (Fincham and Tehrani, 1983), is a complex model that adjusts and the breathing pattern by minimizing the work of breathing through regulation of f (see Equation 9) (Otis et al., 1950). It integrates several peripheral processes and self-adjust the ventilation and breathing pattern at the end of each breath from signals captured by humoral and neurogenic afferent pathways. Simulation of this model, under exercise stimuli, showed a good adjustment of and a proper regulation of arterial gases (P and P). This was the key feature that motivated us to use the neural controller of RS1 in RS3.

The second model analyzed, RS2 (Cheng et al., 2010; Cheng and Khoo, 2012) is a more comprehensive model. It integrates the interaction between the respiratory and cardiovascular systems and allows simulating the dynamic of cardiorespiratory control during wakefulness and sleep. However, unlike RS1, respiratory control is carried out through a proportional controller that does not take into account the work of breathing done by the subject. It adjusts in function of the brain and arterial partial pressures of CO2 and O2 and regulates f and V through lines predefined in Duffin et al. (2000). Moreover, RS2 considers a multiplicative interaction among peripheral chemoreceptors more than a higher sensitive to hypoxia than hypercapnia (Blain et al., 2010; Cui et al., 2012; Kumar and Prabhakar, 2012), such as it is described in RS1 (see Equations 3, 4). It also does not take into account the control neural performed by the central controller to changes in the subject's metabolic rates during exercise (Williamson, 2010; Duffin, 2014; Guyton, 2015). Furthermore, it is not clear why in RS2 both PaO2 and SaO2 were indicated as independent variables since peripheral chemoreceptors are more sensitive to PaO2 rather than SaO2 (Kumar and Prabhakar, 2012). Simulations of this model, using exercise stimuli, did not show a proper regulation of blood gases for both transient and steady-state responses. This fact is understandable if it considers that RS2 is a model designed for sleep-related studies and not for another kind of stimuli such as exercise. One of the challenges of this study was to adapt RS2 to simulate exercise.

The third analyzed model, RS3, was proposed in this study to take advantages of completeness and versatility of RS2 and some properties of RS1. This model was proposed integrating key features of RS1 into RS2. In this sense, the proposed model, RS3, provided to RS2 the ability to adjust the ventilation in function of (a) the brain partial pressure of CO2 and the arterial partial pressure of CO2 and O2 and (b) the tissue metabolic demand. Unlike RS2, RS3 controller allows describing the highest sensitive of such chemoreceptors to hypoxia (see second and third terms of Equations 3, 4). Additionally, this model also supplied to RS2 the capacity to adjust breathing pattern considering an efficiency criterion based on minimization of work of breathing through regulation of f (see Equation 9). These features enable RS3 to accomplish appropriate transient and stationary responses during exercise. Likewise, the model improvement is not only related to prediction error. RS3 showed transient responses faster than RS2 with a better physiological meaning. This was especially important in the dynamics found for P and P (see Figure 4), which were consistent with results published in previous studies (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). On the other hand, although RS3 is more complex and complete than RS1, there is not a higher complexity in RS3 respect to RS2, but there is a substitution of some blocks and equations (from RS1) that allowed getting the improvements mentioned above.

Even though RS3 showed a good regulation of ventilation and blood gases partial pressures and, unlike RS1, it provides information related to cardiac activity such as heart rate, stroke volume, and cardiac output, we are still far to reproduce a real response to this type of stimulus. One of the found handicaps is related to many theories that have been developed so far to describe the underlying mechanism to the cardiorespiratory response during exercise. Particularly, RS3 (like RS1) has a component in its controller that allows it to simulate exercise stimulus from metabolic rate ratio (MRR) and, therefore, from and (see Equations 3, 6). This feature allowed RS3 to reproduce transient responses similar to those reported in the literature, especially at the onset of exercise, and achieve a better performance when experimental data were used. We are aware that this “exercise component” is a straightforward approach to simulate this type of control and additional efforts will be necessary to reproduce detailed changes in ventilation due to, for example, the central command, large vessels chemoreceptors and mechanoreceptors located in exercising muscles. Nonetheless, we consider results obtained here represent one-step beyond simulation of this kind of stimulus because they allowed contrasting the models, taking advantages of their key features, reproducing transient responses more realistic from a physiological point of view and getting a better prediction error.

Another handicap is related to published studies include many variables that influence the subject's response: mode exercise (walking, pedaling), posture, initial conditions, state of the subject (trained, inexperienced, uncomfortable, anxious, anticipating, distracted, tired, etc.), and metabolic rate (Bell, 2006; Fadel, 2013; Duffin, 2014). On the other hand, due to the significant individual variation in the ventilatory responses reported in the literature, a better prediction of real breathing patterns can be achieved only by fitting procedures (i.e., estimation of individual parameters).

Finally, the proposed model RS3 will need a validation in a sample of healthy controlled subjects in different exercise conditions to go in deep in the controller mechanism of the respiratory system during exercise.

Author contributions

All authors made substantial contributions to the conception and design of the paper. LS: Carried out the model simulations, processed the obtained data, and drafted the review; MM, AH, and RR: Revised it critically for content; All authors approved the final submission of the document and agreed to be accountable for all aspects of the work.

Conflict of interest statement

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.