CirculationModels.jl

The case for Julia and reproducible lumped parameter modelling

Harry Saxton, Dr. Torsten Schenkel

Department of Engineering and Mathematics, Sheffield Hallam University

What is Lumped parameter modelling?

Simply put, we represent physiological compartments through electrical analogues. Resistors, capcitors and inductors etc

Key Advantages:

  • Offers the unique ability to examine both cardiac function and global hemodynamics within the context of a single model
  • Being able to quantify both is crucial for the effective prediction and diagnosis of cardiovascular diseases
  • Model computation time is often less due to being ODE/DAE based.
  • The relatively simple structure of the model allows most personalized simulations to be automated. Meaning the ability to embed these lumped parameter models into a clinical workflow could one day become trivial

ODE solvers

4-Element Windkessel

\[\begin{align} \frac{d p_{1}}{d t} = - \frac{R_{c}}{L_{p}} p_{1} + \left( \frac{R_{c}}{L_{p}} - \frac{1}{R_{p} C} \right) p_{2} \\ + R_{c} \frac{d I(t)}{d t} + \frac{I(t)}{C} \end{align}\]

\[\begin{align} \frac{d p_{2}}{d t} = - \frac{1}{R_{p} C} p_{2} + \frac{I(t)}{C} \end{align}\]

def wk4(t, y, I, Rc, Rp, C, Lp, dt):

    dp1dt = (
        -Rc / Lp * y[0]
        + (Rc / Lp - 1 / Rp / C) * y[1]
        + Rc * ddt(I, t, dt)
        + I(t) / C
    )

    dp2dt = -1 / Rp / C * y[1] + I(t) / C

    return [dp1dt, dp2dt]

def ddt(fun, t, dt):
    # Central differencing method
    return (fun(t + dt) - fun(t - dt)) / (2 * dt)

%time sol = sp.integrate.solve_ivp(
        lambda t, y: wk4(t, y, I, Rc, Rp, C, Lp, dt),
        (time_start, time_end),
        initial_cond,
        method="RK45",
        rtol=1e-9,
        vectorized=True,
    )

: CPU times: user 961 ms, sys: 5.47 ms, total: 967 ms
: Wall time: 967 ms
function dP = RHS_defn(t,P,Rc,Rp,C,Lp)

dP    = zeros(2,1);

dP(1) = -Rc / Lp * P(1) ...
        + (Rc / Lp - 1 / Rp / C) * P(2) ...
        + Rc * didt(t) + i(t) / C;

dP(2) = -1 / Rp / C * P(2) + i(t) / C;

end

function didt = didt(t)
dt = 1e-3;
didt = (i(t+dt) - i(t-dt)) / (2 * dt);
end

options        = odeset('Reltol',1e-9);

tic
[t, P] = ode45(@(t,P) RHS_defn(t,P,Rc,Rp,C,Ls,Lp), tspan, P0, options);
toc

Elapsed time is 0.129001 seconds.
using CirculationModels #Load the julia pkg
@named Rc = Resistor(R=R_c) #Characteristic resistance
@named Rp = Resistor(R=R_p) #Peripheral resistance
@named C = Capacitor(C=C_ ) #Capacitance
@named L = Inductance(L=L_) #Inertia
@named source = DrivenCurrent(I=1.0, fun=I sin)
# Flow into arterial network
@named ground = Ground()
circ eqs = [
connect(source.out, Rc.in, L.in)
connect(Rc.out, L.out, C.in, Rp.in)
connect(Rp.out, C.out, source.in, ground.g)]
# Here we connect the pins of the network
circ model = ODESystem(circ eqs, t)
# Declare an ode system dependent on t
circ model = compose( circ model, [L, Rc, Rp, C, source,
ground])
# Specify equations and component names
circ sys = structural simplify(circ model)
# Simplify equations
prob = ODAEProblem(circ sys, u0, (0, 30.0))
# ODAE allows us to only specify u0 the length of the
number of odes
sol = solve(prob)
#Solve the system
plot(sol, idxs=[C.in.p, C.in.q, Rc.q], tspan=(28, 30))

The Main point.

Domain experts can now engage in the modelling process and compartments can be interchanged easily.

Complex Model example

Theodosios Korakianitis and Yubing Shi. Numerical simulation of cardiovascular dynamics with healthy and diseased heart valves. Journal of biomechanics, 39(11):1964–1982, 2006

Compare to other languages

We can compare our CirculationModels.jl to other opensource implementations of the above model.

  • Language Length of code Time to simulate
    Python 410 42.4 (s)
    Matlab 449 12.0 (s)
    CellMLToolkit.jl 3 0.070 (s)
    CirculationModels.jl (23-88) 0.044 (s)
  • CirculationModels.jl CellMLToolkit.jl Matlab Python
    1x 1.6x 272.2x 963.6x

Conclusion

  • We have developed a lumped parameter modelling framework which is efficient, modular, reproducible and inclusive
  • We are the first modelling framework in which domain experts can play an active role within model development. As they do not require an in-depth knowledge of ODE/DAE systems.
  • The speed ups along with the seamless integration to the SciML system means model analysis on such complex systems is now possible and not limited to the number of parameters.
  • The range of applications of these models mean that any domain expert with an interest in any system can combine modelling and simulation with clinical work.

Thanks for listening!

  • Please reach out, would love to bridge the gap between modelling and clinical work!
  • Email: hs1454@hallam.shu.ac.uk
  • Github: H-Sax

CirculationModels.jl