Abstract

Multichannel Blind Deconvolution (MBD) is a powerful tool particularly for the identification and estimation of dynamical systems in which a sensor, for measuring the input, is difficult to place. This paper presents an MBD method, based on the Malliavin calculus MC (stochastic calculus of variations). The arterial network is modeled as a Finite Impulse Response (FIR) filter with unknown coefficients. The source signal central arterial pressure CAP is also unknown. Assuming that many coefficients of the FIR filter are time-varying, we have been able to get accurate estimation results for the source signal, even though the filter order is unknown. The time-varying filter coefficients have been estimated through the proposed Malliavin calculus-based method. We have been able to deconvolve the measurements and obtain both the source signal and the arterial path or filter. The presented examples prove the superiority of the proposed method, as compared to conventional methods.

1. Introduction

In this paper, we present a new approach to monitor central arterial pressure using the Multichannel Blind Deconvolution (MBD) [1, 2]. A multichannel blind deconvolution problem can be considered as natural extension or generalization of instantaneous Blind Source Separation (BSS) problem [3, 4]. The problem of BSS has received wide attention in various fields such as signal analysis and processing of speech, image [5, 6], and biomedical signals, especially, signal extraction, enhancement, denoising, model reduction, and classification problems [79].

The MBD is the technique that allows the estimation of both an unknown input and unknown channel dynamics from only channel outputs. Although one cannot place a sensor [10] to directly measure the input, yet, it may be recovered from the outputs that are simultaneously measured at the multiple branches of the system. The MBD technique distinguishes itself from other techniques that apply a predetermined transfer function [11, 12] to interpret sensor data. The other techniques cannot account for individual differences nor can they account for dynamic changes in the subject’s physiologic state.

The physiologic state of the cardiovascular (CV) system can be most accurately assessed by using the aortic blood pressure or CAP [7, 13, 14] and flow. However, standard measurement of these signals, such as catheter, entails costly and risky surgical procedures. Therefore, most of the practically applicable methods aim to monitor the CV system based on peripheral circulatory signals, for example, arterial blood pressure at a distant site. Various methods have been developed to relate the peripheral signals to the CV state. These include blind deconvolution [1517] methods for recovering the CAP signal from the upper-limb arterial blood pressure, and the estimation of CV parameters such as left ventricular elasticity, end diastolic volume, total peripheral resistance (TPR), and mean aortic flow from arterial BP measurement [18, 19]. A chronic challenge of these previous methods is that the dynamics of the CV system, which relates the aortic and peripheral signals, is unknown and time-varying as well. So this problem turns out to be an ill-posed system identification problem because we are asked to identify both the unknown system dynamics and input signal using the output signal measurement alone.

Because of the practical difficulty in measuring the arterial pressure waveform near the heart [13, 15, 20], several mathematical transformation methods have been developed based on a generalized transfer function approach [21, 22]. Over the years many methods have been suggested for blind deconvolution for estimating central aortic pressure and flow [3, 9, 19, 23, 24]. One of the most popular and effective techniques is to assume an FIR model [1] and IIR model [18, 19, 25] for the modulating channels/paths and to estimate the coefficients [26] of this model. Through the inversion of the FIR filter, we get the original source signals. The principal assumption underlying these methods is that the arterial tree properties are constant over time and between individuals. A few methods have therefore been more recently developed towards “individualizing” the transfer function [1, 12, 21] relating peripheral artery pressure to central aortic pressure. These methods essentially involve (1) modeling the transfer function [11, 23] with physiologic parameters, (2) estimating a subset of the model parameters [7, 15] from the peripheral arterial pressure waveforms and/or other measurements from an individual while assuming values for the remaining parameters, (3) constructing a transfer function based on the estimated and assumed parameter values [26], and (4) applying the transfer function to the measured peripheral arterial pressure waveforms to predict the corresponding central aortic pressure waveform. While these methods attempt to determine a transfer function that is specific to an individual over a particular time period, only a partial individualization is actually obtained. Perhaps, as a result, these methods have found only limited success with results not much, if at all, better than the generalized transfer function approach.

In this paper, we suggest characterizing the channels of the single-input, multioutput system model of the arterial tree by linear and time-variant FIR filters. If we make only one of the FIR filter parameters changing over time, then the problem is handled by the Ito calculus [27, 28], while, if we make more than one of the FIR filter parameters changing over time, then the problem could be handled by the Malliavin calculus [29] as we propose in this paper. This way, the ambiguity in the order of the impulse response is compensated by the time variations of the filter parameters [2931].

In this paper, we introduce a method that could estimate time-varying parameters/coefficients. It is based on the stochastic calculus of variations (Malliavin calculus) [28, 3235]. We derive a closed-form expression for the unknown time varying FIR filter parameters by using the Clark-Ocone formula [30]. This will enable us to find a stochastic differential equation (SDE) for each unknown time-varying parameter of the FIR filter. Each SDE is function of PAP and some other unknown deterministic parameters. The unknown deterministic parameters are estimated through Monte Carlo simulation methods.

The proposed method is then applied to noninvasive monitoring of the cardiovascular system of the swine. The arterial network is modeled as a multichannel system where the CAP is the input and pressure profiles measured at different branches of the artery, for example, radial and femoral arteries, are the outputs. The proposed method would allow us to estimate both the waveform of the input pressure and the arterial channel dynamics from outputs obtained with noninvasive sensors placed at different branches of the arterial network. Numerical examples verify the major theoretical results and the feasibility of the method. In Section 2, we describe the blind deconvolution problem, conventional solution methods. In Section 3 we introduce the proposed method based on the Malliavin calculus. In Section 4, we present the results for the reconstruction of single-input CAP from two distant measures outputs PAP. Finally in Section 5, we provide summary and conclusions. The Appendix contains the technical derivations of the proposed method.

2. Problem Formulation

2.1. Multichannel Dynamic Systems

The cardiovascular system is topologically analogous to a multichannel dynamic system. Pressure wave emanating from a common source, the heart, is broadcast and transmitted through the many vascular pathways. Therefore, noninvasive circulatory measurements taken at different locations (as shown in Figure 3) can be treated as multichannel data and processed with an MBD algorithm. Figure 1 illustrates, in a block diagram form, the relation between the central aortic pressure waveform 𝑢(𝑡) and the peripheral arterial pressure waveforms 𝑦𝑖(𝑡). The arterial channels 𝑖(𝑡) relating the common input to each output represent the vascular dynamic properties of different arterial tree paths as characterized by finite impulse responses (FIRs). The main idea is therefore to determine the absolute central aortic pressure waveform within an arbitrary scale factor by mathematically analyzing two or more PAP waveforms or related signals so as to extract their common features. An ancillary idea is to then determine the parameters of the determined central aortic pressure waveform.

