1D Electrodiffusion Cable Model Compatible with NEURON

Particularly useful chunks of hoc and/or NMODL code. May be pedestrian or stunningly brilliant, may make you gasp or bring tears to your eyes, but always makes you think "I wish I had written that; I'm sure going to steal it."
Post Reply
Keivan
Posts: 127
Joined: Sat Apr 22, 2006 4:28 am

1D Electrodiffusion Cable Model Compatible with NEURON

Post by Keivan » Sat Oct 14, 2017 1:35 pm

A few years ago, I and some of my colleagues started to work on a one dimensional electrodiffusion cable model based on Mori 2007 [A Three-Dimensional Model of Cellular Electrical Activity]. Our mathematical endeavours were successful, but we never tested our equations, since I moved to US to get my PhD in computational neuroscience and it was very hard to find spare time. I want to put our calculations here in this forum in case anybody is interested in testing them, or continue this project.

https://1drv.ms/w/s!Ato0Wsph7iDRjr4t-FYt5WfRkFeeiw

We had also partly implemented a C++ code to test the equations, but it was never completely implemented and properly debugged.

https://1drv.ms/f/s!Ato0Wsph7iDRjr08UXUdh4I3l9JLpA

We had this idea of using predictor corrector numerical methods to increase the accuracy of calculations. I have written the following code for the predictor section of the code.

Code: Select all

#include <iostream>
/*
This is a predictor corrector numerical solver:
here we suppose that dx/dt = f, so f is the slope
Knowing the slope of a function at several points,
it would be possible to predicts its value at a next point
The value of slope (f), and the value of x are stored in circular arrays
this array helps us to access these parameters
at this iteration                 by x[ 0] or f[ 0]
the previous iteration            by x[-1] or f[-1]
the iteration before previous one by x[-2] or f[-2]
and so on
knowing these values it would be possible to predict the value of next x
with different accuracies: 2nd, 3rd, or 4th order accuracies are implemented
in this class.

note initialization of predictor is important.
predictor-corrector solvers are multistep methods.
To achieve a 4th order accuracy, only the initial values of x should be
updated by at least a single step 3rd order accurate (or more) numerical method.
note: Adams-Moulton method is more stable than Adams-Bashforth's
*/
template <typename T>
class predictor
{
public:
    // initialized at first the values of f and x to steady state
    predictor(T initial_value)
    {
        x.initialize(initial_value);
        f.initialize(0);
    }
    ~predictor() {}
    void update(T current_value, float deltaT)
    {
        f.update((current_value - x[0])/deltaT);
        x.update(current_value);
    }
    T predict2ndOrder(float deltaT)
    {
        return x[0] + deltaT/2*(3*f[0] - f[-1]);
    }
    T predict3rdOrder(float deltaT)
    {
        return x[0] + deltaT/12*(23*f[0] - 16*f[-1] + 5*f[-2]);
    }
    // 4th order predictor Adams-Bashforth
    T predict4thOrder(float deltaT)
    {
        return x[0] + deltaT/24*(55*f[0] - 59*f[-1] + 37*f[-2] - 9*f[-3]);
    }
    // 4th order predictor-corrector Adams-Moulton
    T pc4thOrder(float deltaT)
    {
        return x[0] + deltaT/24*(9*(predict4thOrder(deltaT)-x[0])/deltaT + 19*f[0] - 5*f[-1] + f[-2]);
    }
    inline T operator[](int index)
    {
        return x[index];
    }
private:
    // size should be in power of 2 = (2,4,8,16,..,128)
    // since we are going to use high performance bit-wise modulus
    #define size 4
    #define SizeNegOne 3
    template <typename S>
    class carray
    {
    public:
        carray()
        {
            head = SizeNegOne;
        }
        // copy constructor
        carray(const carray& that)
        {
            head = that.head;
            for (int i=0; i<size; i++)
                array[i] = that.array[i];
        }
        // copy assignment operator
        carray& operator=(const carray& that)
        {
            head = that.head;
            for (int i=0; i<size; i++)
                array[i] = that.array[i];
            return *this;
        }
        // fill all the elements of the array with
        // a specific value
        void initialize(S n)
        {
            for (head=0; head<size; head++)
                array[head]=n;
            head = SizeNegOne;
        }
        // add a new value to head of this array
        // and delete one element from its tail
        // we use array for this purpose so,
        // In the other words,
        // if you have reached the end of the array,
        // start again from its head and therefore
        void update(S value)
        {
            // head index = (head+1)%size
            head = (head+1) & SizeNegOne;
            array[head]=value;
        }
        // access elements of this array with negative indexes
        // which shoes how many iterations before they have updated
        inline S operator[](int index)
        {
            return array[(head + index) & SizeNegOne];
        }
    private:
        S array[size];
        signed char head;
    #undef size
    #undef SizeNegOne
    };
    carray<T> f, x; // f=dx/dt
};

int main()
{
    // create our objects and initialize them with a fixed value
    predictor<float> phi(-64), ph(0);
    // suppose that each deltaT = 1 ms, we update the phi value
    phi.update(10, 1);
    phi.update(11, 1);
    phi.update(12, 1);
    phi.update(13, 1);
    phi.update(14, 1);
    // after initialization, we ask the predictor
    // its predictions for the next iteration that would take
    // deltaT = 1 ms
    std::cout<<phi.predict2ndOrder(1)<<std::endl;
    std::cout<<phi.predict3rdOrder(1)<<std::endl;
    std::cout<<phi.predict4thOrder(1)<<std::endl;
    std::cout<<phi.pc4thOrder(1)     <<std::endl;

    // here I show that it is possible to use this object like
    // a regular carray to access previous values of phi
    // we need this feature for the calculation of phi and Concentration
    std::cout<<phi[0] <<std::endl;
    std::cout<<phi[-1]<<std::endl;
    std::cout<<phi[-2]<<std::endl;
    std::cout<<phi[-3]<<std::endl;

    // to make sure that the rule of three is satisfied
    ph=phi;

    // here we test the copied object
    std::cout<<ph.predict2ndOrder(1)<<std::endl;
    std::cout<<ph.predict3rdOrder(1)<<std::endl;
    std::cout<<ph.predict4thOrder(1)<<std::endl;

    std::cout<<ph[0] <<std::endl;
    std::cout<<ph[-1]<<std::endl;
    std::cout<<ph[-2]<<std::endl;
    std::cout<<ph[-3]<<std::endl;

    return 0;
}
I will be in SfN Neuroscience 2017 in DC this year to present my PhD project in case anybody was interested.

Post Reply