The mathematical model was based on the expected physiological response resulting from the interaction between the different variables involved in the medium-term blood pressure regulation. Several mechanisms are known to be involved in the control of blood pressure. However, to simplify the model and ensure identifiability, fundamental inputs that have more significant influence or represent the behavior of a set of variables that affect the circadian variation in blood pressure were used in this work. The chosen input variables included in the model are norepinephrine (NE), physical activity (A), and glycemia (glc), and the dependent variables are systemic vascular resistance (SVR), systolic blood pressure (SBP), diastolic blood pressure (DBP), and heart rate (HR).

Physical activity was chosen since previous works have established that it is one of the determinants of the variation in blood pressure and the variability of heart rate (Kario et al., 1999; Leary et al., 2000; Mansoor, 2002). Furthermore, greater night-time activity has been observed in non-dipper patients, making it one of the influencing factors in altering the circadian profile. However, it has been concluded that it cannot predict a non-dipper pattern by itself (Hermida et al., 2002b). Therefore, other variables related to the alteration of the circadian pattern of blood pressure were also considered. On the one hand, the plasma levels of norepinephrine exhibit a circadian variation, increasing during the day and decreasing during night time rest (Linsell et al., 1985; Candito et al., 1992). It is known that norepinephrine generates vasoconstriction and increases the cardiac frequency and output due to its action on adrenergic receptors. Furthermore, a high correlation has been shown between plasma norepinephrine and the outflow of the sympathetic nervous system (Goldstein et al., 1983), which is part of the most relevant blood pressure regulation mechanisms. For example, renin release has been shown to be influenced by sympathetic outflow (Zanchetti and Stella, 1975) and a direct relationship between norepinephrine levels and plasma renin activity has been observed in both healthy and hypertensive subjects (Beretta-Piccoli et al., 1980; Tuck et al., 1985). Therefore, by including this variable, we include its direct effect and, indirectly, other relevant regulatory mechanisms. Moreover, the variations of this catecholamine have been related to the circadian variation profile of blood pressure and a positive correlation has been observed between urine norepinephrine and the fall in blood pressure at night rest (Arita et al., 1996; Därr et al., 2016). On the other hand, it has been shown that a sharp increase in plasma glucose can decrease peripheral vascular resistance by 13%, increase heart rate by 14%, consequently increase 20% of cardiac output, maintaining unaffected blood pressure (Jern, 1991). Furthermore, blood glucose and the consequent increase in insulin have a positive effect on the sympathetic nervous system without altering plasma norepinephrine (Rowe et al., 1981). Additionally, differences in glucose profiles have been observed between dipper and non-dipper patients with type II diabetes, where non-dipper patients have been shown to have greater insulin resistance (Pistrosch et al., 2007). Differences have also been observed between dipper and non-dipper patients with essential hypertension, where those with a non-dipper pattern show greater insulin resistance (Chen et al., 1998).

The equations were formulated based on the dynamics of the interactions described in the scheme of Figure 1. Direct relationships were modeled by means of first-order kinetics while higher-order interactions involving activating or inhibiting effects were modeled using Hill functions. For the sake of parsimony, a Hill coefficient equal to one was considered when possible. The effects described between variables do not necessarily imply a direct mechanism of activation or inhibition but express the expected physiological response that each of the involved variables exerts on one another. For example, the relationships between physical activity and heart rate and between physical activity and SVR are complex and include several components and regulation systems; however, we have simplified these relationships considering only the initial and final expected effects. Specifically, physical activity generates a greater demand for oxygen and nutrients in the skeletal muscle; this causes the muscle to release vasodilators such as nitric oxide, which reduces vascular resistance and blood pressure in that area. Then, the sympathetic nervous system acts to regulate blood pressure by increasing the vascular resistance in the abdominal viscera and increasing heart rate and stroke volume, thereby increasing cardiac output and, consequently, blood pressure (Nobrega et al., 2014; White and Raven, 2014). Despite the increase in vascular resistance in the visceral abdominal area, its decrease at the peripheral level generates a general reduction in vascular resistance (Casey and Joyner, 2012). The relationships described above are reflected in the SVR and HR equations (Equations 1, 3) by the expressions ki1ki1+A and Aki4+A, which represent the inhibition and activation exerted by physical activity on the SVR and HR variables, respectively. Therefore, the parameters k_i1 and k_i4 are constants that indicate the relevance of the effect of the activity in the changes of SVR and HR, respectively. The larger the value of k_i1, the smaller the regulation (inhibition) of the SVR due to the activity. Analogously, the larger the value of k_i4, the smaller the regulation (activation) of the HR due to the activity.

