Dynamic Fitting Strategy for Physiological Models: A Case Study of Cardiorespiratory Model for Simulation of Incremental Aerobic Exercise

1. Introduction

The increase and variety of chronic diseases, the consequences of the reactive medical approach, and the unsatisfactory results of generic diagnoses and treatments suggest the need for a predictive, preventive, and personalized medical paradigm [1,2]. It is based on integrating healthcare personalization concepts, anticipation and risks assessment, and determination of optimal treatment based on predictive tools [2,3]. Among its applications, those that predict the personalized health status and the risk of presenting pathologies based on physiological information stand out [4,5,6].

The mathematical modeling and computational simulation of physiological systems is a support tool for medical science, whose application has allowed the evaluation of clinical treatments, the development of training tools and medical education, and the design of medical equipment [7,8,9,10]. These applications are possible due to models’ ability to predict clinically critical physiological variables considering different conditions, subjects, and stimuli. The detailed description of the elements, functions, and interactions stand out among its characteristics, differentiating it from other prediction tools and allowing the correct understanding of complex regulations and physiological processes [11,12].

The models of the cardiorespiratory system stand out as one of the most relevant for clinical applications because their predicted variables are the result of complex regulatory processes aimed at maintaining homeostasis, whose analysis is essential for the diagnosis, monitoring, and treatment of different diseases [13,14,15]. Nowadays, there are several studies focused on verifying and using cardiorespiratory models in different situations, including exercise [16,17], cardiovascular surgeries [18], respiratory disorders [19,20], mechanical ventilatory support [21,22], and others. Among these situations, physical exercise is especially important for diagnosis and clinical prediction because it is a natural stimulus that generates significant variations of the related systems and whose simulation results are characterized for healthy human subjects [23,24].

Despite the advantages and variety of validated mathematical models of physiological systems, their usefulness in medical practice, especially under the predictive, preventive, and personalized approach, is restricted by different factors. One of the most important is the mathematical complexity, which is greater as precision and biological relevance increase, and is related to a large number of parameters that are difficult to measure experimentally [11,12]. Although there are reported nominal values of parameters for the model simulation, these were obtained under specific experimental conditions and subject characteristics. Their use for predictions gives a constrained and general idea of the physiological system dynamics, neglecting the properties that vary between people and conditions [11].

There are reported fitting techniques that, based on experimental data such as vital signs, allow the personalization of models. They are related to the identification of parameters that are difficult to measure experimentally [12,25,26]. Even so, these techniques only focus on identifying constant characteristics, so the usefulness of the fitted models is constrained to the static prediction of the physiological response.

In this context, a correct application of physiological models as a predictive, preventive, and personalized medicine tool would involve a strategy to continuously optimize and customize their parameters as a function of time. There are reported works in which a continuous optimization of model parameters is considered to account for the effect of time on the change in the systems’ behavior, obtaining good prediction results for applications related to electronic, control, and image processing systems [27,28,29,30]. This approach, known as time-varying systems modeling, has also been applied to models of physiological systems, but, its application has focused mainly on small models (i.e., a single system model described by a small number of equations), for which the experimental data are directly related to the prediction of the interest variables [31]. Although its application could be extended to more complex models, in which a complete self-regulation of the system is considered, and the responses of different physiological variables are evaluated, this approach assumes that predictions at each time instant depend on the directly related parameters. Therefore, its use neglects the possibility that the changes in the system responses are due to the closed-loop effect of other mechanisms, systems, or controllers, in which the parameter optimizations are related to a better fit of the model.

This work aims to propose a novel strategy for the dynamic fitting of physiological models with complete self-regulation. This strategy will help to predict the future state of the modeled physiological system under the consideration of a large number of parameters and a limited amount of experimental data for fitting. This strategy is composed of the following procedures: a) the dynamic fitting of the physiological model, b) the modeling of the parameters’ time trend, and c) the strategy validation considering the future steady-state response of the modeled physiological system. As a case study, the application of the strategy in a cardiorespiratory model under dynamic aerobic exercise is described in detail. Experimental records obtained from a longitudinal study were used as evaluation scenarios from a population and personalized perspective.

This work is structured as follows. First, the dynamic fitting strategy is presented, including identifying the parameter values at different time instants, modeling their temporal trend, and validating its use to predict future states of the modeled physiological system. Subsequently, the case study is presented, including the qualitative description of the cardiorespiratory model, the experimental data, the computational implementation, and the data analysis. Then, the results of each procedure are presented for the two fitting scenarios, focusing on the strategy validation. Finally, the results of the proposed strategy are discussed.

2. Materials and Methods

2.1. Dynamic Fitting Strategy

The proposed strategy is based on the consecutive application of a static fitting technique of physiological models based on experimental observations of a longitudinal study. The fitting technique should focus on finding the model parameter values that minimize the differences between the predictions and experimental data in each instant. As physiological systems modify their characteristics over time due to natural biological degradation and factors related to different habits, events, and diseases [32,33,34,35,36], the proposed strategy uses the longitudinal results of the static fit of the models. Thus, the variations of the optimized parameters values at different instants expose the temporal changes of the physiological system characteristics, wherein the modeling of their time-trend allows predicting future parameters values.

The strategy includes procedures of dynamic fitting of the physiological model, modeling the time-trend of parameters as a function of their previous results and related factors, and validating the model simulation results using predicted parameter values. These procedures are explained in detail in the following sub-sections.

2.1.1. Dynamic Fitting of the Physiological Model

