Print out values from python / matlab implementations and compare with values from NeuroML2

[slarson]

Moved from a follow on to #13

@pgleeson wrote:

It would be great to get another Python script to print out the corresponding values in the python/matlab impls vs the values used for NeuroML2. I've started a script for this here: https://github.com/openworm/muscle_model/blob/master/BoyleCohen2008/PythonSupport/Main_Version/compareToNeuroML2.py and it would be great if someone could look into this a bit more. Eventually the nml2 can be updated to make sure all the values match.

[slarson]

@pgleeson looking a bit more into your commit, it sounds like you did this, yes?

[pgleeson]

There is still the issue of the factor which changes the calcium current into changes in [Ca+], thiCa. Related to @VahidGh's comment here: https://gitter.im/openworm/muscle_model on 27 Jan. Still looking into this...

[VahidGh]

As for this issue, I guess this line should be replaced with nml2_thiCa = ca_pool.restingConc / nml2_t_ca * get_si_value(densities['ca_boyle_all']) * nml2_Cmem .
But my comment, was relating to another problem with inactivation parameter combination in the original input.csv file (which I think is due to incorrect min/max selection here or maybe the fitness function implementation). my reasons are:

1. According to this study, in which the author has reviewed over 65 studies on l-type calcium channels, the mean half voltage is -38.95 for inactivation and -21 for activation (table 3. page 28), although it is discussed that these values are less in cardiac cells but I couldn't find any half-voltage inactivation value greater than the activation one, where as in the input.csv we have -0.3 mv for activation and 25 mv for inactivation!

2. If you remove VDI/CDI inactivation expressions (f(j,1)*(1 + (h(j,1) - 1)*alphaCa)) from Vclamp and Iclamp implementations, the final output would be the same, but that's not the case for potassium current calculations (although, I guess removing inactivation expressions by the author, intentionally from here, and here could be for generating the same I/V plot for calcium current as the one from this paper)!

So my argument is, if the CDI/VDI in calcium current, should be effective in the final output (by some magnitude), why by removing these equations, nothing happens?! and if they are effective but the combination of values for inactivation parameters are in a way that causing ineffectiveness of the equations (causing the resultant of ~1 for calcium inactivation expressions), how can we correct the model? or better saying, is there any way to validate a model not only based on the closeness of the final output to experiments, but the validity of each part (e.g. inactivation effect here) and the final output in compare with experiments?!

P.S. there are two typos in the Boyle & Cohen paper, one with k_f which is written as +5, that should be -5 mv. (although, as you know, in such studies it's common to use positive values for k, and instead using the minus in the inactivation expression, but as there is no difference between activation/inactivation expressions in this model, the authors decide to use minus for all inactivation k values.)
another typo is in k_h which is -10 µM, but should be -10 nM.
both of these values are correct in the input.csv

[VahidGh]

Sorry, I didn't try to test my suggestion for this issue when I was commenting.
When I replace this line nml2_thiCa = get_si_value(ca_pool.restingConc) / (nml2_t_ca * get_si_value(densities['ca_boyle_all']) * nml2_Cmem * 100) for calculation of thiCa I get this error:
AttributeError: 'FixedFactorConcentrationModel' object has no attribute 'restingConc'
I guess we should add this attribute to the fixedFactorConcentrationModel, at the first place and then replacing 6.41889e-8 mM with 6.1e-6 M in this file.
Assuming I could read this value, running the code by adding this: nml2_thiCa = 6.1e-6 / (nml2_t_ca * get_si_value(densities['ca_boyle_all']) * nml2_Cmem * 100), yields thiCa=33015.6519535

[VahidGh]

@pgleeson, Is there a way that I can solve the error I mentioned by myself? or this requires some change in the NeuroML2 format?
I remember we had another problem for customized HH gates for Ca channels.
I wanted to know how can I solve this kind of issues to not bother you for every little problem related to NeuroML.

[pgleeson]

Regarding the issue of the inactivation curves not having too much of an effect on the ca current due to the specific parameters shown, these plots show the neuron/muscle example with (left) and without (right) inactivation (both of gates f and h removed from the ca_boyle.nml):

It's a small difference, but there are changes in the currents vars and the ca conc. I'd say the best approach is to try to reach the point that we're sure the NeuroML2 version is as close as possible to the original version (using input.csv), and then retune/refine the models (constraining as best we can each individual part) to produce a better version for the network model.

Looking at the nml2_thiCa issue...

[pgleeson]

Regarding the thiCa value, I've updated the rho value in the nml2 file in this commit: a65e674, and now the values for thiCa match in compareToNeuroML2.py.

I believe thiCa should just be the rho factor divided by the surface area, see http://www.neuroml.org/NeuroML2CoreTypes/Cells.html#fixedFactorConcentrationModel. Can you double check this @VahidGh?

restingConc isn't really a factor in this calculation. This refers to the default ca conc the "pool" will tend to when no ca current is flowing. In the original model the resting conc is 0. I just put in a non zero value originally since usually a zero here can lead to infinities, but thankfully not in this case (the reversal potential of ca channels is fixed rather than calcualted from internal/external conc via the nernst eqn).

[VahidGh]

@pgleeson, thanks for your explanations. For the thiCa issue, my suggestion was based on this calculation, in which I guess 6 mM is the Ca concentration in the external solution (and not the resting concentration as you suggested), I guess this value probably coming from this experiment, and if the calculation is true, so dividing by g_ca * t_ca should have the meaning of the dependency not only to the surface area, but also to the conductance of the ion channel and the decaying time. But perhaps it could be calculated in other ways as you suggested, and having the same result :)
For the inactivation issue, I checked this paper again, and the author mentions that for the peak current value (in 10 mv of potential), we have the max CDI (~0.5) and I guess I should check if by removing this expression we'll have the twice of the peak current value for the Calcium, to be sure of the validity of inactivation expressions (of course, in the next phase as you suggested).