Introduction to Modelling in Bioengineering

Abstract

This chapter outlines the basic principles of computational modelling with particular application to physiology and medicine. An overview of the modelling process is provided, along with a description of basic model types, including linear versus non-linear, dynamic versus static, deterministic versus stochastic, continuous versus discrete, and rule-based. Finally, an overview of dimensional analysis and model scaling is also provided. Full Matlab code listings are provided for several example models, including stochastic ion channel kinetics, non-linear passive muscle mechanics, neuronal dendritic branching, and glucose-insulin kinetics. The chapter ends with ten problems covering the fundamentals of modelling, with particular relevance to physiological systems, biology and medicine, with solutions provided at the end of the text.

Problems

1.1

(a) Verify from first principles that the following 1D diffusion equation is linear :

$$\begin{aligned} \frac{\partial c}{\partial t} = D\frac{\partial ^2 c}{\partial x^2} \nonumber \end{aligned}$$

where c is concentration and D is a fixed diffusion coefficient.

(b) Verify using first principles that the following equation governing bacterial population N in a Petri dish is non-linear:

$$\begin{aligned} \frac{\text {d}N}{\text {d}t} = kN\left( 1 - \frac{N}{N_{max}}\right) \nonumber \end{aligned}$$

where k and $N_{max}$ are fixed parameters.

1.2

Determine the dimensions of the following quantities:

• (a) Weight

• (b) Diffusion coefficient (see Problem 1.1a)

• (c) Capacitance

• (d) Resistance

• (e) Cardiac output

• (f) Heart rate

• (g) Radian angle measure

1.3

For the minimal glucose-insulin kinetic model given by

$$\begin{aligned} \frac{\text {d}I}{\text {d}t}= & {} k_2I_p - k_3I \nonumber \\ \frac{\text {d}G}{\text {d}t}= & {} B_0 - (k_1+k_4I+k_5+k_6I)G \end{aligned}$$

(1.25)

the dimensions of I, G and $I_p$ are all $N\,L^{-3}$. Find the corresponding dimensions of $B_0$, $k_1$, $k_2$, $k_3$, $k_4$, $k_5$, and $k_6$.

1.4

(a) In a left ventricular assist pump, the pressure head $\Delta P$ developed across the device is given by the following empirical relationship:

$$\begin{aligned} \Delta P = c_0 + c_1Q_P^3 + c_2\omega ^2 \nonumber \end{aligned}$$

where $\omega$ is the angular velocity of the pump impeller in units of $\mathrm {rad}\,\mathrm {s}^{-1}$ and $Q_P$ is the pump flow rate in $\mathrm {m}^3\,\mathrm {s}^{-1}$. If $\Delta P$ is in units of Pa ($=$N m$^{-2}$), find the dimensions and units of the coefficients $c_0$, $c_1$ and $c_2$.

(b) The opening ($\alpha$) and closing ($\beta$) rate coefficients for potassium ionic channels in a particular excitable cell membrane are given by the following empirical expressions:

$$\begin{aligned} \alpha = \frac{a(V+b)}{1-\mathrm {exp} \left[ \frac{-(V+b)}{c} \right] } \qquad \beta = A\, \mathrm {exp} \left[ \frac{-(V+B)}{C} \right] \nonumber \end{aligned}$$

where V is the transmembrane potential in mV, and $\alpha$, $\beta$ are in units of s$^{-1}$. Find the units of parameters a, b, c, A, B and C.

1.5

(a) The electrical resistance across a block of material of resistivity $\rho$ is given by $R = \rho L/A$, where A is the cross-sectional area and L is the length of the block. Find the dimensions of resistivity.

(b) A monopolar disk stimulating electrode of diameter D is embedded in a medium of resistivity $\rho$. Assuming the medium is an infinite hemisphere centred around the disk, use dimensional analysis to find a formula for the access resistance of the electrode, defined as the resistance between the electrode disk and the hemispherical boundary at infinity.

1.6

The 1D momentum balance equation for the velocity u of a fluid may be written as

$$\begin{aligned} \rho \left( \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} \right) = -\frac{\partial p}{\partial x} + \mu \frac{\partial ^2 u}{\partial x^2} \nonumber \end{aligned}$$

where x is the spatial position, t is time, $\rho$ is fluid density, $\mu$ the viscosity and p is the pressure. Assuming the system has characteristic velocity, length and frequency of V, L, and $\omega$ respectively, scale all variables to dimensionless quantities and obtain the dimensionless form of this equation. How many parameters characterise the dimensionless system?

1.7

The Hodgkin–Huxley [10] equations governing electrical activity in neuronal axons may be written as a coupled system of four differential equations as follows:

$$\begin{aligned} \frac{\text {d}V}{\text {d}t}= & {} -\frac{1}{C} \left[ g_{Na}m^3h(V-V_{Na})+g_Kn^4(V-V_K)+g_{L}(V-V_L) \right] \nonumber \\ \frac{\text {d}n}{\text {d}t}= & {} \alpha _n(1-n) - \beta _n n \nonumber \\ \frac{\text {d}m}{\text {d}t}= & {} \alpha _m(1-m) - \beta _m m \nonumber \\ \frac{\text {d}h}{\text {d}t}= & {} \alpha _h(1-h) - \beta _h h \nonumber \end{aligned}$$

with

