## Abstract

Cardiovascular (CV) variability as a primary vital sign carrying information about CV regulation systems is reviewed by pointing out the role of the main rhythms and the various control and functional systems involved. The high complexity of the addressed phenomena fosters a multimodal approach that relies on data analysis models and deals with the ongoing interactions of many signals at a time. The importance of closed-loop identification and causal analysis is remarked upon and basic properties, application conditions and methods are recalled. The need of further integration of CV signals relevant to peripheral and systemic haemodynamics, respiratory mechanics, neural afferent and efferent pathways is also stressed.

## 1. Introduction

Cardiovascular (CV) variability (CVV) indicates the beat-to-beat changes in CV variables, such as the heart period (HP or equivalently the heart rate) and arterial pressure (AP), which reflect the interactions of many vital regulation systems. In recordings of a few minutes, it displays oscillations related to respiration (3–4 s) and longer ones (10 s or more) related to vascular regulation mechanisms and autonomic modulations (see §2, for details).

A multimodal analysis of a set of CVV signals addresses information carried by their amplitudes, phases, time domain patterns and sensitivity to random events or imposed stimuli, in order to assess the underlying interactions and regulation mechanisms. The causal effects, gains and dynamic relationships represent the added value beyond monovariate analysis. In other words, the multimodal analysis moves from a phenomenological observation to the testing of hypotheses relevant to the origin of the observed rhythms. This approach, mostly based on a variety of linear and nonlinear approaches to CVV closed-loop interactions (Xiao *et al*. 2005), is justified by the complexity of the systems involved (Koepchen 1984) and basically relies on implicit or explicit modelling (Porta *et al*. 2006).

This review paper presents the exploitation of typical techniques for multi-channel signal processing to the field of CVV. Linear multivariate open and closed-loop identification techniques are highlighted in §3. The assessment of causality (Granger 1963) in a generalized sense and a simple extension to the nonlinear case is proposed in §4. In what follows, we will make clear that the main advantages of the multimodal approach to the analysis of CVVs are (i) the frequency and time-domain representation of the relationships among signals in terms of transfer functions (TF) and impulse responses, (ii) the separation of the feedforward actions from feedback ones in closed loops, (iii) the identification of different oscillation sources, either external to the observed loops or related to re-entry mechanisms, (iv) the assessment of causality, (v) the evaluation of indexes of global correlation and correlation along specific causal directions, (vi) the separation of each signal into components causally dependent on the other ones and into residuals, (vii) the decomposition of the power spectral densities into partial spectra dependent on each uncorrelated input (either deterministic or stochastic), and (viii) the easy customization in terms of structure of the modelled interactions, thus allowing the introduction of *a priori* relationships.

So far, most of the multimodal approaches attain to non-invasive or minimally invasive studies basically restricted to HP and AP. However, a physiological interpretation of the extracted features remains an unsolvable problem if the main haemodynamic and vascular variables are overlooked. Progress and limitations in measuring or deriving them will be discussed in §5 in order to foster future developments.

## 2. Main regulation mechanisms affecting CVV

