## Some Visualizations in Mathematical Physiology

“…the equations we deal with are probably more complicated than even most physical scientists are accustomed to. This is because the phenomena we are attempting to describe are generally more complex than most physical systems, although it may reflect our own ineptness in perceiving their underlying simplicity.” -James Murray, Mathematical Biology

Due to the complexity of modelling biological phenomena, analytical and numerical approaches are often used together to give a more thorough understanding of a model than either could provide alone. This is particularly true in the area of mathematical physiology, where there has been a dual development of increasingly sophisticated models, as well as analytical and numerical tools to analyze these models. In this post I want to describe some visualizations of simple models in this area, as well as reference some useful tools to explore these simulations further. These visualizations complement the lecture notes for a corresponding course that I tutored last year, which can be found here for the moment, and a much older version of the notes can be found here.

Mathematically these models are small systems of ordinary differential equations, delay differential equations, and spatial variants consisting of coupled partial differential equations. I will not exhaustively discuss or describe these models but will instead link to relevant literature, and primarily discuss particular example simulations demonstrating relevant behaviours that complement the lecture notes. All of the simulations are performed using different variants of the explicit fourth-order Runge-Kutta method, which can be derived via Taylor series.

1. Hodgkin-Huxley Equations for Transmembrane Ion Transport

The Hodgkin-Huxley Equations were first proposed in a series of papers in the early 1950s concerned with the flow of electric current, mediated by the transport of ions, through the surface membrane of a giant axon from a squid (not to be confused with an axon from a giant squid). They have formed the basis for a wide range of mathematical and scientific work involving ion channels in cellular physiology in particular, and conductance-based excitable media in general. The equations, using the same notation as the lecture notes cited above, can be written as

$\displaystyle - C \frac{d V}{d t} = g_{Na} m^3 h (V-V_{Na})+g_K n^4 (V-V_K)+g_L (V-V_L)-I_{ext}, \ \ \ \ \ (1)$

$\displaystyle \tau_m (V)\frac{d m}{dt} = m_{\infty} (V)-m, \ \ \ \ \ (2)$

$\displaystyle \tau_h(V)\frac{d h}{dt} = h_{\infty} (V)-h, \ \ \ \ \ (3)$

$\displaystyle \tau_n(V)\frac{d n}{dt} = n_{\infty} (V)-n, \ \ \ \ \ (4)$

where ${C}$ is the capacitance of the membrane, ${V}$ is the voltage (potential difference), ${n}$ represents the proportion of open ion channels in the membrane for potassium (K), ${h}$ and ${m}$ are proportions for two different sodium (Na) channels, ${g_{Na}}$, ${g_K}$, ${V_{Na}}$, ${V_K}$ are ion-specific conductances and Nernst-potentials respectively, ${g_L}$ and ${V_L}$ model a leakage current, ${I_{ext}}$ is an externally-imposed current, and finally ${\tau_m (V)}$, ${\tau_h (V)}$, ${\tau_n (V)}$, ${m_{\infty}}$, ${h_{\infty}}$, and ${n_{\infty}}$ represent the timescales and equilibrium values of the gating variables ${m}$, ${h}$ and ${n}$ respectively. Note that these last six parameters are functions of the membrane potential ${V}$, which models the voltage-dependence of the ion gates themselves, although typically these functions are determined experimentally by fitting data to physiologically sensible functional forms.

These equations have many interesting features important in quantitative and qualitative modelling in biophysics. Some of the original motivation came from modelling spiking behaviour in axons, where a sufficiently large external current prompted a large increase in the membrane potential known as an action potential. This corresponds in neurophysiology to individual neurons firing, and has analogues in other biological systems such as the cardiac action potential that gives rise to the periodic contractile behaviour of the heart.

You can observe this behaviour yourself using a JavaScript tool I wrote that solves Equations (1)(4) using Runge-Kutta. The timescales and equilibrium values of the gating variables are given as rational functions involving the membrane potential ${V}$, as done in the paper cited at the beginning of this section. As an interesting historical aside, note that the figures produced in that paper were laboriously computed by hand using a tedious numerical procedure described on page 523. Thankfully, your web browser will allow you to reproduce their results in a fraction of a second, rather than several weeks of tedious arithmetic.