This procedure focuses on obtaining the parameter optimizations that fit the physiological model at different instants regarding experimental observations from a longitudinal study. The number of results will depend on the number of study samples ($N$) since each one ($n$) corresponds to a static fit of the model regarding the data recorded at each instant of time ($\mathrm{Exp}. \mathrm{Data} t_n$). It is composed of different sub-processes that must be applied sequentially in time, considering that the results will reflect the temporal evolution of the physiological system for the same subject or population. The first fit ($n=1$) is based on the nominal values of the model parameters and the subsequent ones on the immediately preceding optimization result (Params. $t_{n-1}$). The procedure also involves the temporary recording of the parameter optimizations’ results (Reg. Params.) and experimental information related to factors that significantly influence the modeled system (Reg. Factors), useful for the subsequent dynamic parameter modeling procedure. Figure 1 shows the general flowchart for the dynamic fitting of the physiological model.

Physiological systems models are highly complex, as demonstrated by their large number of parameters. Attempting to fit them using a limited amount of experimental data can be computationally expensive and results in multiple optimization solutions that lack physiological meaning [12]. For this reason, the proposed procedure is based on the sequential application of a static fitting and validation strategy of physiological models aimed at solving those problems [37]. It is based on techniques of complementary approaches of classification, selection, and optimization of a reduced set of parameters, allowing predictions with adequate precision and accuracy for models under the consideration of a large number of parameters and different stimuli and populations. The effectiveness of this procedure has already been tested in a case study involving a comprehensive cardiorespiratory model.

The main sub-procedures, corresponding to the standardization of the simulation conditions, the classification and selection of parameters, and the parameter optimizations, are described below. Each is based on different techniques and must be applied sequentially for a single fitting procedure. The reader can refer to the paper on the static fitting and validation strategy for full details [37].

Standardization of simulation conditions

It consists of modifying the values of the model parameters by those that can be obtained from the experimental observations (the result is related to the first Update params. process in Figure 1). As a result, the simulation conditions of the model are closer to those used as a reference for the fit, and the computational cost is reduced by not considering the subsequent optimization of the modified parameters. The replaced values are generally related to the characteristics of the subjects, the stimulus evaluated, or the environmental conditions. The parameters unchanged during the standardization, which cannot be associated with the available experimental measures, remain equivalent to the nominal values. Then, they will be subsequently considered in the classification and selection sub-procedures.

Classification and Selection of Parameters

It aims to reduce the complexity of the optimization problem by identifying the reduced set of parameters most representative to be used during the model fitting [37]. It initially consists of classifying parameters based on their role, identifiability, sensitivity to variations, and relationship with the evaluated stimulus. Subsequently, the resulting parameters are selected according to different fitting approaches.

The classification by role consists of grouping the parameters into five sets according to the functions generally found in physiological models, corresponding to time constants, conversion parameters, covariates, initial values, and gain and thresholds. Considering that the implemented static fitting strategy focuses on adjusting the model predictions in steady-state, only the sets of parameters with the role of gain and threshold, whose variations have physiological sense, are selected for the application of the subsequent classification techniques.

The following classification is related to the application of the techniques of subset selection and analysis of the sensitivity. Both generate rankings for the classification of the parameters considering different and complementary criteria. The subset selection classifies the parameters according to how well-conditioned they are to be identified [38]. The sensitivity analysis technique ranks them according to their variation effect in the model predictions [12]. Both are based on measures of finite differences obtained from simulated data. According to the static fitting strategy, their application involves improvements that allow evaluating different levels of stimulus and parameter values and considering the error contribution regarding the experimental data [37].

Finally, selection criteria are applied to define which classified parameters will be considered for optimization. According to the referenced strategy, the selection of the parameters is related to three different optimization processes that must be applied sequentially [37]. The first one corresponds to the base fitting approach of the model, for which parameters with the best results of subset selection and the highest total sensitivity are selected. The second one is the specific fitting approach, for which the parameters with the best relative sensitivities for each variable, which, in turn, have a non-significant effect concerning the others, are selected. The third one is related to the stimulus, wherein the parameters that measure the impact of relevant mechanisms, based on their classification by role, are selected.

Parameters Optimization

It involves identifying the values of previous selected parameters that reduce the difference between the model predictions and the experimental data [31] (the result is related to the second Update params. process in Figure 1). This procedure consists of solving an optimization problem, in which the result is related to the selection and parametrization of an optimization algorithm, the definition of evaluation ranges of parameters, and the choice of a metric to evaluate the prediction error [26]. According to the referenced static fitting and validation strategy, three optimization stages must be applied sequentially, each related to each parameter selection criteria. The first stage corresponds to the base fitting approach and is based on the identification of the most significant parameters of the model. The second stage consists of the specific fitting approach, which aims to reduce the prediction error of each variable without significantly affecting the result of the others. The last stage is the stimulus-related fitting approach, which focuses on bringing the behavior of predictions closer to what is expected by modifying the related mechanisms.

2.1.2. Dynamic Parameter Modeling

Its objective is to model the temporal trend of the parameters so that their simulation results build the prediction of a future fit of the physiological model. There is no specific structure defined for this modeling. Still, applying autoregressive techniques is recommended considering their simplicity and the cumulative effect regarding time and other factors on physiological systems [32,33,34,35,36]. Therefore, the design of a MISO (Multiple-Input Single-Output) model is proposed for each parameter. This latter allows predictions based on their past values and the temporal record of significant factors such as habits, events, or diseases.

Figure 2 shows the general flowchart for the dynamic modeling of one parameter ($P_x$). This procedure focuses on predicting and storing its value (Reg. Pred. Params.) at the next instant of time ($n+1$) regarding the samples used for modeling ($0\to n$). N denotes the total number of time instances available. Only the parameters that modify their value as a result of the dynamic fitting procedure ($Variation\left(P_{x,0}\to P_{x,n}\right)$) are modeled and simulated ($Model/Sim.$), the others remain with previous defined value (the nominal value for the parameters that are not modified as a result of the standardization of the simulation conditions). Similar to the optimization process, it is essential to note that the predictions of the parameters could lack a physiological justification and should be evaluated before being accepted. The goodness-of-fit measures of the models and the rate of change of the values can be used as a criteria to define the use of the predictions or keep the previous values ($Criteria\left(Fit,P_{x,0}\to P_{x,n}\right)$).

