Different Impact of Heart Rate Variability in the Deep Cerebral and Central Hemodynamics at Rest: An in silico Investigation

Background: Heart rate variability (HRV), defined as the variability between consecutive heartbeats, is a surrogate measure of cardiac vagal tone. It is widely accepted that a decreased HRV is associated to several risk factors and cardiovascular diseases. However, a possible association between HRV and altered cerebral hemodynamics is still debated, suffering from HRV short-term measures and the paucity of high-resolution deep cerebral data. We propose a computational approach to evaluate the deep cerebral and central hemodynamics subject to physiological alterations of HRV in an ideal young healthy patient at rest. Methods: The cardiovascular-cerebral model is composed by electrical components able to reproduce the response of the different cardiovascular regions and their features. The model was validated over more than thirty studies and recently exploited to understand the hemodynamic mechanisms between cardiac arrythmia and cognitive deficit. Three configurations (baseline, increased HRV, and decreased HRV) are built based on the standard deviation (SDNN) of RR beats. For each configuration, 5,000 RR beats are simulated to investigate the occurrence of extreme values, alteration of the regular hemodynamics pattern, and variation of mean perfusion/pressure levels. Results: In the cerebral circulation, our results show that HRV has overall a stronger impact on pressure than flow rate mean values but similarly alters pressure and flow rate in terms of extreme events. By comparing reduced and increased HRV, this latter induces a higher probability of altered mean and extreme values, and is therefore more detrimental at distal cerebral level. On the contrary, at central level a decreased HRV induces a higher cardiac effort without improving the mechano-contractile performance, thus overall reducing the heart efficiency. Conclusions: Present results suggest that: (i) the increase of HRV per se does not seem to be sufficient to trigger a better cerebral hemodynamic response; (ii) by accounting for both central and cerebral circulations, the optimal HRV configuration is found at baseline. Given the relation inversely linking HRV and HR, the presence of this optimal condition can contribute to explain why the mean HR of the general population settles around the baseline value (70 bpm).

Background: Heart rate variability (HRV), defined as the variability between consecutive heartbeats, is a surrogate measure of cardiac vagal tone. It is widely accepted that a decreased HRV is associated to several risk factors and cardiovascular diseases. However, a possible association between HRV and altered cerebral hemodynamics is still debated, suffering from HRV short-term measures and the paucity of high-resolution deep cerebral data. We propose a computational approach to evaluate the deep cerebral and central hemodynamics subject to physiological alterations of HRV in an ideal young healthy patient at rest.

Methods:
The cardiovascular-cerebral model is composed by electrical components able to reproduce the response of the different cardiovascular regions and their features. The model was validated over more than thirty studies and recently exploited to understand the hemodynamic mechanisms between cardiac arrythmia and cognitive deficit. Three configurations (baseline, increased HRV, and decreased HRV) are built based on the standard deviation (SDNN) of RR beats. For each configuration, 5,000 RR beats are simulated to investigate the occurrence of extreme values, alteration of the regular hemodynamics pattern, and variation of mean perfusion/pressure levels.
Results: In the cerebral circulation, our results show that HRV has overall a stronger impact on pressure than flow rate mean values but similarly alters pressure and flow rate in terms of extreme events. By comparing reduced and increased HRV, this latter induces a higher probability of altered mean and extreme values, and is therefore more detrimental at distal cerebral level. On the contrary, at central level a decreased HRV induces a higher cardiac effort without improving the mechano-contractile performance, thus overall reducing the heart efficiency.
Conclusions: Present results suggest that: (i) the increase of HRV per se does not seem to be sufficient to trigger a better cerebral hemodynamic response; (ii) by accounting for both central and cerebral circulations, the optimal HRV configuration is found at baseline.

