A stochastic delay differential model of cerebral autoregulation.

Mathematical models of the cardiovascular system and of cerebral autoregulation (CAR) have been employed for several years in order to describe the time course of pressures and flows changes subsequent to postural changes. The assessment of the degree of efficiency of cerebral auto regulation has indeed importance in the prognosis of such conditions as cerebro-vascular accidents or Alzheimer. In the quest for a simple but realistic mathematical description of cardiovascular control, which may be fitted onto non-invasive experimental observations after postural changes, the present work proposes a first version of an empirical Stochastic Delay Differential Equations (SDDEs) model. The model consists of a total of four SDDEs and two ancillary algebraic equations, incorporates four distinct delayed controls from the brain onto different components of the circulation, and is able to accurately capture the time course of mean arterial pressure and cerebral blood flow velocity signals, reproducing observed auto-correlated error around the expected drift.

Introduction

Autoregulation of blood flow denotes the intrinsic ability of an organ or a vascular bed to maintain a constant perfusion in the face of blood pressure changes [1]}. In particular, cerebral autoregulation (CAR) denotes the ability of the circulation to adapt to variations in hydrostatic pressures by means of compensating changes in heart rate, peripheral vascular resistances and venous capacitance, so as to maintain constant and adequate perfusion to the brain.

For the assessment and prognosis of the progression of some diseases, such as cerebrovascular accidents, Alzheimer and others [2-9], the evaluation of the adequacy of cerebral autoregulation provides important information. Compromised cerebral hemodynamics, such as reduced vasodilation, reaction to CO2 and other stimuli, may in fact, be related to reduced post-stenotic perfusion pressure.

Since the early 1990’s several scholars [10-19] have analyzed mathematically the vascular compensation mechanisms responsible for the regulation of cerebral blood flow.

A comprehensive, deterministic mechanistic model of the interplay of cerebral blood flow, cerebral blood volume, intracranial pressures and regulatory mechanisms, had been proposed by Ursino and Lodi [13]. This model, using measured arterial pressure (AP) as driving or input function, proved adequate to reproduce the observed cerebral blood flow velocity (CBFV) profile in human subjects undergoing transition from sitting- to—standing [20].

The Ursino-Lodi model was subsequently used [21] to analyze non-invasive measurements of cerebral blood flow velocity and arterial blood pressure on a healthy young subject undergoing postural changes. The Ursino-Lodi model was also used [21] to study variations in cerebral autoregulation between populations, even though the difficulty in estimating the many model parameters from non-invasive measurements of cerebral blood flow velocity and arterial blood pressure was underscored.

The first major limitation of this approach is that the time course of arterial pressure is simply taken as input function (thereby not explaining the determinants of the variations of arterial pressure itself) and so the time course of arterial pressure is not reproduced from hypotheses on the hemodynamic changes induced by the orthostatism. The second limitation is that the identification of the adequacy of cerebral autoregulation is difficult in a given single patient, because a mechanistic model requires many free parameters to be estimated from data and this limits the clinical applicability of the method. The third limitation is the inadequacy of the deterministic modeling of pressure and flow variables in capturing their irregular variations over time, produced by accidental muscle contractions, station adjustments, hormonal and neurological oscillations etc.

Following the same approach of Ursino and Lodi, Olufsen et al. [22] constructed a model that proposes to explain the morphology of both the arterial pressure signal and that of the corresponding cerebral blood flow velocity. The expected new (apparent) steady state reached after standing was correctly reproduced by this model. However, these authors pointed out that the predicted drop in cerebral blood flow velocity after standing was substantially smaller than observed. This shortcoming was mitigated by the same authors in a subsequent work, in which they detailed the construction of a comprehensive eleven-compartment model, based on a reconstruction of the flows and resistances in idealized, representative segments of the circulation. This improved model captured satisfactorily the behavior of the observed signals, even in the transient phase before establishment of a new equilibrium. The three most important model’s characteristics [18] are: the inclusion of non linear functions describing resistances of the large systemic arteries as functions of pressure; the inclusion of autonomic regulation; and the inclusion of an empirical model describing the dynamics of cerebral vascular resistance.