2.1.3. Validation

This procedure focuses on evaluating the result of the proposed strategy from the approach of predicting the future fit of the physiological model. For this purpose, the experimental values of the variables at each instant are compared against the respective predictions of the physiological model using three sets of parameters: (a) the predicted parameters, (b) the parameters fitted for the used experimental data for each time instant, and (c) the parameters resulting from the traditional single-fit approach, which consists of a population fit in a single instant of time (the first optimization process).

Considering that the sequentially implemented static fitting strategy does not involve the response in a transitory regime, criteria of precision, accuracy, and direction of the physiological model predictions are only evaluated once the steady-state of the variables is reached. This involves the selection and use of a metric that allows the assessment of the model’s performance, usually associated with the calculation of the prediction error.

2.2. Case Study

2.2.1. Cardiorespiratory Model

The mathematical model is a compartmental self-regulated model of the cardiorespiratory system. It comprises 316 parameters and 170 equations distributed in different mechanisms that are integrated to represent the cardiovascular circulation system, the respiratory mechanics system, the gas exchange system, and the respiratory and cardiovascular controllers. It is the result of the adaptation of different previous models reported of systems, controllers, and mechanisms. It has been fitted and validated regarding experimental data from healthy adult subjects under conditions of resting and incremental aerobic exercise [37].

2.2.2. Experimental Data

The experimental data used to evaluate the strategy correspond to signals and measurements from cardiorespiratory variables recorded in an observational and longitudinal study. This was based on an incremental, submaximal, and multistate cardiopulmonary exercise test performed using a programmable cycle ergometer under controlled environmental conditions. The duration of the test was 45 minutes, divided into five stages: rest (2 minutes), warm-up (4 minutes), exercise (25 minutes), recovery (9 minutes), and final rest (5 minutes). Five consecutive steps with load increments of 25 W, each five minutes, and at a constant pedaling speed of 60 RPM were carried out during the exercise stage.

The recorded information was divided into two databases to analyze the usefulness of the dynamic fitting strategy from a population and personalized perspective. The first database (DB1) corresponds to the longitudinal record of twenty-seven male adult volunteers in three different instants with more than six months of difference between them. The second database (DB2) corresponds to the records of three male adult volunteers matched by the degree of consanguinity, medical history, and lifestyle but with different ages, which is proposed as an approximation of the longitudinal record of the same subject at six different stages of his life, and which is consistent with the approach of stratified medicine [39,40,41,42]. These subjects are male members of the same family, corresponding to the father and his two sons, who live together and have similar habits. The first two records correspond to the youngest subject, the next three to his brother, and the last one to his father.

All participants were considered healthy, with a physical fitness that included both sedentary and accustomed to physical training subjects, non-smokers, no history or symptoms of cardiovascular, pulmonary, metabolic, or neurological diseases, and without pacemakers or other implanted electrical stimulators. The records are composed of signals of oxygen uptake (${ \dot{\mathrm{V}}\mathrm{O}}_2$), carbon dioxide output (${ \dot{\mathrm{V}}\mathrm{CO}}_2$), total minute ventilation ($\dot{\mathrm{V}}\mathrm{E}$), inspiratory time (TI), tidal volume (VT), breathing frequency (BF), heart rate (HR), alveolar oxygen partial pressure (${\mathrm{PAO}}_2$) and alveolar carbon dioxide partial pressure (${\mathrm{PACO}}_2$); measurements of systolic blood pressure (PS), diastolic blood pressure (PD), and mean arterial blood pressure (PM). Records comprised of characteristics of the environment, such as the atmospheric pressure (${\mathrm{P}}_{\mathrm{atm}}$), inspired fraction of dry ${\mathrm{CO}}_2$ (${\mathrm{FiCO}}_2$) and inspired fraction of dry ${\mathrm{O}}_2$ (${\mathrm{FiO}}_2$).

Quantifiable information related to factors other than age, reported as significantly influencing the cardiorespiratory system, was also recorded. These data account for habits such as diet [33], physical activity [35], and quality of sleep [34,43], related to measures of the body mass index, average hours of weekly physical activity, and average hours of daily sleep, respectively. Other relevant factors, such as accidents or illnesses, habits such as smoking, or family history of illnesses, are not reported considering that they were considered as exclusion criteria. The use of medication was only reported by the oldest subject of the DB2 and is related to blood pressure control.

Table 1 presents the information recorded for each database related to environmental conditions, characteristics of volunteers, and factors mentioned above. The data presented for DB1 correspond to the mean and standard deviation of the measurements.

All procedures were approved by the Human Research Ethical Committee of the Department of University Research (SIU) of the University of Antioquia (approval certificates 16-59-711, 17-59-711, and 18-59-711). They are conformed to the Declaration of Helsinki. The written consent of all volunteers was obtained after giving them complete information about the experimental protocol and the related risks.

2.2.3. Computational Implementation and Data Analysis

The experimental values of ${ \dot{\mathrm{V}}\mathrm{O}}_2$ and ${ \dot{\mathrm{V}}\mathrm{CO}}_2$ were used as model inputs to simulate the levels of the stimulus of exercise. This implementation has been applied to simulate aerobic exercise in physiological models [18,44,45,46]. It is justified because the physiological response under this stimulus is directly related to metabolic gas rates of ${\mathrm{O}}_2$ consumption and ${\mathrm{CO}}_2$ production, which is directly affected by the workload during exercises (the physiological model does not consider the workload as a model input).

The average of each measurement during the last minute for each subject (except for PS, PM, and PD) was used as a representative sample of each phase of the exercise test at steady-state conditions. All experimental data were restricted using ${ \dot{\mathrm{V}}\mathrm{CO}}_2$ as a reference, starting at 0.3 l/min, considered the state of rest, coinciding with the common minimum value for records of all subjects, and ending at the respective calculated value of AT (Anaerobic threshold) because the model is defined only for aerobic exercise.

