This section provides some applications of differentrial equations to modeling neuroscience problems.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part III of the course APMA0340
Introduction to Linear Algebra with Mathematica


The building block of neurological systems is the neuron that happens to be a biological cell exhibiting physical, electrical, thermal, chemical, mechanical, and physiological properties. The neuron as a living cell performs a multitude of functions, all fundamentally dependent upon the capability of its permeable membrane to control in- and out-flows of ions. The resulting variable voltage across its membrane is the neuron’s most valuable asset, with the action potential. While performing experiments on the living neuron provides a wealth of information on how neurons work, the inclusion of physical and mathematical descriptions of the cell’s behavior may help elucidate phenomena of difficult to understand from the experimental measurements alone. Moreover, developing equivalent physics-based electronic models with both experimental and computational outputs can be invaluable, not only from the scientific point of view but also from the inherently pedagogical opportunities the work offers.

Mathematical models were developed by several people including Jin-Ichi Nagumo (1926--1999), /Richard FitzHugh (1922--2007), and Nobel prize winners in Physiology or Medicine 1963 Alan L. Hodgkin (1914--1998) and Andrew F. Huxley (1917--2012) in the early 1950s. Mathematical modeling of neurons has been largely based on equations describing the physics of action potential generation across the cell membrane and the equivalent electronic circuitry. In the case of the Hodgkin and Huxley results, two defining factors made the realization of their model possible: (a) the selection of the squid giant axon, with a diameter about 0.5 mm, large enough for placing two electrodes inside the cell, and (b) the development of the voltage clamp technique allowing the trans-membrane voltage to be held at a set level during the recording of experimental data. This technique, developed by the biophysicist Kenneth C. Cole (1900-1984), enabled Hodgkin and Huxley to obtain the membrane conductance changes, which happen at different membrane potentials, from the measurement of currents flowing across the membrane. They proposed that the electrical behavior of the membrane could be mimicked with three ion channels (leak, potassium, and sodium) and the corresponding circuitry.

Hodgkin and Huxley conjectured that each of these three ionic currents would be determined by a driving force measurable as a voltage difference and a conductance. In this way, the sodium current, for example, would be given by the sodium channel conductance multiplied by the difference between the membrane potential and the sodium equilibrium potential, and likewise for the potassium and leak currents. They demonstrated that the rising of the action potential was mostly associated with an increase in the conductance of sodium channels, while the dropping of the action potential was mostly associated with an increase in the conductance of potassium channels, consistent with the sodium ions flowing into the cell and the potassium ions flowing out of the cell. An important aspect of Hodgkin and Huxley’s work was that the changes they observed in the cell membrane permeability appeared to be dependent on the membrane potential, not on the membrane current. Based on their measurements, they therefore suggested that the electrical activity across the cell membrane could be represented by the electrical circuit, composed basically of a capacitor in parallel with three dynamic resistors (or conductors) along with their corresponding voltage sources.


Zeeman heartbeat model

Example: The mathematical modeling of biological systems has proven to be a valuable tool by allowing experiments which would otherwise be unfeasible in a real situation. The electrocardiogram (ECG) signal is one of the most obvious effects of the human heart operation. The oscillation between systole and diastole states of the heart is reflected in the heart rate. The surface ECG is the recorded potential difference between two electrodes placed on the surface of the skin at pre-defined points. The largest amplitude of a single cycle of the normal ECG is referred to as the R-wave manifesting the depolarization process of the ventricle. The time between successive R-waves has been widely used as a measure of the heart function, and this helps to identify patients at risk for a cardiovascular event or death. Analysis of variations in this time series is known as heart rate variability (HRV) analysis.

The development of a dynamical model will provide a useful tool to analyse the effects of various physiological conditions on the profiles of the ECG. The model-generated ECG signals with various characteristics can also be used as signal sources for the assessment of diagnostic ECG signal processing devices. The dynamic response of the cardiovascular control system to physiological changes is reflected in HRV and blood pressure. Mathematical models are extremely important for understanding biological processes. Now catastrophe theory is one branch of applied mathematics that was developed in order to describe certain biological processes and has been applied by researchers, especially by Christopher Zeeman, who in 1972 resented an important set of nonlinear dynamical equations for heartbeat modelling:

