## Abstract

Heart rate variability biofeedback intervention involves slow breathing at a rate of ∼6 breaths/min (resonance breathing) to maximize respiratory and baroreflex effects on heart period oscillations. This intervention has wide-ranging clinical benefits and is gaining empirical support as an adjunct therapy for biobehavioral disorders, including asthma and depression. Yet, little is known about the system-level cardiovascular changes that occur during resonance breathing or the extent to which individuals differ in cardiovascular benefit. This study used a computational physiology approach to dynamically model the human cardiovascular system at rest and during resonance breathing. Noninvasive measurements of heart period, beat-to-beat systolic and diastolic blood pressure, and respiration period were obtained from 24 healthy young men and women. A model with respiration as input was parameterized to better understand how the cardiovascular processes that control variability in heart period and blood pressure change from rest to resonance breathing. The cost function used in model calibration corresponded to the difference between the experimental data and model outputs. A good match was observed between the data and model outputs (heart period, blood pressure, and corresponding power spectral densities). Significant improvements in several modeled cardiovascular functions (e.g., blood flow to internal organs, sensitivity of the sympathetic component of the baroreflex, ventricular elastance) were observed during resonance breathing. Individual differences in the magnitude and nature of these dynamic responses suggest that computational physiology may be clinically useful for tailoring heart rate variability biofeedback interventions for the needs of individual patients.

- heart rate variability biofeedback
- respiratory sinus arrhythmia
- adaptability
- baroreflex
- power spectral density
- respiration

the ability to adapt physically and psychologically to the environment depends on bidirectional communication between the body and brain. The efferent path sends information from the brain's autonomic control circuitry to the heart (13, 48, 65). In parallel, the afferent path sends information from the heart to the brain through the baroreceptors (50, 67). This flow of information within what has been termed the baroreflex feedback loop, or baroreflex arc, generates variability in physiological functions, such as heart period and blood pressure (8, 16, 67). This baroreflex feedback loop has been characterized in terms of neurophysiology, functional anatomy, and molecular mechanisms in rodent and primate models as well as human research (9, 11, 16).

In healthy organisms, heart period changes continuously, reacting and modulating arousal to internal and external challenges. This process is usually referred to as heart rate variability (HRV), the variability in consecutive R-wave to R-wave intervals (RRI) of the electrocardiogram (ECG). In humans, indexes of resting-state HRV are robust biomarkers of both physical and psychological health and illness (3, 7, 42, 65, 66). Lower resting HRV has been linked to stress, anxiety, depression, reduced immune response, and all-cause mortality (1, 21, 67, 78, 81). Further, HRV reacts to physical, pharmacological, cognitive, and emotional challenges. Reactive HRV captures integrated brain and cardiovascular control system response in real time and can serve as a gauge of adaptability (5–7, 10, 40, 41, 68, 74, 76).

Behavioral intervention methods, such as HRV biofeedback, build on these naturally occurring changes in HRV to increase the efficacy of the baroreflex feedback loop to improve health and reduce clinical symptoms (19, 24, 28, 31, 55). HRV biofeedback uses slow breathing as a challenge to instigate HRV reactivity. It teaches individuals to induce high-amplitude oscillations in heart period by breathing at a rate of ∼6 respiratory cycles/min (i.e., ∼0.1 Hz), which is considerably slower than typical breathing. This breathing rate has been referred to as resonance breathing (77) because it synchronizes the cardiac oscillations driven by respiratory sinus arrhythmia (i.e., the phenomena of heart acceleration with inhalation and deceleration with exhalation) with the cardiac oscillations driven by the baroreflex (i.e., which links heart acceleration/deceleration to corresponding decreases/increases in blood pressure; 60, 73), thereby maximizing HRV.

The results of clinical efficacy studies support the use of 6–10 weekly sessions of HRV biofeedback training and daily practice to reduce clinical symptoms in disorders characterized by reduced baroreflex sensitivity (34) and difficulty in modulating emotion and arousal, such as asthma (29), fibromyalgia (19), depression (24, 44, 55), and posttraumatic stress disorder (63, 83). HRV biofeedback/resonance breathing also decreased blood pressure and increased baroreflex sensitivity in prehypertensive (33) and hypertensive (54) persons and improved cognitive and emotional regulation in persons with chronic brain damage (26). In addition to these positive cumulative effects, improvements in reaction time and cognitive performance (49), and reductions in performance anxiety (80), have been reported after a single brief HRV biofeedback session. Growing evidence for the clinical benefit of HRV biofeedback is noteworthy because of its public health potential for dissemination through accessible technologies (e.g., portable devices, mobile phone applications) and, unlike most pharmacological interventions, there are no known negative side effects associated with its use. These observations support the need to better understand the dynamic cardiovascular mechanisms of resonance breathing and the extent to which different subgroups of individuals may expect to experience clinical benefit from this intervention.

The purpose of this study was to investigate how cardiovascular processes change when challenged by resonance breathing using a computational physiology model of the human cardiovascular system. Many mathematical models of the human baroreflex and cardiovascular system have been proposed. Existing models describe the mechanics of the blood flow (61) and include a simple beat-to-beat model consisting of a set of difference equations and a time step equal to one heart period (14), as well as several heart-beat nonlinear models based on a system of differential equations (36, 45, 46, 53, 69–71). Lumped parameter models of the circulatory system, which include both respiration and baroreceptor feedback, have also been developed for understanding congenital heart disease (47) and coronary artery disease processes (25, 27). These models significantly advance our understanding of the system-level responses of the cardiovascular system. The current paper focuses on the mathematical description of the human baroreflex and cardiovascular system starting from the validated model developed by Ursino and Magosso (36, 69–71), which includes a pulsating heart and accurate description of blood circulation and the baroreflex. The model consists of a system of delay differential equations and is capable of simulating high-resolution blood pressure and heart period as a function of time.

