Understanding the interplay between baroreflex gain, low frequency oscillations, and pulsatility in the neural baroreflex

The neural baroreflex, which regulates mean arterial pressure (MAP) via the action of the brain, consists of baroreceptors which measure MAP, and actuators that can produce a change in MAP, such as the heart and parts of the peripheral resistance containing innervated smooth muscle. The brain is the controlling unit, maintaining an appropriate MAP in spite of various disturbances. Under certain circumstances, including haemorrhage and other states of distress, the gain of the neural baroreflex can change, causing low frequency (LF) oscillations (sometimes termed Mayer waves) in blood pressure (BP). Though their purpose is unclear, the origins of these LF oscillations has previously been explained via a nonlinear feedback model, though focusing on the peripheral resistance as an MAP actuator only. The present paper now includes analytical and simulation results explaining the LF oscillation phenomenon for the full neural baroreflex, containing both peripheral resistance (PR) and cardiac branches. However, the main contribution of the paper is to examine the effect of blood pulsatility, or a lack of pulsatility, on the neural baroreflex, and how it's effect can manifest in the presence of LF oscillations. This may have importance in cases where pulsatility is reduced (for example where left-ventricular assist devices are present), or completely absent (for example in turbine-based artificial hearts). © 2020 The Author(s). Published by Elsevier B.V. on behalf of Nalecz Institute of Biocybernetics and Biomedical Engineering of the Polish Academy of Sciences. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Available online at www.sciencedirect.com


Introduction
The human body displays a variety of oscillatory phenomena, many related to the necessary natural means by which flow (of both air and fluids) can be generated. In particular, the rhythms of respiration and blood flow are well known and understood, though some rhythms are less so, for example the low frequency (LF) oscillations in blood pressure (BP), often termed Mayer waves [1], occurring at roughly 0.1 Hz in humans. It should be noted that the frequency of this LF oscillation is species dependent, occurring at roughly 0.3 Hz in rabbits and 0.4 Hz in rats [2].
Since we use the same inlet and outlet pathway for breath, it makes sense that a reciprocating (oscillating) flow is employed. However, in the case of the closed-cycle circulatory system, there is a choice between a turbine-type pump and a pulsatile pump. Pulsatility in blood flow is something that we take for granted in healthy humans and animals, and is an inevitable consequence of nature's relative difficulty in synthesising a turbine-type heart, with Moazami et al. [3] contending that ''evolution has favored a pulsatile heart pump to be able to deliver the maximum flow at different levels of systemic vascular resistance, confer kinetic energy to the flow of blood past areas of stenosis and generate low shear stress on blood elements''. Thus, we may imagine that blood flow (and pressure) pulsatility is an unavoidable 'side effect' of an inability to synthesise a constant-flow pump. However, some evidence now suggests that pulsatility may have an important function in maintaining perfusion to ensure that there is an adequate supply of oxygen and nutrients to vital organs. For example [4] demonstrates that when a subject was dependant on a non-pulsatile heart, the perfusion to a number of capillaries was reduced significantly, but not restored to some capillaries when pulsatility was recovered.
In the case where healthy hearts become diseased or lose some of their pumping potential, patients may be fitted with an artificial heart, or left-ventricular assist device (LVAD). Given the artificial nature of these devices, there is a choice as to the use of pulsatile-or continuous-flow designs. In general, turbine-based models are favoured, motivated by the short lifetime and large abdominal cavity required by comparable pulsatile devices [5]. While the benefits of continuous-flow devices are also demonstrated by improvements in patient outcomes [6], there is still some concern regarding the long term physiological effects of diminished blood flow pulsatility.
To some extent, the effect of LVAD insertion on LF oscillations has been considered in the literature, with Cooley et al. [7] show that the presence of an LVAD significantly increases the presence of LF oscillations in R-R interval. While not specifically addressing the issue of LF oscillations, the studies in [8,9] examine broader aspects of non-pulsatile vs. pulsatile blood flow.
Other models for the neural baroreflex, combining both resistance and cardiac branches have been proposed, but employ different analytical tools, and pursue different objectives. For example, [10] employs cross-recurrence analysis to examine the coupling between heart-rate and vascular sympathetic tone, while [11,12] examine heart-rate variability for chaotic behaviour. The model of [13] is used to examine the effect of disease on cardiovascular autonomic regulation, and [14] examines the haemodynamic response to blood volume perturbations.
The studies in [15] and [16] are among the very few to give consideration to the role of pulsatility within the baroreflex loop. The baroreflex model of Ursino [15] includes a comprehensive circulatory model, with both sides of the heart simulated. It also models heart pulsatility, and provides simulation results and comparative experimental results, but conclusions relating baroreflex characteristics to pulsatility are limited to static loading cases. The study in [16] also considers pulsatility, but models pulsatility simply as a ramp signal, and considers the effect of pulsatility on the baroreflex in a qualitative, rather then quantitative way. This paper, through the development of a mathematical baroreflex model, combined with analysis techniques from the control sciences, examines the role of pulsatility in moderating the neural baroreflex and, inter alia, shows a relationship between BP pulsatility and LF (Mayer) waves. In particular, the analysis in the paper shows that, with reduced or absent pulsatility, baroreflex gain may significantly increase with potential consequences for long-term healthy BP dynamics. The analysis in the paper also shows that the increase in baroreflex gain is also generally accompanied by an increase in incidence of LF oscillations [17], which can potentially be used as a surrogate diagnostic tool for baroreflex gain.
This paper extends the work reported in a number of previous publications, notably [2], which examine the conditions under which LF oscillations take place, but only consider the peripheral resistance (PR) baroreflex. [18] considers both PR and cardiac branches, but ignores both BP pulsatility and arterial compliance. [19] articulates the effect of pulsatility on the baroreflex gain, but uses a simplified model of the baroreflex, containing only the PR loop and also ignores arterial compliance. Finally, [20] includes both cardiac and PR branches, along with pulsatility, but omits arterial compliance. The results in [20] are exclusively simulation based. In contrast, the current paper has a comprehensive model for the neural baroreflex, containing both cardiac and PR branches, includes arterial compliance and pulsatility, and develops an analytical solution for the presence of LF oscillations, which is shown to have a solution consistent with results from computer simulation of the model, and also those determined experimentally [21]. Furthermore, we include an analysis which examines the relative importance of individual cardiac (sympathetic and parasympathetic) branches and peripheral resistance branches in mediating LF oscillations. This is performed both through an analytical sensitivity analysis and simulation.

