Windkessel Modelling

1. Windkessel modelling
- 1.1. Matlab to Python
- 1.2. 3 and 4 Element Windkessel Implementation
2. Aortic Windkessel Models - according to cite:Stergiopulos1999a
3. Bibliography

1 Windkessel modelling

1.1 Matlab to Python

Ian has programmed the 3-Element Windkessel in Matlab. Let’s see if we can do that in Python.

This is the Matlab code:

%% Generates 3 element Windkessel forward data
%  Synthesizes t, p(t), p'(t), i(t), di/dt(t) uniformly-spaced
%  samples resembling experimental data

clear all
close all

%% globals

global t_fwd i_fwd didt_fwd

%% Parameters of forward data

% Table 2 Piroli; Line 1 brachiocephalic artery
R1_fwd        = 6.30e7;                % proximal
R2_fwd        = 1.71e9;                % distal
C_fwd         = 10.1e-10;
amp           = 1.3*1.5e-5;                   % current amplitude
dt_fwd        = 0.01;
t_fwd         = [0:dt_fwd:10];
noise         = 0.0;
min_i         = 0.1*amp;V
%% generate noisy p, p' signals

i_fwd          = (amp-min_i) * (1 - cos(2 * 3.1419 / 1 * t_fwd ).*cos(2 * 3.1419 / 1 * t_fwd ))+min_i;
chop           = find ( sin(2 * 3.1419 / 1 * t_fwd ) < 0 );
i_fwd(chop)    = min_i;
didt_fwd       = [ diff(i_fwd),0] / dt_fwd;
IC             = [0 0];                 % initial conditions p, p'
%tspan         = [1 20];                % sets timespan using full vector, not range
tspan          = t_fwd;                 % sets timespan using full vector, not range
options        = odeset('Reltol',1e-3); % options
[t, WM3E_fwd]  = ode45(@(t,WM3E) RHS_defn(t,WM3E,R1_fwd,R2_fwd,C_fwd), tspan, IC, options);
wk             = size(WM3E_fwd(:,1));
WM3E_fwd(:,1)  = WM3E_fwd(:,1)+rand(wk(1),1)*(max(WM3E_fwd(:,1))-min(WM3E_fwd(:,1)))*noise;
WM3E_fwd(:,2)  = WM3E_fwd(:,2)+rand(wk(1),1)*(max(WM3E_fwd(:,2))-min(WM3E_fwd(:,2)))*noise;

%% Plot and store
figure(1)
plot(t,WM3E_fwd(:,1),t,WM3E_fwd(:,2),'LineWidth',2);
xlabel('Time');
ylabel('p'); %labels y-axis
title('Time series plots p and pp ')
legend(['p'],['pp'])

figure(2)
plot(t,i_fwd,'LineWidth',2);
xlabel('Time');
ylabel('i'); %labels y-axis
title('Time series plots i ')
legend(['i'])

clearvars amp chop dt IC noise options t tspan wk IC WM3E;
save("New_forward_data.mat");

%% Functions

function [ d_WM3E ] = RHS_defn(t,WM3E,R1,R2,C)

global t_fwd i_fwd didt_fwd

d_WM3E     = zeros(2,1);
location   = min(find(t_fwd>=t));
i          = i_fwd    (location);
didt       = didt_fwd (location);
d_WM3E(1)  = R1*didt + i/C - WM3E(1)/R2/C;
d_WM3E(2)  =               - WM3E(1)/R2/C + i/C;

end

I’ll split the Python code into chunks, for didactic reasons:.

1.1.1 Input Waveform, generic pulse

First we need to load the necessary libraries:

We need numpy and scipy, which has the functions for ODEs, and since we want to plot, we also need matplotlib.

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

# Seaborn makes the plots look a bit more modern
import seaborn as sns

sns.set_theme()
sns.set(rc={"figure.dpi": 300, "savefig.dpi": 300})
sns.set_style("white")

Then we set the parameters for the WK model and for the input waveform:

# Generic Input Waveform
# half cos^2 wave for systole and flat for diastole

# max value
max_i = 1.0
# min value
min_i = 0.1 * max_i
# systole ratio
systRatio = 2
# end time and time step for array
end_time = 10
dt = 0.01

Define the input waveform as a half cosine wave

The equation for the waveform is:

\[ I = (I_{max}-I_{min}) \left(1- \cos^{2}(2\pi t) \right) + I_{min} \]

