Lumped Parameter Modelling with Ordinary Differential Equations

Abstract

This chapter presents an overview of ordinary differential equations and their use in lumped parameter modelling of physical systems and physiological processes. Methods are described on how to solve some ODEs analytically, followed by a brief overview of numerical solution approaches using Matlab. Detailed examples are presented from models in cardiovascular dynamics and neural activation, followed by sample problems with fully-worked solutions given at the end of the text.

Problems

2.1

In the Hodgkin–Huxley formulation of neural activation, three gating variables n , m, h are employed, satisfying the ODE:

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t} = \alpha _x(V)(1-x) - \beta _x(V)x \nonumber \end{aligned}$$

where $x \equiv n, m, h$ and $\alpha _x(V)$ and $\beta _x(V)$ are known functions of membrane voltage V. Assuming a voltage-clamp experiment is performed, whereby the membrane voltage is stepped suddenly from a value $V_{hold}$ to a new value $V_{clamp}$ and held at this value via a feedback mechanism. $\alpha _x$ and $\beta _x$ are now constant.

(a) Solve this equation analytically for x, with initial value $x(0) = x_0$, stating the homogeneous, particular and general solutions.

(b) What is the steady-state value of x? Hence, what is a reasonable estimate for $x_0$?

2.2

The passive mechanical behaviour of skeletal muscle can be modelled using a simplified lumped parameter representation consisting of a linear spring $k_1$ in series with a parallel linear spring-dashpot combination $k_2$, b, as shown below:

One end of the muscle is fixed, and the displacement of the other end is x. If $x_1$ denotes the change in length from rest in spring $k_1$, then the forces in each element are given by:

$$\begin{aligned} F_1 = k_1x_1 \qquad F_2 = k_2x_2 \qquad F_b = b\frac{\mathrm{d}x_2}{\mathrm{d}t} \nonumber \end{aligned}$$

where $F_1$, $F_2$ and $F_b$ refer to elements $k_1$, $k_2$ and b respectively, and $x_2~=~x-x_1$.

(a) If the length x of the muscle is suddenly stepped and held from 0 to $X_m$ at $t=0$, solve the system analytically for the applied force, F.

(b) If the applied force on the muscle F is suddenly stepped and held from 0 to $F_m$ at $t=0$, find the analytical solution for the change in length, x.

2.3

A simple model of neuronal excitation represents the cell membrane as a resistance R in parallel with a capacitance C. An applied stimulus current I depolarises the membrane to a potential of V, as shown below. If V exceeds a pre-defined threshold $V_{th}$, the neuron will fire.

The currents $i_R$ and $i_C$ flowing through the R and C branches are given by

$$\begin{aligned} i_R = \frac{V}{R} \qquad i_C = C\frac{\mathrm{d}V}{\mathrm{d}t} \nonumber \end{aligned}$$

(a) Assuming the neuron is initially at rest with $V=0$, and a constant stimulus current I is applied at $t=0$, find the time taken T to depolarize the membrane to $V_{th}$. Determine the corresponding stimulus strength-duration characteristic for neuronal activation, i.e. I as a function of T.

(b) Defining the rheobase as the minimum current necessary to activate the neuron, and the chronaxie as the required stimulus duration for an applied current of twice the rheobase, determine both quantities from the above strength-duration characteristic.

2.4

Consider the system below of two coupled masses M, connected to each other and to fixed supports via three linear springs with spring constants k:

(a) Determine the pair of ODEs describing the motion of this system.

(b) Solve this system analytically for the displacements $x_1$ and $x_2$, assuming the masses are initially at rest and displaced by amounts $u_1$ and $u_2$.

Hint: Use the variable substitutions $y_1 = x_1+x_2$, $y_2 = x_1-x_2$.

2.5

A simplified two-compartment model of glucose-insulin kinetics in a human subject proposed by Berman et al. [1]⁴ is shown below. $I_p(t)$ represents insulin injected intravenously into the blood, I is the concentration of insulin in a remote body compartment, and G is the glucose concentration in the blood plasma.

$I_p$, I and G are all in units of mM, with model parameters given below:

 

