thr3ads.net - R help - [R] R and MatLab implementations of the same model differs [Jul 2013]

If this information is useful, please help other people find it:
Share via:

Jannetta Steyn

2013-Jul-04 15:15 UTC

[R] R and MatLab implementations of the same model differs

Hi folks

I have implemented a model of a neuron using Hodgkin Huxley equations in
both R and MatLab. My preference is to work with R but R is not giving me
the correct results. I also can't use ode45 as it just seems to go into an
indefinite loop. However, the MatLab implementation work fine with the same
equations, parameters and initial values and ode45. Below is my R and
MatLab implementations.

I'd really appreciate any help.

R
=
rm(list=ls())

library(deSolve)

ST <-  function(time, init, parms) {
  with(as.list(c(init, parms)),{

    #functions to calculate activation m and inactivation h of the currents
    #Axon
    mNax <- function(v) 1/(1+exp(-(v+24.7)/5.29))
    taumNa <- function(v) 1.32-(1.26/(1+exp(-(v+120)/25)))
    hNax <- function(v) 1/(1+exp((v+48.9)/5.18))
    tauhNa <- function(v)
0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6))));
    mKx <- function(v) 1/(1+exp(-(v+14.2)/11.8))
    taumK <- function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))


    # Currents as product of maximal conducatance(g), activation(m) and
inactivation(h)
    # Driving force (v-E) where E is the reversal potential of the
particular ion

    # AB axon
    iNa_axon_AB <-
gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB)
    iK_axon_AB <- gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB)
    iLeak_axon_AB <- gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)

    # AB Axon equations
    dv_axon_AB <- (0 -iNa_axon_AB -iK_axon_AB -iLeak_axon_AB)/C_axon_AB
    dmNa_axon_AB <- (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB)
    dhNa_axon_AB <- (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB)
    dmK_axon_AB <- (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB)

    list(c(dv_axon_AB,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB
))

  })}
## Set initial state
init = c(v_axon_AB=-55,mNa_axon_AB=1,hNa_axon_AB=0,mK_axon_AB=1
)
## Set parameters
# AB

gNa_axon_AB=300e-3
gK_axon_AB=52.5-3
gLeak_axon_AB=0.0018e-3

# AB Axon Reversal potentials
ENa_axon_AB=50
EK_axon_AB=-80
ELeak_axon_AB=-60

C_axon_AB=1.5e-3


I=0

parms
c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB)
## Set integrations times
times = seq(from=0, to=5000, by = 0.05);

out<-ode(y=init, times=times, func=ST, parms=parms, method="lsoda")
plot(out)
o<-data.frame(out)

plot(o$v_axon_AB,type='l')


MatLab
=====
File: init.m
clear all
close all
clc

global c_axon_AB ...
    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
    ENa_axon_AB EK_axon_AB ELeak_axon_AB I

T_MAX = 5000; % ms
step = 0.05;

gNa_axon_AB=300e-3;
gK_axon_AB=52.5-3;
gLeak_axon_AB=0.0018e-3;

% AB Axon Reversal potentials
ENa_axon_AB=50;
EK_axon_AB=-80;
ELeak_axon_AB=-60;

c_axon_AB=1.5e-3;

I=0;

x0 = [-55, 0, 1, 0];

simulate(T_MAX, step, x0)
% c_axon_AB, gNa_axon_AB, gK_axon_AB, gLeak_axon_AB, ENa_axon_AB,
EK_axon_AB, ELeak_axon_AB, I,



File: simulate.m
===========
function[] = simulate(T_MAX, step,  x0)
%c_axon_AB, gNa_axon_AB, gK_axon_AB, gLeak_axon_AB, ENa_axon_AB,
EK_axon_AB, ELeak_axon_AB, I,
close all
clc
global c_axon_AB ...
    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
    ENa_axon_AB EK_axon_AB ELeak_axon_AB I

tspan=0:step:T_MAX;

%Sx0 = [-55, 0, 1, 0];
% x0 = [v_axon_AB mNa_axon_AB hNa_axon_AB mK_axon_AB]

[t,x] = ode45(@integrate, tspan, x0);