2.
Materials and methods

The neural baroreflex
In essence, the neural baroreflex activates smooth muscle in the peripheral resistance, and modulates heart rate, in response to disturbances in mean arterial pressure (MAP) or metabolic needs. In order to provide a feedback mechanism, BP is measured in the aortic arch and the carotid sinus, with the central nervous system (CNS) using these measurements to provide an appropriate control signal to the BP 'actuators' i.e. the PR and heart, via sympathetic and parasympathetic pathways. In our analysis, we assume that a notional MAP set point is present in the CNS [22], which may be varied, depending on metabolic needs, etc. We also assume that conduction velocity (dromotropy), contraction (inotropy), and relaxation (lusitropy) have more minor effects than cardiac rate (chronotropy), though they are also mediated through sympathetic and parasympathetic neural control [23].
The complete neural baroreflex model [18] is illustrated in Fig. 1, highlighting the baroreceptors, CNS, heart, peripheral resistance and arterial compliance [24] subsystems. The model is parameterised for a New Zealand white rabbit. Note that, while it is beyond the scope of this paper to provide overall b i o c y b e r n e t i c s a n d b i o m e d i c a l e n g i n e e r i n g 4 0 ( 2 0 2 0 ) 1 -1 3 model validation with experimental data, each validated model component has been taken from a reliable source.

Baroreceptors
The baroreceptors are simply modelled as a nerve conduction delay of t a (see Table 1), which also includes the afferent nerve delay [25]. Some researchers [26] have suggested that the baroreceptors may have high-pass frequency characteristics, though it is not entirely clear whether such characteristics are due to the baroreceptors themselves, or the CNS, since typical experiments measure the response in sympathetic nerve activity (SNA) to baroreceptor activation. In any event, the CNS block in our model contains a lead/lag component, articulating both high-pass and low-pass characteristics.

