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

[peter dalgaard]

On Jul 5, 2013, at 00:41 , Jannetta Steyn wrote:

> Hi All
> 
> Thanks for all the comments. I have looked at everything you have pointed
> out and this is the situation at the moment:
> 
>> 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.
> 
> You are right. I have to admit when I referred to equations I referred to
> those in the ST function, and not the initialisation and parameter
> equations. I fixed the line, but look you said it still doesn't produce the
> same results as MatLab.
> 
>> 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.
> 
> The correct equation is:
> 
> iLeak_axon = gLeak_axon_AB.*(v_axon_AB-ELeak_axon_AB);
> 
> I think this error was the result of a variable name change at some point
> to make the to scripts the same. I have fixed the MatLab script but it is
> not making any difference.
> 
>> 
>> Final question for the OP: if the model is supposed to be dynamic, isn't
> it suspicious that Matlab gives a constant result?
> 
> Now this is true ... I'll chat to my MatLab expert about this. I was hoping
> that it is not a constant but that it has small variations in the value
> that I can't see on the graph. However, I did not know how to check the
> values ... another question for my MatLab expert tomorrow.
> 
>> Also, in this block
>> 
>> gNa_axon_AB=300e-3
>> gK_axon_AB=52.5-3
>> gLeak_axon_AB=0.0018e-3
>> 
>> the middle line looks like a typo for 52.5e-3 and the 3rd line looks very
> low -- googling suggests values like (120, 36, 0.3) mS/cm^3, which would be
> more consistent with a value around 1e-3.
> 
> Yes you are right. 52.5-3 should be 52.5e-3. I made the change but still no
> luck.
> 
> With regards to the value of gLeak_axon_AB, the values I'm using for the
> model are from a paper. I tried changing it to what you suggested, but that
> seems to result in a constant value of -60.00006

But isn't -60 = ELeak_axon_AB exactly the right steady-state value for the system with both the K and Na channels closed?

> 
> More late night staring ahead for me ...
> 
> 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.

-- 
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