\[ \begin{split} \varepsilon\,\dot{x} &= T\,x - x^3 - y, \\ \dot{y} &= x - x_d ; \end{split} \qquad\quad x(0) =1, \quad y(0) =0. \]
In this model, the variable x represents the length of a muscle fiber in the heart and the variable y represents an electrical control variable that triggers the electro-chemical wave leading to the heart contraction; positive constant T represents a tension of muscle and is related to blood pressure and ε is a small positive constant which plays a role in fast eigenvalues of the system. Parameter xd is the average muscle length in the diastole (the period of time when the heart refills with blood after contraction).

The following properties of this model are considered as fundamental: (i) The existence of an equilibrium state (fixed point) corresponding to the diastole (relaxed state of the heart). (ii) There must be a threshold for triggering the process whereby the heart contracts from a diastole to a systole (fully contracted state and another equilibrium state). (iii) The model should quickly return to the original equilibrium state after the systole.

First, we find equilibrium solution

\[ x = x_d , \quad y = x_d^3 - T\, x_d , \]
where xd is a positive parameter. The Jacobian of Zeeman system is
\[ {\bf J}(x,y) = \begin{bmatrix} \frac{1}{\varepsilon} \left( T - 3x^2 \right) & -\frac{1}{\varepsilon} \\ 1& 0 \end{bmatrix} . \]
The condition for the real part of the eigenvalue to be negative is \( 3\,x^2 -Y > 0 . \) So the system can be stable if
\[ x \ge \sqrt{T/3} \qquad\mbox{and} \qquad x \le \sqrt{T/3} . \]
These conditions can be satisfied if the value of xd be 1.024. Then the equilibrium point becomes (1.024, -0.0497) in the state space, for parameters T = 1 and ε = 0.2. In this case, all of the trajectories, irrespective to their primary conditions, go to the diastolic equilibrium point.

The system stays at the stable equilibrium point endlessly, until the equilibrium point is stable, except there is an exterior excitation that forces the system move to a new equilibrium point. Based on this new stable system, a control input is added to the Zeeman system as shown below:

\[ \begin{split} \varepsilon\,\dot{x} &= T\,x - x^3 - y, \\ \dot{y} &= \left( x - x_d \right) + \left( x_d - x_s \right) u(t) , \end{split} \]
where, when heart is in systolic state, xs is an additional parameter representing a typical fiber length and u(t) shows the cardiac pacemaker controls mechanism that leads heart into diastolic and systolic states. The feedback controller can be found in the form
\[ u(t) = -k\,x(t- \tau ) , \]
where τ is the time-delay and is nonnegative.    ■


FitzHugh--Nagumo model


Hodgkin and Huxley model


  1. Jafarnia-Dabanloo,N., D.C. McLernon, D.C., Zhang, H., Ayatollahi, A., Johari-Majd, V., , Journal of Theoretical Biology, 2007, 244, pp. 180--189.
  2. Jardóon-Kojakhmetova, H. and Broer, H.W., Polynomial normal forms of constrained differential equations withthree parameters, 2014.
  3. Jones, D.S., Sleeman, B.D., Differential Equations and Mathematical Biology, Chapman & Hall. London. (2003).
  4. McSharry, P.E., Clifford, G.D., Tarassenko, L., Smith, L.A., A Dynamical Model for Generating Synthetic Electrocardiogram Signals, IEEE Transactions on Biomedical Engineering, 2003, 50, pp. 289--294.
  5. Pandiselvi, A. Saravanakumar, R., Sivasankari, M.K., Numerical Solution of Zeeman Heartbeat Systems: A Nonlinear Model, International Journal of Scientific Research and Review,
  6. Rauch, J. and Smoller, J., Qualitative theory of the FitzHugh-Nagumo equations, Advances in Mathematics, 1978, Volume 27, Issue 1, January 1978, Pages 12-44;
  7. Zeeman, E.C., Differential Equations and Nerve Impulse, Towards a Theoretical Biology, 1972, 4, pp. 8-67.


Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions