# data[].attributes.body
1) How many equations do we need for a system?
2) How do we make good/proper assumptions?
3) What assumptions can we make to solve these equations?
4) What phenomena do these equations model?
5) How do we interpret the role of each variable and understand its graphical meaning?
6) How do we formulate models from real world problems?
7) What are the BME applications of differential equations?
8) How do we create a differential equation from experimental data?

# data[].attributes.body
Connecting types of differential equations and BME applications. We can envision building a big table summarizing these connections.

# data[].attributes.body
These equations are used in NMR and MRI to describe the 3-D changes in nuclear magnetization (Mx, My, & Mz).

# data[].attributes.headline
Pharmacokinetic Differential Equation in DCE MRI

# data[].attributes.subject
Pharmacokinetic Differential Equation in DCE MRI

# data[].attributes.body
Common pharmacokinetic model that is used in DCE MRI that describes how a contrast agent perfuses into and away from a tissue of interest (commonly a tumor). This is commonly referred to as: The Extended Tofts Model.

# data[].attributes.body
\rho\frac{d\bar{u}}{dt}=-\nabla p+\mu\nabla ^2\bar{\mu}+\rho\bar{F}

# data[].attributes.body
\dot{u} = D\triangle{u}

# data[].attributes.body
write in latex, then highlight, click on pi symbol

# data[].attributes.headline
The Standard Solid Model (Bone mechanics)

# data[].attributes.subject
The Standard Solid Model (Bone mechanics)

# data[].attributes.body
V=\frac{V_{max} S}{K_m+S}

# data[].attributes.headline
OpenSim Musculoskeletal Model Dynamics Equations

# data[].attributes.subject
OpenSim Musculoskeletal Model Dynamics Equations

# data[].attributes.body
Equations driving muscle contraction and activation dynamics given a specific set of activation's (a) and external forces (F).

# data[].attributes.body
Nonlinear reaction-diffusion equation, where u is the population size wrt both a spatial and a temporal component. D is the diffusion constant, r is the growth rate, and K is the carrying capacity.

# data[].attributes.headline
Analysis of Pressure Distribution in Transfemoral Prosthetic Socket for Prefabrication Evaluation via the Finite Element Method

# data[].attributes.subject
Analysis of Pressure Distribution in Transfemoral Prosthetic Socket for Prefabrication Evaluation via the Finite Element Method

# data[].attributes.body
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6956391/

# data[].attributes.body
These equations model tumor cell growth, diffusion, and traction forces.
https://www.nature.com/articles/s41598-020-68324-4

# data[].attributes.attachment_caption
Useful in a lot of applications (mechanics, optimization, etc) especially when there are restraints. 

# data[].attributes.attachment_link.readable_text
Introduction

The Gray-Scott system is a reaction-diffusion system. This means that it models a process that consists of a reaction and diffusion. In the case of the Gray-Scott model that reaction is a chemical reaction between two substances \(U\) and \(V\), both of which diffuse over time. During the reaction \(U\) gets used up, while \(V\) is produced. The densities of the substances \(u\) and \(v\) are represented in the simulation.

The system is characterised by two parameters: \(F\) is the rate at which \(U\) is replenished, and \(k\) controls the rate at which \(V\) is removed from the system. Varying these parameters leads to a wide range of interesting patterns, some of which look quite familiar.

Overview

The Gray-Scott system models the chemical reaction \begin{equation} U+2V \rightarrow 3V\ . \\ \end{equation} This reaction consumes \(U\) and produces \(V\). Therefore, the amount of both substances needs to be controlled to maintain the reaction. This is done by adding \(U\) at the "feed rate" \(F\) and removing \(V\) at the "kill rate" \(k\). The removal of \(V\) can also be described by another chemical reaction: \begin{equation} V\rightarrow P\ .\\ \end{equation} For this reaction \(P\) is an inert product, meaning it doesn't react and therefore does not contribute to our observations [1]. In this case the Parameter \(k\) controls the rate of the second reaction.

Both substances diffuse over time at the diffusion rates \(D_u\) and \(D_v\) . In this simulation \(D_u\) is double \(D_v\) . This is simply to ensure that patterns can actually form, as they also depend on the diffusion rates.

You can control the feed rate \(F\) and the kill rate \(k\). Most possible combinations of the two will not lead to any patterns, which is why we have given you a number of presets that lead to interesting results. Simply choose one of these presets and observe what happens. Of course, you can also experiment yourself by using the sliders for \(F\) and \(k\).

To the top

Simulation

Select a grid size:

These buttons can be used to change the colour-code:

To the top

Description

What can be controlled are the parameters feed rate \(F\) and kill rate \(k\) to see different patterns.

Feel free to choose from the presets to see some previously defined pattern-formations.

Colours can be changed either to black and white or with the Colour-Sliders R/G/B.

The speed can be controlled (for grid sizes up to 150) by controlling how many times the calculation is done before updating the image. The pattern-formation is changed by this as well.

With the random button some random points are chosen where the concentration of substance \(V\) is initially set to 1 to start the reaction-diffusion.

To the top

Theory

The Gray-Scott system is defined by two equations that describe the behaviour of two reacting substances: \begin{align} \dot{u}&=-uv^2+F(1-u)+D_u\Delta u\\ \dot{v}&=uv^2-(F+k)v+D_v\Delta v \ .\\ \end{align} The variables in these equations \(u\) and \(v\) are the concentrations of the two reacting substances \(U\) and \(V\). On the left-hand side of each equation is the time derivative of one of these concentrations, describing the rate at which it changes. The right-hand sides of the equations both contain three separate terms . The first describes the reaction between the two substances. Since one \(U\) and two \(V\) react, the corresponding term includes \(u\) to the power of one and \(v\) to the power of two: \(uv^2\). As \(U\) gets consumed by the reaction, the term has a negative sign in the first equation. In the second equation it has a positive sign, as \(V\) is produced in the reaction.

The second term of the first equation describes the rate at which \(U\) is replenished externally. This is necessary, as \(U\) would otherwise simply be used up. The feed rate is given by the parameter \(F\). \(F\) is multiplied by \((1-u)\) to ensure that \(u\) is replenished at a rate dependent on the current concentration, which never exceeds \(1\). \(V\) does not need to be replenished, since it is produced in the reaction. Instead, it needs to be removed in order maintain the reaction. The rate of removal, the "kill rate", is controlled by the parameter \(k\). To remove \(V\) faster than \(U\) is added, \(k\) is added to \(F\) and multiplied by \(v\), since the removal of \(V\) is also supposed to be dependent on its concentration.

The last term in both equations describes the diffusion of \(U\) and \(V\), respectively.

Phase diagram of the Gray-Scott system.[1]

As you can see in the simulation above, there is a variety of patterns that can form for different values for \(F\) and \(k\). However, most possible combinations of values for those parameters lead to rather boring results: For combinations of \(F\) and \(k\) on one side of a certain threshold the screen is completely white, indicating that it is completely filled with substance \(U\). There is no substance \(V\) left to react with, and so the system is in equilibrium. The system will always return to this state.

For combinations of \(F\) and \(k\) on the other side of this threshold the screen completely fills with black. All of the added substance \(U\) reacts immediately, so the area is almost completely filled with substance \(V\). In this case the system is also in equilibrium, and will always return to this state.

The region of interest of the simulation is close to this threshold. The system changes from one stable state to the other. The system potentially has three fixed points: a trivial fixed point at \(u=1\) and \(v=0\), and two non-trivial fixed points . Which ones of these actually exists depends on the parameters \(F\) and \(k\). The trivial fixed point always exists, since the reaction cannot take place without \(v\). The two non-trivial fixed points exist beyond a bifurcation line at \(F=4(F+k)^2\) (the solid line in the figure). One of these is always a saddle, meaning it is stable in one direction and unstable in the other. The other one changes its stability at the dotted line. [2] These changes in number and stability of fixed points are what creates the patterns.

To the top

Example Pictures

Some snapshots of the patterns after they filled the whole simulation-frame.

Flower-pattern: \(F\)=0.055 ; \(k\)=0.062

Mazes-pattern: \(F\)=0.029 ; \(k\)=0.057

Mitosis-pattern: \(F\)=0.028 ; \(k\)=0.062

Solitons-pattern: \(F\)=0.03 ; \(k\)=0.06

References

[1] John E. Pearson: "Complex Patterns in a Simple System", Science Vol. 261

[2] Claudius Gros: "Adaptive and Dynamical Systems", Springer 4th Edition

To the top

# data[].attributes.attachment_link.title
Gray-Scott Model of a Reaction-Diffusion System

# data[].attributes.attachment_link.description
The Gray-Scott system is a reaction-diffusion system. This means that it models a process that consists of a reaction and diffusion. In the case of the Gray-Scott model that reaction is a chemical reaction between two substances \(U\) and \(V\), both of which diffuse over time.

# data[].attributes.attachment_link.readable_text
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

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).

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

Generic solution for the two-element Windkessel model. $$ \begin{align} p(t) = p_0 \, e^{-\frac{(t-t_{0})}{\tau}} + \frac{1}{C} \, \int_{t_0}^{t} Q(t') \, e^{\frac{t'-t}{\tau}} \, dt' \tag{6.23} \end{align} $$

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

Impedance of two element Windkessel. $$ \begin{align} Z_n^{WK} &= \frac{R}{1+j\omega_n R C} \tag{6.31} \end{align} $$

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}

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

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.

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

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}

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.

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).

# data[].attributes.attachment_link.keywords[]
principal invariants of the stress tensor

# data[].attributes.attachment_link.keywords[]
characteristic equation of the stress tensor