Computational physiological modeling affords the unique ability to analyze the interdependent nature of system changes in cardiovascular and cardiorespiratory states, flows, and pressures, which is typically unobservable in the human laboratory, at least when noninvasive measures are used. The present study experimentally assessed cardiovascular control systems in humans at rest and in reaction to a slow breathing manipulation that served as a proxy for HRV biofeedback by maximizing cardiac oscillations. We then assessed how closely the computational model approximated the experimentally assessed data, namely, heart period and blood pressure at rest and during resonance breathing, and examined the importance of the derived parameters within and between individuals. Our primary aim was to understand the changes that arise from interactions among dynamic cardiovascular processes to identify key parameters (i.e., the processes) that change significantly during resonance breathing compared with a restful state. The predictive validity of these changes can be tested in future studies to determine whether acute response to resonance breathing is a prognostic indicator of the effectiveness of HRV biofeedback at the level of the individual or subgroups of individuals.

## METHODS

#### Participants.

Data for model validation were provided by 12 men and 12 women who were healthy college students between 21 and 24 yr of age (mean age = 21.7 ± 0.9 yr). They were participants in an experiment that was designed *1*) to develop a computational physiology approach to model how cardiovascular processes change when the baroreflex mechanism is challenged and then *2*) to determine how alcohol adversely affects cardiovascular processes at rest and during challenge. The present study addressed the first aim using a breathing challenge paced at 0.1 Hz (i.e., 6 breaths/min).

The racial and ethnic composition of the sample was 67% White, 16.5% Asian, 12.5% Black, and 4% Other or Mixed, with 12.5% Hispanic origin. Exclusion criteria were history of a psychiatric, neurological, or substance use disorder or treatment, weekly use of illicit or prescribed drugs, prenatal exposure to drugs of abuse, suicidal ideation, a medical condition or medication that can compromise interpretation of physiological assessments (e.g., asthma, diabetes, cardiac conditions), pregnancy (based on urine analysis), and body weight 20% more than or less than average according to Metropolitan Life Insurance norms (1983). In addition, if abnormal cardiac rhythm was detected during baseline recording, we excluded their data from analysis. All participants provided written informed consent and were reimbursed for their time. This study was approved by the Rutgers University Institutional Review Board for the protection of human subjects involved in research.

#### Experimental procedure.

Each participant completed one laboratory session that started between 10 a.m. and 2 p.m. to minimize biological circadian variations. They were instructed to fast for 4 h following a low-fat meal prior to reporting to the experiment. Physiological assessment was completed in a sound-attenuated, dimly-lit testing room. Participants were seated in a comfortable chair located 1 m in front of a LCD TV screen, and physiological sensors and electrodes were attached.

Participants completed two, 5-min low-demand baseline tasks (B1 and B2), in which a rectangle, presented in the center of the TV screen, changed color every 10 s. They were asked to silently count the number of blue rectangles. This baseline task provided a more standardized and reproducible baseline condition than an uncontrolled resting state (“Plain Vanilla” task, 22). A juice beverage was administered to participants between B1 and B2. This manipulation served as a control for the other experimental groups (not included in this study) that consumed an alcohol or placebo beverage. Data from both B1 and B2 were modeled because participants may become more relaxed (i.e., vagally influenced) during B2 after they have acclimated to the testing situation (7). Due to an error in data collection, B1 data were not available for one participant. Participants next completed a 5-min resonance breathing task (6P) (30, 33, 73, 75), during which they breathed at a rate of ∼6 breaths/min following a visual pacer (Easy Air, Biofeedback Foundation of Europe, Montreal, Canada) presented on the TV screen. They inhaled as the pacer moved vertically up and exhaled as it moved down. A brief practice preceded the 6P task wherein participants were instructed how to avoid hyperventilation by not breathing too deeply.

#### Computational model.

The present model was developed from Ursino and Magosso's mathematical model of the human cardiovascular system (36, 69–71). This model includes mathematical descriptions of a pulsating heart, as well as the mechanics of blood flow (61) and baroreflex activation (36, 69–71). It consists of a system of delay differential equations and was adapted in this study to be applicable for repeated measurements on the same person under different conditions (resting state, resonance breathing).

The model includes more than 90 parameters (see appendix a, Tables A1 and A2) and 21 states (pressures, flows, volumes, resistances, and elastances; see appendix a Table A3). Twenty-one delay differential equations reflect conservation of mass and balance of forces at arteries and veins, as well as delayed physiological responses to vagal and sympathetic neural activity. This allows for simulation of high-resolution blood pressure and heart period as a function of time. We modified the Ursino and Magosso model (36, 69–71) to use experimentally derived respiration period as a model input. In addition, we set external noise (parameter *A*; see appendix a, Table A2) from the Ursino and Magosso model to zero because of the noise in the respiratory input used in our model.

Although human respiration can be simulated (23), it strongly depends on the environment and is uncertain (i.e., difficult to predict). Thus the present model represents respiration as a dynamic input in the cardiovascular control processes that generate HRV. By using respiration rather than blood pressure as a model input, the closed-loop circuit between heart period and blood pressure was captured in its entirety. This approach (36, 69–71) coupled with our use of a respiratory input extends prior approaches that modeled heart period dynamics as a function of blood pressure, without feedback from heart period to blood pressure.