The parameters values used as reference in each parameter temporal fit of the model corresponded to the immediately preceding optimization result, except for the first fit of each database, in which the nominal values reported in [47] were used (traditional single-fit approach).

The simulation time for each stimulus level was 2000 s, which was enough to guarantee the reported settling time measure of all variables evaluated [47]. The model responses in steady-state corresponded to the mean values of each variable at the last minute of each simulated stimulus level.

The model simulations and data processing were implemented in SIMULINK/MATLAB®. The computational characteristics for the simulation were the same as reported for the model validation, corresponding to the numerical solver ODE23 (Bogacki and Shampine BS23 algorithm for solution of ordinary non-rigid differential equations) and variable step size between $1x{10}^{-2}$ s and $1x{10}^{-3}$ s. The simulations were run once due to all the model equations are determinists. The specific implementation characteristics for each procedure are described below.

Dynamic Model Fitting

Standardization of simulation conditions

The parameters of the physiological model that were assigned according to experimental data corresponded to ${\mathrm{FiO}}_2$, ${\mathrm{FiCO}}_2$, ${\mathrm{P}}_{\mathrm{atm}}$, basal respiratory tidal volume (${\mathrm{VT}}_{\mathrm{n}}$), anaerobic threshold ($\mathrm{AT}$), total blood volume (${\mathrm{V}}_{\mathrm{tot}}$) and constant values of unstressed blood volumes. ${\mathrm{FiO}}_2$, ${\mathrm{FiCO}}_2$ and ${\mathrm{P}}_{\mathrm{atm}}$ account for the environmental conditions, were considered constant in all records and equal to 640 mmHg, 0.0421 %, and 21.0379 %, respectively. Values of ${\mathrm{VT}}_{\mathrm{n}}$, $\mathrm{AT}$, ${\mathrm{V}}_{\mathrm{tot}}$ and the unstressed blood volumes are related to specific characteristics of the subjects. The value of ${\mathrm{VT}}_{\mathrm{n}}$ was determined from VT measurements; $\mathrm{AT}$ was calculated using the noninvasive technique v-slope [20,48]. ${\mathrm{V}}_{\mathrm{tot}}$ in ml was also estimated for each record as a function of the body surface area (BSA) in ${\mathrm{m}}^2$ according to the Equations (1) and (2) [49,50]. Finally, the parameters related to the unstressed blood volumes were calculated using the proportion of their reported nominal values regarding ${\mathrm{V}}_{\mathrm{tot}}$ [47].

$$\mathrm{BSA}=\sqrt{\frac{\left(\mathrm{Weigth}·\mathrm{High}\right)}{3600}},$$

(1)

$${\mathrm{V}}_{\mathrm{tot}}=1000 \left(3.29·\mathrm{BSA} -1.29\right),$$

(2)

The values of ${\mathrm{VT}}_{\mathrm{n}}$, $\mathrm{AT}$ and ${\mathrm{V}}_{\mathrm{tot}}$ used for the standardization of each fit are presented in Table 2. The data presented for DB1 correspond to the mean and standard deviation, except for $\mathrm{AT}$ whose value corresponds to the median and interquartile range considering the data distribution.

Classification and Selection of Parameters

Steady-state simulation results for variations of stimulus level and parameter values were used for classification techniques. The number, order, and values of stimulus levels and parameter variations were the same as reported for the static fitting strategy previously applied to the same model [37]. Three simulated exercise stimulus levels (using ${ \dot{\mathrm{V}}\mathrm{O}}_2$ and ${ \dot{\mathrm{V}}\mathrm{CO}}_2$ as model inputs) were evaluated as incremental steps, corresponding to states of rest, an intermediate level of exercise, and $AT$. Five uniformly distributed percentage variations in a range of ± 5% of the reference value were evaluated for the selected parameters.

Regarding the selection, the fitting approach related to the stimulus was not applied because its result does not represent a significant decrease in the prediction error, and the exercise mechanisms had already been fitted for the model [37]. The selection criteria and the number of parameters defined were also the same as those reported for the model according to the static fitting strategy.

Parameters Optimization

The optimization algorithm implemented corresponds to the evolutionary strategy with covariance matrix and adaptation (CMA-ES). It is a stochastic global optimization algorithm based on adaptive and evolutionary strategies [51]. It was selected considering the reported results of convergence speed, precision, and accuracy for the fitting of the case study model [37] and other related ones [26]. The parameter values used for implementing CMA-ES algorithm are the same as those reported for the model according to the static fitting strategy [37].

The evaluation ranges for the parameter optimizations were in accordance with the static fitting strategy [37]. The parameter values reported as nominal for the model were used as variation reference [47]. It was established as a general evaluation range of ± 30% regarding the nominal value, expanded to ± 50% only for weighting parameters that are unrelated to direct physiological measures. These ranges were modified according to the information reported about optimization results or physiological sense (values with physiological justification).

The evaluated model predictions corresponded to the simulation results in steady-state under eight consecutive, equidistant, and incremental steps of the values of ${ \dot{\mathrm{V}}\mathrm{O}}_2$ and ${ \dot{\mathrm{V}}\mathrm{CO}}_2$ between rest and AT. The experimental data for comparison were the same as the parameter’s classification procedure. The cost function implemented corresponds to Equation (3), the same RMSE-based metric used for the reported fit of the case study model [37].

$$\mathrm{CF}=\frac{1}{\mathrm{I}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{I}}\sqrt{\frac{1}{\mathrm{K}}{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{K}}{\left(\frac{{\mathrm{y}}_{{\exp }_{\mathrm{i},\mathrm{k}}}-{\mathrm{y}}_{{\mathrm{sim}}_{\mathrm{i},\mathrm{k}}}}{{\mathrm{y}}_{{\exp }_{\mathrm{i},\mathrm{k}}}}\right)}^2} ,$$

