sg-s / xolotl

A MATLAB neuron simulator. Very fast (written in C++). Flexible (fully object oriented). Immediate (live manipulation in MATLAB). Comes with a powerful parameter optimizer. Get started ➡️
https://go.brandeis.edu/xolotl
GNU General Public License v3.0
45 stars 8 forks source link

Documentation should talk about how to plot activation, inactivation variables over time (ChannelProbe) #544

Closed JRS92 closed 3 years ago

JRS92 commented 3 years ago

These are not included in any of the output_type(s).

For example, I would like to see how m and h change over time. Or any other state variables for a conductance.

sg-s commented 3 years ago

unfortunately, there is no easy way to plot them. a decision was made a long time ago that this use case was marginal, and therefore we dropped it (which greatly sped up the simulation and used less memory)

I'm guessing you want to plot the activation & inactivation variables of some custom channel? In principle, one can write a "probe" mechanism that reads it out and returns it -- now that I think of it, I think I will do exactly this. Give me a minute...

JRS92 commented 3 years ago

Yes, exactly! Thank you!

sg-s commented 3 years ago

OK I wrote a quick and dirty probe. It's called ChannelProbe

% example usage

% example neuron
x = xolotl.examples.neurons.BurstingNeuron;

% let's probe the NaV channel
x.AB.NaV.add('ChannelProbe');

% read out the m, h from NaV
% remember the channel probe is a mechanism, so the 3rd output is what we want
[~,~,M]=x.integrate;

plot(M)

and we see something like this

Screen Shot 2020-12-10 at 5 22 00 PM

hope this helps

JRS92 commented 3 years ago

Ah I see. I am trying to do this for a C - O - I markov model and thus I don't think this would work, since it uses a custom integration routine. Maybe I should write the channel as a mechanism instead of a conductance?

sg-s commented 3 years ago

No, don't write a channel as a mechanism. That will break a lot of things.

The way ChannelProbe is simply to read out m,h at every time step. If your custom channel has some custom variables, insert whatever you want to read out in this line and you should be good to go.

JRS92 commented 3 years ago

The conductance has the custom variables C and O but they are not being read.

In file included from /Users/j/Documents/MATLAB/xolotl/X_a539ec013670e1a1101cd0ae8a0ee9d9.cpp:26:
/Users/j/Library/Application Support/MathWorks/MATLAB Add-Ons/Toolboxes/xolotl/c++/mechanisms/ChannelProbe.hpp:66:32: error: no member named 'O' in
'conductance'
    cont_state[idx] = channel->O;
                      ~~~~~~~  ^
/Users/j/Library/Application Support/MathWorks/MATLAB Add-Ons/Toolboxes/xolotl/c++/mechanisms/ChannelProbe.hpp:68:32: error: no member named 'C' in
'conductance'
    cont_state[idx] = channel->C;
                      ~~~~~~~  ^
2 errors generated.

Error in xolotl/compile (line 59)
    mex('-silent',ipath,mexBridge_name,'-outdir',filelib.cachePath('xolotl'))

Error in xolotl/integrate (line 107)
            self.compile;

Error in untitled4 (line 12)
[V,~,M]=x.integrate;
sg-s commented 3 years ago

i can only echo what the compiler says: it looks like there is no member called "O" or "C" in your channel. it should be public. can you send me the channel you're trying to probe?

JRS92 commented 3 years ago

The channel is there along with the ChannelProbe mechanism.

 markovNA object (3e2f4cd) 

            kd : -2.19298807971014
             O : 0
            kb : 7.7542293131695
            rc : 0.00324397221993918
            kc : 8.17784833805925
             E : 40
            sc : 0
            sd : 14.3831303102262
          gbar : 22.6086956521739
            sa : 48.2375264040895
            rd : 2.81426630434783
            r4 : 0
            sb : 44.4485864202081
            r2 : 0.0180255968638314
            ra : 44.0520551839643
            ka : -5.34381014886591
            rb : 54.827282669052
             C : 1
  ChannelProbe : ChannelProbe object
sg-s commented 3 years ago

ah, i see the confusion.

what cpplab shows is the constructor signature: every field here is an argument to the constructor. however, you can have things in the constructor that aren't members of that channel.

For example, in this mechanism, you see how "Target" is defined as a member (under public)? Only these fields will be readable by anything else.

