6.2 Lumped models
Figure 58: The heart ejects blood intermittently into the compliant aorta, like water pumped in to a Windkessel with a compliant air chamber (adopted from [13]).
Currently, the term lumped models allude to a family of mathematical models for the load to the heart, in which the physics of the entire systemic arterial tree is represented by a few, lumped parameters.
Stephen Hales (1733) presented a basic, and conceptual model of the arterial tree, based on the observation that the blood flow in the peripheral arteries is relatively smooth, despite the pulsatile action of the heart. He envisioned that the interaction between the heart and the arteries has similarities with the working principle of a fire hose, in which the pulsatile action of the pump is damped by an air chamber (Windkessel in German). In the cardiovascular system, large arteries play the role of the air chamber of a Windkessel (fig. 58)
In 1899 Otto Frank formulated these ideas mathematically, by introducing the two element Windkessel which consisted of two building blocks: the peripheral resistance, \( R \), (representing the peripheral arterioles and capillaries), and the total arterial compliance \( C \), which accounts for the compliance (or elasticity) of the larger conduit vessels.
A more accurate description of the pressure-flow relation may be obtained by including more lumped elements in a lumped model of the cardiovascular system. For example the characteristic impedance \( Z_c \) of the aorta is added for the 3-element Windkessel model, resulting in a better description of the high frequency modes of the aortic input impedance and accordingly a better prediction of the aortic pressure pulse [14]. An overview of the most commonly used lumped models for the vascular tree is given by [15] [15]. Increasing complexity of the lumped models provides better prediction of the pressure-flow relation. However, as the number of lumped elements increase, the physical meaning and hence the determination of each individual parameter becomes less clear [5]. It has been shown that the three-element Windkessel systematically overestimates arterial compliance [16] [17]. As pointed out in [7], whenever more and more complexity must be introduced into a scientific theory to make its predictions match reality, it is usually time to look for a new hypothesis.
However, in spite of the severity of the underlying assumptions, lumped parameter models can be used for various purposes:
- Estimation of the mechanical properties of the arterial circulation.
- Representation of a part of the circulatory system (e.g., the arterioles and capillaries or the part of the vascular tree not accounted for in greater detail), for more complex models such as 1D network models accounting for wave propagation phenomena or 3D FSI-models and CFD models.
- A tool to get more insight into the physical problem.
6.2.1 The Windkessel model
Figure 59: The rate of change of arterial volume equals the difference between aortic inflow (\( Q \)) and outflow (\( Q_p \)) towards the periphery.
The mathematical representation of the Windkessel model attributed to Otto Frank, is obtained by the requirement of mass conservation for a model of the vascular tree as seen in figure 59. Let \( Q \) represent the inflow to the aorta, and \( Q_p \) the flow towards the periphery, whereas \( Q_a \) represents the stored volume per time unit in the aorta. The ordinary differential equation, representing the Windkessel model is then obtained by the requirement of mass conservation: $$ \begin{align} Q &= Q_a + Q_p = \partd{V}{p} \partd{p}{t} + \frac{p}{R} \tag{6.14}\\ & = C \partd{p}{t} + \frac{p}{R} \tag{6.15} \end{align} $$
where \( C=\partial V/\partial p \) is the total arterial compliance, and \( Q_p = p/R \) and \( R \) is the total peripheral resistance. We may rewrite equation (6.15) on a conventional differential equation form:
Two element Windkessel model $$ \begin{align} \partd{p}{t} + \frac{1}{RC} \, p & = \frac{Q(t)}{C} \tag{6.16} \end{align} $$ which is the differential equation for the Windkessel model in figure 60. Mathematically we realize that equation (6.16) represents a first order, linear, temporally driven differential equation. The homogenous solution \( p_h \), to equation (6.16), i.e., a vanishing right hand side of equation (6.16), is given by: $$ \begin{align} p_h = p_0 \, e^{-\frac{(t-t_{0})}{\tau}} \tag{6.17} \end{align} $$ where \( p_0 \) is a constant depending on the initial conditions at time \( t_0 \) and \( \tau=RC \) represents a time constant. The solution \( p_h \) is sometimes referred to as the autonomous solution, as in this case there is no driving force on the right hand side of the differential equation (6.16).
Figure 60: Aortic flow (lower pane) from measurements and measured and predicted aortic pressures (upper panel). The measured pressure is represented with a solid dark blue line.
The solution of equation (6.16) with a right hand side unequal to zero, normally referred to as the inhomogenous solution or the particular solution , may found by assuming that \( p_0 = p_0(t) \), i.e., that \( p_0 \) is a function of time: $$ \begin{align} p = p_0(t) \, e^{-\frac{t}{\tau}} \tag{6.18} \end{align} $$
Differentiation of equation (6.18) yields: $$ \begin{align} \partd{p}{t} = \partd{p_0}{t} \, e^{-\frac{t}{\tau}} - \frac{1}{\tau} \,p_0(t) \, e^{-\frac{t}{\tau}} \tag{6.19} \end{align} $$ which by substitution in equation (6.16) cancels the second term and we are left with: $$ \begin{align} \partd{p_0}{t} & = \frac{Q(t)}{C} \, e^{\frac{t}{\tau}} \tag{6.20} \end{align} $$ Equation (6.20) may then readily be integrated to give: $$ \begin{align} p_0(t) = \frac{1}{C} \, \int_{t_0}^{t} Q(t') \, e^{\frac{t'}{\tau}} \, dt' \tag{6.21} \end{align} $$ and by subsequent substitution of equation (6.21) into equation (6.18) we obtain a particular solution: $$ \begin{align} p_p(t) = \frac{1}{C} \, \int_{t_0}^{t} Q(t') \, e^{\frac{t'-t}{\tau}} \, dt' \tag{6.22} \end{align} $$
Finally, we obtain the generic solution for the two-element Windkessel model, represented by equation (6.16), as the sum of the homogenous solution \( p_h \) in equation (6.17) and the particular solution \( p_p \) in equation (6.22):
We observe that the flow \( Q(t) \) occurs as a driving function in the second term of equation (6.23), and that this term has the form of a convolution integral. Therefore, for a given initial pressure \( p_0 \) at time \( t_0 \), and flow as a function of time \( Q(t) \), on may compute the pressure \( p(t) \) from equation (6.23). The integral may readily computed numerically, e.g. with trapezoidal numerical integration, whenever an analytical solution is not tractable.
However, the numerical integration in equation (6.23) may be avoided, based on the following reasoning. Normally, pressure and flow measurements are sampled at discrete time instances \( t_i = i \, \Delta t \), where \( \Delta t \) is the sample interval and consequently \( t_1 = \Delta t \). Further, we introduce the standard notation: \( p_0 = p(t_0), p_1=p(t_1), \ldots, p_i=p(t_i) \) and similarly for samples of the flow \( Q_i = Q(t_i) \). The integral in equation (6.23) may then be evaluated for \( t=t_1 \), i.e., the interval for the time integration is \( t_1 - t_0 = \Delta t \). A fair assumption to make is that \( \Delta t/\tau \ll 1 \), as the sampling rate \( \Delta t \), normally is small (typically \( 5\cdot 10^{-3} \, s \), and \( \tau \) is a relaxation scale which should be a fraction of the heart rate (approximately the period of diastole, \( \tau = RC\propto 0.5 \,s \)). Subject to these assumptions we get from equation (6.23): $$ \begin{align} p(t_1) & = p_0 \, e^{-\frac{\Delta t}{\tau}} + \frac{1}{C} \, \int_{t_0}^{t_1} Q(t') \, e^{\frac{t'-t_1}{\tau}} \, dt' \tag{6.24}\\ & = p_0 \, e^{-\frac{\Delta t}{\tau}} + \frac{Q(t^*_1)}{C} \, \int_{t_0}^{t_1} e^{\frac{t'-t_1}{\tau}} \, dt', \qquad t_0 \leq t^*_1 \leq t_1 \tag{6.25}\\ & = p_0 \, e^{-\frac{\Delta t}{\tau}} + \left . \frac{Q(t^*_1) \, \tau }{C} \, e^{\frac{t'-t_1}{\tau}} \right |^{t_1}_{t_0} = p_0 \, e^{-\frac{\Delta t}{\tau}} + \frac{Q(t^*_1) \, \tau }{C} \, \left (1- e^{\frac{-\Delta t}{\tau}} \right ) \tag{6.26}\\ & \approx p_0 \, \left (1 -\frac{\Delta t}{\tau} \right ) + \frac{Q(t^*_1) \, \Delta t}{C} \, \tag{6.27} \end{align} $$ The first identity in equation (6.24) follows from the mean value theorem (8.33), where \( Q(t^*_1) \) denote a mean value of the flow \( Q(t) \) in the sampling interval \( [t_0, t_1] \). The two following identities are analytical evaluations of the remaining integral. Finally, a Taylor expansion yields an approximation in the last expression in the equation (6.24). From, the estimation of the integral in equation (6.24), a generalized expression for \( p(t_i) \) follows naturally: $$ \begin{equation} p_i \approx p_{i-1} \, \left (1 -\frac{\Delta t}{\tau} \right ) + \frac{Q(t^*_i) \, \Delta t}{C} \, \tag{6.28} \end{equation} $$
The expression in equation (6.23), works remarkably well in correlating experimental data on pressure \( p \) and flow \( Q \), in particular during diastole (see [18] (p. 23) and cite[pp. 11, 210, 423]{mcdonald73}). Otto Frank reasoned that the decay of the diastolic pressure in the ascending aorta, when flow is zero, can be expressed by and exponential curve [12] (p. 122). Observe that this observation is in agreement with the homogenous solution of equation (6.16) given by equation (6.17). The time constant \( \tau = RC \), expresses the time needed for the pressure to decrease to 37\% of the starting pressure. The larger the resistance, the slower the blood will leave the compliant arteries. Conversely, the larger the compliance, the more blood is stored in the compliant arteries, and the longer time is needed for this blood to leave these vessels.
6.2.1.1 Non-invasive estimation of mean flow
The cross-sectional area compliance \( C_A \), may be estimated if the pulse wave velocity is measured (see equation (7.79)) along with the averaged cross-sectional area [12]. Subsequently, the volume compliance \( C \) may be estimated as \( C \approx C_A L_a \), if the length \( L_a \), of the aorta is measured. The peripheral resistance may then be calculated as: $$ \begin{equation} R = \frac{\tau}{C} \tag{6.29} \end{equation} $$ Finally, the mean flow \( \bar{Q} \) may then be found from mean pressure \( \bar{p} \) by: $$ \begin{equation} \bar{Q} = \frac{\bar{p}}{R} \tag{6.30} \end{equation} $$ The assumption that all compliance is lumped in the aorta, and implicitly assuming infinite wave speed in the aorta, introduces an error. Only after pulsatile flows could be measured, and the arterial input impedance could be determined, the shortcomings of the two-element Windkessel became clear [12] (chap 24).6.2.1.2 The impedance of the two-element Windkessel
By proposing that the pressure and the flow may be expressed as Fourier series (see equation (6.6) ) (and thus implicitly assume periodic solutions) the impedance for a two element Windkessel may be found by substation in equation (6.16):
which also may be conveniently represented on polar form (see [19]): $$ \begin{align} Z_n^{WK} &= \frac{R}{\sqrt{1+ \left ( \omega_n \tau \right )^2}} \; e^{j \theta} \quad \text{with} \quad \theta = -\arctan (\omega \tau) \tag{6.32} \end{align} $$
where \( \tau=RC \) has been introduced as a relaxation time for the Windkessel.
The impedance for the Windkessel is a complex quantity with a modulus and phase, which has been plotted for each harmonics (frequency) in figure 61. Observe that for the two element Windkessel, both modulus and phase decrease monotonous as a function of the frequency.
\subsubsection{The Windkessel model as a pure resistor or a pure capacitor}
Figure 61: A comparison of the input impedances of the two-element, three-element and four element Windkessel models (adopted from [12]).
We observe for later use, that the impedance of the Windkessel will approach a pure resistor if: $$ \begin{equation} Z_n^{WK} = \frac{R}{1+j\omega_n R C} \approx R \tag{6.33} \end{equation} $$ which is satisfied if \( RC \ll 1 \), i.e., \( R \ll 1/C \). At the same time we require that \( R \gg C \), as the impedance is supposed to be more of a resistor type than a compliant chamber type. Consequently, for a Windkessel to approach a pure resistor we require: $$ \begin{equation} Z_n^{WK} \approx R \quad \text{if} \quad C\ll R \ll 1/C \tag{6.34} \end{equation} $$
On the contrary, the Windkessel will approach a pure capacitor if: $$ \begin{equation} Z_n^{WK} = \frac{1}{1/R+j\omega_n C} \approx \frac{1}{j\omega C} \quad \text{if} \quad C \gg 1/R \tag{6.35} \end{equation} $$
6.2.2 The three-element Windkessel model
Figure 62: The three-element Windkessel model.
The three-element Windkessel model (see figure 62) is based on the two-element Windkessel model of Otto Frank, with an additional characteristic impedance \( Z_c \) [20]. The aortic characteristic impedance \( Z_c \) account for the combined effects of compliance and inertance of the very proximal aorta, and thus forms a link with the transmission line models. In order to derive the mathematical representation of the threee-element Windkessel, let \( p \) denote the pressure at the aortic root and \( p_d \) an aortic pressure at an undetermined but more distal position in the aorta (see figure 62). The pressure drop \( p-p_d \) between these two location is then given by: $$ \begin{equation} p-p_d = Z_c Q \tag{6.36} \end{equation} $$ whereas the relation between \( p_d \) and \( Q \) is given by the previously derived differential equation for the two-element Windkessel (equation (6.16)), when \( p_d \) is substituted for \( p \). Thus, $$ \begin{equation} \partd{p_d}{t} + \frac{p_d}{R C } = \frac{Q}{C} \tag{6.37} \end{equation} $$ The distal pressure \( p_d \) may then be eliminated by substitution of equation (6.36) into equation (6.37) which results in the following differential equation for the three-element Windkessel model: $$ \begin{equation} \partd{p}{t} + \frac{p}{RC} = \frac{Q}{C} \left (1 + Z_c/R \right ) + Z_c \partd{Q}{t} \tag{6.38} \end{equation} $$ The input impedance of the three-element Windkessel may be shown to be: $$ \begin{equation} Z_i = \frac{R+Zc + j \omega R Z_c C}{1 + j \omega R C} \tag{6.39} \end{equation} $$
From measurements for pressure and flow, the total peripheral resistance \( R \) is calculated from the time averaged pressure \( p_m \) and flow \( Q_m \). One the has the option to set \( R=p_m/Q_m \) or \( R+Z_c = p_m/Q_m \) [16]. Based on the input impedance from the measurements, the characteristic impedance may be estimated by averaging the high frequencies [21]. Leaving only one unknown model parameter, namely the compliance \( C \).
The remaining model parameter \( C \) is finally determined by using the measured flow as input to the three-element Windkessel model and fitting the (complete) estimated pressure response to the measured pressure.
Note one can also estimate the model parameters by fitting input impedance of the model to the measured input impedance.
The procedure can be applied in a number of different ways, depending on the calculation of \( R \) and \( Z_c \) from the measurements. The various approaches have been assessed in [16] applying a non-linear 1D network-model for the vascular tree. The overall conclusion is that methods based on the three-element Windkessel model overestimate compliance (from 15-40\%). Further, the calculation of the pressure response yields the best results, whereas the integral method is the inferior approach of the ones tested.
Windkessel models are frequently used to terminate more advanced models such as 1D network models and FSI/CFD-models. Such models are frequently discretize in the time-domain, and thus it is advantageous to discretize the Windkessels in the time-domain too. The discretization of the three-element Windkessel model is conveniently done by using and implicit time discretization of the differential equation (6.38): $$ \begin{equation} \frac{p^{n+1}-p}{\Delta t} + \frac{p^{n+1}}{RC} = \frac{f(Q)}{C} \tag{6.40} \end{equation} $$ where we for simplicity introduce the source on the right hand side as a function of the flow: $$ \begin{equation} f(Q) = Q \left (1 + Z_c/R \right ) + Z_c C \partd{Q}{t} \tag{6.41} \end{equation} $$ which yields: $$ \begin{equation} p^{n+1} = \frac{p + \Delta t/C \, f(Q)}{1 + \Delta t/\tau} \tag{6.42} \end{equation} $$ where \( \tau = RC \). This procedure holds for both two-element and three-element Windkessel models, as the only difference the \( f(Q)=Q \) in for the two-element Windkessel (see equations (6.16) and (6.38)).
6.2.3 Methods for estimation of total arterial compliance
\subsubsection{The time decay method (TDM) for estimation of total arterial compliance} The time decay method (TDM) is based on the assumption (see fig. 60) that the aortic flow is approximately zero during diastole (\( Q \approx 0 \)). From equation (6.23) we then see that for a Windkessel model the pressure reduces to: $$ \begin{equation} p(t) = p_{t=t_d} \, e^{-\frac{(t-t_{d})}{\tau}} \tag{6.43} \end{equation} $$ where \( t_d \) the time for the onset of diastole. By nature, this method is tailored to approximate the compliance of the aorta, and has thus limited applicability to other parts of the cardiovascular network due to the assumption of zero flow. To avoid high frequency disturbances, the onset of the diastole is normally not taken to be the time for the incisura in the aortic pressure (see upper panel of figure 60) , but rather the last \( 2/3^{\text{rd}} \) of the diastole.
Figure 63: Impedance absolute value \( |Z| \) and phase angle \( \angle Z \) for measurements and various methods for estimation of total arterial compliance.
From the expression for the impedance of the two element Windkessel in equation (6.32) we see that: $$ \begin{equation} Z_0^{WK} \equiv \frac{p_0}{q_0} = R \tag{6.44} \end{equation} $$ The TDM method requires the knowledge of some measure of the flow rate in the aorta (i.e., the mean flow rate \( q_0 \) with units \( [q_0] = m^3/s \) or more commonly the cardiac output \( CO \) which is expressed in \( l/min \)). From equation (6.44) we get: $$ \begin{equation} R = \frac{60 \; p_0}{1000 \; CO } \tag{6.45} \end{equation} $$ i.e., the peripheral resistance \( R \) is obtained form the averaged pressure \( p_0 \) and the cardiac output \( CO \). The only remaining unknown in the pressure given by equation (6.17) is the total arterial compliance, which may be obtained by a least square procedure (curve fitting), where the difference between the measured pressure and the predicted pressure by equation (6.17) is minimized.
6.2.4 Example 18: Estimation of total arterial compliance with the TDM method
Figure 64: The TDM method for the complete diastole (+) and for the last \( 2/3^{\text{rd}} \) (o) of the diastole (taken from [5]).
To see how the TDM method may be applied, consider some synthetic pressure and flow data, simulated with a linear transmission line model [22] [5]. In the transmission line model, non-linear convective terms and viscous losses are omitted from the governing mass and momentum equations for pulsatile flow in compliant vessels (see equations (7.24) and (7.25)). Further, viscoelasticity is not accounted for in the transmission line model.
The input for the simulation is a physiological flow wave at a heart rate of 60 beats/min. The material parameters for the model are chosen such that the total peripheral resistance \( R_{pm}=1 \text{ mmHg/(ml/s)} \) and the total compliance \( C_m=0.86 \text{ ml/mmHg} \). The simulated pressure resulting from the physiological flow wave input, is given by the solid line in figure 64.
An estimated total arterial compliance \( C_e \approx 0.84 \text{ ml/mmHg} \) $(\pm 2.3 \%)$ is obtained, using the TDM on the complete diastole. When the TDM is applied on the last \( 2/3^{\text{rd}} \) of the diastole, the estimate of the total arterial compliance reduces to \( C_e \approx 0.81 \text{ ml/mmHg} \) $(\pm 5.8 \%)$. The resulting pressures for the two TDM approaches are illustrated in figure 64.
\subsubsection{The area method (AM) for estimation of total arterial compliance}
Figure 65: Illustration of the area method (taken from [5]).
To motivate the area method[23], let \( t_1 \) and \( t_2 \) denote two moments in time during diastole, in which we assume the flow \( Q \approx 0 \). From the differential equation (6.16) representing the Windkessel model we get: $$ \begin{align} \int_{t_1}^{t_2} \partd{p}{t} + \frac{1}{RC} \, p \; dt = 0 \tag{6.46} \end{align} $$ which after integration and rearrangement under the assumption of constant \( R \) and \( C \) may be presented as: $$ \begin{align} RC= \frac{\int_{t_1}^{t_2} p \; dt}{p(t_2) - p(t_1)} \tag{6.47} \end{align} $$ We observe that the integral in the numerator of equation (6.47) is the area under the pressure curve between the two moments in time \( t_1 \) and \( t_2 \), respectively. Naturally, The name of the method is inspired from this property of the method.
As for the TDM, the peripheral resistance \( R \) is estimated from the mean pressure \( p_0 \) and the \( CO \) by equation (6.45). Thus, the total arterial compliance may be computed directly from equation (6.47). An advantage of the AM in comparison with the TDM, is that the exponentially decaying pressure curve does not have to be fitted. For the same conditions as in Example 6.2.4 Example 18: Estimation of total arterial compliance with the TDM method, the AM provides an estimate of the total arterial compliance of \( C_{AM} \approx 0.87 \text{ ml/mmHg} \) (\( \pm 1.2 \%) \). The method is also applicable when the pressure whenever the diastolic pressure is no true exponential function and the method can, according to the authors, be extended with a non-linear pressure volume relation. However, the compliance, obtained at mean pressure, corresponds well with the value obtained with a constant compliance [23]. Note that different diastolic time intervals may be chosen, and the resulting estimate will be sensitive to the specific interval. The authors suggest that the time interval between the dicrotic notch and the end diastole is chosen (as in figure 65).
\subsubsection{The pulse pressure method (PPM) for estimation of total arterial compliance}
The Pulse Pressure Method (PPM), like the previous methods in this section, provides an estimate of the total arterial compliance based on measurements of pressure and flow as a function of time at the ascending aorta.
Figure 66: Pressure from a 2WK with a total arterial compliance from the PPM (dots), and measured pressure (solid line). Illustration of the area method (taken from [5]).
The PPM originally suggested by [16], predicts pressure by an imposed flow from the two element Windkessel model (see e.g., equation (6.28)), where the resistance is taken as the ratio of mean pressure to mean flow (see equation (6.44)). However, the PPM differs from the other methods by the fact that it is only the pulse pressure (i.e., the difference between systolic(maximum pressure and diastolic pressure (minimal pressure)) which is used in the in the least square estimation procedure for estimation of the total arterial compliance. This results in a simple and stable iteration scheme, which is independent of the initial values of the parameter to be estimated. Neither does the method depend on a given time interval in the diastole. And further, the method applies for both normal and pathological cases.
Application of the PPM (approximation in the frequency domain) on synthetic data yields a total arterial compliance of \( 0.82 \text{ ml/mmHg} \), which is an underestimation of \( 4.2\% \) [5].
Exercise 4: Windkessel model
a) The two element Windkessel model (see (6.16)) is given by the con olution integral in (6.23). An approximated expression for the convolution integral is provided in {eq:394}. Show that an identical approximation may be obtained by discretizing equation (6.16) with a backward Euler scheme.
b) Outline a strategy for an evalutation of the various methods to estimate the total arterial compliance in section (6.2.1 The Windkessel model).
A backward Euler discretization of equation (6.16) is: $$ \begin{equation} \frac{p-p_{i-1}}{\Delta t} + \frac{1}{\tau} p = \frac{Q}{C} \tag{6.48} \end{equation} $$ where we have dropped subscipts for time level \( i \) for convenience, \ie \( p=p_i \)). By collecting terms of \( p \) in equation (6.48) we get: $$ \begin{equation} p \left ( 1 + \frac{\Delta t}{\tau} \right ) = p_{i-1} + \frac{Q \, \Delta t}{C} \tag{6.49} \end{equation} $$ and by division by the expression in the paranthesis on the lhs of equation (6.49) we get: $$ \begin{equation} p = p_{i-1} \frac{1}{1 + \frac{\Delta t}{\tau}} + \frac{Q \, \Delta t}{C} \frac{1}{1 + \frac{\Delta t}{\tau}} \tag{6.50} \end{equation} $$
which by Taylor expansion and neglection of higher order termes simplifies to: $$ \begin{equation} p \approx p_{i-1} \left ( 1 -\frac{\Delta t}{\tau} \right ) + \frac{Q \, \Delta t}{C} \tag{6.51} \end{equation} $$
which is identical with the approximation obtained by using the convolution integral in (6.23).