Implemented as:

# create linear matrix for time
t = np.linspace(0, 10, num=int(end_time / dt))

# basic cos^2 wave
i = (max_i - min_i) * (1 - pow(np.cos(2 * np.pi / (2 / systRatio) * t), 2)) + min_i

# plot curve
plt.xlabel("time [s]")
plt.ylabel("")
plt.xlim(8, 10)
plt.plot(t, i)
plt.grid()

Figure 1: \(cos^{2}\) waveform is the basis for the generic pulse waveform

We want the waveform to be more realistic and have a systolic and a diastolic phase. So we cut every second wave off. For this we need to replace the values for every second wave with the minimum flow.

The mask is defined by the modulo of \(t \% 1\) and all values for which the modulo is greater than 0.5, will be replaced with the minimum flow rate.

# replace all values where t % 1 is outside the systolic time with min_i
i = np.where( t % 1 < (1/systRatio),
        (max_i - min_i) * (1 - pow(np.cos(2 * np.pi / (2/systRatio) * t), 2)) + min_i,
        min_i
        )

# plot curves
plt.xlabel("time [s]")
plt.ylabel("")
plt.xlim(8,10)
plt.plot(t,i);
plt.grid()

Figure 2: \(cos^{2}\) waveform with diastolic flatline

This looks a lot more like a blood flow waveform with distinct systole and diastole. We’ll use this as our generic input waveform.

1.1.2 Physiological maximum and minimum flow rates

Based on our CFD model realistic values for the common carotid artery (CCA) are:

\(V_{max}\ [\mathrm{m/s}]\)	\(V_{min}\ [\mathrm{m/s}]\)	\(A\ [\mathrm{m^2}]\)
1	0.2	3e-5

So we’ll use those for max_i and min_i:

# Generic Input Waveform
# half cos^2 wave for systole and flat for diastole

# max volume flow in m^3/s
max_i = 1.0 * 3e-5

# min volume flow in m^3/s
min_i = 0.2 * 3e-5

# minimum current value
# end time and time step for array
end_time = 10
dt = 0.01

# create linear matrix for time
t = np.linspace(0, 10, num=int(end_time / dt))

# i_fwd waveform
i = np.where( t % 1 < (1/systRatio),
        (max_i - min_i) * (1 - pow(np.cos(2 * np.pi / (2/systRatio) * t), 2)) + min_i,
        min_i
        )

didt = np.where( t % 1 < (1/systRatio),
            4
            * np.pi / (2/systRatio)
            * (max_i - min_i)
            * np.sin(2 * np.pi / (2/systRatio) * t)
            * np.cos(2 * np.pi / (2/systRatio) * t),
            0
        )
# plot curves
plt.xlabel("time [s]")
plt.ylabel("volume flow [$\\mathrm{m^{3}/s}$]")
plt.xlim(8, 9)
plt.plot(t, i)
plt.plot(t, didt)
plt.grid()

1.1.3 Windkessel model implementation

The 3-element windkessel model can be written as on ODE as follows (this seems to be the widely accepted form, but I have not yet found the original reference):

\[ \left(1+\frac{R_{1}}{R_{2}}\right)I(t) + C R_{1} \frac{dI(t)}{dt} = \frac{P(t)}{R_{2}} + C \frac{dP(t)}{dt} \]

We want to solve the ODE for \(\frac{dP}{dt}\), so re-arrange this equation:

\[ \frac{dP(t)}{dt} = \frac{1}{C} \left(1+\frac{R_{1}}{R_{2}}\right)I(t) + R_{1} \frac{dI(t)}{dt} - \frac{P(t)}{C R_{2}} \]

Note:

In the end I did not use the arrays for I as defined above, but used the equations to define the functions. solve_inv seems to work better with functions than with arrays.

We can define the functions for \(I(t)\) and \(dI(t)/dt\), and then the whole ODE, parameters from [Pirola2017]

# Table 2 Pirola; Line 2 Left Common Carotid Artery
# R1 = 17.60e7
# R2 = 41.7e8
# C = 4.1e-10
# Table 2 Pirola; Line 1 BrachioCephalic Artery
R1 = 6.30e7
R2 = 17.1e8
C  = 10.1e-10
L = 1e7