Several mechanisms regulate vascular resistance. Some regulatory molecules are capable of generating a short-term response; for example, the action of catecholamines, like norepinephrine, on alpha-adrenergic receptors causes vasoconstriction (Takeda and Maemura, 2011; Tank and Lee Wong, 2015). Moreover, as mentioned above, norepinephrine has a high correlation with the outflow of the sympathetic nervous system (Goldstein et al., 1983), which is related to vasoconstriction. On the other hand, insulin secretion, triggered by increased blood glucose, generates vasodilation (Muniyappa et al., 2007). Hence, the influence of both norepinephrine and glucose was included in the model on the expression NEglcn+NE. In a previous version of the model, the effect of glucose and norepinephrine was modeled by additive inhibition and activation expressions, respectively. However, a better fit to the data was obtained by relating them in a single expression, suggesting a synergistic effect. Therefore, the parameter n is a constant that represents the impact of glucose on vasodilatation and k₁ represents the basal rate of vasoconstriction (without the influence of the activity and the glucose). Finally, the expression SVRki2+SVR*SVR, that represents the regulation mechanism of vascular resistance of blood vessels produced when blood flow to an organ decreases generating metabolic-induced vasodilatation (Nobrega et al., 2014), was also included in Equation (1). Consequently, the k_i2 value is related to the self-regulatory capacity of SVR.

Equations (2) and (4) describe the change of blood pressure over time, which depends directly on SVR because, for a given blood flow, an increase of vascular resistance generates an increase in the pressure on the blood vessels' walls (SBP and DBP). In addition, the influence of BP over HR is captured by the expression HRki3+HR; a higher HR (together with a higher stroke volume) generates an increase in cardiac output, increasing blood flow, which also induces an increase in SBP and DBP. Due to this, the parameter k_i3 is related to the ability of HR to modify the value of BP. Besides, heart rate is regulated by the autonomic nervous system, which activity at this level is regulated, in part, by the activation of baroreceptors, that are able to detect changes in BP (Guyenet, 2006). This relationship between BP and HR is known as baroreflex regulation, and it was included in Equations (2)–(4). Baroreflex regulation allows regulating the values of HR and BP, increasing the value of HR when BP decreases and vice versa. Finally, a constant rate of increase of HR or baseline is expressed by k₂. Moreover, parameters k₃ and k₄ were incorporated into the model to ensure dimensional homogeneity.

Thus, the following system of four differential equations and seven parameters was obtained:

Some variables, as changes in fluid volume, were not considered due to unavailability of data; however, these variables were enough to capture the general pattern of the circadian rhythm of blood pressure.

Experimental data on physical activity, heart rate, and blood pressure (systolic and diastolic) for dipper and non-dipper subjects measured every 30 min over 24-h were extracted using the Engauge Digitizer software from the work of Hermida et al. (2002b). The data correspond to the average of 130 male patients (53.7 ± 14.0 years of age).

The norepinephrine data for dipper patients was extracted from the work of Linsell et al. (1985) and van Dijk and van Loon (2015) and the values for non-dipper patients were simulated by increasing by 15% during the night rest the values obtained from dipper patients, according to the work of Därr et al. (2016), who reported a significant difference between day-night urinary norepinephrine excretion ratio between dippers and non-dippers.

