Hi, I have a question regarding the spatial neuron.

I cannot understand how you derived the equation related to the coupling system.

What does u_plus and u_minus indicate and what does P matrix refer to?

Thanks.

Hi, I have a question regarding the spatial neuron.

I cannot understand how you derived the equation related to the coupling system.

What does u_plus and u_minus indicate and what does P matrix refer to?

Thanks.

I get it nowâ€¦ u_plus and u_minus refer to the current flowing through the left and right side of each section.

Also, you can derive the equation for the matrix P assuming that the total current flowing into each node is zero.

Hi @hunjunlee, we should definitely document this in detail somewhere, it has been on my list for quite a whileâ€¦

Our approach basically follows this one (in particular the â€śMethods for PDEsâ€ť part):

Mascagni, M. V., & Sherman, A. S. (1998). Numerical methods for neuronal modeling. Methods in Neuronal Modeling, 2.

https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=bc950e1d71da936934d9f5ff650d493bfd82e9fb

I seeâ€¦ Can I ask another question?

The spatial neuron update is costly due to the implicit Euler method which requires evaluating the inverse of the tridiagonal matrix.

Also, the explicit Euler method does not require such complex equations and potentially reduce the compute overhead.

Soâ€¦ can we adopt the semi-implicit Euler method instead?

For example,

Existing Brian2 code requires solving the following equation (for the compartment i).

Cm * (V i,t+1 - V i,t) / dt + Iion i,t+1 = invr * (V i-1,t+1 - V i,t+1) + invr * (V i+1,t+1 - V i,t+1)

(I made some simplifications)

On the other hand, if you replace the equation as follows,

Cm * (V i,t+1 - V i,t) / dt + Iion i,t = invr * (V i-1,t+1 - V i,t+1) + invr * (V i+1,t+1 - V i,t+1)

(I replaced Iion i,t+1 to Iion i,t)

In this case, the tridiagonal matrix would be constant at all time steps and the simulator no longer requires a complex equation to evaluate the inverse of the matrix in real time.

Can you give me some comments on this type of approximation?

Thanks.

Hi again. I am a bit confused here. We already do a staggered update, i.e. we first calculate the current from the ionic channels, and then the current flow between compartments. I am not sure how the matrix could be â€śconstant at all time stepsâ€ť, given that the membrane potential (and the resulting current flow) changes dynamically. Or maybe I misunderstood your suggestion?

Thanks for the reply.

Maybe I did not make myself clear.

In the current implementation, Brian2 code solves the following equation:

Cm * (V i,t+1 - V i,t) / dt + g_tot * Vi,t+1 - I0 = invr * (V i-1,t+1 - V i,t+1) + invr * (V i+1,t+1 - V i,t+1).

For that, it coverts the equation into the following form:

invr * Vi-1,t+1 + (-invr - invr - g_tot - Cm/dt) * Vi,t+1 + invr * Vi+1,t+1 = Cm / dt* Vi,t - I0.

Then, it uses the matrix to solve the problem:

G * V_vector,t+1 = Cm / dt* V_vector,t - I0.

where G is a tridiagonal matrix with the supra- and sub-diagonal elements being invr and the diagonal elements being (-invr - invr - g_tot - Cm/dt).

However, the problem is that the diagonal elements are not constant as g_tot changes, and thus, Brian2 should solve a complex equation for each time step.

So, I suggest adopting a mixed Euler method and solve the following equation instead:

Cm * (V i,t+1 - V i,t) / dt + g_tot * Vi,t - I0 = invr * (V i-1,t+1 - V i,t+1) + invr * (V i+1,t+1 - V i,t+1).

(I used g_tot * Vi,t instead of g_tot * Vi,t+1)

In this case, the simulator should solve the following equation.

G * V_vector,t+1 = Cm / dt* V_vector,t + g_tot * V_vector,t - I0.

where G is a tridiagonal matrix with the supra- and sub-diagonal elements being invr and the diagonal elements being (-invr - invr - Cm/dt).

Therefore, the G matrix becomes constant for each time step and most of the computations can be done in the â€śbefore_runâ€ť function.

So, Iâ€™m wondering if this kind of implementation is also a plausible option.

Hi again @mstimberg.

I implemented the above code in cython and observe that it achieves around 10% speedup in a typical random balanced network with 8-compartment / 7-section morphology with AdEX-like neuron model.

Also, I assume that the such implementation would greatly simplify the hardware implementation.

So, I wonder what you guys think.

Hi @hunjunlee, this is just to let you know that I did not forget about this issue! I just came back from holidays and will need a few days before I can have an in-depth look into this. Thanks for your patience