Numerical methods in NEURON

[zyc]

Dear all
I am reading the NEURON book and try to understand the numerical methods in NEURON. In my understanding, when NEURON integrate one step, it first setup a matrix GV=I, where G is the Hines matrix representing the conductance, I is the current vector (including axial current, injected current and current of ion channels) and V is the vector of v_new (membrane voltages in next time step). I think V is the vector of v_new because when calculating current, the equation i=g*(v-e) is used, and v in this equation is membrane voltage.
However, in impelementation of NEURON, it seems that after we solve the matrix equation GV=I, we get a vector of dv (i.e. v_new-v_old) rather than v_new. Could you please tell me what's wrong with my understanding and why we get dv not v_new after solving GV=I in NEURON? Thank you so much.

[ted]

In my understanding, when NEURON integrate one step, it first setup a matrix GV=I, where G is the Hines matrix representing the conductance, I is the current vector (including axial current, injected current and current of ion channels) and V is the vector of v_new (membrane voltages in next time step).

Where did you find that? If you have the book, please read 3.3 Cables and note that spatial discretization of the cable equation produces a set of ODEs, each of which involves dv/dt at a particular location in space.

[zyc]

Thank you for replying. I watched the source code of NEURON, and I found that when NEURON computes the ion channels (e.g. hh.c), there are two statements:

Code: Select all

NODERHS(_nd) -= _rhs;
NODED(_nd) += _g;

so I think the matrix represent the conductance and rhs represent current. When current is computed by i=g*(v-e), and chapter 4 of NEURON book introduces how to compute v_new from v_old (numerical integration methods), so I guess the solution of the matrix equation are voltages rather than dv.

[ted]

OK, we're talking about different things.

The model is described by the set of ODEs
y' = f(y,t)
(bold terms are vectors).
The simplest method for calculating a(n approximate) numerical solution to this is the explicit Euler method
y(t+dt) = y(t) + dt * f(y,t)
where the first term on the right hand side is the value of y at the current time t, and the second term dt * f(y,t) is the "delta y vector" i.e. the amount by which y changes from time t to time t + dt.

NEURON never calculates a "delta y vector." Instead, by default NEURON uses the implicit Euler method, which requires solving
y(t+dt) = y(t) + dt * f(y(t+dt),t)
Gathering all terms at time t+dt onto one side of the equation, and all terms at time t onto the other, we have
y(t+dt) * (I - dt * f(y(t+dt),t)) = y(t)
(the I here is the identity matrix). Multiplying both sides by the inverse of (I - dt * f(y(t+dt),t)) gives us
y(t+dt) = y(t) * (I - dt * f(y(t+dt),t))^-1

[zyc]

Thank you, I'm more clear now. I see that in "nrn_fixed_step_thread" function, mainly three steps are done:
1) setup_tree_matrix
2) solve the tree matrix
3) update.
So step 1) is to setup the matrix (I - dt * f(y(t+dt),t)), step 2) is to get the inverse of matrix, is that true? If it's true, why in update step, the operation is vec_v(i) += vec_rhs(i), not vec_v = M*vec_v (M is the inverse of matrix)? That's also why I thought NEURON calculates a "delta y vector"

[hines]

In nrn/src/nrnoc/treeset.c there is a lengthy comment (beginning line 81) that discusses how NEURON is doing the calculation.

[zyc]

Thank you so much!