We estimate the FIR filter coefficients by the conventional methods in Section 2.2. We introduce our proposed method that is based on Malliavin calculus in Section 3. We will be working in the probability space (Ω,,𝒫). To simplify the exposure, we shall assume that we have only two measurements outputs of a modulated version of the source signal that are given as𝑦1(𝑡)=1(𝑡)𝑢(𝑡)+𝜀1𝑦(𝑡),(2.1)2(𝑡)=2(𝑡)𝑢(𝑡)+𝜀2(𝑡),(2.2) where 𝑢(𝑡) is the unknown source signal (central AP), 1(𝑡) and 2(𝑡) are unknown filters (hemodynamic response at time 𝑡) or arterial paths, “” is the convolution operation, 𝑦1(𝑡) (femoral AP) and 𝑦2(𝑡) (radial AP) are the observed measurements, and 𝜀1(𝑡) and 𝜀2(𝑡) are the measurements noise. The objective is to deconvolve 𝑦1(𝑡) and 𝑦2(𝑡) to estimate 𝑢(𝑡). If we convolve 𝑦1(𝑡) with 2(𝑡), we will get2(𝑡)𝑦1(𝑡)=2(𝑡)1(𝑡)𝑢(𝑡)+2(𝑡)𝜀1(𝑡),(2.3) since the convolution is a commutative operation, then exchanging 1(𝑡) and 2(𝑡) and on the right-hand side we get: 2(𝑡)𝑦1(𝑡)=1(𝑡)2(𝑡)𝑢(𝑡)+2(𝑡)𝜀1(𝑡)=1(𝑡)𝑦2(𝑡)1(𝑡)𝜀2(𝑡)+2(𝑡)𝜀2(𝑡).(2.4) Thus, 2(𝑡)𝑦1(𝑡)=1(𝑡)𝑦2(𝑡)1(𝑡)𝜀2(𝑡)+2(𝑡)𝜀1(𝑡).(2.5) Note that this equation does not include the unknown input 𝑢(𝑡). It represents the constraints among the channel dynamics or filters and observed output. Substituting a measured time series of output data for 𝑦1(𝑡) and 𝑦2(𝑡), the above equation can be solved for the unknown parameters involved in 1(𝑡) and 2(𝑡). Once the filters are obtained, we will use their inverses to find an estimate for the source signal. To simplify the exposure further, assume that the modulating filters, that represent the signal paths or channel dynamics, are second-order linear time invariant and have the 𝑍 transforms: 1(𝑧)=1+𝛽1𝑧1+𝛽2𝑧2,(2.6) that is, 𝑦1(𝑘)=𝑢(𝑘)+𝛽1𝑢(𝑘1)+𝛽2𝑢(𝑘2)+𝜀1(𝑘),(2.7) and in matrix format for 𝑁 data points:𝑦1=𝛽(2)𝑦(3)𝑦(𝑁1)2𝛽11000𝛽2𝛽11000𝛽2𝛽11+𝜀𝑢(0)𝑢(1)𝑢(𝑁1)1𝜀(2)1(𝜀3)1(𝑁1),(2.8) that is, 𝑌1=𝐻1𝑈+𝜀1,(2.9) where𝑌1=𝑦1(2)𝑦1(3)𝑦1(𝑁1)𝑇,𝑈(2.10)=[]𝑢(0)𝑢(1)𝑢(𝑁1)𝑇𝜀,(2.11)1=𝜀1(2)𝜀1(3)𝜀1(𝑁1)𝑇𝐻,(2.12)1=𝛽2𝛽11000𝛽2𝛽11000𝛽2𝛽11,(2.13) and “𝑇” stands for transpose.

Similarly, 2(𝑧)=1+𝛼1𝑧1+𝛼2𝑧2,(2.14) that is, 𝑦2(𝑘)=𝑢(𝑘)+𝛼1𝑢(𝑘1)+𝛼2𝑢(𝑘2)+𝜀2(𝑘),(2.15) where 1(𝑧) and 2(𝑧) are the 𝑍 transforms of the discrete versions of 1(𝑡) and 2(𝑡), respectively. Since, 2(𝑡)𝑦1(𝑡)=1(𝑡)𝑦2(𝑡)1(𝑡)𝜀2(𝑡)+2(𝑡)𝜀1(𝑡),(2.16) taking the 𝑍 transform of the discrete version of both sides, we get2(𝑧)𝑦1(𝑧)=1(𝑧)𝑦2(𝑧)1(𝑧)𝜀2(𝑧)+2(𝑧)𝜀1(𝑧).(2.17) In the time domain, we get the equation: 𝑦1(𝑘)+𝛼1𝑦1(𝑘1)+𝛼2𝑦1(𝑘2)=𝑦2(𝑘)+𝛽1𝑦2(𝑘1)+𝛽2𝑦2+(𝑘2)𝜀2(𝑘)𝛽1𝜀2(𝑘1)𝛽2𝜀2(𝑘2)+𝜀1(𝑘)+𝛼1𝜀1(𝑘1)+𝛼2𝜀1.(𝑘2)(2.18)

2.2. A Conventional Method for the Estimation of the System

The familiar scalar regression format as shown in (2.18) is the shape of a regression equation with a correlated error (colored noise). Unless we take this into consideration, the ordinary least square (OLS) method will yield biased estimates for the unknown coefficients 𝛼1,𝛼2,𝛽1,and𝛽2.

Rearrange (2.18), we gets 𝑎𝑦1(𝑘)𝑦2(𝑘)=𝛼1𝑦1(𝑘1)𝛼2𝑦1(𝑘2)+𝛽1𝑦2(𝑘1)+𝛽2𝑦2+(𝑘2)𝜀2(𝑘)𝛽1𝜀2(𝑘1)𝛽2𝜀2(𝑘2)+𝜀1(𝑘)+𝛼1𝜀1(𝑘1)+𝛼2𝜀1(.𝑘2)(2.19) For the general case where the order of the FIR filters are 𝐼 and 𝐽, we get 𝑦1(𝑘)𝑦2(𝑘)=𝐼𝑖=1𝛼𝑖𝑦1(𝑘𝑖)+𝐽𝑗=1𝛽𝑗𝑦2+(𝑘𝑗)𝐽𝑗=1𝛽𝑗𝜀2(𝑘𝑗)+𝐼𝑖=1𝛼𝑖𝜀1+𝜀(𝑘𝑖)1(𝑘)𝜀2.(𝑘)(2.20) This could be approximated as𝑦1(𝑘)𝑦2(𝑘)𝐼𝑖=1𝛼𝑖𝑦1(𝑘𝑖)+𝐽𝑗=1𝛽𝑗𝑦2(𝑘𝑗)+𝜀(𝑘).(2.21) Since the filter orders are unknown, one could use the corrected Akaike information criterion AIC𝑐 to determine both “𝐼” and “𝐽”. Assume that the error term, 𝜀(𝑘), is zero mean Gaussian with variance 𝜎2, the AIC𝑐 is defined as [36]AIC𝑐=𝑛ln𝜎2++12𝑛(𝑞+1)(𝑛𝑞2),(2.22) where 𝑛 is the number of observations, 𝑞=𝐼+𝐽 is a number of unknown, and 𝜎2 is an estimate of the error variance. We choose the order 𝑞 such that AIC𝑐 is minimized.

Once the coefficients of the FIR filter are estimated, we use inverse filtering to find an estimate for the source signal 𝑈 as follows: 𝑈=𝐻𝑇1𝐻1𝐻𝑇11𝑌1,(2.23) where 𝑌1=𝑦1(2)𝑦1(3)𝑦1(𝑁1)𝑇,𝑈(2.24)=[]̂𝑢(0)̂𝑢(1)̂𝑢(𝑁1)𝑇𝐻,(2.25)1=̂𝛽2̂𝛽10̂𝛽1002̂𝛽1̂𝛽10002̂𝛽11.(2.26) The symbol “” on top of the variable refers to estimation. For example, ̂𝛽2 and ̂𝛽1 are the estimates for 𝛽2 and 𝛽1, respectively. Due to its simplicity, the above mentioned method is the one that is commonly used [19].

3. The Malliavin Calculus and the Generalized Clark-Ocone Formula for the Estimation of the Unknown Time-Varying Coefficients: The Proposed Approach

