Targil 3

1. Finding Roots Algorithms

2. Ordinary Differential Equations

3. Using the Runge-Kutta Method to solve the Hodgkin-Huxley Model
The Problem
I choose to describe a problem discussed in an undergraduate course I have taken in the EE department - Introduction to Biological Signals and Systems.
In physiology, an action potential is a short-lasting event in which the electrical membrane potential of a cell rapidly rises and falls, following a consistent trajectory. The shape of a typical action potential is depicted in the following figure. The membrane potential remains near a baseline level until at some point in time it abruptly spikes upward and then rapidly falls.
action potential
The Hodgkin–Huxley model, or conductance-based model, is a mathematical model that describes how action potentials in neurons are initiated and propagated, through a system of four first-order non-linear differential equations:
HH model equations
where V is the membrane potential, I is the externally applied current. VNa,VK,VCl are the ions' resting potential. The current of each ions type is proportional to the potential difference between the membrane potential and the rest potential. The parameters gikk and the rest potentials are dependent on the neuron type and are found experimentally. Following is an example set of parameters:
HH parameters
We wish to solve the above set of equations for V(t), given a set of initial conditions and a given externally applied current.
I would chose to solve this problem using the fourth order Runge-Kutta method. This method is chosen because it is more accurate than the Euler method and since the equations are not linear, the implicit methods can be difficult. The steps for solving would be the following:
  1. Substitute the initial V(1) to calculate α(V(1)),β(V(1)) and substitute to get m(1),n(1),h(1).
  2. For each time step i: we shall define ni=dn/dt(i) and the same for mi and hi and calculate
    • k1=dt*HH(i,V(i),n(i),m(i),h(i)) where HH is the differential equation for dV/dt.
    • k2=dt*HH(i+dt/2,V(i)+k1/2,n(i)+n1/2,m(i)+m1/2,h(i)+h1/2)
    • k3=dt*HH(i+dt/2),V(i)+k2/2,n(i)+n2/2,m(i)+m2,h(i)+h2)
    • k4=dt*HH(i+dt),V(i)+k3,n(i)+n3,m(i)+m3,h(i)+h3)
  3. Create the next step for each variable:
    • V(i+1)=V(i)+1/6*(k1+2*k2+2*k3+k4)
    • n(i+1)=n(i)+1/6*(n1+2*n2+2*n3+n4)
    • m(i+1)=m(i)+1/6*(m1+2*m2+2*m3+m4)
    • h(i+1)=h(i)+1/6*(h1+2*h2+2*h3+h4)
I would use MATLAB since it is the easiest and the language I am used to work with, but it can be easily implemented in any language.
For a given set of initial conditions and not too long a time vector, this should not be too long of a calculation, and I would plot V(t) at the end. But, if we wish to examine the response to changing external currents for a long period of time, I might choose to plot the progress of V(t) at every iteration.
Accuracy and Validation
The accuracy is the order of dt4, so the smaller dt the better. But, as we've learned, this can be a problem if dt is smaller than the precision of the arithmetic in the computer.
Another problem might occur if the initial conditions and external current are such that will create chaos. Reduction of the above equations to two dimensions under certain assumptions can be examined in state-space representations and there isn't always a steady state. In real life this cannot happen (since n and h bring V down) but if the given conditions are not realistic one might receive chaos. This should not happen.