As detailed in the following section, we selected as model output the cost function that takes into account power spectral densities and time averages of several observables, such as heart period. Instead of doing brute-force optimization on the cost function with over 90 parameters in the model, we used the following procedure: *1*) an initial sensitivity analysis was performed to select the most important parameters to tune, and *2*) optimization of these most important parameters was performed to minimize the cost function. A great deal of economy in computing efficiency was gained by using this procedure. Furthermore, it provided additional insight into most sensitive “knobs” of the system.

#### Model parameterization.

Parameterizing models with oscillatory outputs is difficult due to their potentially strong sensitivity to phase. Even small changes in initial conditions or external disturbances can play a big role in the phase of the system response (37). For this reason, it is not advisable to attempt to precisely match phases of multiple outputs, such as modeled and experimental heart period. Here, we adopted an approach to model parameterization and validation, suggested in Ref. 39, in which phase-independent quantities are used as outputs. Specifically, we used the time averages and power-spectral densities (PSDs) of physiologically relevant outputs. The PSD of heart period was used to determine the contributions of sympathetic and parasympathetic (vagal) nervous system activities to heart period dynamics. The power spectrum is classically divided into three spectral components: a very-low-frequency component, which may be associated with humoral and thermal control or vasomotion; a low-frequency component, potentially associated with both sympathetic and vagal activity; and a high-frequency component, associated mostly with vagal activity (64). High-frequency HRV is often considered to be synonymous with cardiorespiratory influences such as respiratory sinus arrhythmia (e.g., 18); however, respiration also plays a role in the genesis of low-frequency HRV in the heart period power spectrum, particularly when respiration period is slowed to ∼0.1 Hz (30, 76). Greater low-frequency power in respiratory signals is linked to greater low-frequency power in model-predicted heart period (82).

Figure 1 shows a schematic overview of the model, wherein parameters were estimated for each participant individually, at rest and during the breathing challenge, by finding the sensitive model parameters that provided the closest match between the model outputs and the experimental data collected during B1, B2, and 6P. This was accomplished by using several runs of the mesh adaptive direct search algorithm that accommodated discontinuous cost functions (4). We used GoSUM software (2) to minimize the cost function with respect to sensitive model parameters. GoSUM is a software program for global sensitivity analysis, optimization, uncertainty analysis, and model reduction. It employs kernel-based algorithms to accurately learn various models either directly from data or from a model executable file. Once an input-output model is learned in GoSUM using support vector regression (56) and support vector machines (for continuous and binary outputs, respectively), GoSUM provides tools for the identification of parameters responsible for model variability, quantification of global uncertainty in model outputs, and optimization of an objective function in the presence of equality and inequality constraints. GoSUM also has routines for performing direct optimization, without using the learned input-output model.

The parameterization or calibration process included the following steps. First, the initial estimate for blood volume was calculated for each participant based on his/her height (in m) and weight (in kg) (43). For female and male participants, blood volume (V_{blood}, in liters) was calculated by using the following formulas:

and

respectively. The following parameters (see appendix a for all definitions) in the Ursino and Magosso model were rescaled for humans by multiplying them by V_{blood}/(5.3 liters): *C*_{pa}, *C*_{pp}, *C*_{pv}, *C*_{sa}, *C*_{sp}, *C*_{ep}, *C*_{ev}, *C*_{sv}, V_{u,sp}, V_{u,ep}, V_{u,pp}, V_{u,pv}, V_{u,sv,max}, V_{u,sv,min}, V_{u,ev,max}, and V_{u,ev,min}. The following parameters were rescaled for humans by dividing them by V_{blood}/(5.3 liters): *R*_{pa}, *R*_{pp}, *R*_{pv}, *R*_{sa}, *R*_{ev}, *R*_{sv}, *R*_{sp,max}, *R*_{sp,min}, *R*_{ep,max}, *R*_{ep,min}, *L*_{pa}, and *L*_{sa}.

Second, the duration of each breath (i.e., respiration period) was extracted from experimental respiration data. Thoracic and abdominal pressures for each breath (P_{thor} and P_{abd}) were modeled as described in the appendix b. Even though the amplitude of each breath was not constant, we found that the constant-amplitude approximation was sufficient for model parameterization. In contrast, the irregularity of respiration period significantly contributed to both high- and low-frequency oscillations in heart period and blood pressure (model outputs). The SD of experimentally derived respiration period for each participant in each task is presented in Fig. 2, *top*. Note the significant variability in respiration period across participants (horizontal) and tasks (vertical). The relative uncertainty of the experimentally derived respiration period (SD divided by the mean) is shown in Fig. 2, *bottom*. As expected, respiration period for the resonance breathing task (6P) had very small relative uncertainty (<10%).

Third, we derived cost functions for each task (B1, B2, and 6P). Initial cost functions were built based on conceptual knowledge of how respiration influences cardiovascular dynamics under different breathing conditions (18, 30, 76, 82). Cost function terms and coefficients were then refined iteratively. Due to processes not included in the current model (e.g., metabolic, vascular, neurological), it was not always possible to closely match the experimental and modeled time trace of blood pressure (i.e., the model peak occasionally appeared slightly before or after the actual data peak), indicating the phase sensitivity of the process as discussed above. Therefore, the cost function included PSD of the heart period and systolic and diastolic blood pressure, which were derived from the model output (blood pressure) using the physiological formulas shown in appendix c. Heart period (HP_{model}) was calculated as the length of time between the minima of the blood pressure (Psa_{model}) curve. Systolic pressure (SAP_{model}) was defined from the maxima of the Psa_{model} curve, and diastolic pressure (DAP_{model}) was defined from the minima of the Psa_{model} curve. HP_{exp}, SAP_{exp}, and DAP_{exp} were obtained in an analogous way from the experimental data (see appendix c). The cost function was designed to match model output and experimental data in terms of heart period and the power in the very-low-frequency, low-frequency, and high-frequency bands, as well as systolic and diastolic pressure and their PSDs. Coefficient weighting of heart period spectral power in different frequency bands was varied by task because respiration has a major influence on high-frequency oscillations in heart period at rest (B1, B2), whereas it has a greater influence on low-frequency oscillations in heart period during the resonance breathing task (6P). The terms that were dependent on the PSD of systolic and diastolic blood pressure had a small effect on the optimized parameters, and thus the weights of these terms were an order of magnitude smaller than the heart period-based terms.