The plausible hypotheses for autonomic and cerebrovascular regulation on which the model is based, allow a good quantitative match with physiological observations, even if this match is obtained through a rather complex set of equations and parameters.

Chiu et al. [23] and Liau et al. [24] focused their studies on diabetic patients and applied time-domain cross-correlation as a technique to assess the relationship between blood pressure and cerebral blood flow velocity signals after postural changes. Liau et al. [25] extended the application of the same techniques to stroke patients and showed that mean arterial blood pressure changes in response to postural challenges were reduced in stroke patients.

The urgency of devising protocols to validate the reproducibility and ranges of the dynamic parameters extracted from these models, and the importance of developing multivariate models that take into account time-varying parameters was stressed by Panerai et al. [26].

Very recently a new empirical model has been presented in [27]. The model aims at representing the viscoelastic response of tissues (blood vessels, cerebral arteries…), which exhibit a continuous relaxation in response to stress and in the specific case to postural change from sitting to standing. The presented model is a mechanical analog model for predicting the CBFV in response to AP changes, incorporated in the model as an input function. While the strength of the model derives from its simplicity (it has only four parameters to be estimated), its major shortcoming lies in the lack of representing the physiological mechanisms intervening during the experiment.

With the above considerations in mind, a simpler model (if it could correctly reproduce observations) would clearly be of interest in the quest for practical assessment of the efficiency of CAR. At the same time, evident autocorrelated departures from the expected, smooth signal, indicate that stochastic elements are involved (stemming in all likelihood from moment-to-moment variations in muscle activity, hormonal concentrations, sympathetic/parasympathetic tone and other influences). Finally, it would therefore also be of interest to consider more prolonged variations, such as the likely continuing, slow increase in arterial pressure, after fast compensation, instead of assuming attainment of the Steady State within a few tens of seconds after standing.

The goal of the present work is therefore to describe a simple stochastic delay differential equations model of cardiovascular regulation during toilt-table or sitting-to-standing maneuvers. We aim to show that this model is able to reliably reproduce the time course of mean AP and CBFV time courses after changes in posture, inclusive of some auto-correlated oscillations around the expected signals.

Materials and Methods

Model description

In the present section a stochastic delay differential equations (SDDE’s) model for cerebral autoregulation is proposed. The model is composed of four compartments and each one of the four equations (1–4) represents one of the main components of overall cerebral regulation. The stochastic component of the model appears in Equation (1).

The model equations are as follows: with the initial conditions C(0) = C0, A(0) = A0, H(0) = H0, F(0) = F0. C (cmH20) is central venous pressure, A (mmHg) is arterial blood pressure, H (Hz) represents heart rate (as an index of sympathetic activity), and F (mL/sec) is cerebral blood flow velocity.

Moreover, two algebraically defined variables are introduced: the variable B (mmHg), representing brain arterial pressure, and variable R (mmHg/mL/sec), representing Peripheral Vascular Resistances as a function of brain arterial pressure . R is modeled as a decreasing function of in order to represent the mechanism by which, if brain arterial pressure increases, peripheral arterioles dilate and peripheral resistance falls:

Therefore, in the present formulation, Peripheral Resistance represents one of the control mechanisms of the cerebral autoregulation system, driven, with delay, by cerebral pressure.

Brain arterial perfusion pressure B differs from arterial blood pressure A by a pressure delta determined by the distance from the head to the heart (dist), multiplied by the sine of the tilt angle, that is the radians upright from the supine position: and with where ρ(t) is the tilt angle at time t. In the present formulation we have set:

Notice that the model, as formulated, can be used to fit data from generic tilt-table maneuvers (with parameter ρ in eq. 1 taking values in [0, π/2], depending on the final tilt-angle) and also sitting-to-standing maneuvers (with the parameter ρ = π /2 or, which is the same, with sin(ρ) = 1, expressing a simple Heaviside step from one position to the other). In the following we will describe the (continuous) geometry relative to tilt-table experiments, since sitting-to standing experiments can be simply represented as theoretically instantaneous transitions between two positions, arbitrarily indexed by sin(ρ) = 0 and sin(ρ) = 1.

The delta pressure in cmH2O from the head to heart is converted to mmHg dividing by the conversion factor 1.36.

The first equation (Equation 1) describes the variation of central venous pressure (CVP, in the model indicated with C) over time. In order to allow the heart to pump sufficient blood to maintain cerebral perfusion pressure, and hence cerebral blood flow, it is supposed that compensating mechanisms act to maintain a sufficient level of CVP and ventricular filling. CVP is at steady state when the tendencies to increase CVP, necessary to maintain cerebral perfusion pressure, and to decrease CVP, when a sufficient level of brain perfusion pressure is obtained, equilibrate.

For ρ = 0 the subject lies supine and the ‘elastic’ CVP target value is c , but when the subject undergoes a head-up tilt experiment, ρ increases proportionally to the tilt angle and the central venous pressure drops by a hydrostatic factor equal to . After the maneuver, CVP increases again: the gain in CVP is due to both reduced loss and concurrent drop of suppressing (delayed) brain perfusion pressure (the multiplicative second term). After postural change the control mechanism restores CVP, allowing cardiac filling to be brought back to normal again. The system which drives CVP trend is subject to random fluctuations due to different mechanisms such as muscular contraction, respiration, deglutition and other uncontrolled factors, represented in the model by the random system noise σW(t).

The second equation (Equation 2) describes variations of arterial blood pressure. These variations depend linearly on heart rate H, on CVP (these two multiplied together determine cardiac output) and on peripheral resistances R.

The third equation (Equation 3) describes the equilibrium between spontaneous increase and (delayed) brain perfusion pressure-dependent suppression of heart rate H. In fact, when the brain enjoys abundant perfusion, sympathetic tone may decrease, and heart rate may also correspondingly decreases [28].

The fourth equation (Equation 4) describes variations of cerebral blood flow velocity F as depending on current brain perfusion pressure B according to a Michaelis-Menten relationship: cerebral flow increases with increasing B, reaching a maximum level (in order to represent self-protecting cerebral vasoconstriction mechanisms, which may be missing in some pathological conditions). The parameter B 50 represents instead the brain perfusion pressure level at which a half-maximal rate of flow increase is obtained. The spontaneous decay of the flow is represented by the second term on the right hand of the equation, where the variable F is delayed, .

All the delayed variables in the model appear with a tilde hat on the top of the variable. In the present model all the considered delays are distributed: the influence of the variable on the dynamics of the system at a certain time t is based on previous values assumed by the variable until time t, weighted with a delay kernel function according to the following definition (Y indicates a generic variable): with

The exponential delay kernel ω(s) is parameterized by the parameter α, which is greater than zero: for α large, recent values carry more weight. For each of the delayed variables considered in the model (,,,) a different parameter α of the kernel function is considered.

Fig. 1 shows the schematic representation of the model.

Fig 1

Block diagram of the ESDDeCAR-02 model.

Blocks in continuous lines represent the model state variables; blocks in dashed lines represent the two defined algebraic variables.

The system admits a continuum of different equilibria, indexed by ρ, when the model drift (that is the deterministic component only) is considered. Each of the equations have two steady state conditions, before the maneuver and at infinite time after its completion:

While H 0 and R 0 are supposed to be assigned (fixed).

From the steady state conditions the following parameters are also determined:

Table 1 reports the model parameters along with their description, units of measurement and values used in the simulations.

Table 1

List and description of parameters values used in the simulation.