v_axon_AB = x(1);

vars = getVariables(t,x0);


figure
set(gcf,'numbertitle','off','name','AB - Membrane
potential')
plot(t,v_axon_AB)
title('v axon')


function out = mNax(v)
out = 1/(1+exp(-(v+24.7)/5.29));

function out = taumNa(v)
out = 1.32-(1.26/(1+exp(-(v+120)/25)));

function out = hNax(v)
out = 1/(1+exp((v+48.9)/5.18));

function out = tauhNa(v)
out = 0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6))));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function out = mKx(v)
out = 1/(1+exp(-(v+14.2)/11.8));

function out = taumK(v)
out = 7.2-(6.4/(1+exp(-(v+28.3)/19.2)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xdot = integrate(t,x)
global c_axon_AB ...
    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
    ENa_axon_AB EK_axon_AB ELeak_axon_AB I

v_axon_AB = x(1);
mNa_axon_AB = x(2);
hNa_axon_AB = x(3);
mK_axon_AB = x(4);

iNa_axon = gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB);
iK_axon = gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB);
iLeak_axon = ELeak_axon_AB*(v_axon_AB-ELeak_axon_AB);

dv_axon_AB = (0 -iNa_axon -iK_axon -iLeak_axon)/c_axon_AB;

dmNa_axon_AB = (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB);
dhNa_axon_AB = (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB);
dmK_axon_AB = (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB);

xdot = [dv_axon_AB dmNa_axon_AB dhNa_axon_AB dmK_axon_AB]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function out = getVariables(t,x)
global c_axon_AB ...
    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
    ENa_axon_AB EK_axon_AB ELeak_axon_AB I


v_axon_AB = x(1);
mNa_axon_AB = x(2);
hNa_axon_AB = x(3);
mK_axon_AB = x(4);

iNa_axon gNa_axon_AB.*mNa_axon_AB.^3.*hNa_axon_AB.*(v_axon_AB-ENa_axon_AB);
iK_axon = gK_axon_AB.*mK_axon_AB.^4.*(v_axon_AB-EK_axon_AB);
iLeak_axon = gLeak_axon_AB.*(v_axon_AB-ELeak_axon_AB);

out = [iNa_axon iK_axon iLeak_axon ];



Regards
Jannetta

-- 

==================================Web site: http://www.jannetta.com
Email: jannetta@henning.org
==================================
	[[alternative HTML version deleted]]

Berend Hasselman

2013-Jul-04 15:48 UTC

head link

[R] R and MatLab implementations of the same model differs

On 04-07-2013, at 17:15, Jannetta Steyn <jannetta at henning.org> wrote:
> Hi folks
> 
> I have implemented a model of a neuron using Hodgkin Huxley equations in
> both R and MatLab. My preference is to work with R but R is not giving me
> the correct results. I also can't use ode45 as it just seems to go into
an
> indefinite loop. However, the MatLab implementation work fine with the same
> equations, parameters and initial values and ode45. Below is my R and
> MatLab implementations.
> 
No problem in running your R file. Have plot.
(Mac mini Core2Duo 2.66Ghz running R 3.0.1 Patched (2013-06-19 r62992) on Mac OS
X 10.8.4:  16.5 seconds.)

Trying to run your Matlab file in Octave. No success yet because of unavailable
ode45.

