Understanding NeuroML2 and LEMS

[Robby Simpson]

Hi,

I'm a newbie and wanted to get a better understanding of the underlying physics so that I can contribute to OpenWorm.  To do so, I extracted the NeuroML2 data from the muscle model (https://github.com/openworm/muscle_model/tree/master/NeuroML2) and the LEMS files for the schema (https://github.com/NeuroML/NeuroML2) and created a simple C program to update the voltage for a single muscle cell.

Unfortunately, the voltage varies widely before going to infinity.

I'm hoping that someone can take a look at the attached C code and let me know if I am missing something.

Your help is most appreciated.

Sincerely,
Robby

[Padraig Gleeson]

Hi Robby,

Thanks for your interest in the muscle model and sorry for the delay in getting back to you. The first thing I see from the code is that there is no external current injection, so the cells shouldn't do anything beyond gradually approach resting membrane potential. I'd recommend adding a current clamp input, and testing the behaviour of the cell with just a leak conductance and gradually add the other conductances. You could make an equivalent model in NML2 (based on this) and test the simpler C version against that (running in jNeuroML). I've added your code here: https://github.com/openworm/muscle_model/tree/master/NeuroML2/C. Feel free to fork that and work there. Add a makefile so it's easier for others to test.

If you want to start on an easier to tackle problem, there are 2 implementations of the Hodgkin Huxley model here in Python & NML2: https://github.com/openworm/hodgkin_huxley_tutorial/tree/master/Tutorial/Source, along with a tutorial http://hodgkin-huxley-tutorial.readthedocs.org/en/latest/_static/Tutorial.html. You could make an implementation of this model in C which could eventually be added to this repository. The C muscle cell could be updated afterwards based on your code there.

Regards,
Padraig

--
Visit us online at:
http://openworm.org
http://blog.openworm.org
http://github.com/openworm
http://twitter.com/openworm
https://plus.google.com/s/openworm
---
You received this message because you are subscribed to the Google Groups "OpenWorm-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openworm-discu...@googlegroups.com.
To post to this group, send email to openworm...@googlegroups.com.
Visit this group at https://groups.google.com/group/openworm-discuss.
To view this discussion on the web visit https://groups.google.com/d/msgid/openworm-discuss/2300bc2c-bf00-4e64-8769-7983ec0ac0fa%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

-- 

-----------------------------------------------------
Padraig Gleeson
Room 321, Anatomy Building
Department of Neuroscience, Physiology&  Pharmacology
University College London
Gower Street
London WC1E 6BT
United Kingdom

+44 207 679 3214
p.gl...@ucl.ac.uk
-----------------------------------------------------

[Robby Simpson]

Hi Padraig, All,

Thanks for taking a look and no worries on the delay!

In my code, to make things as simplistic as possible until it is working, I did not include any external currents.  My goal was to make certain the cell should just go to the resting membrane potential before attempting anything else.  Unfortunately, the membrane potential pretty quickly goes to infinity.  Can you tell what I am missing or doing wrong?

For source material, I used:
https://github.com/openworm/muscle_model/blob/master/NeuroML2/CaPool.nml -- this gave me the values for the concentration model (Lines 96-100)
https://github.com/openworm/muscle_model/blob/master/NeuroML2/Leak.channel.nml -- this ultimately goes unused
https://github.com/openworm/muscle_model/blob/master/NeuroML2/ca_boyle.channel.nml -- this gave me the values for the calcium channel (Lines 133 - 153) as well as the dynamics for the customHGate (Lines 176 - 178 and 254 - 256).  Note, the customHGate does not seem to be time-dependent other than depending on the changing calcium concentration.
https://github.com/openworm/muscle_model/blob/master/NeuroML2/k_fast.channel.nml -- this gave me the values for the k_fast channel (Lines 106 - 121)
https://github.com/openworm/muscle_model/blob/master/NeuroML2/k_slow.channel.nml -- this gave me the values for the k_slow channel (Lines 123 - 131)
https://github.com/openworm/muscle_model/blob/master/NeuroML2/SingleCompMuscle.cell.nml -- this gave me the values for the surface area (Line 88), the various condDensity's and erev's, the spike threshold, specificCapacitance, initMembPotential, and the initial calcium concentration (Line 100).

I then used data from the LEMS files from https://github.com/NeuroML/NeuroML2/tree/master/NeuroML2CoreTypes to update derived and time-dependent variables.  This primarily includes:
- inf for the various GateHHTauInf:   value="rate / (1 + exp(0 - (v - midpoint)/scale))"
- iCa
- fopen for the various GateHHTauInf:   value="q^instances"
- deltaConcentration:    value="(iCa/surfaceArea) * rho - ((concentration - restingConc) / decayConstant)"
- deltaQ:    value="(inf - q) / tau"
- deltaV:     value="(iChannels + iSyn) / totCap"

I must admit, I'm at a bit of a loss here...

As a note, the C code I posted has no external dependencies and is standalone.  You should be able to compile it with your favorite C compiler.  For gcc, this would be:
gcc testMuscleOpenworm.c -o testMuscleOpenworm

You should then be able to run the program "testMuscleOpenworm"

Please let me know if you have any questions as I'd love to get to the bottom of this issue.

Thanks,
Robby

[Stephen Larson]

What if you try the simulation with only one channel at a time enabled?  Does the membrane potential always go nuts, or is it one channel, or maybe just a handful?  What about pairwise combinations of channels?  Knowing this will enable you to isolate if it is a problem with your numerical integration process, or if it is about the implementation of a specific channel.

Best,

  Stephen

ᐧ

Project coordinator

http://openworm.org