INTRODUCTION
Heart rate variability (HRV), defined as the variability between successive RR heartbeats, has grown as hot topic, with a simple "Heart rate variability" topic search currently listing more than 44,000 results on Web of Science. The great interest elicited within the scientific community ranges from psychophysiology (Laborde et al., 2017) to exercise (Michael et al., 2017), cardiovascular risk factors (Tsuji et al., 1996), and chronic fatigue syndrome (Meeus et al., 2013). This increasing and wide interest is also due to the fact that HRV measurements are easy to perform, quite reliable, and non-invasive. In particular, HRV can be considered a surrogate measure of the cardiac vagal tone (Laborde et al., 2017), although there are controversial aspects in using HRV to determine cardiac vagal tone independently from other factors, such as age, physical condition, and the presence of cardiac pathologies (Boyett et al., 2019). In this paper, we refer to the parasympathetic activity within cardiac regulation as cardiac vagal tone. HRV parameters in time and frequency domain able to reflect cardiac vagal tone (Malik et al., 1996;Berntson et al., 1997;Laborde et al., 2017) show that an increased cardiac vagal tone is associated to a higher HRV, while a reduction of the cardiac vagal tone response is linked to a HRV reduction (Hayano et al., 1991;Singh et al., 2018). Moreover, it has recently being recognized that HRV (as measured by means of different metrics, such as SDNN, RMSSD, pvRSA, and HF), inversely correlates with the heart rate (HR), and this relation should be taken into account when dealing with increased/reduced HRV (Monfredi et al., 2014;Shaffer and Ginsberg, 2017;Boyett et al., 2019;de Geus et al., 2019). The inverse relation between HR and HRV can be explained considering that a faster HR reduces the interval between successive beats and the opportunity for the RR beating to vary, lowering the HRV. On the contrary, a slower HR increases the cardiac interval and enhances the chance for RR to vary, raising HRV (Shaffer and Ginsberg, 2017).
From an overall cardiovascular point of view, a HRV increase in the physiological range seems beneficial at rest, while a reduced HRV is symptomatic of a stress action, can be predictor of cardiovascular risk factors and is correlated to higher morbidity and mortality (Tsuji et al., 1996;Dekker et al., 2000;Meeus et al., 2013). Thus, the clinical interest has been mainly focused on the prognostic significance of HRV related to the risk factors of cardiovascular pathologies. There is substantial evidence that a decreased HRV enhances several risk factors and is associated to cardiovascular diseases, such as left ventricle hypertrophy (Acharya et al., 2006), sudden cardiac death (Sessa et al., 2018), and myocardial infarction (Buccelletti et al., 2009). On the contrary, lower risk factor profiles are associated with increased HRV (Thayer et al., 2009b). The HRV scenario is opposite under physical effort: for increasing moderate-to-vigorous exercise intensity HRV decreases (Tulppo et al., 1998;Michael et al., 2017). During exercise a HRV decrease is beneficial and forced by an increased HR, however there is an immediate post-exercise recovery of HRV (Michael et al., 2017) and exercise training can increase resting HRV (Billman and Kukielka, 2006).
Although the overall HRV impact on the central hemodynamics is quite clear, especially related to cardiac disease conditions, its influence on the cerebral hemodynamics is barely known so far. In fact, definitive association between HRV and cerebral circulation is still missing and mainly involve pathologic conditions, such as chronic fatigue syndrome (Meeus et al., 2013;Boissoneault et al., 2019) and stroke (Fyfe-Johnson et al., 2016). Attention has recently grown about the possible association between increased HRV and improved cognitive performance. Results are conflicting, as some studies have demonstrated the association (Kim et al., 2006;Shah et al., 2011;Schaich et al., 2020), while some others observed that a reduced HRV does not contribute to cognitive impairment or even dementia (Allan et al., 2005;Britton et al., 2008;Mahinrad et al., 2016;Hazzouri et al., 2017). To understand this controversy it should be kept in mind that, among many influencing factors, such as the sample and the objectives of the study, most of the literature is based on short term measures, which have lower prognostic value than 24-h HRV (Malik et al., 1996;Nunan et al., 2010;Shaffer and Ginsberg, 2017).
Though a possible association is widely debated, the underlying hemodynamic mechanisms are so far only hypothesized and mostly unclear (Schaich et al., 2020). It is presently unknown whether HRV is able to influence the normal pressure and flow rate pattern at distal-capillary level, to alter the mean perfusion and pressure levels, and to enhance the occurrence of extreme values. Detailed information are lacking since current non-invasive techniques, such as transcranial doppler ultrasonography, are not able to offer high-resolution data in terms of pressure and flow rate beyond the circle of Willis. Therefore, in the deep cerebral circulation hemodynamic measurements are unreliable or almost absent. Thus, it can be of great significance understanding and predicting from a computational point of view how deep cerebral hemodynamics is affected by HRV at rest. Despite the intrinsic limits, computational hemodynamics and cardiovascular modeling are becoming increasingly important to isolate specific mechanisms and investigate processes where clinical data are not yet feasible, accurate, or easily measurable.
We propose an in-silico study where, through a validated combined cardiovascular-cerebral model, we evaluated the deep cerebral hemodynamics subject to physiological alterations of HRV in an ideal young healthy patient at rest, thereby comparing the cerebral to central hemodynamics responses to HRV changes. With deep cerebral circulation we refer to the microcirculation beyond the circle of Willis (i.e., distal and capillary-venous circulation), while the set of parameters defined through cardiac variables and central arterial pressure characterizes the central hemodynamics. The complete lumped-parameter model is composed by a suitable combination of electrical counterparts (resistances, compliances/elastances, inertances), accounting for the arterial and venous circuits of both systemic and pulmonary circulations, an active representation of the four cardiac chambers, an accurate valve motion description, and a short-term baroreceptor mechanism. The cerebral circulation is divided into three main regions: large arteries, distal arterial circulation, and capillary/venous circulation. Cerebrovascular control mechanisms of autoregulation and CO 2 reactivity are taken into account. The present computational approach has been validated and recently exploited in a wide area of applications, such as the hemodynamic response to exercise in atrial fibrillation (Anselmino et al., 2017), the linking mechanisms between cardiac arrythmia and cognitive impairment (Anselmino et al., 2016;Scarsoglio et al., 2017a,b;Saglietto et al., 2019), and the cardiovascular deconditioning emerging in altered gravity conditions (Gallo et al., 2020).
We considered three HRV configurations (baseline, increased HRV, decreased HRV), each of them analyzed over 5,000 RR beats, so that our results are not affected by transient behavior and are statistically significant and stable. HRV variations were assessed through SDNN, which is the standard deviation of all RR beats and allowed us to define in a straightforward and univocal way the RR stochastic beating extraction. The model is able to provide the whole pressure and flow rate time-series from proximal to distal circulation as well as beat-to-beat values. The focus was on the possible occurrence of extreme values, such as hypertensive events, alteration of the regular hemodynamics pattern, and variation of mean perfusion/pressure levels. All these aspects are to date mostly unexplored, but extremely useful to understand the role of HRV on the cerebral circulation and how, in turn, the hemodynamic alteration can impact the cognitive sphere. The present work can offer precious insights by: (i) isolating the net basic mechanisms induced by HRV changes into the cerebral circulation in healthy resting conditions, and (ii) quantifying how differently HRV acts on the cerebral hemodynamics with respect to the central hemodynamics, thus fostering necessary future clinical measurements within this topic.