Parameter	Unit of Measurements	Description	Value
k c	/sec	Rate of central venous pressure increase	1
C tgt	cmH20	Central venous pressure target value	determined
CΔmax	cmH20	Proportional constant drop of central venous pressure target following the tilt experiment	determined
k xca	/sec/mmHg	Central venous pressure decrease rate per mmHg of brain arterial pressure	0.01
ρ	radians	tilt angle	0,π/2
dist	cm	Distance between head and heart	30
C 0	cmH20	Central venous pressure value before the tilt experiment	8
C end	cmH20	Central venous pressure value at equilibrium after the tilt experiment	4
k arch	mL/cmH2O/bpm	Third-order arterial pressure increase rate	0.2
k xa	/sec	Arterial pressure decrease rate	0.5
A 0	mmHg	Arterial pressure value before the tilt experiment	determined
A end	mmHg	Arterial pressure value at equilibrium after the tilt experiment	determined
k ha	/sec/mmHg	Heart rate decrease rate per mmHg of brain arterial pressure	determined
k h	bpm/sec	Spontaneous increase of heart rate	0.05
H 0	bpm	Heart rate value before the tilt experiment	60
H end	bpm	Heart rate value at equilibrium after the tilt experiment	determined
kfmax	mL/sec/sec	Maximal increase in cerebral blood flow velocity	28.9
k xf	/sec	Cerebral blood flow velocity decrease rate	0.5
F 0	mL/sec	Cerebral blood flow velocity value before the tilt experiment	determined
F end	mL/sec	Cerebral blood flow velocity value at equilibrium after the tilt experiment	determined
ν	#	Rapidity with which cerebral blood flow velocity increase reaches its maximum with increasing B	1
B 50	mmHg	Brain arterial pressure level at which an half-maximal increment of cerebral blood flow velocity is obtained	40
B 0	mmHg	Brain arterial pressure value before the tilt experiment	determined
B end	mmHg	Brain arterial pressure value at equilibrium after the tilt experiment	determined
R inf	mmHg/mL/sec	Minimum value obtained for peripheral vascular resistance	0
R max	mmHg/mL/sec	Maximum value obtained for peripheral vascular resistance	determined
λ	/mmHg	Rate of decay of peripheral vascular resistance with increasing brain arterial pressure	0.035
R 0	mmHg/mL/sec	Peripheral vascular resistance value before the tilt experiment	0.47
R end	mmHg/mL/sec	Peripheral vascular resistance value at equilibrium after the tilt experiment	determined
α i	/sec	Delay kernel rate constant for delay kernel i	0.01,0.1,0.2,0.6
t 0	sec	Time of tilt maneuver	30
σ	cmH20/sec1/2	CVP volatility	0.5

The model was implemented in Matlab using a fixed-step, fourth-order Runge-Kutta numerical integration scheme [29].

Qualitative analysis of the model solutions

The present deterministic qualitative analysis refers to the drift part of model (1–8) only (without the additional stochastic term σdW in (Eq. 1)).

The following subsections investigate the conditions for which the solutions are positive and the model is persistent, report the proof that a unique positive equilibrium point exists and provide asymptotic local stability analysis.

Positivity of the model solutions and model persistence

Theorem 1. The system (1), (2), (3) and (4) admits positive solutions for any positive initial condition.

Let C(0)>0. According to the continuity of the solution of a differential equation, C(t) would become non-positive if there existed a t*>0 such that C(t*) = 0, C(t)>0 for any 0 ≤t

This proves that C(t)>0 (it never vanishes and it is always positive). Similarly it can be proven that, if H(0)>0, also H(t) never vanishes and is always positive: let H(0)>0 and assume that ∃t*>0 such that H(t*) = 0 and H(t)>0 for any 0≤ t

For (equation 2) let A(0)>0. A(t) would become non-positive if there existed a t*>0 such that A(t*)>0, for any 0 ≤ t

This proves that A(t) > 0.

For (equation 3) let F(0)>0. Again,F(t) would become non-positive if there existed a t*>0 such that F(t*) = 0, F(t)>0 for any 0 ≤ t

At equilibrium, , from which it follows that since A(t*)–dist/1.36 > 0 due to physiological considerations (A(t) is larger than 29 mmHg in all humans still alive). This means that at equilibrium F(t) is always strictly positive.

In the transient, the behavior of F(t) depends on . Since the analytic investigation of the conditions, for which the positivity of F(t) always holds, is very cumbersome, a numerical exploration of the behavior of F(t), as depending on a range of values of the parameters α and R o, has been performed.

Fig. 2 show the behavior of the minimum of the cerebral blood flow velocity time course over the duration of the experiment. A region of the plain is explored, determined by a range of values of α from 0.1 to 1.5 and by a range of R 0 from 0.1 to 1.5 (mmHg/mL/sec). The color coding shows with darker color smaller but positive values of the min(F), while white denotes negative values of CBFV, which are of course physiologically unacceptable.