The glycemia data for dipper patients was extracted from the work of van Dijk and for non-dipper patients, these values were increased by 4.7, 13.3, and 8.1% 2 h after each of the three meals: breakfast, lunch, and dinner, respectively, according to what was reported by Pistrosch et al. (2007).

The systemic vascular resistance (SVR) data to fit the model was obtained using the following equation (Coats et al., 1989):

assuming an average stroke volume (SV) for subjects with essential hypertension of 79.5 ml (Messerli et al., 1983) and replacing the systolic blood pressure (SBP), diastolic blood pressure (DBP) and heart rate (HR) of dipper and non-dipper subjects from the work of Hermida et al. (2002b). The value of the mean right atrial pressure for essential hypertension patients (5 mmHg) (Ferlinz, 1980) was used instead of the central venous pressure (CVP), since, as previously reported, both values are similar given the low resistance of large vessels (Magder, 2005). The standard deviation of SVR was calculated from the standard deviation of the variables SBP, DBP and HR.

The data were organized in such a way that the simulation begins at 0 a.m. and wake-up time is 8 a.m. for all the data.

The ODE-based model was developed using previous knowledge about the main regulatory mechanisms of blood pressure, detailed in the Mathematical Model and Experimental Data sections, and was implemented in MATLAB 2018a. Each model's parameters were fitted to the experimental data employing global optimization using a scatter search method (MEIGO) (Egea et al., 2014). Scatter Search is an evolutionary meta-heuristic able to solve combinatorial and non-linear optimization problems, which allows fitting the model to the experimental data in order to obtain the values of the parameters that minimize an objective function, i.e., given a certain range for each parameter, MEIGO iteratively searches for the combination of the parameters that minimize the error between the experimental data and the model predictions. To integrate the model, the solver ode23s of MATLAB was chosen and the objective function is expressed as a sum of squares:

where N is the number of experimental data per variable, k is the number of variables, and y_{exp i} is the experimental value of the i variable, and y_{model i}. is the value predicted by the model. Twenty-four hour data were used to adjust the model.

Data from both types of subjects were used independently to estimate the parameters. Then, according to the values obtained individually, some parameters were established in common for both types of subjects and the remaining parameters were re-fitted.

The analysis of residuals was carried out using Equation (7), where the residuals, R(t_n), are calculated using the experimental values (HR, SBP and DBP), y_expNCi(t_n), and the values predicted by the model, y_modeli(t_n), at each time, t_n.

To evaluate the residuals, a QQ plot and a histogram were plotted using MATLAB's qqplot and hist functions. Furthermore, to corroborate these results, the Shapiro Wilks test of normality was used using MATLAB's swtest function with an alpha value of 0.05.

Local sensitivity analysis allows quantifying the impact of the local variation of a parameter on the value of each of the output variables (in contrast with global sensitivities that explore the entire parameter space). In this work, the local sensitivity of each variable y_i with respect to each parameter p_j evaluated for the best parameter set was computed using the SENS_SYS third-party MATLAB function (Garcia Molla, 2019):

Then, normalization was performed using the following equation:

Srel_ij(t_n) was calculated for each parameter (j), variable (i), and time point (t_n) obtaining relative sensitivity trajectories for each parameter and each variable (SVR, SBP, HR, and DBP). The analysis of these results is complex, due to the large number of sensitivity values (for each time and variable). Therefore, a sensitivity index for each parameter was calculated using the following equation (Brun et al., 2001), which transforms the relative sensitivity trajectories corresponding to each parameter into positive values by squaring and averaging them across the different time points and variables:

where N is the number of experimental data per variable, k is the number of variables, and y_{model in} is the value n of the i variable predicted by the model.

To carry out an identifiability analysis, the Fisher information matrix (FIM) was obtained from the local sensitivity matrix S(t_n) and the covariance matrix Q(t_n) (Ljung, 1999):

where Q(t_n) is a 4-by-4 diagonal matrix calculated using the standard deviation of the data for SBP, DBP, HR, and the estimated standard deviation for SVR at each time point and S(t_n) is a 7-by-4 matrix obtained from the local sensitivity values of each of the four variables with respect to each of the seven parameters at each time point (Equation 8).