Hales in the earliest observations of CV function through AP already noted slight beat-to-beat changes of the pulse wave, thus marking the very beginning of CVV studies (for an historical review see Koepchen 1984). Respiratory-related (second order) waves were first catalogued, while the slower (third order) Mayer waves were later recognized as an independent phenomenon related to vasomotor control (Koepchen 1984; Cohen & Taylor 2002). The spectral analysis of heart rate variability (Hyndmann *et al*. 1971; Akselrod *et al*. 1981; Pagani *et al*. 1986) opened a new window on CV regulation by simple analyses of the HP series. HP rhythms are congruent with those observed in AP variability: the respiratory sinus arrhythmia at high frequency (HF, 0.25 Hz in humans, depending on respiratory rate) and Mayer waves related oscillations at low frequency (LF, 0.1 Hz in humans, hence the 10 s rhythm). Longer recordings (up to 24 hours or more) revealed that very low-frequency oscillations display a 1/*f* power law (Goldberger & West 1987) related to the complexity of long-term regulations.

Concerning the LF and HF oscillations, the relationship with sympathetic and parasympathetic (alias vagal) tones and modulations was investigated (Pomeranz *et al*. 1985; Eckberg 1997; Malliani *et al*. 1998, Malliani 2000). General consensus was achieved that the HF power is almost exclusively sensitive to parasympathetic modulation, while the LF one, especially when expressed in normalized units, is linked to both sympathetic and parasympathetic modulations. Increased normalized LF power was shown to mark a number of manoeuvres eliciting sympathetic arousal (Malliani 2000). Sympathetic outflow to the vessels was indeed recognized as a key element in rhythm genesis (Pagani *et al*. 1997; Malpas 2002; Julien 2006). Conversely, disruption of the LF rhythm was found in physiological, e.g. intense orthostatic exercise (Perini *et al*. 1993), and pathological, e.g. severe heart failure (van de Borne *et al*. 1997), conditions eliciting a high sympathetic tone.

As to multi-channel analyses, the variability in HP and AP is the most studied pair of signals considered, thanks to the non-invasive (via plethysmography for AP) or minimally invasive availability. HP–AP closed-loop interaction involves the arterial baroreflex, thus allowing the estimation of baroreflex gain from spontaneous CVV. Spontaneous baroreflex gain estimates mimic the well-known indices of baroreflex sensitivity derived from provocative manoeuvres, e.g. pharmacological vasoconstriction (Smyth *et al*. 1969). Several baroreflex gain estimates were proposed: selection of short (e.g. three beats) sequences in time domain (Bertinieri *et al*. 1988), open-loop gain in LF band assumed to be driven by AP (Robbe *et al*. 1987; Pinna *et al*. 2002), empirical ratio of HP to AP effective root mean square (RMS) amplitude in HF and LF bands (Pagani *et al*. 1988). Causal relationships in closed loop (see §3) were parametrized to distinguish the reflex (feedback) from the HP to AP feed-forward pathway, with or without respiration as reference exogenous signal (Porta *et al*. 2000*a*). Many comparative studies of baroreflex gain indexes were performed (Pinna *et al*. 2002; Laude *et al*. 2004); however, questions relevant to their different sensitivity to confounding factors, e.g. cardiopulmonary reflexes (Lucini *et al*. 2000), are still open, fostering improvements based on multimodal approaches. Several modelling studies indicate the baroreflex loop marginal stability as a primary cause of the LF oscillations or at least as an amplifying mechanism (Guyton & Harris 1951; Kitney 1979; Boer de *et al*. 1985; Malpas 2002; Julien 2006). However, other neural and vascular mechanisms can be considered as alternative or concurrent sources as reviewed by Koepchen (1984). Also, the role of synchronization of the distributed vasomotor activity is not completely understood (Baselli *et al*. 2006).

Respiration is a primary reference in CVV studies and it is necessary to disentangle the CVV component related to breathing. Respiration perturbs CV regulation by way of cardiopulmonary reflexes to the heart (Kappagoda *et al*. 1975) and vessels (Chen *et al*. 2008), central respiratory modulations, haemodynamic effects (mainly on venous return and stroke volume) and related baroreflex responses. Accordingly, both HP and AP respiratory modulations were recounted in modelling studies (Boer de *et al*. 1985; Saul *et al*. 1989).

The imposed respiratory pattern strictly influences CVV. In normal respiration, the respiratory peak is in the HF band; however, effects are clearly in broadband, due to breathing irregularities and can extend towards the LF band (Kitney *et al*. 1985) partially contributing to LF oscillations (Baselli *et al*. 1997). Metronome respiration (experimentally imposed or spontaneous in some instances as deep sleep) elicits a narrow HF peak decreasing in power above a corner frequency (Hirsch & Bishop 1981). Slow metronome respiration close to the LF rhythm entrains it (Kitney 1979). Random metronome breathing (Saul *et al*. 1989) was also proposed as a broadband input in order to enhance system identification. Sleep apnoea or periodic breathing largely entrains CVV at frequencies below the usual LF rhythm (Jo *et al*. 2007).

Numerous studies have stressed the importance of AP control variables such as total peripheral resistance (Scher & Young 1963), cardiac output (Shoukas & Sagawa 1973) and myocardial contractility (Martin *et al*. 1969). An integration of negative feedback (Guyton & Harris 1951; Katona *et al*. 1968; Sagawa *et al*. 1973) and positive feedback (Malliani *et al*. 1986) was also considered as the expression of baroreceptive control complexity. Chemoreflex mechanisms as well as responses to muscular workload and metabolic demand (Schmidt *et al*. 2005; Iellamo *et al*. 2006) strictly interact with baroreflex regulation. Interestingly, also chemoceptive mechanisms were intensively studied by closed-loop identification methods (Ghazanshahi & Khoo 1997; Topor *et al*. 2004).

## 3. Closed-loop identification of CV interactions

The ongoing activity of CV regulation working in closed loop can be analysed from a comparison of CVV signals by means of proper identification methods. This approach avoids invasive manoeuvres necessary to either open the regulatory loops or highlight specific reflexes by overwhelming stimuli that may cause unwanted reset or adaptation. Furthermore, the observation time window is not limited by re-entry mechanisms that occur when the feedback shows up its action. Here, some basic properties of the usual open-loop identification are recalled in order to stress the working hypotheses necessary for a closed-loop extension and to discuss the pitfalls relevant to improper application.

### (a) Open-loop condition

System identification (Ljung 1999) provides a black-box description of input–output relationships from data series. Considering the block diagram of figure 1*a*, an estimate of the unknown TF *G*(*z*) from input *u* to output *y* is derived, by choosing a parametric family of models *G*(*z*|** ϑ**) and computing the parameter vector producing the best fit of

*N*output samples given the input. Here, and in what follows, the notation of linear time invariant discrete time systems is maintained, by which TFs are

*z*-transforms of impulse responses and

*z*

^{−1}is the one-step lag operator. Hence, the TF is considered an estimate of

*G*(

*z*). Usually, identification is performed by minimizing the differences between data and predictions generated according to a prediction error model, often by the least-squares (LS) procedure. While estimating the deterministic input–output relationship, knowledge is also gained about the stochastic part of the model

*H*(

*z*). This is the spectral factor that describes the dynamics (memory and oscillations) of residual

*ν*by colouring the white noise

*ϵ*. Better predictions are obtained if the stochastic part is described by a correct estimate leading to the white prediction error series.

LS identification based on the prediction error model shows several noticeable properties.

*Consistency.* If the model family includes the true model, i.e. if a exists such that *G*(*z*)=*G*(*z*|** ϑ**°) and

*H*(

*z*)=

*H*(

*z*|

**°); than for**

*ϑ**N→∞*.

*Accuracy.* Estimate error covariance at a given angular frequency *ω* (i.e. for *z*=e^{jω}) is proportional to the ratio between parameter number *n* and data number *N* and to the ratio between the spectrum of residuals and that of the input(3.1)the last item states that the system must be perturbed by the input at each frequency, i.e. the input must be persistently exciting, such as a white noise or at least a broadband signal.

*Lack of bias.* If the stochastic part is forced to an incorrect shape *H*(*z*|** ϑ**)=

*H*

^{*}(

*z*) (typically due to a poor parametrization), then the parameters of the deterministic TF tend to an estimate that minimizes the following functional,(3.2)clearly, if the parametrization includes the true model of

*G*, the functional minimum is there, i.e. consistency holds even if the stochastic part is not estimated and errors are not whitened.

Note that in modern identification theory (Ljung 1999), the above properties are accredited to the TF, rather than to the parameter vector, thus rendering them more readily interpretable.

Popular software packages such as the System Identification Toolbox of Matlab permit easy estimation parameters of all parametric structures that are particular cases of the general parametric model(3.3)where *B* is a polynomial in *z*^{−1}, while *A*,*C*,*D*,*F* are monic polynomials in *z*^{−1}, with unit zero order coefficient. Their order determines the number of past samples used in regressions (alias model order) to be optimized by means of well known model order criteria, e.g. final prediction error, Akaike's information criterion or minimum description length.

### (b) Closed-loop extension

Above open-loop results require that the deterministic input *u* is not correlated with the stochastic one *ν*; i.e. no feedback may exist from the output to the input. This is a severe condition rarely fulfilled in intact physiological systems devoted to maintain homeostasis by means of closed loops (Guyton & Hall 1996; Khoo 2000).

Nonetheless, conditions, specific for closed loops, permit the use of open-loop identification algorithms. The basic problem structure for two signals is exemplified in figure 1*b*: two deterministic blocks, *G*_{12} (feed-forward) and *G*_{21} (feedback) fed by two coloured and uncorrelated processes, *ν*_{1} and *ν*_{2}. (The usual minus sign of negative feedbacks is not evidenced here, for symmetry reasons; it is included in one of the loop blocks, when required.)

A central role in the loop description and in its identifiability is played by the sensitivity function(3.4)which describes the amplification (or attenuation) of any input immediately after the entering node, due to the loop recurrence. Note that feedback working well, with high loop gain—*G*_{12}(*z*) *G*_{21}(*z*), is stiff versus external noise with a low sensitivity and poor identifiability. Fortunately, physiological loops are rarely stiff and their continuous reaction to perturbations (i.e. innovation) can be exploited.

Depending on the measured data and on the *a priori* knowledge, either an indirect or a direct identification approach can be adopted (Ljung 1999).

### (c) Indirect approach

As shown in figure 1*c*, the indirect approach permits use of the open-loop identification by looking over the closed loop from an external input *u*=*ν*_{2} to the output *y*_{1} (*step-1*) and next disentangling the loop by introducing the *a priori* knowledge of a part of it (*step-2*).

*Step 1* is based on the necessary condition that an exogenous input *u* is measured immediately outside the loop and enters it at a specific point. Hence, the closed loop can be encapsulated in TF *G*(*z*), thus satisfying the open-loop condition. Given the *a priori* knowledge of either *G*_{12} or *G*_{21}, *step 2* computes the unknown part by solving . For example in problems of automatic control tuning, *u* can be obtained by perturbing the reference signal, while the loop is solved by exploiting the knowledge of the controller.

These conditions are rarely fulfilled in physiology; nonetheless, interesting applications were proposed. Saul *et al*. (1989) used randomly paced respiration as non-invasive broadband input to probe the sinus node and the baroreflex system. Though the authors did not explicitly present their method as an indirect closed-loop identification, nonetheless their approach parallels *step 1* above, since it creates a condition in which the measured input is taken outside the loop. The second step was overlooked thanks to the important perturbation imposed and to selective parasympathetic and sympathetic activation/blockades.

Mukkamala *et al*. (2006*a*) explicitly applied an indirect approach for the identification of the vasomotor cardiopulmonary reflex and found more consistent results than by a direct approach. A deterministic input, representing a surrogate of right atrial pressure, was obtained by means of a comparison between cardiac output and stroke volume changes. This signal was exogenous to the arterial control of AP; hence, the relation from right atrial pressure to AP satisfied the open-loop condition (*step-1*). Finally, *a priori* knowledge of vascular Windkessel dynamics and of the arterial baroreflex was used to disentangle the pathway from right atrial pressure to AP from the AP to AP re-entry (*step-2*).

### (d) Direct approach and its parametrization

Often, only data from the inside of the loop are measurable (*y*_{1} and *y*_{2} in figure 1*d*) and none of *G*_{12} and *G*_{21} is known. A relevant example is the HP–AP pair interacting in the baroreflex loop. In this case, it is necessary to use a direct method, which virtually opens the loop while identifying one of the pathways (e.g. feedforward *G*_{12}) simply by ignoring the opposite one (e.g. feedback *G*_{21}, accordingly dashed in figure 1*d*). Clearly, the open-loop uncorrelation between *y*_{2}, playing the role of input and *ν*_{1}, playing the role of output noise, does not hold. Hence, different closed-loop conditions must be verified: (i) causality of feedforward and feedback and (ii) uncorrelation of stochastic inputs *ν*_{1} and *ν*_{2}. Causality is a reasonable hypothesis when physical events (e.g. reflexes and control actions) are considered and parametric modelling permits this feature to be described. Strict causality is necessary for at least one of the blocks in each loop (e.g. HP can display haemodynamic effects on the next beat), to avoid algebraic loops. Non-strict causality can be admitted in the other blocks, where zero-lag parameters describe effects within a sampling period (e.g. AP can affect the concurrent HP by way of fast baroreflex). Uncorrelation of inputs is by far a more severe and problematic hypothesis. *G*_{12} and *G*_{21} over-parametrization can force *ν*_{1} and *ν*_{2} uncorrelation confounding external correlations with loop effects. Nothing but experimental validation or physiological knowledge can solve this ambiguity.

Further conditions concerning correct parametrization of both the causal and the stochastic parts and sufficient innovation govern identification properties and impose more severe constraints compared with the open-loop case.

*Consistency*—This fundamental property holds if the true model of both the deterministic *G*_{12} and the stochastic *H*_{1} part are correctly modelled.

*Accuracy*—If the spectrum of *y*_{2} is divided into the sum of the two uncorrelated components responding to *ν*_{1} and *ν*_{2},(3.5)then, different from (3.1), only the latter contributes to parameter convergence, while feedback noise is useless(3.6)Since , note that a low sensitivity amplifies identification errors.

*Bias*—An incorrect identification of the noise spectral factor forcing *H*_{1}(*z*|** ϑ**)=

*H*

_{1}

^{*}(

*z*) causes a bias of , which converges to , with(3.7)where

*λ*

_{1}is the variance of

*ϵ*

_{1}. Note that the bias is proportional to the feedback noise , hence null in the open-loop case, as previously stated.

In general, equation (3.7) shows that direct identification tries to attribute noise dynamics to the loop, if an insufficient parametrization of the stochastic part does not allow it to locate externally. In other words, noise oscillations may be misinterpreted as loop resonances. The consequent confusion of oscillation sources is particularly severe in CVV multimodal analysis, especially when assessing the origin of CVV rhythms (Koepchen 1984; Malpas 2002; Julien 2006). Thus, specific model structures are required (Baselli *et al*. 1994).

Moreover, Granger's causality analysis is misleading. Let us consider the joint process(3.8)where is the matrix spectral factor that indicates the absence of either causal direction when upper or lower triangular, i.e. with *G*_{21}=0 (no feedback) or *G*_{12}=0 (no feedforward).

E.g. if *G*_{12}=0, if the output noise *ν*_{1} and the feedback *G*_{21} are not null, and if fits a wrong *H*_{1}^{*}, then ; since all the factors in (3.7) are strictly positive. So, a false feedforward causality is detected due to the bias. In addition to a correct modelling of *H*_{1}, hypothesis testing by surrogate series (see §4) is always recommended, in order to separate possible bias effects. Interestingly, in the further hypothesis of low input noise, *ν*_{2}→0, then (with *S*=1, for *G*_{12}=0)) and (3.7) tends to(3.9)thus showing that the inverse of the feedback is misinterpreted for the feedforward and this error is weighted by the per cent error of the stochastic part.

Two kinds of parametrizations were practically proposed for the direct identification of CVV loops.

(i) Bivariate AR, formed by a feedforward and a feedback ARX model, as in figure 2

*a*. Bivariate AR is derived from the two ARX by setting*A*_{12}(*z*)=−*B*_{12}(*z*) and*A*_{21}(*z*)=−*B*_{21}(*z*).(ii) Bivariate dynamic adjustment, formed by a feedforward and a feedback ARXAR model, as in figure 2

*b*.

In ARX identification LS parameters are obtained in closed form. Adaptive (time variant) extensions are available (Mainardi *et al*. 1997) as well. However, poles of the feedforward (feedback) deterministic part are constrained to be the same as of the stochastic part *H*_{1} (*H*_{2}). Equivalently, the joint process displays only the poles of its determinant. Wold's theorem still assures identifiability asymptotically with the model order. This solution can be applicable for a causal analysis or for the extraction of block gains if properly tested; but it contrasts with a classification of oscillation sources by means of a limited number of poles.

An ARXAR identification requires iterative algorithms (e.g. generalized LS); however, actual software packages do not require the user to deal with algorithmic details. Unlike ARX, ARXAR allows the residuals (*ν*_{1} and *ν*_{2}) with poles different from those created by the feedback loops to be described thus disentangling the respective spectral components (Baselli *et al*. 1997). This structure was proposed to analyse the interactions between HP and AP (namely, SAP), including a HP–AP loop and an AP–AP one (Baselli *et al*. 1988); respiration was also considered as further exogenous input affecting both HP and AP. The hypothesis of a further short-term (1 or 2 heart beats) HP–HP loop, describing sinus node dynamics was later successfully tested (Porta *et al*. 2003).

Identification results can be displayed as TFs, partial spectra, gains and other significant parameters. Also, time domain plots of impulse responses can be readily interpreted if regularized by proper kernel expansions (Xiao *et al*. 2004). A trivial plot of the partial signals, derived by filtering the measured signals by the identified TFs, displays information relevant to the identified linear features and also enhances transient deviations from the model, due to nonlinearities or non-stationarities (Baselli *et al*. 2001).

## 4. Causal CV interactions

Traditional multi-channel spectral and cross-spectral analysis (Kay & Marple 1981) enhances correlation in frequency bands but loses causal directions. Conversely, methods based on Granger's (1963) causality do accept a directional interaction from a signal *u* to a signal *y* if the additional knowledge of *u* past values improves the prediction of *y*. These approaches exploit primarily cross-predictability (Schiff *et al*. 1996; Le Van Quyen *et al*. 1998; Wiesenfeldt *et al*. 2001; Feldmann & Bhattacharya 2004; Faes *et al*. 2006) or transfer entropy (Porta *et al*. 1999; Palus *et al*. 2001; Kaiser & Schreiber 2002) or partial or causal coherence (Baccalà & Sameshima 2001; Kaminski *et al*. 2001; Porta *et al*. 2002).

Examples of the application of recent methods (Porta *et al*. 2000*b*, 2002; Nollo *et al*. 2002, 2005) over short temporal sequences (i.e. with *N*=300 heart beats) will be given. Beat-to-beat variability series of HP measured as the temporal distance between two consecutive R peaks (RR) and SAP at rest and during standing is shown in figure 3*a*,*b*,*d*,*e*. In the following, RR and SAP, normalized in order to have zero mean and unit variance, will be indicated as and *sap*(*i*), .

### (a) Generalized causal analysis

Let us define as and the ordered sequences of *n rr* and *m sap* samples. A generalized (linear or nonlinear) relationship between *rr* and *sap* along the causal direction from *sap* to *rr* can be denoted as(4.1)where describes the dependence of *rr*(*i*) on *n* past values of *rr* and on the current and *m*-1 past values of *sap* and *ϵ*_{rr} is a white noise independent of . Under this hypothesis the best prediction of *rr* (i) is(4.2)where is the best estimate of .

If is linear, equation (4.1) reduces to the ARX model (§3*d*) and reference can be made to the feedback part of figure 2*a*: *y*_{2}=*rr* (single output), *y*_{1}=*sap* (exogenous input), *ϵ*_{2}=*ϵ*_{rr} (white noise on *rr*), *A*_{22}(*z*)=*A*_{rr}=(*z*) (stochastic part) and (deterministic part). If is nonlinear it can be approximated globally by NARX models (Chen & Billings 1989) or Volterra series (Marmarelis 1993) or locally by *k*-nearest-neighbour strategies (Abarbanel *et al*. 1994). Both in the linear and the nonlinear cases, the agreement between the predicted and original series can be evaluated as the squared correlation coefficient between *rr*(*i*) and as(4.3)where *k* is the maximum between *n* and *m* (Porta *et al*. 2007*a*). is bounded between 1 (perfect correlation, i.e. *rr* is perfectly predictable given *sap*) and 0 (full uncorrelation, i.e. *rr* cannot be predicted given *sap*).

The causal direction (here ) can be assessed according to three approaches. The first and the second ones test a significant difference from 0 of either (i.e. past samples of *rr* are not taken into account for setting the best prediction; Schiff *et al*. 1996; Le Van Quyen *et al*. 1998) or (Wiesenfeldt *et al*. 2001; Faes *et al*. 2006), respectively. The third approach is based on the concept of predictability improvement: it compares (i.e. past samples of *sap* are not taken into account) to (Feldmann & Bhattacharya 2004). If the latter is significantly larger than the former, then the two series are coupled with *sap*→*rr*.

All considerations expressed above for *sap*→*rr* do symmetrically apply to *sap*←*rr*. If *sap*→*rr* coexists with *sap*←*rr* a closed-loop interaction is demonstrated.

### (b) Causal analysis based on cross-conditional entropy

Instead of using (4.3), predictability of *rr* given previous *sap* values can be evaluated using cross-conditional entropy (CCE; Porta *et al*. 1999, 2000*b*; Nollo *et al*. 2002). measures the information carried by *rr* when (*m*−1) *sap* values are given as(4.4)where is the probability of and is the conditional probability of *rr* given the present and *m*−2 past *sap* samples. In (4.4), the first sum is extended to all different 's found in the set of data for and the second sum to all different *rr*'s following a specific *sap* sequence . Probabilities were evaluated on six bins after uniform quantization. If is significantly lower than (i.e. *rr* entropy), then *sap* knowledge lowers *rr* uncertainty and *sap*→*rr* holds. is strongly biased towards 0 due to poor statistics as the increase of *m* is accompanied by an exponential growth of possible patterns. Bias correction strategies have been accordingly proposed (Porta *et al*. 2007*b*). Corrected (Porta *et al*. 1999) adds a corrective term proportional to patterns appearing only once, which give a null contribution to and artefactually lower it. It was shown that corrected (i) remains high for all *m* values when *rr* and *sap* are uncoupled, (ii) decreases to 0 if *rr* deterministically depends on *sap and* (iii) exhibits a minimum if the knowledge of *sap* is helpful to reduce the uncertainty of *rr* prediction.

Figure 3*c*,*f* shows corrected calculated over *rr* and *sap* series (bold line) at rest and during standing, respectively. Corrected exhibits a well-defined minimum in both conditions. The lower minimum during standing indicates that the increased sympathetic activity and augmented LF rhythm enforce causality from *sap* to *rr* (Porta *et al*. 2000*b*; Nollo *et al*. 2002). In figure 3*c*,*f*, corrected *CCE*s derived from the original data (bold line) are compared with those relevant to surrogate series generated according to the null hypothesis of uncoupling. Surrogate pairs were obtained by applying two-independent phase-randomization procedures preserving both the power spectrum and the sample distribution of the original *rr* and *sap* series (Schreiber & Schmitz 1996; Palus 1997). Since the minima of the corrected CCEs derived from the original series are significantly lower than the surrogate range, the coupling is significant.

### (c) Frequency domain causal analysis by causal coherence

The above methods provide global indicators of causality; by contrast, frequency domain linear methods permit specific CVV rhythms to be dealt with. Squared coherence function evaluates linear correlation (between 0 and 1) in the frequency domain as(4.5)

The cross spectrum, and the power spectra of *rr* () and *sap* () are derived from a bivariate AR process (figure 2, with , , , )(4.6)yielding(4.7)where the TF matrix of the joint process , is its Hermitian, and is the variance matrix of the white uncorrelated noises and *ϵ*_{rr}.

In figure 4*a*,*d,* the squared coherence at rest and during standing are, respectively, shown. Again curves derived from the original series (bold lines) are compared with those obtained from surrogates (light lines). The presence of coherence peaks at LF and HF frequencies above the distribution of surrogates demonstrates significant linear coupling. Note that also the surrogates display artefactual coherence peaks due to biases relevant to power peaks; hence, a comparison with surrogates is mandatory in order to prevent false positives (Porta *et al*. 2002).

does not assess causal direction, since it relies both on feedforward and feedback . Fixing to zero either block can therefore provide a causal coherence relevant to the opposite pathway alone (Porta *et al*. 2002). Accordingly, *rr*→*sap* is tested by(4.8)and *sap*→*rr* by(4.9)

At rest both (figure 4*b*) and (figure 4*c*) exhibit significant peaks at LF (); conversely, the HF peak is significant (i.e. above surrogates) only in the latter direction (*sap*→*rr*), in relation to a prevalent baroreceptive feedback and vagal drive. During standing (figure 4*e*,*f*) both in LF and HF bands is detected. This behaviour was also observed during sympathetic activation produced by head-up tilt (Nollo *et al.* 2005).

In conclusion, both linear parametric and entropic methods allow the assessments of Granger's causality by means of the significance of either a regression or a probabilistic dependence. The latter approach comprises the former and extends the analysis to nonlinear features. The lack of modelling constraints of entropic methods, however, penalizes them by an explosion of patterns, as more and more past values are considered, and limits the length of patterns that can be explored.

Linear causal analysis is strictly connected to direct closed-loop identification; although, emphasis is given to the relevance of causal effects rather than to the features of the causal blocks and of the oscillation sources. Identifiability conditions are the same as for closed-loop identification (see §3*d*), but comparison with surrogates can effectively separate estimation biases.

## 5. Towards the integration of haemodynamic variables and CVV modelling

Most CVV studies are for invasiveness reasons limited to AP, lacking cardiac output and flow measures, which are necessary to assess the general status of vasomotor control, peripheral resistances and the metabolic demand. For instance, even a gross drop in cardiac output due to haemorrhage is masked by peripheral resistance compensation if AP regulation only is observed (Guyton & Harris 1951).

While AP dynamics can be measured in various ways, either invasively or non-invasively (see Allen 2007 for a review), the access to cardiac output beat-to-beat variability is strictly invasive. Non-invasive or moderately invasive gold standards such as the Fick technique, thermodilution, dye dilution and acetylene rebreathing (Hoeper *et al*. 1999) are limited to averages over a time duration encompassing several beats. Thus, over the last few years, innovative studies have aimed at reconstructing cardiac output changes from peripheral measurements of AP. Most of them are based on lumped Windkessel modelling or on TF estimate (Lu & Mukkamala 2006; Mukkamala *et al*. 2006*b*; Swamy *et al*. 2007). These techniques are now sufficiently robust to be exploited in a multimodal approach including cardiac output into the set of signals for CV monitoring.

As to peripheral vasomotion, laser Doppler flowmetry allows a non-invasive, continuous measure of microcirculatory blood flow and the direct assessment of peripheral vasomotion. Stefanovska *et al*. (1999) showed that the laser Doppler signal contains significant information. The following five characteristic bands were labelled: (i) heart pulse at heart rate (0.4–1.6 Hz), (ii) a respiratory band (0.15–0.4 Hz), (iii) a band (0.06–0.15 Hz) largely overlapped by the LF one and designated as myogenic, (iv) a very LF band (0.02–0.06 Hz) hypothesized to reflect metabolic regulation and (v) an ultra LF band (0.009–0.02 Hz) related to endothelial activity and nitric oxide release. Although several studies addressed the link between vasomotion and systemic CVV (Bernardi *et al*. 1997; Stauss *et al*. 1998), the possibility of integrating signals describing peripheral vasomotion with other CVVs remains largely unexplored. The study of the CV system might take advantage from the introduction of the information content related to vascular regulation activity. A multimodal approach including vasomotion into the set of CVV signals might allow the assessment of the origin of the slow oscillations clearly present in diastolic AP series and, sometimes, evident in SAP series and the study of synchronization among peripheral districts responsible for AP changes at systemic level via synchronous vasoconstrictions (Baselli *et al*. 2006).

## 6. Conclusion

Multimodal approaches might provide the framework to integrate physiological and clinical knowledge, biomedical instrumentation, non-invasive techniques, modelling methodologies and multivariate data analysis, thus allowing a more insightful representation and description of a highly integrated and complex system such as the CV one.

## Acknowledgements

Partially supported by ASI, DCMC project and by EU Heart Cycle project.

## Footnotes

One contribution of 13 to a Theme Issue ‘Signal processing in vital rhythms and signs’.

- © 2008 The Royal Society