The linear filter assumption is just an approximation to reality. Sometimes the media (arterial tree) is nonlinear, time varying, random, or all. Moreover and above the measured signals are usually noisy. The filter order is usually unknown. All the factors suggest that the FIR filter model is just an approximation. To compensate for these approximations, we suggest making some of the unknown filter coefficients varying with time, for example, (𝛼1(𝑡),𝛼2(𝑡),𝛽1(𝑡),𝛽2(𝑡)). Notice that in this case, the filter coefficients do not have any physical meaning and we are only concerned with the estimation of their values. Now the problem becomes that of the estimation of the unknown time-varying coefficients. The details of the estimation procedure are, briefly, given in this section. More details could be found in the appendix.

We now recast the problem in the format that could be handled by the Malliavin calculus. Using (2.21), the observed signal/regressor, 𝑝𝑦(𝑘)=𝑦1(𝑘)𝑦2(𝑘), is the sum of stochastic processes, 𝑦1(𝑘1),𝑦1(𝑘2),,𝑦2(𝑘1),𝑦2(𝑘2), weighted by the unknown time-varying parameters, 𝛼1(𝑘),𝛼2(𝑘),,𝛽1(𝑘),𝛽2(𝑘),. The objective is to estimate, from the observations, the unknown time-varying coefficients 𝛼1(𝑘),𝛼2(𝑘),,𝛽1(𝑘),𝛽2(𝑘), .

3.1. The Estimates for the Time-Varying Filter Parameters

To solve the estimation of the time-varying parameters problem, we imbed the sum of the processes (the observations) into another signal, 𝑝(𝑡), which we call the augmented observed signal. The augmented observations are the original observations plus a deterministic component. The addition of the known deterministic component is needed to facilitate the analysis.

The augmented observed signal, 𝑝(𝑡), with (𝑑) components, could be modeled, in the continuous time, as follows:𝑝(𝑡)=𝜂1(𝑡)𝑆1(𝑡)+𝑑𝑖=2𝜂𝑖(𝑡)𝑆𝑖[](𝑡)𝑡0,𝑇=𝜂1(𝑡)𝑆1(𝑡)+𝑝𝑦(𝑡),(3.1) where 𝑝𝑦(𝑡)=[𝑦1(𝑡)𝑦2(𝑡)]; see (2.21). 𝜂𝑖(𝑡),𝑖>1, is the 𝑖th unknown time-varying coefficient and are defined as follows for second-order filters:

𝜂2(𝑡)=𝛼1(𝑡),𝜂3(𝑡)=𝛼2(𝑡),𝜂4(𝑡)=𝛽1(𝑡),𝜂5(𝑡)=𝛽2(𝑡),𝑆2(𝑡)=𝑦1(𝑡Δ),𝑆3(𝑡)=𝑦1(𝑡2Δ),𝑆4𝑦(𝑡)=2(𝑡Δ),𝑆5(𝑡)=𝑦2(𝑡2Δ),Δ is the sampling interval, 𝜂1(𝑡) = known constant, with 𝜂1(𝑡)=knownconstant,(3.2)𝑑𝑆1(𝑡)=𝜆1𝑆1(𝑡)𝑑𝑡,(3.3) where 𝜆1 and 𝑆1(𝑡) are known.

In this analysis, the presence of a deterministic component will make the exposition easier. In this case, It is 𝑆1(𝑡). This component acts as a reference signal or a numeraire. Since the observed signal does not usually come with a known deterministic component, we should add a known deterministic component and proceed with the analysis. We will then try to find an expression for the SDE that describes the evolution of 𝑝(𝑡). This expression will be a function of the unknown stochastic coefficients 𝜂𝑖(𝑡). Assuming a shape for the final value 𝑝(𝑇), as a function of the Wiener processes 𝑊(𝑡), and using the generalized Clark-Ocone formula [30, 37, 38], we will be able to find a closed-form expression for the coefficients 𝜂𝑖(𝑡) as a function of the Wiener processes 𝑊(𝑡). The details of the derivations are given in the appendix and are also given in [29].

Assume that the unknown coefficients, 𝜂𝑖(𝑡), are varying slowly with time, that is𝑑𝜂𝑖(𝑡)0,(3.4) or more precisely, 𝑑𝜂𝑖(𝑡)𝑑𝑆𝑖(𝑡).

There are several commonly models that represent different physical situations. In our case, the signals, 𝑆𝑖(𝑡), each has an SDE of the form of an Ornstein-Uhlenbeck (OU) process with trend𝑑𝑆𝑖(𝑡)=𝑐𝑖𝑂(𝑡)𝑖(𝑡)𝑆𝑖(𝑡)𝑑𝑡+𝑒𝑖(𝑡)𝑑𝑊𝑖(𝑡),𝑖>1.(3.5) This form represents an Ornstein-Uhlenbeck (OU) process and has the solution [12, 39]𝑆𝑖(𝑡)=𝑒𝐶𝑖(𝑡)𝑆𝑖(0)+𝑡0𝑐𝑖(𝑢)𝑂𝑖(𝑢)𝑒𝐶𝑖(𝑢)𝑑𝑢+𝑒𝐶𝑖(𝑡)𝑡0𝑒𝑖(𝑢)𝑒𝐶𝑖(𝑢)𝑑𝑊𝑖(𝑢),𝑖>1,(3.6) where 𝐶𝑖(𝑡)=𝑡0𝑐𝑖(𝑢)𝑑𝑢.