The following cost functions for tasks B1, B2, and 6P, respectively, were used:

where

is the relative difference (error) between the PSDs of the model (calculated) and experimental heart periods in the very low-frequency band;

is the relative PSD error for the low-frequency heart period band;

is the relative PSD error for the high-frequency heart period band;

is the relative error in mean heart periods;

is the relative error in mean systolic blood pressure;

is the relative error in mean diastolic blood pressure;

is the relative PSD error for systolic arterial pressure;

is the relative PSD error for diastolic arterial pressure; *F*_{L1} = 0.002 is the lower limit of the very low-frequency band; *F*_{L2} = 0.05 is the upper limit of the very low-frequency band and the lower limit of the low-frequency band; *F*_{H1} = 0.15 is the upper limit of the low-frequency band and lower limit of the high-frequency band; *F*_{H2} = 0.50 is the upper limit of the high-frequency band; and *T*_{0} = 300 s is the duration of each task.

To set initial conditions for each of the 21 states in the model, we considered the largest delay in the HRV baroreflex closed-loop system, which is ∼5 s and corresponds to the time between a shift in blood pressure and the corresponding heart period change. Thus ∼5 s of experimental blood pressure (P_{sa}) time series data (for each of the three tasks) was supplied to the solver as historic values. The initial conditions for the remaining 20 states (see appendix a, Table A3) for each of the three tasks were found as follows. First, a reasonable guess was taken for the initial conditions. Then 5 min of model outputs were simulated, and the second half of the simulation was analyzed to find the smallest difference between the minima of high-resolution blood pressure from the model and from the experiment. For each task, the time moment corresponding to the smallest difference was chosen, and the remaining 20 states at that time moment were saved as initial conditions. Finally, each simulation used in the cost function minimization was started from a time moment during the first heart period cycle (after a 5 s history) when the experimental blood pressure was at its minimum. The above procedure to find initial conditions was found to be robust and resulted in a good match between the model and experiment from the beginning of the simulation. We assume that, without noise, the model in our paper is globally stable to an attractor (52). This attractor may not be simple in nature (i.e., not periodic); rather, its spectrum may be broad and continuous (38). However, for such systems, it is not the exact initial conditions per se that are important for computations of the spectrum, but the amount of time the system needs to converge to the attractor. Our computations of the spectrum showed that it did not substantially change when different initial conditions or even longer initial times were considered.

The estimated initial values for parameters in the model were set as follows: *T*_{min} and *T*_{max} were equal to the minimum and maximum heart period, respectively, during the corresponding task. P_{san} was set equal to the average experimental blood pressure for the corresponding task. V_{blood} was set to the estimated blood volume and the initial estimates for all remaining parameters were set to the values given in Refs. 69 and 71 and rescaled according to V_{blood} as described earlier.

Parameter values were then optimized to allow individualization of the model. As described in detail in results, we used a novel two-part optimization procedure. First, we started by using global sensitivity analysis to identify the most important input parameters that affected the outputs. This was achieved by sampling the space of all parameters, learning an analytic model, and performing sensitivity analysis on that model. The full set of input parameters is specified in appendixes a and b. Second, we optimized the reduced set of parameters with respect to the cost function described above, keeping the rest of the input parameters at their initial estimated values. For each participant, this process yielded a unique set of key parameters to the model. The reduced set of key parameters is specified below in *Model optimization*.

## RESULTS

#### Global sensitivity analysis.

Global sensitivity analysis was performed for all subjects and all tasks in order to select the parameters for model optimization. Specifically, we used variance-based global sensitivity analysis based on Sobol global sensitivity estimates and obtained analytical model representation (20, 32, 57) (see appendix d). All 80 nonconstant model parameters (see appendix a, Table A3) were assumed to be uncertain with a uniform distribution and a range of ±10% from the initial value. A similar pattern of relative parameter importance was observed for all 24 participants for all tasks. The following 12 parameters that were among the upper quartile of parameters across the B1, B2, and 6P tasks and that were conceptually relevant to baroreflex functioning were chosen as the optimization parameters: *k*_{E,lv}, V_{blood}, *G*_{aTv}, *G*_{pTv}, *G*_{aTs}, V_{Ln}, *C*_{sp}, P_{0,lv}, *E*_{lv,max}, *E*_{lv,min}, *R*_{ep,max}, and *R*_{ep,min}.

#### Model optimization.

The optimization was performed over the 12 parameters that were selected based on the global sensitivity analysis. To optimize this black-box model, we used the NOMAD optimizer, which is a C++ implementation of the mesh adaptive direct search algorithm and is directly incorporated in GoSUM. The chosen set of optimization parameters (and the range around the initial value used) were V_{blood} (±40%), *V*_{Ln} [±50% (B1, B2), ±100% (6P)], *C*_{sp}, P_{0,lv}, *E*_{lv,max}, *E*_{lv,min}, *R*_{ep,max}, *R*_{ep,min}, *G*_{aTv}, *G*_{aTs} (±100%), and *k*_{E,lv}, and *G*_{pTv} (±300%). All of these parameters were constrained to be greater than zero during this step. The rest of the parameters were kept at their initially estimated values. For each participant, the optimization was repeated five times for each task without any modifications between runs to account for the stochastic nature of the optimization algorithm (doing more than five times was computationally prohibitive). The maximum number of model evaluations for each step was set to 5,000. The solution with the smallest cost was selected as the final solution. The personalized parameters for each participant are given in appendix e.

