Home | About | Journals | Submit | Contact Us | Français |

**|**Europe PMC Author Manuscripts**|**PMC2933828

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. THE METHODOLGY
- III. BAYESIAN INFERENCE OF STOCHASTIC NONLINEAR DYNAMICAL MODELS
- IV. ESTIMATION OF PARAMETERS OF CARDIORESPIRATORY INTERACTION FROM UNIVARIATE TIME-SERIES DATA
- V. VALIDATION OF THE METHOD USING MODEL-GENERATED TIME-SERIES DATA
- VI. DISCUSSION
- VII. CONCLUSION
- References

Authors

Related links

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.021905PMCID: PMC2933828

EMSID: UKMS31797

D.G. Luchinsky,^{1,}^{2,}^{3} M.M. Millonas,^{2} V.N. Smelyanskiy,^{2} A. Pershakova,^{3} A. Stefanovska,^{4} and P.V.E. McClintock^{3}

See other articles in PMC that cite the published article.

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.

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.

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.

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 *f _{r}* ≈ 0. 2 Hz and

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 *f _{c}* and

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

$${\stackrel{.}{x}}_{r}={a}_{1}{x}_{r}+{b}_{1}{y}_{r},\phantom{\rule{1em}{0ex}}{\stackrel{.}{y}}_{r}=\sum _{i=1}^{N}{\alpha}_{i}{\varphi}_{i}(\mathbf{x},\mathbf{y})+\sum _{j=1}^{2}{\sigma}_{1j}{\xi}_{j},$$

(1)

$${\stackrel{.}{x}}_{c}={a}_{2}{x}_{c}+{b}_{2}{y}_{c},\phantom{\rule{1em}{0ex}}{\stackrel{.}{y}}_{c}=\sum _{i=1}^{N}{\beta}_{i}{\varphi}_{i}(\mathbf{x},\mathbf{y})+\sum _{j=1}^{2}{\sigma}_{2j}{\xi}_{j},$$

(2)

$$\langle {\xi}_{i}\left(t\right)\rangle =0,\phantom{\rule{1em}{0ex}}\langle {\xi}_{i}\left(t\right){\xi}_{j}(t\prime )\rangle ={\delta}_{ij}\delta (t-t\prime ).$$

