How the diffusion current is calculated?

Support
,
1

Description of problem

I’m wondering how the diffusion current is calculated in the multi-compartmental neurons. There’s not much in the documentation about it. I ask this because I assumed that diffusion current shouldn’t impact the voltage trace for SpatialNeurons with n=1 segment. Similarly, the value of Ri argument must be irrelevant. But the simple code below shows otherwise.

Minimal code to reproduce problem

The following script defines and simulates a neuron governed by a leakage current and an external current I. In addition, two auxiliary functions are defined:

  • make_n_run(scale = 1, shape='sphere', L=1, Ri_scale=1): Makes a morphology of type Soma if shape= 'sphere' otherwise, a single compartmental Cylinder. By default, the area of both objects is 1 * meter**2. Yet, it is possible to scale the diameter of both objects by a factor of scale. In addition, Ri_scale can scale the axial resistance. Yet, I guess that it should not be relevant for the single-compartmental object created here.
  • plotter(neuron, mon): Plots, v, I, Im during the course of simulation.
import brian2 as b2
from brian2 import Soma, Cylinder, Section
from brian2.units import *
import matplotlib.pyplot as plt
from brian2tools import *
import numpy as np
from brian2.units.constants import faraday_constant as F
from brian2.units.constants import gas_constant as R


El = -84*mV
gl = 30* nsiemens/(1*meter**2) 
Ri = 0.4*metre/siemens # axial resistance
Cm = 1.4*pF/meter**2

test_eq = """
Im = gl * (El - v) + I: amp/meter**2
I : amp/meter**2
"""

def make_n_run(scale = 1, shape='sphere', L=1, Ri_scale=1):
    """
    Scales diamter by a factor of `scale`. Also one can select the shape.
    """
    b2.start_scope()

    if shape=='sphere':
        morpho = Soma(diameter=scale*np.sqrt(1/np.pi)*meter) # area scaled by sqrt(scale)
    else:
        morpho = Cylinder(diameter=scale/(np.pi*L)*meter, length=L*meter, n=1) # area scaled by scale
    
    neuron = b2.SpatialNeuron(morphology=morpho, 
                          model= test_eq,
                          Cm=Cm,         
                          Ri=Ri*Ri_scale, # irrelevent for single compartments or not?
                          method='exponential_euler',
                          )
    neuron.v = El
    neuron.I = 0*mA/meter**2
    mon = b2.StateMonitor(neuron, ['v', 'I', 'Im'], record=True)
    
    b2.run(2*ms)
    neuron.I = 0.3*nA/meter**2
    b2.run(2*ms)
    neuron.I = 0*amp/meter**2
    b2.run(3*ms)
    
    return neuron, mon


def plotter(neuron, mon):
    fig, axs = plt.subplots(3,1, figsize=(8,8))

    for i in range(0,mon.v.shape[0]):
        axs[0].plot(mon.t,mon.Im[i], 'k', 
                    alpha=(i+1.)/mon.v.shape[0],label=i, )
        axs[1].plot(mon.t,mon.I[i], 'r', alpha=(i+1.)/mon.v.shape[0],label=i, )
        axs[2].plot(mon.t, mon.v[i], 'b', alpha=(i+1.)/mon.v.shape[0],label=i, )

    axs[0].set_title('Im [Amp/m^2]')
    axs[1].set_title('I [Amp/m^2]')
    axs[2].set_title('v [volt]')

    # originally was desigend for a REAL multicompartmental neuron
    for ax in axs:
        ax.legend(ncol=6)
    
    print(neuron.Im[0], neuron.Im[-1]) # for reference

Using these functions we can look at the following cases:

  1. large soma (scale = 1e6)
  2. small soma (scale = 1e-4)
  3. thick axon (scale = 1e7, shape='cylinder')
  4. thin axon (scale = 1e-7, shape='cylinder')
  5. large soma with high axial resistivity (scale = 1e6, Ri_scale=1e8)
  6. small soma with high axial resistivity (scale = 1e-4, Ri_scale=1e8)
  7. thick axon with high axial resistivity(scale = 1e7, shape='cylinder', Ri_scale=1e4)
  8. thin axon with high axial resistivity(scale = 1e-7, shape='cylinder', Ri_scale=1e4)

with the following commands:

neuron, mon= make_n_run(*suitable_args)
plotter(neuron, mon)

What you have aready tried

Turned out that once the axial resistance is not high (case 1-4), the results are sensitive to the object’s diameter, yet in different manners. While the voltage trace of larger somas (or tiny axons) are independent of the diameter, small somas (or huge axons) exhibit some diameter-dependent response:

small soma

soma_small
soma_small464×509 44 KB

large axon

axon_large
axon_large475×511 46.6 KB

whereas for large soma, small axon, or high enough axial resistivity, all configurations converge to the following output:

large axon with high axial resistivity

axon_large_high_R
axon_large_high_R523×506 47.5 KB

Note that the figure above coincides with a point neuron implementation – where there’s no propagation.

So my question is:

Why do the voltage traces (blue lines on the lower panels) depend on the diameter? Given that the equations are expressed in terms of current density, I expect all morphologies, independent of size or shape, to exhibit the behavior of a point neuron.

Thanks for sharing your thoughts. A reference link would also suffice.

Arash

2

Hi @arashgmn. Unfortunately, the multicompartmental simulations are lacking a bit of documentation for the numerical methods that are used. As a reference for the basic approach, please have a look at:

Mascagni, Michael V., and Arthur S. Sherman. “Numerical Methods for Neuronal Modeling.” Methods in
Neuronal Modeling 2 (1998).
https://www.cs.fsu.edu/~mascagni/papers/RCEV1996_1.pdf

I have run some tests with similar code as yours, and I don’t think there is anything wrong in particular. Varying the “irrelevant” parameters within reasonable bounds does not change anything. I think what you are seeing are simply numerical issues when you use extreme values (which you can somewhat compensate for by using extreme values elsewhere). The diffusion current simulation does quite a few things, including solving systems of equations for consistency between the currents flowing in and out of compartments. Almost all of this is of course unnecessary for a single compartment, but this wasn’t a relevant use case that we handled, i.e. even for a single compartment all these computations are performed. This will introduce minor numerical artifacts, that will become relevant when multiplied by huge areas (fun fact: in your “large axon” example, the axon has a membrane surface approximately equal to the surface area of Paris :blush: ), or divided by tiny values. Basically, I think the equations that are solved contain some v = \dots + a - b term at some point, where a=b for a single compartment – if we had a special case for a single compartment, these terms wouldn’t appear, but our equations take only the general case into account. With extremely big/small values for areas/resistances, small numerical discrepancies in a and b will lead to a significant difference.

But all this is not based on a thorough analysis of the calculations for the single compartment case – if you find something that is different from a point-neuron model for reasonable parameter values please let us know!

3

Thanks again, @mstimberg, for your detailed explanation! After extra modifications, I came to the same conclusion–numerical issues rather than something fundamental. But honestly, wasn’t sure. Your response was a relief. And great point about the area of my gigantic neuron :sweat_smile:!

Btw, the link for the book mentioned isn’t working for me (even after appending the pdf to). Maybe you can recheck it for future reference.

Thanks again for your support!
Arash

4

Sorry about that, try https://www.cs.fsu.edu/~mascagni/papers/RCEV1996_1.pdf (I also edited my comment above for anyone who stumbles across this post in the future).

1 Like