R and MatLab implementations of the same model differs

On 04-07-2013, at 19:56, peter dalgaard <pdalgd at gmail.com> wrote:

On Jul 4, 2013, at 19:11 , Berend Hasselman wrote:

On 04-07-2013, at 18:42, Jannetta Steyn <jannetta at henning.org> wrote:

Hi Ben and others

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 you run the two scripts it is easy to see the
difference. MatLab produces a line on x= -55. This is what I expect - a
more or less straight line. R on the other hand the result drops from -55
(the initial value) to -80 and then goes up to -71.37092.

The two scripts have exactly the same equations.

I don't think so.
In the R script you have

init = c(v_axon_AB=-55,mNa_axon_AB=1,hNa_axon_AB=0,mK_axon_AB=1)

That is not the same as in your Matlab script. To make them the same you should replace the line with

init = c(v_axon_AB=-55,mNa_axon_AB=0,hNa_axon_AB=1,mK_axon_AB=0)

Using this line gives quite different results. But not the same as Matlab.

It seems that you ought to carefully check all your equations.
I don't know enough about Matlab/Octave syntax to determine if the results of the two systems  are identical.

Also, the line

iLeak_axon_AB <- gLeak_axon_AB*(v_axon_AB-ELeak_axon_AB)

differs from the Matlab code. Changing it gave me all-NA results, though. And the with() construct assumes that init and parms are named vectors; parms is not, so the values are taken from the global environment. It is not clear whether that makes a difference.

Same NA results here. Making parms a named vector makes no difference but is is certainly necessary to do that.

In that case there is probably an inconsistency in the file simulate.m near the end:
The last non empty line before the line out = ?  reads

iLeak_axon = gLeak_axon_AB.*(v_axon_AB-ELeak_axon_AB);

and this doesn't agree with the line

iLeak_axon = ELeak_axon_AB*(v_axon_AB-ELeak_axon_AB);

in the function xdot.

Final question for the OP: if the model is supposed to be dynamic, isn't it suspicious that Matlab gives a constant result?

Berend

Berend

I have even named the
variables the same in the two scripts and copied the equations across to
make sure I haven't made any typos.

Can one attach images to posts? I'll try. The flatline image is the plot
from MatLab and the other is the plot from R.

Thanks
Jannetta


On 4 July 2013 16:52, Ben Bolker <bbolker at gmail.com> wrote:

Berend Hasselman <bhh <at> xs4all.nl> writes:

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.

I'm impressed that you (BH) went to the trouble of checking
on this vague "doesn't work" question.  If you want to go farther
I think you can get ode45 for octave by installing the odepkg
package:

http://octave.sourceforge.net/odepkg/index.html

-- 

===================================
Web site: http://www.jannetta.com
Email: jannetta at henning.org
===================================
<Rplot01.png><MatPlot01.png>______________________________________________
R-help at r-project.org mailing list

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com