To summarize, you should declare "O" and "C" in your channel under public. I bet you're not doing it

this is wrong, see below

JRS92 commented 3 years ago

Here is the conductance. It is from this paper. I am optimizing the parameters to recordings from neurons I am studying. The default values are old/work in progress. I will update them and submit a pull request with my folder full of conductances once the paper is published :)

#ifndef MARKOVNA
#define MARKOVNA
#include "conductance.hpp"

//inherit conductance class spec
class markovNA: public conductance {

public:

    double r2 = 0.2;
    double r4 = 0.05;
    double ra = 55;
    double sa = 6.4;
    double ka = -15.9;
    double rb = 60;
    double sb = 32;
    double kb = 10;
    double rc = 30;
    double sc = 77.5 ;
    double kc = 12; 
    double rd = 60;
    double sd = 32;
    double kd = 10; 
    double C = 1.0;
    double O = 0.96; 

    // specify parameters + initial conditions
    markovNA(double gbar_,double E_, double O_, double C_, double r2_, double r4_, double ra_, double sa_, double ka_, double rb_, double sb_, double kb_, double rc_, double sc_, double kc_, double rd_, double sd_, double kd_)
    {
        gbar = gbar_;
        O = O_;
        C = C_; 
        E = E_;

        // activation parameters
        r2 = r2_;
        r4 = r4_;
        ra = ra_;
        sa = sa_;
        ka = ka_;
        rb = rb_;
        sb = sb_;
        kb = kb_;
        rc = rc_;
        sc = sc_;
        kc = kc_;
        rd = rd_;
        sd = sd_;
        kd = kd_; 

        // defaults 
        if (isnan(gbar)) { gbar = 0; }

        if (isnan (E)) { E = 40; }

        UseMInfApproximation = 0;
        UseHInfApproximation = 0; 

        // do not allow this channel to be approximated

        E = 40; 
    }

    double A_inf(double, double);
    double B_inf(double, double);
    double r3(double, double);
    double r1(double, double);
    void integrate(double,double);
    string getClass(void);   
};

string markovNA::getClass(){return "markovNA";}

double markovNA::A_inf(double V, double Ca) {return ra/(1+exp((V+sa)/ka));}
double markovNA::B_inf(double V, double Ca) {return rb/(1+exp((V+sb)/kb));}
double markovNA::r3(double V, double Ca) {return rc/(1+exp((V+sc)/kc));}
double markovNA::r1(double V, double Ca) {return rd/(1+exp((V+sd)/kd));}

void markovNA::integrate(double V, double Ca) {   

   O = O + dt*(A_inf(V,Ca)*C + r2*(1-C-O) - (B_inf(V,Ca) + r1(V,Ca))*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); 

}

#endif
sg-s commented 3 years ago

Sorry, I now understand your issue.

The reason this doesn't work is because while your channel defines "O" and "C", it is actually stored as a reference to a type "conductance" (which it inherits from). The problem is that conductance doesn't know what "O" and "C" are

What you want to do is read out a member of a class externally, but in principle without knowning about that class. This is impossible in C++ (I think)

This is why the getState method exists. You call getState with some integer argument, and it returns whatever you want it to return. It's a workaround, but I've used it extensively in xolotl. I've modified your channel to include this method so you can see how it works. Add this to your channel:


double markovNA::getState(int idx) {
    if (idx == 1) {return O;}  //ideally this should be a switch
    else if (idx == 2) {return C;}
    else {return std::numeric_limits<double>::quiet_NaN();}
}

and modify ChannelProbe so that it looks like:

int ChannelProbe::getFullState(double *cont_state, int idx) {
    // read out O and C
    cont_state[idx] = channel->getState(1);
    idx++;
    cont_state[idx] = channel->getState(2);
    idx++;
    return idx;
}

I believe this solves your issue

JRS92 commented 3 years ago

I added getState to the conductance but now I am just getting

error: no member named 'getState' in 'conductance'
    cont_state[idx] = channel->getState(1);
JRS92 commented 3 years ago

I solved this by adding getState as a virtual function in conductance along with the changes you suggested to ChannelProbe and the markovNA conductance. Thank you! Will leave open since it looks like you want to update the documentation.

sg-s commented 3 years ago

I solved this by adding getState as a virtual function in conductance

Holy shit you're right, there is no virtual getState in conductance, and there should be. I'll add this. Thanks!