Difference between revisions of "NumHH"

From Sean_Carver
Jump to: navigation, search
Line 28: Line 28:
  
 
== Numerical Solution of Differential Equations ==
 
== Numerical Solution of Differential Equations ==
 +
 +
 +
{|
 +
| valign="top" |<source lang="matlab">[X,Y] = meshgrid(-10:0.25:10,-10:0.25:10);
 +
f = sinc(sqrt((X/pi).^2+(Y/pi).^2));
 +
mesh(X,Y,f);
 +
axis([-10 10 -10 10 -0.3 1])
 +
xlabel('{\bfx}')
 +
ylabel('{\bfy}')
 +
zlabel('{\bfsinc} ({\bfR})')
 +
hidden off
 +
</source>
 +
| &nbsp;&nbsp;&nbsp;
 +
| valign="top" |<source lang="matlab">
 +
[X,Y] = meshgrid(-10:0.25:10,-10:0.25:10);
 +
f = sinc(sqrt((X/pi).^2+(Y/pi).^2));
 +
surf(X,Y,f);
 +
axis([-10 10 -10 10 -0.3 1])
 +
xlabel('{\bfx}')
 +
ylabel('{\bfy}')
 +
zlabel('{\bfsinc} ({\bfR})')
 +
</source>
 +
|-
 +
| This code produces a '''[[wire frame model|wireframe]]''' 3D plot of the two-dimensional unnormalized [[sinc function]]:
 +
| &nbsp;&nbsp;&nbsp;
 +
| This code produces a '''surface''' 3D plot of the two-dimensional unnormalized [[sinc function]]:
 +
|-
 +
| align="center" |[[Image:MATLAB mesh sinc3D.svg]]
 +
| &nbsp;&nbsp;&nbsp;
 +
| align="center" |[[Image:MATLAB surf sinc3D.svg]]
 +
|}
  
 
Remember the equation for the cell with only leak channels.
 
Remember the equation for the cell with only leak channels.

Revision as of 20:08, 23 February 2009

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): 

<source lang="matlab">
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)];
</source>


Numerical Solution of Differential Equations

<source lang="matlab">[X,Y] = meshgrid(-10:0.25:10,-10:0.25:10);

f = sinc(sqrt((X/pi).^2+(Y/pi).^2)); mesh(X,Y,f); axis([-10 10 -10 10 -0.3 1]) xlabel('{\bfx}') ylabel('{\bfy}') zlabel('{\bfsinc} ({\bfR})') hidden off </source>

    <source lang="matlab">

[X,Y] = meshgrid(-10:0.25:10,-10:0.25:10); f = sinc(sqrt((X/pi).^2+(Y/pi).^2)); surf(X,Y,f); axis([-10 10 -10 10 -0.3 1]) xlabel('{\bfx}') ylabel('{\bfy}') zlabel('{\bfsinc} ({\bfR})') </source>

This code produces a wireframe 3D plot of the two-dimensional unnormalized sinc function:     This code produces a surface 3D plot of the two-dimensional unnormalized sinc function:
File:MATLAB mesh sinc3D.svg     File:MATLAB surf sinc3D.svg

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.

Retrieved from "https://seancarver.org/index.php?title=NumHH&oldid=450"