Defining integration method for conductance

[JRS92]

This is not a bug but more likely a lack of understanding on my part. I've tried to get a better grip on how channel integration works by looking at conductance, conductance2, and compartment but still not quite getting it.

I am defining a new conductance and want to integrate the equations below using Euler method , but is not working. I am able to get the expected result (inward current with inactivation) using the same equations using scipy odeint(), which uses a different solver. markovNA.txt

O = O + dt*(A_inf(V,Ca)*C + r2*(1-C-O) - (B_inf(V,Ca) + r1)*O); C = C + dt*(B_inf(V,Ca)*O + r3(V,Ca)*(1-C-O) - (A_inf(V,Ca) + r4)*C); g = gbar * fast_pow(O,3);

[sg-s]

Sure, let me see if I can make this clearer. Using default options, conductance objects are integrated as follows:

For each time step:

1. network.integrate calls the integrateChannels method in every compartment
2. the integrateChannels method then calls integrate on every conductance object in it
3. The conductance.integrate method does the actual (default) integration.

By default, the exponential Euler method is used. see this for the equation and how it's solved

If you want, you can define your own integration method (for anything) -- simply create a method called integrate() in that component and that will be used instead of the default since integerate is a virtual method -- if you define one in your object, it will be preferentially used.

[sg-s]

I can't see anything very obviously wrong in what you've done in your channel, let me know if it still doesn't work and I'll see if I can write it for you

where's your channel from? I don't think the link in your source file matches what you want

[JRS92]

Ah, I simply used another conductance file as a template which is why the url does not point to the right place.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4860362/

The model is shown in figure 1B and parameters are in the figure legend.

The integration works fine using odeint in python, which uses LSODA. I brought the same equations into matlab (changing syntax for c++) so I'm not sure why its not working. Below equations are how they appear there

As they appear in python:

      r1 = 1.
      r2 = .2
      r4 = .05
      ra =55.
      sa = 33.
      ka = -7.
      rb = 60.
      sb = 32.
      kb = 10.
      rc = 30.
      sc = 77.5
      kc = 12.
      gmax =50.

     ana =  ra/(1+np.exp((V+sa)/ka))
     bna =  rb/(1+np.exp((V+sb)/kb))
     r3 =     rc/(1+np.exp((V+sc)/kc))

    dOdt = (ana*C + r2*(1-C-O) - (bna + r1)*O)
    dCdt = (bna*O + r3*(1-C-O) - (ana + r4)*C)

    INa =gmax*O*O*O*(vi-ENa)

[JRS92]

Here are -100 to +40 mV voltage clamp steps from python showing the intended shape.

[JRS92]

OK, I checked and voltage-clamping this channel works fine. I was getting strange results when simulating Vm because the step size was too small. I changed sim_dt from 0.05 to 0.01 and now it seems to be working. I would prefer using the default integration method but I don't know how to use exponential euler for these equations. Do you have any recommendations?

[sg-s]

Vm because the step size was too small.

Do you mean the step size was too large? Euler integration will explode if the step size is too large -- which is why exponential Euler is preferred. Exponential Euler is only possible for equations of the form indicated in the docs -- so you're probably stuck with Euler (which is fine if the step size is sufficiently small)

Are you able to reproduce the original paper using a smaller step size?

[JRS92]

Ah yes, typo. Step size was too large!

I also didn't think it was possible which is why I went with Euler. Slower but not unbearable.

[sg-s]

excellent. feel free to send me a PR with your component if you want to!