Using Mathematica to Model Cardiovascular Function

Tim Shine
Mayo Clinic
(with assistance from Roger Germundsson and Andrew Moylan, Mathematica Summer School 2009)

Abstract

The aim of this project is to develop a model of the arterial pressure wave that contains information user to patients and physicians in understanding the state of the arterial system that the wave traverses.

Articles

M.E. Green, J.W. Clark Jun., D.N. Mohr and H.M. Bourland, A Mathematical Model of Left-Ventricular Function, Medical and Biological Engineering, March 1973, pp. 126-134 [web]

M.E. Green and J.W. Clarke, The Innervated Left Ventricle: A Mathematical Model of Function, Medical and Biological Engineering, July 1973, pp. 464-468 [web]

J.W. Clark, R.C. Pruett, D.L. Baldridge, R. Srinivasan, H.M. Bourland, J.S. Cole, R.W. Brower, Functional Model for the Characterisation of the Ventricular Mechanics of the Human Subject, Medical and Biological Engineering & Computing, July, 1977, pp. 335-348 [web]

G. Avanzolini and A. Capello, Estimation of time-varying systolic properties of left ventricular mechanics, Medical and Biological Engineering and Computing, Volume 24, Number 3 / May, 1986, pp. 261-266 [web] (Refers to the Green et.al March 1973 article, but may be a simpler model.)

Models

A brief video reviewing cardiovascular function.

Heart model, left-ventricle: (Subscript[R, M]) mitral-valve resistance, (Subscript[C, V]) left-ventricle capacitance, (Subscript[R, A]) aortic valve resistance, (Subscript[P, V]) left-ventricle pressure, (Subscript[P, A]) aortic root pressure

Arterial model: L lumped inductance of the blood, (Subscript[R, B]) resistance between systems, (Subscript[R, L]) cerebral and coronary resistance, (Subscript[R, 1])(Subscript[C, L]) models visco elastic properties of ascending aorta and aortic arch, (Subscript[R, R]) systemic peripheral resistance, and (Subscript[R, 2])(Subscript[C, R]) models the visco elastic properties the descending and abdominal aorta.

Heart-model and arterial model:

Simplified models during the different phases: (a) Filling phase. (b) Isovolumic phase. (c) Ejection phase  [Note: Diodes (i.e. valves) are considered ideal ...]

Equations

The variables in the models satisfy these differential equations

Equations for Filling Phase ((Subscript[P, V])<(Subscript[P, AS])): Ep is the passive elastance
FillingPhase={
x1'[t]==((-x1[t]+(Subscript[R, R])x4[t])/((Subscript[C, R])((Subscript[R, 2])+(Subscript[R, R])))),
x2'[t]==((-(Subscript[E, V]) x2[t])/(Subscript[R, M]))+((Subscript[P, AS])/(Subscript[R, M])),
x3'[t]==((-x3[t]-(Subscript[R, L])x4[t])/((Subscript[C, L])((Subscript[R, 1])+(Subscript[R, L])))),
x4'[t]==((-x4[t])/L)((Subscript[R, B])+(((Subscript[R, 1])(Subscript[R, L]))/((Subscript[R, 1])+(Subscript[R, L])))+(((Subscript[R, 2])(Subscript[R, R]))/((Subscript[R, 2])+(Subscript[R, R]))))+(((Subscript[R, L])x3[t])/(L((Subscript[R, 1])+(Subscript[R, L]))))-(((Subscript[R, R])x1[t])/(L((Subscript[R, 2])+(Subscript[R, R]))))
};
Equations for Isovolumic Contraction and Relaxation Phase ((Subscript[P, AS])<(Subscript[P, V])<(Subscript[P, A])): 
IsovolumicPhase={
x1'[t]==((-x1[t]+(Subscript[R, R])x4[t])/((Subscript[C, R])((Subscript[R, 2])+(Subscript[R, R])))),
x2'[t]==0,
x3'[t]==((-x3[t]-(Subscript[R, L])x4[t])/((Subscript[C, L])((Subscript[R, 1])+(Subscript[R, L])))),
x4'[t]==((-x4[t])/L)((Subscript[R, B])+(((Subscript[R, 1])(Subscript[R, L]))/((Subscript[R, 1])+(Subscript[R, L])))+(((Subscript[R, 2])(Subscript[R, R]))/((Subscript[R, 2])+(Subscript[R, R]))))+(((Subscript[R, L])x3[t])/(L((Subscript[R, 1])+(Subscript[R, L]))))-(((Subscript[R, R])x1[t])/(L((Subscript[R, 2])+(Subscript[R, R]))))
};
Equations for Ejection Phase ((Subscript[P, V])>(Subscript[P, A])):  Eh is the active elastance ((Subscript[R, T]) and (Subscript[R, G]) are defined in equations 15 and 16 on page 130)
EjectionPhase={
x1'[t]==((-x1[t]+(Subscript[R, R])x4[t])/((Subscript[C, R])((Subscript[R, 2])+(Subscript[R, R])))),
x2'[t]==(-((Subscript[R, 1])+(Subscript[R, L]))(Subscript[E, V])x2[t]-(Subscript[R, L])x3[t]-(Subscript[R, 1])(Subscript[R, L])x4[t])/(Subscript[R, T]),
x3'[t]==((Subscript[R, L])(Subscript[E, V]) x2[t]-((Subscript[R, L])+(Subscript[R, A]))x3[t]-(Subscript[R, L])(Subscript[R, A])x4[t])/((Subscript[C, L])(Subscript[R, T])),
x4'[t]==((-(Subscript[R, R])x1[t])/(L((Subscript[R, 2])+(Subscript[R, R]))))+((Subscript[R, 1])(Subscript[R, L])(Subscript[E, V])x2[t]+(Subscript[R, L])(Subscript[R, A])x3[t]-((Subscript[R, 1])(Subscript[R, L])(Subscript[R, A])+(Subscript[R, G])(Subscript[R, T]))x4[t])/((Subscript[R, T])L)
};