The coefficients 𝑐𝑖(𝑡),𝑒𝑖(𝑡), and the trend 𝑂𝑖(𝑡) are estimated from the observed data. For constant quantities 𝑐𝑖,𝑂𝑖,𝑒𝑖, and after some manipulations we end up with the estimate of the unknown time-varying parameters (see the appendix):𝜂𝑖𝑒(𝑡)=𝜆1𝑡𝑒𝑖𝐸𝑄𝐷𝑖𝑡𝐹(𝑇)𝐹(𝑇)𝑇𝑡𝐷𝑖𝑡𝜃𝑖𝑠,𝑊𝑖𝑑𝑊(𝑠)𝑖(𝑠)𝑡,(3.7) where 𝜃𝑖𝑢,𝑊𝑖=𝑐(𝑢)𝑖𝜆1𝑒𝑖𝑆𝑖(0)𝑒𝑐𝑖𝑢+𝑢0𝑒𝑖𝑒𝑐𝑖(𝑢𝑟)𝑑𝑊𝑖(𝑟),(3.8)𝐹(𝑇)=𝑒𝜆1𝑇𝑝(𝑇),(3.9)𝐸𝑄{𝑋/𝑡} is the conditional expectation of the random variable X under the probability measure 𝑄 [37] and 𝑋=[𝐷𝑖𝑡𝐹(𝑇)𝐹(𝑇)𝑇𝑡𝐷𝑖𝑡𝜃𝑖(𝑠,𝑊𝑖𝑊(𝑠))𝑑𝑖(𝑠)]. This measure 𝑄 is related to the old probability measure 𝒫 through a function of the stochastic process 𝜃𝑖(𝑢,𝑊𝑖(𝑢)).𝑡 is the filtration generated by the Wiener vector 𝑊(𝑡). 𝑡[0,𝑇]. 𝑡 contains exactly the information learned by observing the Wiener vector processes up to time 𝑡. 𝐷𝑖𝑡𝜃𝑖(𝑠,𝑊𝑖(𝑠))is the Malliavin derivative of 𝜃𝑖(𝑠) with respect to the 𝑖th Wiener process 𝑊𝑖(𝑠)𝐹(𝑇) is, in general, a function of the Wiener vector 𝑊(𝑇). 𝐸𝑄{𝑋/𝑡} is usually difficult to calculate. 𝜂𝑖(𝑡) is in general a function of 𝑊(𝑡), and it could be approximated by the SDE:𝑑𝜂𝑖(𝑡)=𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗(𝑡),(3.10) and 𝑒𝑖𝜂𝑖(𝑡) could be represented by the SDE𝑑𝑒𝑖𝜂𝑖(𝑡)=𝑒𝑖𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑒𝑖𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗(𝑡).(3.11) Notice that the drift and the diffusion are related to each other through 𝜃𝑖(𝑡,𝑊𝑖(𝑡)) [29, 32] see the appendix. Let 𝑒𝑖𝜇𝑖𝑊(=𝑡),𝑡𝑗𝜙𝑖𝑗𝑡,𝑊(𝜃𝑡)𝑗𝑡,𝑊(𝑡),(3.12) then 𝑒𝑖𝜅𝑖𝑗𝑊(𝑡),𝑡=𝜙𝑖𝑗𝑡,𝑊(𝑡).(3.13) The last equality comes from the solution of the CO formula 𝜃𝑖(𝑡,𝑊𝑖(𝑡)) (see the appendix). An expression for 𝜙𝑖𝑗(𝑡,𝑊(𝑡)) is needed to solve for the unknowns 𝜂𝑖(𝑡). An exact expression, however, is difficult to obtain. Instead, we assume a polynomial shape with unknown parameters (𝛿𝑖𝑗,𝛾𝑘𝑖𝑗,𝜆𝑖𝑗𝑘𝑙), namely,𝜙𝑖𝑗𝑡,𝑊(𝑡)=𝛿𝑖𝑗+𝑘𝛾𝑘𝑖𝑗𝑊𝑘(𝑡)+𝑘𝑙𝜆𝑖𝑗𝑘𝑙𝑊𝑘(𝑡)𝑊𝑙(𝑡)+.(3.14) The unknown parameters are estimated by minimizing the squared error, 𝑒2(𝑇), between the observed, 𝑝(𝑡), and the estimated observations ̂𝑝(𝑡), wherê𝑝(𝑡)=𝑖̂𝜂𝑖(𝑡)𝑆𝑖𝑑𝑒(𝑡),𝑖1,(3.15)𝑖̂𝜂𝑖(𝑡)=𝑒𝑖𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑒𝑖𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗𝑒(𝑡),(3.16)𝑖𝜇𝑖𝑊=(𝑡),𝑡𝑗𝜙𝑖𝑗𝑡,𝑊𝜃(𝑡)𝑗𝑡,𝑊,𝑒(𝑡)(3.17)𝑖𝜅𝑖𝑗𝑊=𝜙(𝑡),𝑡𝑖𝑗𝑡,𝑊,𝜙(𝑡)(3.18)𝑖𝑗𝑡,𝑊(𝑡)=𝛿𝑖𝑗+𝑘𝛾𝑘𝑖𝑗𝑊𝑘(𝑡)+𝑘𝑙𝜆𝑖𝑗𝑘𝑙𝑊𝑘(𝑡)𝑊𝑙𝑒(𝑡),(3.19)2(𝑇)=𝑇0[]𝑝(𝑡)̂𝑝(𝑡)2𝑑𝑡.(3.20)

Once the parameters 𝜂𝑖(𝑡) are estimated, we use inverse filtering to find 𝑈 (the aortic pressure waveform). In the case that the filter order is two, that is, only two time-varying parameters, ̂𝜂2(𝑡)=𝛼1(𝑡),̂𝜂3(𝑡)=𝛼2(𝑡),̂𝜂4̂𝛽(𝑡)=1(𝑡), and ̂𝜂5̂𝛽(𝑡)=2(𝑡). In this case, the aortic pressure 𝑈 is estimated as 𝑈=𝐻𝑇1𝐻1𝐻𝑇11𝑌1,(3.21) where 𝑌1=𝑦1(2)𝑦1(3)𝑦1(𝑁1)𝑇,𝑈=[]̂𝑢(0)̂𝑢(1)̂𝑢(𝑁1)𝑇,𝐻1=̂𝛽2𝑡1̂𝛽1𝑡20̂𝛽1002𝑡2̂𝛽1𝑡3̂𝛽10002𝑡𝑁2̂𝛽1𝑡𝑁11.(3.22) We note from matrix that the parameters ̂𝛽1(𝑡) and ̂𝛽2(𝑡) have values that are changing across the sample.

In order to further reduce the mathematical complexity, one could restrict the analysis to the estimation of only two coefficients 𝛼1(𝑡) and 𝛽1(𝑡). The rest of filters coefficients are considered constants and the OLS estimates of them are used in the equations.

3.2. Summary of the Proposed Algorithm and its Assumptions

Our MBD technique is based on a set of assumptions as follows.(1)The common input is obtained at the output of two sensors. The channels relating the common input to each distinct output, as in Figure 1, are linear time variant. (2)The system is represented by FIR filters. The filters 1(𝑡) and 2(𝑡) are second-order and contain many parameters some of which are time-varying.(3)The FIR filter parameters 𝛼1(𝑡),𝛼2(𝑡),𝛽1(𝑡),𝑎𝑛𝑑𝛽2(𝑡) are varying slowly with time.(4)An OU model is assumed to describe the signals/regressors 𝑆𝑖(𝑡) which represent the processes 𝑦1(𝑡Δ),𝑦1(𝑡2Δ),𝑦2(𝑡Δ), and 𝑦2(𝑡2Δ)); see (3.1) and (3.6).(5)The time-varying FIR filter parameters 𝜂𝑖(𝑡) are estimated by using Malliavin calculus and the generalization Clark-Ocone formula; see (3.7). (6)Equataion (3.21) is used with the time-varying parameters/coefficients ̂𝜂𝑖(𝑡) to find an estimate for the input AP.

An absolute central aortic pressure waveform is determined by scaling the estimated input based on the measured waveforms. In certain embodiments, the reconstructed waveform is calibrated to absolute pressure based on the measured peripheral artery pressure waveforms. For example, the reconstructed waveform is scaled to have a mean value similar to the mean value of the measured waveforms. This scaling step is well justified, since the paths from the central aorta to peripheral arteries offer very little resistance to blood flow due to Poiseuille's law [1]. Certain embodiments scale the reconstructed waveform to have a mean value specifically equal to that of one of the measured peripheral artery pressure waveforms or the waveform with the largest mean value plus a constant (whose value may be between, e.g., 0 and 3 mmHg). Certain alternative embodiments scale the reconstructed waveform to have a mean value equal to the mean (or medium) of the mean values of some or all of the measured peripheral artery pressure waveforms plus a constant (whose value may be between, e.g., 0 and 3 mmHg).

The reconstructed waveform 𝑈 is calibrated to absolute pressure by scaling it to have the same mean value as the measured peripheral artery pressures as follows:𝑈𝑆(𝑡)=𝑈(𝑡)𝑁1𝑡=0𝑌1(𝑡)+𝑌2(𝑡)/2𝑁1𝑡=0[]𝑈(𝑡),𝑡0,𝑁1,(3.23) where 𝑈𝑆(𝑡) is the final absolute (scaled) estimated aortic pressure, that is, 𝑈𝑆(𝑡)=𝑈(𝑡) × Scale factor, Scalefactor=averagefemoralpressure/2+averageradialpressure/2averageestimatedaorticpressure.(3.24)

3.3. The Proposed Algorithm