def I(t):
    # This needs to return a scalar, not an array.
    # Otherwise it won't work in the windkessel solver
    if t % 1 < (1/systRatio):
        I = (max_i - min_i) * (1 - pow(np.cos(2 * np.pi / (2/systRatio) * t), 2)) + min_i
    else:
        I = min_i

    return I


def didt(t):
    if t % 1 < (1/systRatio):
        didt = (
            4
            * np.pi / (2/systRatio)
            * (max_i - min_i)
            * np.sin(2 * np.pi / (2/systRatio) * t)
            * np.cos(2 * np.pi / (2/systRatio) * t)
        )
    else:
        didt = 0

    return didt


def WK3(t, p):
    dp = 1 / C * (1 + R1 / R2) * I(t) + R1 * didt(t) - p / (C * R2)
    return dp

Using the solve_ivp function in scipy we can solve this (note the from ... import construct, which allows to import functions to be used without their namespace):

from scipy.integrate import solve_ivp

# Solve WK3 for time 0-30 seconds, initial value 0
# for values along linear time array with 1000 values
pressure = solve_ivp(WK3, (0, 30), [0], t_eval=np.linspace(0, 30, 10000), method="RK45")
# the solution is returned as an array with labels 't' and 'y'
# where 'y' is a vector array (in this case only 1-d)

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(0, 30)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure["t"], pressure["y"][0])
plt.grid()

You see that the Runge-Kutta algorithm in 5^th(4^th) order formulation (RK45) shows an instability. Let’s try another method (RK23):

# Solve WK3 for time 0-30 seconds, initial value 0
# for values along linear time array with 1000 values
# using the 3rd(2nd) order Runge-Kutta scheme
pressure = solve_ivp(WK3, (0,30), [12000], t_eval=np.linspace(0,30,10000), method='RK23')
# the solution is returned as an array with labels 't' and 'y'
# where 'y' is a vector array (in this case only 1-d)

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(0,30)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure['t'],pressure['y'][0]);
plt.grid()

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(28, 30)
plt.ylim(19000, 24000)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure["t"], pressure["y"][0])
plt.grid()

Save the pressure for comparison with alternative implementation:

pressure_gen_orig = pressure

1.2 3 and 4 Element Windkessel Implementation

In this I’m going to include the 4-Element Windkessel model as well as the \(p, p'\) formulation.

The 4-Element WK can be described using the equation (solved for \(\frac{dP(t)}{dt}, \frac{dP'(t)}{dt}\)):

\[ \frac{dP(t)}{dt} = \frac{1}{C} \left(1+\frac{R_{1}}{R_{2}}\right)I(t) + R_{1} \frac{dI(t)}{dt} - \frac{P(t)}{C R_{2}} + \frac{L}{R_{2}C} \frac{dI(t)}{dt} + L \frac{d^{2}I(t)}{dt^{2}} \]

1.2.1 TODO Check the following equation with Ian

\[ \frac{dP'(t)}{dt} = \frac{1}{C} \left(1+\frac{R_{1}}{R_{2}}\right)I(t) - \frac{P(t)}{C R_{2}} + \frac{L}{R_{2}C} \frac{dI(t)}{dt} + L \frac{d^{2}I(t)}{dt^{2}} \]

1.2.2 Pedrizzetti Waveform

This is a simple Fourier series approximation of a realistic wave form. Parameters and formulation from [Pedrizzetti1996].

# Parameters
C0 = 0.4355
C2 = 0.05
C4 = -0.13
C6 = -0.1
C8 = -0.01

S2 = 0.25
S4 = 0.13
S6 = -0.02
S8 = -0.03

# Area of the CCA
# area = 3e-5
area = 600


def I(t):
    return (
        C0
        + C2 * np.cos(2 * np.pi * t)
        + S2 * np.sin(2 * np.pi * t)
        + C4 * np.cos(4 * np.pi * t)
        + S4 * np.sin(4 * np.pi * t)
        + C6 * np.cos(6 * np.pi * t)
        + S6 * np.sin(6 * np.pi * t)
        + C8 * np.cos(8 * np.pi * t)
        + S8 * np.sin(8 * np.pi * t)
    ) * area -100


# def didt(t):
#     return (
#         0
#         - C2 * 2 * np.pi * np.sin(2 * np.pi * t)
#         + S2 * 2 * np.pi * np.cos(2 * np.pi * t)
#         - C4 * 2 * np.pi * np.sin(4 * np.pi * t)
#         + S4 * 2 * np.pi * np.cos(4 * np.pi * t)
#         - C6 * 2 * np.pi * np.sin(6 * np.pi * t)
#         + S6 * 2 * np.pi * np.cos(6 * np.pi * t)
#         - C8 * 2 * np.pi * np.sin(8 * np.pi * t)
#         + S8 * 2 * np.pi * np.cos(8 * np.pi * t)
#     ) * area


# def d2idt2(t):
#     return (
#         0
#         - C2 * 4 * np.pi ** 2 * np.cos(2 * np.pi * t)
#         - S2 * 4 * np.pi ** 2 * np.sin(2 * np.pi * t)
#         - C4 * 4 * np.pi ** 2 * np.cos(4 * np.pi * t)
#         - S4 * 4 * np.pi ** 2 * np.sin(4 * np.pi * t)
#         - C6 * 4 * np.pi ** 2 * np.cos(6 * np.pi * t)
#         - S6 * 4 * np.pi ** 2 * np.sin(6 * np.pi * t)
#         - C8 * 4 * np.pi ** 2 * np.cos(8 * np.pi * t)
#         - S8 * 4 * np.pi ** 2 * np.sin(8 * np.pi * t)
#     ) * area

The inlet wave form then looks like this:

# plot curves
plt.xlabel("time [s]")
plt.ylabel("volume flow [$\\mathrm{m^{3}/s}$]")
plt.xlim(28, 30)

t = np.linspace(0, 30, 10000)
plt.plot(t, I(t))
plt.grid()

And solving the windkessel model. See the SciPy Documentation for details.

def WK3(t, p):
    dp = np.zeros(2)

    # The usual WK equation
    dp[0] = R1 * didt(t) + (1 + R1 / R2) * I(t) / C - p[0] / (R2 * C)
    dp[1] = -p[0] / (R2 * C) + (1 + R1 / R2) * I(t) / C

    return dp


def WK4(t, p):
    dp = np.zeros(2)
    # The 4-Element WK
    dp[0] = R1 * didt(t) + (1 + R1 / R2) * I(t) / C - p[0] / (R2 * C) + L * 1 / ( R2 * C ) * didt(t) + L * d2idt2(t)
    # This is most likely incorrect:
    dp[1] = -p[0] / (R2 * C) + (1 + R1 / R2) * I(t) / C + L * 1 / ( R2 * C ) * didt(t)
    return dp

# Solve WK3 for time 0-30 seconds, initial value 0
# for values along linear time array with 1000 values

R1 = 6.30e7
R2 = 17.1e8
C = 10.1e-10
L = 5e6


pressure3 = solve_ivp(
    WK3, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
)
pressure4 = solve_ivp(
    WK4, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
)
# the solution is returned as an array with labels 't' and 'y'
# where 'y' is a vector array (in this case only 1-d)

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(0, 30)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure3["t"], pressure3["y"][0], label="p, WK3")
plt.plot(pressure3["t"], pressure3["y"][1], label="p', WK3")
plt.plot(pressure4["t"], pressure4["y"][0], label="p, WK4")
plt.plot(pressure4["t"], pressure4["y"][1], label="p', WK4")
plt.legend()
plt.grid()


NameErrorTraceback (most recent call last)
<ipython-input-54-0692a19adbeb> in <module>
     11     WK3, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
     12 )
---> 13 pressure4 = solve_ivp(
     14     WK4, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
     15 )

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/ivp.py in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, **options)
    475         method = METHODS[method]
    476 
--> 477     solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
    478 
    479     if t_eval is None:

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/rk.py in __init__(self, fun, t0, y0, t_bound, max_step, rtol, atol, vectorized, first_step, **extraneous)
     98         self.max_step = validate_max_step(max_step)
     99         self.rtol, self.atol = validate_tol(rtol, atol, self.n)
--> 100         self.f = self.fun(self.t, self.y)
    101         if first_step is None:
    102             self.h_abs = select_initial_step(

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/base.py in fun(t, y)
    137         def fun(t, y):
    138             self.nfev += 1
--> 139             return self.fun_single(t, y)
    140 
    141         self.fun = fun

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/base.py in fun_wrapped(t, y)
     19 
     20     def fun_wrapped(t, y):
---> 21         return np.asarray(fun(t, y), dtype=dtype)
     22 
     23     return fun_wrapped, y0

<ipython-input-53-034bd14175b9> in WK4(t, p)
     12     dp = np.zeros(2)
     13     # The 4-Element WK
---> 14     dp[0] = R1 * didt(t) + (1 + R1 / R2) * I(t) / C - p[0] / (R2 * C) + L * 1 / ( R2 * C ) * didt(t) + L * d2idt2(t)
     15     # This is most likely incorrect:
     16     dp[1] = -p[0] / (R2 * C) + (1 + R1 / R2) * I(t) / C + L * 1 / ( R2 * C ) * didt(t)

NameError: name 'd2idt2' is not defined

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(28, 30)
plt.ylim(18500, 25500)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure3["t"], pressure3["y"][0], label="p, WK3")
plt.plot(pressure3["t"], pressure3["y"][1], label="p', WK3")
plt.plot(pressure4["t"], pressure4["y"][0], label="p, WK4")
plt.plot(pressure4["t"], pressure4["y"][1], label="p', WK4")

plt.legend()
plt.grid()


NameErrorTraceback (most recent call last)
<ipython-input-55-499b7fb98c6b> in <module>
      9 plt.plot(pressure3["t"], pressure3["y"][0], label="p, WK3")
     10 plt.plot(pressure3["t"], pressure3["y"][1], label="p', WK3")
---> 11 plt.plot(pressure4["t"], pressure4["y"][0], label="p, WK4")
     12 plt.plot(pressure4["t"], pressure4["y"][1], label="p', WK4")
     13 

NameError: name 'pressure4' is not defined

1.2.3 Generic Waveform

Comparison of the alternative implementation for the generic waveform.

systRatio = 2
def I(t):
    # This needs to return a scalar, not an array.
    # Otherwise it won't work in the windkessel solver
    if t % 1 < (1/systRatio):
        I = (max_i - min_i) * (1 - pow(np.cos(2 * np.pi / (2/systRatio) * t), 2)) + min_i
    else:
        I = min_i

    return I


def didt(t):
    if t % 1 < (1/systRatio):
        didt = (
            4
            * np.pi / (2/systRatio)
            * (max_i - min_i)
            * np.sin(2 * np.pi / (2/systRatio) * t)
            * np.cos(2 * np.pi / (2/systRatio) * t)
        )
    else:
        didt = 0

    return didt

# Solve WK3 for time 0-30 seconds, initial value 0
# for values along linear time array with 1000 values
pressure3 = sp.integrate.solve_ivp(
    WK3, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
)
pressure4 = sp.integrate.solve_ivp(
    WK4, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
)
# the solution is returned as an array with labels 't' and 'y'

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(0, 30)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure3["t"], pressure3["y"][0], label="p WK3")
plt.plot(pressure3["t"], pressure3["y"][1], label="p' WK3")
plt.plot(pressure4["t"], pressure4["y"][0], label="p, WK4")
plt.plot(pressure4["t"], pressure4["y"][1], label="p', WK4")
plt.legend()
plt.grid()


NameErrorTraceback (most recent call last)
<ipython-input-57-df6e0a3ac019> in <module>
      4     WK3, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
      5 )
----> 6 pressure4 = sp.integrate.solve_ivp(
      7     WK4, (0, 30), [0, 0], t_eval=np.linspace(0, 30, 10000), method="RK23"
      8 )

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/ivp.py in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, **options)
    475         method = METHODS[method]
    476 
