The algorithm that we have developed for the calculation of the reservoir pressure from the measured pressure waveform is now available as a Matlab script:
(reservoir pressure algorithm matlab script)
The algorithm is based upon the theory published in Aguado-Sierra J, Alastruey J, Wang J-J, Hadjiloizou N, Davies J, Parker KH (2008) Separation of the reservoir and wave pressure and velocity from measurements at an arbitrary location in arteries. Proc Inst Mech Eng [H], 222, 403-416, ISSN: 0954-4119. [pdf]. It has undergone extensive revision and testing since the publication of that paper and we felt it would be useful to make it available to other researchers so that they can try it themselves.
The algorithm is being made freely available under the conditions of the GNU project copyright licence (GNU General Public License). The software can be modified for your own use but please do not distribute copies of it except in the original form.
The variables and use of the algorithm are described in the header to the file. Here we reiterate and slightly amplify those notes.
The program is designed to calculate the reservoir pressure from the pressure measured over a single cardiac cycle. It is assumed that the first point (t=0) in the input pressure corresponds to the diastolic point, i.e. the point of minimum pressure just before the rapid increase in pressure at the start of systole. Inputting a pressure waveform starting at any other time (e.g. the peak of the R-wave on the ECG) will produce erroneous results. The input pressure should terminate during late diastole, preferably at the diastolic point for the next beat.
The format of the function call is:
[Pr,A,B,Pinf,Tn,Pn]=kreservoir_v8(P,Tb,Tn)
% inputs P - pressure starting at diastolic pressure
% Tb - duration of beat
% Tn - time of the dicrotic notch (if missing it is calculated
% found by the function; if ==0 it is done interactively
% outputs Pr - reservoir pressure waveform
% A - rate constant relating Pw to U
% B - rate constant = 1/tau
% Pinf - pressure asymptote
% Tn - detected time of the dichrotic notch
% Pn - pressure at dichrotic notch
In slightly more detail:
P is the pressure waveform for a single cardiac cycle. It should start at the diastolic point (i.e. the minimum pressure before the rapid rise of pressure at the start of systole).
Inputs:
Tb is the time of the cardiac cycle that is being analysed (in s).
Tn is the time of the dicrotic notch (in s). There are three options for this calling parameter:
1. If it is omitted, the program will try to find the dicrotic notch automatically. If the dicrotic notch is clearly discernible, the program will find the minimum of the notch. If there is no minimum it will look for a suitable point of inflection. If the automatic selection of this time, which is critical to the performance of the algorithm since it is used to define the diastolic period, one of the other options should be used.
2. If Tn = 0, the program will display a graph of the input pressure waveform and ask the operator to click on the dicrotic notch using the ginput function in Matlab. This interactive mode is convenient if only a few pressure waveforms are being analysed. Its disadvantage is that it is very difficult to repeat the calculation exactly because of the uncertainties in positioning the cross-hairs at precisely the same position.
3. If Tn is a number other than zero, the program will take this as the time of the dicrotic notch (in s). This option is convenient if the time of the dicrotic notch is already known to the user.
Outputs:
Pr is the calculated reservoir pressure. It is a vector of the same length as the input pressure.
A is a rate constant relating the wave pressure to the velocity at the root of the aorta. In the paper it is the variable a. Theoretically it is related to the input impedance of the aorta and other parameters of the system. However, it is probably best to think of it simply as an unknown constant of proportionality that is found by iteration to minimise the square error between measured pressure and the pressure calculated by the theory.
B is another rate constant that is simply the inverse of the reservoir time constant which is given by RC, the product of the total resistance to outflow from the arteries and the total compliance of the arteries. B is determined by the program by fitting an exponential curve to the pressure during diastole.
Pinf is the asymptote of the exponential decay of the pressure during diastole. Theoretically, it is the pressure at which flow through the microcirculation ceases. As discussed in the paper, this pressure is not necessarily the venous pressure but is related to the interstitial pressure surrounding the capillaries. This parameter is also determined by the program by fitting an exponential curve to the pressure during diastole.
Tn is the time of the dicrotic notch (in s). If Tn is specified as in input parameter, it is simply copied to the output. Otherwise it is the time of the dicrotic notch as determined either automatically by the program or interactively by the operator.
Pn is the pressure at the time of the dicrotic notch. It is included to that it is easy to plot this point as a check to see if the results of the program are reasonable.
The algorithm is written as a function to be imbedded into a Matlab program by the user. Its performance, therefore, is dependent, to some degree, by the calling program. We have tested it extensively on a wide variety of arterial pressure waveforms measured in many different ways. Invasive measurements at different sites in the aorta, femoral, brachial and radial arteries, as well as measurements at different locations in the coronary arteries have been analysed using the program. Non-invasive measurements, mainly from the radial artery, by applanation tonometry have also been extensively tested. We have tested the program in a limited number of cases where simultaneous measurements of the volume flow rate from the left ventricle were available, and found that the reservoir pressure calculated using the flow measurements were very similar to the results of the program using only the pressure measurements. [get some examples from Justin??] In a recent study (paper in preparation) the algorithm was used to calculate the reservoir pressure from more applanation tonometry measurements in than ?? patients and exhibited excellent robustness; it failed in only a dozen of cases. [ditto??]
The algorithm could be simplified somewhat by using explicit functions. However, these are not compatible with earlier versions of Matlab and so we have retained the rather cumbersome method of using global variables to communicate between the different functions that are called in the minimisation routine to be compatible with earlier versions of Matlab.