Figure 3*A* shows the distributions of the optimized parameter values. The considerable spread of parameter values across participants highlights the importance of model personalization. Figure 3*B* shows the distribution of every term in the optimized cost function (without the weighting coefficients). As a standard of reference, values close to zero show excellent match of the modeled outputs to the experimental data. The values of the heart period (HP) and systolic and diastolic blood pressure (SAP and DAP), and the power of the dominant heart period spectral bands, high frequency (HPH) for tasks B1 and B2, and low frequency (HPL) for 6P task, fitted well for all participants. Values close to one for SAPT and DAPT were not unexpected as the weights for these terms in the cost function were very small. Note that the distributions for the cost function terms across tasks were quite similar, demonstrating the high quality of the optimization method. Figure 3*C* shows the distribution of the total cost function minima for all three tasks. The relatively small values and similar distributions demonstrate good personalization of the model.

#### Comparison of experimental and modeled outputs.

Experimental and modeled blood pressure and heart period for one participant are shown in Figs. 4–6 (also see appendix f). Figure 4 shows modeled and experimental blood pressure changes that were in good agreement quantitatively. However, modeled and experimental data often had a phase difference, as described previously. Phase differences such as shown in appendix g were likely due to the fact that the current model did not include all relevant biological processes (e.g., metabolic, vascular, neurological) that are known to affect the cardiovascular system. This was the rationale for our decision to match derived quantities rather than using the difference between the model and experimental blood pressure curves. The use of derived spectral quantities was supported by the good agreement between the model and the experimental data in terms of heart period (Fig. 5, top panels). In the return maps (Fig. 5, bottom panels), the heart period state (at discrete time *i, X*-axis) is plotted against its state at time *i*+1 (*Y*-axis). In the language of nonlinear dynamics, this represents a two-dimensional embedding. In other words, in a seemingly random time series of the heart period signal, the structure of the loops in the embedding plot shows that there is an underlying partial determinism in the signal.

Finally, Fig. 6 demonstrates physiological validity by demonstrating good agreement between the PSDs of the experimental and modeled heart period smoothed using the Welch method (79) in the high-frequency range for B1 and B2 and in the low-frequency range for 6P. This is conceptually interesting in that, at rest, high-frequency oscillations in heart period are dominated by the respiration input (i.e., reflecting respiratory sinus arrhythmia). However, when breathing is slowed to 6 breaths/min, as in task 6P, respiration exerted a major influence on low-frequency oscillations in heart period, whereas it had very little influence in the high-frequency range. The excellent fit between the PSDs in these task-specific, respiration-driven frequency ranges supports the use of constant-amplitude respiration in the model (i.e., only respiration period varied).

#### Model validation.

To illustrate the validity of the model, the following steps were performed. The model parameterization was performed using only the first half of data. The model output for the second half of each 5-min interval was obtained by supplying the found parameters and the second half of the experimental data. The heart period and its power spectral density, modeled and experimental, are presented in Fig. 7. Note that when the model was parameterized by using only the first half of the data and then predicted for the second half of the data, the results were very similar to those obtained using the whole set of data.

#### Resonance breathing effects on model parameters.

Paired *t*-tests (*P* < 0.05) were conducted to examine whether statistically significant changes occurred within individuals in several important model parameters between the baseline tasks (B1, B2) and resonance breathing task (6P). As expected, we identified significant changes in *T*_{min} between 6P and both baseline tasks. We also found a significant change in *E*_{lv,max} between B1 and 6P, but not between B2 and 6P. In addition, Fig. 8 shows that splanchnic peripheral compliance (*C*_{sp}, top panel), arterial receptor gain on sympathetic control of heart period (*G*_{aTs}, middle panel), and minimum left ventricular elastance (*E*_{lv,min}, bottom panel) significantly increased during the resonance breathing task compared with both baselines. We further examined these three dynamic parameters for each individual. Consistent with the group-level paired *t*-test results, a specific positive reaction to resonance breathing was observed as the modal (most frequent) single response pattern for each parameter. The left panels of Fig. 9 show participants who demonstrated a specific positive response to 6P compared with B1 and B2. There were also several less frequently observed patterns of change between baselines and the resonance breathing task (e.g., high stable values, low stable values, large changes between B1 and B2, but not 6P, and other patterns). The right panel of Fig. 9 shows these less dynamic and other nonmodal patterns of cardiovascular responses to B1, B2, and 6P. Fourteen individuals showed specific positive C_{sp} reactions to 6P, nine showed specific positive *G*_{aTs} reactions to 6P, and 11 showed specific positive *E*_{lv,min} reactions to 6P. Five individuals showed these specific 6P reactions to all three parameters, and six others showed specific 6P reactions to two parameters. An additional four participants demonstrated substantial elevations in *E*_{lv,min} to both B2 and 6P compared with B1. Only one individual showed no positive reaction to 6P across all three parameters.

## DISCUSSION