(3)

Where $\mathrm{CF}$ denotes de cost function for the model; ${\mathrm{y}}_{\exp }$ and ${\mathrm{y}}_{\mathrm{sim}}$ represent the experimental and simulated values of the variable, respectively; $\mathrm{I}$ and K denote the number of variables and stimuli levels; the subscripts i and k indicate each variable and stimulus level.

Dynamic Parameter Modeling

The parameter values corresponding to five records were predicted for the case study. For DB1 the optimized parameter values of the record three were predicted, and for DB2 the optimized parameter values of the records from three to six were predicted. Only the parameters that presented variations in their value due to the previous fit procedures were modeled (dynamic parameters). The others remain with their same value (nominal value or standardization result of the simulation conditions).

Deterministic and parametric MISO models were implemented for each dynamic parameter. The implemented model structure corresponded to an ARX model (autoregressive with exogenous input), selected to consider the effect of the past values of the outputs and inputs on the current value [31]. Equation (4) corresponds to the representation of the model structure.

$$A\left(q\right)P\left(t\right)={\displaystyle \sum }_{i=1}^I{B}_i\left(q\right)F_i\left(t-n{k}_i\right)+e\left(t\right),$$

(4)

Where $P$ represents the value of the parameter at each instant of time; $F_i$ is the $i$th input, which is associated with the factors evaluated; A and B are polynomials expressed in the time-shift operator $q^{-1}$ and are related to the parameter and the inputs, respectively; $I$ is the total number of inputs; $n{k}_i$ is the input delay; and $e\left(t\right)$ is considered random noise.

As inputs of each model, the time between records, the average weekly physical activity, the average daily sleep hours (Table 1), and AT (Table 2) were used. The time was considered as the accumulated difference between the ages in each record, assigning the value zero to the first fit of the model for each database. The AT values were used because they account for the metabolic response of the subject to cardiorespiratory changes [48,52]. Except for time, the values of the inputs of the previous record were considered as the equivalent values in the record to be predicted.

Different model orders associated with the polynomials were evaluated for each parameter, and the model with the best fit was selected. The past values of the parameter and the inputs were interpolated, considering that the design of the models requires constant sampling times. Sampling times between 0.1 and 0.4 years were evaluated in this work considering the temporal difference between records for each fit. The orders of polynomials A and B were evaluated from zero to the maximum number of samples available. $n{k}_i$ was considered zero regarding all inputs due to the limited number of available records. The coefficients of the polynomials for each model were estimated using the least-squares method. The best-order selection was based on applying Akaike’s Information Criterion (AIC) according to Equation (5).

$$\mathrm{AIC}=\log \left(\mathrm{V}\right)+\frac{2\mathrm{d}}{\mathrm{N}} ,$$

(5)

Where $\mathrm{V}$ is the loss function (normalized sum of squared prediction errors), $\mathrm{d}$ is the total number of parameters in the structure in question, and $\mathrm{N}$ is the number of data points used for the estimation.

Criteria based on the goodness of fit of the models and on the variation of the predicted value regarding a bounded range were applied to evaluate the physiological justification of the predictions. Model with goodness of fit less than 60%, with the sign change, or variations greater than 50% regarding the range of recorded values were discarded, and the values corresponding to the immediately previous record were used instead.

Validation

The validation consisted of comparing the simulation results of the case study model using the optimized parameter values for the records (from the dynamic fitting of the physiological model), the predicted values (from the dynamic parameter modeling), and the nominal values. The experimental values of variables for the respective record were used as reference. These comparisons were aimed at demonstrating the usefulness of the dynamic fitting strategy from the perspective of predicting the future state of the modeled system against the traditional approach of using population fit results in a single instant (single-fit).

The model performances were evaluated by the measure of the prediction error (PE) regarding the steady-state experimental measures of the cardiorespiratory variables. The metric used to calculate the PE corresponds to Equation (6), which is a modification of the Mean Absolute Error (MAE) that allows an interpretation of the differences as a proportion of the experimental data. It consider the error for each subject, variable and level of stimulus; and it is the same reported in the validation work of the case study model [37].

$$\mathrm{PE}(\%)=100·\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{N}}\left({\mathrm{Median}}_{\mathrm{i},\mathrm{k}}\left(\left|\frac{{\mathrm{y}}_{{\exp }_{\mathrm{i},\mathrm{j},\mathrm{k}}}-{\mathrm{y}}_{{\mathrm{sim}}_{\mathrm{i},\mathrm{j},\mathrm{k}}}}{{\mathrm{y}}_{{\exp }_{\mathrm{i},\mathrm{j},\mathrm{k}}}}\right|\right)\right),$$

(6)

Where ${\mathrm{y}}_{\exp }$ is the experimental variable value; ${\mathrm{y}}_{\mathrm{sim}}$ is the simulation prediction value; $\mathrm{N}$ denotes the number of variables; the subscripts i, j and k are indexes that represent each subject, variable, and stimulus level, respectively.

3. Results

3.1. Dynamic Fitting of the Physiological Model

Table 3 shows the PE results for the model fitted to each record of the databases. The results regarding the variables correspond to the median and the interquartile range.

For both databases, it is observed that the highest PE values are generally related to the respiratory variables and the lowest to the gas exchange variables. More specifically, DB2’s fittings generally show results with lower variabilities and lower overall PEs. It is also highlighted that the highest values of overall PE in DB2 are related to the most significant changes in the cardiorespiratory system according to the proposed approach (section 2.2.2. Experimental Data).

3.2. Dynamic Parameter Modeling

Table 4 and Table 5 compare the parameter predictions (Pred.) with the optimized results (Fitted) for the records evaluated of each database. The results for each parameter are grouped according to the directly related system or controller. The definition and description of each parameter are not presented in this section for brevity reasons. The interested reader can refer to the work of the case study model for complete details [47]. The symbol (*) indicates that the parameter was not predicted for the record because there was no modification regarding previous values, but it was optimized for the evaluated record. The symbol (**) indicates that the predicted value did not satisfy the physiological justification criteria and, therefore, the value from the previous record was used instead.