?EvolveGreene1973Model
\!\(\*"\""\)\!\(\*"EvolveGreene1973Model"\)\!\(\*"["\)\!\(\*"{"\)\!\(\*"x1"\)\!\(\*"["\)\!\(\*"0"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*"x2"\)\!\(\*"["\)\!\(\*"0"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*"x3"\)\!\(\*"["\)\!\(\*"0"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*"x4"\)\!\(\*"["\)\!\(\*"0"\)\!\(\*"]"\)\!\(\*"}"\)\!\(\*","\)\!\(\*" "\)\!\(\*"steps"\)\!\(\*"]"\)\!\(\*" "\)\!\(\*"evolves"\)\!\(\*" "\)\!\(\*"the"\)\!\(\*" "\)\!\(\*"heart"\)\!\(\*" "\)\!\(\*"through"\)\!\(\*" "\)\!\(\*"steps"\)\!\(\*"/"\)\!\(\*"4"\)\!\(\*" "\)\!\(\*"heartbeats"\)\!\(\*","\)\!\(\*" "\)\!\(\*"using"\)\!\(\*" "\)\!\(\*"the"\)\!\(\*" "\)\!\(\*"given"\)\!\(\*" "\)\!\(\*"initial"\)\!\(\*" "\)\!\(\*"conditions"\)\!\(\*"."\)\!\(\*" "\)\!\(\*"x2"\)\!\(\*"["\)\!\(\*"0"\)\!\(\*"]"\)\!\(\*" "\)\!\(\*"is"\)\!\(\*" "\)\!\(\*"effectively"\)\!\(\*" "\)\!\(\*"the"\)\!\(\*" "\)\!\(\*"end"\)\!\(\*"-"\)\!\(\*"diastolic"\)\!\(\*" "\)\!\(\*"volume"\)\!\(\*"."\)\!\(\*"\\n"\)\!\(\*"{"\)\!\(\*"endtime"\)\!\(\*","\)\!\(\*" "\)\!\(\*"r"\)\!\(\*","\)\!\(\*" "\)\!\(\*"alltimes"\)\!\(\*"}"\)\!\(\*" "\)\!\(\*"is"\)\!\(\*" "\)\!\(\*"returned"\)\!\(\*","\)\!\(\*" "\)\!\(\*"where"\)\!\(\*" "\)\!\(\*"endtime"\)\!\(\*" "\)\!\(\*"is"\)\!\(\*" "\)\!\(\*"the"\)\!\(\*" "\)\!\(\*"total"\)\!\(\*" "\)\!\(\*"length"\)\!\(\*" "\)\!\(\*"of"\)\!\(\*" "\)\!\(\*"time"\)\!\(\*" "\)\!\(\*"taken"\)\!\(\*","\)\!\(\*" "\)\!\(\*"r"\)\!\(\*" "\)\!\(\*"are"\)\!\(\*" "\)\!\(\*"rules"\)\!\(\*" "\)\!\(\*"for"\)\!\(\*" "\)\!\(\*"the"\)\!\(\*" "\)\!\(\*"functions"\)\!\(\*" "\)\!\(\*"{"\)\!\(\*"x1"\)\!\(\*"["\)\!\(\*"t"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*"x2"\)\!\(\*"["\)\!\(\*"t"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*"x3"\)\!\(\*"["\)\!\(\*"t"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*"x4"\)\!\(\*"["\)\!\(\*"t"\)\!\(\*"]"\)\!\(\*","\)\!\(\*" "\)\!\(\*SubscriptBox["P", 
RowBox[{"AS", ","}]]\)\!\(\*SubscriptBox["P", "V"]\)\!\(\*","\)\!\(\*" "\)\!\(\*SubscriptBox["P", "A"]\)\!\(\*","\)\!\(\*" "\)\!\(\*SubscriptBox["E", "V"]\)\!\(\*"}"\)\!\(\*","\)\!\(\*" "\)\!\(\*"and"\)\!\(\*" "\)\!\(\*"times"\)\!\(\*" "\)\!\(\*"is"\)\!\(\*" "\)\!\(\*"a"\)\!\(\*" "\)\!\(\*"list"\)\!\(\*" "\)\!\(\*"of"\)\!\(\*" "\)\!\(\*"all"\)\!\(\*" "\)\!\(\*"transitions"\)\!\(\*" "\)\!\(\*"between"\)\!\(\*" "\)\!\(\*"heart"\)\!\(\*" "\)\!\(\*"phases"\)\!\(\*"."\)\!\(\*"\""\)