Here the noise matrix *σ* mixes zero-mean white Gaussian noises *ξ _{j}*(

$$\varphi =\{1,{x}_{r},{x}_{c},{y}_{r},{y}_{c},{x}_{r}^{2},{x}_{c}^{2},{y}_{r}^{2},{y}_{c}^{2},{x}_{r}{y}_{r},{x}_{c}{y}_{c},{x}_{r}^{3},{x}_{c}^{3},{x}_{r}^{2}{y}_{r},{x}_{c}^{2}{y}_{c},{x}_{r}{y}_{r}^{2},{x}_{c}{y}_{c}^{2},{y}_{r}^{3},{y}_{c}^{3},{x}_{r}{x}_{c},{x}_{r}^{2}{x}_{c},{x}_{r}{x}_{c}^{2}\}.$$

(3)

Signals *x _{r}* and

Following the logic of the inverse modelling approach, we must then identify the parameters $\mathcal{M}=\{\mathbf{a},\mathbf{b},\alpha ,\beta ,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.

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

$$\stackrel{.}{\mathbf{x}}\left(t\right)=\mathbf{f}\left(\mathbf{x}\right)+\epsilon \left(t\right)=\mathbf{f}\left(\mathbf{x}\right)+\sigma \xi \left(t\right),$$

(4)

where (*t*) is an additive stationary white, Gaussian vector noise process characterized by

$$\langle \xi \left(t\right)\rangle =0,\phantom{\rule{1em}{0ex}}\langle \xi \left(t\right){\xi}^{T}(t\prime )\rangle =\widehat{\mathbf{D}}\delta (t-t\prime )$$

(5)

where $\widehat{\mathbf{D}}$ is a diffusion matrix.

It is assumed that the trajectory *x*(*t*) of this system is observed at sequential time instants {*t _{k}*;

The *a priori* expert knowledge about the model parameters, if any, is summarized in the so-called *prior* PDF ${p}_{\mathrm{pr}}\left(\mathcal{M}\right)$. 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 $\mathcal{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 ${p}_{\text{post}}(\mathcal{M}\mid \mathcal{S})$, which is related to the prior PDF via Bayes’ theorem:

$${p}_{\text{post}}(\mathcal{M}\mid \mathcal{S})=\frac{\ell (\mathcal{S}\mid \mathcal{M}){p}_{\mathrm{pr}}\left(\mathcal{M}\right)}{\int \ell (\mathcal{S}\mid \mathcal{M}){p}_{\mathrm{pr}}\left(\mathcal{M}\right)\mathrm{d}\mathcal{M}}.$$

(6)

Here $\ell (\mathcal{S}\mid \mathcal{M})$, usually termed the *likelihood*, is the conditional PDF that relates measurements $\mathcal{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 $\mathcal{S}$, $\mathcal{S}\prime $, etc. The posterior PDF computed from block $\mathcal{S}$ serves as the prior PDF for the next block $\mathcal{S}\prime $, etc. For a sufficiently large number of observations, ${p}_{\text{post}}(\mathcal{M}\mid \mathcal{S},\mathcal{S}\prime ,\dots )$ becomes sharply peaked at a certain most probable model $\mathcal{M}={\mathcal{M}}^{\ast}$.

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

$$\mathbf{f}\left(\mathbf{x}\right)=\widehat{\mathbf{U}}\left(\mathbf{x}\right)\mathbf{c}\equiv \mathbf{f}(\mathbf{x};\mathbf{c}),$$

(7)

where $\widehat{\mathbf{U}}\left(\mathrm{x}\right)$ is an *N* × *M* matrix of suitably chosen basis functions {*U _{nm}*(

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]):

$$-\frac{2}{K}\mathrm{log}\ell (\mathbf{y}\mid \mathcal{M})=\mathrm{ln}\phantom{\rule{thinmathspace}{0ex}}\mathrm{det}\phantom{\rule{thinmathspace}{0ex}}\widehat{\mathbf{D}}+\frac{h}{K}\sum _{k=0}^{K-1}[\mathbf{v}\left({\mathbf{y}}_{k}\right)\mathbf{c}+{({\stackrel{.}{\mathbf{y}}}_{k}-{\widehat{\mathbf{U}}}_{k}\mathbf{c})}^{T}{\widehat{\mathbf{D}}}^{-1}({\stackrel{.}{\mathbf{y}}}_{k}-{\widehat{\mathbf{U}}}_{i}\mathbf{c})]+N\phantom{\rule{thinmathspace}{0ex}}\mathrm{ln}\left(2\pi h\right),$$

(8)

which relates the dynamical variables **x**(*t*) of the system (4) to the observations **s**(*t*). Here we introduce the following notation: ${\widehat{\mathbf{U}}}_{k}\equiv \widehat{\mathbf{U}}\left({\mathbf{y}}_{k}\right)$;${\stackrel{.}{\mathbf{y}}}_{k}\equiv {h}^{-1}({\mathbf{y}}_{k+1}-{\mathbf{y}}_{k})$; and vector **v**(**x**) with components

$${\mathrm{v}}_{m}\left(\mathbf{x}\right)=\sum _{n=1}^{N}\frac{\partial {U}_{nm}\left(\mathbf{x}\right)}{\partial {x}_{n}},\phantom{\rule{1em}{0ex}}m=1:M.$$

The vector elements {*c _{m}*} and the matrix elements {

Choosing the prior PDF in the form of Gaussian distribution

$${p}_{\mathrm{pr}}\left(\mathcal{M}\right)=\sqrt{\frac{\mathrm{det}\left({\widehat{\Sigma}}_{\mathrm{pr}}^{-1}\right)}{{\left(2\pi \right)}^{M}}}\mathrm{exp}\left[\frac{-1}{2}{(\mathbf{c}-{\mathbf{c}}_{\mathrm{pr}})}^{T}{\widehat{\Sigma}}_{\mathrm{pr}}^{-1}(\mathbf{c}-{\mathbf{c}}_{\mathrm{pr}})\right].$$

(9)

and substituting ${p}_{\mathrm{pr}}\left(\mathcal{M}\right)$ and the likelihood $\ell (\mathcal{S}\mid \mathcal{M})$ into (6) yields the posterior PDF ${p}_{\text{post}}(\mathcal{M}\mid \mathcal{S})=\mathrm{const}\times \mathrm{exp}[-L(\mathcal{M}\mid \mathcal{S})]$, where

$$L(\mathcal{M}\mid \mathcal{S})\equiv {L}_{s}(\mathbf{c},\widehat{\mathbf{D}})=\frac{1}{2}{\rho}_{s}\left(\widehat{\mathbf{D}}\right)-{\mathbf{c}}^{T}{\mathbf{w}}_{\mathrm{s}}\left(\widehat{D}\right)+\frac{1}{2}{\mathbf{c}}^{T}{\widehat{\Xi}}_{\mathrm{s}}\left(\widehat{\mathbf{D}}\right)\mathbf{c}.$$

(10)

Here, use was made of the definitions

$${\rho}_{\mathrm{s}}\left(\widehat{\mathbf{D}}\right)=h\sum _{k=0}^{K-1}{\stackrel{.}{\mathbf{s}}}_{k}^{T}{\widehat{\mathbf{D}}}^{-1}{\stackrel{.}{\mathbf{s}}}_{k}+K\phantom{\rule{thinmathspace}{0ex}}\mathrm{ln}\left(\mathrm{det}\phantom{\rule{thinmathspace}{0ex}}\widehat{\mathbf{D}}\right),$$

(11)

$${\mathbf{w}}_{\mathrm{s}}\left(\widehat{\mathbf{D}}\right)={\widehat{\Sigma}}_{\mathrm{pr}}^{-1}{\mathbf{c}}_{\mathrm{pr}}+h\sum _{k=0}^{K-1}\left[{\widehat{\mathbf{U}}}_{k}^{T}{\widehat{\mathbf{D}}}^{-1}{\stackrel{.}{\mathbf{s}}}_{k}-\frac{\mathbf{v}\left({\mathbf{s}}_{k}\right)}{2}\right],$$

(12)

$$\widehat{\Xi}\left(\widehat{\mathbf{D}}\right)={\widehat{\Sigma}}_{\mathrm{pr}}^{-1}+h\sum _{k=0}^{K-1}{\widehat{\mathbf{U}}}_{k}^{T}{\widehat{\mathbf{D}}}^{-1}{\widehat{\mathbf{U}}}_{k}.$$

(13)

The mean values of **c** and $\widehat{\mathbf{D}}$ in the posterior distribution give the best estimates for the model parameters for a given block of data $\mathcal{S}$ of length *K* and provide the global minimum of ${L}_{s}(\mathbf{c},\widehat{\mathbf{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 $\widehat{\mathbf{D}}$ has a mean ${\widehat{\mathbf{D}}}_{\text{post}}^{\prime}={\widehat{\u03f4}}_{\mathrm{s}}\left(\mathbf{c}\right)$ that provides a minimum of ${S}_{\mathrm{s}}(\mathbf{c},\widehat{\mathbf{D}})$ with respect to $\widehat{\mathbf{D}}={\widehat{\mathbf{D}}}^{T}$. Its matrix elements are

$${\widehat{\u03f4}}_{\mathrm{s}}^{nn\prime}\left(\mathbf{c}\right)\equiv \frac{1}{K}\sum _{k=0}^{K-1}{[{\stackrel{.}{\mathbf{s}}}_{k}-\widehat{\mathbf{U}}\left({\mathbf{s}}_{k}\right)\mathbf{c}]}_{n}{[{\stackrel{.}{\mathbf{s}}}_{k}-\widehat{\mathbf{U}}\left({\mathbf{y}}_{k}\right)\mathbf{c}]}_{n\prime}^{T}.$$

(14)

Alternatively, assume next that $\widehat{\mathbf{D}}$ is known, and note from (10) that in this case the posterior distribution over **c** is Gaussian. Its covariance is given by ${\widehat{\Xi}}_{\mathrm{s}}\left(\widehat{\mathbf{D}}\right)$ and the mean ${\mathbf{c}}_{\text{post}}^{\prime}$ minimizes ${L}_{\mathrm{s}}(\mathbf{c},\widehat{\mathbf{D}})$ with respect to **c**

$${\mathbf{c}}_{\text{post}}^{\prime}={\widehat{\Xi}}_{\mathrm{s}}^{-1}\left(\widehat{\mathbf{D}}\right){\mathbf{w}}_{\mathrm{s}}\left(\widehat{\mathbf{D}}\right).$$

(15)

We repeat this two-step optimization procedure iteratively, starting from some prior values **c**_{pr} and ${\widehat{\Sigma}}_{\mathrm{pr}}$.

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 *s _{r}*(

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, *x _{r}*(

$${b}_{n}{y}_{n}\left({t}_{k}\right)=\frac{{s}_{n}({t}_{k}+h)-{s}_{n}\left({t}_{k}\right)}{h}+{a}_{n}{s}_{n}\left({t}_{k}\right),$$

(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**(*t _{k}*)}. 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

$$\stackrel{.}{\mathbf{y}}=\widehat{\mathbf{U}}(\mathbf{s},\mathbf{y})\mathbf{c}+\xi \left(t\right),$$

(17)

where *ξ*(*t*) is a two-dimensional Gaussian white noise with correlation matrix $\widehat{\mathbf{D}}$, and the matrix $\widehat{\mathbf{U}}$ will have the following block structure

$$\widehat{\mathbf{U}}=\left[\left[\begin{array}{cc}\hfill {\varphi}_{1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\varphi}_{1}\hfill \end{array}\right]\cdots \left[\begin{array}{cc}\hfill {\varphi}_{2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\varphi}_{2}\hfill \end{array}\right]\cdots \left[\begin{array}{cc}\hfill {\varphi}_{B}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\varphi}_{B}\hfill \end{array}\right]\right].$$

(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* = 2*B*.

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 $\langle {\xi}_{c}^{2}\left(t\right)\rangle =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 *D*_{11} = 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.

(Color online) (a) Time series of the cardiac oscillations *x*_{c}(*t*_{n}) = *s*_{c}(*t*_{n}) 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 **...**

(Color online) (a) Time series of the respiratory oscillations *x*_{r}(*t*_{n}) = *s*_{r}(*t*_{n}) 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.

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*) = *x _{r}* (

First we verify that the decomposition of the input signal *s*(*t*) into low-frequency ${\stackrel{~}{s}}_{r}$ and high-frequency ${\stackrel{~}{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 *y _{r}*(

(Color online) (a) Velocity of the respiratory oscillations of the original signal *y*_{r}(*t*) (green line) is compared with the signal ${\stackrel{~}{y}}_{r}\left({t}_{k}\right)$ (black dashed line) obtained as a result of filtration of *s*(*t*) followed by the embedding ${b}_{1}{\stackrel{~}{y}}_{r}\left({t}_{k}\right)=({s}_{r}({t}_{k}+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.

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.

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 {*x _{r}*(

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.

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

$${D}_{i}={S}_{i-1}\mathrm{exp}[(-2\u22153){I}_{i-1}\u2215T],$$

(19)

$${S}_{i}={D}_{i}+\gamma {I}_{i-1}+{C}_{1}+A\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\left(2\pi ft\right),$$

(20)

$${I}_{i}={G}_{v}{S}_{i-{\tau}_{v}}^{\prime}+{G}_{\beta}F(S\prime ,{\tau}_{\beta})+{C}_{2},$$

(21)

Here *C*_{1}, *C*_{2}, and *C*_{3} are constants and the sigmoidal nature of the baroreceptor sensitivity is accounted for by defining an effective systolic pressure (*S′*) [57]

$${S}_{i}^{\prime}={S}_{0}+18\phantom{\rule{thinmathspace}{0ex}}\mathrm{arctan}[(S-{S}_{0})\u221518].$$

(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 ${G}_{v}{S}_{i-{\tau}_{v}}^{\prime}$; and (ii) the slower *β*-sympathetic control of the heart rate *G _{β}F*(

$$F(S\prime ,\tau )=\sum _{k=-M}^{M}{a}_{k}{S}_{i-\tau +k}^{\prime}=({S}_{i-\tau -2}^{\prime}+2{S}_{i-\tau -1}^{\prime})+3{S}_{i-\tau}^{\prime}+2{S}_{i-\tau +1}^{\prime}+{S}_{i-\tau +2}^{\prime})\u22159$$

Further we assume for simplicity that the pressure oscillations do not deviate far from the working point *S*_{0} in (22), i.e. *S′* ≈ *S*.

(a) The BP signal in the frequency range of cardiac oscillations (black line). Systolic pressure (*S*_{n}), diastolic pressure (*D*_{n}), RR intervals (*I*_{n}), 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 *I _{n}* 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

$$\{\begin{array}{c}\stackrel{.}{x}=\u220a(y-\beta x),\hfill \\ \stackrel{.}{y}=\alpha y+\gamma {y}^{2}+\delta {y}^{3}-x+C,\hfill \end{array}\phantom{\}}$$

(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 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

$${x}_{0}\left(t\right)=\{\begin{array}{c}({S}_{n-1}+\frac{\mid {a}_{-}\mid}{\beta}){e}^{-\beta t}-\frac{\mid {a}_{-}\mid}{\beta},\phantom{\rule{thinmathspace}{0ex}}\text{for}\phantom{\rule{thinmathspace}{0ex}}0<t<\frac{2}{3}{I}_{n};\hfill \\ ({D}_{n}-\frac{\mid {a}_{+}\mid}{\beta}){e}^{-\beta t}+\frac{\mid {a}_{+}\mid}{\beta},\phantom{\rule{thinmathspace}{0ex}}\text{for}\phantom{\rule{thinmathspace}{0ex}}\frac{2}{3}{I}_{n}<t<{I}_{n}.\hfill \end{array}\phantom{\}}$$

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

Time evolution of the dynamical variables *x* (solid line) and *y* (dashed line) of the FHN system with parameters: = 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.

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.

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

[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] Scafer 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |