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 at henning.org
===================================

	[[alternative HTML version deleted]]