HRV Configurations
We considered three HRV configurations: baseline condition (reported in the following sections with blue color), increased HRV (in red color) and decreased HRV (in black color). To assess HRV variations we adopted the simplest HRV metric, i.e., SDNN-which is the standard deviation of all RR beats. SDNN belongs to the time-domain analyses, which are computationally simpler and easier to apply than frequency-domain analyses (Michael et al., 2017). Since we focused on resting configurations lasting hours, we chose SDNN because it is the gold standard for long-term measurements (Malik et al., 1996;Shaffer and Ginsberg, 2017). Moreover, the SDNN adoption allows us to define in a straightforward and univocal way the RR beating extraction. In terms of SDNN, HRV was found to be dependent on the circadian day/night cycle (Shaffer and Ginsberg, 2017), e.g., 93.06 (ms) day, 121.31 (ms) night, 101.71 (ms) 24-h (Talib et al., 2005)-and to decrease with age in adulthood (van den Berg et al., 2018). In general, SDNN recording was lower in resting supine condition [65 (ms) (Lehavi et al., 2019); 49 (ms) (Nunan et al., 2010)] than in non-resting or upright position [e.g., 93.06 (ms) (Talib et al., 2005); 127 (ms) (Genovesi et al., 2007)]. It was recently observed that HRV inversely correlates with HR through an exponential function (Monfredi et al., 2014;Kazmi et al., 2016;Shaffer and Ginsberg, 2017;van den Berg et al., 2018;de Geus et al., 2019), and this relation should be taken into account when dealing with increased/reduced HRV.
To consider SDNN for a healthy young person in resting supine condition during day (awake) and account for the interplay between SDNN and HR, we adopted the exponential relation obtained by de Geus et al. (2019) in leisure baseline condition: The baseline configuration was set at HR = 70 bpm and SDNN = 71.14 (ms) was obtained from the above relation. Recalling that the beating period is RR = 60/HR (s), we define the coefficient of variation as c v = SDNN/RR. The resulting c v is Frontiers in Neuroscience | www.frontiersin.org 0.08 which well represents normal sinus rhythm daily values in resting supine condition [e.g., c v = 0.05 (Pikkujämsä et al., 2001); c v = 0.07 (Lehavi et al., 2019); c v = 0.06 (Nussinovitch et al., 2011)]. Moreover, c v = 0.08 as baseline value is close to the one adopted by our group (c v = 0.07) to simulate resting supine sinus rhythm (Scarsoglio et al., 2014). Increased and decreased HRV configurations were achieved by changing SDNN by +20 and −20%, respectively. This threshold represents a significant variation within a range of physiological values. The corresponding HRs were individuated again using Equation (1) as proposed by de Geus et al. (2019). Figure 1A reports the SDNN(HR) curve in resting condition introduced by de Geus et al. (2019) with the values here chosen and below a summarizing table of the three cases.
For each configuration, we evaluated 5,000 RR beats to have statistically significant and stationary results. RR beatings were extracted from an in silico pink-correlated Gaussian distribution, which well reproduces the typical beating features of normal sinus rhythm recorded in vivo (Hayano et al., 1997;Pikkujämsä et al., 2001;Hennig et al., 2006;Scarsoglio et al., 2014), having mean and standard deviation values identified as above. In the range SDNN ∈ [30, 100] ms, we also evaluated the RMSSD (i.e., the root mean square of successive RR interval differences) for the RR series extracted with the adopted time-correlation features and composed by 10 8 beats. A high linear correlation value (R 2 = 0.99) was found between RMSSD and SDNN, with fitting law RMSSD = 0.65 SDNN + 0.54, as reported in Figure 1B. Although RMSSD and SDNN don't have the same physiological origin and only RMSSD can be assumed to fully reflect cardiac vagal activity (Malik et al., 1996;Berntson et al., 1997;Laborde et al., 2017), SDNN is a reliable proxy of RMSSD for the considered RR time series. Bottom panels of Figure 1 display the RR series (panel C) and the corresponding probability density functions (PDFs) for the three HRV cases here investigated (panel D).