The present study used a computational physiology approach to better understand, in vivo, how cardiovascular processes change when the baroreflex mechanism is challenged by resonance breathing. We modeled cardiovascular function during a resonance breathing challenge because it is the basis of HRV biofeedback, a promising adjunct intervention for prevalent biobehavioral disorders associated with reduced baroreflex sensitivity and affective dysregulation (33). The present research built on the model originally proposed by Ursino and Magosso (36, 69–71) by using respiration period as the input to the model. In doing so, it extends previous mathematical models of the human cardiovascular system that considered only the mechanical effect of respiration (14), assumed respiration period was constant (36, 53, 69–71), or neglected the effect of respiration altogether (45, 46). In these latter models (45, 46), high-resolution blood pressure was used as the input into the model and thus feedback from heart period to blood pressure could not be modeled. Future computational modeling approaches should also consider adding a spatially distributed component to the model (e.g., 25, 27) in conjunction with a calibration method similar to the one performed here.

The computational results provided support for the elaborated model, which adequately approximated the experimental data. It accurately described human cardiovascular system activities at rest and in response to the resonance breathing challenge. We note that the number of important parameters selected during parameterization was considerably smaller than the total number of parameters in the model. This is consistent with current theories that describe reasons for sensitivity of large systems with a small number of outputs to only a small number of input parameters (e.g., 35). We further suggest that a system can be equally sensitive to a large number of parameters, provided the output is a symmetric function of its inputs. For example, in the linear case where the output equals

where *x*_{i} are input variables, *p*_{i} are parameters, and all *p*_{i} have the same order of magnitude. Thus, the finding that only a small number of parameters was identified as important is noteworthy, given the predominantly linear nature of the model.

The use of statistical inference testing in combination with the computational model was also heuristic in terms of conceptualizing internal cardiovascular processes that are enhanced during resonance breathing compared with baseline (B1 and B2). *T*_{min}, a mechanical parameter related to heart period, was identified as expected based on the high-amplitude cardiovascular oscillations that occurred in response to breathing paced at the heart rate baroreflex resonance frequency. Prior studies have demonstrated that resonance breathing increases the dynamic range of heart and vessel responses to promote adaptability of the cardiovascular system to perturbations (e.g., 30, 33) and, through its actions on baroreflex sensitivity, is likely to promote inhibitory neural processing to support integrated behavioral regulation (for review, see 15).

Three dynamic parameters (splanchnic peripheral compliance, arterial receptor gain on sympathetic control of heart period, and minimum left ventricular elastance) also demonstrated large reactions to resonance breathing, at least in a modal subset of healthy study participants. These parameters capture elements of autonomic nervous system receptor activity (arterial receptor gain), heart muscle condition (minimum left ventricular elastance), and peripheral blood flow (splanchnic peripheral compliance). The results suggest that a machine-learned model of cardiovascular reactions to slow breathing paced at the resonance frequency of the heart rate baroreflex maps well onto what is known from experimental data about cardiorespiratory phenomena. Based on prior literature and our conceptual model of resonance breathing effects, we offer the following explanations for the model results.

First, increases in splanchnic peripheral compliance in response to resonance breathing (seen in 14 participants) are consistent with resonance breathing-induced improvements in blood flow to internal organs. This may occur as the result of greater oscillation amplitudes in vascular tone, which should act to increase compliance. We note that peripheral circulation to the brain's vascular beds is not included in the current model; thus, the impact of resonance breathing on neural circulation is not addressed in this study. Second, increases in arterial receptor gain on sympathetic control of heart period in response to resonance breathing (seen in 9 participants) likely reflect enhanced heart rate baroreflex sensitivity. This increase implies that, during resonance breathing, the baroreceptors respond to smaller changes in blood pressure and thus initiate a more fine-grained (i.e., adaptive) neural response. Increased activation of the baroreceptors is generally observed in response to the high oscillations in blood pressure that are elicited by resonance breathing. This activation also has been shown to promote reductions in stress (49, 62, 63, 80). Finally, elevations in minimum left ventricular elastance in response to resonance breathing (seen in 11 participants) suggest increased heart interaction with systemic vasculature (12). This parameter reflects the passive, unactuated stiffness of the left ventricle. It can change under increased load (51, 72), such as the load created by resonance breathing. Each of these parameters is a potential physiological mechanism that may underlie clinical benefit of resonance breathing and HRV biofeedback interventions, and is consistent with the high-amplitude cardiovascular oscillations elicited by resonance breathing. Moreover, we note that all three of these parameters are sensitive to loading and demonstrate nonlinear responses to high loads (17, 51, 72). The resonance breathing task used in this study causes a substantial cardiorespiratory load and thus may explain one mechanism by which resonance breathing elicits these effects.

It is noteworthy that, on the one hand, resonance breathing appeared to push minimum left ventricular elastance to its maximum level for the majority of participants. On the other hand, there was greater individual variability in the amount of change during resonance breathing observed for splanchnic peripheral compliance and arterial receptor gain on sympathetic control of heart period (see Fig. 8). Variation in these parameters may be informative about individual differences in the clinical response to HRV biofeedback. The distinctive patterns of cardiovascular responding observed across individuals further support a computational physiology approach to parameterize the baroreflex feedback loop model at the individual level.

Although promising, the interpretation of these results is subject to limitations of the study design. The participants in the present study were relatively young and healthy; thus it is unknown whether the important parameters identified here will generalize to older persons and specific clinical populations. Replication and extension to patient populations, such as those with congenital heart disease or other conditions that disrupt venous flow, is an important goal for future study. We also note, consistent with our use of respiration period as the only model input, that this computational model was most accurate in capturing baroreflex activity mediated by respiratory influences in the high-frequency range in the baseline tasks, and in the low-frequency range during resonance breathing (e.g., see Fig. 7). The lack of fit in the very low- and low-frequency ranges of B1 and B2 also likely reflects processes not included in the present model. Further model elaboration is needed to include more of these processes, and to examine the influence of a broader range of internal and external challenges on cardiovascular response. In addition, we note that initial estimates were calculated for each participant based on his/her height (in m) and weight (in kg); however, holding V_{blood} constant across the tasks did not always lead to a good match with the experimental data. Many factors can influence blood volume, including fluid consumption and normal hormonal variations. Thus, when V_{blood} was identified as a sensitive parameter, we chose to include it in the optimization procedure. After optimization, however, the blood volume of a person did not significantly differ across the three tasks. In some sense, this limited variability provided additional validation for our procedure. Moreover, this method for identifying volume of blood is nonintrusive and may be a useful strategy for future studies.