$$\begin{aligned} \alpha _n&= \frac{A_n(V+V_{an})}{1-\mathrm {exp} \left[ \frac{-(V+V_{an})}{s_{an}} \right] } \qquad \beta _n&= B_n \mathrm {exp} \left[ \frac{-(V+V_{bn})}{s_{bn}} \right] \nonumber \\ \alpha _m&= \frac{A_m(V+V_{am})}{1-\mathrm {exp} \left[ \frac{-(V+V_{am})}{s_{am}} \right] } \qquad \beta _m&= B_m \mathrm {exp} \left[ \frac{-(V+V_{bm})}{s_{bm}} \right] \nonumber \\ \alpha _h&= A_h \mathrm {exp} \left[ \frac{-(V+V_{ah})}{s_{ah}} \right] \quad \beta _h&= \frac{B_h}{1+\mathrm {exp} \left[ \frac{-(V+V_{bh})}{s_{bh}} \right] } \nonumber \end{aligned}$$

where V is the transmembrane potential (typically in units of mV), n, m, h are dimensionless gating variables , and $i_{\mathrm {stim}}(t)$ is the input stimulus current. The remaining terms represent 25 parameters, namely: C, $g_{Na}$, $V_{Na}$, $g_K$, $V_K$, $g_L$, $V_L$, $A_n$, $V_{an}$, $s_{an}$, $B_n$, $V_{bn}$, $s_{bn}$, $A_m$, $V_{am}$, $s_{am}$, $B_m$, $V_{bm}$, $s_{bm}$, $A_h$, $V_{ah}$, $s_{ah}$, $B_h$, $V_{bh}$ and $s_{bh}$.

Variables V and t may be transformed into corresponding dimensionless variables using

$$\begin{aligned} V^* = \frac{V-V_K}{V_{Na}-V_K} \qquad t^* = B_n t \nonumber \end{aligned}$$

Using these, find the dimensionless form of the above Hodgkin–Huxley equations. How many parameters are needed to characterise this dimensionless system?

1.8

A hydrogel-based optical sensor responds to concentrations of an analyte (c) through changes in the mean wavelength ($\overline{\lambda }$) of its emitted reflectance. The wavelength – log-concentration profile consists of three sequential levels of saturation, as shown below. Formulate a possible model of this behaviour using a system of differential equations.

 

1.9

Steady-state passive electrical behaviour of an unmyelinated nerve cell axon can be approximated by a cylinder of length 10 mm, radius 25 µm, filled with axoplasmic medium of resistivity 0.2 $\Omega$ m, and surrounded by a membrane of resistance 0.1 $\Omega$ m². If a current of 1 mA is injected into one end of the axon, and assuming the extracellular potential is set everywhere to ground, this system can be represented as N discrete circuit elements shown below:

 

where the value of each $R_i$ and $R_m$ will depend on N. The membrane voltage across the $R_m$’s will monotonically decrease from the current source, and the distance at which it has decayed to $\mathrm{e}^{-1} \approx 0.368$ of its maximum value is known as the length constant of the axon. Determine the length constant for $N = 5, 10, 20, 40, 80$ and 160. If the length constant for $N =160$ is assumed to be the exact value, what minimum value of N in this list is sufficient to guarantee a length constant accuracy of 1 %?

1.10

A simplified model of cardiac ventricular electrical activity, based on the cellular automata model of Mitchell et al. [15], consists of a $50\times 50$ square grid, each square representing a small $2\times 2$ mm² region of electrical activity in the heart. The grid is wrapped around into a cylinder, so that the left and right edges are assumed to be in contact. Every region has exactly eight neighbours, with the exception of those squares on the top and bottom edges alone. Each square in the grid can be in one of four states: 0 (quiescent), 1 (relative refractory), 2 (absolute refractory) and 3 (excited). Once a region is excited, it will move through states $3 \rightarrow 2 \rightarrow 1 \rightarrow 0$ in that order, unless it is prematurely excited again.

State	Description	Behaviour
0	quiescent	the region is excited on the next time step if at least one if its eight neighbours is currently excited
3	excited	the region is excited, and can excite neighbouring regions on the next time step
2	absolute refractory	the region cannot be excited, nor can it excite any of its neighbours
1	relative refractory	the region is excited on the next time step if more than one if its eight neighbours is currently excited. The number of excited neighbours required is dependent on the time elapsed since entering this state, according to the following:
		Time elapsed (in ms)	Excited neighbours required
		$\le$2	8
		$\le$4	7
		$\le$6	6
		$\le$8	5
		$\le$12	4
		$\le$20	3
		$\le$50	2

The duration of the excited state is fixed at $T_{ex}$ and the duration of the relative refractory period is fixed at 50 ms. To account for the inhomogeneous refractory properties of the tissue, the duration of the total refractory period (absolute + relative) is a random variable for each region, drawn from a normal distribution of mean $T_{rp}$ and standard deviation $\sigma _{rp}$. The ventricles are excited by a periodic impulse from the atrioventricular (AV) node, of period T, located in a small square on the top border of the grid one quarter of the way in from the left edge. Whenever applied, this impulse can only excite the small square if its current state is not absolute refractory. All state behaviours and parameter values are given on the following page.

Parameter	Description	Value
$T_{ex}$	Duration of excited state	70 ms
$T_{rp}$	Mean total refractory period	250 ms
$\sigma _{rp}$	Standard deviation of total refractory period	100 ms
T	AV node pacing period	0.8 s
$d_T$	Time step increment	2 ms

(a) Implement this model in Matlab, and plot a snapshot of the heart’s electrical activity at $t=1.65$ s.

(b) Decrease the pacing period to 0.2 s, and plot the electrical activity again at $t=1.7$ s. You should observe chaotic activation akin to ventricular fibrillation.