Options[EvolveGreene1973Model]
{\!\(\*"\""\)\!\(\*"PhaseEquations"\)\!\(\*"\""\):>$PhaseEquations,\!\(\*"\""\)\!\(\*"PhaseConditions"\)\!\(\*"\""\):>$PhaseConditions,\!\(\*"\""\)\!\(\*"HeartParameters"\)\!\(\*"\""\):>$HeartParameters,\!\(\*"\""\)\!\(\*"ArterialParameters"\)\!\(\*"\""\):>$ArterialParameters}
{endtime,rules,times}=EvolveGreene1973Model[{30,0.037,30,0.004},8,\!\(\*"\""\)\!\(\*"PhaseConditions"\)\!\(\*"\""\)->{(Subscript[P, V])<1.05(Subscript[P, A]),(Subscript[P, A])≤1.05(Subscript[P, V]),And[(Subscript[P, AS])<(Subscript[P, V])<=(Subscript[P, A]),t<0.32],t<0.34},
\!\(\*"\""\)\!\(\*"HeartParameters"\)\!\(\*"\""\)->{(Subscript[R, A])->60,(Subscript[R, M])->22,(Subscript[P, AS])->10},
\!\(\*"\""\)\!\(\*"ArterialParameters"\)\!\(\*"\""\)->{(Subscript[C, L])->0.000235`,(Subscript[R, 1])->60,(Subscript[C, R])->0.000133`,(Subscript[R, 2])->300,(Subscript[R, L])->12640,(Subscript[R, B])->60,(Subscript[R, R])->8550,L->21.9`}
];
times
{0,0.020352548724302505`,0.06299089825113122`,0.06828502498062514`,0.34000000294699295`,0.36182283570008844`,0.4024695771535979`,0.40772842961971617`,0.6800000029766966`}
phaselines={Lighter[Red,3/4],Table[Line[{{x,-1},{x,1}}],{x,times}],Table[Line[{{x,-1000},{x,1000}}],{x,times}]};
Table[Plot[var[t]/.rules//Evaluate,{t,0,endtime},PlotLabel->HeartVariableInterpretation[var],ImageSize->400,PlotRange->All,Epilog->phaselines],
{var,{x1,x2,x3,x4}}
]
{,,,}

<<PlotLegends`
Plot[{(Subscript[P, V]),(Subscript[P, A]),(Subscript[P, AS])}/.rules//Evaluate,{t,0,endtime},PlotLegend->{\!\(\*"\""\)\!\(\*"Ventricular"\)\!\(\*" "\)\!\(\*"pressure"\)\!\(\*"\""\),\!\(\*"\""\)\!\(\*"Aortic"\)\!\(\*" "\)\!\(\*"pressure"\)\!\(\*"\""\),\!\(\*"\""\)\!\(\*"Atrial"\)\!\(\*" "\)\!\(\*"source"\)\!\(\*" "\)\!\(\*"pressure"\)\!\(\*"\""\)},LegendPosition->{1,0},ImageSize->Large,PlotRange->{0,All},Epilog->phaselines]

Plot[(Subscript[E, V])/.rules//Evaluate,{t,0,endtime},PlotRange->All,PlotLabel->\!\(\*"\""\)\!\(\*"Elastance"\)\!\(\*" "\)\!\(\*SubscriptBox["E", "V"]\)\!\(\*"\""\),Epilog->phaselines]

The Cardiac cycle

Summary

The cardiac cycle passes through 4 distinct phases determined by which valves are opened or closed and pressure differences between the heart chambers.

We solved the first set of equations until the pressure condition was reached and then moved on to the next phase and continued on to the next phase .

The pressure conditions were:
• Filling phase: Ventricular pressure is greater than atrial pressure. 
• Isolvolumic contraction and relaxation: The pressure in the atrium was less than the pressure in the ventricle  which was less than the pressure in the Aorta.
• Ejection phase: The pressure in the ventricle was greater than pressure in the Aorta.

After solving the equations we found that relying on pressure condition may not be realistic to cause the valves to open and close.

To switch phases we had to use realistic conditions to switch to the next phase ,We added a timing parameter to make the model switch phases we will have to decide on better conditions to switch phases.