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

[Berend Hasselman]

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
>>>> 
>>>> ______________________________________________
>>>> 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.
>>>> 
>>> 
>>> 
>>> 
>>> -- 
>>> 
>>> ===================================
>>> Web site: http://www.jannetta.com
>>> Email: jannetta at henning.org
>>> ===================================
>>> <Rplot01.png><MatPlot01.png>______________________________________________
>>> 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.
>> 
>> ______________________________________________
>> 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.
> 
> -- 
> 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
> 
> 
> 
> 
> 
> 
> 
>