"Talk is cheap, show me the code" - L. Torvalds

To view this discussion on the web visit https://groups.google.com/d/msgid/openworm-discuss/03e9e5ea-e5fc-4678-ae31-2b4f5777114e%40googlegroups.com.

[Robby Simpson]

Hi Stephen, All,

Thanks for the suggestion!

The issue only appears with the "e" gate of the calcium ion channel (ca_boyle at https://github.com/openworm/muscle_model/blob/master/NeuroML2/ca_boyle.channel.nml).  All others seem to perform as desired, which has me even more puzzled as the behavior/code is the same.

I have attached a very simple C file that only calculates the membrane potential based on this single GateHHtauInf (https://www.neuroml.org/NeuroML2CoreTypes/Channels.html#gateHHtauInf).  I've also attached the same in Java, in case some are more familiar with that language.

For reference, here is what happens with my original code with the other ion channels and gates:

• With no ion channels, the membrane potential stays a constant -75 mV
• With only the leak ion channel, the membrane potential goes to 10 mV (erev of that channel)
• With only k_fast, the membrane potential stays a constant -75 mV (which does strike me as odd as erev of k_fast is -54.9998 mV)
• With only k_slow, the membrane potential goes to -64.346100 mV (erev of that channel)
• With both leak and k_fast, the membrane potential goes to -8.299032 mV
• With both leak and k_slow, the membrane potential goes to -28.034009 mV
• With both k_slow and k_fast, the membrane potential goes to -64.346100 mV
• With leak, k_fast, and k_slow, the membrane potential goes to -28.035636 mV
• Any combinations including the "e" gate of ca_boyle quickly (~30 ms) goes to infinity.
• Any combinations of ca_boyle excluding the "e" gate (including the "custom" gate) do not have this issue and the membrane potential goes to 49.11 mV (erev of ca_boyle)
• With everything (leak, k_fast, k_slow, ca_boyle.f_gate, ca_boyle.h_gate) except ca_boyle.e_gate, the membrane potential goes to 0.725609 mV

Any help is, of course, greatly appreciated.

One theory I've had is that perhaps I'm over- or under-flowing in the calculations.  However, I've inserted quite a lot of debugging on my end and do not see this to be the case.

Any thoughts?

Thanks,

Robby

[Bóris Marin]

Hi Robby

You don't seem to be multiplying dq/dt (the "deltas" for gating variable dynamics) by the integration timestep (akin to to the 1e-3 factor in the v update).

It might be helpful to study https://en.wikipedia.org/wiki/Euler_method , which is the numerical integration method you are trying to use.

Best,

Boris

[Robby Simpson]

Hi Boris,

Thanks so much for your help!  Adding additional scaling certainly does stop the runaway to infinity.

However, I'm hoping you can help me to better understand the desired units (the reason the 1.0e-3 factor was in the v update).

For dq/dt, I have:
dq/dt = (inf - q) / tau
As inf and q are dimensionless, this leaves dq/dt with units of 1 / ms  (as tau is specified in ms).

For dv/dt, I have:
dv/dt = condDensity * fopen * (erev - v) / specificCapacitance
This gives units of:
(mS_per_cm2) * (mV) / (mF_per_cm2) = mS * mV / mF = mA / F = mV / s
As my voltages are all in mV, I think the numerator is correct.

So, should these units be per "s" or per "ms"?  And should the integration factor (1 ms or 1.0e-3) then be multiplied in?

I genuinely appreciate everyone's patience and hope I'm nearing the solution.  I will certainly update the code on github once correct in the hopes that others can learn from my trials.

Thanks again,
Robby

[Chris Jensen]

It's not obvious to me from glancing at the code the exact form of the equation being solved, but generally you want the time step to be smaller than the time constant (tau). For the e channel tau is roughly 0.1 and you are running stepsize = 1.

I can't help with the units.

[Bóris Marin]

Robby,

Once again, you need to familiarise yourself with numerical integration of ODEs. You are trying to use Euler's method, a fixed stepsize method, which involves choosing a discretization stepsize dt, and updating state as x_new = x_old + dx_dt * dt
That needs to be done for each one of the dynamical variables in your model.

Best,
Boris

[Bóris Marin]

Robby,

Once again, you need to familiarise yourself with numerical integration of ODE. You are trying to use Euler's method, a fixed stepsize method, which involves choosing a discretization stepsize dt, and updating state as x_new = x_old + dx_dt * dt

That needs to be done for each one of the dynamical variables in your model.

Best,
Boris

--

Visit us online at:
http://openworm.org
http://blog.openworm.org
http://github.com/openworm
http://twitter.com/openworm
https://plus.google.com/s/openworm
---
You received this message because you are subscribed to the Google Groups "OpenWorm-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openworm-discu...@googlegroups.com.
To post to this group, send email to openworm...@googlegroups.com.
Visit this group at https://groups.google.com/group/openworm-discuss.

To view this discussion on the web visit https://groups.google.com/d/msgid/openworm-discuss/bb2e15b1-e26f-4876-a0c3-bc7ad22d7a13%40googlegroups.com.

[Robby Simpson]

Eureka!!!  I think this (a stepsize of 1 ms when tau is roughly 0.1 ms) was the issue Chris!

I'm actually quite familiar with Euler's Method and so was scratching my head regarding the timestep issue - I believe all of my time derivatives were already using millisecond timesteps.  Though, admittedly, not very obviously.

However, moving to a smaller timestep has solved the problem (even just a timestep of 0.1 ms)!

Thanks again everyone for your help.  I will post updated code to github shortly that will hopefully help others in the future.

- Robby