Berend
> I'd really appreciate any help.
> 
> R
> => 
> rm(list=ls())
> 
> library(deSolve)
> 
> ST <-  function(time, init, parms) {
>  with(as.list(c(init, parms)),{
> 
>    #functions to calculate activation m and inactivation h of the currents
>    #Axon
>    mNax <- function(v) 1/(1+exp(-(v+24.7)/5.29))
>    taumNa <- function(v) 1.32-(1.26/(1+exp(-(v+120)/25)))
>    hNax <- function(v) 1/(1+exp((v+48.9)/5.18))
>    tauhNa <- function(v)
> 0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6))));
>    mKx <- function(v) 1/(1+exp(-(v+14.2)/11.8))
>    taumK <- function(v) 7.2-(6.4/(1+exp(-(v+28.3)/19.2)))
> 
> 
>    # Currents as product of maximal conducatance(g), activation(m) and
> inactivation(h)
>    # Driving force (v-E) where E is the reversal potential of the
> particular ion
> 
>    # AB axon
>    iNa_axon_AB <-
> gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB)
>    iK_axon_AB <- gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB)
>    iLeak_axon_AB <- gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)
> 
>    # AB Axon equations
>    dv_axon_AB <- (0 -iNa_axon_AB -iK_axon_AB -iLeak_axon_AB)/C_axon_AB
>    dmNa_axon_AB <- (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB)
>    dhNa_axon_AB <- (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB)
>    dmK_axon_AB <- (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB)
> 
>    list(c(dv_axon_AB,dmNa_axon_AB, dhNa_axon_AB, dmK_axon_AB
> ))
> 
>  })}
> ## Set initial state
> init = c(v_axon_AB=-55,mNa_axon_AB=1,hNa_axon_AB=0,mK_axon_AB=1
> )
> ## Set parameters
> # AB
> 
> gNa_axon_AB=300e-3
> gK_axon_AB=52.5-3
> gLeak_axon_AB=0.0018e-3
> 
> # AB Axon Reversal potentials
> ENa_axon_AB=50
> EK_axon_AB=-80
> ELeak_axon_AB=-60
> 
> C_axon_AB=1.5e-3
> 
> 
> I=0
> 
> parms >
c(ENa_axon_AB,EK_axon_AB,ELeak_axon_AB,gNa_axon_AB,gK_axon_AB,gLeak_axon_AB,C_axon_AB)
> ## Set integrations times
> times = seq(from=0, to=5000, by = 0.05);
> 
> out<-ode(y=init, times=times, func=ST, parms=parms,
method="lsoda")
> plot(out)
> o<-data.frame(out)
> 
> plot(o$v_axon_AB,type='l')
> 
> 
> MatLab
> =====> 
> File: init.m
> clear all
> close all
> clc
> 
> global c_axon_AB ...
>    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
>    ENa_axon_AB EK_axon_AB ELeak_axon_AB I
> 
> T_MAX = 5000; % ms
> step = 0.05;
> 
> gNa_axon_AB=300e-3;
> gK_axon_AB=52.5-3;
> gLeak_axon_AB=0.0018e-3;
> 
> % AB Axon Reversal potentials
> ENa_axon_AB=50;
> EK_axon_AB=-80;
> ELeak_axon_AB=-60;
> 
> c_axon_AB=1.5e-3;
> 
> I=0;
> 
> x0 = [-55, 0, 1, 0];
> 
> simulate(T_MAX, step, x0)
> % c_axon_AB, gNa_axon_AB, gK_axon_AB, gLeak_axon_AB, ENa_axon_AB,
> EK_axon_AB, ELeak_axon_AB, I,
> 
> 
> 
> File: simulate.m
> ===========> 
> function[] = simulate(T_MAX, step,  x0)
> %c_axon_AB, gNa_axon_AB, gK_axon_AB, gLeak_axon_AB, ENa_axon_AB,
> EK_axon_AB, ELeak_axon_AB, I,
> close all
> clc
> global c_axon_AB ...
>    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
>    ENa_axon_AB EK_axon_AB ELeak_axon_AB I
> 
> tspan=0:step:T_MAX;
> 
> %Sx0 = [-55, 0, 1, 0];
> % x0 = [v_axon_AB mNa_axon_AB hNa_axon_AB mK_axon_AB]
> 
> [t,x] = ode45(@integrate, tspan, x0);
> 
> v_axon_AB = x(1);
> 
> vars = getVariables(t,x0);
> 
> 
> figure
> set(gcf,'numbertitle','off','name','AB -
Membrane potential')
> plot(t,v_axon_AB)
> title('v axon')
> 
> 
> function out = mNax(v)
> out = 1/(1+exp(-(v+24.7)/5.29));
> 
> function out = taumNa(v)
> out = 1.32-(1.26/(1+exp(-(v+120)/25)));
> 
> function out = hNax(v)
> out = 1/(1+exp((v+48.9)/5.18));
> 
> function out = tauhNa(v)
> out = 0.67/(1+exp(-(v+62.9)/10))*(1.5+(1/(1+exp((v+34.9)/3.6))));
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> 
> function out = mKx(v)
> out = 1/(1+exp(-(v+14.2)/11.8));
> 
> function out = taumK(v)
> out = 7.2-(6.4/(1+exp(-(v+28.3)/19.2)));
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> 
> function xdot = integrate(t,x)
> global c_axon_AB ...
>    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
>    ENa_axon_AB EK_axon_AB ELeak_axon_AB I
> 
> v_axon_AB = x(1);
> mNa_axon_AB = x(2);
> hNa_axon_AB = x(3);
> mK_axon_AB = x(4);
> 
> iNa_axon = gNa_axon_AB*mNa_axon_AB^3*hNa_axon_AB*(v_axon_AB-ENa_axon_AB);
> iK_axon = gK_axon_AB*mK_axon_AB^4*(v_axon_AB-EK_axon_AB);
> iLeak_axon = ELeak_axon_AB*(v_axon_AB-ELeak_axon_AB);
> 
> dv_axon_AB = (0 -iNa_axon -iK_axon -iLeak_axon)/c_axon_AB;
> 
> dmNa_axon_AB = (mNax(v_axon_AB)-mNa_axon_AB)/taumNa(v_axon_AB);
> dhNa_axon_AB = (hNax(v_axon_AB)-hNa_axon_AB)/tauhNa(v_axon_AB);
> dmK_axon_AB = (mKx(v_axon_AB)-mK_axon_AB)/taumK(v_axon_AB);
> 
> xdot = [dv_axon_AB dmNa_axon_AB dhNa_axon_AB dmK_axon_AB]';
> 
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> 
> function out = getVariables(t,x)
> global c_axon_AB ...
>    gNa_axon_AB gK_axon_AB gLeak_axon_AB ...
>    ENa_axon_AB EK_axon_AB ELeak_axon_AB I
> 
> 
> v_axon_AB = x(1);
> mNa_axon_AB = x(2);
> hNa_axon_AB = x(3);
> mK_axon_AB = x(4);
> 
> iNa_axon >
gNa_axon_AB.*mNa_axon_AB.^3.*hNa_axon_AB.*(v_axon_AB-ENa_axon_AB);
> iK_axon = gK_axon_AB.*mK_axon_AB.^4.*(v_axon_AB-EK_axon_AB);
> iLeak_axon = gLeak_axon_AB.*(v_axon_AB-ELeak_axon_AB);
> 
> out = [iNa_axon iK_axon iLeak_axon ];
> 
> 
> 
> Regards
> Jannetta
> 
> -- 
> 
> ==================================> Web site: http://www.jannetta.com
> Email: jannetta at henning.org
> ==================================> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Berend Hasselman

2013-Jul-05 08:44 UTC

head link

[R] R and MatLab implementations of the same model differs

On 05-07-2013, at 09:53, Jannetta Steyn <jannetta at henning.org> wrote:
> 
> 
> >
> > I don't quite know how to explain the "doesn't work"
in more detail without
> > any visual aid.
> 
> You said that R got into an indefinite loop, whatever that maybe.
> 
> 
> > When I change the solver to ode45 the script never stops running. I
have even left it over night at one stage but the next day it was still busy. Is
there a way to see what it is doing and to determine why it seems to be in an
> > "infinite loop"?
> 
> The script just ran using ode45!! For the first time ever.
> 
Please keep the conversation on R-help.
Don't reply to me personally.

Which script?
With the corrections suggested by me and Peter? (gK_axon_AB=52.5-3 ==>
gK_axon_AB=52.5e-3)

For what it's worth: lsdoda seems quicker.
Variable v_acon_AB now converges to -60 (the equilibrium state?)

Berend

R help - Jul 2013 - R and MatLab implementations of the same model differs

[R] R and MatLab implementations of the same model differs

[R] R and MatLab implementations of the same model differs

[R] R and MatLab implementations of the same model differs