Parameter	Value	Parameter	Value
$k_1$	0.015 min$^{-1}$	$k_5$	0.035 min$^{-1}$
$k_2$	1 min$^{-1}$	$k_6$	0.02 mM$^{-1}$ min$^{-1}$
$k_3$	0.09 min$^{-1}$	$B_0$	0.5 mM min$^{-1}$
$k_4$	0.01 mM$^{-1}$ min$^{-1}$

An intravenous dose of insulin is administered as a square-pulse waveform according to:

$$\begin{aligned} I_p(t) = \left\{ \begin{array}{ll} 200~\mathrm {mM} \quad &{} 0 \le t < 0.1~\mathrm {min} \\ 0 &{} t \ge 0.1~\mathrm {min} \end{array} \right. \nonumber \end{aligned}$$

The initial values of I and G at $t=0$ are 0 and 10 mM respectively.

Solve for I and G using Matlab over the time interval $0 \le t \le 60~\mathrm {min}$.

2.6

A simplified lumped parameter electric-analogue model of the heart and systemic circulation⁵ is shown below:

where $P_s$ is the systemic pressure and $L_o$ represents the blood inertance within the aortic root such that the pressure drop across this element is given $L_o\frac{\mathrm{d}Q_L}{\mathrm{d}t}$, where $Q_L$ is the flow (in cm³ s$^{-1}$) through it. $P_v(t)$ represents the developed ventricular pressure as a function of time, given by the following simplified periodic square-wave profile:

$$\begin{aligned} P_v(t+T)= & {} P_v(t) \nonumber \\ P_v(t)= & {} \left\{ \begin{array}{ll} P \quad &{} 0 \le t < t_c \\ 0 &{} t_c \le t \le T \end{array} \right. \nonumber \end{aligned}$$

All other elements are similar to those given earlier in the example of Sect. 2.3.1. Remaining model parameters and descriptions are given below:

 

Parameter	Description	Value
$R_o$	Aortic root resistance	0.06 mmHg s cm$^{-3}$
$L_o$	Aortic root blood inertance	0.2 mmHg s2 cm$^{-3}$
$C_s$	Systemic compliance	1 cm$^{3}$ mmHg$^{-1}$
$R_s$	Systemic resistance	1.4 mmHg s $^{-3}$
P	Peak ventricular pressure	120 mmHg
T	Heart period	1 s
$t_c$	Active contraction interval (systole)	0.35 s

(a) Determine the ODEs describing this system.

(b) Using an appropriate choice of initial values, solve this system using Matlab for the steady-state oscillations in systemic pressure $P_s$.

2.7

The following set of ODEs modified from McSharry et al. [8] reproduce a synthetic electrocardiogram (ECG) waveform in variable z:

$$\begin{aligned} \frac{\mathrm{d}x}{\mathrm{d}t}= & {} \alpha x - \omega y \nonumber \\ \frac{\mathrm{d}y}{\mathrm{d}t}= & {} \alpha y + \omega x \nonumber \\ \frac{\mathrm{d}z}{\mathrm{d}t}= & {} -\sum _{i=1}^5 a_i(\theta -\theta _i)\exp \left( -\frac{(\theta -\theta _i)^2}{2b_i^2}\right) - (z-z_0) \nonumber \end{aligned}$$

where $\alpha = 1 -\sqrt{x^2 + y^2}$, $\omega = 2\pi$ rad s$^{-1}$, $z_0 = 0$ and $\theta = \mathrm {atan2}(y,x)$, where atan2 represents the four-quadrant inverse tangent implemented by the Matlab function atan2. Remaining model parameters are given below:

 

Index i	1	2	3	4	5
$\theta _i$ (rad)	$-\frac{1}{3}\pi$	$-\frac{1}{12}\pi$	0	$\frac{1}{12}\pi$	$\frac{1}{2}\pi$
$a_i$ (mV s$^{-1}$ rad $^{-1}$)	1.2	−5	30	−7.5	0.75
$b_i$ (rad)	0.25	0.1	0.1	0.1	0.4

Given the initial values, $x(0) = -1$, $y(0) = 0$, $z(0) = 0$, numerically solve this ECG model in Matlab from $t=0$ to 1 s, plotting z against t, where z, t are in units of mV and s respectively.

2.8

The Frankenhaeuser-Huxley neural action potential model [3] consists of the following ODE system:

$$\begin{aligned} \frac{\mathrm{d}V}{\mathrm{d}t}= & {} -\frac{1}{C_m}\left[ i_{Na}+i_K+i_P+i_L-i_{stim}\right] \nonumber \\ \frac{\mathrm{d}m}{\mathrm{d}t}= & {} \alpha _m\left( 1-m\right) - \beta _mm \nonumber \\ \frac{\mathrm{d}h}{\mathrm{d}t}= & {} \alpha _h\left( 1-h\right) - \beta _hh \nonumber \\ \frac{\mathrm{d}n}{\mathrm{d}t}= & {} \alpha _n\left( 1-n\right) - \beta _nn \nonumber \\ \frac{\mathrm{d}p}{\mathrm{d}t}= & {} \alpha _p\left( 1-p\right) - \beta _pp \nonumber \end{aligned}$$

with membrane ionic currents given by

$$\begin{aligned} i_{Na}= & {} m^2h\bar{P}_{Na}\left( \frac{EF^2}{RT}\right) \left[ \frac{[\mathrm {Na}]_{\mathrm {o}}-[\mathrm {Na}]_{\mathrm {i}}\exp \left( \frac{EF}{RT}\right) }{1-\exp \left( \frac{EF}{RT}\right) } \right] \nonumber \\ i_K= & {} n^2 \bar{P}_K\left( \frac{EF^2}{RT}\right) \left[ \frac{[\mathrm {K}]_{\mathrm {o}}-[\mathrm {K}]_{\mathrm {i}}\exp \left( \frac{EF}{RT}\right) }{1-\exp \left( \frac{EF}{RT}\right) } \right] \nonumber \\ i_P= & {} p^2\bar{P}_{P}\left( \frac{EF^2}{RT}\right) \left[ \frac{[\mathrm {Na}]_{\mathrm {o}}-[\mathrm {Na}]_{\mathrm {i}}\exp \left( \frac{EF}{RT}\right) }{1-\exp \left( \frac{EF}{RT}\right) } \right] \nonumber \\ i_L= & {} g_L(V-V_L) \nonumber \end{aligned}$$

where E is the transmembrane potential, V is the membrane potential displacement from its resting value $E_r$ ($V = E-E_r$), F is Faraday’s constant, R is the gas constant, and T is the absolute temperature. [Na]$_\text {o}$, [Na]$_\text {i}$, [K]$_\text {o}$ and [K]$_\text {i}$ represent outside (extracellular) and intracellular Na$^{+}$ and K$^{+}$ concentrations, whilst $i_{Na}$, $i_K$, $i_P$ and $i_L$ represent the membrane Na$^{+}$, K$^{+}$, non-specific (mainly Na$^{+}$), and leakage currents respectively. The membrane permeabilities of $i_{Na}$, $i_K$ and $i_P$ are $\bar{P}_{Na}$, $\bar{P}_K$ and $\bar{P}_P$ respectively, with the kinetics of these currents determined from the m, h, n and p gating variables. $i_{stim}$ is an applied intracellular stimulus current given by

$$\begin{aligned} i_{stim} = \left\{ \begin{array}{ll} I_s &{} t_{on} \le t < t_{on} + t_{dur} \\ 0 &{} \mathrm {otherwise} \end{array} \right. \ \nonumber \end{aligned}$$

where $I_s$, $t_{on}$ and $t_{dur}$ represent the stimulus current amplitude, onset time and duration respectively. The voltage-dependent forward ($\alpha$) and reverse ($\beta$) rates (in ms$^{-1}$) are given by:

$$\begin{aligned} \alpha _m&= \frac{0.36(V-22)}{1-\exp \left[ \frac{(22-V)}{3}\right] } \qquad \beta _m&= \frac{0.4(13-V)}{1-\exp \left[ \frac{(V-13)}{20}\right] } \nonumber \\ \alpha _h&= \frac{0.1(-10-V)}{1-\exp \left[ \frac{(V+10)}{6}\right] } \qquad \beta _h&= \frac{4.5}{1+\exp \left[ \frac{(45-V)}{10}\right] } \nonumber \\ \alpha _n&= \frac{0.02(V-35)}{1-\exp \left[ \frac{(35-V)}{10}\right] } \qquad \beta _n&= \frac{0.05(10-V)}{1-\exp \left[ \frac{(V-10)}{10}\right] } \nonumber \\ \alpha _p&= \frac{0.006(V-40)}{1-\exp \left[ \frac{(40-V)}{10}\right] } \qquad \beta _p&= \frac{0.09(-25-V)}{1-\exp \left[ \frac{(V+25)}{20}\right] } \nonumber \end{aligned}$$

where V is in units of mV. Initial values are $V=0\,\mathrm {mV}$, $m=0.0005$, $h = 0.8249$, $n=0.0268$ and $p=0.0049$.

The drug tetrodotoxin (TTX) is known to selectively block the membrane $i_{Na}$ current. Assuming that at one given dosage, TTX reduces parameter $\bar{P}_{Na}$ to 20% of its original value. Solve this model using Matlab over the time interval $t=0$ to 5 ms, plotting the transmembrane potentials E (in mV) on the same graph for the following two cases:

(1) control (i.e. no TTX) and

(2) TTX applied.

All model parameters are given in the following table:

 

Parameter	Value	Parameter	Value
$C_m$	2 µF cm$^{-2}$	$[\mathrm {K}]_\mathrm {o}$	2.5 mM
$\bar{P}_{Na}$	0.008 cm s$^{-1}$	$[\mathrm {K}]_\mathrm {i}$	120 mM
$\bar{P}_K$	0.0012 cm s$^{-1}$	$I_s$	1 mA cm$^{-2}$
$\bar{P}_P$	0.00054 cm s$^{-1}$	$t_{on}$	1 ms
$g_L$	30.3 ms cm$^{-2}$	$t_{dur}$	0.12 ms
$V_L$	0.026 mV	F	96.49 C mmol$^{-1}$
$[\mathrm {Na}]_\mathrm {o}$	114.5 mM	R	8.31 J mol$^{-1}$ K$^{-1}$
$[\mathrm {Na}]_\mathrm {i}$	13.74 mM	T	310 K
$E_r$	−70 mV

2.9

A simple three-element model of active cardiac muscle contraction consists of a passive non-linear spring in parallel with a contractile and passive series element, as shown in the diagram. The model structure and equations have been modified from Fung [4].

The total tension T in the muscle is given by the sum of tensions in the parallel and series elements:

$$\begin{aligned} T = T_p+T_s \nonumber \end{aligned}$$

where $T_p$ and $T_s$, the tensions in the parallel and series elements respectively, are given by:

$$\begin{aligned} T_p = \beta \left( \mathrm{e}^{\alpha (L-L_0)}-1\right) , \qquad T_s = \beta \left( \mathrm{e}^{\alpha L_s}-1\right) \nonumber \end{aligned}$$

where $\alpha$, $\beta$ are parameters and $L_0$ denotes the resting length of the muscle. For the contractile element, the velocity of its shortening is described by:

$$\begin{aligned} \frac{\mathrm{d}L_c}{\mathrm{d}t} = \frac{a\left[ T_s-S_0f(t)\right] }{T_s + \gamma S_0} \nonumber \end{aligned}$$

where a, $\gamma$, $S_0$ are muscle parameters and f(t) is the muscle activation function given by

$$\begin{aligned} f(t) = \left\{ \begin{array}{ll} \sin \left( \frac{\pi }{2} \left[ \frac{t+t_0}{t_{ip}+t_0} \right] \right) \qquad &{} 0 \le t < 2t_{ip} + t_0 \\ 0 &{} t \ge 2t_{ip} + t_0 \end{array} \right. \ \nonumber \end{aligned}$$

with $t_0$ and $t_{ip}$ constants defining activation phase offset and the time to peak isometric contraction respectively.

Model parameters for a cardiac papillary muscle specimen are given in the table below:

 

Parameter	Value	Parameter	Value
$L_0$	10 mm	$S_0$	4 mN
$\alpha$	15 mm	$\gamma$	0.45
$\beta$	5 mN	$t_0$	0.05 s
a	0.66 mm $^{-1}$	$t_{ip}$	0.2 s

Note that in the relaxed state, the length of the series element $L_s$ is 0.

(a) Using Matlab, solve for and plot the total tension T against time during an isometric contraction in which the muscle is clamped at its resting length $L_0$.

(b) Solve for and plot muscle length L against time during an isotonic contraction in which the muscle is allowed to freely contract with no imposed load (i.e. $T = 0$).