A novel strategy to fit and validate physiological models: a case study of acardiorespiratory model for simulation of incremental aerobic exercise

Citation: Sarmiento, C.A.; Serna, L.Y.; Hernández, A.M.; Mañanas, M.Á. A Novel Strategy to Fit and Validate Physiological Models: A Case Study of a Cardiorespiratory Model for Simulation of Incremental Aerobic Exercise. Diagnostics 2023,13, 908. https://doi.org/10.3390/ diagnostics13050908 Academic Editor: Fleur T. Tehrani Received: 24 January 2023 Revised: 17 February 2023 Accepted: 23 February 2023 Published: 27 February 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). diagnostics Article A Novel Strategy to Fit and Validate Physiological Models: A Case Study of a Cardiorespiratory Model for Simulation of Incremental Aerobic Exercise Carlos A. Sarmiento 1,* , Leidy Y. Serna 2,3 , Alher M. Hernández 1and Miguel Á. Mañanas 2,3 1Bioinstrumentation and Clinical Engineering Research Group, Bioengineering Department, Engineering Faculty, Universidad de Antioquia UdeA, Calle 70 # 52-51, Medellin 050016, Colombia 2Departament d’Enginyeria de Sistemes, Automàtica i Informàtica Industrial (ESAII), Universitat Politècnica de Catalunya, 08028 Barcelona, Spain 3CIBER de Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN), 28040 Madrid, Spain *Correspondence: [email protected] Abstract: Applying complex mathematical models of physiological systems is challenging due to the large number of parameters. Identifying these parameters through experimentation is difficult, and although procedures for fitting and validating models are reported, no integrated strategy exists. Additionally, the complexity of optimization is generally neglected when the number of experimental observations is restricted, obtaining multiple solutions or results without physiological justification. This work proposes a fitting and validation strategy for physiological models with many parameters under various populations, stimuli, and experimental conditions. A cardiorespiratory system model is used as a case study, and the strategy, model, computational implementation, and data analysis are described. Using optimized parameter values, model simulations are compared to those obtained using nominal values, with experimental data as a reference. Overall, a reduction in prediction error is achieved compared to that reported for model building. Furthermore, the behavior and accuracy of all the predictions in the steady state were improved. The results validate the fitted model and provide evidence of the proposed strategy’s usefulness. Keywords: fitting; mathematical modeling; sensitivity analysis; parameter estimation; cardiorespiratory system; aerobic exercise 1. Introduction Mathematical modeling is an interdisciplinary field that, applied to medicine, has allowed a better understanding of physiological functions and relationships. Its support in this field has been related to the education and training of clinical staff, the identification, monitoring, and treatment of diseases, and equipment development [ 1 – 5 ]. Mathematical models of physiological systems constitute approximations fitted to bounded populations and specific sets of physiological events or stimuli, so their practical and conceptual utility depends on their ability to predict experimental measurements of physiological variables [ 2 , 3 , 6 ]. Its mathematical complexity increases with accuracy and physiological relevance, so the number of parameters and equations is considerably high for describing the regulation mechanisms of the most complex physiological systems [7,8]. The cardiorespiratory system stands out as one of the most relevant physiological systems regarding the diagnosis, treatment, monitoring, and prevention of diseases, mainly due to the importance and usefulness of their related variables for the identification of the correct state and functioning of the organism [ 9 – 11 ]. The cardiorespiratory variables result from different and complex regulatory mechanisms focused on maintaining correct physiological functioning, even under different external conditions, such as diseases [ 12 , 13 ]. Different mathematical models have focused on the cardiorespiratory system, highlighting Diagnostics 2023,13, 908. https://doi.org/10.3390/diagnostics13050908 https://www.mdpi.com/journal/diagnostics Diagnostics 2023,13, 908 2 of 33 its support for personalized diagnosis applications due to its ability to be fitted to a subject or population under a specific stimulus [1,8,14]. Physical exercise is a natural stimulus that generates significant cardiorespiratory variations that cannot be evidenced under resting conditions and constitute a well-defined characteristic pattern in healthy human subjects [ 15 , 16 ]. Several mathematical models of the cardiorespiratory system have been proposed, but few of them can correctly represent the human response to exercise. The handicaps of those reported models to simulate the cardiorespiratory response to this stimulus are diverse. For instance, although some physiological models [ 8 , 17 ] detail some behaviors, they do not consider important mechanisms of neuronal regulation and omit the prediction of essential variables related to the exercise response. In contrast, other more complete models [ 18 , 19 ], despite including a complete description of the systems and allowing the prediction of variables useful for exercise analysis, are not adequate for predicting the human response under the stimulus due to the lack of related physiological mechanisms. For this aspect, models of the respiratory [ 20 , 21 ] and cardiovascular [ 15 ] controller, and even more recent works such as the model by [ 22 ], involve mechanisms and dynamics related to the response under aerobic exercise. However, the application of these models in a complex study is limited by their specificity and specialization in only one of the systems involved. Although mathematical models can describe the physiological response under a specific stimulus, their computational simulation using parameters’ nominal values constrains their use in populations with the same characteristics and conditions used during their development [ 14 , 23 ]. Fitting processes allow the personalization of physiological models to a specific subject or population by estimating the parameter values that minimize the difference between predictions and experimental observations [ 6 , 8 ]. Applying these techniques to physiological models could involve high complexity due to the large number of parameters that need to be identified. In addition, the restricted number of experimental observations could imply multiple solutions, and several parameter values could not have a physiological sense [8]. Different techniques to solve the fitting problem of models with many parameters have been proposed. They all focused on identifying the minimum number of parameters that must be fitted to correctly predict the experimental response. In [ 6 ], a deterministic sensitivity analysis focused on those parameters that are resolvable in the presence of noise, based on the effect of parameter value variations on the interest variables, was described. In a different approach, Ref. [ 24 ] validated a subset selection method that focused on finding the well-conditioned independent parameters for reliable identification using experimental measurements. Refs. [ 6 , 24 ] do not consider the evaluation of different levels of stimulus. In contrast, Ref. [ 8 ] proposed a classification technique based on the overall sensitivity of the model regarding parameter variations and considering different stimulus levels [ 23 ]. However, this proposal underestimates sensitivity measures by considering the sum of changes in different directions and assuming that all variables contribute equally to the error. Although the subsequent fitting could be designed for a steady state, its sensitivity computation includes time lags that are not consistent. The validation procedure aims to verify whether the model predictions agree with the experimental data for the defined conditions and disturbances. For this purpose, measures that consider the transient and steady-state response features of the analyzed physiological variables are usually implemented. Although there is no consensus on the proper methodology for validating mathematical models of physiological systems, different works coincide with the criteria proposed by [ 25 ]. Such criteria allow validating a model when the dynamic of its responses is consistent with the expected behavior and the values in the steady state are accurate. This work aims to propose a fitting and validation strategy for physiological models with many parameters that consider different populations and stimuli. This strategy includes several techniques and approaches to classify, select, and optimize model parameters under steady-state conditions. A model validation methodology is also presented. Diagnostics 2023,13, 908 3 of 33 Some described methods are new proposals or correspond to improvements to previously reported methods. The strategy application in a cardiorespiratory model during dynamic aerobic exercise is presented and described as a case study. This paper is structured as follows. First, the fitting and validation strategies are presented. It describes in detail the procedures for the classification, selection, and optimization of parameters and methodology validation. Second, a case study is presented. It includes a qualitative description of the cardiorespiratory model and the simulation, fitting, and validation. Third, predictions at different levels of aerobic exercise, both with the fitted model and the nominal values, are compared with experimental data obtained from healthy subjects under a cardiopulmonary exercise test [ 26 ]. Finally, the proposed strategy and the model performance results are discussed. 2. Materials and Methods 2.1. Fitting and Validation Strategy The strategy proposed in this work consists of applying three main procedures or steps that must be done sequentially: (a) classification and selection of parameters, (b) model fitting, and (c) model validation. The first procedure involves classifying the model parameters according to some predefined roles in the model (gain, threshold, initial conditions, etc.) and selecting the most relevant ones according to different model-fitting approaches. The second one consists of four sequential stages of parameter value identification related to the model fitting regarding the available experimental information, the overall accuracy, the specific prediction of each variable, and the predictions’ behavior concerning the evaluated stimulus. The final procedure focuses on validating the fitted model regarding its predictions’ accuracy, behavior, and transient and steady-state regimes. Each procedure is presented in Sections 2.1.1–2.1.3, and its application is described according to the case study. Figure 1shows a schematic summary of the strategy, in which the main procedures and their respective stages are shown. Diagnostics 2023, 13, x FOR PEER REVIEW 3 of 35 This work aims to propose a fitting and validation strategy for physiological models with many parameters that consider different populations and stimuli. This strategy includes several techniques and approaches to classify, select, and optimize model parameters under steady-state conditions. A model validation methodology is also presented. Some described methods are new proposals or correspond to improvements to previously reported methods. The strategy application in a cardiorespiratory model during dynamic aerobic exercise is presented and described as a case study. This paper is structured as follows. First, the fitting and validation strategies are presented. It describes in detail the procedures for the classification, selection, and optimization of parameters and methodology validation. Second, a case study is presented. It includes a qualitative description of the cardiorespiratory model and the simulation, fitting, and validation. Third, predictions at different levels of aerobic exercise, both with the fitted model and the nominal values, are compared with experimental data obtained from healthy subjects under a cardiopulmonary exercise test [26]. Finally, the proposed strategy and the model performance results are discussed. 2. Materials and Methods 2.1. Fitting and Validation Strategy The strategy proposed in this work consists of applying three main procedures or steps that must be done sequentially: (a) classification and selection of parameters, (b) model fitting, and (c) model validation. The first procedure involves classifying the model parameters according to some predefined roles in the model (gain, threshold, initial conditions, etc.) and selecting the most relevant ones according to different model-fitting approaches. The second one consists of four sequential stages of parameter value identification related to the model fitting regarding the available experimental information, the overall accuracy, the specific prediction of each variable, and the predictions’ behavior concerning the evaluated stimulus. The final procedure focuses on validating the fitted model regarding its predictions’ accuracy, behavior, and transient and steady-state regimes. Each procedure is presented in Sections 2.1.1–2.1.3, and its application is described according to the case study. Figure 1 shows a schematic summary of the strategy, in which the main procedures and their respective stages are shown. Figure 1. Schematic diagram of the fitting and validation strategy. The dotted lines divide the diagram into the main procedures and their respective stages: the first procedure is at the top and corresponds to the classification and selection of parameters; the second is in the middle and corresponds to model fitting, and the third one is at the bottom and corresponds to model validation. Figure 1. Schematic diagram of the fitting and validation strategy. The dotted lines divide the diagram into the main procedures and their respective stages: the first procedure is at the top and corresponds to the classification and selection of parameters; the second is in the middle and corresponds to model fitting, and the third one is at the bottom and corresponds to model validation. 2.1.1. Classification and Selection of Parameters This procedure aims to highlight and select the parameters that can be considered best-fit candidates. The selection of parameters is justified according to their relationship to available experimental data, role in the model, identifiability, sensitivity to variations, and Diagnostics 2023,13, 908 4 of 33 their relationship to the stimulus evaluated. Initially, it reduces the number of parameters depending on their role in the model and the objective of fitting predictions in steady-state conditions. Subsequently, the resulting parameters are classified and selected according to two reported techniques and established criteria regarding the model fitting approaches. Classification and Selection of Parameters by Role It consists of an initial classification of the model parameters according to their roles. It comprises five parameter classes generally found in structured models of physiological systems. They correspond to (i) time constants, i.e., parameters related to transient response; (ii) conversion parameters, which are constant values related to equivalences among measurement units; (iii) covariates, corresponding to values that allow defining simulation conditions regarding external disturbances, environment conditions, and features of the population to be simulated; (iv) initial values, corresponding to initial conditions of the model variables, usually required as initial states of integrators and whose action mainly affect the temporal characteristics of the responses before reaching the model steady state; and (v) gain and thresholds, that either module or saturate variables related to model mechanisms, i.e., the weighting of the chemoreceptors response to set the parasympathetic and sympathetic activity regarding the regulation of peripheral resistances in cardiovascular control models. The selection of parameters comprises choosing those parameters that mainly define the steady-state behavior of the cardiorespiratory system response and, therefore, can be used to fit the model response in such a condition. The time constant parameters and initial values are discarded, considering that the proposed fitting procedure focuses on the response once the steady state is reached. Conversion parameters do not need fitting because of their nature and meaning. Therefore, only gain and threshold parameters, whose variations have physiological sense, and are correlated with experimental conditions, and covariates will be considered for the subsequent selection procedures and fitting process. Covariates are only used for the standardization of the simulation conditions. Parameter Classification Techniques The proposed strategy involves modifying and implementing two of the most widely implemented parameter classification techniques for fitting physiological models [ 8 , 23 ]. A detailed description of each is presented below. Subset Selection of Parameters It is a technique based on the classification of model parameters according to how well-conditioned or ill-conditioned they can be identified [ 27 ]. Well-conditioned parameters correspond to those that can be reliably estimated from the constrained experimental data, while ill-conditioned ones are those for which there are multiple fitting solutions. In this study, the subset selection of parameters is based on QR Factorization with the Column Pivoting method [ 28 ]. It was selected because it has been widely implemented in different physiological models [ 2 , 8 , 23 , 29 ], presented one of the best results regarding cardiorespiratory models among the other methods reported [ 28 ], and considers the difference between model predictions and experimental data as a reference [ 2 ]. This technique is based on the solution of optimization problems using gradient techniques. It establishes a ranking of the well-conditioned parameters for identification by analyzing the interdependencies in the Jacobian. This technique is complemented by integrating the different stimulus exercise levels (according to the experimental data) and variations in the parameters previously selected in the appropriate physiological range (around their nominal values). These additions would allow the selection of parameters that are more consistent with the desired fit. According to the above, the Jacobian is calculated according to Equation (1). r0 klujl= ∂riklujl ∂ujl , (1) Diagnostics 2023,13, 908 5 of 33 where the indexes i, j, and k represent the variable, parameter, and stimulus level to be evaluated, respectively; l is an index that identifies the parameter variation regarding its nominal value. r0 kl is the Jacobian for a specific stimulus level and parameter value variation, formed as the matrix of variations of the differences for each variable rikl regarding the change of each parameter for a specific variation level ujl ; r corresponds to the difference between each experimental data and the model prediction given each parameter’s change for a particular level of variation. The singular value decomposition is calculated for r0 kl , according to Equation (2). r0 klujl=UjklΣjklVT jkl, (2) where Ujkl is the matrix of left singular vectors, Σjkl is the diagonal matrix of singular values of r0 klujl in decreasing order, Vjkl is the matrix of the right singular vectors, and T denotes the matrix transpose. Matrix Vjkl must be partitioned, as expressed in Equation (3). Vjkl =hVjkl,ρ(j,k,l)Vjkl,W−ρ(j,k,l)i, (3) where W is the total number of parameters analyzed, and ρ(j,k,l) is a numerical rank that indicates the number of maximally independent columns of r0 klujl . ρ(j,k,l) is equivalent to the number of parameters that can be identified given the model output and can be determined by the selection of the smallest allowed singular value according to the relation expressed in Equation (4) [23]. σjkl,ρ(j,k,l) σjkl,1 =σNjkl,ρ(j,k,l)>ε, (4) where σjkl,ρ(j,k,l) is the singular value for ρ(j,k,l) , σjkl,1 is the largest singular value, σNjkl,ρ(j,k,l) is the normalized singular value for ρ(j,k,l) , and ε is a tolerance value that allows differentiating the most significant eigenvalues. The parameters associated with the ρ(j,k,l) highest singular values are found using QR decomposition with column pivoting, according to Equation (5). VT jkl,ρ(j,k,l)Pjkl =QjklRjkl, (5) where Pjkl is a permutation matrix, Qjkl is an orthogonal matrix and the first ρ(j,k,l) columns of Rjkl form an upper triangular matrix with diagonal elements in decreasing order. Pjkl is then used to reorder the parameters according to Equation (6). ˆ ujkl =PT jkluI, (6) where uI is an identification vector of the parameters and ˆ ujkl is the vector of the parameters reordered, which is partitioned regarding ρ(j,k,l), as presented in Equation (7). ˆ ujkl =hˆ ujkl,ρ(j,k,l)ˆ ujkl,W−ρ(j,k,l)i, (7) where ˆ ujkl,ρ(jk,l) is a vector containing the ρ(j,k,l) estimable parameters and ˆ ujkl,W−ρ(k,l) are the parameters that would be fixed at nominal values. The standard 2-norm of the normalized singular value for each parameter in each ranking obtained for each parameter variation must be determined to obtain a general ranking independent of the stimulus level, as presented in Equation (8). Zjl =kZjlkk2=r1 K∑K k=1Zjkl2, (8) Diagnostics 2023,13, 908 6 of 33 where Zjl corresponds to the standard 2-norm of the normalized singular value of the parameter j for the variation level l, Zjkl is the normalized singular value of the parameter j for the stimulus level k and the parameter variation level l, and K is the number of stimulus levels. Finally, the result of this technique corresponds to the presented in Equation (9). Zj=kZjlk2=r1 L∑L l=1Zjl2, (9) where Zj is the standard 2-norm of Zjl of the parameter j, and L is the total number of parameter variations. This technique is usually implemented with sensitivity analysis techniques because their approaches complement each other regarding the problem solution. QR decomposition identifies the parameters to which the model is sensitive as a group, while sensitivity analysis finds the parameters to which the model is individually sensitive [2]. Sensitivity Analysis Sensitivity analysis involves the integration and improvement of the techniques introduced by [ 6 , 8 , 23 ]. It comprises a deterministic analysis that evaluates the global model and variables’ sensitivity to variations in model parameters and stimulus levels. It also involves an experimental data dependency term that weighs such sensitivity by the error reached by the model at each parameter variation and stimulus level. The model variable’s sensitivity regarding each parameter variation is based on the standard local differential equation described by [ 30 ] and used by [ 8 , 19 , 23 , 31 , 32 ] to calculate time-dependent sensitivities in cardiovascular and respiratory models. Equation (10) shows the computation of this sensitivity. sij(t, u)u=un=∂Yit,uj ∂uj uj Yi(t,u)u=un ;uj,Yi(t,u)6=0 , (10) where sij represents the relative sensitivity of the variable Yi to parameter uj , which is dimensionless by the ratio between the parameter uj and the variable Yi values at nominal conditions (no parameter variation). n refers to the parameter’s nominal value. This work proposes modifying the measurement of relative sensitivity to make it independent of the time and dependent on the stimulus. Therefore, an approach based on steady-state conditions at different stimulus levels is used. Relative sensitivity is evaluated by varying each parameter over a range around its nominal value, while the others are kept at their nominal values. Equation (11) shows the proposed relative sensitivity measure. Time independence avoids significant differences between sensitivity measures due to time lags, considering that the fitting model procedure is focused on minimizing the steadystate differences between the experimental measurements and the model predictions. The evaluation of the stimulus is included by considering the influence of its variations on sensitivity measures. sijk = ∂Yikuj ∂uj ujn Yik (un) , (11) where, i, j, and k are indexes that refer to the analyzed variable, parameter, and stimulus level. sijk is a scalar value of the relative sensitivity; Y is the variable value in the steady state; u is the parameter value; n refers to the parameter’s nominal value. Therefore, Equation (11) evaluates changes in the variable Yi regarding variations in the parameter uj at the stimulus level k (first quotient), which is dimensionless by the ratio between the parameter’s nominal values ujn and the variable Yiat the stimulus level k (second quotient). Deriving sensitivity equations can be tedious and error-prone for large systems, mainly when they involve nonlinear features, such as those analyzed in this case study. Alternatively, Equation (11) can be solved using a computational approach consisting of a simple Diagnostics 2023,13, 908 7 of 33 finite difference method. A numerical approximation of the derivatives, which also considers the variation of the parameter relative to its nominal value, is expressed in Equation (12). sijlk ≈ Yikujn +hlujn−Yik(un) ujn +hlujn−ujn ujn Yik (un) , (12) where sijkl is a scalar value of the sensitivity for the variable Yi regarding the variation hl in the parameter uj at the stimulus level k . h is the vector of change proportions of the parameter regarding its nominal value. This expression can be reduced and organized as expressed in Equation (13). sijlk ≈ Yik(hl+1)ujn−Yik(un) Yik (un) 1 hl =DYijlk 1 hl , (13) where the first quotient of the equation, also identified as DYijkl , is a dimensionless term representing the rate of change of the variable Yi when a percentual variation hl in the parameter uj is applied at the stimulus level k. The second quotient is a dimensionless term representing a weighting factor, giving heavier importance to slight parameter variations. Therefore, the sensitivity sijlk corresponds to a dimensionless scalar value that measures the weighted relative variations of the analyzed variables according to four degrees of freedom. Standard 2-norm is applied to sijlk to obtain a measure of the relative sensitivity independent of the parameter variations, as expressed in Equation (14). As a result, a positive dimensionless scalar value related to the relative sensitivity’s mean trend is obtained for all parameters. sijk =ksijlkk2=s1 L∑L l=1DYijlk 1 hl2 , (14) where L is the total number of parameter variations or elements of the vector h. The stimulus level independency relative sensitivity is calculated according to Equation (15). sij = sP,ijk max(sik)2 =s1 K∑K k=1Pik·sijk max(sik)2 , (15) where K is the total number of stimulus levels; max(sik) is the maximum relative sensitivity obtaining for the variable i at the stimulus level k considering all parameters; sP,ijk is the weighted relative sensitivity; and Pik is a weighting factor representing the variable error at the stimulus level k. Sij corresponds to the standard 2-norm of sijk normalized by its maximum value at each stimulus level k and weighted by Pik . Normalization of sijk allows comparing the sensitivity measures between different stimulus levels, whereas the inclusion of the weighting factor Pik prioritizes the sensitivity obtained at the stimulus levels for which the variable’s error is more significant (i.e., errors obtained from predictions resulting from using parameter’s nominal values). As a result, a ranking of the model parameters concerning their variation impact for each variable is obtained. Pik is calculated according to Equation (16). Pik =Eik K·EiT , (16) where Eik is the error of the variable Yi at the stimulus level k, and EiT is the total error for the variable Yi (considering all stimulus levels). Equations (17) and (18) are proposed to calculate the mentioned errors. Diagnostics 2023,13, 908 8 of 33 Eik = Yexp,ik −Yik(un) Yexp,ik  , (17) EiT =1 K K ∑ k=1 Eik, (18) where Yexp,ik represent the experimental value of the variable Yiat the stimulus level k. The model’s total sensitivity to each parameter is calculated according to Equation (19). sj= sP,ij max(si)2 =s1 I∑I i=1Pi·sij max(si)2 , (19) where I is the total number of variables evaluated; max(si) is the maximum value of sensitivity among all the parameters for the variable Yi ; sP,ij is the weighted sensitivity of each parameter ujfor the variable Yi; Piis the weighting factor related to the variable and it is calculated according to the Equation (20). Pi=Ei I·ET, (20) where Ei is the error of the variable i, and ET is the total model error. Equations (21) and (22) are proposed to calculate the mentioned errors. Ei=v u u t1 K∑K k=1 Yexp,ik −Yik(un) Yexp,ik !2 , (21) ET=1 I K ∑ k=1 Ei, (22) Sj corresponds to the standard 2-norm of Sij normalized by the maximum sensitivity obtained in each variable and weighted by Pi . Normalization of sij allows comparing the sensitivity measures of the parameters among different variables, whereas Pi prioritizes the sensitivities of the parameters for which the variable’s error is more significant. As a result , a ranking of the model parameters representing the effect of their variation on the model’s whole output is obtained. Parameter Selection Criteria This procedure focuses on applying the selection criteria of the classified parameters concerning four different fitting approaches. Selection for the Standardization of Simulation Conditions It consists of selecting the model’s parameters, which can be determined from the available experimental data. It involves those related to the subjects’ characteristics, the stimulus evaluated, or the environmental conditions. They can be established by direct equivalence or by applying previously validated equations. Selection for the Base Fitting Approach It comprises selecting the set of parameters for which the model has the highest global sensitivity. This selection aims to reduce the model’s total prediction error by considering the parameter’s overall effect on model behavior. They correspond to the union of the set of parameters obtained in subset selection (i.e., those parameters that have been identified as well-conditioned to be adjusted, Equation (9)) and total sensitivity techniques ( i.e., those parameters whose variations generate significant changes in the model output, Equation (19) ). In this work, this selection is defined according to the following criteria. Diagnostics 2023,13, 908 9 of 33 1. Parameters from the subset selection ranking are chosen according to the criterion defined in Equation (4), considering ε as the square root of the termination tolerance on the function evaluation defined for the fitting procedure (Table 1). 2. Parameters from the total sensitivity ranking are chosen in descending order until at least one is obtained for each system and controller of the model, including those already selected in the previous step. Table 1. Parameter values used for the CMA-ES algorithm. Name Definition Value Fitness limit Value to reach Infinite TolFun Termination tolerance in the function evaluation 10−12 TolX Termination tolerance on x 10−3 MaxIter Maximal number of iterations 100 ×N2* MaxFunEval Maximal number of function evaluations 500 MaxRestart Number of restarts 10 PopSize Population size 4+3ln(N)2 σInitial coordinate-wise standard deviation(s) 0.2(UB −LB)** * N is defined by the number of the objective variables. ** UB and LB are the upper and lower bounds of objective variables (related to the model parameters’ evaluation ranges). Selection for a Specific Fitting Approach It corresponds to the set of parameters for each variable’s highest sensitivity at the individual level. Its optimization aims to modify each variable’s predictability to bring it closer to the respective experimental data without significantly affecting the other predictions. It is based on relative sensitivity measures independent of the stimulus level (Equation (15)). Based on the following criteria, only one parameter is selected by each ranking obtained for the evaluated physiological variables. 1. Remove the parameters selected for the base fitting approach from each variable ranking . 2. Remove the parameters of systems and controllers that are not directly related to regulating the variable of interest. 3. Select only those parameters whose sensitivity is high for the variable of interest and low for the remaining ones. Parameters with high sensitivity for other variables could be selected if those variables belong to or are dependent on the same system or controller of the variable of interest. Selection for the Stimulus-Related Fitting Approach This corresponds to the parameter selection criteria that relate the stimulus to the regulation mechanisms addressed in the model under study. Their selection is based on the parameters’ role concerning the mechanisms mainly associated with the stimulus, highlighting the weighting factors or gains that link them to regulatory activities. Eight mechanisms directly related to the exercise stimulus were evaluated in this work. Only one parameter per mechanism was selected. 2.1.2. Model Fitting Fitting a parametric model involves solving an optimization problem in which the values of the set of parameters minimizing the differences between experimental data and model predictions are identified [ 6 ]. The identification of the parameter values in this strategy results in applying an optimization algorithm in three stages that must be carried out sequentially: first, a base optimization; second, a specific optimization; and third, a stimulus-related optimization. Each stage focuses on the value identification of a specific and reduced number of parameters. It is proposed to apply the following procedures before the mentioned optimization stages to obtain correct, fast, and physiological meaning results: the standardization of simulation conditions, the selection and parameterization of the optimization algorithm, the definition of the parameter evaluation ranges, and the choice Diagnostics 2023,13, 908 16 of 33 Diagnostics 2023, 13, x FOR PEER REVIEW 16 of 35 Figure 2. Top ten positions in the rankings of classification techniques: (a) correspond to the results for subset selection and (b) to the results for sensitivity analysis techniques. The light gray bars correspond to the selected base setting parameters, and the dark gray bars correspond to the parameters that will remain fixed at the nominal values. Eight base fitting parameters were selected from the rankings, three from the subset selection technique under consideration of the tolerance range of 1×10 (Table 1 and Equation (4)), and five from the total sensitivity analysis corresponding to those whose sensitivity was greater than 46% of the maximum sensitivity found (Section 2.1.1). Regarding the subset-selection results, it is highlighted that there is a single significant difference between the parameter in the first position and the others presented, evidencing the difficulty of identifying a single optimal solution for more than one parameter using gradient-based techniques [27]. Regarding the total sensitivity analysis, all the selected parameters directly affect the model’s controllers, which are the ones that have the most significant influence on the model’s base behavior, as expected. This approach evidences the impact of the modifications on selection methods, mainly those related to the error and stimulus-level evaluation. In this sense, the distribution of the parameters regarding the systems and controllers (5 parameters of the respiratory controller, 2 of the cardiovascular system and controller, and 1 of the gas exchange system) is associated with the error contribution of the related variables [26]. After applying the selection criteria, the number of parameters considered for this fitting procedure was reduced by 97.4%. The selected parameters, the corresponding model equations and mechanisms, and the possible variation range information are presented below. 𝐾𝑅 is a parameter that describes the dependence of the left ventricle resistance (R) on the isometric pressure (P,) in the cardiovascular system [47]. Its nominal value results from scaling data extracted from animal experiments to reflect changes in ventricle volume in human beings. It has been used in different works without presenting variations [15,18,19,22,47,48], even in studies focusing on personalized cardiovascular models [35]. Considering the simulation results under nominal conditions of R and P,, and applying the rescaling approach of vascular resistance implemented in [35] concerning the mean total blood volume from experimental data, variations close to 6% of the nominal value could be expected. Figure 2. Top ten positions in the rankings of classification techniques: ( a ) correspond to the results for subset selection and ( b ) to the results for sensitivity analysis techniques. The light gray bars correspond to the selected base setting parameters, and the dark gray bars correspond to the parameters that will remain fixed at the nominal values. This approach evidences the impact of the modifications on selection methods, mainly those related to the error and stimulus-level evaluation. In this sense, the distribution of the parameters regarding the systems and controllers (5 parameters of the respiratory controller, 2 of the cardiovascular system and controller, and 1 of the gas exchange system) is associated with the error contribution of the related variables [26]. After applying the selection criteria, the number of parameters considered for this fitting procedure was reduced by 97.4%. The selected parameters, the corresponding model equations and mechanisms, and the possible variation range information are presented below. KRlv is a parameter that describes the dependence of the left ventricle resistance ( Rlv ) on the isometric pressure ( Pmax,lv ) in the cardiovascular system [ 47 ]. Its nominal value results from scaling data extracted from animal experiments to reflect changes in ventricle volume in human beings. It has been used in different works without presenting variations [ 15 , 18 , 19 , 22 , 47 , 48 ], even in studies focusing on personalized cardiovascular models [ 35 ]. Considering the simulation results under nominal conditions of Rlv and Pmax,lv , and applying the rescaling approach of vascular resistance implemented in [ 35 ] concerning the mean total blood volume from experimental data, variations close to 6% of the nominal value could be expected. MRBCO2 denotes the metabolic production rate for CO2 in the brain tissue [ 21 ] and allows us to relate PaCO2 with PvbCO2 , the CO2 brain venous blood pressure. Its reported nominal value is 0.0009 L/s STPD and has not been fitted in the different validation works in which the associated model has been used [ 14 , 21 , 22 , 45 ]. Works related to other validated models report values for parameters with the same physiological sense ranging from 0.0007 to 0.00104 L/s STPD [ 8 , 19 , 49 – 51 ]. Therefore, fitted values between − 22.2% and 15.6% of the nominal value could be expected. T0 is a cardiovascular controller parameter representing the heart period (HP) in the absence of cardiac innervation [ 48 ]. It is the offset term in the equation that relates to the changes in HP induced by sympathetic and parasympathetic stimulation ( ∆Ts and Diagnostics 2023,13, 908 17 of 33 ∆Tv , respectively). Its reported nominal value results from a fitting process to mimic animal experimental data [ 47 , 52 ]. It has been used without any modification in different publications related to cardiovascular control models [ 15 , 18 , 19 , 22 , 47 , 48 ]. Considering this value as a proportion of the reported HP at rest (0.833 s), variations close to –7% of the nominal value could be expected whether the median experimental value of HP at rest (0.775 s) is considered a reference. Pmax is a parameter of the respiratory controller that denotes the maximum inspiratory pressure [ 20 ]. It relates to the inspiratory muscles’ capacity to minimize the work of breathing. In [ 20 ], a nominal value of 150 cmH2O is proposed, and variations of around ± 66% of this value are evaluated. Subsequent work has used a value of 50 cmH2O on fitting and validating such a model with healthy subjects [ 53 ]. Following the above, although the optimization procedure can identify a value near the nominal value, variations higher than 100% could be expected. Ers is a parameter that denotes the overall elastance of the ventilatory system. It is used in the model’s respiratory controller and lung mechanics to represent the motion equation of the respiratory system. Its nominal value is 21.9 cmH2O /l and agrees with that reported in [ 20 , 53 ] for healthy adult subjects. However, other authors report lower values of 10, 8.55, and 8.52 cmH2O/l . According to the above, variations of around –66% of the nominal value could be expected [14,21,50,54]. KpO2 , KcCO2 , and Kbg are the parameters of the respiratory controller associated with the control of ventilation. KpO2 and KcCO2 correspond to constant weighting factors related to the peripheral chemoreceptors for O2 and the central chemoreceptors for CO2 contributions, respectively. Kbg is an offset term that relates to the blood gas dissociation constant for the controller [ 26 ]. These parameters correspond to fitted values to mimic the change in alveolar ventilation from rest. They have not been modified in studies in which experimental data on healthy adult subjects are used [14,21,22,55]. 3.1.4. Parameters Selected for the Specific Fitting Approach Figure 3shows each variable’s relative sensitivity ranking results after removing the parameters selected for the base fitting approach. They are presented in descending order from left to right. The parameters selected for the specific fitting approach are highlighted in the ranking of each variable. The results highlighted that not all the parameters selected for each variable corresponded to the first ranking positions, which shows the significant influence of the regulation results between the model’s different systems and controllers. This fact was presented mainly for variables not directly related to the controllers (PM, PD, and PACO2 ), evidencing that the parameters of the main regulatory mechanisms have the most significant influence on the model. The parameter selected for PACO2 , corresponding to VLCO2 , is not even among the positions presented in the respective ranking, showing the high sensitivity of the variable to the results of other systems and controllers. The number of parameters considered for this fitting procedure was reduced by 97.2%. The selected parameters, the related model equations and mechanisms, and the information reported regarding their possible variation ranges are presented below. GT,v corresponds to the weighting factor that relates parasympathetic activity with heart rate regulation in the cardiovascular controller [ 48 ]. This parameter was initially adjusted to mimic the humans’ cardiac period’s response in [ 47 ]. Its nominal value is not related to direct physiological measurement and has not been modified in subsequent applications of the model [15,18,19]. Rsa is a parameter of the cardiovascular system that represents systemic arterial hydraulic resistance [ 48 ] and relates Psa to systemic arterial blood flow ( Qsa ). Its nominal value was initially computed according to experimental cardiac output and cardiovascular pressure measurements [ 56 ]. It has been modified in different cardiovascular system versions due to the inclusion of new vascular beds [ 47 ]. The latest reported nominal value corresponds to that defined for the case study model [ 26 ]. Although it has not been Diagnostics 2023,13, 908 18 of 33 modified in other validation works [ 15 , 19 , 48 , 57 ], a rescaling approach was proposed in [ 35 ] for personalized fitting based on experimental data on total blood volume. Following the above, using the average experimental value of total blood volume reported in [ 26 ] as a reference, a variation close to 5.42% of the nominal value could be expected. Diagnostics 2023, 13, x FOR PEER REVIEW 18 of 35 Figure 3. Relative sensitivity rankings after removing the parameters selected for the base fitting approach. The light gray bars correspond to the selected parameters for the specific fitting approach, and the dark gray bars correspond to the parameters that remain fixed at nominal values. The subfigures from (a) to (i) correspond to the result for each of the variables evaluated. The results highlighted that not all the parameters selected for each variable corresponded to the first ranking positions, which shows the significant influence of the regulation results between the model’s different systems and controllers. This fact was presented mainly for variables not directly related to the controllers (PM, PD, and PACO), evidencing that the parameters of the main regulatory mechanisms have the most significant influence on the model. The parameter selected for PACO, corresponding to VLCO, is not even among the positions presented in the respective ranking, showing the high sensitivity of the variable to the results of other systems and controllers. The number of parameters considered for this fitting procedure was reduced by 97.2%. The selected parameters, the related model equations and mechanisms, and the information reported regarding their possible variation ranges are presented below. 𝐺, corresponds to the weighting factor that relates parasympathetic activity with heart rate regulation in the cardiovascular controller [48]. This parameter was initially adjusted to mimic the humans’ cardiac period’s response in [47]. Its nominal value is not related to direct physiological measurement and has not been modified in subsequent applications of the model [15,18,19]. 𝑅 is a parameter of the cardiovascular system that represents systemic arterial hydraulic resistance [48] and relates P to systemic arterial blood flow (Q). Its nominal value was initially computed according to experimental cardiac output and cardiovascular pressure measurements [56]. It has been modified in different cardiovascular system versions due to the inclusion of new vascular beds [47]. The latest reported nominal value Figure 3. Relative sensitivity rankings after removing the parameters selected for the base fitting approach. The light gray bars correspond to the selected parameters for the specific fitting approach, and the dark gray bars correspond to the parameters that remain fixed at nominal values. The subfigures from (a) to (i) correspond to the result for each of the variables evaluated. KE,lv describes the left ventricle’s function at the end of the diastole based on the pressure–volume relationship [ 47 ]. Its nominal value was initially identified to fit the exponential relationship proposed in [ 52 ] to healthy humans’ experimental measurements. According to the experimental measurements reported in [ 52 , 58 , 59 ], variations of less than − 35.7% and greater than 7.1% of the nominal value could be associated with symptoms of cardiac pathologies. Φmax is a cardiovascular controller parameter associated with the vasodilation of the peripheral resistance of the active muscles ( Ramp ) during exercise [ 26 ]. It corresponds to the upper saturation bound of a sigmoid function used to describe the vasodilator effect independent of tissue hypoxia during exercise [ 15 , 60 ]. The parameter’s nominal value is not related to a direct physiological measure but results from optimization procedures to imitate experimental data. λ1 is a parameter related to the breathing pattern optimizer and corresponds to a weighting factor that relates the mechanical work of the inspiratory phase ( . WI ) with the average square magnitude of volume acceleration [ 26 ]. Its value was fitted in [ 14 , 33 ] to the experimental data of healthy subjects. Considering their results implies possible variations in the nominal value of around −43%. Diagnostics 2023,13, 908 19 of 33 V 0 dead is a parameter of the respiratory controller related to the regulation of . VE [ 26 ]. This corresponds to the offset term of the empirical equation used in [ 21 ] to calculate the dead space volume as a function of alveolar ventilation. It has not been modified in subsequent validation studies on healthy subjects [14,22]. n is a parameter related to the breathing pattern optimizer, used as a power index of efficiency factors that relate . WI with Pmusc(t) [ 26 ]. Following the values fitted in the validation work by [ 20 , 33 ], a variation ranging from − 9.17% to 81.7% of the nominal value could be expected. C1 is a parameter related to the dissociation of oxygen in the blood and denotes the maximum concentration of hemoglobin-bound oxygen [ 50 ]. Its nominal value was taken from [ 61 ] and calculated from predefined pressure, temperature, and the amount of hemoglobin conditions. Although this value has not been modified in subsequent studies of the same blood gas dissociation model, other values are reported for this physiological measure [62], allowing a variation of approximately −3.6% of the nominal value. VLCO2 is a parameter that relates PACO2 with the blood concentrations of CO2 and Qpp [ 26 ]. It denotes the lungs’ storage volume for CO2 , and can be understood as a fraction of the functional alveolar volume of the lung [ 22 , 50 , 51 , 63 ]. Therefore, using the alveolar volume values reported in [ 63 ], a study based on subjects under exercise as a reference, variations around 50% of the nominal value could be expected. 3.1.5. Stimulus-Related Fitting Parameters For this approach, the selected parameters relate to each of the eight mechanisms reported for the case study model [ 26 ]. Only one parameter for each mechanism’s regulatory activity, which was not selected in the previous fitting approaches, was selected. The mechanisms evaluated were: the central command action (I-EP) and the central respiratory neuromuscular drive (NT) on regulation control activities; the central vasodilatory action on active muscles due to central command (IRamp ); the independent description of venous vascular beds from active muscles (VRamv ); the muscle (MP) and the respiratory (RP) pumps; the neural driving ventilation related to metabolism (MRV); and the respiratory control action based on mechanical work of breathing minimizing (minWOB). I-EP, IRamp , VRamv , MP and RP are mechanisms that relate I (exercise intensity) to the regulatory activities of the cardiovascular controller and system. Their parameters were defined and optimized from human and animal experimental data [ 15 ]. No parameter was selected concerning blood flow signals in the gas exchange system because it is not considered a mechanism that explicitly relates to the exercise stimulus. The selected parameters, their corresponding mechanisms, and regulatory activity are presented in Table 3. The number of parameters considered for this fitting procedure was reduced by 95.3% with respect to the total after applying the selection criteria. The related mechanisms and the information reported on the possible variation range of each selected parameter are presented below. γsh,max , γsp,max , γsv,max and γv,max are the upper saturation values of sigmoid functions that relate I to sympathetic and vagal efferent activity. In IRamp , gM is a static gain that relates I to the effect of tissue hypoxia concerning the regulation of the active skeletal muscle’s peripheral resistance ( Ramp) . In VRamv , kr,am is a constant parameter that characterizes the inversely proportional behavior of the active muscles’ venous vascular beds ( Ramv ) concerning the total volume of blood it contains ( VTamv ) during exercise. In MP, Aim is a parameter that denotes the peak value of Pim , which affects the vascular venous pressure during exercise ( Pamv ). Finally, in RP, gabd and gthor are constant gain factors that relate changes in the tidal volume with the maximum and minimum values of Pabd and Pthor, respectively. Wt,sh , Wt,sp , Wt,sv and Wt,v are weighting factors that relate Nt to sympathetic and parasympathetic efferent activity. Their values were optimized in [ 50 ] to ensure that the dynamic behavior of the model under various conditions remains realistic, and they have been used in later works without any modification [14,22]. Diagnostics 2023,13, 908 20 of 33 Table 3. Selected parameters for stimulus-related fitting. Mechanism Parameter Regulatory Activity I-EP γsh,max Sympathetic activity to heart γsp,max Sympathetic activity to peripheral resistances γsv,max Sympathetic activity to veins volumes γv,max Vagal activity NT Wt,sh Sympathetic activity to heart Wt,sp Sympathetic activity to peripheral resistances Wt,sv Sympathetic activity to veins volumes Wt,vVagal activity I−Ramp gMEffect of hypoxia on vascular vasodilation of active muscles during exercise V−Ramv kr,am Changes in venous resistance of active muscles during exercise MP Aim Venous return of active muscles RP gabd Abdominal pressure signal gthor Thoracic pressure signal MRV KcMRV Change in alveolar flow from its resting value minWOB λ2Total work of breathing KcMRV is a parameter of the respiratory controller associated with ventilation. It corresponds to a constant gain that relates MRV to . VA . It was initially defined as equal to 1 in [ 21 ] under the consideration of a direct action of MRV. However, in the case study model, it was defined as a constant parameter to consider a weighting factor for the related mechanism. λ2 is a factor that weighs the expiratory mechanical work ( . WE ) in the total mechanical work of breathing ( . WT ). According to [ 14 , 20 , 33 ], a variation ranging between − 28.4% and 170% of its nominal value could be considered. 3.2. Model Fitting A general evaluation range of ± 30% regarding the nominal value was defined for those selected parameters whose variation could not be established or constrained to reported values in the literature. For that, the following considerations were taken into account: (a) it is suitable for the constrictions of each associated mechanism of the model, and (b) it agrees with the expected closeness of the results for the nominal values. It was also considered that these parameter values are related to subjects with physiological characteristics similar to the subject’s experimental data (healthy adult males). The evaluation ranges for the other parameters are presented in Table 4, following what was previously described in the selection parameter results. Table 4. Evaluation ranges of parameters according to the variations and constraints reported. Parameter Lower Bound (%) Upper Bound (%) Ers −70 30 KE,lv −30 5 λ1−50 30 λ2−30 200 n−30 90 Pmax −30 200 VLCO2−30 50 The values are shown for upper and lower bounds correspond to the percentage of variation regarding the parameter’s nominal value. It should be noted that most of the selected parameters were adapted to the defined general evaluation range; the remaining ones were mainly related to the increase of either lower or upper bounds. Only KE,lv required a decrease in the range due to its relationship with cardiovascular diseases. Diagnostics 2023,13, 908 21 of 33 Table 5compares the nominal parameter values against the best optimization results for each fitting stage. Most optimizations were between ± 15% of the nominal value, confirming the expected closeness, considering that the nominal values and the experimental data are related to subjects with similar physiological characteristics. The most significant changes that occurred in the second and third fitting stages related to the specific and stimulus-related fitting approaches. The following main modifications were obtained regarding the specific fitting approach: (a) a decreased effect of volume on left ventricular pressure due to the decrease in KE,lv ; (b) an increase in VD , and therefore an increase in . VE not related to . VA , due to the increase in V 0 dead ; (c) a decrease in the PACO2 due to the decrease in VLCO2. Table 5. Comparison of the parameter nominal values and optimization results for each fitting stage. Parameter Nominal Value Fitted Value Units Base fitting approach KpO24.7200 ×10−94.3473 ×10−9mm Hg−4.9 KRlv 3.7500 ×10−43.9120 ×10−4s/mm Hg MRBCO29.0000 ×10−49.3366 ×10−4L/s STPD Kbg 17.4000 16.6734 Dimensionless KcCO20.2332 0.2395 mm Hg−1 T00.5800 0.5877 s Pmax 50.0000 42.4337 cm H2O Ers 21.9000 22.2940 cm H2O/L Specific fitting approach GT,v0.0900 0.0951 Dimensionless Rsa 0.0600 0.0522 mm Hg·s/mL KE,lv 0.0140 0.0105 mL−1 Φmax 20.0000 22.0349 Dimensionless λ10.8600 0.8901 Dimensionless V0dead 0.1587 0.2059 L n1.1010 1.0157 Dimensionless C19.0000 9.5133 mmol/L VLCO23.0000 2.1479 L Stimulus-related fitting approach γsh,max 9.0 6.3 spikes/s γsp,max 5.50 3.85 spikes/s γsv,max 64.90 84.37 spikes/s γv,max 1.9000 2.0708 spikes/s Wt,sh 0.4000 0.5043 Dimensionless Wt,sp 0.4000 0.4016 Dimensionless Wt,sv 0.4000 0.4268 Dimensionless Wt,v0.4000 0.4341 Dimensionless gM40 28 Dimensionless kr,am 24.1700 27.7488 s/mL Aim 50.0000 59.6992 mm Hg gabd 3.3900 4.0826 mm Hg/L gthor 6.800 6.818 mm Hg/L KcMRV 1.000 0.895 Dimensionless λ20.4890 0.3423 Dimensionless Where STPD is standard temperature and pressure, dry. For the Stimulus-related fitting approach, the following results were obtained: (a) a decrease of the sympathetic activity related to the regulation of peripheral resistances and heart elastances and an increase for the sympathetic activity for venous volumes due to the modifications in γsh,max , γsp,max and γsv,max ; (b) a decrease in the tissues’ hypoxic effect on the vasodilation action of Ramp due to the decrease in gM ; (c) an increase in the muscular and respiratory pump activity on venous return due to increased Aim and gabd ; (d) a decreased weighting of . WE regarding . WT due to decreased λ2 . The results obtained Diagnostics 2023,13, 908 22 of 33 from NT are mainly related to increased sympathetic activity regarding heart elastances, but the I-EP results overshadowed this effect. 3.3. Validation 3.3.1. Steady-State Response Figure 4presents the steady-state model predictions under nominal conditions and at each fitting stage. Eight equidistant step inputs of . VCO2 , from 0.3 L/min to 1.0 L/min were used. Diagnostics 2023, 13, x FOR PEER REVIEW 23 of 35 Figure 4. Steady-state model predictions for each cardiorespiratory variable evaluated. Results are shown as a function of 𝑉󰇗𝐶𝑂 values. Gray dots denote the experimental data limited at each subject’s AT; the black dot-line the experimental data average; the square marker the model predictions under nominal conditions (Nominal); the cross marker the predictions of the model at the first fitting stage (Base); the plus sign marker the model predictions at the second fitting stage (Specific); and the diamond marker the model predictions at the third fitting stage (Stimulus). The subfigures from (a) to (j) correspond to the result for each of the variables evaluated. The steady-state results confirm the overall improvement of the model predictions after the fitting stages. The base fitting stage improves respiratory variables and gas exchange predictions during rest and exercise. Similar behaviors as nominal results are observed, but with an offset change, wherein the model results are closer to the experimental data’s mean trend. The predictions of cardiovascular variables showed only an improvement in PS. The model’s specific fitting improves most respiratory and cardiovascular predictions, evidencing behavior modifications under the stimulus’s increase concerning the base fitting results. V󰇗E, VT, BF, PS, PM, and PD showed improvements mainly regarding rest and moderate stimulus levels, while HR improved his behavior for high levels of exercise. The stimulus fitting does not significantly improve the model accuracy but primarily benefits the systematic blood pressure predictions. Figure 5 shows the PE values obtained from the model predictions under nominal conditions and at each fitting stage. The mean, median, interquartile range, overall values, and statistically significant differences among the predictions are presented. Figure 4. Steady-state model predictions for each cardiorespiratory variable evaluated. Results are shown as a function of . VCO2 values. Gray dots denote the experimental data limited at each subject’s AT; the black dot-line the experimental data average; the square marker the model predictions under nominal conditions (Nominal); the cross marker the predictions of the model at the first fitting stage (Base); the plus sign marker the model predictions at the second fitting stage (Specific); and the diamond marker the model predictions at the third fitting stage (Stimulus). The subfigures from ( a ) to (j) correspond to the result for each of the variables evaluated. The steady-state results confirm the overall improvement of the model predictions after the fitting stages. The base fitting stage improves respiratory variables and gas exchange predictions during rest and exercise. Similar behaviors as nominal results are observed, but with an offset change, wherein the model results are closer to the experimental data’s mean trend. The predictions of cardiovascular variables showed only an improvement in PS. The model’s specific fitting improves most respiratory and cardiovascular predictions, evidencing behavior modifications under the stimulus’s increase concerning the base fitting results. . VE , VT, BF, PS, PM, and PD showed improvements mainly regarding rest and moderate stimulus levels, while HR improved his behavior for high levels of exercise. The stimulus fitting does not significantly improve the model accuracy but primarily benefits the systematic blood pressure predictions. Figure 5shows the PE values obtained from the model predictions under nominal conditions and at each fitting stage. The mean, median, interquartile range, overall values, and statistically significant differences among the predictions are presented. The PE results confirm the steady-state predictions’ observations. The highest errors are obtained under nominal conditions, and the last fitting stage reduces the overall PE by 21.9%. The most significant PE value changes are presented for the predictions of respiratory variables and systemic blood pressure measurements, with statistically significant differences for . VE , VT, PS, PM, and PD. Although the results for PACO2 and Diagnostics 2023,13, 908 23 of 33 PAO2 do not show significant differences between consecutive fitting stages, a lower dispersion for PACO2 and an improvement of 0.07% and 39.30% was obtained for PACO2 and PAO2 predictions, respectively. The most significant decrease in the overall PE was obtained for the specific fitting approach, which also had the most negligible negative effect on the previous prediction results, followed, in order, by the base fitting approach and the stimulus-related fitting approach. Diagnostics 2023, 13, x FOR PEER REVIEW 24 of 35 Figure 5. PE results for steady-state predictions under nominal conditions and at each fitting stage. The bar graph represents the errors’ median values, and the whiskers represent the interquartile distance. The gray bars represent the PE under nominal conditions (Nominal), the line pattern bars the results at the first fitting stage (Base), the dots pattern bars the results at the second fitting stage (Specific), and the white bars represent the results at the third fitting stage (Stimulus). Symbols (*) and (**) highlight the statistically significant differences found between the obtained PE values (ρ < 0.05 and ρ < 0.01, respectively). The PE results confirm the steady-state predictions’ observations. The highest errors are obtained under nominal conditions, and the last fitting stage reduces the overall PE by 21.9%. The most significant PE value changes are presented for the predictions of respiratory variables and systemic blood pressure measurements, with statistically significant differences for V󰇗E, VT, PS, PM, and PD. Although the results for PACO and PAO do not show significant differences between consecutive fitting stages, a lower dispersion for PACO and an improvement of 0.07% and 39.30% was obtained for PACO and PAO predictions, respectively. The most significant decrease in the overall PE was obtained for the specific fitting approach, which also had the most negligible negative effect on the previous prediction results, followed, in order, by the base fitting approach and the stimulus-related fitting approach. Table 6 presents the PE mean and standard deviation values for the model predictions under nominal conditions and at each fitting stage for the related subsystems. Table 6. Prediction error results (%) at each fitting stage for the model subsystems. Fitting Stage Cardiovascular Respiratory Mechanics Gas Exchange Nominal 8.33 ± 1.73 15.54 ± 3.87 4.03 ± 0.97 Base 8.71 ± 1.62 13.25 ± 2.02 3.27 ± 1.67 Specific 8.22 ± 1.46 10.83 ± 2.13 3.24 ± 1.78 Stimulus-related 7.41 ± 1.53 11.16 ± 1.35 3.28 ± 1.78 Table 6 shows the model PE under nominal conditions and at each fitting stage. The results are presented in function to the model subsystems: cardiovascular, respiratory mechanics, and gas exchange. In general, a significant decrease in PE was found at each fitting stage. Slight adverse effects were evidenced at the base and stimulus-related fitting approaches. Improvements associated with cardiovascular variables were mainly obtained in the last fitting stage. They were related to the effects of exercise mechanisms on systemic arterial pressures. Figure 5. PE results for steady-state predictions under nominal conditions and at each fitting stage. The bar graph represents the errors’ median values, and the whiskers represent the interquartile distance. The gray bars represent the PE under nominal conditions (Nominal), the line pattern bars the results at the first fitting stage (Base), the dots pattern bars the results at the second fitting stage (Specific), and the white bars represent the results at the third fitting stage (Stimulus). Symbols (*) and (**) highlight the statistically significant differences found between the obtained PE values ( ρ< 0.05 and ρ< 0.01, respectively). Table 6presents the PE mean and standard deviation values for the model predictions under nominal conditions and at each fitting stage for the related subsystems. Table 6. Prediction error results (%) at each fitting stage for the model subsystems. Fitting Stage Cardiovascular Respiratory Mechanics Gas Exchange Nominal 8.33 ±1.73 15.54 ±3.87 4.03 ±0.97 Base 8.71 ±1.62 13.25 ±2.02 3.27 ±1.67 Specific 8.22 ±1.46 10.83 ±2.13 3.24 ±1.78 Stimulus-related 7.41 ±1.53 11.16 ±1.35 3.28 ±1.78 Table 6shows the model PE under nominal conditions and at each fitting stage. The results are presented in function to the model subsystems: cardiovascular, respiratory mechanics, and gas exchange. In general, a significant decrease in PE was found at each fitting stage. Slight adverse effects were evidenced at the base and stimulus-related fitting approaches. Improvements associated with cardiovascular variables were mainly obtained in the last fitting stage. They were related to the effects of exercise mechanisms on systemic arterial pressures. The most significant results for gas exchange predictions were obtained after the model’s base fitting. The other stages do not present significant variations in PE, as Figures 4and 5show. The results related to respiratory mechanics present the most significant variations between stages. This subsystem has the highest contribution to PE under nominal conditions and was reduced by around 30% after the specific fitting approach. Diagnostics 2023,13, 908 24 of 33 3.3.2. Transient Response Figure 6presents the model transient results under nominal conditions and at each fitting stage for a single-step load simulation. Experimental and model-predicted data are depicted as proportional changes to the variable’s initial values. Comparisons for every variable are presented, even though the experimental record length was shorter than the model simulation time. Diagnostics 2023, 13, x FOR PEER REVIEW 25 of 35 The most significant results for gas exchange predictions were obtained after the model’s base fitting. The other stages do not present significant variations in PE, as Figures 4 and 5 show. The results related to respiratory mechanics present the most significant variations between stages. This subsystem has the highest contribution to PE under nominal conditions and was reduced by around 30% after the specific fitting approach. 3.3.2. Transient Response Figure 6 presents the model transient results under nominal conditions and at each fitting stage for a single-step load simulation. Experimental and model-predicted data are depicted as proportional changes to the variable’s initial values. Comparisons for every variable are presented, even though the experimental record length was shorter than the model simulation time. Figure 6. Transient results for a single load step. Model inputs correspond to variations of 𝑉󰇗𝑂 and 𝑉󰇗𝐶𝑂 from 0.64 to 0.82 L/min and 0.58 to 0.76 L/min, respectively. They were obtained from the experimental mean values at the beginning and end of the first exercise load step. The dotted lines are the simulation results under nominal conditions (Nominal); the dashed lines are the simulation results at the first fitting stage (Base); the dash-dot lines are the simulation results at the second fitting stage (Specific); and the dotted lines with cross marker are the simulation results at the third fitting stage (Stimulus). The simulation results are compared with the corresponding experimental data. The solid gray and black lines represent the subject-by-subject experimental data restricted at AT, and their total mean value. The subfigures from (a) to (i) correspond to the result for each of the variables evaluated. The simulation results do not present significant variations for the different fitting stages, consistent with the proposed fitting strategies focused on predictions in the steady state. As part of the validation strategy, it is shown that the predictions had consistent behaviors regarding the dynamic response of the experimental data, mainly highlighting the similarity for HR and V󰇗E. The rest of the predictions do not entirely mimic the dynamics of the experimental data. Table 7 shows the model predictions’ settling times for a single load step simulation under nominal conditions and at each fitting stage. The settlement time results showed changes among the fitting stages. The obtained time values show adverse effects on the respiratory variables’ predictions and significant improvements concerning the variables of the gas exchange system. No significant effect was obtained for HR. Figure 6. Transient results for a single load step. Model inputs correspond to variations of . VO2and . VCO2 from 0.64 to 0.82 L/min and 0.58 to 0.76 L/min, respectively. They were obtained from the experimental mean values at the beginning and end of the first exercise load step. The dotted lines are the simulation results under nominal conditions (Nominal); the dashed lines are the simulation results at the first fitting stage (Base); the dash-dot lines are the simulation results at the second fitting stage (Specific); and the dotted lines with cross marker are the simulation results at the third fitting stage (Stimulus). The simulation results are compared with the corresponding experimental data. The solid gray and black lines represent the subject-by-subject experimental data restricted at AT, and their total mean value. The subfigures from ( a ) to ( i ) correspond to the result for each of the variables evaluated. The simulation results do not present significant variations for the different fitting stages, consistent with the proposed fitting strategies focused on predictions in the steady state. As part of the validation strategy, it is shown that the predictions had consistent behaviors regarding the dynamic response of the experimental data, mainly highlighting the similarity for HR and . VE . The rest of the predictions do not entirely mimic the dynamics of the experimental data. Table 7shows the model predictions’ settling times for a single load step simulation under nominal conditions and at each fitting stage. The settlement time results showed changes among the fitting stages. The obtained time values show adverse effects on the respiratory variables’ predictions and significant improvements concerning the variables of the gas exchange system. No significant effect was obtained for HR. Table 7. Settling time (seconds) of the model predictions for each fitting stage. Fitting Stage . VE VT BF TI PACO2PAO2HR Nominal 280.8 277.9 290.4 277.9 2544.6 2831.0 290.6 Base 314.5 314.4 328.2 314.5 2136.1 2795.4 281.9 Specific 295.1 294.0 276.8 266.0 1506.6 2407.5 281.6 Stimulus-related 351.6 351.5 351.6 351.6 1800.6 1622.8 284.3 Diagnostics 2023,13, 908 25 of 33 4. Discussion 4.1. Parameter Classification This work presents a methodology focused on selecting reduced sets of model parameters that can be reliably fitted in steady-state conditions following several classification approaches. The methodology is based on applying complementary and sequential techniques of parameter classification that evaluate criteria, such as the role in the model, identifiability, and sensitivity to its variation. It involves modifications to reported techniques that evaluate different stimulus levels, imply variations in the parameter values, and consider the experimental data by evaluating the error contribution. The parameters’ role classification initially constrained those parameters that should be considered for subsequent selection and fitting procedures. Five role classification groups of the parameters most commonly found in mathematical models of physiological systems were proposed. They were based on time constants, conversion parameters, covariates, initial values, and gain and thresholds that facilitated the case study model’s parameter selection. The parameter identifiability classification was based on Jacobian matrix analysis, which focused on determining which parameters can be reliably fit from available experimental data. Different stimulus levels were evaluated, considering their effects on the model’s prediction behavior. The first three parameters, classified as most identifiable for the case study, were selected (Figure 2a). Each one belongs to one of the main subsystems of the evaluated model (cardiovascular, respiratory, and gas exchange), and their well-conditionality for reliable estimation could be related to the order of magnitude of their nominal values (Table 5). A single significant difference between the parameters in the first position concerning the others presented shows the difficulty of finding a unique identification solution considering more than one parameter (Figure 2a) [ 27 ]. This fact is related to the model’s limitation to better mimic the experimental data regarding the results under nominal conditions, which could already be considered sufficiently close. The above can be related to the small decreases in overall PE at each fitting stage (Figure 5) and was involved in the value selection of termination tolerance in the function evaluation for the optimization algorithm (Table 1). Parameter sensitivity analysis is based on finite difference measurements. It evaluates model variable predictions in steady-state conditions under parameter variations and different stimulus levels. In particular, steady-state conditions allowed considering variable magnitude changes instead of changes to temporary lags among simulated variables, mainly observable in variables with oscillating dynamic behaviors such as Psa , BF, PAO2 and PACO2 . Normalization measures were applied to the obtained relative sensitivities for each stimulus level and variable to obtain unbiased results (Equations (15) and (19)). Finally, given that the sensitivity analysis does not consider the closeness between the predicted and experimental data, PE for each case (variable and stimulus level) was used as a weighting factor to highlight cases where the error was exceptionally high. This consideration results in classifications that are more consistent with the selecting and fitting procedures (Equations (15)–(22)). These modifications are related to obtaining the most respiratory control and mechanics parameters in the first positions of the total sensitivity ranking (Figure 2b), considering that their predictions provide the most significant error under nominal conditions at different stimulus levels (Figure 5and Table 6). For the specific rankings by variable, obtaining parameters related to the own systems and controllers in the first positions is related to the equitable evaluation of the sensitivities, except for PAO2 and PACO2 due to their high sensitivity concerning variables from other systems (Figure 3h,i). 4.2. Parameter Selection Parameter selection based on the parameter’s classification according to their role in the model allowed a significant reduction of the complexity associated with subsequent selection and fitting procedures, mainly evidenced by the computational cost involved in Diagnostics 2023,13, 908 32 of 33 13. Poon, C.S. Ventilatory Control in Hypercapnia and Exercise: Optimization Hypothesis. J. Appl. Physiol. 1987 ,62, 2447–2459. [CrossRef] [PubMed] 14. Hernández, A.M. Sistema de Control Respiratorio Ante Estímulos y Patologías: Análisis, Modelado y Simulación; OmniScript; Publicia: Paris, France, 2014; ISBN 978-3-639-55914-9. 15. Magosso, E.; Ursino, M. Cardiovascular Response to Dynamic Aerobic Exercise: A Mathematical Model. Med. Biol. Eng. Comput. 2002,40, 660–674. [CrossRef] 16. Silva, S.C.D.; Monteiro, W.D.; Farinatti, P.D.T.V. Exercise Maximum Capacity Assessment: A Review on the Traditional Protocols and the Evolution to Individualized Models. Rev. Bras. Med. Esporte 2011,17, 363–369. [CrossRef] 17. Lu, K.; Clark, J.W.; Ghorbel, F.H.; Ware, D.L.; Bidani, A. A Human Cardiopulmonary System Model Applied to the Analysis of the Valsalva Maneuver. Am. J. Physiol. Circ. Physiol. 2001,281, H2661–H2679. [CrossRef] 18. Cheng, L.; Khoo, M.C.K. Modeling the Autonomic and Metabolic Effects of Obstructive Sleep Apnea: A Simulation Study. Front. Physiol. 2011,2, 111. [CrossRef] 19. Albanese, A.; Cheng, L.; Ursino, M.; Chbat, N.W. An Integrated Mathematical Model of the Human Cardiopulmonary System: Model Development. Am. J. Physiol. Circ. Physiol. 2016,310, H899–H921. [CrossRef] 20. Poon, C.S.; Lin, S.L.; Knudson, O.B. Optimization Character of Inspiratory Neural Drive. J. Appl. Physiol. 1992 ,72, 2005–2017. [CrossRef] 21. Fincham, W.F.; Tehrani, F.T. A Mathematical Model of the Human Respiratory System. J. Biomed. Eng. 1983 ,5, 125–133. [CrossRef] [PubMed] 22. Serna, L.Y.; Mañanas, M.A.; Hernández, A.M.; Rabinovich, R.A. An Improved Dynamic Model for the Respiratory Response to Exercise. Front. Physiol. 2018,9, 69. [CrossRef] [PubMed] 23. Ellwein, L.M. Parameter Identifiability of a Respiratory Mechanics Model in an Idealized Preterm Infant. arXiv Prepr. 2018 , arXiv:1808.00998. 24. Ipsen, I.C.F.; Kelley, C.T.; Pope, S.R. Rank-Deficient Nonlinear Least Squares Problems and Subset Selection. SIAM J. Numer. Anal. 2011,49, 1244–1266. [CrossRef] 25. Summers, R.L.; Ward, K.R.; Witten, T.; Convertino, V.A.; Ryan, K.L.; Coleman, T.G.; Hester, R.L. Validation of a Computational Platform for the Analysis of the Physiologic Mechanisms of a Human Experimental Model of Hemorrhage. Resuscitation 2009 ,80, 1405–1410. [CrossRef] 26. Sarmiento, C.A.; Hernandez, A.M.; Serna Higuita, L.Y.; Mañanas, M.Á. An Integrated Mathematical Model of the Cardiovascular and Respiratory Response to Exercise: Model-Building and Comparison with Reported Models. Am. J. Physiol. Circ. Physiol. 2021 , 320, H1235–H1260. [CrossRef] 27. Burth, M.; Verghese, G.C.; Velez-Reyes, M. Subset Selection for Improved Parameter Estimation in On-Line Identification of a Synchronous Generator. IEEE Trans. Power Syst. 1999,14, 218–225. [CrossRef] 28. Pope, S.R. Parameter Identification in Lumped Compartment Cardiorespiratory Models; North Carolina State University: Raleigh, NC, USA, 2009. 29. Aoi, M.C.; Kelley, C.T.; Novak, V.; Olufsen, M.S. Optimization of a Mathematical Model of Cerebral Autoregulation Using Patient Data. IFAC Proc. Vol. 2009,42, 181–186. [CrossRef] 30. Eslami, M. Theory of Sensitivity in Dynamic Systems; Springer: Berlin/Heidelberg, Germany, 1994; ISBN 978-3-662-01634-3. 31. Ellwein, L.M.; Tran, H.T.; Zapata, C.; Novak, V.; Olufsen, M.S. Sensitivity Analysis and Model Assessment: Mathematical Models for Arterial Blood Flow and Blood Pressure. Cardiovasc. Eng. 2008,8, 94–108. [CrossRef] 32. Pope, S.R.; Ellwein, L.M.; Zapata, C.L.; Novak, V.; Kelley, C.T.; Olufsen, M.S. Estimation and Identification of Parameters in a Lumped Cerebrovascular Model. Math. Biosci. Eng. 2009,6, 93–115. [CrossRef] [PubMed] 33. Serna, L.Y.; Mañanas, M.Á.; Marín, J.; Hernández, A.M.; Benito, S. Optimization Techniques in Respiratory Control System Models. Appl. Soft Comput. 2016,48, 431–443. [CrossRef] 34. Fix, L.E.; Khoury, J.; Moores, R.R.; Linkous, L.; Brandes, M.; Rozycki, H.J. Theoretical Open-Loop Model of Respiratory Mechanics in the Extremely Preterm Infant. PLoS ONE 2018,13, e0198425. [CrossRef] 35. Fonoberova, M.; Mezi´c, I.; Buckman, J.F.; Fonoberov, V.A.; Mezi´c, A.; Vaschillo, E.G.; Mun, E.Y.; Vaschillo, B.; Bates, M.E. A Computational Physiology Approach to Personalized Treatment Models: The Beneficial Effects of Slow Breathing on the Human Cardiovascular System. Am. J. Physiol. Hear. Circ. Physiol. 2014,307, H1073–H1091. [CrossRef] 36. Koziel, S.; Yang, X.-S. Computational Optimization Methods and Algorithms, 1st ed.; Koziel, S., Yang, X.-S., Eds.; Springer: Berlin/Heidelberg, Germany, 2011; ISBN 9783642175565. 37. Hansen, N.; Ostermeier, A. Completely Derandomized Self-Adaptation in Evolution Strategies. Evol. Comput. 2001 ,9, 159–195. [CrossRef] 38. Magosso, E.; Cavalcanti, S.; Ursino, M. Theoretical Analysis of Rest and Exercise Hemodynamics in Patients with Total Cavopulmonary Connection. Am. J. Physiol. Heart Circ. Physiol. 2002,282, H1018–H1034. [CrossRef] 39. Chai, T.; Draxler, R.R. Root Mean Square Error (RMSE) or Mean Absolute Error (MAE)?—Arguments against Avoiding RMSE in the Literature. Geosci. Model Dev. 2014,7, 1247–1250. [CrossRef] 40. Ait-Amir, B.; Pougnet, P.; El Hami, A. Meta-Model Development. In Embedded Mechatronic Systems 2; El Hami, A., Pougnet, P., Eds.; Elsevier: London, UK; Oxford, UK, 2015; Volume 2, pp. 151–179, ISBN 9780081004692. Diagnostics 2023,13, 908 33 of 33 41. American Thoracic Society; American College of Chest Physicians. ATS/ACCP Statement on Cardiopulmonary Exercise Testing. Am. J. Respir. Crit. Care Med. 2003,167, 211–277. [CrossRef] 42. Heyward, V.; Gibson, A. Advance Fitness Assessment and Exercise Prescription, 7th ed.; Human Kinetics: Champaign, IL, USA, 2014; ISBN 9781450466004. 43. Latin, R.W.; Berg, K.E.; Smith, P.; Tolle, R.; Woodby-Brown, S. Validation of a Cycle Ergometry Equation for Predicting Steady-Rate VO2. Med. Sci. Sports Exerc. 1993,25, 970–974. [CrossRef] 44. Harada, T.; Kubo, H.; Mori, T.; Sato, T. Pulmonary and Cardiovascular Integrated Model Controlled with Oxygen Consumption. In Proceedings of the 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference, Shanghai, China, 17–18 January 2006 ; Volume 7, pp. 304–307. 45. Mananas, M.A.; Hernandez, A.M.; Romero, S.; Grino, R.; Rabinovich, R.; Benito, S.; Caminal, P. Analysis of Respiratory Models at Different Levels of Exercise, Hypercapnia and Hypoxia. In Proceedings of the 25th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (IEEE Cat. No.03CH37439), Cancun, Mexico, 17–21 September 2003; Volume 3, pp. 2754–2757. 46. Cooper, C.B.; Storer, T.W. Exercise Testing and Interpretation, 1st ed.; Cambridge University Press: Cambridge, UK, 2001; ISBN 9780521648424. 47. Ursino, M. Interaction between Carotid Baroregulation and the Pulsating Heart: A Mathematical Model. Am. J. Physiol. 1998 ,275, H1733–H1747. [CrossRef] 48. Ursino, M.; Magosso, E. Acute Cardiovascular Response to Isocapnic Hypoxia. I. A Mathematical Model. Am. J. Physiol. Heart Circ. Physiol. 2000,279, H149–H165. [CrossRef] 49. Batzel, J.J.; Kappel, F.; Timischl-Teschl, S. A Cardiovascular-Respiratory Control System Model Including State Delay with Application to Congestive Heart Failure in Humans. J. Math. Biol. 2005,50, 293–335. [CrossRef] [PubMed] 50. Cheng, L.; Ivanova, O.; Fan, H.-H.; Khoo, M.C.K. An Integrative Model of Respiratory and Cardiovascular Control in SleepDisordered Breathing. Respir. Physiol. Neurobiol. 2010,174, 4–28. [CrossRef] [PubMed] 51. Grodins, F.S.; Buell, J.; Bart, A.J. Mathematical Analysis and Digital Simulation of the Respiratory Control System. J. Appl. Physiol. 1967,22, 260–276. [CrossRef] [PubMed] 52. Levy, M.N.; Zieske, H. Autonomic Control of Cardiac Pacemaker Activity and Atrioventricular Transmission. J. Appl. Physiol. 1969,27, 465–470. [CrossRef] 53. Serna Higuita, L.Y.; Mananas, M.A.; Mauricio Hernandez, A.; Marina Sanchez, J.; Benito, S. Novel Modeling of Work of Breathing for Its Optimization During Increased Respiratory Efforts. IEEE Syst. J. 2016,10, 1003–1013. [CrossRef] 54. Otis, A.B.; Fenn, W.O.; Rahn, H. Mechanics of Breathing in Man. J. Appl. Physiol. 1950,2, 592–607. [CrossRef] 55. Tehrani, F.T. Mathematical Analysis and Computer Simulation of the Respiratory System in the Newborn Infant. IEEE Trans. Biomed. Eng. 1993,40, 475–481. [CrossRef] 56. Ursino, M.; Antonucci, M.; Belardinelli, E. Role of Active Changes in Venous Capacity by the Carotid Baroreflex: Analysis with a Mathematical Model. Am. J. Physiol. 1994,267, H2531–H2546. [CrossRef] 57. Ursino, M.; Magosso, E. A Theoretical Analysis of the Carotid Body Chemoreceptor Response to O 2 and CO 2 Pressure Changes. Respir. Physiol. Neurobiol. 2002,130, 99–110. [CrossRef] [PubMed] 58. Sinning, D.; Kasner, M.; Westermann, D.; Schulze, K.; Schultheiss, H.-P.; Tschöpe, C. Increased Left Ventricular Stiffness Impairs Exercise Capacity in Patients with Heart Failure Symptoms Despite Normal Left Ventricular Ejection Fraction. Cardiol. Res. Pract. 2011,2011, 1–10. [CrossRef] 59. Westermann, D.; Kasner, M.; Steendijk, P.; Spillmann, F.; Riad, A.; Weitmann, K.; Hoffmann, W.; Poller, W.; Pauschinger, M.; Schultheiss, H.P.; et al. Role of Left Ventricular Stiffness in Heart Failure with Normal Ejection Fraction. Circulation 2008 ,117, 2051–2060. [CrossRef] [PubMed] 60. Pawelczyk, J.A.; Hanel, B.; Pawelczyk, R.A.; Warberg, J.; Secher, N.H. Leg Vasoconstriction during Dynamic Exercise with Reduced Cardiac Output. J. Appl. Physiol. 1992,73, 1838–1846. [CrossRef] [PubMed] 61. Spencer, J.L.; Firouztale, E.; Mellins, R.B. Computational Expressions for Blood Oxygen and Carbon Dioxide Concentrations. Ann. Biomed. Eng. 1979,7, 59–66. [CrossRef] 62. Angleys, H.; Jespersen, S.N.; Østergaard, L. The Effects of Capillary Transit Time Heterogeneity on the BOLD Signal. Hum. Brain Mapp. 2018,39, 2329–2352. [CrossRef] 63. Edwards, A.D.; Jennings, S.J.; Newstead, C.G.; Wolff, C.B. The Effect of Increased Lung Volume on the Expiratory Rate of Rise of Alveolar Carbon Dioxide Tension in Normal Man. J. Physiol. 1983,344, 81–88. [CrossRef] 64. Wasserman, K.; Hansen, J.; Sietsema, K.; Sue, D.Y.; Stringer, W.W.; Sun, X.-G.; Whipp, B.J. Principles of Exercise Testing and Interpretation: Including Pathophysiology and Clinical Applications, 5th ed.; Lippincott Williams & Wilkins, Ed.; Lippincott Williams & Wilkins: Philadelphia, PA, USA, 2012; ISBN 978-1-60913-899-8. 65. Hester, R.L.; Brown, A.J.; Husband, L.; Iliescu, R.; Pruett, D.; Summers, R.; Coleman, T.G. HumMod: A Modeling Environment for the Simulation of Integrative Human Physiology. Front. Physiol. 2011,2, 12. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.