In summary, the computational physiology approach used in this study has potential applications for clinical intervention and, if replicated in clinical samples, may be useful for matching persons to the more effective or targeted interventions. Prior experimental evidence suggests that HRV biofeedback increases baroreflex sensitivity (30). The present study further suggests that increased blood flow to internal organs, sensitivity of the sympathetic component of baroreflex, and ventricular elastance may also be active mechanisms of breathing interventions, such as HRV biofeedback. A computational physiological approach may be useful for distinguishing subgroups of individuals who may benefit from breathing-based interventions compared with those who may not because they deviate from the modal positive cardiovascular responses to resonance breathing in distinctive ways. Thus continued development of this approach appears to have potential to provide information about the relative expected utility of interventions for enhancing specific cardiovascular functions in different persons.

## GRANTS

This study was supported by the National Institute on Alcohol Abuse and Alcoholism through Grants HHSN275201000003C, K01-AA-017473, K02-AA-00325, K24-AA-021778, R01-AA-019511, and AFOSR FA9550-12-1–0230.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s). The AIMdyn's software GoSUM was used for the parameterization and global sensitivity analysis.

## AUTHOR CONTRIBUTIONS

Author contributions: M.F., I.M., J.F.B., V.A.F., A.M., E.V., E.-Y.M., B.V., and M.E.B. conception and design of research; M.F., V.A.F., E.V., and E.-Y.M. analyzed data; M.F., I.M., J.F.B., V.A.F., and M.E.B. interpreted results of experiments; M.F. and E.-Y.M. prepared figures; M.F., I.M., J.F.B., V.A.F., and M.E.B. drafted manuscript; M.F., I.M., J.F.B., E.V., E.-Y.M., B.V., and M.E.B. edited and revised manuscript; M.F., I.M., J.F.B., V.A.F., A.M., E.V., E.-Y.M., B.V., and M.E.B. approved final version of manuscript; B.V. performed experiments.

## Appendix A

#### Model Parameters and States

appendix a provides a detailed list of model parameters (Table A1), model parameters treated as constants (Table A2), and model states (Table A3). The abbreviation, definition, value, and units are given.

## Appendix B

#### Quantitative Model Description

Below we list all equations of the mathematical model. Unless noted otherwise, the equations are taken from Refs. 69 and 71.

The thoracic and abdominal pressures (in mmHg) for each respiratory cycle are calculated by using the following formulas:

where *T*_{resp} is the respiration period (which can be different for each respiratory cycle), *T*_{i} = 0.4 × *T*_{resp} and *T*_{e} = 0.35 × *T*_{resp} denote the duration of inspiration and expiration, respectively. α is a dimensionless variable, ranging between 0 and 1, which represents the fraction of the respiratory cycle. α = 0 corresponds to the beginning of inspiration.

The instantaneous heart period is computed as

where

The heart period changes induced by vagal and sympathetic stimulations satisfy the following equations (1) (2)

where

P_{sa}(*t*) is systemic arterial pressure (one of the states of the model) and the lung volume V_{L}(*t*) is computed from intrathoracic pressure as follows

where the volume is in milliliters and the pressure is in millimeters of Hg.

The sympathetic regulation mechanisms (peripheral resistances, venous unstressed volumes, end-systolic elastances) include a static sigmoidal relationship in series with a dynamic characteristic. The latter incorporates a pure delay and a first-order low-pass filter. The input to the sigmoid is the weighted sum of the information coming from the two groups of receptors. Arterial baroreceptors are sensitive to systemic arterial pressure, whereas lung stretch receptors respond to changes in lung volume. The following equations hold: (3) (4) (5)

(6) (7) (8)

where

where θ can be *E*_{lv}, *E*_{rv}, V_{u,sv}, V_{u,ev}, *R*_{ep}, or *R*_{sp}.

In the following, pressures P_{pa}, P_{pp}, P_{pv}, P_{la}, and P_{ra} were shifted by value P_{thor}(*t*), and pressures P_{sp} and P_{sv} were shifted by value P_{abd}(*t*). This accounts for the fact that extravascular pressure is different from atmospheric pressure in the thoracic and abdominal cavities. Thus vessel transmural pressure inside the thoracic cavity was computed as the difference between the intravascular and intrathoracic pressure, whereas transmural pressure at splanchnic vessels was computed as intravascular pressure minus abdominal pressure,

Balance of forces at pulmonary arteries (*L*_{pa}):
(9)

Balance of forces at systemic arteries (*L*_{sa}):
(10)

Conservation of mass at pulmonary arteries (*C*_{pa}):
(11)

Conservation of mass at pulmonary peripheral circulation (*C*_{pp}):
(12)

Conservation of mass at pulmonary veins (*C*_{pv}):
(13)

Conservation of mass at systemic arteries (*C*_{sa}):
(14)

where

Conservation of mass at peripheral systemic circulation (*C*_{sp} and *C*_{ep}):
(15)

Conservation of mass at extrasplanchnic venous circulation (*C*_{ev}):
(16)

Conservation of mass at right atrium: (17)

Conservation of mass at left atrium: (18)

Conservation of mass at right ventricle: (19)

Conservation of mass at left ventricle: (20)