Most of the parameters optimized for DB1 show slight differences regarding the predicted ones and those that were not modified (*) (variations less than 10%). Therefore, similar simulation results could be expected. In this case, all predictions were in accordance with the defined physiological justification criteria, which suggests good fits for the restricted number of samples, and predictions with consistent variations considering the prediction time and the input data.

From the results related to DB2, it was observed that the number of non-modeled parameters (*) decreases with the number of samples, which is related to a greater coverage of prediction of the physiological characteristics related to cardiorespiratory variations. Even so, the predictions that did not satisfy the physiological justification criteria increased in the last two records (**), related to lower goodness-of-fit measures and most significant variations, mainly for record six, considering the longer prediction time. The most significant differences were obtained for the respiratory controller parameters due to the results of the static fit in each record. They could be related to higher reported errors of the physiological model regarding the respiratory variables [37].

3.3. Validation

The validation results correspond to the comparison of simulations in steady-state and the PE values in each instant evaluated. The compared simulations correspond to those obtained by using the set of parameters resulting from the predictive approach (Predicted) against those of the traditional single-fit approach (Single-Fit). The model-specific fit results for each record in the databases are presented as reference (Fitted). For the PE results, in addition to the overall value, mean, median, interquartile range, and statistically significant differences among predictions are presented.

Figure 3 shows the steady-state results of the physiological model simulations regarding the third DB1 record. It was shown that the simulations of the predictive approach were the most similar to those of the fitted model for this record, mainly for the variables $\dot{\mathrm{V}}\mathrm{E}$, $\mathrm{VT}$, $\mathrm{TI}$, ${\mathrm{PACO}}_2$ and $\mathrm{HR}$, both in magnitude and behavior. Although $\mathrm{PM}$ and $\mathrm{PD}$ show some of the most significant differences in magnitude, their behaviors concerning the stimulus are similar, mainly for $\mathrm{PD}$, compared to the results of the traditional single-fit approach. Some variables that present the most significant differences for the predictive approach are related to the parameter values with the highest variations regarding those optimized for the record. Thus, the results of PM and PD can be attributed mainly to the value of ${\mathrm{I}}_{0,\mathrm{met}}$ and ${\mathrm{PaCO}}_{2,\mathrm{n}}$, and $\mathrm{BF}$ to the value of ${\mathsf{\lambda }}_2$ and ${\mathrm{P}}_{\max }$ (Table 4). For ${\mathrm{PAO}}_2$, the differences in the parameter values are not so significant, so the results can be mainly related to the effect of the systemic arterial pressures on the blood flows that enter the gas exchange system (Equations of Gas Exchange and Mixing, and Gas transport in the Supplementary Material).

Figure 4 shows the PE values obtained for the physiological model predictions regarding the third DB1 record. The obtained overall results indicate that the simulations related to the predictive approach are more suitable than those of the traditional single-fit approach. The highest errors associated with the traditional single-fit approach are for $\dot{\mathrm{V}}\mathrm{E}$ and $\mathrm{VT}$, in which the most statistically significant differences are obtained regarding the results of the model fitted to the record. It is also evidenced that the most significant PE’s contribution to the predictive approach is mainly related to the variables $\mathrm{PM}$ and $\mathrm{PD}$.

Figure 5 shows the comparison results of the steady-state model simulations against the third to the sixth DB2 records. The results obtained for each record are described below.

The simulations of the predictive approach related to the third DB2 record show less similarity to the model-fitted results than for DB1, despite the fact that they match the same record number. The most significant differences are presented for the variables $\mathrm{BF}$, $\mathrm{TI}$, $\mathrm{PS}$, $\mathrm{PM}$, $\mathrm{PD}$, ${\mathrm{PACO}}_2$ and ${\mathrm{PAO}}_2$, and are related to the parameter values used for the simulations (Table 5). The value of ${\mathsf{\lambda }}_2$, which is a weighting factor that relates the work of breathing during expiration, significantly affects the regulation of $\mathrm{BF}$ and $\mathrm{TI}$ (Equations of the Breathing Pattern Optimizer in the Supplementary Material). The lower magnitudes of the systemic arterial pressures are related to a lower predicted value of ${\mathrm{KE}}_{\mathrm{lv}}$, due to its effect on the left ventricular pressure signal (Equations of The Heart in the Supplementary Material). For ${\mathrm{PACO}}_2$ and ${\mathrm{PAO}}_2$, although the differences could also be related to the parameter values, especially ${\mathrm{K}}_2$ for ${\mathrm{PACO}}_2$, considering the high sensitivity of the gas exchange system to cardiovascular feedback (Equations of the Gas Exchange and Mixing in the Supplementary Material), it could be mainly attributed to the lower systemic blood pressures.

In the case of the fourth DB2 record, the predictive approach results showed improvements in the behavior and magnitude of the variables regarding the previous one, mainly for $\mathrm{TI}$, ${\mathrm{PACO}}_2$ and ${\mathrm{PAO}}_2$. The difference shown for the respiratory variables could be related to the high variation between the predicted and optimized parameters regarding the respiratory controller (Table 5). The difference between the $\mathrm{HR}$ slopes regarding the stimulus is mainly due to the prediction of higher sympathetic activity (${\mathrm{G}}_{\mathrm{T},\mathrm{s}}$) and lower parasympathetic activity (${\mathrm{G}}_{\mathrm{T},\mathrm{v}}$) (Equations of the Effectors for Reflex Control in the Supplementary Material), a result that is coherent considering the time trend of the previous records of these parameters. The improvement in the simulations of systemic blood pressure is due to the increased similarity between the predicted and optimized values of ${\mathrm{KE}}_{\mathrm{lv}}$. Regarding ${\mathrm{PACO}}_2$ and ${\mathrm{PAO}}_2$, considering that the differences in the related parameters are not significant, the improvement can be attributed mainly to the best predicted systemic blood pressure values.