Cardiovascular and Cerebral Circulation Modeling
Following RR extraction, the cardiovascular-cerebral model was run to obtain the hemodynamic cerebral signals. The complete lumped model is composed by electrical components able to reproduce the response of the different (cardiac and vascular) regions and their features, in terms of resistance R (diffusive effects), inertance L (inertial effects), and compliance C or elastance E (elasticity/contractility effects). The cardiovascular dynamics includes the arterial and venous circuits of both systemic and pulmonary circulations, an active representation of the four cardiac chambers, and an accurate valve motion FIGURE 3 | (A,B) Examples of percentile evaluation for P dm,left (p(P dm,left ) is the probability density function). Ten and 90% blue dashed lines individuate the 10th and 90th percentiles in the baseline configuration (A,B, blue areas), while they correspond to the 2nd and 85th percentiles in the decreased HRV configuration (A, black areas), and to the 22nd and 94th percentiles in the increased HRV configuration (B, red areas). (C,D) P dm,left examples of hypertensive event lasting two beats in the decreased HRV configuration (average pressure per beat, P, is represented by gray horizontal lines and hypertensive events are marked with bold gray lines) and hypotensive event in the increased HRV configuration (P is represented by red horizontal lines and hypothensive events are marked with bold red lines). The 10th and 90th percentile baseline thresholds are displayed through the dashed blue horizontal lines, while full circles indicate the heartbeat extremes. In (D) examples of P max (△), P min (2), and P pv (↔) values per beat are evidenced for P dm,left .
Frontiers in Neuroscience | www.frontiersin.org description. A short-term baroreceptor mechanism is also modeled, accounting for the inotropic effect of both ventricles, as well as the control of the systemic vasculature (peripheral arterial resistances, unstressed volume of the venous system, and venous compliance). The chronotropic effects due to the heart rate regulation are intrinsically taken into account by the RR extraction. The cerebral circulation is divided into three principal regions: large arteries, distal arterial circulation, and capillary/venous circulation. Cerebrovascular control mechanisms of autoregulation and CO 2 reactivity are taken into account. The model is expressed in terms of pressures, P, flow rates, Q, volumes, V, and valve opening angles, θ .
The cardiovascular model has been proposed and validated over more than thirty studies to check the consistency of the hemodynamic response during AF: extensive details of this evaluation are reported in Scarsoglio et al. (2014). Then, the model has been exploited to study the impact of atrial fibrillation on the cardiovascular system (Anselmino et al., 2015(Anselmino et al., , 2017Scarsoglio et al., 2016a,b). The complete cardiovascular-cerebral model has been used to understand the hemodynamic mechanisms between atrial fibrillation and cognitive deficit (Anselmino et al., 2016;Scarsoglio et al., 2017a,b;Saglietto et al., 2019). A full description of the governing equations and model parameters is given in the Supplementary Material.
We here focused-in terms of pressure P and flow rate Q-on the left internal carotid artery-middle cerebral artery (ICA-MCA) pathway, which already turned out to be representative of the hemodynamics from proximal to distal cerebral districts (Anselmino et al., 2016;Scarsoglio et al., 2017a,b;Saglietto et al., 2019). The left ICA-MCA path starts at the internal carotid level (P a : systemic arterial pressure; Q ICA,left : left internal carotid flow rate), goes through the middle cerebral artery (P MCA,left : left middle cerebral artery pressure; Q MCA,left : left middle cerebral artery flow rate), includes the middle distal regions (P dm,left : left middle distal pressure; Q dm,left : left middle distal flow rate), and ends with capillary-venous districts (P c : cerebral capillary pressure; Q pv : proximal venous flow rate). A schematic representation of the systemic-pulmonary circulation and ICA-MCA pathway is reported in Figure 2, together with examples of left ventricle PV loops, aortic flow rate (Q ao ), systemic arterial pressure (P a ), and cerebral capillary pressure (P c ) time-series for different HRV configurations.

Variable Definition and Data Analysis
We recall the definition of mechano-energetic and oxygen consumption indexes, which will be used to describe the central hemodynamics. End-systolic left ventricular volume, V lves [ml], is the left ventricle volume at the closure of the aortic valve, while end-diastolic left ventricular volume, V lved [ml], corresponds to the closure of the mitral valve. Stroke volume is defined as SV = (V lved − V lves ) [ml], ejection fraction is EF = V/V lved · 100 [%]. Cardiac output is CO = SV · HR [l/min], stroke work per minute, SW/min [J/min], is measured as the area within the left ventricle pressure-volume loop per beat. Oxygen consumption is evaluated through three estimates (Westerhof et al., 2010): (i) the rate pressure product, RPP = P a,syst · HR [mmHg/min], where P a,syst is the aortic systolic pressure; (ii) the tension time index per minute, TTI/min = P lv · RR · HR [mmHg s/min], where the symbol f indicates the mean value of the generic hemodynamic variable f averaged over a RR beat (in this case f = P lv ); and (iii) the pressure volume area per minute, PVA/min = (PE+SW)·HR, where PE = [P lves · (V lves − V lv,un )/2 − P lved · (V lved − V lv,un )/4] is the elastic potential energy (V lv,un = 5 ml is the unstressed left ventricle volume), SW is the stroke work, P lves and P lved are the end-systolic and end-diastolic left ventricle pressures. The left ventricular efficiency LVE is defined by the ratio SW/PVA.
The cerebral hemodynamics is assessed by means of two main approaches: (i) analysis of the continuous time-series, where the signal is continuous and defined by all the temporal instants of the whole time-series, and (ii) beat-to-beat analysis, where the signal was discretized and one-per-beat data were obtained. In so doing, beat-to-beat signals are composed by 5,000 values, corresponding to the 5,000 RR beats simulated.
For the analysis of the continuous time-series, extremely high or low cerebral hemodynamic values were evaluated through the percentile analysis of the hemodynamic signals. We adapted the definition used for studying atrial fibrillation (Anselmino et al., 2016;Saglietto et al., 2019), to assess the possibility of extreme events related to HRV changes. For the increased and decreased HRV cases, we estimated to which percentile the baseline (10th and 90th) reference thresholds  are able to make extreme values more frequently reached with respect to the baseline case (e.g., in Figure 3B the 10th baseline percentile corresponds to the 22nd percentile for increased HRV). For the beat-to-beat analysis, the i-th element of the discretized series may contain the average value (Q and P), the maximum (Q max and P max ) and minimum (Q min and P min ) values, as well as the pulse value (Q pv and P pv , defined as the difference between maximum and minimum values) of the related hemodynamic variables computed over the i-th beat. In Figure 3D examples of P, P max , P min , and P pv are represented for P dm,left .
The average values per beat, P and Q, were still exploited to enrich the beat-to-beat analysis, by accounting for the persistence of extreme values over the whole beat (Anselmino et al., 2016;Saglietto et al., 2019), and not only the occurrence of instantaneous peak values. We introduce hypoperfusions (or hypotensive events), which take place when the mean flow rate per beat Q (or mean pressure per beat P) is below the 10th percentile referred to the whole baseline signal. On the contrary, hyperperfusions (or hypertensive events) emerge when Q (or P) is above the 90th percentile of the baseline condition. Note that, by definition, extreme events cannot emerge in the baseline configuration. Bottom panels of Figures 3C,D show examples of hypertensive and hypotensive events for P dm,left . Table 1 reports the main statistics (mean µ, standard deviation σ , and coefficient of variation c v ) of the continuous time-series of the hemodynamic variables along the ICA-MCA pathway. Figure 4 shows probability density functions (PDFs) of the whole hemodynamic time-series for the three configurations. Considering the continuous hemodynamic signals, HRV influenced much more pressure than flow rate series. In fact, mean values were overall maintained (see Table 1) and PDFs revealed self-similar features (see Figure 4) for flow rates Q, while HRV significantly affected mean values and PDFs of pressures P (only P c at the end of the pathway regained a self-similar shape for the three configurations). At a given district, HRV impacted σ relative variations with respect to baseline similarly for both pressure and flow rate, while c v relative variations induced by HRV with respect to baseline were more evident toward the distal-capillary circulation.

Continuous P and Q Time-Series Analysis
It is then useful to evaluate how extreme values of the continuous time-series P(t) and Q(t) distribute as HRV changes, by considering the percentile analysis (as represented in Figures 3A,B). By assessing to which percentile the baseline reference thresholds correspond in the configurations with altered HRV, we can quantify whether HRV is able to modify the probability of reaching extremely high/low values (see Figure 5).  Distal-capillary regions were the most affected by HRV: the 10th percentile for baseline P dm,left corresponds to over the 20th percentile for increased HRV, while the 90th percentile for baseline P c corresponds to about the 82th percentile for the reduced HRV configuration (see top panels of Figure 5). It is worth noting that a similar scenario was found for flow rates as well. The increased HRV enhanced extreme values of the continuous time-series toward the venous circulation: the 10th percentile corresponds to the 17th, while the 90th to the 82th (see bottom panels of Figure 5).

Beat-to-Beat Analysis
The scenario of higher impact of HRV on cerebral pressures than flow rates was also confirmed in the beat-to-beat analysis. Table 2 displays mean and standard deviation values for the minimum (P min and Q min ), average (P and Q), maximum (P max and Q max ), and pulse (P pv and Q pv ) values per beat of pressures P and flow rates Q as computed over 5,000 RR beats. PDFs of P and Q along the ICA-MCA pathway are depicted in Figure 6. The altered HRV did not substantially modify Q values, with relative variations of increased and decreased HRV <0.1% with respect to baseline, see the Q column in Table 2 and the corresponding PDFs (see Figure 6, right panels). Relative variations of minimum (Q min ) and maximum (Q max ) values per beat were within 8%. P values as well as the corresponding PDFs were instead much more affected by HRV: only at the capillary level, similar P c values were recovered for the three HRV configurations (with relative variations <0.1%, see bottom left panel of Figure 6 and the mean value column of Table 2), while in all the other regions relative P variations were within 3%. Relative variations of maximum (P max ) and minimum (P min ) values per beat did not exceed 5% with respect to the baseline configuration (see Table 2). The different impact of HRV change on pressure and flow rate can be explained noting that cerebral control mechanisms act on the continuous flow rate values at the distal level only, by modifying pial arterial-arteriolar compliances and resistances (Ursino and Giannessi, 2010). Thanks to the autoregulation effects, mean values of the continuous flow rate time-series ( Table 1) and beat-to-beat Q recordings ( Table 2) were preserved unvaried in the three configurations, while the same was not true for cerebral pressure. Despite cerebral autoregulation mechanisms are concentrated in the distal district solely, their effects involve both upstream and downstream regions. In fact, the mean flow rate is guaranteed as constant for the three configurations throughout the ICA-MCA pathway (i.e., from Q ICA,left to Q pv , see Q mean values in Table 1 and Q in Table 2).
Based on what emerged so far, we expect a higher occurrence of extreme events in the beat-averaged pressure, P, than beataveraged flow rate, Q. Table 3 reports the total number of rare one-beats hypotensive, hypertensive, hypoperfusion, and hyperperfusion events. We recall that these events, by definition, cannot occur in the baseline configuration, which is taken to set the reference thresholds (10th and 90th percentiles). For the other configurations, one of these events happens at a specific district if the mean value per beat for P or Q reaches extremely high (above the 90th percentile of the baseline configuration) or low (below the 10th percentile of the baseline configuration) values. In particular, out of 5,000 RR beats we counted 42 hypotensive events at distal level for the increased HRV configuration, and 69 hypertensive events for reduced HRV. No hypoperfusions were found, while other more occasional events, such as nine distal hyperperfusion events, six capillary hypertensive events, and four distal hypertensive events, were induced by the higher HRV configuration. The reason of the low occurrence of extreme perfusive events can be understood by observing that the widening of the PDFs for the continuous flow rate time-series caused by increased HRV (see the red curves in right panels of Figure 4) is related to instantaneous peak values rather than enduring events in time. Thus, increased HRV was not able to trigger significant extreme perfusive events over a whole beat. In fact, PDFs of Q values were quite coincident for the three 3 | Total number of one-beat extreme events (out of 5,000 RR beats) along the ICA-MCA pathway for the two decreased and increased HRV configurations.

Configuration
Pressure

Hypotensive events Hypoperfusion events
Decreased configurations (right panels of Figure 6), leading to very sporadic hyperperfusion events and no hypoperfusion events at all. It should be also recalled that HRV variations were taken in a physiological range (±20% of the baseline value), therefore it is reasonable not expecting a large number of extreme events based on averaged-beat variables.

DISCUSSION
The present study aims at computationally investigating the role of HRV on the deep cerebral and central hemodynamics. SDNN has been adopted as a proxy measure for HRV and, although it does not have the same physiological origin reflecting cardiac vagal activity as RMSSD, it is highly correlated with the latter and more suited to the proposed computational approach. HRV seems to have overall a stronger impact on P than Q values, due to autoregulation effects acting to maintain an adequate average cerebral perfusion. However, in terms of extreme (minimum and maximum) values and percentile distributions, HRV similarly altered both pressure and flow rate. By comparing reduced and increased HRV, this latter induced a higher probability of altered mean values and extreme values (the only exception was the P dm,left case with reduced HRV and hypertensive events), and is therefore more detrimental at distal cerebral level.
As a further confirmation of the worsening linked to a higher HRV, it is observed that starting with an increased HRV, the variability of hemodynamic variables increased from proximal to distal regions along the left ICA-MCA pathway. On the contrary, starting with a reduced HRV, the hemodynamic variability decreased toward the deep cerebral circulation. This behavior is evident for the continuous signals in terms of c v ratio (as reported in Figure 7), as well as at beat-tobeat level, through relative variations of pulse values, P pv and Q pv (see Figure 8). It should be noted that the relative variations of P pv and Q pv showed this increasing/decreasing trend, although P min , P max , Q min , and Q max did not singularly present, as mentioned before, a similar behavior: relative variations of P min and P max (Q min and Q max ) with respect to baseline did not exceed 5% (8%) along the ICA-MCA pathway (please refer to Table 2, columns with minimum and maximum values). Moreover, as displayed by Figures 7, 8, the behavior was analogously found for both pressure and flow rate.
Even though all results fell within a range of physiological values, increased HRV seems more deleterious than reduced HRV along the proximal-to-distal cerebral pathway. An unsafe increased pulsatility of the hemodynamic signals was found especially at capillary-venous levels, and can promote extremely high/low pressure and flow rate values. This result can be explained recalling what observed for the cerebral hemodynamics during atrial fibrillation (Scarsoglio et al., 2017b). The intrinsic structural latency revealed by the cerebral microcirculation and combined with a sequence of inseries irregular RR beats produced an alteration of the cerebral circulation. Indeed, due to the delay related to mechanical and structural properties, the inertia of the system increased toward the microvasculature.
Atrial fibrillation represents a succession of disturbances in terms of RR beating, which are then transmitted from the carotid level to the deep cerebral circulation. Each of these perturbations singularly produced an alteration of the cerebral circulation. The continuous sequence of transient perturbations at the carotid entrance represented by the fibrillated beating did not allow the system to recover the physiological state before another disturbance arrives. When this disturbance chain spread throughout the cerebral vessel network, the distal and capillary regions remained altered for longer. The behavior is similar to a system of springs in series and parallel, which is externally forced at one end: the stiffness of each spring combines with the others and the oscillation survives even if the external perturbation is ceased. It is important to note that this mechanism occurred regardless of whether the beat in atrial fibrillation is accelerated or not. In fact, the comparison between NSR and atrial fibrillation (Scarsoglio et al., 2017b) was proposed at the same mean heart rate, highlighting the role of RR variability in promoting the observed hemodynamic alteration.
With due proportions, the basic mechanism detected in atrial fibrillation seems to be similar in the case of increased HRV. Clearly, fibrillated outcomes were much more exacerbated and pathological since variability was much higher, and the RR beating was uncorrelated too. However, also here the component of increased variability seems to plausibly combine with the intrinsic latency of the system, causing the observed altered scenario.
The hemodynamic framework is inverted if we focus at central level, in terms of cardiac contractility, efficiency end energetic indexes (see Table 4). The decreased HRV configuration showed an increase of oxygen consumption indexes (+16% RPP, +6% TTI/min, +11% PVA/min) and cardiac work (+9% EW/min), in front of a reduction of the ejection fraction (−5%) and a limited increase of cardiac output (+6%). The reduction of the left ventricle efficiency LVE, though slight (−2%), confirms that a higher cardiac effort did not translate into a better mechanical and contractile cardiac performance. Moreover, variations of mean central flow rate (in terms of CO) did not entail average changes of the flow rate in the deep cerebral circulation (see cerebral mean flow rates in Tables 1, 2), due to the autoregulation mechanisms at the distal cerebral level. In other words, the greater central energy expenditure did not produce any gain in terms of mean flow rate in the cerebral microcirculation.
To extend the discussion, we considered a configuration of heart failure as representative of a very different condition from the healthy state and often present as a co-morbidity in the elderly population. With respect to the baseline heart failure condition, we analyzed two HRV configurations also here, increased (+20% SDNN) and decreased (−20%  Figure 9. Although the values achieved are different from the case healthy (see the blue thick and thin curves), the HRV variations with respect to the baseline are quite close to the healthy case (see the red/black thick and thin curves), especially in the distal regions. Table 5 shows the values of the central hemodynamic parameters in the heart failure case. Compared to the healthy case, the values achieved are very different and show a rather compromised hemodynamic framework (e.g., at baseline CO = 3.55 l/min, EF = 16.41%, LVE = 0.31). However, the HRV variations compared to baseline are in agreement with the healthy case (decreased HRV: +12% RPP, +3% TTI/min, +10% PVA/min, +4% EW/min, −8% EF, +3% CO, −5% LVE).
For both healthy and heart failure cases, at central level there is a qualitative agreement with what has been widely observed in literature, that is an increased HRV is beneficial for the cardiac performance and efficiency. At cerebral level, where instead results are still debated, an increase of HRV only-regardless any other variation, compensatory mechanism or pathological condition-seems to be associated to a worse hemodynamic response. A plausible interpretation of these findings with respect to the possible association between HRV and cognitive performance observed in some literature studies is that we only considered the net contribution of HRV alterations, without any complementary and further variation of the sympathetic system. In fact, changes of the sympathetic activity are hardly quantifiable and out of the objectives of the present study. According to the neurovisceral integration model vagally-mediated HRV metrics, such as RMSSD and HF, are expected to be linked to cognitive performance, but not SDNN (Thayer et al., 2009a). Based on our outcomes the increase of HRV, measured by means of SDNN, does not seem to be sufficient per se to trigger a better cerebral hemodynamic response.

Limitations
The present computational approach has some limiting aspects related to the modeling hypotheses. We considered the same young healthy resting subject as forced through three different HRV configurations. In particular, the three configurations only differ by the entrance inputs, the beating sequence RR, while the remaining hemodynamic framework is set as in baseline condition, regardless of other compensatory mechanisms that may occur or altered sympathetic activity. Cerebral autoregulation and autonomic regulation are not fully coupled, but this interplay is present in one-way direction only: through the baroreceptor mechanisms acting on the systemic arterial pressure, autonomic regulation influences the cerebral autoregulation, but there are no feedbacks from cerebral blood flow to autonomic control. Moreover, coronary circulation was not included into the model. In the end, we evaluated a monitoring temporal window of about 1 h (5,000 RR beats), thus no long-term remodeling effects related to increased or decreased HRV were taken into account.

CONCLUSIONS
In conclusion, by recalling the relation inversely linking SDNN and HR (Monfredi et al., 2014;Shaffer and Ginsberg, 2017;de Geus et al., 2019), at deep cerebral level a higher SDNN is worse than a higher HR. On the contrary, at central level a higher HR is more negatively impacting than a higher SDNN. From an overall point of view which contemporarily accounts for both central and deep cerebral circulations, present results suggested that the optimal HRV configuration is found at baseline. In fact, an increased HRV (lower HR) is detrimental for the cerebral circulation causing possible hypertensive events and extreme pressure and flow rate values, while a decreased HRV (higher HR) induces a higher cardiac effort without improving the mechano-contractile performance, thus overall reducing the heart efficiency. The presence of this optimal condition is the main finding and can contribute to explain why the mean HR of the general population settles around 70 bpm: the baseline value is a good compromise able to guarantee an adequate hemodynamic response for both central and deep cerebral regions, preserving from extreme cerebral hemodynamic events and at the same time maintaining satisfactory cardiac efficiency.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.