Importing the required packages

using CirculatorySystemModels
using CirculatorySystemModels.ModelingToolkit
using CirculatorySystemModels.DifferentialEquations
using Plots
using DisplayAs

A simple single-chamber model

Single chamber, closed-loop, lumped parameter model of the systemic circulation and the left ventricle. The circuit equivalent formulation of the model is depicted, with the pressures of each compartment, as well as most of the mechanical parameters. The model describes three compartments: the left ventricular, arterial and venous compartments. 𝑃𝑡ℎ is the intrathoracic pressure, 𝑃𝑙𝑣 is the left ventricular pressure and 𝐸𝑙𝑣(𝑡) indicates the left ventricular elastance function.

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.

Jupyter Notebook

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
0.006

Inital Pressure (mean cardiac filling pressure)

MCFP = 7.0
7.0

Calculating the additional k parameter

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} 0 =& - out_{+}p\left( t \right) + in_{+}p\left( t \right) \\ p\left( t \right) =& in_{+}p\left( t \right) \\ 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) \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) \\ 0 = LV_{+}in_{+}p\left( t \right) - LV_{+}out_{+}p\left( t \right) \\ LV_{+}p\left( t \right) = LV_{+}in_{+}p\left( t \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) \\ 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:

states(circ_sys)
3-element Vector{Any}:
 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) \\ LV_{+}out_{+}p\left( t \right) =& LV_{+}p\left( t \right) \\ MV_+{\Delta}p\left( t \right) =& LV_{+}p\left( t \right) - Csv_{+}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_{+}in_{+}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) \\ LV_{+}in_{+}p\left( t \right) =& LV_{+}out_{+}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) \\ Rs_{+}q\left( t \right) =& \frac{ - Rs_+{\Delta}p\left( t \right)}{Rs_{+}R} \\ Csa_{+}out_{+}p\left( t \right) =& Csa_{+}in_{+}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 non-existant states (like LV.p in this case), which allows for exchanging the compliances or the ventricle for one that's defined in terms of $dp/dt$).

u0 = [
        LV.p => MCFP
        LV.V => MCFP/Eₘᵢₙ
        Csa.p => MCFP
        Csa.V => MCFP*C_sa
        Csv.p => MCFP
        Csv.V => MCFP*C_sv
        ]
6-element Vector{Pair{Symbolics.Num, Float64}}:
  LV₊p(t) => 7.0
  LV₊V(t) => 233.33333333333334
 Csa₊p(t) => 7.0
 Csa₊V(t) => 7.909999999999999
 Csv₊p(t) => 7.0
 Csv₊V(t) => 77.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}:
 233.33333333333334
   7.909999999999999
   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
Example block output

This page was generated using Literate.jl.