Time Domain Estimation of Arterial Parameters using the Windkessel Model and the Monte Carlo Method

Numerous parameter estimation techniques exist for characterizing the arterial system using electrical circuit analogs. However, they are often limited by their requirements and usually high computational burdain. Therefore, a new method for estimating arterial parameters based on Monte Carlo simulation is proposed. A three element Windkessel model was used to represent the arterial system. The approach was to reduce the error between the calculated and physiological aortic pressure by randomly generating arterial parameter values, while keeping constant the arterial resistance. This last value was obtained for each subject using the arterial flow, and was a necessary consideration in order to obtain a unique set of values for the arterial compliance and peripheral resistance. The estimation technique was applied to in vivo data containing steady beats in mongrel dogs, and it reliably estimated Windkessel arterial parameters. Further, this method appears to be computationally efficient for on-line time-domain estimation of these parameters.


Introduction
With the intent of practical use, this work aims to estimate the parameters of the arterial system, at a low calculation cost, prioritizing speed and simplicity. We chose the Windkessel model, a lumped model [1,2], over tube models [3][4][5] and anatomically based distributed models [6][7][8], because it provides a quantitative point of view of the main parameters of a vascular system. Therefore it can be used for the systemic arterial system and the pulmonary arterial bed of all mammals.
The three-element Windkessel models the arterial system represented as analogous to an electric circuit, as shown in Figure 1. From this model we deduced the two parameters that we sought to estimate, the peripheral resistance R p and the peripheral compliance C p . There are numerous widely used and accepted methods for the estimation of these parameters [9]. The most common ones are: the decay time method [10], the stroke volume over pulse pressure method [11][12][13], the area method [14], the two-area method [15], the pulse pressure method [16,17], the transient method [18].

Materials and Methods
Measurements of pressures from the left ventricle and the aorta arterial were taken to a group of mongrel dogs. Aortic (P ao ) and left ventricular (P lv ) pressures of the Windkessel model were registered with a sampling rate of 100 Hz during 4 seconds per measurement. Simulations were performed in Scipy [31] using a computer with an Intel Core i5 2.53 GHz processor, with 4 GB DDR3 RAM.

Windkessel Model
The three-element Windkessel model is a way to characterize part of the arterial system as an electrical circuit, the electrical diagram is shown in figure 1. The three-element Windkessel employs equations derived from the analysis of the electrical circuit. Below are the final forms for the calculation of aortic pressure. These equations require three parameters and the left ventricular pressure as input data in order to simulate an aortic pressure curve. The equation for the simulation of aortic pressure P ao during systole is: where: and And the equation for the simulation of aortic pressure P ao during diastole is:

Monte Carlo Parameter Estimation and Curve fitting
The Windkessel model used in this case generates a P ao with three inputs, which are the model's parameters. Generating a P ao vector requires only the three input parameters, and considering that the flow is not taken into account in these calculations, there are infinite solutions for the three impedances (C p , R a and R p ) that would yield the same P ao , each corresponding to a different flow.
In order to have a determined system, the aortic flow was taken into account calculating R a as a function of it. To find this relation, R a can be computed by means of equation (7), obtaining the difference between mean left ventricular pressure and mean aortic pressure during systole, divided by the measured aortic flow. This could also be approximated by using the Hagen-Poiseuille equation. This calculation, though hindering an optimal calculation of R a , allows a much more precise approximation to the real values of the peripheral impedances, C p and R p . Moreover, another alternative allows R p to remain constant, and computing R p as the difference between mean arterial pressure and mean right atrial pressure, divided by the mean aortic flow. In this case, R a could be calculated more precisely, if that is desired.
The aortic flow during systole results in: then R a is: The Monte Carlo method was used to generate pseudo-random values for the peripheral parameters R p and C p . These pseudo-random values are generated within a fixed interval limited by the physiological boundaries for these parameters [15,32]. The initial intervals we set for R p and for C p were [1,20] and [0. 1,5] respectively, selected in a way that the physiological values were included,and with some additional range in case of any anomaly.
An ammount of N sets were randomly generated and input into equations 1 and 4 along with experimental data for P lv . This resulted in a set of N simulated aortic pressure curves as a function of time. Then, the curve which best fitted the experimental aortic pressure curve was selected employing the least squares method. Also, the sixth best curve was selected and compared with the best curve using an algorithm that shortened the interval that limits the pseudo-random variables generated.
During several iterations, the original parameter intervals were subsequently shortened to obtain a determined parameter value. To do so, the following algorithm was used: The value of the parameters of the best fit curve bv were set as midpoints. Those corresponding to the sixth best fit curve 6bv determined the distance from the midpoint to the upper limit and from the midpoint to the lower limit. If the new upper limit surpassed the upper physiological boundary, the former was set as the upper limit. The same was done for the lower limit.
The process of generating N curves and narrowing the parameter intervals was repeated until the top and bottom values of the intervals were within range of desired precision, we set it to ten decimal points. The estimated optimal values for the parameters thus correspond to the last iteration of the process.

Results
The estimated parameter values obtained with a precision of 10 −4 for three different data sets are presented in Table 3. Here, the error is computed as the sum of the squares of the errors in each point and N represent the number of curves generated that allowed achieving the desired precision. Figures 2, 3, and 4 present the experimental data plotted against the best fit curve obtained, corresponding to the optimum values for R p and C p .  Figure 2. Estimated P ao (red), experimental P ao (blue) and P lv (green) for Data Set 1. Table 2 presents the computation time for all data sets varying the number of initial curves, N . Observing the calculation time for each data set, it was observed they had similar duration, with variations attributed to signal characteristics of each data set, specially the signal noise. It should be noted that, considering three heartbeats, the desired precision is achieved faster for  Figure 3. Estimated P ao (red), experimental P ao (blue) and P lv (green) for Data Set 2.

Discussion
The aortic pressure curves were approximated with a relatively low error, implying that the parameters obtained should be close to the real values. However, we were not able to take any physical measurements of the arterial parameters to compare to our results. An alternative would involve estimating this parameters with other methods, but it is important to note that most of them are useful for estimating only C p . These physical measurements also have errors themselves which would be something to be compared with much attention. Our data also lacked the mean flow, which forced us to calculate the arterial resistance, R a , using numerical methods. By creating an error function, which resulted being unimodal, we found the R a which provided the minimum of the error function for several trials of our method.
Using the Monte Carlo method implies that results are subject to probability and could provide varying results. By increasing the number of generated pseudo-random values, N , this variation is reduced. Moreover, P ao curves with noise and higher frequency components required a higher N value, and therefore, a higher computation time to reach a certain precision, as evidenced by the fact that the third data set (Data Set 3) required higher N values to reach the same precision than the first two sets (Data Sets 1 and 2).
It was also observed that when N is too large, precision is lost, because of an inherent complication of our interval limiting algorithm. When too many random aortic pressure curves are generated, the best and sixth best curve end up being too close and thus drastically shortening the interval limits for the next iteration. Sometimes this shortening is so drastic that it cuts the ideal curve out and the algorithm converges to parameters outside logical boundaries.
The desired precision was achieved for different number of starting sets depending on the quality of the signal. We choose N = 1000 to avoid any extreme variations of the Monte Carlo method and the issues that arise with large N . In order to improve the performance of our method when facing noisy signals, filtering the noise frequency component could be useful, and therefore a smaller N should achieve the desired precision.
For Data Sets 1 and 2, an approximate time of 7.5 s was required for estimating the parameters of a single heart-beat, which makes this method a good alternative for estimating these parameters on a clinical environment, and using these results as a diagnostic tool.

Author Contributions
Analyzed the data: VG IP GRP GVD. Contributed reagents/materials/analysis tools: VG IP GRP GVD HMV MR. Wrote the paper: VG IP GRP GVD HMMV MR.