As an external current is added (at ${t=25}$), a small deviation from the resting potential is observed. Increasing the value of ${I_{ext}}$ above ${3}$ causes a sharp rise in the membrane potential. At ${I_{ext}=6}$ two spikes are observed. Finally, for ${I_{ext}\geq 7}$, a continuous spiking behaviour is observed. This can be understood via Hopf bifurcations as the external current is varied, although the calculations for the full four-dimensional system become quite cumbersome, and so the use of sophisticated mathematical and computational tools is typically used to analyze these bifurcations. Another common approach is to use typical values of the timescales and equilibrium values of the gating variables to asymptotically reduce Equations (1)(4) to a system of two equations, where classical phase-plane analysis (as well as our intuition for two-dimensional systems) can be used. This is done in the lecture notes mentioned at the beginning, as well as in the book by Keener and Sneyd.

Some intuition for both the action potential itself, and sustained oscillations, can be gained from analyzing a caricature of this reduced two-dimensional model. However, there are qualitative behaviours in the full system that cannot be captured by a two-dimensional reduction or a caricature thereof. These include bifurcations not found in planar dynamical systems, such as the twisting of homoclinic orbits reported in this paper. Additionally, for small subsets of the parameter space, numerical evidence for intermittent deterministic chaos can be found, which cannot exist in any two-dimensional vector field (unless there is a non-autonomous forcing, such as a time-dependent external current).