We now give a brief summary in an algorithmic form to describe the methodology as follows.(1)Insert data in [1] as shown in Figure 2.(2)Describe the form of the stochastic process OU using the formula in (3.5).(3)Estimate the parameters 𝜂𝑖(𝑡) by the following steps for each unknown time-varying coefficient 𝜂𝑖(𝑡).(i)Initiate an estimate for the unknown parameters 𝛿𝑖𝑗,𝛾𝑘𝑖𝑗,𝜆𝑖𝑗𝑘𝑙.(ii)Use Monte Carlo simulation method to generate Wiener processes 𝑊(𝑡).(iii)Use (3.19) to an estimate for 𝜙𝑖𝑗(𝑡,𝑊(𝑡)).(iv)Use (3.16) to find an estimate for ̂𝜂𝑖(𝑡).(v)Calculate the summed squared error of (3.20).(vi)Change the estimate for the unknown parameters 𝛿𝑖𝑗,𝛾𝑘𝑖𝑗,and𝜆𝑖𝑗𝑘𝑙 using the least square method, for example, such that the error is reduced.(vii)Repeat steps (ii)–(vi) and stop when the error is minimum or the number of iterations is exceeded.(4)Calculate the matrix 𝐻1 in (3.25).(5)Estimate the source signal by (3.21).(6)Calculate the estimated aortic pressure within scale factor by (3.23).

4. Computed Examples

To test the proposed approach, first, we took real data from [1] by inserting the graph that contains the measured data on software CURVESCAN and extracting points (see Figure 2). Second, we simulated a set of 300 data points on computer measured data.

Multichannel blind deconvolution was experimentally evaluated with respect to measured data in which femoral artery pressure, radial artery pressure waveforms, and aortic pressure waveform were simultaneously measured (see Figure 3). Third, we demonstrate the ability of the proposed approach to extract aorta AP waveform from multichannel. All the calculations in the algorithm were performed under Excel by using Risk Solver Platform V9.0 [40, 41]. Fourth, we evaluated the proposed method by the performance measure SNRE. The signal-to-noise ratio of the estimates (SNRE) was taken as the measure of performance for this evaluation. It is defined asSNRE=10log𝑘𝑢2(𝑘)𝑘[]𝑢(𝑖)̂𝑢(𝑖)2𝑑𝑏,(4.1) where ̂𝑢(𝑖) is the estimated value of the pressure at instant “𝑖” (see Table 1).

In the first example, it was assumed that the data were noise free but the order of each of the filters was unknown. In the second example, the data was assumed noisy and the order of each of the filters was unknown. The filter order was estimated using a corrected Akaike information criterion (AIC) [36]. In both examples, by using the proposed method, all coefficients of the two filters were assumed to be unknown slowly time-varying, 𝜂2(𝑡),,𝜂𝑑(𝑡). It was assumed that the signals/regressors 𝑆𝑖(𝑡) follow an Ornstein-Uhlenbeck (OU) process [31] with some unknown constant parameters see; (3.5).

In both cases, the proposed approach outperformed the conventional method. The pressure at the root of the aorta (central AP) was successfully estimated.

4.1. Example  1 Noise-Free Measurements

We have two measurements, 𝑦1(𝑘) (femoral AP) and 𝑦2(𝑘) (radial AP) that are assumed to be represented by the equations 𝑦1(𝑘)=𝑢(𝑘)+𝛽1𝑢(𝑘1)+𝛽2𝑦𝑢(𝑘2),(4.2)2(𝑘)=𝑢(𝑘)+𝛼1𝑢(𝑘1)+𝛼2𝑢(𝑘2)+𝛼3𝑢(𝑘3)+𝛼4𝑢(𝑘4).(4.3) In the conventional method, the order of the filters was assumed to be just two, that is, 𝛼3=𝛼3=0. In the proposed method, we made the same order assumption of two but all coefficients are slowly time-varying. Specifically, we assumed that the SDE of the signals, 𝑆𝑖(𝑡), of each has the form of the Ornstein-Uhlenbeck (OU) process 𝑑𝑆𝑖(𝑡)=𝑐𝑖𝑂𝑖(𝑡)𝑆𝑖(𝑡)𝑑𝑡+𝑒𝑖𝑑𝑊𝑖(𝑡),𝑖>1,(4.4) where 𝑐𝑖,𝑒𝑖 are unknown constants, 𝑂𝑖(𝑡) is a polynomial in time that represents the trend or the base data, and 𝑆𝑖(𝑡) bounces around 𝑂𝑖(𝑡). Remember that in the conventional method 𝑦1(𝑘)𝑦2(𝑘)=𝛼1𝑦1(𝑘1)𝛼2𝑦1(𝑘2)+𝛽1𝑦2(𝑘1)+𝛽2𝑦2(𝑘2).(4.5) We estimate the unknown coefficients 𝛼1,𝛼2,𝛽1,and 𝛽2, through regression analysis. In the proposed method, however, the slowly time-varying coefficients ̂𝜂𝑖(𝑡), 𝑖=2,3,4,5; are estimated by the (3.16) in Figure 5 and assumed ̂𝜂1(𝑡)=1.

Figure 4 shows the estimated pressure at the root of the aorta AP using the proposed method based on the Malliavin calculus (Section 3.2) and measured aorta or central AP, (SNRE = 25.111 db) (see Table 1). In Figure 7, we compared the estimated central arterial pressure using the conventional method (Section 2.2) with second-order FIR-2 (SNRE = 15.125 db) and the proposed Malliavin calculus-based method (SNRE = 25.111 db). So applying the conventional method with fourth-order FIR-4, we obtained a better result than by using FIR-2. But still the Malliavin calculus (MC) proved to be the best (see Table 1 and Figure 7). We note the difference between three methods in Figures 6 and 7. Therefore, typical estimates for ̂𝜂2(𝑡),̂𝜂3(𝑡),̂𝜂4(𝑡),̂𝜂5(𝑡) using the Malliavin calculus are shown in Figure 5, where ̂𝜂2(𝑡)=𝛼1(𝑡),̂𝜂3(𝑡)=𝛼2(𝑡),̂𝜂4̂𝛽(𝑡)=1(𝑡), and ̂𝜂5̂𝛽(𝑡)=2(𝑡).

4.2. Example  2 Noisy Measurements

In this example, white Gaussian noise was added to the original values of the measured pressures. It is assumed that the channels are represented by a second-order FIR filter and fourth order FIR filter (4.2) and (4.3). However, in this case, for conventional method, the order of the filters was estimated by using the corrected AIC. The estimated orders turned out to be 𝐼=4, 𝐽=2. Table 2 shows the estimated values of central AP using both methods.

Figure 9 shows the measured pressure at the root of the aorta (with added high noise) and the estimated pressure using the proposed method that is based on the Malliavin calculus. By using the proposed method based on Malliavin calculus (stochastic calculus of variations), the SNRE (as shown in Table 2) at noise variance 70 = 18.808 db, where the signal-to-noise ratio (SNR = 19.833 db). While, in Figure 8, the estimated central arterial pressure (with added the same noise variance) using the conventional method with second-order FIR-2 (SNRE = 12.748 db) and with four order FIR-4 (SNRE = 16.548 db). The performance measure (as shown in Figure 6) of the FIR filter-based methods and the proposed method are compared in Table 2 for different noise levels. The results show the superior performance of the proposed method.

5. Discussion and Conclusions

