Page 1 of 1

### 1D Electrodiffusion Cable Model Compatible with NEURON

Posted: 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)/deltaT);
x.update(current_value);
}
T predict2ndOrder(float deltaT)
{
return x + deltaT/2*(3*f - f[-1]);
}
T predict3rdOrder(float deltaT)
{
return x + deltaT/12*(23*f - 16*f[-1] + 5*f[-2]);
}
// 4th order predictor Adams-Bashforth
T predict4thOrder(float deltaT)
{
return x + deltaT/24*(55*f - 59*f[-1] + 37*f[-2] - 9*f[-3]);
}
// 4th order predictor-corrector Adams-Moulton
T pc4thOrder(float deltaT)
{
return x + deltaT/24*(9*(predict4thOrder(deltaT)-x)/deltaT + 19*f - 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 <<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 <<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.