PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of wtpaEurope PMCEurope PMC Funders GroupSubmit a Manuscript
 
Phys Rev E Stat Nonlin Soft Matter Phys. Author manuscript; available in PMC 2010 September 6.
Published in final edited form as:
Phys Rev E Stat Nonlin Soft Matter Phys. 2005 August; 72(2 Pt 1): 021905.
Published online 2005 August 19. doi:  10.1103/PhysRevE.72.021905
PMCID: PMC2933828
EMSID: UKMS31797

Nonlinear Statistical Modelling and Model Discovery for Cardiorespiratory Data

Abstract

We present a Bayesian dynamical inference method for characterizing cardiorespiratory (CR) dynamics in humans by inverse modelling from blood pressure time-series data. The technique is applicable to a broad range of stochastic dynamical models, and can be implemented without severe computational demands. A simple nonlinear dynamical model is found that describes a measured blood pressure time-series in the primary frequency band of the CR dynamics. The accuracy of the method is investigated using model-generated data with parameters close to the parameters inferred in the experiment. The connection of the inferred model to a well-known beat-to-beat model of the baroreflex is discussed.

Keywords: Bayesian inference, nonlinear time-series analysis, cardiorespiratory interaction, respiratory-sinus arrhythmia

I. INTRODUCTION

Model identification is important for both fundamental and applied research [1-8] on the human cardiovascular system (CVS). Because of the complexity of CVS dynamics, and the multiplicity of its mechanisms, it is inherently difficult or impossible to isolate and study individual response mechanisms in the intact organism [9]. In such cases mathematical models of cardiovascular control that are consistent with the experimental data can provide valuable insights [10, 11]. Altered dynamics of the cardiovascular system is associated with a range of cardiovascular diseases and with increased mortality, and it is hoped that dynamical metrics will provide a means of evaluating autonomic activity, and eventually form the basis for diagnostic tests for many conditions [12-16].

Notwithstanding that most cardiovascular controls are demonstrably nonlinear [8, 10, 13, 17-20], and that they are perturbed by stochastic inputs [8, 21, 22], assumptions of model linearity [1, 3-6, 23, 24] and/or determinism [10, 17] are often made in an attempt to facilitate progress in cardiovascular system identification. Such choices are often influenced more by the ready availability of particular statistical tools and methodologies than by biophysical or medical considerations. It is therefore highly desirable to develop reliable methods of system identification that are free of such limitations, and are capable of treating more realistic models. The latter could be used to relate difficult-to-access parameters to noninvasively-measured data [11].

Thus, while a number of numerical schemes have been proposed recently to deal with different aspects of the inverse problem using linear approximations [1, 5, 25, 26], or estimation of either the strength of some of the nonlinear terms [27, 28] or the directionality of coupling [29, 30], inverse cardiovascular problems remain difficult because of the complexity and nonlinearity of the cardiovascular interactions. The stochasticity of many dynamical inputs to the system presents an additional complication. The problem of nonlinear cardiovascular system identification has been addressed in a number of publications [7, 8, 35]. Nonlinearities generally require the use of more complex and involved numerical techniques [36-42], while the presence of dynamical noise in continuous systems can introduce systematic errors in the estimation of the model parameters [43, 44]. Analogous difficulties arise in a broad range of scientific disciplines, including problems in lasers [45] and molecular motors [46], in epidemiology [47], and in coupled matter–radiation systems in astrophysics [48]. An obstacle to progress in these fields is the lack of general methods of dynamical inference for stochastic nonlinear systems. Accordingly the methods described in this paper should be of broad interdisciplinary interest.

Following a short report published elsewhere [49], we now provide a full description of our method, applying it to the analysis of cardiorespiratory dynamics and doing so within a general nonlinear Bayesian framework for the inference of stochastic dynamical systems [44]. The basic methodology is described in Sec. II. The Bayesian inference scheme is outlined in Sec. III, and the technique is applied to the analysis of a univariate blood pressure time-series in Sec. IV. In this way, a simple nonlinear dynamical model based on coupled nonlinear oscillators [13, 50, 51] is found, able to describe time-series data in the relevant frequency range. In Sec. V, the accuracy of the method is investigated using model time series data with parameters inferred from the experimental data. The results are considered in Sec. VI, including a discussion of the connections between this model and the well-known beat-to-beat model of the baroreflex. Finally, Sec. VII summarises the conclusions that can be drawn.

II. THE METHODOLGY

Our methodological framework involves three essential steps: (i) the input data is prepared; (ii) a parameterized class of models is chosen, with the hope that one of them may describe the data; and (iii) the parameters of this model class are inferred from the time-series data.

A. Data

We analyse a particular recording of the central venous blood pressure (BP), a sample of which is shown in Fig. 1(a). A feature of this BP time-series is the presence of the two oscillatory components at frequencies approximately fr ≈ 0. 2 Hz and fc ≈ 1.7 Hz (see spectrum of this signal shown in the Fig. 1) corresponding respectively to respiratory and cardiac oscillations. The effect of nonlinear terms, including those related to the nonlinear cardiorespiratory interaction (corresponding to the side peaks) are clearly evident. We note that the relative intensities and positions of the cardiac and respiratory components vary from subject to subject, with the average frequency of respiration being around 0.3 Hz and that of the heart beat being near 1.1 Hz.

FIG. 1
Data derived from time series record 24 of the MGH/MF Waveform Database available at www.physionet.org, before and after preprocessing. (a) Original time series of the central venous blood pressure and (b) its power spectrum. (c) Respiratory component ...

In preprocessing cardiovascular data for model identification one has to bear in mind that CVS power spectra reflect a variety of complex cardiovascular interactions which give rise to peaks and other features over a broad frequency range [21, 52-54]. In order to make sense of these multi-scale phenomena, parametric modelling is usually restricted to a specific part of the power spectrum. It is clear that, in order to model the cardiorespiratory interaction, the frequency range considered must include at least the basic frequencies of cardiac and respiratory oscillations fc and fr and their combinational frequencies. Moreover, as was pointed out earlier [55, 56] locally measured blood pressure signals resemble a steady-state oscillation and the sum of the first three harmonics contains more then 70% of the total signal variance. Therefore, it is desirable that at least the 2-nd and the 3-rd harmonics (see also discussion below), besides the basic frequencies of the respiratory and cardiac oscillations, are included into the frequency range of modelling.

B. Models

When one considers modelling the cardiovascular system, one usually envisages the construction of a model based on biophysical principles, one that is capable of generating solutions that reproduce, to some degree, the measured data. This approach tackles the forward modelling problem [10, 17, 22, 57]. One may also consider the inverse modelling problem, in which models are built specifically to describe measured data [1, 3-5, 35]. Both approaches have proven useful in the context of cardiovascular research, with the forward approach providing valuable insight into the system and its causal relationships, and the inverse approach providing a useful means of intelligent monitoring of cardiovascular function in patients.

As a third alternative one may try to bridge the two approaches by building a model that accurately reproduces the experimental observations while at the same time being based on the physiological principles of circulation. In such a case the form of the mathematical model is taken from physiological principles, with its component parts corresponding to a greater or lesser degree to specific physiological mechanisms, while the values of some or all of the parameters of the mathematical model are inferred directly from the data. In such a case it is to be hoped that information with direct physiological significance, more than mere mathematical or statistical characterization, can be inferred from the data.

Many studies have been carried out to explore the physiological mechanisms underlying cardiorespiratory interactions [58-60]. Those shown to be involved are the modulation of cardiac filling pressure by respiratory movements [61], the direct respiratory modulation of parasympathetic and sympathetic neural activity in the brain stem [62], and the respiratory modulation of baroreceptor feedback control [63]. A common feature that these mechanisms is that they are nonlinear, have a dynamical (or memory) component, and are subject to exogenous fluctuations [7, 13, 50, 64-67].

A simple beat-to-beat model was introduced by DeBoer et al. [57, 68] to describe the cardio-respiratory system. The model has further been elaborated recently in [10, 11, 22]. Insight into cardio-respiratory dynamics can also be gained through inverse modelling, treating the cardiac and respiratory cycles as coupled nonlinear oscillators [13, 21, 50-52]. In this approach spectral and synchronization features [32-34] observed in the time-series data are interpreted physiologically, and related to the model parameters [13]. However, the model parameters could not be identified directly from the time-series data. Instead, they were deduced on the basis of physiological assumptions and then evaluated through extensive computer simulations [21].

The simplest model that can reproduce steady-state oscillations of the blood pressure signal at two fundamental frequencies is a system of two coupled limit cycles on a plane. The Poincaré-Bendixson theory of planar dynamical systems implies that, for a system to have limit cycle in a simply connected region, the divergence of the vector field must change sign within this region [69]. We conclude, therefore, that the simplest system that can reproduce the BP signal features considered is a planar one with limit cycles whose vector field contains polynomials of order 3. Accordingly, we model the time-series data as a system of two coupled oscillators with vector fields including nonlinearities (including nonlinearities in coupling terms) up to the 3rd order in the form

x.r=a1xr+b1yr,y.r=i=1Nαiϕi(x,y)+j=12σ1jξj,
(1)

x.c=a2xc+b2yc,y.c=i=1Nβiϕi(x,y)+j=12σ2jξj,
(2)

ξi(t)=0,ξi(t)ξj(t)=δijδ(tt).

Here the noise matrix σ mixes zero-mean white Gaussian noises ξj(t), which is related to the diffusion matrix D = σσT. The base functions are chosen in the form

ϕ={1,xr,xc,yr,yc,xr2,xc2,yr2,yc2,xryr,xcyc,xr3,xc3,xr2yr,xc2yc,xryr2,xcyc2,yr3,yc3,xrxc,xr2xc,xrxc2}.
(3)