We have developed a new technique to estimate the aortic pressure waveform from multiple measured peripheral artery pressure waveforms. The technique is based on multichannel blind deconvolution in which two or more measured outputs (peripheral artery pressure waveforms) of a single input, multioutput system (arterial tree) are mathematically analyzed so as to reconstruct the common unobserved input (aortic pressure waveform). Each channel is modeled as an FIR filter. We assumed that all of the FIR filter parameters/coefficients are varying slowly with time. Their values were estimated using methods based on the Malliavin calculus. By this assumption, time-varying parameters, we were able to compensate for the wrong FIR filter order and the possible time variations/nonlinearities of the channels. The blind deconvolution problem was reformulated as a regression problem with unknown time-varying regression coefficients. We introduced a new method for the estimation of these slowly time-varying regression coefficients. The method relies heavily on the Malliavin calculus (stochastic calculus of variations) and the generalized Clark-Ocone formula to find a closed-form expression for the estimates of the unknown time-varying coefficients. Some approximations were needed to find mathematically tractable estimates. While the approach is quite general and can be applied to any signal model, we present only one signal model, the Ornstein Uhlenbeck process. Other models could have been used as well.

We tested the proposed technique in swine experiments, and our results showed superior performance for our proposed approach compared to conventional methods. Our way to reconstructed AP is simple and straightforward. Our method needs only the calculation of pressure wave components in the time domain and does not need calculations in the frequency domain and no need to large computer time. Because of this simplicity it is quite possible to implement this method in monitoring central pressure AP on-line. In the future, we suggest expanding this method by applying it to real data taken from human cardiovascular simulator. We suggest studying the case where the time-varying parameters/coefficients are rapidly changing over time. Another extension is the estimation of the flow at the root of the Aorta from peripheral measurements of the flow/pressure.

Appendix

A. Derivation of the Estimates [37, 38, 42]

Let the augmented signal 𝑝(𝑡) be defined as𝑝(𝑡)=𝑑𝑖=1𝜂𝑖(𝑡)𝑆𝑖(𝑡).(A.1) It is assumed that we have slowly varying unknown coefficients 𝜂𝑖(𝑡), that is;𝑑𝜂𝑖(𝑡)0,(A.2) or more precisely 𝑑𝜂𝑖(𝑡)𝑑𝑆𝑖(𝑡).(A.3) The signals, 𝑆𝑖(𝑡), each has an SDE of the form:𝑑𝑆𝑖(𝑡)=𝑎𝑖𝑡,𝑆𝑖(𝑡)𝑑𝑡+𝑏𝑖𝑡,𝑆𝑖(𝑡)𝑑𝑊𝑖(𝑡),𝑖>1,(A.4) which has a solution:𝑆𝑖(𝑡)=𝑆𝑖(0)+𝑡0𝑎𝑖𝑢,𝑆𝑖(𝑢)𝑑𝑢+𝑡0𝑏𝑖𝑢,𝑆𝑖(𝑢)𝑑𝑊𝑖(𝑢),𝑖>1,(A.5) and assume that 𝑑𝑆1(𝑡)=𝑎1𝑡,𝑆1(𝑡)𝑑𝑡=𝜆1𝑆1(𝑡)𝑑𝑡.(A.6) with known 𝑆1(0) and known 𝜆1.

Under the slowly varying-coefficient assumptions we have: 𝑑𝑝(𝑡)=𝑖𝜂𝑖(𝑡)𝑑𝑆𝑖(𝑡).(A.7) Substituting for 𝑑𝑆𝑖(𝑡), we get an expression for the SDE of 𝑝(𝑡) as𝑑𝑝(𝑡)=𝑖𝜂𝑖(𝑡)𝑎𝑖𝑡,𝑆𝑖(𝑡)𝑑𝑡+𝑖𝜂𝑖(𝑡)𝑏𝑖𝑡,𝑆𝑖(𝑡)𝑑𝑊𝑖(𝑡),(A.8) Rearrange (A.1), we get𝜂1(𝑡)=𝑝(𝑡)𝑖>1𝜂𝑖(𝑡)𝑆𝑖(𝑡)𝑆1(𝑡),(A.9) then, 𝑑𝑝(𝑡)=𝑖>1𝜂𝑖𝑎(𝑡)𝑖𝑡,𝑆𝑖𝑎(𝑡)1𝑡,𝑆1(𝑡)𝑆1𝑆(𝑡)𝑖+(𝑡)𝑑𝑡𝑝(𝑡)𝑎1𝑡,𝑆1(𝑡)𝑆1(𝑡)𝑑𝑡+𝑖>1𝜂𝑖(𝑡)𝑏𝑖𝑡,𝑆𝑖(𝑡)𝑑𝑊𝑖(𝑡)(A.10) which is reduced to:𝑑𝑝(𝑡)=𝜆1𝑝(𝑡)𝑑𝑡+𝑖>1𝜂𝑖(𝑎𝑡)𝑖𝑡,𝑆𝑖(𝑡)𝜆1𝑆𝑖(𝑡)𝑑𝑡+𝑖>1𝜂𝑖(𝑡)𝑏𝑖𝑡,𝑆𝑖(𝑡)𝑑𝑊𝑖(𝑡).(A.11) We need to find an expression for 𝑝(𝑡) that is function only of 𝑝(𝑡) and 𝑆𝑖(𝑡); more specifically as function of 𝑝(𝑡) and the parameters of 𝑆𝑖(𝑡). The first step is to find 𝜂𝑖(𝑡) as function of 𝑝(𝑡) and 𝑆𝑖(𝑡) or 𝑊𝑖(𝑡). This will be obtained through the generalized Clark-Ocone formula. To do this, we need to do some mathematical manipulations such as the change of the probability measure. This will enable us to find an expression for 𝜂𝑖(𝑡) as function of 𝑝(𝑡) and 𝑆𝑖(𝑡).

Let 𝑣𝑖(𝜂𝑖,𝑆𝑖,𝑡)=𝜂𝑖(𝑡)𝑏𝑖(𝑡,𝑆𝑖(𝑡)), 𝑆(𝑡)=[𝑆2(𝑡),,𝑆𝑑(𝑡)]𝑇, 𝜂(𝑡)=[𝜂2(𝑡),,𝜂𝑑(𝑡)]𝑇, 𝑣(𝑡)=[𝑣2(𝑡),,𝑣𝑑(𝑡)]𝑇, 𝑊(𝑡)=[𝑊2(𝑡),,𝑊𝑑(𝑡)]𝑇, a (𝑑1)-dimensional Brownian motion, 𝑡, is the accompanying filtration, 𝜃(𝑡)=[𝜃2(𝑡),,𝜃𝑑(𝑡)]𝑇 is a (𝑑1)-dimensional adapted process and 0𝑡𝑇.

Define 𝑊𝑗(𝑡)=𝑡0𝜃𝑗(𝑢)𝑑𝑢+𝑊𝑗(𝑡),𝑗=2,,𝑑,(A.12) that is, 𝑑𝑊𝑗(𝑡)=𝜃𝑗(𝑡)𝑑𝑡+𝑑𝑊𝑗(𝑡),(A.13)𝑍(𝑡)=exp𝑡0𝜃𝑇(𝑢)𝑑𝑊(1𝑢)2𝑡0𝜃𝑇(𝑢)𝜃(,𝑢)𝑑𝑢(A.14) and define the new probability measure𝑑𝑄=𝑍(𝑇)𝑑𝒫,(A.15a) where 𝒫 is the old probability measure and 𝑄 is the new probability measure.

Notice that, for a random variable 𝑥, 𝐸𝑄𝑥𝑡=1𝐸𝑍(𝑡)𝒫𝑍(𝑇)𝑥𝑡=1𝐸𝑍(𝑡)𝑍(𝑇)𝑥𝑡,(A.15b)where 𝐸𝑄[] the expected value with respect to the new probability measure Q, and 𝐸[] is the expected value with respect to the old probability measure 𝒫.