--> 477     solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
    478 
    479     if t_eval is None:

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/rk.py in __init__(self, fun, t0, y0, t_bound, max_step, rtol, atol, vectorized, first_step, **extraneous)
     98         self.max_step = validate_max_step(max_step)
     99         self.rtol, self.atol = validate_tol(rtol, atol, self.n)
--> 100         self.f = self.fun(self.t, self.y)
    101         if first_step is None:
    102             self.h_abs = select_initial_step(

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/base.py in fun(t, y)
    137         def fun(t, y):
    138             self.nfev += 1
--> 139             return self.fun_single(t, y)
    140 
    141         self.fun = fun

/usr/lib/python3/dist-packages/scipy/integrate/_ivp/base.py in fun_wrapped(t, y)
     19 
     20     def fun_wrapped(t, y):
---> 21         return np.asarray(fun(t, y), dtype=dtype)
     22 
     23     return fun_wrapped, y0

<ipython-input-53-034bd14175b9> in WK4(t, p)
     12     dp = np.zeros(2)
     13     # The 4-Element WK
---> 14     dp[0] = R1 * didt(t) + (1 + R1 / R2) * I(t) / C - p[0] / (R2 * C) + L * 1 / ( R2 * C ) * didt(t) + L * d2idt2(t)
     15     # This is most likely incorrect:
     16     dp[1] = -p[0] / (R2 * C) + (1 + R1 / R2) * I(t) / C + L * 1 / ( R2 * C ) * didt(t)

NameError: name 'd2idt2' is not defined

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{Pa}$]")
plt.xlim(28, 30)
plt.ylim(19000, 25000)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure3["t"], pressure3["y"][0], label="p WK3")
plt.plot(pressure3["t"], pressure3["y"][1], label="p' WK3")
plt.plot(pressure4.t, pressure4.y[0], label="p, WK4")
plt.plot(pressure4.t, pressure4.y[1], label="p', WK4")