Signals xr and xc can be directly related to the respiratory and cardiac components of the blood pressure signal shown in the Fig. 1 (c) and Fig. 1 (e) (see also Sec. IV below for a more detailed discussion). The restrictions imposed on the equations for x.r and x.c in (1) and (2) are determined mainly by the fact that we have to infer four unknown dynamical variables using univariate time-series data (see Sec. III for further details). The parametric representation (1) and (2) covers a wide range of models with limit cycles in the plane. In particular, with appropriate choices of model parameters it can describe either a van der Pol or a FitzHugh-Nagumo oscillator systems, both of which have been widely used in the context of cardiovascular modelling. Furthermore, the choice of the parametric model in the form (1), (2) allows one to relate them to physiological parameters characterizing the autonomous nervous system (note discussion below in Sec. VI). See [13] for an alternative choice and corresponding physiological reasoning, taking account of more complex interactions and considering the influence of blood flow and pressure propagation through the closed system of vessels.

C. Parameters

Following the logic of the inverse modelling approach, we must then identify the parameters M={a,b,α,β,D} of the model (1), (2) that reproduce the dynamical and spectral features of the BP signal shown in the Fig. 1. Terms representing nonlinear cardiorespiratory interactions are described by the last three base functions in (3). The correspondence of these terms to the experimentally observed combinational frequencies in the BP signal is summarized in the Fig. 2. It can be seen from the figure that the same combinational frequencies correspond to the nonlinear coupling terms in both limit cycle systems in the model. A nonlinear time-series analysis is therefor a requirement for the identification of such a model.

FIG. 2
Summary of the main harmonics of the cardiac and respiratory components observed in the BP signal. The correspondences between the nonlinear terms of the model (1), (2) and the frequencies observed in the time-series data are shown by arrows.

III. BAYESIAN INFERENCE OF STOCHASTIC NONLINEAR DYNAMICAL MODELS

Details of the Bayesian technique can be found elsewhere [44] but, for completeness, we now provide a brief description of the main steps of the algorithm.

Stochastic nonlinear dynamical models of the type (1), (2) can be expressed as a multi-dimensional nonlinear Langevin equation

x.(t)=f(x)+ε(t)=f(x)+σξ(t),
(4)

where [sm epsilon](t) is an additive stationary white, Gaussian vector noise process characterized by

ξ(t)=0,ξ(t)ξT(t)=D^δ(tt)
(5)

where D^ is a diffusion matrix.

It is assumed that the trajectory x(t) of this system is observed at sequential time instants {tk; k = 0, 1, …, K} so that the time series S={sks(tk)} thus obtained is related to the (unknown) “true” system states X={xkx(tk)} through some conditional probability density function (PDF) po(SX).

The a priori expert knowledge about the model parameters, if any, is summarized in the so-called prior PDF ppr(M). In our case we chose the prior PDF in the form of a zero-mean Gaussian distribution for the model parameters and uniform distributions for the coefficients of the diffusion matrix.

If experimental time-series data S are available they can be used to improve the estimation of the model parameters. The improved knowledge of the models parameters is summarized in the posterior conditional PDF ppost(MS), which is related to the prior PDF via Bayes’ theorem:

ppost(MS)=(SM)ppr(M)(SM)ppr(M)dM.
(6)

Here (SM), usually termed the likelihood, is the conditional PDF that relates measurements S to the dynamical model. The denominator on the right hand side of (6) is a normalizing factor. In practice, (6) can be applied iteratively using a sequence of data blocks S, S, etc. The posterior PDF computed from block S serves as the prior PDF for the next block S, etc. For a sufficiently large number of observations, ppost(MS,S,) becomes sharply peaked at a certain most probable model M=M.

The main efforts in research on stochastic nonlinear dynamical inference are focused on construction of the likelihood function, which compensates noise induced errors, and on the introduction of efficient algorithms for optimization of the likelihood function and integration of the normalization factor (cf. [36, 37, 43, 70]).

In our earlier work [44] a novel technique of nonlinear dynamical inference of stochastic systems was introduced that solves both problems. To avoid extensive numerical methods of optimization of the likelihood function and integration of the normalization factor we suggested parametrization of the vector field of (4) in the form

f(x)=U^(x)cf(x;c),
(7)

where U^(x) is an N × M matrix of suitably chosen basis functions {Unm(x); n = 1 : N, m = 1 : M}, and c is an M-dimensional coefficient vector. An important feature of (7) is that, while possibly highly nonlinear in x, f(x; c) is strictly linear in c.

The computation of the likelihood function can be cast in the form of a path integral over the random trajectories of the system [71, 72]. Using the uniform sampling scheme introduced above we can write the logarithm of the likelihood function in the following form for a sufficiently small time step h (cf. [71, 72]):

2Klog(yM)=lndetD^+hKk=0K1[v(yk)c+(y.kU^kc)TD^1(y.kU^ic)]+Nln(2πh),
(8)

which relates the dynamical variables x(t) of the system (4) to the observations s(t). Here we introduce the following notation: U^kU^(yk);y.kh1(yk+1yk); and vector v(x) with components

vm(x)=n=1NUnm(x)xn,m=1:M.

The vector elements {cm} and the matrix elements {Dnn′} together constitute a set M={c,D^} of unknown parameters to be inferred from the measurements S.

Choosing the prior PDF in the form of Gaussian distribution

ppr(M)=det(Σ^pr1)(2π)Mexp[12(ccpr)TΣ^pr1(ccpr)].
(9)