Girsanov's theorem states that 𝑊𝑗(𝑡) is a Wiener process with respect to the probability measure 𝑄. In addition 𝑊𝑗(𝑡) is a 𝑡 martingale with respect to 𝑄 [37, 43]. We choose 𝜃𝑖(𝑡) such that the drift coefficient is function only of 𝑝(𝑡). Other options or shapes could also be used. Thus, let𝜃𝑖𝑎(𝑡)=𝑖𝑡,𝑆𝑖(𝑡)𝜆1𝑆𝑖(𝑡)𝑏𝑖𝑡,𝑆𝑖(𝑡),𝑖>1.(A.16a) Notice that 𝜕𝜃𝑖(𝑡)𝜕𝑆𝑗(𝑡)=0,𝑖𝑗,𝑖,𝑗>1.(A.16b)Substitute (A.13), (A.16a), and (A.16b) into (A.11), we get:𝑑𝑝(𝑡)=𝜆1𝑝(𝑡)𝑑𝑡+𝑣𝑇(𝑊𝑡)𝑑(𝑡),(A.17) this has the form𝑑𝑝(𝑡)=𝜌(𝑡)𝑝(𝑡)𝑑𝑡+𝑣𝑇(𝑊𝑡)𝑑(𝑡),(A.18) where 𝜌(𝑡)=𝜆1.(A.19)

A.1. The Estimates for the Stochastic Time-Varying Coefficients 𝜂𝑖(𝑡) [30, 34]

The generalized version of the Clark-Ocone formula will now be used to find an estimate for the unknown time-varying coefficients 𝜂𝑖(𝑡) as function of 𝑝(𝑡) and 𝑆𝑖(𝑡).

Let 𝑈(𝑡)=𝑒𝑡0𝜌(𝑠)𝑑𝑠𝑝(𝑡)=𝑒𝜆1𝑡𝑝(𝑡), applying Ito lemma, we get𝑑𝑈(𝑡)=𝑒𝜆1𝑡𝑣𝑇(𝑊𝑡)𝑑(𝑡).(A.20) Or, equivalently, 𝑒𝜆1𝑡𝑝(𝑇)=𝑝(0)+𝑇0𝑒𝜆1𝑡𝑣𝑇𝑊(𝑡)𝑑(𝑡).(A.21) Define the random variable 𝐹=𝐹(𝑇)=𝑒𝜆1𝑇𝑝(𝑇)=𝑝(0)+𝑇0𝑒𝜆1𝑡𝑣𝑇𝑊(𝑡)𝑑(𝑡).(A.22) We shall now use the generalized Clark-Ocone formula [30] that represents the random variable 𝐹 as:𝐹=𝐸𝑄[𝐹]+𝑇0𝐸𝑄𝐷𝑡𝐹𝐹𝑇𝑡𝐷𝑡𝜃𝑊(𝑠)𝑑(𝑠)𝑡𝑑𝑊(𝑡),(A.23) where 𝐷𝑡𝐹 is the Malliavin derivative of 𝐹 [34] and it is a (𝑑1)×1 vector, 𝐷𝑡𝜃(𝑠) is the Malliavin derivative of 𝜃(𝑠) and it is a (𝑑1)×(𝑑1) diagonal matrix since (𝜕𝜃𝑖(𝑡))/(𝜕𝑆𝑗(𝑡))=0,𝑖𝑗. Equating (A.22) and (A.23) we get𝑇0𝑒𝜆1𝑡𝑣𝑇𝜂,𝑆𝑑𝑊,𝑡(𝑡)=𝑇0𝐸𝑄𝐷𝑡𝐹𝐹𝑇𝑡𝐷𝑡𝜃𝑊(𝑠)𝑑(𝑠)𝑡𝑑𝑊(𝑡).(A.24) Equating the coefficients of 𝑑𝑊(𝑡) on both sides, we get𝑒𝜆1𝑡𝑣𝑇𝜂,𝑆,𝑡=𝐸𝑄𝐷𝑡𝐹𝐹𝑇𝑡𝐷𝑡𝜃𝑊(𝑠)𝑑(𝑠)𝑡,(A.25a) hence,𝑣𝑇𝜂,𝑆,𝑡=𝑒𝜆1𝑡𝐸𝑄𝐷𝑡𝐹𝐹𝑇𝑡𝐷𝑡𝜃𝑊(𝑠)𝑑(𝑠)𝑡.(A.25b)Remember that 𝑣𝑖(𝜂𝑖,𝑆𝑖,𝑡)=𝜂𝑖(𝑡)𝑏𝑖(𝑡,𝑆𝑖(𝑡)), 𝜂(𝑡)=[𝜂2(𝑡),,𝜂𝑑(𝑡)]𝑇, and 𝑣(𝑡)=[𝑣2(𝑡),,𝑣𝑑(𝑡)]𝑇. After some manipulations, we get (𝑑1) equations in the (𝑑1) unknowns [𝜂2(𝑡)𝜂𝑑(𝑡)]. Using the probability measure 𝑄, each equation is of the form: 𝜂𝑖𝑒(𝑡)=𝜆1𝑡𝑏𝑖𝑡,𝑆𝑖𝐸(𝑡)𝑄𝐷𝑖𝑡𝐹𝐹𝑇𝑡𝜕𝜃𝑖/(𝑠)𝜕𝑆𝑖𝐷(𝑠)𝑡𝑆𝑖𝑊(𝑠)𝑑𝑖(𝑠)𝑡.(A.26) Or, using the probability measure 𝒫𝜂𝑖𝑒(𝑡)=𝜆1𝑡𝑏𝑖𝑡,𝑆𝑖𝐸𝐷(𝑡)𝑍(𝑡)𝑍(𝑇)𝑖𝑡𝐹𝐹𝑇𝑡𝜕𝜃𝑖/(𝑠)𝜕𝑆𝑖𝐷(𝑠)𝑡𝑆𝑖𝑊(𝑠)𝑑𝑖(𝑠)𝑡,(A.27) where 𝜃𝑖(𝑡)=(𝑎𝑖(𝑡,𝑆𝑖(𝑡))𝜆1𝑆𝑖(𝑡))/(𝑏𝑖(𝑡,𝑆𝑖(𝑡))).

In our case, 𝑎𝑖(𝑡,𝑆𝑖(𝑡))=𝑐𝑖(𝑂𝑖(𝑡)𝑆𝑖(𝑡)), and 𝑏𝑖(𝑡,𝑆𝑖(𝑡))=𝑒𝑖. Substituting these values and setting 𝜆1=0, we get𝑒𝑖𝜂𝑖(𝑡)=𝐸𝑄𝐷𝑖𝑡𝐹𝐹𝑇𝑡𝜕𝜃𝑖/(𝑠)𝜕𝑆𝑖𝐷(𝑠)𝑡𝑆𝑖𝑊(𝑠)𝑑𝑖(𝑠)𝑡.(A.28)

Equation (A.28) when solved yields an estimate for the unknown time-varying coefficient 𝜂𝑖(𝑡). This requires a knowledge of the shapes of 𝐷𝑡𝐹, 𝑎𝑖(𝑡,𝑆𝑖(𝑡)), 𝑏𝑖(𝑡,𝑆𝑖(𝑡)), 𝐷𝑡𝑆𝑖(𝑠), and 𝑆𝑖(𝑡). Since a mathematical model for 𝑆𝑖(𝑡) is available, it is easy to get the desired values.

A.2. An Approximate Expression for 𝐸𝑄{[𝐷𝑖𝑡𝐹(𝑇)𝐹(𝑇)𝑇𝑡𝐷𝑖𝑡𝜃𝑖𝑊(𝑠)𝑑𝑖(𝑠)]/𝑡}

