BME 355: Physiological Systems Modelling
Estimated study time: 10 minutes
Table of contents
Sources and References
Primary texts — Keener and Sneyd, Mathematical Physiology (2 vols.), 2nd ed. (Springer). Khoo, Physiological Control Systems: Analysis, Simulation, and Estimation, 2nd ed. (Wiley-IEEE).
Supplementary texts — Hoppensteadt and Peskin, Modeling and Simulation in Medicine and the Life Sciences, 2nd ed. (Springer). Ljung, System Identification: Theory for the User, 2nd ed. (Prentice Hall). Enderle, Introduction to Biomedical Engineering, 3rd ed. (Academic Press).
Online resources — MIT OCW HST.508 Biological Engineering II and HST.582J Biomedical Signal and Image Processing. OpenSim open-source musculoskeletal modelling platform (Stanford). Physiome Project models (open CellML repository). Hodgkin–Huxley ModelDB entries. MATLAB/Simulink physiological modelling examples.
Chapter 1: The Modelling Mindset
1.1 What a Model Is
A physiological model is a mathematical statement of how a quantity of interest evolves in time or space under identified inputs and parameters. A good model is fit for purpose: complex enough to capture the phenomenon at the scale of interest, simple enough that its behaviour can be interrogated analytically or computationally.
1.2 Systems Theory
Regardless of the biological substrate, models share a common vocabulary: state \( \mathbf{x} \), input \( \mathbf{u} \), output \( \mathbf{y} \), dynamics \( \dot{\mathbf{x}} = f(\mathbf{x}, \mathbf{u}, t; \boldsymbol{\theta}) \), output map \( \mathbf{y} = h(\mathbf{x}, \mathbf{u}; \boldsymbol{\theta}) \). Linear time-invariant models with \( f = A\mathbf{x} + B\mathbf{u} \) admit transfer-function and state-space analysis; nonlinear models require numerical integration and care about stability.
Chapter 2: Compartmental Models
2.1 Single Compartment
A single well-mixed compartment with volume \( V \), concentration \( C \), inflow \( q_{in} C_{in} \), outflow \( q_{out} C \), and elimination \( kVC \) gives
\[ V\frac{dC}{dt} = q_{in}C_{in} - q_{out}C - kVC . \]For pharmacokinetics this becomes a first-order decay \( C(t) = C_0 e^{-kt} \) with half-life \( t_{1/2} = \ln 2/k \). Clearance \( Cl = kV \) and volume of distribution \( V_d \) are the classical descriptive parameters.
2.2 Multi-Compartment Systems
Transfer rates between compartments assemble into a matrix \( A \) with column sums set by conservation and diagonal entries equal to negative sums. Eigenvalues determine characteristic time scales; eigenvectors reveal modes of exchange. Two-compartment models capture the fast distribution phase and the slow elimination phase observed after drug bolus injection.
Chapter 3: Systems of the Body
3.1 Skeletal and Muscular Systems
Hill-type muscle models combine an active contractile element with series and parallel elastic elements:
\[ F_{\text{muscle}} = F_0\, f_L(\ell)\, f_V(\dot{\ell})\, a(t) + F_{\text{PE}}(\ell) , \]with length-tension \( f_L \), force-velocity \( f_V \), activation dynamics \( a(t) \), and passive elasticity \( F_{\text{PE}} \). OpenSim and similar platforms assemble musculoskeletal multibody models where inverse dynamics computes joint torques from kinematics, and static optimization distributes torques among redundant muscles.
3.2 Neuromuscular Control
Motor unit recruitment follows Henneman’s size principle: small, slow, fatigue-resistant units recruited first. The relationship between EMG envelope and muscle force is nonlinear and subject-specific. Activation dynamics are modelled as low-pass first-order systems
\[ \dot{a} = \frac{u - a}{\tau(u,a)} , \]with time constant depending on whether the muscle is activating or deactivating.
3.3 Central Nervous System
Neuron models span scales: single-compartment integrate-and-fire, Hodgkin–Huxley ionic, multi-compartment cable. The cable equation
\[ \lambda^2 \frac{\partial^2 V}{\partial x^2} = \tau \frac{\partial V}{\partial t} + V - V_{\text{rest}} \]describes passive membrane with space constant \( \lambda = \sqrt{r_m/r_i} \) and time constant \( \tau = r_m c_m \). Network models (integrate-and-fire, rate-based) describe population dynamics underlying reflexes and movement planning.
3.4 Cardiovascular System
Windkessel models represent the arterial system by a capacitor (compliance) and resistor:
\[ C \frac{dP}{dt} = Q_{\text{in}}(t) - \frac{P}{R} . \]Three- and four-element extensions add characteristic impedance and inertance. Lumped-parameter heart models couple atrial and ventricular chambers through time-varying elastance \( E(t) \); pressure–volume loops quantify stroke work.
3.5 Respiratory System
A linear respiratory mechanics model
\[ P_{\text{airway}} = R\dot{V} + E V + P_0 \]identifies resistance \( R \) and elastance \( E \) from ventilator data. Gas-exchange models combine alveolar ventilation with capillary blood flow and \( \mathrm{O_2} \) and \( \mathrm{CO_2} \) dissociation curves.
Chapter 4: Parameter Identification
4.1 Least Squares
Given model \( y = h(\mathbf{x}(t, \boldsymbol{\theta})) \) and measurements \( \{y_k\} \), parameter estimation minimizes
\[ J(\boldsymbol{\theta}) = \sum_k w_k\,(y_k - h(\mathbf{x}(t_k, \boldsymbol{\theta})))^2 . \]Gauss–Newton, Levenberg–Marquardt, and interior-point algorithms solve the resulting nonlinear optimization. Initial guesses and parameter bounds are critical in physiological problems, where likelihood landscapes have many local minima.
4.2 Identifiability
Structural identifiability asks whether parameters are uniquely determined by the model structure and input–output map given noise-free data. Practical identifiability addresses noisy, finite-sample data. Profile likelihood and Fisher-information analysis identify parameters that are informed, sloppy, or unidentifiable. Reparametrizing to physically meaningful combinations often restores identifiability.
4.3 Experimental Design
The Fisher information matrix \( F = \sum_k \nabla h_k \nabla h_k^T / \sigma_k^2 \) measures the information content of a protocol. D-optimal designs maximize \( \det F \); input signals and sampling times are chosen to maximize identifiability given experimental constraints.
Chapter 5: Simulation and Sensitivity
5.1 Numerical Integration
Explicit methods (RK4) suit non-stiff problems; implicit methods (BDF, trapezoidal) handle stiff systems where time scales span many orders of magnitude, typical in ionic models. Adaptive step size is essential; event detection (e.g., action-potential upstroke crossing) requires dense-output interpolation.
5.2 Sensitivity Analysis
Local sensitivity \( S_{ij} = \partial y_i / \partial \theta_j \) is computed by finite difference or forward sensitivity equations. Global sensitivity (Sobol indices, Morris screening) captures nonlinear, interaction effects across full parameter ranges. Sensitivity results guide which parameters to identify, which to fix, and where additional experiments would pay off.
5.3 Model Validation
Validation compares model predictions to independent data not used for identification. Metrics include normalized RMSE, concordance correlation, and residual diagnostics. Falsification — finding where the model fails — is as valuable as successful fitting. A validated model is always provisionally validated: valid within the conditions tested.
Chapter 6: Applications
6.1 Diagnostic Model-Based Metrics
Model-derived quantities (ejection fraction estimated from image segmentation, insulin sensitivity from glucose tracer experiments) extend direct measurement. Their reliability depends on model fidelity and identifiability in the clinical protocol.
6.2 Model-Based Control
Artificial pancreas, closed-loop ventilators, and neurostimulator adaptation use physiological models as internal plants for model-predictive control. Safety constraints dominate design: robust performance across patient variability outweighs nominal optimality.
6.3 Digital Twins and Population Models
Subject-specific models calibrated from wearable and imaging data form patient “digital twins” supporting personalized therapy planning. Population models pool information across subjects via nonlinear mixed effects, estimating population-typical and individual parameters simultaneously. These tools blur the line between physiology, statistics, and engineering — which is exactly where the biomedical modeller should stand.