and substituting ppr(M) and the likelihood (SM) into (6) yields the posterior PDF ppost(MS)=const×exp[L(MS)], where

L(MS)Ls(c,D^)=12ρs(D^)cTws(D^)+12cTΞ^s(D^)c.
(10)

Here, use was made of the definitions

ρs(D^)=hk=0K1s.kTD^1s.k+Kln(detD^),
(11)

ws(D^)=Σ^pr1cpr+hk=0K1[U^kTD^1s.kv(sk)2],
(12)

Ξ^(D^)=Σ^pr1+hk=0K1U^kTD^1U^k.
(13)

The mean values of c and D^ in the posterior distribution give the best estimates for the model parameters for a given block of data S of length K and provide the global minimum of Ls(c,D^). We handle this optimization problem in the following way. Assume for the moment that c is known in (10). Then the posterior distribution over D^ has a mean D^post=ϴ^s(c) that provides a minimum of Ss(c,D^) with respect to D^=D^T. Its matrix elements are

ϴ^snn(c)1Kk=0K1[s.kU^(sk)c]n[s.kU^(yk)c]nT.
(14)

Alternatively, assume next that D^ is known, and note from (10) that in this case the posterior distribution over c is Gaussian. Its covariance is given by Ξ^s(D^) and the mean cpost minimizes Ls(c,D^) with respect to c

cpost=Ξ^s1(D^)ws(D^).
(15)

We repeat this two-step optimization procedure iteratively, starting from some prior values cpr and Σ^pr.

IV. ESTIMATION OF PARAMETERS OF CARDIORESPIRATORY INTERACTION FROM UNIVARIATE TIME-SERIES DATA

In order to apply algorithm (12)-(15) for the identification of the model of nonlinear cardio-respiratory dynamics (1), (2) from the univariate BP time-series of the type shown in Fig. 1(a) we have to extract time-series data corresponding to the four dynamical variables in the model. Accordingly we divide the total spectrum into a low-frequency respiratory component sr(t) and high-frequency cardiac component sc(t) as is shown in Fig. 1(c) and (e).

A discussion of the physiological relevance of this spectral separation can be found in [13, 51]). However, it is perfectly correct to consider this separation to be a mathematical ansatz. The filter parameters (see Fig. 1 caption) were chosen to preserve the 2-nd and 3-rd harmonics of these signals. The two dynamical variables of the model, xr(t) and xc(t) given by (1), (2), can be identified with the two-dimensional time-series of observations s(t) = {sr(t), sc(t)} introduced above and can be interpreted as the contribution to blood pressure from cardiac and respiratory activity. The remaining two dynamical variables y(t) = {yr(t), yc(t)} can be related to the observations {s (tk)} as follows

bnyn(tk)=sn(tk+h)sn(tk)h+ansn(tk),
(16)

where n = r, c. In this way we obtain the velocity of blood pressure changes contributed by the cardiac and respiratory components. The relation (16) is a special form of embedding. As mentioned above, it allows one to infer a wide class of dynamical models of the cardiorespiratory interactions including FitzHugh-Nagumo oscillators. The reason that we introduced restrictions on the form of the first equations in (1), (2) is now clear: it is to reduce the number of embedding parameters that must be selected to minimize the cost (10) and provide the best fit to the measured time series {s(tk)}. The corresponding simplified model of the nonlinear interaction between the cardiac and respiratory limit cycles can now be written in a form corresponding to the parametrization (7), as follows

y.=U^(s,y)c+ξ(t),
(17)

where ξ(t) is a two-dimensional Gaussian white noise with correlation matrix D^, and the matrix U^ will have the following block structure

U^=[[ϕ100ϕ1][ϕ200ϕ2][ϕB00ϕB]].
(18)

Here B = 22, so that there are 22 2 × 2 diagonal blocks formed by the basis functions given in (3), and the vector of unknown parameters c is of length M = 2B.

Finally, the model (17), (18) has to be inferred using method described in the previous section. A comparison between the time series of the inferred and actual cardiac oscillations is shown in Fig. 3. Similar results are obtained for the respiratory oscillator as shown in the Fig. 4. In each case, the level of agreement obtained is encouraging. The nonlinear coupling parameters and noise intensity of the cardiac oscillations have been estimated to have the following values: β20 = 2.2, β21 = 0.27, β22 = −8.67, and ξc2(t)=8.13. The parameters characterising coupling of respiratory oscillations to the cardiac oscillations were estimated as: α20 = 0.12, α21 = 0.048, α22 = −0.066, and D11 = 0.18. Consistent with expectations, in all experiments the parameters of the nonlinear coupling are more than an order of magnitude higher for the cardiac oscillations as compared to their values for the respiratory oscillations, reflecting the fact that respiration strongly modulates cardiac oscillations, while the effect of the cardiac oscillations on respiration is relatively weak.