The total unstressed volume:

Pressure in the splanchnic venous circulation:

Ventricle activation function is a squared half-sine wave

where *T* is the heart period and the duration of systole *T*_{sys} decreases linearly with heart period

The fraction of the cardiac cycle *u*(*t*) is obtained by using an additional state variable ξ(*t*):
(21)

with

where the function “fractional part” frac(ξ) resets the variable *u*(*t*) to zero as soon as it reaches the value +1.

Isometric pressure in left ventricle:

where φ = 1 at maximum contraction, and φ = 0 at complete relaxation.

Isometric pressure in right ventricle:

The viscous resistance of the left ventricle (*R*_{lv}) is assumed to be proportional to the isometric pressure, and so it increases significantly during contraction:

The instantaneous left ventricle pressure can be obtained as the difference between the isometric pressure and viscous losses:

Viscous resistance of right ventricle:

Instantaneous right ventricle pressure:

## Appendix C

#### Calculation of Heart Period HP_{model}, Systolic Pressure SAP_{model}, and Diastolic Pressure DAP_{model}

Modeled heart period HP_{model} was calculated as the length of the time intervals between the minima (see circles in Fig. C1) of the blood pressure Psa_{model} curve multiplied by a constant. If (*t*1,*y*1), (*t*2,*y*2), and (*t*3,*y*3) are three consecutive minima of the Psa_{model} curve, then [*t*1,(*t*2 − *t*1) × const] and [*t*2,(*t*3 − *t*2) × const] are two corresponding points for the HP_{model}. The experimental and modeled HP were then spline interpolated at the same points.

Systolic pressure SAP_{model} (see crosses in Fig. C2) is the maxima of the blood pressure Psa_{model} curve. Diastolic pressure DAP_{model} (see circles in Fig. 11) is the minima of the blood pressure Psa_{model} curve. The experimental and calculated SAP and DAP were then spline interpolated at the same points.

## Appendix D

#### Variance-Based Global Sensitivity

Let *f*(*x*_{1},*x*_{2},…*x*_{n}) be a square integrable function defined in the domain *R*^{n}. The inputs are treated as random variables and their probability density functions represent the associated uncertainty. The impact of the multiple input variables on the output can be independent as well as cooperative, and the analysis of variance (ANOVA) expresses the model output *f*(*x*) as a finite hierarchical cooperative function expansion in terms of its input variables. In order to express the input-output relationship of complex models with a large number of input variables, the mapping between the input variables *x*_{1},…*x*_{n} and the output variables *f*(*x*) = *f*(*x*_{1},…,*x*_{n}), in the domain *R*^{n} can be written in the following form (57):

where *f*_{0} is the constant mean effect (zeroth order), function *f*_{i}(*x*_{i}) is a first-order term describing the effect of variable *x*_{i} acting independently upon the output *f*(*x*), and function *f*_{i},_{j}(*x*_{i},*x*_{j}) is a second-order term describing the cooperative effects of variables *x*_{i} and *x*_{j} upon the output *f*(*x*). The higher-order terms reflect the cooperative effects of increasing numbers of input variables acting together to influence the output *f*(*x*). The last term *f*_{1,2,…n}(*x*_{1},*x*_{2},…*x*_{n}) contains any residual *n*th order cooperative contribution of all input variables. All terms in ANOVA decomposition are orthogonal to each other.

The total variance *D* is computed as follows:

where ρ(*x*) is the probability density of distribution of input variables. Partial variances are defined as follows:

The total partial variances *D*_{i}^{tot} for each parameter *x*_{i}, *i* = 1,2,…,*n*, can be obtained as:

where <*i*> means summation over all that contain index *i*. After the total partial variances are determined, the total sensitivity indexes can be calculated as follows:

By definition, the total partial variance D_{i}^{tot} for each parameter *x*_{i} is

where *f*|*x*_{−i} is a function of *x*_{i} with all other parameters fixed. For independent parameters, the last equation can be rewritten in the following way (58, 59):

## Appendix E

#### Personalized Results of the Parameterization for All Three Experimental Tasks

appendix e provides the personalized parameters for each of the 24 participants for each of the tasks. The parameterization results for tasks B1, B2, and 6P are shown in Tables E1, E2, and E3, respectively.

## Appendix F

#### Parameter Importance in the Model

To determine parameter importance, variance sensitivity analysis was again performed. Parameters were again assumed to be uncertain with uniform distribution. For this analysis, a range of ±10% from the value found during the parameterization was used. Figure F1 shows a representative heat map of variance sensitivity for one participant, with lighter colored cells depicting greater importance. The column numbers match the parameter numbers provided in appendix a, Table A1. V_{Ln} (*column 44*), V_{blood} (*column 37*), P_{san} (*column 43*), *k*_{E,lv} (*column 33*), *T*_{max} (*column 46*), *R*_{ep,max} (*column 62*), *T*_{min} (*column 45*), G_{aTs} (*column 40*), V_{u,ev,max} (*column 70*), V_{u,sv,max} (*column 66*), *C*_{sa} (*column 5*), and *R*_{ep,min} (*column 63*) were identified as important for this individual, and a similar pattern of relative importance was observed for all other participants. This suggests that, once again, only a subset of parameters was important for all terms in the cost function. This provides validation for the selection of the 12 parameters used in model parameterization.

## Appendix G

#### Experimental and Modeled Systolic and Diastolic Blood Pressure During Baseline Task B1, and the Corresponding PSD

Figure G1 shows experimental and modeled systolic and diastolic blood pressure and the corresponding PSDs for task B1 for one of the participants.

- Copyright © 2014 the American Physiological Society

## REFERENCES

- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