multisplit compatibility with extracellular

[cccohen]

The following is a request to allow cells split by splitcell or, preferably, multisplit, to be run with the extracellular mechanism inserted. At present, this is not possible (see warning in documentation for multisplit and splitcell).

Thank you.

[nrnhines]

That is a very complex change and is unlikely to be accomplished in the forseeable future.
If you are using extracellular in order only to calcualate i_membrane, eg for calculation of the
local field potential, one can do this without inserting extracellular and instead use cvode.use_fast_imem(1)
whiich also allows use of multisplit.

[cccohen]

Unfortunately, it seems multisplit (or splitcell) are incompatible with cvode. This is specified in the same documentation (linked above).

[nrnhines]

My apologies. Cvode ( when the model is a set of ODE), multisplit, and parallel are allowed together and the reference needs to be updated. The presence of extracellular generates a DAE model which cannot
use either multisplit or parallel. The Cvode notation is confusing as the interpreter wrapper is called
Cvode and originally wrapped the sundials implementation of Cvode. This predated the ability
to simulate networks. Now the interpreter level wraps an implementiton level class called NetCvode
which in turn uses an implementation of Sundials Cvodes for ODE models and sundials IDA for
DAE models. ODE models can use parallel threads, mpi , and local and global variable time steps
managed by Cvodes.
DAE models are presenly only globally solvable on a single thread machine.

[cccohen]

Understood. Going forward, in terms of solutions excluding cvode, it seems DAE models could be implemented in multisplit using Sundials IDA. I suppose another possibility would be to convert the extracellular DAE to an ODE model, but that is likely non-trivial.

[cccohen]

Hi Michael, I would like to devote some renewed attention to this issue. I think it is critical for the future solvability of myelinated axons. If there are no better options, perhaps you can instruct me on what needs to be done? -Charles.

[nrnhines]

At a sufficiently conceptual level the tree solver is the same for fixed step , CVODE, and IDA without extracellular or LinearMechanism instances and without multisplit. Addition of multisplit into the fixed step and CVODE mix involves
a modification of the tree solver in two ways. All the details are in https://www.neuron.yale.edu/neuron/publications/fully-implicit-parallel-simulation-single-neurons but the primary constraint is that pieces of a tree must have two or fewer distinct connection points to other pieces. Pieces with only one connection point are handled very similarly to trees. I.e the connection point is the root, triangularization proceeds normally from the leaves, but backsubtitution needs to wait until the node that the connection point connects to is fully solved. Pieces with two connection points can be triangularized from the leaves to nodes that define the "backbone" between the two connection points. This backbone tridiagonal can be transformed into an "N" matrix but further solving needs to wait until the two nodes that the two connection points connect to are fully solved. The nodes that are in the set of nodes that need to be "fully solved" forms
a reduced tree topology that can be solved by a normal tree solver. At that point all the pieces can be back substituted
to finish the whole tree solve.

Now, the addition of extracellular does not much modify this conceptually. Ie each node needs to be considered not as a single state but as a short vector of size equal to the number 1 (internal) + number of extracellular layers. (default extracellular with two layers means a size of 3). The processing is that each multiplication, addition, subtraction, and division is not for scalar matrix elements but 3x3 matrix matrix elements. (at least for connection points and perhaps for the entire backbones and reduced tree). Not clear what the precise performance impact would be.

Now I have to throw another software maintenance technicality into the mix. Ie I'm stuck at the end of a long languishing
upgrade to using SUNDIALS with NEURON (the cvode3 branch). So long languishing that although it is called cvode3
it is now for SUNDIALS version 5. The problem Is that the NVECTOR module for SUNDIALS has changed its API for threads and mpi. And I have not figured out how to use it properly. Another, and more vexing, issue is that cvode3 has a subsidiary branch called idainit which attempts to move from my heuristic IDA initialization methods (small initial step) to the IDA supported initialization methods. This latter requires that one mark each state as either differential or algebraic. I can do that by introducing new states with the algebraic equation vmij = vi - vj, marking the left hand side states as
algebraic, and the right hand side as differential, and replacing all c*(dvi/dt - dvi/dt) terms with c*dvmij/dt. I've gotten it mostly to work but the extra states and equations are highly unsatisfying and needed only for initialization.

[cccohen]

I'm happy to read some progress has been made. It even sounds like it may be possible under the current architecture. I guess, for one, I am not sure how I can help at this point in the process. For another, I am not quite sure I understood what is preventing completion. You mention, "the extra states and equations are highly unsatisfying". Maybe this is something I can help with?

[nrnhines]

extra states and equations are highly unsatisfying

I should mention that for multisplit solving they also add n-1 to the rank of the matrix matrix elements so those elements are then 5x5 matrices instead of just 3x3 (for two extracellular layers)

You are welcome to jump into the swamp with me:) The extra states really did resolve the initialization issues for LinearMechanism (extra state for two nodes connected by a capacitor or inductor) and I don't think there is any substantive problem for the extracellular case (always an extra state due to membrane capacitance and occasionally an extra state if extracellular capacitors between extracellular layers are non-zero). I expect the problem is a bug that I can't see but will be obvious when I do see it. These long-languishing branches are a pain since when one starts up with them again it is necessary to merge a dramatically changed master into them. Also need to do the transformation from autotools to cmake so the sundials part can be a submodule (actually, I can't remember if I did that already). Anyway, the best way to start is to have conversations where I give you a tour of the code (and relearn what is there myself). Then you will be better able to assess whether the effort is reasonable.