According to the Cramér–Rao theorem, the inverse of the FIM represents an approximation of the covariance matrix of the error of the objective estimator of minimum variance (Ljung, 1999). In this case, the inverse of the FIM is a 7-by-7 matrix where each element represents an approximation of the parameter estimation error covariance between the parameters j and h (σjh2):

Therefore, the diagonal of the inverse of the FIM (σjj2) is an approximation of the variance of the parameters. Thus, assuming a normal distribution, the 95% confidence intervals (CI) of a parameter can be approximated by p ± 1.96 σ_jj (Walter and Pronzato, 1997).

Moreover, the correlation between parameters (κ_jh) can be calculated using the following equation (Ljung, 1999):

The correlation matrix measures the interrelationship between the parameters and gives an idea of the compensation effects of changes in the parameter values on the output variables. If two parameters are highly correlated, a change in the model output caused by a change in a model parameter can be (nearly) compensated by an appropriate change in the other parameter value. This prevents the parameters from being uniquely identifiable even if the model output is very sensitive to changes in the individual parameters. A singular FIM indicates the presence of unidentifiable parameters, and correlations between parameters that are >0.95 may lead to singular FIM.

After several iterations in the search for equations representing the data, it was possible to adjust the same model to both dipper and non-dipper subjects. The experimental data and the model predictions for the dipper and non-dipper subjects can be seen in Figures 2, 3, respectively. The estimated values and 95% CI of the parameters are shown in Table 1. The parameters k₃ and k₄ were fixed to one due to its high correlation with k_i3 and k₂, respectively; they were kept in the model to ensure dimensional homogeneity. The table shows the common parameters between dipper and non-dipper subjects (k₁, n, k_i2, k₃, and k₄) and the parameters with different values for each type of subject (k_i1, k_i3, k_i4, k₂).

The parameters k₁, n, k_i1 and k_i2 are part of Equation (1), of which only k_i1 presents a different value between both types of subjects, suggesting a difference in the SVR response to physical activity in non-dipper subjects. On the other hand, in Equation (2), a greater value for parameter k_i3 was obtained for dipper subjects. This can be due to the range of HR values that is slightly smaller in non-dipper subjects (minimum value is lower in dipper subjects). In Equation (3), parameters k_i4, k₂ have higher values for dipper subjects, which can be explained by a higher activation baseline that requires a smaller increase in each time interval and a greater influence of physical activity on HR. The 95% confidence interval values indicate that all parameters are statistically significant. The value of the objective function was 0.0015; namely, the mean of the normalized error for each data point is 0.15%.

The QQ plot and histogram of the residuals (Supplementary Figures 1, 2) show a distribution close to normal. This was verified using the Shapiro Wilks test, which indicates that, with a 95% confidence, the residuals come from a normal distribution (p-value = 0.762).

The results obtained for the sensitivity analysis are summarized in Figure 4. Parameters n and k_i1 have higher sensitivity for non-dipper subjects than for dipper subjects, parameters k₁, k_i2, k_i3, and k₂ have similar sensitivity for dipper and non-dipper subjects, and parameter k_i4 has higher sensitivity for dipper subjects than for non-dipper subjects. Parameters k_i1 and k_i4 have different sensitivity and they are related to physical activity, which indicates a different response to physical activity for non-dipper subjects. Additionally, the higher sensitivity of parameter n can be due to higher glycemia (glc) values for non-dipper subjects.

The results obtained for the identifiability analysis are summarized in Figure 5. The FIM matrix is not singular; therefore, it was possible to calculate its inverse and obtain an approximation of the correlation matrix for parameters of dipper subjects (Figure 5A) and non-dipper subjects (Figure 5B). Both types of subjects showed a high correlation between k₁ and k_i2, but only dippers subjects had a correlation value >0.95. The correlation between the rest of the parameters was <0.95. Therefore, the confidence intervals of the estimated parameters indicate that the model is capable of representing the data of dipper and non-dipper subjects with essential hypertension.