FIG. 3
(Color online) (a) Time series of the cardiac oscillations xc(tn) = sc(tn) in arbitrary units (black line) obtained from measurements of the central venous blood pressure. The sampling rate was 90 Hz after resampling of the original signal. Inferred time ...
FIG. 4
(Color online) (a) Time series of the respiratory oscillations xr(tn) = sr(tn) in arbitrary units (black line) obtained from measurements of central venous blood pressure. The sampling rate was 90 Hz after resampling of the original signal. Inferred time ...

We have thus shown that this method enables one to infer simultaneously the relevant coupling strengths and noise parameters directly from a non-invasively measured time series. We view this demonstration of principle as a first step towards the practical use of this technique for cardio-respiratory modelling and its potential for clinical applications. A number of very important physiological and mathematical issues arise in relation to the application of the technique to specific problems, and we hope to address some of these in future publications. Next, however, we consider explicitly the problem of estimating the accuracy of the method, and initiate a discussion of the connections between the inferred parameters and the indices of autonomous cardiovascular control.

V. VALIDATION OF THE METHOD USING MODEL-GENERATED TIME-SERIES DATA

It is desirable to check performance of the method on synthesized time-series data obtained by numerically simulating the model (1), (2) using parameters inferred from the CVS data. To this end we consider a model-generated signal x(t) = xr (t) + xc(t) as the time-series data input s(t) for the inference. Here xr(t), xc(t) are obtained from numerical simulations of the model (1), (2) with parameters inferred from the experimental BP signal described in the previous section.

First we verify that the decomposition of the input signal s(t) into low-frequency s~r and high-frequency s~c components allows one to reconstruct the original signal. The decomposition is effected using two band-pass Butterworth filters, which were followed by application of the embedding procedure (16). In Fig. 5 we compare the velocity of the respiratory component of the original signal yr(t) with the reconstructed velocity y~r(t). The agreement is excellent. Similar results are obtained for the reconstruction of the high-frequency component. We notice in particular that the noise introduced by embedding can be neglected since it is more than an order of magnitude smaller then the dynamical noise in the signal.

FIG. 5
(Color online) (a) Velocity of the respiratory oscillations of the original signal yr(t) (green line) is compared with the signal y~r(tk) (black dashed line) obtained as a result of filtration of s(t) followed by the embedding b1y~r(tk)=(sr(tk+h) ...

Now we can apply the inference procedure described in the previous section to estimate nonlinear coupling parameters of the model from the synthesized univariate time-series data. The results of the estimation are summarized in the Table I. It can be seen from the table that the method allows one to estimate the nonlinear coupling parameter, at least to the correct order of magnitude. For some parameters the estimation accuracy is much better, but in practice the correct values are not known.

TABLE I
Absolute values of the coefficients of nonlinear cardiorespiratory interactions corresponding to the last three base functions in (3), {xr,xc, xr2xc,xrxc2}. The coefficient {αi} corresponds to respiration coupling to the cardiac rhythm. Coefficients ...

Similar results are obtained for the estimation of other parameters of the model. Using parameter values estimated from the univariate model-generated data, one can reconstruct very closely the dynamical and spectral features of the original system as shown in the Fig. 6. The largest estimation errors are for the noise intensity, as shown in the last two columns of Table I. This result can be easily understood in that filtration of the signals has the strongest effect on the noise spectrum of the system. However, the filter-induced errors are systematic and so can be corrected for based on tests with model-generated data.

FIG. 6
(Color online) (a) Synthetically generated cardiac time series xc(tn) in arbitrary units (black line) obtained from model (1), (2), compared with inferred time series of the cardiac oscillator (green line). (b) Power spectrum of the model-generated cardiac ...

The main source of error is related to the spectral decomposition of the univariate data and it is therefore systematic. To illustrate this point we use the original synthesized time-series data {xr(t), xc(t), yr(t), yc(t)} for two coupled oscillators to infer parameters of the model (1), (2). The results of inference of the coupling parameters are summarized in the Table II. It can be seen that the values of the parameters can be estimated with relative error < 10%. In particular, the relative error of estimation of the noise intensity is now < 4%. The accuracy of the estimation can be further improved by increasing the total time of observation of the system dynamics as explained in [44].

TABLE II
Absolute values of the coefficients of nonlinear cardiorespiratory interactions corresponding to the last three base functions in (3), {xr,xc, xr2xc,xrxc2}. Coefficients {αi} correspond to the respiration coupling to cardiac rhythm. Coefficients ...

These results should be compared with estimates of either relative strengths of some of the nonlinear terms [27], or the directionality of coupling [29, 30] from bivariate time-series data. It becomes clear that our algorithm provides an alternative effective approach to the analysis of cardiovascular coupling. In particular, the results of this section validate the application of the method to the measured cardiovascular data and demonstrate that it is indeed possible to estimate simultaneously the strength, directionality and the noise of nonlinear cardiorespiratory coupling form the univariate blood pressure signal.

VI. DISCUSSION

It is important to establish a relationship between the model parameters and physiological parameters of the cardiovascular system. A beat-to-beat model describing the relationships between blood pressure and respiration in simple but physiologically meaningful terms is that due to DeBoer and his collaborators [57, 68]. The DeBoer model incorporates several well-known physiological laws of the cardio-respiratory system based on static relationships. Recent extensions and modifications of this model have included [10, 11, 22]. The problem of inverse modelling was not addressed in this earlier work and it is therefore very desirable to connect it to the approach presented here.

The DeBoer model describes the beat-to-beat evolution of the state variables shown in the Fig. 7 (a): systolic pressure (S), diastolic pressure (D), RR intervals (I), and arterial decay time (T = R × C = peripheral resistance × arterial compliance). Following the brief account of the DeBoer model given in [17], and neglecting for the sake of simplicity variations in the peripheral resistance, we can write the corresponding set of difference equations as

Di=Si1exp[(23)Ii1T],
(19)

Si=Di+γIi1+C1+Asin(2πft),
(20)

Ii=GvSiτv+GβF(S,τβ)+C2,
(21)

Here C1, C2, and C3 are constants and the sigmoidal nature of the baroreceptor sensitivity is accounted for by defining an effective systolic pressure (S′) [57]

Si=S0+18arctan[(SS0)18].
(22)

The first equation (19) follows from the Windkessel model of the circulation, while the second equation (20) expresses the contractile properties of the myocardium in accordance with Starling’s law, which takes into account the mechanical effect of circulation on the BP [55, 57]. The last equation (21) includes explicitly two mechanisms of cardiovascular control defined by their respective gain (G) and delay (τ ): (i) the fast vagal control of the heart rate GvSiτv; and (ii) the slower β-sympathetic control of the heart rate GβF(S′, τβ). Here F(S′, τ) is a linear weighted sum of the form

F(S,τ)=k=MMakSiτ+k=(Siτ2+2Siτ1)+3Siτ+2Siτ+1+Siτ+2)9

