Group Abstract Group Abstract

Message Boards Message Boards

Hodgkin-Huxley cable equation or propogation equation

GROUPS:
I am trying to find the equations and the solution to an equation that seems to have several names: Hodgkin-Huxley Cable Equation, Propogation Equation, Wave Equation.  I am looking for the equation and solution to the action potential as it goes down an axon in Mathematica code.

I am stuck with being in a wheel chair and tied to oxygen thus I cannot get to a library.  I no longer work as a professor so I have no one to ask.  So I am stuck. 

I would appreciate any help
Thank you
Jake
POSTED BY: Jake Trexel
Answer
1 year ago
When I Google for Mathematica Hodgkin Huxley Cable Equation I see several results that might help you.

This includes this link: Modelling Hodgkin-Huxley which provides some explanation and an appendix with a Mathematica implementation.
POSTED BY: Bill Simpson
Answer
1 year ago
There is also a very extensive application that simulates this equation at the Wolfram Demonstrations Project:

The Hodgkin-Huxley Equations for Transmission of Electrical Impulses along an Axon

It of course is written in Mathematica and the source code can be downloaded there via a link on the right called "Download Author Code".

POSTED BY: Darya Aleinikava
Answer
1 year ago
Thank you Darya
That is an excellent program that you wrote.  Bill helped me find the equations that I was looking for.

Thanks
Jake
POSTED BY: Jake Trexel
Answer
1 year ago
Thank you.  I looked at the code to what you posted, but I don't see any equations.  Maybe I am missing the point of your code or how it is written, and if so I am sorry.  I exepcted to see a partial DE.  Can you tell me more about what you did?  Its is an excellent program, I am just trying to understand it.

Thanks
Jake
POSTED BY: Jake Trexel
Answer
1 year ago
If you were looking at the Wolfram Demonstrations Project page he mentioned and
you clicked on the "Download Author Code" on the right side of the screen on that page and
you opened the notebook that it downloads and
you position your mouse exactly on the outermost thin blue bracket to the far right side of "Initialization Code" and
give it a quick left mouse button double click then you should see that this reveals the code they have hidden.

Most of the way to the bottom of that cell in the notebook is the hidden NDSolve[] I think you are looking for.
POSTED BY: Bill Simpson
Answer
1 year ago
Thank you Bill.  Thats what I was looking for.  Thanks for your help.

Jake
POSTED BY: Jake Trexel
Answer
1 year ago
Jake, you are most welcome, I am glad to help. Just a note - it was not me who wrote the program, - the author's name is Pedro Faria. That site you downloaded the program from is a publication site hosted by Wolfram Research. There are many people who publish their programs for free there. You can search that site for more relevant publications. Here is what I found - click on this link: Demonstrations Search on Hodgkin-Huxley

Bill thaks for explaining how to get the source code.

Darya
POSTED BY: Darya Aleinikava
Answer
1 year ago
I would like to use the code that was used in the demonstrations in my book.  Is this allowed??  I would be more than happy to give the author the credit for the program.  My book is titled "The Dynamic Neuron" and it will be a free book.  It is written using Mathematica so the code will fit right in.  I just dont have what I need to write my own code for this part of my book.  

Please let me know if I have the right to do this or not. 

Thanks
Jake
POSTED BY: Jake Trexel
Answer
11 months ago
I need more help with the program, if you dont mind.  I am no expert in Mathematica and I having problems.  The example you sent me was excellent, but could break it down to something much more simple.  All I would like it to do is plot the membrane potential over time the very first plot or the one one in blue.  I am not even interested in having a manipulation of the program.  Just the solution to the DE and a plot of it    I am not seperate out what I need from the example.  I got lost.  I sould like a school kid, but I am out of my field and not know Mathematica, my mental processor is over loaded with too much information.

I sure hope you are ok with me asking more from you.

Thanks
Jake
POSTED BY: Jake Trexel
Answer
1 year ago
This strips out the Manipulate and some of the code for plotting other than the membrane potential from the demonstration.
What remains is initialization of variables defining the initial conditions,
defining some utility functions to make other parts of the code a little more compact,
initialization of variables controlling what stimulus to apply and what range to plot,
numerically solving the DE (so there is no closed form algebraic solution you can see)
and finally plotting the potential versus time.

It might be possible to remove a little more from this, but I hesitate to take any more risk breaking something.

Hopefully you can "scrape" the contents of the box below, paste into your Mathematuca notebook, evaluate and see the resulting plot.

If there is anything else anyone can do to try to make this more useful for you then feel free to ask.

 (*Values for the physical quantities relevant in the experiment performed by Hodgkin and Huxley with the squid giant axon. \
 Conductances are given in mho/cm^2, voltages in mV*)
 GNa = 120.0*10^-3; (*Sodium Maximum Conductance*)
 GK = 36*10^-3; (*Potassium Maximum Conductance*)
 GL = 0.3*10^-3; (*Leakage Current Maximum Conductance*)
 VNa = 115; (*Sodium Diffusion Potential*)
 VK = -12; (*Potassium Diffusion Potential*)
 VL = 11; (*Leakage Current Diffusion Potential*)
 a = 238*10^-4; (*Radius of the axon, given in cm*)
Cap = 10^-3; (*Capacitance per unit area of the axon membrane, given in mF/cm^2*)
r = 2*10^4; (*Longitudinal resistance per unit length of the axon, given in ohm/cm*)
m0 = 0.053;
h0 = 0.59;
n0 = 0.318;
V0 = 0;

(* Define utility functions *)

(*Voltage-dependent coefficients to be used in the equations for the dynamics of gate-variables m,h and n*)
am[V_] := (0.1 (25 - V))/(E^((25 - V)/10) - 1);
bm[V_] := 4 E^(-(V/18));
ah[V_] := 0.07 E^(-(V/20));
bh[V_] := 1/(E^((30 - V)/10) + 1);
an[V_] := (0.01 (10 - V))/(E^((10 - V)/10) - 1);
bn[V_] := 0.125 E^(-(V/80));
(*Equations for the dynamics of gate-variables m,h and n. dm, dh and dn stand for dm/dt, dh/dt and dn/dt, respectively*)
dm[t_] := am[V[t]] (1 - m[t]) - bm[V[t]] m[t];
dh[t_] := ah[V[t]] (1 - h[t]) - bh[V[t]] h[t];
dn[t_] := an[V[t]] (1 - n[t]) - bn[V[t]] n[t];
(*Aperture of the sodium and potassium channels*)
ChannelsNa[u_] := (m[u])^3 h[u];
ChannelsK[u_] := (n[u])^4;
(*Current density of sodium, potassium and leakage current, together with the resultant ion current density*)
JNa[u_] := GNa ChannelsNa[u] (V[u] - VNa);
JK[u_] := GK ChannelsK[u] (V[u] - VK);
JL[u_] := GL (V[u] - VL);
Jion[u_] := JNa[u] + JK[u] + JL[u];
(*UPS stands for 'Unique Point Solution', since it finds and plots a numerical solution for the behaviour of voltage, currents and gates in a chosen point along the axon. J0 is the stimulus current,dV is the function dV/dt and s is the numerical solution returned by the Mathematica function NDSolve*)

(* Initialize stimulus and plot range values *)

stimulus = 0.02;(* stimulus strength (mA/cm^2) *)
start = 10;(* start of stimulus (ms) *)
duration = 40;(* stimulus duration (ms) *)
domaini = 0;(* initial time *)
domainf = 50;(* final time *)
rangel = -20;(* minimum range *)
rangeh = 120;(* maximum range *)

(* Numerically solve the DE and plot the potential versus time *)

Block[{J0, dV, s},
J0 = Function[t, If[start <= t <= duration, stimulus, 0]];
dV = Function[t, (J0[t] - Jion[t])/Cap];
s = NDSolve[{
    V'[t] == dV[t],
    m'[t] == dm[t],
    h'[t] == dh[t],
    n'[t] == dn[t],
    V[0] == V0,
    m[0] == m0,
    h[0] == h0,
    n[0] == n0},
   {V, m, h, n},
   {t, domainf}
   ];
Plot[
  Evaluate[{V[t]} /. s],
  {t, domaini, domainf},
  PlotRange -> {rangel, rangeh},
  AxesLabel -> {"time (ms)", "membrane voltage (mV)"},
  AxesOrigin -> {domaini, rangel},
  GridLines -> Automatic, GridLinesStyle -> Dashed,
  ImageSize -> 450, PlotStyle -> {Blue}]
]

POSTED BY: Bill Simpson
Answer
1 year ago
Dear Bill
That reduced it a lot for me.  I really appreciate you taking the time to reduce it for me. 

thanks
Jake
POSTED BY: Jake Trexel
Answer
1 year ago
Hey there, I have added a thread about this simulation using Wolfram SystemModeler. Please take a look and let me know what I should add.
POSTED BY: Shenghui Yang
Answer
1 year ago
I am very impressed with what you have done Shenghui.  The best part I liked about the program is the what is showed the plots.  I dont have the program, so is was quite a shock to see it, with this problem.

Thanks
Jake
POSTED BY: Jake Trexel
Answer
1 year ago
Jake, I am glad that you like it. The SystemModeler is very powerful. Yet I still did not show everything in the thread that how you call this model from within Mathemtica. This would be another thread or topic. 
POSTED BY: Shenghui Yang
Answer
1 year ago
Shenghui
Is it easier to make models using the modler, or doing it with just Mathematica?  Or do you need both?

Jake
POSTED BY: Jake Trexel
Answer
1 year ago
I always pick the right tool to do my job.

You can see that Mathematica's Manipuate is very handy to demonstrate the basic principle of the system. But you may have to reprogram many things if you are asked to show how two neurons respond to the signal if they are connected in parallel or in series. Moreover, You cannot change the input source without NOT looking into the Mathemtica code.

Now comes the unique strength of modeling language that you can find in SystemModeler - details are in this post:

A Hands-on Guide to Simulating a Neuron Cable with Wolfram SystemModeler

I can change input by set the parameter in the "pulse" via GUI also I can directly link the voltage output in the model to a magnifying circult compoent or a monitor, just like you put a brain on plate and connect wires to it in a lab (you won't fry the brain here if you accidentally tune the wrong button).  For instance, In the attached picture, I simply drag and drop a build-in low-pass filter and link it to cell from the voltage port. I did not touch any thing in the axon model yet I had added a bunch of new stuff to my model. The extendability is seen here clearly. 


POSTED BY: Shenghui Yang
Answer
1 year ago
To respond to the question raised by Jake Trexel, to the best of my knowledge the noncommercial usage you describe is permitted. To check for full details I will direct you to two relevant links below. There may be provisions regarding correct attribution, notice of copyright, and the like so I recommend you read through everything carefully.
I hope this provides useful guidance.
POSTED BY: Daniel Lichtblau
Answer
11 months ago
Why is the current represented as a function vs. a given value? I am trying to solve the H-H equations that he did in this simulation, and have no problems execpt for understanding why Current is not a given value. I want someone to go to the demonstrations page a look up the the program at the Demonstrations site page and look at the code. I also just tired putting in the code. Here is the code, I finally found a way to post it, I hope. The line that is I do not understand is starts with this Istim . Why is not a consant. The way I look at it, is the neuron fired a current (just one current) to the axon.
 GK = 36; GNa = 120; VL = 10.6;
 alpham[V_] := (0.1 (25 - V))/(Exp[(25 - V)/10] - 1);
 betam[V_] := 4 Exp[-V/18];
 alphan[V_] := (0.01 (10 - V))/(Exp[(10 - V)/10] - 1);
 betan[V_] := 0.125 Exp[-V/80];
 alphah[V_] := 0.07 Exp[-V/20];
 betah[V_] := 1/(Exp[(30 - V)/10] + 1);
 gNa[t_] := ((m[t])^3) h[t];
 gK[t_] := (n[t])^4;
Istim[t_, duration_, strength_] := If[0 < t < 100, strength, 0];


Manipulate[
Show[
  Plot[
   Evaluate[
    V[t] /. NDSolve[{V'[t] ==
        1/Cm (Istim[t, duration,
            strength] - ((GNa) gNa[t] (V[t] - VNa)) - ((GK) gK[
              t] (V[t] - VK)) - ((GL) (V[t] - VL))),
       m'[t] == (alpham[V[t]] (1 - m[t])) - (betam[V[t]] m[t]),
       n'[t] == (alphan[V[t]] (1 - n[t])) - (betan[V[t]] n[t]),
       h'[t] == (alphah[V[t]] (1 - h[t])) - (betah[V[t]] h[t]),
       V[0] == -0.0741, m[0] == 0.0525, n[0] == 0.3183,
       h[0] == 0.588}, {V, m, n, h}, {t, 0, 80},
      MaxSteps -> 100000]],
   {t, 0, 60}, PlotRange -> {-20, (VNa + 10)}, AspectRatio -> 1,
   PlotLabel -> "membrane voltage versus time",
   AxesLabel -> {"time (msec)", "membrane voltage (mV)"},
   ImageSize -> {550, 300}, PlotStyle -> {Thickness[.015], Blue}],
  Graphics[{Thickness[0.008], Dashed, Red,
    Line[{{0, VNa}, {60, VNa}}], Line[{{0, VK}, {60, VK}}]}]],
{{strength, 6.3, "applied current"}, 0, 18,
  Appearance -> "Labeled"},
{{VNa, 115, "sodium equilibrium potential"}, 110, 150,
  Appearance -> "Labeled"},
{{VK, -12, "potassium equilibrium potential"}, -15, 5,
  Appearance -> "Labeled"},
Delimiter,
{{Cm, 1, "membrane capacitance"}, {1, 10^-2.99}, Checkbox},
{{GL, 0.3, "leakage channel conductance"}, {0.3, 0}, Checkbox},
SaveDefinitions -> True]
POSTED BY: Jake Trexel
Answer
11 months ago