As shown in Figure 5, the results from the predictive approach for the fifth DB2 record show more similar behaviors to those from the fitted model, mainly for the variables BF, TI, PS, PM, and PD. This can be attributed to a smaller difference between the predicted and optimized parameter values compared to previous records due to better predictions and the use of values of earlier optimizations that are not modified for the last fit (Table 5). Regarding the differences in magnitude between the predictive approach results and the fitted model, a higher $\dot{\mathrm{V}}\mathrm{E}$ value could be mainly related to a high dead space volume ($\mathrm{V}{0}_{\mathrm{dead}}$) (Equations of the Ventilation Controller in the Supplementary Material); the lower TI and the higher BF are related to a higher ${\mathrm{P}}_{\max }$ value [46]; a lower HR is due to a higher ${\mathrm{T}}_0$ (Equations of the Effectors for Reflex Control in the Supplementary Material); and the differences in ${\mathrm{PACO}}_2$ and ${\mathrm{PAO}}_2$ could be attributed to the variation in the concentration of hemoglobin (${\mathrm{C}}_1$) and to the difference in the blood and respiratory flows (related to $\dot{\mathrm{V}}\mathrm{E}$ and the systemic arterial pressures).

The predictive approach results for the sixth DB2 record present magnitudes more similar to those of the fitted model than those of the traditional single-fit approach. Regarding behavior, there was an increase in the variation regarding the results in previous records, which coincides with the increase in the differences between the predicted and optimized parameter values. The main differences occur for systemic arterial pressures and are due to the predicted value of ${\mathrm{PaCO}}_{2,\mathrm{n}}$ is higher than the optimized (Table 5). The predicted value for this parameter is similar to the optimized for the previous record, and its effect on many of the cardiovascular controller mechanisms is related to an inflection behavior of systemic arterial pressures at mid-stimulus levels (Equations of the Afferent Pathways, Blood Flow Local Control and CNS Ischemic Response in the Supplementary Material). Other major differences are related to the lower slopes and offset of BF and TI and could be related to a higher value of $\mathrm{n}$ and ${\mathrm{P}}_{\max }$ [46].

Figure 6 shows the PE results obtained for the model predictions regarding the DB2 records from third to sixth. The results for each record are described below.

The results for the third record corroborate the significant difference in the simulations related to the predictive approach compared to those of the fitted model. The above evidence, similar as for dynamic fitting results (Table 3), a high sensitivity of the predictive approach regarding significant changes in the variables evaluated. For this record, the significant change mentioned evidence a limitation of the strategy proposed for DB2 regarding the personalized approach (section 2.2.2. Experimental Data), considering a small number of samples used for the prediction (only two), that corresponds to a different subject from the one trying to predict.

The values presented for the fourth record evidence an improvement of the predictive approach, even presenting better overall results than the traditional single-fit approach. Although the distribution of the PE values is similar to the previous record, the improvement in the behavior of the variables allows results closer to those of the fitted model, mainly for the gas exchange variables.

The values for the predictive approach in the fifth record also show better results than the traditional single-fit approach, even for more variables than in the previous record. This consecutive improvement is related to the greater similarity regarding the magnitudes and behaviors of the variables (Figure 5), a consequence of the smaller overall difference between the predicted and optimized parameter values (Table 5).

The results for the sixth record show better results of the predictive approach than the traditional one for most of the variables evaluated. The differences obtained regarding the fitted model are related to that described for the steady-state results and were expected considering the long prediction time for this record (see Table 1).

4. Discussion

4.1. Dynamic Model Fitting

It is highlighted that the highest deviations and errors are mainly related to the respiratory variables, contrary to what is presented for the gas exchange variables (Table 3). These results are similar to those reported in the building and validation works of the case study model [37,47] and, therefore, could be mainly related to the predictability characteristics of the model rather than as a result of the method presented.

It was also evidenced that the variabilities and overall PE values were lower for DB2 (Table 3). These results are related to the characteristic dispersion of the DB1 records and show the better performance of the model under a personalized fitting perspective.

Regarding DB2, there were continuous decreases in the overall PE value between the first and second records and from the third to the fifth (Table 3). These results could be related to a sensitivity characteristic of the dynamic fitting strategy, which shows the records with significant changes in the cardiorespiratory response. Thus, considering the approximation proposed for DB2 (section 2.2.2. Experimental Data), the model presents better results between successive records of the same subject due to the use and accumulation of previously optimized parameter values (which do not change regarding the next fit). Still, when there is a significant change in the response of the system, which in this case is due to the variation of the recorded subject, the prediction error of the variables increases. This can not be observed for DB1, possibly due to the variability of the experimental data in its records.

4.2. Dynamic Parameter Modeling

The number of parameters that were not modeled regarding the third record was higher for DB2 (see Table 4 and Table 5). This evidence a significant change in the cardiorespiratory characteristics, which is related to the change of subject recorded according to the approximation implemented (section 2.2.2. Experimental Data). In the case of DB1, despite the variability of a popu.ation record, this effect is attenuated considering the fit regarding the mean trend of the data.

It was also observed that the number of parameters that did not satisfy the physiological justification criteria was higher for DB2, and the highest number is presented for its sixth record (see Table 4 and Table 5). The above also shows the sensitivity of the strategy in relation to significant changes in the records, which for DB1 are attenuated considering the fit regarding the mean trend of the experimental data. For DB2, these results are related to greater difficulty in modeling the temporal behavior of the parameters regarding the available factors and optimization records. This complexity increases, considering a longer prediction time for the sixth record (see Table 1).

It should be noted that as the number of records increases, so does the number of parameters to be modeled. This is because the optimized parameters are not necessarily the same between consecutive model fits. According to the incorporated static fitting strategy, the selection of the parameters contemplates the evaluation of the differences between the experimental data and the model predictions, so that the variations between characteristics or magnitudes promote the optimization of the most related parameter.

Although a greater number of model parameters is consistent, considering that all physiological characteristics change with time, this also implies a greater computational load, which may not be justified regarding the temporal trends of the optimizations. ${ \dot{\mathrm{V}}}_{\max }$ and ${\mathrm{VTCO}}_2$ regarding DB2 are examples of the above, and although their values are modified in some of the records, they do not change again for subsequent ones, or their variations are not representative (see Table 5).

4.3. Validation

Regarding the simulation comparisons for the third record of DB1, the higher overall PE value of the traditional single-fit approach evidence the change in the cardiorespiratory system over time and, therefore, the importance of a dynamic fit of the model to obtain more adequate predictions of the variables (Figure 4). Additionally, the greater similarity of the results related to the predictive approach with those of the fitted model, both in magnitude and behavior (Figure 3), indicates the usefulness of considering the temporal trend of associated factors and parameter optimizations to predict the future state of the modeled system.

The predictive approach results evidenced a greater sensitivity of the strategy regarding the personalized records (DB2) than with population records (DB1) (Figure 3 to Figure 6). This sensitivity is related to the significant variation of the variables to be predicted regarding their previous records and the few samples used for the prediction. Although the predictive approach considered the same number of samples regarding the third record of both databases, the variation in the temporal trend of the records of DB1 is less than DB2 due to the mean trend of the registered subjects. For DB2, the mentioned variation evidences a limitation of the approximation used for the personalized records (section 2.2.2. Experimental Data), since the records used for the prediction correspond to a different subject from the one attempted to predict.

As the number of temporary records evaluated for DB2 increases, the predictive approach exhibits better overall results regarding accuracy and behavior than the traditional single-fit approach (¡Error! No se encuentra el origen de la referencia. Figure 5 and Figure 6). This can be mainly attributed to two characteristics of the proposed strategy. Firstly, using a more significant number of samples allows us to obtain model parameters with more precise predictions regarding the optimized results (see Table 5). Secondly, using previously fitted parameters, some of which do not change between optimizations, allows us to get predictions with behaviors and magnitudes more similar to those of the model fitted to each record. The above confirms the advantage of dynamic fit regarding the traditional application of a fitted model in a single instant, in addition to evidence that a personalized approach is more appropriate than a population one.

Regarding the last DB2 record, despite the good overall result of the personalized approach, the difference in the behavior of some variables regarding the fitted model simulations was highlighted (Figure 5). Mainly, this difference was because using previous values of optimized parameters implies results with characteristics of previous fittings (see Table 5). In this sense, when evaluating a longer prediction time and, therefore, significant variations of the variables (due to the expected change of the modeled system in relation to age), the memory characteristic of the proposed strategy could be evidenced.

4.4. Application, Limitations, and Future Work

In this work, a dynamic fitting strategy of physiological models was proposed and evaluated regarding a model of the cardiorespiratory system. Experimental data from a longitudinal study of healthy subjects under aerobic exercise were used. The validation of the strategy consisted of comparing the model prediction results using the proposed strategy against those of the traditional approach of the population fit in a single instant of time. The comparison was made from a population and individualized perspective, and its objective was to demonstrate the importance of dynamic fitting of the physiological model.

In accordance with what is presented in this work, the proposed strategy can be considered more proper to represent the cardiorespiratory system. This is due to generally better results in terms of accuracy and behavior of the predicted variables, which validate the strategy both for the fitting of the model at each instant of time and for the prediction of the future state of the physiological system. The above evidence the importance of considering the temporal evolution of the physiological systems and associated factors such as habits, accidents, illnesses, clinical history, among others.

The results obtained also show the usefulness of the proposed strategy in relation to medical practice, mainly under the predictive, preventive, and personalized approach. Thus, it can be considered a tool that allows mathematical models to be suitable for applications of analysis and simulation of the temporal response of the physiological systems, both from a population and personalized perspective. In addition, it provides the ability to anticipate physiological behavior by predicting the future state of physiological systems.

This work also identified aspects that may limit the application of the proposed strategy. The sensitivity of the results regarding significant changes of variables to be predicted between records is highlighted, both for the dynamic fitting of the physiological model and for the prediction of the future state of the system. This aspect can be related to large time differences between consecutive records and specifically to the predictive approach with a short number of samples. The memory characteristic is also standing out, so the quality of the experimental data records can significantly influence the results. Finally, the availability of sufficient and adequate factors related to the temporal response of the system can also affect the results of the predictive approach.

The proposed strategy can have various future applications due to the good validation results. Considering that in this study, the strategy and the procedures that comprise it were generically described, its application regarding other physiological systems, stimuli, characteristics of subjects, and experimental recording conditions can also be evaluated. Works related to the temporal analysis and the prediction of subjects with different pathologies, clinical treatments, medication, or physical conditioning could be considered for the case study model.

Although the proposed methodology was adequate to validate the strategy, some aspects could be considered for future studies. These include the evaluation of the effect of time between records; the evaluation of different strategies for modeling the time trend of parameters; the identification, analysis, and evaluation of different factors related to the temporal response of the system to the predictive approach; and the definition and evaluation of criteria that allow defining which parameters justify being modeled based on the time trend of their optimizations.

5. Conclusions

This paper presents the description, application, and validation of a dynamic fitting strategy of physiological models under the consideration of a large number of parameters and a limited amount of experimental data. The strategy consists of different procedures that allow the dynamic fit of the model, the modeling of the temporal trend of their parameters, and the validation considering its usefulness to predict the future state of the modeled system. The results allow us to conclude that the proposed strategy provides the physiological models with the ability to be used in applications with a predictive, preventive, and personalized approach.