CNS
As discussed in Section 2.1.1, the CNS contains a lead/lag component G CNS (s) [18], shown in Eq. (1), along with a variable gain K c . K c may vary, depending on physiological condition (e.g. stress, haemorrhage, etc.), leading to the onset, or disappearance of LF oscillations in MAP.
The CNS also contains the activation characteristics, f pc , f sc , and f sr which articulate the nonlinear steady-state response of parasympathetic cardiac (pc), sympathetic cardiac (sc) and sympathetic peripheral resistance (sr) nerve signals, respectively, to deviations in MAP away from the set point, BP set . f pc , f sc , and f sr are all sigmoidal in form (see Fig. 2), each parameterised with a tan À1 ( ) function [27,28]: with parameters as shown in Table 2. Note that, in Table 2, no x* value is given. Given that we utilise a set point, the base (homeostatic) condition will be that the BP error will be zero (i.e. x * =0), with the y* values therefore setting the baseline tone of sympathetic and parasympathetic activity.

Heart
The inputs to the heart block are the sympathetic and parasympathetic cardiac signals from the CNS, with the output being blood flow rate. Separate dynamics operate on the sympathetic (G sc ) and parasympathetic (G pc ) channels [18], via: G sc ðsÞ ¼ k sc 1 1:29s 2 þ 1:29s þ 1 ; where k pc = 0.148 and k sc = 0.181, while there are also nerve conduction delays of t sc and t pc in the sympathetic and parasympathetic channels, respectively (see Table 1). The heart rate has been determined as a nonlinear function of both sympathetic and parasympathetic signals [29,30] from data recorded by Kawada et al. [31] and has been linearised [18] to the following affine function: where k s = 0.88, k p = À1.2 [18], and f o = 220 [32] denotes the base heart rate, with value f h = 215 beats/min for a rabbit. Finally, blood flow, q b , is determined from where V h is the heart stroke volume, given as 0.51 cc per kg of body weight, m r [33], where m r = 2.9 kg [34].

Peripheral resistance
The modelled peripheral resistance subsystem focuses exclusively on the neural control of PR, above a base resistance of r* (taken as r * =0.16 mmHg/(ml/min) [35]). Hormonal, metabolic, myogenic and paracrinal effects are ignored, given the short timescale of PR effectors considered in this study, and a single branch is used to represent the aggregated effects of innervated resistance in the kidney, muscles, gut, and skin. A sympathetic nerve conduction delay of t sr (see Table 1) is included in the PR branch, as well as a dynamic block, G sr (s) [34] relating the response in PR to SNA: where k r = 0.0005.

Arterial compliance
Arterial compliance is modelled using a first-order Windkessel model [24] containing resistive (R p denoting total PR) and capacitive effects, as shown in Fig. 3, and effectively dampening the pulsatile effects of blood flow, due to its low-pass frequency characteristics. A value of C = 0.000171 g À1 cm 4 s 2 [24] was employed. Considering that R p is a variable, but appears as a parameter of the Windkessel model, it is potentially problematic from an analytical viewpoint, so an approximation will be employed in Section 2.2.1.

Model simplification
One of the main aims of this study is to develop a set of analytical conditions for LF oscillations in the baroreflex. To this end, it is important to have a compact and systemoriented description of the system, and one which, ideally, places the system firmly within the framework of linear timeinvariant (LTI) analysis. This, inter alia, requires that the external pulsatile signal be absorbed into a system-type description, and that an effective linearised, but amplitude dependent (describing function), approximation be employed for representation of the nonlinear characteristics f pc , f sc , and f sr . Furthermore, the parametric variation in the Windkessel model is also simplified, to aid analysis.

Arterial compliance
A significant difficulty with the arterial compliance model, from an analysis perspective, is that the subsystem currently belongs to the class of linear parameter-varying (LPV) systems.
In particular, the relationship between q b and BP is: where R p (t) is a time varying parameter of the system. It is possible to retain the essence of the effect of a varying R p by using the true value of R p in the numerator of (8), while using a constant (mean) value, R p in the denominator of (8): where p b is the MAP signal derived from the classic Ohm's law relationship, excluding compliance, as: as demonstrated in Fig. 4. This is justified by the fact that the main effect is a gain change, due to the numerator instance of R p (note that the dc gain of (8) is R p (t)), while a variation in the dynamic response around the mean value of R p in the denominator (R p ) will be minimal. Though the system on the righthand side of Fig. 4 is still parameter varying, it is significantly easier to analyse. A steady-state analysis yields:

