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 can be very hard to make both 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 function, based on a key function, hh.m.

Before we discuss the key function the nextState equation requires explanation:

nextState = previousState + oneDT * f(previousState)

This is a famous equation known as Euler's method to approximate solutions to the differential equation:

d(state)/dt = f(state)

Lets look at Euler's method next.

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 NEURON code for visualizing the numerical solution of differential equations. Find the file icon in the SaveFilesHere directory and double click on it.

Creating Models From The Key and Template Files

Download:

• HH_step_template.m: This is the function that steps from one time to the next.

• HH_meas_template.m: This is the function that models the measurement of the membrane potential with noise.

• HH.m This is the key file.

• CreateCell.m: This it the function that converts the key and template files into the MATLAB code files.

Now type

createCell(hh)

This creates five files:

• HH_step_gen.m

• HH_step_fit.m

• HH_meas_gen.m

• HH_meas_fit.m

• HH_ics.m

If you change the template or the key you get a different set of 5 functions. That's the whole point. The "gen" functions are for generating data. The "fit" functions are for fitting the model to data. The point is the model used to generate data can be different from the one used for fitting. The "ics" function returns the initial conditions (parameters and state).