[R] R and MatLab implementations of the same model differs
Berend Hasselman
bhh at xs4all.nl
Thu Jul 4 17:48:50 CEST 2013
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.
More information about the R-help
mailing list