Pulsatility effects
The pulsatile component, P, of the BP signal is absorbed into the static nonlinear characteristics, f pc , f sc , and f sr to produce equivalent nonlinearities f e pc , f e sc , and f e sr . This provides a significant simplification by observing the frequency separation between the variations in MAP (which include LF oscillations) and pulsatility. In the 'equivalent nonlinearity' (EQNL) procedure, a system with an input x(t), containing a relatively LF input r(t) and a high frequency 'dither' signal d(t), can be represented as a modified nonlinear system [36], as shown in Fig. 5. Here, the procedure for the determination of an EQNL is outlined, developing a new nonlinearity (Nonlinear system B ( f e (x))), subject only to a LF signal r(t), from an original nonlinear system A ( f(x)) which is subject to LF (r(t)) and HF (d (t)) signals. For further details, see [36] and [19]. Since the EQNL is specific to the nature of the dither signal d(t) employed, a typical pulsatile BP signal from a rabbit is employed [37], as shown in Fig. 6.
With an original nonlinear system The output of the EQNL corresponding to the nonlinear system is: f ðrðtÞ þ qÞ pðqÞ dq (13) where p(q)dq is the probability that, for any time t, the dither signal d(t) lies in the range q to q + dq, and d(q) is the probability density function for the dither signal. The probability density function for the dither signal is: The parameters for the dither signal were obtained from a piecewise linear approximation of the clinical measurement of BP in the abdominal aorta of a rabbit [37]. One cycle of the data   was divided into 4 segments, three triangles and one flat line, with the fit, as shown in Fig. 6, determined manually. The specific parameters for each segment can be found in Table 3. It is found that the EQNL is not especially sensitive to the finer details of the dither signal, with the amplitude extremes and proportion of positive/negative transition the dominant influences.
Segments 1, 2 and 4 are represented by triangles with varying heights, A, and their respective EQNL outputs can be obtained by inserting the appropriate A into Eq. (15): Segment 3 is approximated by a flat line, resulting in a constant probability function and giving an EQNL output of: The total EQNL output is the sum of all four segment components, weighted by factors a i , corresponding to the relative portions of the total period t p À t o , from Fig. 6: where a 1 = 0.252, a 2 = 0.385, a 3 = 0.167, a 4 = 0.196. Note that Sa i = 1. In order to get the f ! f e for each of f pc , f sc , and f sr , Eq. (18) is used in conjunction with (16) and (17), with the appropriate b and y* values for y 1 , y 2 , y 3 , and y 4 taken from Table 2, as appropriate.
Note that, since the EQNL depends only on the relative areas of the piecewise approximations to the pulsatile signal segments, and the proportions of the period occupied by those segments, but not on the period of the dither signal itself (so long as d(t) is significantly HF, relative to r(t), for example by an order of magnitude, which would be a typical ratio between heart-rate and LF oscillation, across different species), a fixed frequency for d(t) (effectively the heart rate) can be used in subsequent analysis and simulation.

Describing function approximation
The static nonlinear characteristics, f pc , f sc , and f sr , are represented by their describing functions, simplifying the analytical manipulation of the functions f pc , f sc , and f sr . The describing function (DF) approximation essentially assumes a sinusoidal input to a nonlinearity (useful for the analysis of LF oscillations) and calculates the effective 'gain' of the system (which is dependent on amplitude of the input) with respect to the fundamental component of the output [38]. For the tan À1 ( ) function of Eq. (2) and Fig. 2, the DF has been calculated [27] as: where the input to the nonlinearity is x(t) = X sin(vt). The expression in (19) may be further simplified [27] for asymptotic values of X to: for large X and, for small X. The use of the DF in this application is justified, since the major assumption that the DF relies on, i.e. that output harmonics are significantly attenuated, is satisfied through G pc , G sc and G sr all having low-pass frequency response characteristics. Furthermore, we are primarily interested in the propagation of LF (fundamental) oscillations, rather than the existence of harmonic components. In any case, the dominance of the fundamental will be demonstrated in Section 3.1.2.