plt.legend()
plt.grid()


NameErrorTraceback (most recent call last)
<ipython-input-58-52c212eca7cc> in <module>
      9 plt.plot(pressure3["t"], pressure3["y"][0], label="p WK3")
     10 plt.plot(pressure3["t"], pressure3["y"][1], label="p' WK3")
---> 11 plt.plot(pressure4.t, pressure4.y[0], label="p, WK4")
     12 plt.plot(pressure4.t, pressure4.y[1], label="p', WK4")
     13 

NameError: name 'pressure4' is not defined

pressure3

2 Aortic Windkessel Models - according to [Stergiopulos1999a]

We are also interested in the aortic windkessel behaviour. Key differences will be:

No minimum flow rate in diastole. Diastolic flow rate will be zero.
Larger influence of proximal values (valve resistance) and compliance.

Values here are taken from [Stergiopulos1999a], Fig 4, table 2. Non SI units: mm_Hg, ml/s.

# Generic Input Waveform
# max volume flow in ml/s - type C, fig 2, Stergiopulos 1999
max_i = 425

# min volume flow in m^3/s
min_i = 0.0

# Period time in s
T = 0.8

# Syst. Time in s
systTime = 2 / 5 * T

# Dicrotic notch time
dicrTime = 0.02

# minimum current value
# end time and time step for array
end_time = 10
dt = 0.01

Define the waveform as functions for the WK models:

def I(t):
    # implicit conditional using boolean multiplicator
    # sin waveform
    I = (max_i - min_i) * np.sin(np.pi / systTime * (t % T)) * (
        t % T < (systTime + dicrTime)
    ) + min_i
    # sin^2 waveform
    # I = (max_i - min_i) * np.sin(np.pi / systTime * (t % T)) ** 2 * (
    #     t % T < systTime

    return I


# time derivative needed for WK model
# def didt(t):
#     didt = (
#         (np.pi / systTime)
#         * (max_i - min_i)
#         * np.cos(np.pi / systTime * (t % T))
#         * (t % T < (systTime + dicrTime))
#     )

#     return didt

def didt(t):
    return (I(t + 1e-6) - I(t - 1e-6)) / 2e-6


def d2idt2(t):
    return (didt(t + 0.5e-5) - didt(t - 0.5e-5)) / 1e-5



# def d2idt2(t):
#     d2idt2 = (
#         (np.pi ** 2 * (max_i - min_i) / (systTime ** 2))
#         ,* (- np.sin(np.pi / systTime * (t % T)))
#         ,* (t % T < (systTime + dicrTime))
#     )

t = 0
ta = np.linspace(0, 30, 30000)

plt.plot(ta, I(ta), label="$\dot{V}$")
plt.plot(ta, 0.1 * didt(ta), label="$d\dot{V}/dt$")
plt.plot(ta, 0.01 * d2idt2(ta), label="$d^{2}\dot{V}/dt^{2}$")
plt.xlabel("time [$\\mathrm{s}$]")
plt.ylabel("$\dot{V}\ [\mathrm{\\frac{ml}{s}}],\ \\frac{d\dot{V}}{dt}\ [10\ \mathrm{\\frac{ml}{s^{2}}}],\ \\frac{d^{2}\dot{V}}{dt^{2}}\ [100\  \mathrm{\\frac{ml}{s^{3}}}]$ ")
plt.xlim(0, 3)
plt.ylim(-800,800)
plt.grid()
plt.legend(loc='lower center', ncol=3)

