NumHH

We have looked at the action potential with the Hodgkin-Huxley model using the textbook software, Neurons in Action. Today we are going to use MATLAB. We will generate data, then change parameters, then calculate the likelihood that the model (with different parameters) produced the data.

Step Function

Here is the function that simulates the model (one step at a time):

function nextState = hh_step_template(state,input,noise,param)

alpha_m = (2.5 - 0.1*stateV)./(exp(2.5 - 0.1*stateV) - 1);
alpha_n = (0.1 - 0.01*stateV)./(exp(1 - 0.1*stateV) - 1);
alpha_h = 0.07*exp(-stateV./20);
beta_m = 4*exp(-stateV./18);
beta_n = 0.125*exp(-stateV./80);
beta_h = 1./(exp(3-0.1*stateV)+1); 

ionicCurrent = gNa*stateM^3*stateH*(stateV-ENa) ...
     + gK*stateN^4*(stateV-EK) ...
     + gL*(stateV-EL);

nextState = [stateV + oneDT./CAP*(injectedCurrent - ionicCurrent) ...
                    + sqrtDT*PNoiseLevel*PNoise; ...
             stateM + oneDT*(alpha_m*(1-stateM) - beta_m*stateM); ...
             stateN + oneDT*(alpha_n*(1-stateN) - beta_n*stateN); ...
             stateH + oneDT*(alpha_h*(1-stateH) - beta_h*stateH)];

This function returns the nextState (state at current time plus oneDT) as a function of current state (stateV,stateM,stateN,stateH) input (injectedCurrent) noise (PNoise) and parameters (param).

Note that this function won't run because a lot of things (like injectedCurrent and ENa) aren't defined. A function like this one is very hard to both make human readable and run quickly. So I have written preprocessor createCell.m which takes the human readable template function above and creates machine readable MATLAB Code.

Numerical Solution of Differential Equations

Remember the equation for the cell with only leak channels.

$C \frac{dV}{dt} = I(t) - g_L(V - E_L)$

Let's simplify: suppose there is no injected current and that the reversal potential for the leak channels is $E_L = 0$. Then our equation is

$\frac{dV}{dt} = - \frac{g_L}{C} V$

Using different letters for the variables (because this is done in the software linked below):

$\frac{dy}{dt} = - k y$

Here k is the rate constant, 1/k is the time constant, 1/k is $\frac{C}{g_L} = RC$ in the notation above. A leaky cell is what is called an RC circuit -- a resistor and capacitor together in a circuit. The time constant of an RC circuit is RC. The bigger k, the higher the rate of convergence, and the smaller the time constant 1/k. The time constant is the time it takes the solution to decay to 1/e of its value.

Solution of differential equations happens at discrete times: $y_k$ , separated by small time intervals dt.

The simplest way of solving this equation is with Euler's method:

$y_k = y_{k-1} + \Delta t \ (-k y_{k-1})$

This is a special case of the general formula for Euler's method applied to the (vector) differential equation

$\frac{dy}{dt} = f(y)$

$y_k = y_{k-1} + \Delta t \ f(y_{k-1})$

Euler's equation is the simplest way to solve a differential equation numerically. However it is often not the preferred method: often you need to take much smaller time steps with Euler than with some other methods, so it takes longer to get as good a solution. Still if you are doing something complicated, like solving an equation with noise, or Bayesian filtering (to compute likelihood), an argument can be made that a simpler method is desirable -- at least as a first step.

Click here for code for visualizing the numerical solution of differential equations.