LF oscillations
Using the model of Section 2.1, LF oscillations can be analysed in two ways. The first, and most direct, approach is to simulate the model of Fig. 1 and examine for conditions under which oscillations in BP occur. While such simulation studies are useful, they provide little indication of the bigger picture, outside the set of specific tests performed. It would therefore be beneficial if an alternative analytical route could be pursued, employing the simplifications of Section 2.2.
In [2], it was possible to employ a DF approximation in combination with classical (linear) Nyquist stability analysis to determine conditions under which LF oscillations would occur. In that analysis, the gain of the CNS was seen to be instrumental in the modulation of LF oscillations, as will be shown in the current case. However, for the model of Fig. 1, such an approach is not possible, given the multiple branches and diverse dynamics in each path. Rather, an approach, inspired by [18], is employed, where a notional LF (sinusoidal) oscillatory signal is injected at x(t), and propagated around the baroreflex, to determine a set of conditions under which the original signal might arrive back at x 0 (t) i.e. the conditions for sustained oscillations, as shown in Fig. 7.
If sustained oscillations occur, assuming the fundamental component of x 0 (t) to be x 00 ðtÞ ¼ A 00 sinðv 00 t þ f 00 Þ þ X 00 , then the following must hold: f 00 ¼ 0; X 00 ¼ 0; Parameter corresponding to conditions on frequency, amplitude, phase and offset, respectively. We assume that x(t) has no dc component, since the system of Fig. 1 is assumed to be in homeostasis (see more on this in Section 3). Using the configuration in Fig. 7, the signal x(t) first passes through the three nonlinearities, f pc , f sc and f sr , resulting in the DF outputs: and The heart output flowrate, f 0 h , is now: where, N pc ¼ k p k pc jG pc ðv 0 ÞjDF pc A 0 ; (33) N sc ¼ k s k sc jG sc ðv 0 ÞjDF sc A 0 ; f sc ¼ ffG sc ðv 0 Þ À v 0 t sc ; f h ¼ k p k pc y Ã pc þ k s k sc y Ã sc þ f o : (37) or j ¼ N pc sinðf pc Þ þ N sc sinðf sc Þ: From (6), where Now propagating z sr (t) through the peripheral resistance branch, get where, The cardiac (43) and resistance (45) branches are combined via the simplified arterial compliance model as shown in Fig. 4: Note that, in (53), the amplitude of the harmonics of v 0 are assumed to be negligible compared with the fundamental, so just the fundamental component is retained. This is justified, as in Section 2.2.3, by the fact that the dynamics of G pc , G sc , and G sr are low-pass and therefore attenuate relatively high frequency harmonics, with further validation in Section 3.1.2.
The BP 0 signal is now fed back through the baroreceptor delay to the CNS dynamics: Since we have assumed, in (53), that the harmonics of BP 0 are negligible, x 00 in (57) clearly contains only the fundamental component of v 0 and the condition of (22) becomes redundant. The amplitude condition of (23), in residual form, from (57), becomes: while the phase and offset conditions of (24) and (25) become, respectively: and K c jG CNS ðv 0 ÞjfPB set À r q 0 b À kðv 0 Þg ¼ 0: Eqs. (58)-(60) describe the set of conditions of sustained oscillations with a fundamental frequency v 0 .

Solution to equations
From Eqs. (58)-(60), we can identify the cost function to minimise as: where where the weighting factors g i reflect the relative magnitude units of the r i .
The Levenberg-Marquardt algorithm [39] is used to solve the non-linear least squares problem represented by minimising (61). In order to investigate the sensitivity of the solution to the optimisation algorithm used, trust-region [40] and simplex [41] algorithms were also tested, giving consistent results.

Homeostatic analysis
In order to prove that the system, from a static balance point of view (homeostasis), is consistent, a steady-state analysis is performed. Under such conditions, the BP error D BP should be zero and all dynamic blocks reduce to their dc gain, as shown in Fig. 8. Also, the arterial compliance subsystem reverts to the 'Ohm's law' relationship of (10). Since D BP = 0, the outputs of f pc , f sc , and f sr are y pc *, y sc *, and y sr *, respectively, as given in Table 2.
Under these conditions, BP is 80 mmHg, which is consistent with the condition DBP = BP set À BP specified in Fig. 1, confirming homeostasis. While the specific value of Bp set = 80 mmHg is not crucial in confirming homeostasis, it has been chosen as typical of a normotensive rabbit (see, for example, the experimental data in Fig. 6).

EQNL calculations
Following the procedure in Section 2.2.2, the equivalent nonlinearities f e pc , f e sc , and f e sr are now determined for the nonlinear characteristics f pc , f sc , and f sr , shown in Fig. 9.
Note that a tan À1 ( ) fit has been generated to the EQNL in each case, so that a specific change in the linearised 'gain' at (x * , y *), due to the presence of pulsatility, may be calculated for each case. The results are shown in Table 4. It may be noted, from Fig. 9 that, while the fit of the tan À1 ( ) function to the EQNL is generally good in the region of the function origin, it is less good at the extremities. This is deliberate, in trading off error in central values for extreme values, within the limits of approximation of the fitted tan À1 ( ) functions. Specifically, the EQNL is precisely matched at the origin, so that a representative b value can be determined for the EQNL. If desired, a more precise parametric approximation can be determined using a combination of parametric functions, including an tan À1 ( ) and a saturation characteristic [19].

Determining a value for K c
The value for K c is not well defined, since any experimental evaluation would require that all nerve fibres at the outputs of f pc , f sc , and f sr are fully recruited, which is impossible to achieve in experimental practice. However, K c can be estimated by recognising that K c is primarily responsible for mediating LF oscillations [2], and it therefore seems likely that the brain modulates K c in response to different adverse physiological conditions e.g. hypoxia [42], haemorrhage [43], etc. In particular, given that LF oscillations are instigated/ accentuated during adverse physiological conditions, one possibility is to establish a value of K c in the 'normal' physiological state (including pulsatility) which is just below that required to initiate LF oscillations. This value is set at K nom c ¼ 7 for the (normal) pulsatile case, where a 10% increase results in the onset of LF oscillations. Therefore, a reasonably substantial effective increase in overall baroreflex gain (due to a lack of pulsatility or otherwise) is required to instigate LF oscillations.

Sensitivity calculations
Given that, in the model of Fig. 1, there are multiple paths via which the neural control of BP is effected, and BP oscillations may propagate, a reasonable question might arise as to the predominance of any of these branches. To examine the relative strength/gain of each of the branches of the model in Fig. 1, we perform a simple sensitivity study, which examines the raw (dc) gain from D BP to BP, from the simplified model of Fig. 8. While a dc analysis does not include the dynamic effects of G pc , G sc and G sr , their effects will be relatively consistent between paths, and their low-pass cutoff frequency is relatively high, compared to the LF oscillation's frequency of 0.33 Hz, as are the arterial compliance dynamics.
From Fig. 8, the following may be easily calculated: which articulates the sensitivity of BP to D BP along the peripheral resistance path, and which articulates the sensitivity of BP to D BP along the complete cardiac path (pc & sc), where R p is given from Eq. (11) and The individual sensitivities along the parasympathetic and sympathetic cardiac branches can also be evaluated, respectively, as: and

LF oscillations
LF oscillations are now initiated by an increase in K c by 25% over the nominal value calculated in Section 2.4.3, consistent with the data in [27] showing a 25% increase in b in the change from normoxia to hypoxia (10% O 2 + 3%CO 2 ). In our simulation, this will be effected by an increase in K c from its nominal value of K c = 7 to K 0 c ¼ 1:25K nom c ¼ 8:75. This is confirmed by both the analytical solution from Eqs. (58)-(60), and computer simulation of the model in Fig. 1, using the Matlab/Simulink platform. The most useful solution parameters from the minimisation of (61) are v 0 and A 0 , which give the frequency and amplitude of the LF oscillations, respectively. The examination of A 0 may be used as binary indication of oscillations, where a value of A 0 = 0 indicates no oscillations. Fig. 11 shows the change in oscillation amplitude (A 0 ) with variation in K c . Note the expected threshold at K c = 7.7 for the normal (pulsatile) case.  Table 5, showing good agreement (deviation is 5.6% for the analytical solution and 2.3% for the simulation) with the experimental LF BP oscillation frequency of 0.3 Hz reported in [21].

3.2.
Removing pulsatility Under certain medical interventions, blood pulsatility may be reduced or eliminated, for example by insertion of an artificial (turbine-type) heart, or the inclusion of a LVAD. The effect of pulsatility on the neural baroreflex can be examined both via manipulation of the baroreflex functions f pc , f sc , and f sr , or by inclusion/exclusion of the pulsatility component P in the simulation.   Some initial indication can be had from examination of Fig. 9. In particular, with reference to Table 4, it is clear that the addition of pulsatility causes a significant reduction in overall baroreflex gain, demonstrated by the decrease in the b xx from their 'original' values to their EQNL values, with a gain reduction of up to 55% (for b pc ). The change in baroreflex gain, resulting from an absence of pulsatility, is confirmed by the simulation, where Fig. 14 shows the appearance of LF oscillations, following the removal of pulsatility (Case 2), consistent with experimental evidence [17]. Pulsatility can therefore be represented, for simulation purposes, either by including a pulsatile component P in the BP signal, as shown in Fig. 1 (Case 1), or by omitting P but using the EQNL values for the b xx (Case 3). Table 6 shows 3 cases examined, including both pulsatility, using both an explicit P component or the EQNL values for the b xx , and a non-pulsatile scenario, where original b xx values are used and P is omitted. Fig. 11 shows the variation in oscillation amplitude (A 0 ) with variations in K c for the non-pulsatile case (red trace). Note the significant reduction in the threshold value of K c . In Fig. 14, the trace for Case 1 (see Table 6) contains the pulsatile (green line) BP which is not completely dampened by arterial compliance, but confirms consistence with the trace for Case 3, i.e. the absence of LF oscillations.