def WK3(t, p):
    dp = np.zeros(2)

    # The usual WK equation
    dp[0] = R1 * didt(t) + (1 + R1 / R2) * I(t) / C - p[0] / (R2 * C)
    dp[1] = -p[1] / (R2 * C) + I(t) / C

    return dp


def WK4(t, p):
    dp = np.zeros(2)
    # The 4-Element WK
    dp[0] = (
        R1 * didt(t)
        + (1 + R1 / R2) * I(t) / C
        - p[0] / (R2 * C)
        + L / (R2 * C) * didt(t)
        + L * d2idt2(t)
    )
    # This is most likely incorrect:
    dp[1] = -p[1] / (R2 * C) + I(t) / C
    return dp

def WK3v(t, p):
    dp = np.zeros(2)

    # The usual WK equation
    dp[0] = R1 * didt(t) + (1 + R1 / R2) * I(t) / C - p[0] / (R2 * C)
    dp[1] = -p[0] / (R2 * C) + (1 + R1 / R2) * I(t) / C

    return dp


def WK4v(t, p):
    dp = np.zeros(2)
    # The 4-Element WK
    dp[0] = (
        R1 * didt(t)
        + (1 + R1 / R2) * I(t) / C
        - p[0] / (R2 * C)
        + L / (R2 * C) * didt(t)
        + L * d2idt2(t)
    )
    # This is most likely incorrect:
    dp[1] = -p[0] / (R2 * C) + (1 + R1 / R2) * I(t) / C + L / (R2 * C) * didt(t)
    return dp

# Solve WK3 for time 0-30 seconds, initial value 0
# for values along linear time array with 1000 values


# R1_3 = R1_4
# R2_3 = R2_4
# C_3 = C_4

# R1 = 0.03
# R2 = 0.63
# C = 5.16

R1 = 0.05
R2 = 0.9
C = 1.0666


pressure3 = solve_ivp(
    WK3,
    (0, 30),
    [0, 0],
    t_eval=np.linspace(0, 30, 30000),
    method="RK45",
    rtol=1e-6,
    vectorized=True,
)

R1 = 0.125
C = 2
L = 0.001

pressure4 = solve_ivp(
    WK4,
    (0, 30),
    [0, 0],
    t_eval=np.linspace(0, 30, 30000),
    method="RK45",
    rtol=1e-6,
    vectorized=True,
)
# the solution is returned as an array with labels 't' and 'y'
# where 'y' is a vector array (in this case only 1-d)

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{mm_{Hg}}$]")
plt.xlim(0, 30)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure3["t"], pressure3["y"][0], label="p, WK3")
plt.plot(pressure3["t"], pressure3["y"][1], label="p', WK3")
plt.plot(pressure4["t"], pressure4["y"][0], label="p, WK4")
plt.plot(pressure4["t"], pressure4["y"][1], label="p', WK4")
plt.legend()
plt.grid()

# plot curves
plt.xlabel("time [s]")
plt.ylabel("pressure [$\\mathrm{mm_{Hg}}$]")
plt.xlim(28, 30)
plt.ylim(60, 140)

# plot first component (index 0) of array 'y'
# over array 't':
plt.plot(pressure3["t"], pressure3["y"][0], label="p WK3")
plt.plot(pressure3["t"], pressure3["y"][1], label="p' WK3")
plt.plot(pressure4.t, pressure4.y[0], label="p, WK4")
plt.plot(pressure4.t, pressure4.y[1], label="p', WK4")

plt.legend()
plt.grid()

3 Bibliography

[Pedrizzetti1996] Pedrizzetti, Unsteady Tube Flow over an Expansion, Journal of Fluid Mechanics 310, 89-111 (1996).
[Pirola2017] Pirola, Cheng, Jarral, O'Regan, Pepper, Athanasiou & Xu, On the Choice of Outlet Boundary Conditions for Patient-Specific Analysis of Aortic Flow Using Computational Fluid Dynamics, Journal of Biomechanics 60, 15-21 (2017).
[Stergiopulos1999a] Stergiopulos, Westerhof & Westerhof, Total Arterial Inertance as the Fourth Element of the Windkessel Model, American Journal of Physiology-Heart and Circulatory Physiology 276(1), H81-H88 (1999).

(1 > 2) * 2