# # Importing the required packages
using CirculatorySystemModels
using ModelingToolkit
using OrdinaryDiffEq
using Plots
using DisplayAs
A simple single-chamber model
This follows Bjørdalsbakke et al.
Bjørdalsbakke, N.L., Sturdy, J.T., Hose, D.R., Hellevik, L.R., 2022. Parameter estimation for closed-loop lumped parameter models of the systemic circulation using synthetic data. Mathematical Biosciences 343, 108731. https://doi.org/10.1016/j.mbs.2021.108731
Changes from the published version above:
- Capacitors are replaced by compliances. These are identical to capacitors, but have an additional parameter, the unstrained volume $V_0$, which allows for realistic blood volume modelling. Compliances have an inlet and an oulet in line with the flow, rather than the ground connector of the dangling capacitor.
- The aortic resistor is combined with the valve (diode) in the
ResistorDiode
element.
Define the parameters
All the parameters are taken from table 1 of [Bjørdalsbakke2022].
Cycle time in seconds
#
τ = 0.85
0.85
Double Hill parameters for the ventricle
Eₘᵢₙ = 0.03
Eₘₐₓ = 1.5
n1LV = 1.32;
n2LV = 21.9;
Tau1fLV = 0.303 * τ;
Tau2fLV = 0.508 * τ
#
0.4318
Resistances and Compliances
R_s = 1.11
C_sa = 1.13
C_sv = 11.0
11.0
Valve parameters
Aortic valve basic
Zao = 0.033
0.033
Mitral valve basic
Rmv = 0.006
# Inital Pressure (mean cardiac filling pressure)
MCFP = 7.0
### Calculating the additional `k` parameter
7.0
The ventricle elastance is modelled as:
\[E_{l v}(t)=\left(E_{\max }-E_{\min }\right) e(t)+E_{\min }\]
where $e$ is a double-Hill function, i.e., two Hill-functions, which are multiplied by each other:
\[e(\tau)= k \times \frac{\left(\tau / \tau_1\right)^{n_1}}{1+\left(\tau / \tau_1\right)^{n_1}} \times \frac{1}{1+\left(\tau / \tau_2\right)^{n_2}}\]
and $k$ is a scaling factor to assure that $e(t)$ has a maximum of $e(t)_{max} = 1$:
\[k = \max \left(\frac{\left(\tau / \tau_1\right)^{n_1}}{1+\left(\tau / \tau_1\right)^{n_1}} \times \frac{1}{1+\left(\tau / \tau_2\right)^{n_2}} \right)^{-1}\]
nstep = 1000
t = LinRange(0, τ, nstep)
kLV = 1 / maximum((t ./ Tau1fLV).^n1LV ./ (1 .+ (t ./ Tau1fLV).^n1LV) .* 1 ./ (1 .+ (t ./ Tau2fLV).^n2LV))
1.6721792928965973
Set up the model elements
Set up time as a variable t
@variables t
\[ \begin{equation} \left[ \begin{array}{c} t \\ \end{array} \right] \end{equation} \]
Heart is modelled as a single chamber (we call it LV
for "Left Ventricle" so the model can be extended later, if required):
@named LV = DHChamber(V₀ = 0.0, Eₘₐₓ=Eₘₐₓ, Eₘᵢₙ=Eₘᵢₙ, n₁=n1LV, n₂=n2LV, τ = τ, τ₁=Tau1fLV, τ₂=Tau2fLV, k = kLV, Eshift=0.0)
\[ \begin{align} p\left( t \right) &= p_0 + \left( E_{m i n} + \frac{\left( \frac{\mathrm{rem}\left( t + \left( 1 - Eshift \right) \tau, \tau \right)}{\tau_1} \right)^{n_1} \left( - E_{m i n} + E_{m a x} \right) k}{\left( 1 + \left( \frac{\mathrm{rem}\left( t + \left( 1 - Eshift \right) \tau, \tau \right)}{\tau_1} \right)^{n_1} \right) \left( 1 + \left( \frac{\mathrm{rem}\left( t + \left( 1 - Eshift \right) \tau, \tau \right)}{\tau_2} \right)^{n_2} \right)} \right) \left( - V_0 + V\left( t \right) \right) \\ \frac{\mathrm{d} V\left( t \right)}{\mathrm{d}t} &= out_{+}q\left( t \right) + in_{+}q\left( t \right) \\ 0 &= - out_{+}p\left( t \right) + in_{+}p\left( t \right) \\ p\left( t \right) &= in_{+}p\left( t \right) \end{align} \]
The two valves are simple diodes with a small resistance (resistance is needed, since perfect diodes would connect two elastances/compliances, which will lead to unstable oscillations):
@named AV = ResistorDiode(R=Zao)
@named MV = ResistorDiode(R=Rmv)
\[ \begin{align} {\Delta}p\left( t \right) &= out_{+}p\left( t \right) - in_{+}p\left( t \right) \\ 0 &= out_{+}q\left( t \right) + in_{+}q\left( t \right) \\ q\left( t \right) &= in_{+}q\left( t \right) \\ q\left( t \right) &= \frac{ - {\Delta}p\left( t \right) \left( {\Delta}p\left( t \right) < 0 \right)}{R} \end{align} \]
The main components of the circuit are 1 resistor Rs
and two compliances for systemic arteries Csa
, and systemic veins Csv
(names are arbitrary). Note: one of the compliances is defined in terms of $dV/dt$ using the option inV = true
. The other without that option is in $dp/dt$.
@named Rs = Resistor(R=R_s)
@named Csa = Compliance(C=C_sa)
@named Csv = Compliance(C=C_sv, inP=true)
\[ \begin{align} 0 &= - out_{+}p\left( t \right) + in_{+}p\left( t \right) \\ p\left( t \right) &= in_{+}p\left( t \right) \\ V\left( t \right) &= V_0 + C p\left( t \right) \\ \frac{\mathrm{d} p\left( t \right)}{\mathrm{d}t} &= \frac{out_{+}q\left( t \right) + in_{+}q\left( t \right)}{C} \end{align} \]
Build the system
Connections
The system is built using the connect
function. connect
sets up the Kirchhoff laws:
- pressures are the same in all connected branches on a connector
- sum of all flow rates at a connector is zero
The resulting set of Kirchhoff equations is stored in circ_eqs
:
circ_eqs = [
connect(LV.out, AV.in)
connect(AV.out, Csa.in)
connect(Csa.out, Rs.in)
connect(Rs.out, Csv.in)
connect(Csv.out, MV.in)
connect(MV.out, LV.in)
]
\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( LV_{+}out, AV_{+}in \right) \\ \mathrm{connect}\left( AV_{+}out, Csa_{+}in \right) \\ \mathrm{connect}\left( Csa_{+}out, Rs_{+}in \right) \\ \mathrm{connect}\left( Rs_{+}out, Csv_{+}in \right) \\ \mathrm{connect}\left( Csv_{+}out, MV_{+}in \right) \\ \mathrm{connect}\left( MV_{+}out, LV_{+}in \right) \\ \end{array} \right] \end{equation} \]
Add the component equations
In a second step, the system of Kirchhoff equations is completed by the component equations (both ODEs and AEs), resulting in the full, overdefined ODE set circ_model
.
Note: we do this in two steps.
@named _circ_model = ODESystem(circ_eqs, t)
@named circ_model = compose(_circ_model,
[LV, AV, MV, Rs, Csa, Csv])
\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( LV_{+}out, AV_{+}in \right) \\ \mathrm{connect}\left( AV_{+}out, Csa_{+}in \right) \\ \mathrm{connect}\left( Csa_{+}out, Rs_{+}in \right) \\ \mathrm{connect}\left( Rs_{+}out, Csv_{+}in \right) \\ \mathrm{connect}\left( Csv_{+}out, MV_{+}in \right) \\ \mathrm{connect}\left( MV_{+}out, LV_{+}in \right) \\ LV_{+}p\left( t \right) = LV_{+}p_0 + \left( LV_{+}E_{m i n} + \frac{\left( \frac{\mathrm{rem}\left( t + \left( 1 - LV_{+}Eshift \right) LV_+\tau, LV_+\tau \right)}{LV_+\tau_1} \right)^{LV_{+}n_1} \left( - LV_{+}E_{m i n} + LV_{+}E_{m a x} \right) LV_{+}k}{\left( 1 + \left( \frac{\mathrm{rem}\left( t + \left( 1 - LV_{+}Eshift \right) LV_+\tau, LV_+\tau \right)}{LV_+\tau_1} \right)^{LV_{+}n_1} \right) \left( 1 + \left( \frac{\mathrm{rem}\left( t + \left( 1 - LV_{+}Eshift \right) LV_+\tau, LV_+\tau \right)}{LV_+\tau_2} \right)^{LV_{+}n_2} \right)} \right) \left( - LV_{+}V_0 + LV_{+}V\left( t \right) \right) \\ \frac{\mathrm{d} LV_{+}V\left( t \right)}{\mathrm{d}t} = LV_{+}in_{+}q\left( t \right) + LV_{+}out_{+}q\left( t \right) \\ 0 = LV_{+}in_{+}p\left( t \right) - LV_{+}out_{+}p\left( t \right) \\ LV_{+}p\left( t \right) = LV_{+}in_{+}p\left( t \right) \\ AV_+{\Delta}p\left( t \right) = - AV_{+}in_{+}p\left( t \right) + AV_{+}out_{+}p\left( t \right) \\ 0 = AV_{+}in_{+}q\left( t \right) + AV_{+}out_{+}q\left( t \right) \\ AV_{+}q\left( t \right) = AV_{+}in_{+}q\left( t \right) \\ AV_{+}q\left( t \right) = \frac{ - \left( AV_+{\Delta}p\left( t \right) < 0 \right) AV_+{\Delta}p\left( t \right)}{AV_{+}R} \\ MV_+{\Delta}p\left( t \right) = MV_{+}out_{+}p\left( t \right) - MV_{+}in_{+}p\left( t \right) \\ 0 = MV_{+}in_{+}q\left( t \right) + MV_{+}out_{+}q\left( t \right) \\ MV_{+}q\left( t \right) = MV_{+}in_{+}q\left( t \right) \\ MV_{+}q\left( t \right) = \frac{ - \left( MV_+{\Delta}p\left( t \right) < 0 \right) MV_+{\Delta}p\left( t \right)}{MV_{+}R} \\ Rs_+{\Delta}p\left( t \right) = Rs_{+}out_{+}p\left( t \right) - Rs_{+}in_{+}p\left( t \right) \\ 0 = Rs_{+}in_{+}q\left( t \right) + Rs_{+}out_{+}q\left( t \right) \\ Rs_{+}q\left( t \right) = Rs_{+}in_{+}q\left( t \right) \\ Rs_+{\Delta}p\left( t \right) = - Rs_{+}R Rs_{+}q\left( t \right) \\ 0 = - Csa_{+}out_{+}p\left( t \right) + Csa_{+}in_{+}p\left( t \right) \\ Csa_{+}p\left( t \right) = Csa_{+}in_{+}p\left( t \right) \\ Csa_{+}p\left( t \right) = \frac{ - Csa_{+}V_0 + Csa_{+}V\left( t \right)}{Csa_{+}C} \\ \frac{\mathrm{d} Csa_{+}V\left( t \right)}{\mathrm{d}t} = Csa_{+}out_{+}q\left( t \right) + Csa_{+}in_{+}q\left( t \right) \\ 0 = - Csv_{+}out_{+}p\left( t \right) + Csv_{+}in_{+}p\left( t \right) \\ Csv_{+}p\left( t \right) = Csv_{+}in_{+}p\left( t \right) \\ Csv_{+}V\left( t \right) = Csv_{+}V_0 + Csv_{+}C Csv_{+}p\left( t \right) \\ \frac{\mathrm{d} Csv_{+}p\left( t \right)}{\mathrm{d}t} = \frac{Csv_{+}out_{+}q\left( t \right) + Csv_{+}in_{+}q\left( t \right)}{Csv_{+}C} \\ \end{array} \right] \end{equation} \]
Simplify the ODE system
The crucial step in any acausal modelling is the sympification and reduction of the OD(A)E system to the minimal set of equations. ModelingToolkit.jl does this in the structural_simplify
function.
circ_sys = structural_simplify(circ_model)
\[ \begin{align} \frac{\mathrm{d} LV_{+}V\left( t \right)}{\mathrm{d}t} &= - Csv_{+}out_{+}q\left( t \right) - Csa_{+}in_{+}q\left( t \right) \\ \frac{\mathrm{d} Csa_{+}V\left( t \right)}{\mathrm{d}t} &= - Csv_{+}in_{+}q\left( t \right) + Csa_{+}in_{+}q\left( t \right) \\ \frac{\mathrm{d} Csv_{+}p\left( t \right)}{\mathrm{d}t} &= \frac{Csv_{+}out_{+}q\left( t \right) + Csv_{+}in_{+}q\left( t \right)}{Csv_{+}C} \end{align} \]
circ_sys
is now the minimal system of equations. In this case it consists of 3 ODEs for the ventricular volume and the systemic and venous pressures.
Note: this reduces and optimises the ODE system. It is, therefore, not always obvious, which states it will use and which it will drop. We can use the states
and observed
function to check this. It is recommended to do this, since small changes can reorder states, observables, and parameters.
States in the system are now:
unknowns(circ_sys)
3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
LV₊V(t)
Csa₊V(t)
Csv₊p(t)
Observed variables - the system will drop these from the ODE system that is solved, but it keeps all the algebraic equations needed to calculate them in the system object, as well as the ODEProblem
and solution object - are:
observed(circ_sys)
\[ \begin{align} LV_{+}p\left( t \right) &= LV_{+}p_0 + \left( LV_{+}E_{m i n} + \frac{\left( \frac{\mathrm{rem}\left( t + \left( 1 - LV_{+}Eshift \right) LV_+\tau, LV_+\tau \right)}{LV_+\tau_1} \right)^{LV_{+}n_1} \left( - LV_{+}E_{m i n} + LV_{+}E_{m a x} \right) LV_{+}k}{\left( 1 + \left( \frac{\mathrm{rem}\left( t + \left( 1 - LV_{+}Eshift \right) LV_+\tau, LV_+\tau \right)}{LV_+\tau_1} \right)^{LV_{+}n_1} \right) \left( 1 + \left( \frac{\mathrm{rem}\left( t + \left( 1 - LV_{+}Eshift \right) LV_+\tau, LV_+\tau \right)}{LV_+\tau_2} \right)^{LV_{+}n_2} \right)} \right) \left( - LV_{+}V_0 + LV_{+}V\left( t \right) \right) \\ Csa_{+}p\left( t \right) &= \frac{ - Csa_{+}V_0 + Csa_{+}V\left( t \right)}{Csa_{+}C} \\ Csv_{+}in_{+}p\left( t \right) &= Csv_{+}p\left( t \right) \\ Csv_{+}V\left( t \right) &= Csv_{+}V_0 + Csv_{+}C Csv_{+}p\left( t \right) \\ MV_+{\Delta}p\left( t \right) &= LV_{+}p\left( t \right) - Csv_{+}p\left( t \right) \\ LV_{+}out_{+}p\left( t \right) &= LV_{+}p\left( t \right) \\ AV_{+}in_{+}p\left( t \right) &= LV_{+}p\left( t \right) \\ MV_{+}out_{+}p\left( t \right) &= LV_{+}p\left( t \right) \\ Rs_+{\Delta}p\left( t \right) &= - Csa_{+}p\left( t \right) + Csv_{+}p\left( t \right) \\ Csa_{+}out_{+}p\left( t \right) &= Csa_{+}p\left( t \right) \\ AV_{+}out_{+}p\left( t \right) &= Csa_{+}p\left( t \right) \\ Rs_{+}in_{+}p\left( t \right) &= Csa_{+}p\left( t \right) \\ Csv_{+}out_{+}p\left( t \right) &= Csv_{+}in_{+}p\left( t \right) \\ Rs_{+}out_{+}p\left( t \right) &= Csv_{+}in_{+}p\left( t \right) \\ MV_{+}in_{+}p\left( t \right) &= Csv_{+}in_{+}p\left( t \right) \\ MV_{+}q\left( t \right) &= \frac{ - \left( MV_+{\Delta}p\left( t \right) < 0 \right) MV_+{\Delta}p\left( t \right)}{MV_{+}R} \\ AV_+{\Delta}p\left( t \right) &= Csa_{+}p\left( t \right) - MV_+{\Delta}p\left( t \right) - Csv_{+}p\left( t \right) \\ LV_{+}in_{+}p\left( t \right) &= LV_{+}out_{+}p\left( t \right) \\ Rs_{+}q\left( t \right) &= \frac{ - Rs_+{\Delta}p\left( t \right)}{Rs_{+}R} \\ Csa_{+}in_{+}p\left( t \right) &= Csa_{+}out_{+}p\left( t \right) \\ Csv_{+}out_{+}q\left( t \right) &= - MV_{+}q\left( t \right) \\ MV_{+}out_{+}q\left( t \right) &= - MV_{+}q\left( t \right) \\ AV_{+}q\left( t \right) &= \frac{ - \left( AV_+{\Delta}p\left( t \right) < 0 \right) AV_+{\Delta}p\left( t \right)}{AV_{+}R} \\ Rs_{+}out_{+}q\left( t \right) &= - Rs_{+}q\left( t \right) \\ Csa_{+}out_{+}q\left( t \right) &= - Rs_{+}q\left( t \right) \\ LV_{+}in_{+}q\left( t \right) &= - Csv_{+}out_{+}q\left( t \right) \\ MV_{+}in_{+}q\left( t \right) &= - MV_{+}out_{+}q\left( t \right) \\ Csa_{+}in_{+}q\left( t \right) &= AV_{+}q\left( t \right) \\ AV_{+}out_{+}q\left( t \right) &= - AV_{+}q\left( t \right) \\ LV_{+}out_{+}q\left( t \right) &= - AV_{+}q\left( t \right) \\ Rs_{+}in_{+}q\left( t \right) &= - Rs_{+}out_{+}q\left( t \right) \\ Csv_{+}in_{+}q\left( t \right) &= - Csa_{+}out_{+}q\left( t \right) \\ AV_{+}in_{+}q\left( t \right) &= - AV_{+}out_{+}q\left( t \right) \end{align} \]
And the parameters (these could be reordered, so check these, too):
parameters(circ_sys)
18-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
LV₊V₀
LV₊p₀
LV₊Eₘᵢₙ
LV₊Eₘₐₓ
LV₊n₁
LV₊n₂
LV₊τ
LV₊τ₁
LV₊τ₂
LV₊k
LV₊Eshift
AV₊R
MV₊R
Rs₊R
Csa₊V₀
Csa₊C
Csv₊V₀
Csv₊C
Define the ODE problem
First defined initial conditions u0
and the time span for simulation:
Note: the initial conditions are defined as a parameter map, rather than a vector, since the parameter map allows for changes in order. This map can include not only state-variables, but also observables (like LV.p
in this case), which makes setting up the initial conditions easier, if the known initial conditions are not included in the resulting model states, which is often the case in MTK models. Since MTK9 these cannot be redundant, however, so we need to decide if we want to define pressures or volumes in this case.
u0 = [
LV.p => MCFP
Csa.p => MCFP
Csv.p => MCFP
]
3-element Vector{Pair{Symbolics.Num, Float64}}:
LV₊p(t) => 7.0
Csa₊p(t) => 7.0
Csv₊p(t) => 7.0
tspan = (0, 20)
(0, 20)
in this case we use the mean cardiac filling pressure as initial condition, and simulate 20 seconds.
Then we can define the problem:
prob = ODEProblem(circ_sys, u0, tspan)
ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 20)
u0: 3-element Vector{Float64}:
0.0
0.0
7.0
Simulate
The ODE problem is now in the MTK/DifferentialEquations.jl format and we can use any DifferentialEquations.jl solver to solve it:
sol = solve(prob, Vern7(), reltol=1e-12, abstol=1e-12);
Results
p1 = plot(sol, idxs=[LV.p, Csa.in.p], tspan=(16 * τ, 17 * τ), xlabel = "Time [s]", ylabel = "Pressure [mmHg]", hidexaxis = nothing) # Make a line plot
p2 = plot(sol, idxs=[LV.V], tspan=(16 * τ, 17 * τ),xlabel = "Time [s]", ylabel = "Volume [ml]", linkaxes = :all)
p3 = plot(sol, idxs=[Csa.in.q,Csv.in.q], tspan=(16 * τ, 17 * τ),xlabel = "Time [s]", ylabel = "Flow rate [ml/s]", linkaxes = :all)
p4 = plot(sol, idxs=(LV.V, LV.p), tspan=(16 * τ, 17 * τ),xlabel = "Volume [ml]", ylabel = "Pressure [mmHg]", linkaxes = :all)
img = plot(p1, p2, p3, p4; layout=@layout([a b; c d]), legend = true)
img = DisplayAs.Text(DisplayAs.PNG(img))
img
#

This page was generated using Literate.jl.