Predominance of resistance or cardiac branches
Following the sensitivity functions developed in Section 2.5, Table 7 shows the enumerated values for the sensitivity functions in (62)-(66). The sensitivity values suggest a predominance of the cardiac branch over the resistance branch, while the parasympathetic cardiac is favoured over the sympathetic cardiac. Simulation studies were undertaken to examine the relative importance of various branches in mediating LF oscillations and Table 7 also shows the achieved oscillation amplitude (OA) in BP and the associated oscillation frequency, for removal of various branches (for example, in the pc column, only the pc branch is active). The simulation results are quite revealing. For a value of K c ¼ K 0 c (with pulsatility included), no single branch is capable of supporting LF    Table 6.  b i o c y b e r n e t i c s a n d b i o m e d i c a l e n g i n e e r i n g 4 0 ( 2 0 2 0 ) 1 -1 3 oscillations. As K c is increased (to 2K 0 c , and subsequently to 3K 0 c ), individual branches become capable of sustaining LF oscillations but, crucially, none of the achieved LF oscillations frequency values correspond with experimental evidence [21]. This is primarily due to the variety of conduction and phase delays in each individual branch, while the full model contains a mixture of these dynamical components that appears to correctly articulate the oscillations observed experimentally. Overall, there is a strong suggestion that both cardiac and resistance branches are crucial in the support of LF BP oscillations. Interestingly, Liu et al. [44] conclude that the cardiac branch is relatively unimportant in the dynamic regulation of BP and the mediation of LF oscillations in MAP. However, their results do show an 8 dB (equivalent to a factor of 2.5) drop in gain between 0.1 and 0.3 Hz due to vagotomy and b 1 -receptor blockade (effectively inhibiting neural control of cardiac output), which this study has shown to be significant in mediating LF oscillations, given that a 10% increase in baroreflex gain could, potentially, instigate LF oscillations.

Discussion
The paper presents a model for the neural baroreflex focused on the components that mediate LF oscillations in BP. A number of analytical tools are employed to render the model amenable to an algebraic solution for LF oscillations conditions, which gives some insight into the role that baroreflex gain plays in mediating LF oscillations. In particular, changes in the 'gain' of the central nervous system are shown to be a key factor, which can result from external stimulus or internal stress (e.g. haemorrhage, hypoxia, etc.). In particular, the use of the 'equivalent nonlinearity', in representing pulsatility, reveals a significant increase in overall baroreflex gain as a result of a loss of pulsatility, which may have important implications for patients in receipt of artificial (turbine-based) hearts, or LVADs. Such an increase in baroreflex gain is shown to precipitate an elevation in LF oscillations, which is borne out by experimental evidence [17]. In addition, the increase in experimentally measured LF oscillations following LVAD insertion documented in [7] in R-R is consistent with the analytical results developed in this paper relating to the role of the heart in mediating LF oscillations and an increase in baroreflex gain following a drop in pulsatility. However, it should also be noted that the study of Cooley et al. [7] did not report an increase in LF MAP oscillations, opening the possibility that the steady flow from the autonomous LVAD may have dominated the native cardiac output. While experimental evidence therefore exists which ties pulsatility decreases (LVAD insertion) to increases in baroreflex gain, and LVAD insertion to the increased presence of LF oscillations in R-R interval, the mechanisms by which these are linked has, to date, remained unclear. We hope that the unifying framework in this paper provides a plausible explanation for these experimentally observed phenomena. Overall, the model, and the analysis contained in the paper, provides a toolbox for the analysis of the interplay between baroreflex gain, LF oscillations and blood pulsatility. One of the important conclusions is that LF oscillations are mediated by a combination of cardiac and peripheral resistance branches, while establishing the key role of baroreflex gain as the mechanism by which LF oscillations are modulated. In turn, the baroreflex gain is seen to be substantially influenced by pulsatility and physiological stress.

Declarations of interest
None declared.