Most of the time, it is difficult to find a closed-form expression for this equation. Instead we shall try to put it as a summation of three terms: (1) constant, (2) integration w.r.t. time, and (3) Integration w.r.t. a Wiener process. This way we will be able, in our case, to find a closed-form solution to this formula. We shall present an approximate solution to the above formula.

Since 𝜂𝑖(𝑡) is in general function of 𝑊(𝑡), it could be represented by the SDE 𝑑𝜂𝑖(𝑡)=𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗(𝑡).(A.29) In our analysis we are more concerned with 𝑒𝑖𝜂𝑖(𝑡). Thus we need an SDE for 𝑒𝑖𝜂𝑖(𝑡). This has the form𝑑𝑒𝑖𝜂𝑖(𝑡)=𝑒𝑖𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑒𝑖𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗(𝑡).(A.30) Let 𝑓𝑖𝑊,𝑡=𝐸𝑄𝐷𝑖𝑡𝐹(𝑇)𝐹(𝑇)𝑇𝑡𝐷𝑖𝑡𝜃𝑖𝑊(𝑠)𝑑𝑖(𝑠)𝑡.(A.31) Using the generalized Clark-Ocone formula for the vector case and the fact that “𝑑𝑊𝑗(𝑡)=𝜃𝑗(𝑡)𝑑𝑡+𝑑𝑊𝑗(𝑡),” one obtains an expression for the random variable 𝑓𝑖(𝑊,𝑡) as: 𝑓𝑖𝑊,𝑡=𝐸𝑄𝑓𝑖𝑊+,𝑡𝑡0𝑗𝜙𝑖𝑗𝑠,𝑊(𝑑𝑊𝑠)𝑗(𝑠)=𝐸𝑄𝑓𝑖𝑊+,𝑡𝑡0𝑗𝜙𝑖𝑗𝑠,𝑊𝜃(𝑠)𝑗𝑠,𝑊(𝑠)𝑑𝑠+𝑡0𝑗𝜙𝑖𝑗𝑠,𝑊(𝑠)𝑑𝑊𝑗(𝑠),(A.32) where, from the generalized Clark-Ocone formula,𝜙𝑖𝑗𝑢,𝑊(𝑢)=𝐸𝑄𝐷𝑗𝑢𝑓𝑖𝑊,𝑡𝑓𝑖𝑊,𝑡𝑡𝑢𝐷𝑗𝑢𝜃𝑗𝑢,𝑊(𝑑𝑊𝑢)𝑗(𝑢)𝑢.(A.33) In (A.32), the term “𝐸𝑄{𝑓𝑖(𝑊,𝑡)}” could be approximated as a constant value. Thus, (A.32) or the random variable 𝑓𝑖(𝑊,𝑡) could be approximated as𝑓𝑖𝑊,𝑡Constant+Integrationw.r.t.time+Integrationw.r.t.theWienerprocesses.(A.34) An approximate SDE for 𝑓𝑖(𝑊,𝑡) could thus be obtained as𝑑𝑓𝑖𝑊,𝑡𝑗𝜙𝑖𝑗𝑡,𝑊(𝜃𝑡)𝑗𝑡,𝑊(𝑡)𝑑𝑡+𝑗𝜙𝑖𝑗𝑡,𝑊(𝑡)𝑑𝑊𝑗(𝑡).(A.35) Since 𝑒𝑖𝜂𝑖(𝑡)=𝐸𝑄𝐷𝑖𝑡𝐹𝐹𝑇𝑡𝜕𝜃𝑖/(𝑠)𝜕𝑆𝑖𝐷(𝑠)𝑡𝑆𝑖𝑊(𝑠)𝑑𝑖(𝑠)𝑡=𝑓𝑖𝑊,𝑡.(A.36) Then 𝑑𝑒𝑖𝜂𝑖(𝑡)=𝑑𝑓𝑖𝑊,𝑡.(A.37) Substituting the different expressions for 𝑑[𝑒𝑖𝜂𝑖(𝑡)] and for 𝑑𝑓𝑖(𝑊,𝑡) and equating them, we get𝑒𝑖𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑒𝑖𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗(𝑡)=𝑗𝜙𝑖𝑗𝑡,𝑊(𝜃𝑡)𝑗𝑡,𝑊(+𝑡)𝑑𝑡𝑗𝜙𝑖𝑗𝑡,𝑊(𝑡)𝑑𝑊𝑗(𝑡),(A.38) that is, 𝑒𝑖𝜇𝑖𝑊(=𝑡),𝑡𝑗𝜙𝑖𝑗𝑡,𝑊(𝜃𝑡)𝑗𝑡,𝑊(,𝑒𝑡)𝑖𝜅𝑖𝑗𝑊(𝑡),𝑡=𝜙𝑖𝑗𝑡,𝑊.(𝑡)(A.39)

Notice that the drift and the diffusion coefficients of 𝑑[𝑒𝑖𝜂𝑖(𝑡)] are related to each other. This came from the approximate solution of the CO formula.

Without an expression for 𝜙𝑖𝑗(𝑡,𝑊(𝑡)), we would not be able to solve for the unknowns 𝜂𝑖(𝑡). Instead of an exact expression, which is difficult to obtain, we shall assume a polynomial shape with unknown parameters 𝛿𝑖𝑗, 𝛾𝑘𝑖𝑗, and 𝜆𝑙𝑖𝑗, namely,𝜙𝑖𝑗𝑡,𝑊(𝑡)=𝛿𝑖𝑗+𝑘𝛾𝑘𝑖𝑗𝑊𝑘(𝑡)+𝑘𝑙𝜆𝑖𝑗𝑘𝑙𝑊𝑘(𝑡)𝑊𝑙(𝑡)+.(A.40) The unknown parameters are estimated by minimizing the squared error, 𝑒2(𝑇), between the observed, 𝑝(𝑡), and estimated observations ̂𝑝(𝑡), wherê𝑝(𝑡)=𝑖̂𝜂𝑖(𝑡)𝑆𝑖(𝑑𝑒𝑡),(A.41)𝑖̂𝜂𝑖(𝑡)=𝑒𝑖𝜇𝑖𝑊(𝑡),𝑡𝑑𝑡+𝑒𝑖𝑗𝜅𝑖𝑗𝑊(𝑡),𝑡𝑑𝑊𝑗(𝑒𝑡),(A.42)𝑖𝜇𝑖𝑊=(𝑡),𝑡𝑗𝜙𝑖𝑗𝑡,𝑊𝜃(𝑡)𝑗𝑡,𝑊,𝑒(𝑡)(A.43)𝑖𝜅𝑖𝑗𝑊(=𝜙𝑡),𝑡𝑖𝑗𝑡,𝑊(,𝜙𝑡)(A.44)𝑖𝑗𝑡,𝑊=̂𝛿(𝑡)𝑖𝑗+𝑘̂𝛾𝑘𝑖𝑗𝑊𝑘(𝑡)+𝑘𝑙̂𝜆𝑖𝑗𝑘𝑙𝑊𝑘(𝑡)𝑊𝑙𝑒(𝑡),(A.45)2(𝑇)=𝑇0[]𝑝(𝑡)̂𝑝(𝑡)2𝑑𝑡.(A.46)

Abbreviations

AIC:Akaike information criterion
AP:Arterial pressure
BP:Blood pressure
CO:Clark-Ocone
CV:Cardiovascular
FIR:Finite impulse response
OLS:Ordinary least square
PAP:Peripheral arterial pressure
SDE:Stochastic differential equation
SNR:Signal-to-noise ratio
SNRE:Signal-to-noise ratio of the estimate
TPR:Total peripheral resistance.