Fig 2

Behavior of the minimum of the cerebral blood flow velocity (CBFV) trend over time varying parameters α and R 0.

White region represents negative (unacceptable) values of the minimum, graded colored region represents the (positive) value of the minimum of the CBFV over the perturbation time interval.

Theorem 2. The system of equations (

Recall that a model is persistent if there exists a pair of positive real numbers (m, M) such that: for each component X of the state vector. Denote:

The proof is achieved by proving the following eight statements:

Step 1. In order to show the boundedness of the evolution of central venous pressure, assume that C = +∞

Since C is by definition differentiable and since C = limsupt→+∞, we can define a sequence of time points {t } such that:

t →+∞ as n→+∞

C(t )→+∞ as n→+∞

If we consider (equation 1) at each of time points t , we find that which is a contradiction, so that C < + ∞.

Step 2. In order to show the boundedness of the evolution of arterial blood pressure, heart rate, cerebral blood flow velocity, assume that A = H = F = + ∞. From the same considerations as in (27) it should be:

A(t )→ + ∞,H(t )→ + ∞,F(t )→ + ∞ and .

However: which are contradictions, so that A < +∞, H < + ∞,F < + ∞.

Step 3. From Step 1, it follows that C ≤ C < + ∞. Since C is by definition differentiable and since C liminf C(t), we can define a sequence of time points {t } such that: which means:

t → +∞ as n→ +∞

C(t )→C as n→ +∞

According to inequality (28), Cm>0 because

Step 4. From Step 2, it follows that H ≤ H <+∞, A ≤ A <+∞ and. F ≤ F <+∞ From the same considerations as in step 3 it follows that, for each state variable H, A and F, it is possible to define a sequence of time points {t }{t } and {t } such that: and which means:

Existence of positive equilibrium points

It is shown here that each equation of the model (1–4) admits positive equilibrium points over a set indexed by ρ.

Theorem 3. System (1–4) has a unique positive equilibrium point for each ρ.

As derived above, for t = 0 (a(ρ(t)) = 0 the equilibrium point is given by (C 0,A 0,H 0,F 0), and for t = t (a(ρ ) = 1) the equilibrium point is given by (C ,A ,H ,F ).

For t≠0 and t≠t the equilibrium points of (1–4) satisfy the following algebraic equations: from which it follows that:

Substituting (33–35) and (37–38) in (36): Notice that and

This condition is satisfied ∀ρ(t), in fact

and due again to physiological considerations it is always A≫29mmHg.

From the above it follows that φ (·) is a decreasing function for positive argument, starting from a positive value at zero, and hence that it has at most one positive root. The existence of such a root is ensured by the limit:

Asymptotic local stability analysis

In order to study the local stability of the equilibrium points (indexed by ρ), the system (1)-(4) is linearized around its generic equilibrium point (C A H F )

Adopting a change of variable, the delayed variable can be written as:

Denoting:

The following 6-dimentional ordinary differential system is obtained:

Taking into account the extended system (39), the local asymptotic stability of the equilibrium point is achieved if and only if all the roots of the characteristic polynomial associated with the Jacobian computed in X have negative real part. The elements in (40) are defined as follows:

After computation, the characteristic polynomial is: where p 1 = a, p = c, p = e. The polynomial (41) is zero if and only if (p – λ) = 0∀i and (λ3 -(g+h)λ2 + gh λ-f) = 0.

It follows that , where a, c and e are negative values.

For the study of the sign of the eigenvalues associated with the polynomial (λ3 -(g+h)λ2 + gh λ-f) the Routh-Hurwitz criterion is used. The first column of the Routh-Hurwitz table is:

and we have to determine for which conditions on parameters g and f the elements of (42) assume positive values. These conditions assure that the roots of the polynomial all have negative real part:

-2g is always positive, since g<0.

which is true if

-f is always positive, since f<0.

The condition in ii) translates into the standard requirement for delayed systems that, in order to guaranteed model stability, a sufficiently small (average) delay on the considered delayed state variable F is necessary: in fact, the larger the parameter α , the smaller the associated average delay.

Results

Fig. 3 shows predicted time courses of arterial pressure (panel A) and of cerebral blood flow velocity (panel B), as derived by the stochastic (dashed lines) and deterministic (continuous line) model (1–8). These time courses are directly comparable (both in timing and in amplitudes) with non-invasive observation data for sitting-to-standing maneuvers reported in the literature (e.g. Fig. 2 of [20]). Predictions in Fig. 3 and Fig. 4 were in fact obtained by calibrating parameter values for a sitting-to-standing experiment: parameter values are reported in Table 1.

Fig 3

Simulated model behavior in normal physiological conditions.

Panel A: time course of model-simulated arterial pressure (mmHg) (dashed line: stochastic model, continuous line: model drift). Pressure drops of about 40 mmHg after standing, with a nadir at nearly 5–8 seconds post-maneuver, climbs back within another 10 to 15 seconds with an overshoot of approx 20 mmHg before stabilizing. Panel B: time course of model-simulated cerebral blood flow velocity (CBFV,mL/sec) (dashed line: stochastic model, continuous line: model drift). CBFV drops of about 30 mL/sec after standing, with a nadir at 5 to 8 seconds post-maneuver, climbs back within another 10–15 seconds with an overshoot of approx 5 mL/sec before stabilizing.

Fig 4

Simulated model behavior in abnormal physiological conditions obtained by varying R 0.

Panel A reports model-simulated Resistance (R, mmHg/mL/sec) as an algebraic function of brain arterial pressure (mmHg) for different values of R . Moreover for each value of R the equilibrium points of brain arterial pressure (and corresponding resistance) are shown. Panel B and C report the time course of model-simulated delayed brain arterial pressure (mmHg) and of cerebral blood flow velocity (mL/sec) respectively, in correspondence of the different values of R 0. Panel D reports delta CBFV (difference between the lowest value of CBFV over time and F 0) as a function of R 0.

In order to explore the behaviour of the model in medically meaningful situations, a series of simulations have been performed (Fig. 4, Fig. 5), making some parameters vary, one at the time, from normal to abnormal values. These simulations associate the variation of model parameters, as they could be physiologically postulated in relation to some dysfunction, with the time courses of pressures and blood flow velocity predicted by the model. Some of parameters affect more than others the behaviour of brain arterial pressure (BAP) or arterial pressure (AP) and cerebral blood flow velocity, as shown in Fig. 4 and Fig. 5.

Fig 5

Simulated model behavior in abnormal physiological conditions obtained by varying some parameters.

Panel A reports the Time course of model-simulated arterial pressure (mmHg), on the left, and the model-simulated cerebral blood flow velocity (CBFV, mL/sec), on the right, varying the parameter k . The continuous line describes the trend in the physiological situation (k = 0.5), the dashed black line describes the trend with k = 0.7 whereas the dashed gray line describes the trend with k = 0.3. Panel B reports the time course of model-simulated arterial pressure (mmHg), on the left, and model-simulated cerebral blood flow velocity (CBFV, mL/sec), on the right, varying the parameter k . The continuous line describes the trend in the physiological situations (k = 0.05), the dashed black line describes the trend with k = 0.07 and the dashed gray line describes the trend with k = 0.02. Panel C reports the time course of model-simulated arterial pressure (mmHg), on the left, and model-simulated cerebral blood flow velocity (CBFV, mL/sec), on the right, varying the parameter k . The continuous line describes the trend in the physiological situations (k = 0.01), the dashed black line describes the trend with k = 0.06 whereas the dashed gray line describes the trend for k = 0.005.

For instance, Fig. 4 displays the time courses of BAPP and CBFV after sudden orthostatism, causing intravascular volume depletion, in correspondence of different values of the parameter R 0 (equilibrium resistance before the maneuvre). Panel A shows that decreasing R 0, i.e. decreasing the equilibrium Peripheral Vascular Resistance, determines lower equilibria for brain arterial pressure. Panel B reports the time course of delayed brain arterial pressure for different values of R 0 and shows that, while the pattern remains substantially unchanged, the lowest R 0 value determines the curve at the bottom of the graph, which fails to rise up to the considered reference values. Panel C reports the trend of CBFV: despite lower values attained during the transient period for smaller values of R 0, once the equilibrium is achieved the trend is not much altered, given the concurrent action of the other control mechanisms.

In panel D the greatest attainable delta in CBFV (the largest CBFV drop occurring in the transient phase), is reported as a function of R 0, showing that cerebral blood flow velocity undergoes greater decrements in correspondence with smaller resistances.

Fig. 5 describes the behavior of the model corresponding to variations of the parameters k ,k and k in panels A, B and C respectively. The term represents the local vasoconstriction, which is the fastest considered control mechanism: an increase in k determines a faster decrease, but a more rapid increase (due to the introduction of the delay) in CBFV with respect to the standard situation. Further, the achievement of the equilibrium state (panel A) is slower (more oscillating). A decrease in the parameter, of course, produces a smaller and slower fall of CBVF along with a slower achievement of the physiological values, but a more rapid establishment of the equilibrium, with fewer and shallower oscillations.

The other two parameters represent the other two control mechanisms: on heart rate (k ) and on venous capacitance (k ). Variations of these two parameters (that is variations of the efficiency of the sympathetic system and variations of central venous pressure recovery) influence very little both AP and CBFV (see panels B and C in Fig. 5).

Discussion

Since cerebral autoregulation (CAR) indicates the ability of the circulation to adapt to variations in hydrostatic pressures in order to maintain constant and adequate perfusion to the brain, the study of the mechanisms and of the overall adequacy of cerebral autoregulation has meaningful implications for the assessment of the clinical severity of degenerative and cardiovascular diseases, such as the clinical states associated with chronic hemodynamic compromise (e.g. obstructive carotid artery disease). In these cases detection of impaired cerebral autoregulation might help to identify patients at risk, as already shown for cerebrovascular reserve capacity [4,30].

Assessing the adequacy of cerebral autoregulation is however not immediate. Empirical indices of CAR efficiency such as ARI, autoregulatory index and ARMA-ARI, autoregressive-moving average, [31] or indices derived from spectral analysis of the oscillatory pulse signal, for instance G / G , the ratio of respiratory gain and the gain of the first cardiac harmonic [32], have been proposed in the past to offer semi-quantitative indications of a generically better or worse clinical situation. These indices have however a limited interpretability and their connection with well established physiological quantities is conceptually rather labile. Mathematical models of the response of the cerebral and systemic circulations to perturbation maneuvers could offer the opportunity to directly identify and quantitatively describe the key regulatory steps in the overall circulatory dynamics.

Mathematical models of the cardiovascular system and of cerebral autoregulation [10-19] have been developed over several years in order to describe the time courses of pressure and flow subsequent to postural changes. In general, many such models suffer from either one (or both) of the following shortcomings: they may be very large, consisting of many state variables bound to each other by interacting feedback loops; or, they may fail to fit available observations well. While the second shortcoming clearly determines major problems when foreseeing a possible clinical application of the model in question, the first shortcoming also gives rise to practical difficulties, connected with the difficult identifiability of overparameterized models from routine clinical data sets. In fact, the need exists, according to Panerai et al. [26], to develop multivariate models that take into account time-varying parameters with the aim to studying protocols to validate the reproducibility and ranges of the dynamic parameters of cerebral autoregulation.

The combination of the above considerations makes it so that it would be desirable to have a mathematical model, mechanistically representing accepted physiological phenomena, with as simple a structure as possible (but no simpler [33]) hence with as few free parameters to be estimated from data as possible, and still able to fit data well.

Attempting to respond to these requirements, the present work details a simple stochastic delay differential equations model of cardiovascular regulation during the postural perturbation maneuver (ESDeCAR-02). The model does in fact provide a realistic mathematical description of cardiovascular control, which may be fitted onto non-invasive experimental observations after postural changes, and is shown to replicate well observations performed during a sitting-to-standing manoeuvre.

The ESDeCAR-02 model attempts first of all to overcome the essential limitation of the first major model in this field, the Ursino-Lodi model [13], which considered the time course of arterial pressure as a forcing function, without explaining arterial pressure variations as a result of the hemodynamic changes induced by sudden orthostatism. In this way the new model uses the information content carried by the arterial pressure tracing in order to provide indications about the values of the autoregulation dynamics parameters.

While previous versions of pressure-flow-capacitance models [22] had significant problems explaining observations, the most recent Olufsen model [18] does an excellent job of reproducing observed arterial pressure and cerebral blood flow velocity tracings, and it is furthermore based on mechanical elements directly expressing relevant physiologic quantities (thereby making interpretation immediate). However, the Olufsen model is rather large and complex, making its use in routine clinical applications difficult.

The ESDeCAR-02 model is relatively simple, with only four differential equations (in place of 11) and only 7 free parameters to be estimated instead of the 62 of the Olufsen model.

An additional, relevant feature of the ESDeCAR-02 model is the stochastic modeling of pressure and flow variables, deriving from the representation of central venous pressure with a stochastic differential equation. This formalism is more appropriate to represent irregular variations of pressures and flows over time, produced by accidental muscle contractions, station adjustments, hormonal and neurological oscillations, respiration and swallowing, etc. The inclusion of system noise may be judged irrelevant (and hence the physiological system may be judged to be essentially deterministic), only after the volatility σ has been assessed, and it may very well be that for some subjects, but not for others, the inclusion of system noise could be indispensable to represent the evolution of the whole pressure-flow dynamics over time.

Finally, the model introduces explicitly distributed delays with respect to the four main controls (fastest to slowest: on local vasoconstriction, on heart rate, on peripheral vascular resistances, on venous capacitance), which may be presumed to depend on brain arterial perfusion pressure, consistently with generally accepted neurophysiological notions.

In the present work it has been shown that the ESDeCAR-02 model is well able to reproduce the time course of mean AP and CBFV after a sitting-to-standing change in posture, inclusive of some auto-correlated oscillations around the expected signals. In particular the model-simulated arterial pressure drops by approx 40 mmHg after standing, with a nadir at nearly 5–8 seconds post-maneuver, climbs back within another 10 to 15 seconds with an overshoot of approx 20 mmHg before stabilizing (Fig. 3, Panel A) and CBFV drops by approx 30 mL/sec after standing, with a nadir at 5 to 8 seconds post-maneuver, climbs back within another 10–15 seconds with an overshoot of approx 5 mL/sec before stabilizing (Fig. 3, Panel B). These features correspond point-by-point to what has been reported in the literature [20].

In order to explore the behaviour of the model in medically meaningful situations, a series of simulations has been produced, making single parameters vary in turn from normal to abnormal values (Fig. 4, Fig. 5). These simulations associate the variation of model parameters, such as could be physiologically postulated in relation to some dysfunction, with the time courses of pressures and flows which the model predicts would be observed in the presence of those dysfunctions.

All in all, the simulations appear very plausible: what is more important, they allow future validation or criticism of the model by comparison with experimental data sets, thereby paving the way to a future progressive refinement and calibration of CAR models.

While the present model’s structure appears reasonable, and its predictions seem consistent with observations, much work remains in fact to be done in assessing the possibility to use the model in order to discriminate between observation error and system volatility, and in determining the applicability of (simplifications of) this model to reliably estimate cerebral autoregulation in individual subjects under standard clinical conditions, which still remains an unsatisfied goal of biomathematical research in the neurosciences.

Conclusion

The proposed empirical SDDE model of cardiovascular regulation after postural change, ESDeCAR-02, is able to reproduce observed average arterial pressure and cerebral blood flow velocity profiles, including autocorrelated random oscillations around the expected time courses. The model introduces explicitly distributed delays with respect to the four main controls (fastest to slowest: on local vasoconstriction, on heart rate, on peripheral vascular resistances, on venous capacitance), which may be presumed to depend on brain arterial perfusion pressure. Much work remains to be done in assessing the possibility to use the model in order to discriminate between observation error and system volatility, and in determining the applicability of (simplifications of) this model to reliably estimate cerebral autoregulation in individual subjects.