You can use the web application linked above to explore some of these dynamics further. One way to do this is to introduce a noise term into the external current. the variable ${I_{rand}}$ does this by adding white noise to this current. Without discussing the technicalities involved in defining this noise term, you can think of it as for every time step ${\delta t}$, the external current is modified by a normally distributed random variable with zero mean and variance ${\delta t}$. Adding a deterministic constant shifts the mean of this process by the same amount, and multiplying by a constant multiplies the variance at each step. If you set the deterministic part of the external current ${I_{ext}=6}$, and the random part ${I_{rand}=3}$ you can see that while the total external current rarely moves far from the mean, it can induce action potentials in random intervals of the time series. I encourage you to explore the effects this noise has on inducing (or collapsing) sustained oscillations for varying these two values. If the variable ${RandSeed}$ is set, then the simulation uses the same `realization’ of the noise each time the simulation is computed. If it is not set then a different realization is used every time the simulation is done.

Another interesting effect can be seen by setting the resting potential ${V_{rest}=32}$. Increasing the resting potential from the physiological value by a factor of about three causes the membrane potential to spontaneously oscillate on its own. If you now introduce an external current, by slowly increasing ${I_{ext}}$ using the sliders, you will see that for some values of the current, the oscillations appear to collapse after the current is turned off. To see this more clearly, you can change the range of the time series by setting ${tStop=500}$. Now, for the case when ${I_{ext}=4}$, if you increase the initial condition on the membrane potential (which corresponds to shifting the phase of the oscillation while the external current is present) for ${0 \leq V(0) \lesssim 12}$, the membrane potential stops oscillating when the current is turned off after ${t=175}$. However for ${V(0) = 13}$, the oscillation becomes self-sustained. This is due to the limit cycle and the steady state both being locally asymptotically stable, and changing the initial condition changes which basin of attraction the solution falls into as time goes on. This is an example of multistability, as the long-time dynamics of the system are dependent upon the initial conditions. Quasi-periodic behaviours can also be seen by increasing ${V_{rest}}$ and ${I_{ext}}$.

While the above asymptotic behaviours can occur in systems of ODEs with fewer equations, it is important to remember that planar (two-dimensional) dynamical systems can only exhibit steady state behaviour or oscillations for asymptotically long periods of time. Adding non-autonomous effects (such as the explicitly time-dependent noise term ${I_{rand}}$) does increase the effective dimension of the dynamics as well, but unless this non-autonomous term induces a kind of resonant behaviour as ${t \rightarrow \infty}$, the long-time dynamics will still be either a limit cycle or a steady state.

2. FitzHugh-Nagumo Equations and Phase Plane Analysis

Next, I want to discuss a few tools to explore planar dynamical systems numerically, in order to help develop some intuition for these simplified systems. Here is a generic solver for two-dimensional systems of ODEs that I wrote to illustrate certain bifurcation behaviours for the FitzHugh-Nagumo equations, which are a caricature of Equations (1)(4). The FHN equations can be written as,

$\displaystyle \epsilon \frac{d v}{d \tau} = Av(v-a)(1-v)-w+I_{ext}^*, \ \ \ \ \ (5)$

$\displaystyle \frac{d w}{d \tau} = -w + bv, \ \ \ \ \ (6)$

where ${0 < \epsilon \ll 1}$, ${A>0}$, ${a \in (0,1)}$, and ${b>0}$ is sufficiently large to guarantee that there is only one equilibrium point. Using the application above, you can see the emergence of oscillations by increasing the rightmost term in the first equation (corresponding to ${I_{ext}}$) from ${0}$ to ${0.3}$ through ${0.4}$.

There is also this excellent JavaScript solver written by Daryl Nester that shows a direction field and plots trajectories of the system starting from a given initial point. Eventually I hope to extend this application to include plots of the nullclines, but you can find these easily for the FHN system. For now, you can get an idea of the role that these nullclines and their intersections (corresponding to steady states) have on the dynamics by varying the parameters and visualizing various orbits in the plane. Performing the same increase in ${I_{ext}}$, you can see the emergence of a small amplitude limit cycle for ${I_{ext}=3.2}$, which becomes a large-amplitude limit cycle for larger external currents.

If you look under, “Numerical solution, timeplots, and tables,” you will see a time series plot for the FHN system. Hovering over various sections of the time series will display where that point corresponds to in the phase plane diagram above, and so this can be used to see the “fast” dynamics along the horizontal part of the orbit, and the “slow” dynamics in the vertical parts of the trajectory, which correspond to the limit cycle crawling along a nullcline. These timescale arguments presented in the lecture notes (and elsewhere) can be generalized beyond the two-dimensional setting and applied to the original four-dimensional system above, but the details are more complicated, and it is not generically true that oscillations in these systems always separate into fast and slow components that have a simple structure as in the two-equation case. See this paper for some of the details in this direction.

Alongside this system are several other two-dimensional models, such as the Two-pool model of cellular calcium dynamics and the Lotka-Volterra Predator-Prey equations that can be explored using phase-plane analysis. Murray’s first volume of Mathematical Biology goes through this analysis in detail for many of these biological systems, and analysis of the Two-Pool model can be found in the book by Keener and Sneyd mentioned above.

3. Delay and Partial Differential Equations

Finally I want to discuss two other common kinds of equations used to model physiological phenomena, as well as other things in mathematical biology. The first of these are partial differential equations, often in the form of reaction-diffusion equations. The second kind of model that has seen recent interest in the past few decades are delay differential equations.

The former can sometimes be derived using tools from statistical mechanics or by analogy with well-understood physicochemical processes, such as diffusion in a stationary fluid. Often qualitative insight can be gained by phenomenologically adding a diffusion term to systems of ordinary differential equations, such as equations (1)(4) or (5)(6). Analytical and numerical analysis is typically more difficult to do for these systems, partly due to the infinite-dimensionality of PDEs, and the complicating effects that boundaries have on the dynamics. Nevertheless, several interesting phenomena can be found in these models, such as spiral waves, and their higher-dimensional analogues, scroll waves. Figure 1 of the article on scroll waves gives an interesting visualization of how these higher-dimensional waves can propagate and interact in non-intuitive ways.

Lastly I wanted to briefly mention delayed dynamical systems. These kinds of equations are less well-studied than ODEs and PDEs, but there is a growing literature analyzing them from both a theoretical point of view, and using them to model a variety of phenomena. One of these is the Mackey-Glass equation, which was originally proposed to demonstrate bifurcations that may be relevant to physiological systems, such as the change from a periodic heart-rate to an aperiodic arrythmia, or the transition from periodic breathing to Cheyne Stokes respiration. Note that there are several non-equivalent Mackey-Glass equations in the literature. As with PDEs, these equations are not as easy to simulate as ordinary differential equations. Partly this is due to the complicated phase space (which is also infinite-dimensional), and partly because these equations have not had as much attention as the other kinds of models historically. There are good numerical routines in Matlab, Octave, and Maple, however, so I would encourage you to explore some of these models. Here is a nice illustration of some of these delay models in Maple that I found useful.