Further we assume for simplicity that the pressure oscillations do not deviate far from the working point S0 in (22), i.e. S′S.

FIG. 7
(a) The BP signal in the frequency range of cardiac oscillations (black line). Systolic pressure (Sn), diastolic pressure (Dn), RR intervals (In), and arterial decay time (T = RC =??) are shown for the n-th heartbeat. (b) Comparison of the BP signal (thin ...

To establish the connection between DeBoer (19) - (21) model and the model (1), (2) introduced in this paper we note that the DeBoer model is a piece-wise approximation of the actual BP signal. In particular, it describes the BP signal as an exponential decay during 2/3 of the RR interval and a linear increase during 1/3 of the In as shown in Fig. 7 (b). We note also that this model of cardiac oscillations (2) resembles the FitzHugh-Nagumo (FHN) model of the system

{x.=(yβx),y.=αy+γy2+δy3x+C,}
(23)

where we have neglected for a moment the cardiorespiratory interaction. The approximation of the BP signal by the output of the FHN system is also shown in Fig. 7 (b). It can be seen already from a comparison between the two approximations that there is a close connection between the DeBoer model and the model of coupled oscillators considered in this paper. This can be further illustrated by noticing that for small [sm epsilon] the limit cycle in the FHN system consists of fast motion with practically constant values of y, when y jumps between negative and positive values, and slow motion, when y changes very little (see Fig. 8). Assuming a constant value of y at the top |a+| and at the bottom −|a| of the dashed curve corresponding to slow motion along the limit cycle, we can integrate the first equation in (23) to obtain

x0(t)={(Sn1+aβ)eβtaβ,for0<t<23In;(Dna+β)eβt+a+β,for23In<t<In.}

This solution closely resembles Eqs. (19) and (20) of the DeBoer model.

FIG. 8
Time evolution of the dynamical variables x (solid line) and y (dashed line) of the FHN system with parameters: [sm epsilon] = 0.01, β = −0.05, C = −0.125, α = 0.5, γ = −1, δ = 1.

It can be seen even from this simplified discussion that the parameters of the model (1), (2) found in the present paper can be related directly to the physiological parameters of the autonomous control of circulation. Furthermore, this discussion suggests that it should be possible at least in principle to bridge inverse and forward modelling and to infer parameters of the autonomous nervous control of the cardiovascular system directly from the time-series data.

We emphasize, however, that the results obtained represent only a first step in this direction. In particular, the DeBoer model itself has to be modified in various ways, including more realistic functional form of the feedback terms and specifically so as to take into account the fact that the baroreflex control is a closed loop [18, 20]. In fact it was shown [74] that a multi-compartment closed loop model of the cardiovascular responses can simulate well the experimentally observed variations in the time-series. On the other hand, this comparison suggests that the inference scheme used in this paper has to be modified in a various ways to facilitate convergence and guarantee deeper physiological meaning of the model parameters as will be discussed in more details elsewhere. It is also important to emphasize that dynamical inference of more sophisticated multi-dimensional models of the type [74], as well as coupled oscillatory models [13], can now be addressed within the framework of full Bayesian inference of the unknown dynamical variables.

VII. CONCLUSION

In the present paper we have introduced a technique for nonlinear dynamical inference of cardiovascular interactions from blood pressure time-series data. The method is applied to the simultaneous estimation of the dynamical couplings and noise strengths in a model of the nonlinear cardio-respiratory interaction. We have identified a simple nonlinear stochastic dynamical model of the cardio-respiratory interaction that describes, within the framework of inverse modelling, the time-series data in a particular frequency band. The method was validated by use of synthesized data obtained by numerically integrating the inferred model itself. We have shown that main source of error in the method is the decomposition of the blood pressure signal into two oscillatory components. We illustrate in the discussion that the dynamical model of the cardiorespiratory interaction identified in the present research can be related to the well-known beat-to-beat model of cardiovascular control introduced by DeBoer and co-workers [68]. The method developed in this paper can be used to infer the parameters of stochastic nonlinear dynamical models from observed phenomena, and is applicable across many scientific disciplines.

Footnotes

PACS numbers: 87.19.Hh, 02.50.Tt, 05.45.Tp, 05.10.Gg

References

[1] Berger RD, Saul JP, Cohen RJ. Am. J. Physiol.: Heart. Circ. Physiol. 1989;256:H142. [PubMed]
[2] Faucheux LP, Bourdieu LS, Kaplan PD, Libchaber AJ. Phys. Rev. Lett. 1995;74:1504. [PubMed]
[3] Mullen TJ, et al. Am J Physiol Heart Circ Physiol. 1997;272:H448. [PubMed]
[4] Mukkamala R, et al. Am J Physiol Regul Integr Comp Physiol. 1999;276:R905.
[5] Mukkamala R, Cohen RJ. Am J Physiol Heart Circ Physiol. 2001;281:H2714. [PubMed]
[6] Nollo G, et al. Am J Physiol Heart Circ Physiol. 2001;280:H1830. [PubMed]
[7] Chon KH, Mullen TJ, Cohen RJ. IEEE Trans. Biomed. Eng. 1996;43:530. [PubMed]
[8] Lu S, Chon KH. IEEE Trans on Sig. Proc. 2003;51:3020.
[9] Jordan D. In: Cardiovascular regulation. Jordan D, Marshall J, editors. Portland Press; Cambridge: 1995.
[10] Seidel H, Herzel H. In: Modeling the Dynamics of Biological Systems. Mosekilde E, Mouritsen OG, editors. Springer; Berlin: 1995. pp. 205–229.
[11] Seidel H, Herzel H. Physica D. 1998;115:145.
[12] Berntson GG, et al. Psychophysiology. 1997;34:623. [PubMed]
[13] Stefanovska A, Bračič M. Contemporary Physics. 1999;40:31.
[14] Malpas SC. Am. J. Physiol.: Heart. Circ. Physiol. 2002;282:H6. [PubMed]
[15] Majercak I. Bratisl Lek Listy. 2002;103:368. [PubMed]
[16] McClintock PVE, Stefanovska A. In: Fluctuations and Noise in Biological, Biophysical and Biomedical Systems II. Abbott D, Bezrukov SM, Der A, Sanchez A, editors. SPIE; Bellingham: 2004. pp. 54–68.
[17] Eyal S, Akselrod S. Meth. of Inform. in Medicine. 2000;39:118. [PubMed]
[18] Sato T, et al. Am J Physiol Heart Circ Physiol. 1999;276:H2251. [PubMed]
[19] van Leeuwen P, Bettermann H. Herzschr Elektrophys. 2000;11:127.
[20] Ringwood JV, Malpas SC. Am. J. Physiol.: Reg. Integr. and Compar. Physiol. 2001;280:R1105. [PubMed]
[21] Stefanovska A, Luchinsky DG, McClintock PVE. Physiol. Meas. 2001;22:551. [PubMed]
[22] Kotani K, Takamasu K, Ashkenazy Y, Stanley HE, Yamamoto Y. Phys. Rev. E. 2002;65:051923. [PubMed]
[23] Taylor JA, et al. Am J Physiol Heart Circ Physiol. 2001;280:H2804. [PubMed]
[24] Zou R, Chon KH. IEEE Trans. Biomed. Engin. 2004;51:219. [PubMed]
[25] Taylor JA, et al. Am J Physiol Heart Circ Physiol. 2001;280:H2804. [PubMed]
[26] Chon KH. IEEE Trans. Biomed. Engin. 2001;48:622. [PubMed]
[27] Jamšek J, Stefanovska A, McClintock PVE, Khovanov IA. Phys. Rev. E. 2003;68:016201. [PubMed]
[28] Jamšek J, Stefanovska A, McClintock PVE. Phys. Med. Biol. 2004;49:4407. [PubMed]
[29] Rosenblum MG, Cimponeriu L, Bezerianos A, Patzak A, Mrowka R. Phys. Rev. E. 2002;65:041909. [PubMed]
[30] Paluš M, Komárek V, Hrnčíř Z, Štěbrová K. Phys. Rev. E. 2001;63:046211. [PubMed]
[31] Winfree AT. Springer-Verlag. YEAR; New York: 1980.
[32] Sc[h with umlaut]afer C, Rosenblum MG, Kurths J, Abel HH. Nature. 1998;392:239. [PubMed]
[33] Janson NB, Balanov AG, Anishchenko VS, McClintock PVE. Phys. Rev. Lett. 2001;86:1749. [PubMed]
[34] Janson NB, Balanov AG, Anishchenko VS, McClintock PVE. Phys. Rev. E. 2002;65:036212/1. [PubMed]
[35] Lu S, Ju H, Chon KH. IEEE Trans. Biomed. Engin. 2001;48:1116. [PubMed]
[36] Meyer R, Christensen N. Phys. Rev. E. 2001;65:016206. [PubMed]
[37] McSharry PE, Smith LA. Physical Review Letters. 1999;83:4285.
[38] Heald JPM, Stark J. Phys. Rev. Lett. 2000;84:2366. [PubMed]
[39] Meyer R, Christensen N. Physical Review E. 2000;62:3535. [PubMed]
[40] Fullana J-M, Rossi M. Physical Review E. 2002;65:031107. [PubMed]
[41] Siegert S, Friedrich R, Peinke J. Phys. Lett. A. 1998;253:275.
[42] Siefert M, Kittel A, Friedrich R, Peinke J. Europhys. Lett. 2003;61:466.
[43] Fullana J-M, Rossi M. Phys. Rev. E. 2002;65:031107. [PubMed]
[44] Smelyanskiy VN, Timucin DA, Bandrivskiy A, Luchinsky DG. Model reconstruction of nonlinear dynamical systems driven by noise. physics/0310062 http://arxiv.org/
[45] Willemsen MB, van Exter MP, Woerdman JP. Phys. Rev. Lett. 2000;84:4337. [PubMed]
[46] Visscher K, Schnitzer MJ, Block SM. Nature. 1999;400:184. [PubMed]
[47] Earn DJD, Levin SA, Rohani P. Science. 2000;290:1360. [PubMed]
[48] Christensen-Dalsgaard J. Rev. Mod. Phys. 2002;74:1073.
[49] Smelyanskiy VN, Luchinsky DG, Stefanovska A, McClintock PVE. Phys. Rev. Lett. 2005;94:098101. [PubMed]
[50] Saul PJ, Kaplan DT, Kitney RI. Computers in Cardiology. IEEE Comput. Soc. Press; Washington: 1989. pp. 299–302. 1989.
[51] Stefanovska A, Bračič M, Strle S, Haken H. Physiol. Meas. 2001;22:535. [PubMed]
[52] Stefanovska A, Krošelj P. Open Syst. and Inf. Dyn. 1997;4:457.
[53] T. F. of the European Society of Cardiology. the North American Society of Pacing, and Electrophysiology Circulation. 1996;93:1043. [PubMed]
[54] Taylor AJ, Carr DL, Myers CW, Eckberg DL. Circulation. 1998;98:547. [PubMed]
[55] R MW. Hemodynamics. Williams & Wilkins; Baltimore: 1989.
[56] Javorka I, ans Zila M, Javorka K, Calkovska A. Physol. Res. 2002;51:227. [PubMed]
[57] DeBoer RW, Karemaker JM, Strackee J. Am. J. Physiol. 1987;253:H680. [PubMed]
[58] Taylor EW, Jordan D, Coote JH. Physiol. Rev. 1999;79:855. [PubMed]
[59] Koepchen HP. In: Mechanisms of Blood Pressure Waves. Miyakawa K, Koepchen HP, Polosa C, editors. Springer; Berlin: 1984.
[60] Melcher A. Acta Physiol. Scand. Suppl. 1976;435:1. [PubMed]
[61] Abel FL, Waldhausen JA. Am. heart J. 1969;78:266. [PubMed]
[62] Gilbey MP, Jordan D, Richter DW, Spyer KM. J. Physiol. (Lond.) 1984;365:67. [PubMed]
[63] Glass L, Mackey MC. From Clocks to Chaos. Princeton University Press; Princeton: 1988.
[64] Braun C, et al. Am J Physiol Heart Circ Physiol. 1998;275:H1577. [PubMed]
[65] Suder K, Drepper FR, Schiek M, Abel H-H. Am. J. Physiol.: Heart. Circ. Physiol. 1998;275:H1092. [PubMed]
[66] Novak V, et al. J. Appl. Physiol. 1993;74:617. [PubMed]
[67] Kanters JK, Hojgaard MV, Agner E, Holstein-Rathlou NH. Am. J. Physiol. 1997;272:R1149. [PubMed]
[68] de Boer RW, Karemaker JM, Strackee J. Psychophysiology. 1985;22:147. [PubMed]
[69] Arrowsmith DK, Place CM. Ordinary Differential Equations. Chapman and Hall; London: 1982.
[70] Meyer R, Christensen N. Phys. Rev. E. 2000;62:3535. [PubMed]
[71] Graham R. Z. Phys. B. 1977;26:281.
[72] Dykman MI. Phys. Rev. A. 1990;42:2020. [PubMed]
[73] Stan G-B, Sepulchre R. 42nd IEEE Conference on Decision and Control; PUBLISHER, Maui, Hawaii, USA. 2003.pp. 4169–4173.
[74] Heldt T, Shim EB, Kamm RD, Mark RG. J Appl Physiol. 2002;92:1239. [PubMed]
[75] Kaplan DT, Bremer CL. Physica D. 2001;160:116.