Surface area of compartmental models

Support
1

Hi,
I have a question about how the surface area of a compartmental model is computed.
My understanding from the docstring of a Cylinder is that the Length argument specifies its total length, independent of how many compartments (n) it has.
Yet the input resistance of a cylinder appears to change when the number of compartments changes in a way that cannot be explained by the increased accuracy of the spatial discretization.
I suppose I’m getting the concept of the “compartments” wrong - I’m thinking of it as the equivalent of “nseg” in the NEURON simulation environment.

Here is an example to reproduce the issue:

import numpy as np
import brian2 as b2  # Sorry I'm obsessed with separating namespaces

gL = 1e-4*b2.siemens/b2.cm**2  # Leak conductance
EL = -80*b2.mV
tbase = 50.0*b2.ms
tpulse = 200.0*b2.ms 
duration = 400.0*b2.ms
Iamp = 10e-12*b2.amp  # Current injection amplitude (10 pA)
nseg = 11 # Number of compartments in the dendrite; change this to see the effects on input resistance 

# A simple passive membrane
eqs = '''
Im = gL*(EL-v) : amp/meter**2
I : amp (point current)
'''

morpho = b2.Soma(diameter=20*b2.um)
morpho.dendrite = b2.Cylinder(length=1000*b2.um, diameter=10*b2.um, n=nseg)

# Axial resistance is negligibly low - effectively this is a single compartment
neuron = b2.SpatialNeuron(morphology=morpho, model=eqs, Cm=1*b2.uF/b2.cm**2, Ri=1e-3*b2.ohm*b2.cm)

neuron.v = EL
neuron.I = 0*b2.amp

# Monitors
mon = b2.StateMonitor(neuron, 'v', record=True)
b2.run(tbase)
neuron.I = Iamp
b2.run(tpulse)
neuron.I = 0*b2.nA
b2.run(duration-tpulse-tbase)

# Compute input resistance from simulation values
vbase = mon.v[0][mon.t<tbase][-1]
vmax = mon.v[0][mon.t<tbase+tpulse][-1]
delta_v = vmax-vbase
rin_sim = delta_v / Iamp

# Compute total membrane area according to brian2, and compare with theory
total_area_sim = np.sum(morpho.area)  + np.sum(morpho.dendrite.area)
total_area_theory = \
    4.0*np.pi*(morpho.diameter/2.0)**2 + \
    np.pi*morpho.dendrite.diameter[0]*np.sum(morpho.dendrite.length)

# Compute theoretical input resistance
rin_theory = 1.0/(gL * total_area_sim)

print("Total surface area (brian2): ", total_area_sim)
print("Total surface area (theory): ", total_area_theory)
print("Input resistance (brian2): ", rin_sim)
print("Input resistance (theory): ", rin_theory)

Output for nseg=11:
Total surface area (brian2): 0.03267256 mm^2
Total surface area (theory): [0.03267256] mm^2
Input resistance (brian2): 0.36725744 Gohm
Input resistance (theory): 30.60671983 Mohm

I.e. the input resistance is nseg+2 times higher than what I would have expected.

I’m running brian2 v2.4.2 with Python 3.8.

What am I getting wrong here?

Thank you for your help.

Christoph

1 Like
2

Hi Christoph. Your understanding is correct, the number of compartments is essentially identical to nseg in NEURON, except that in Brian it is a fixed property of the morphology, i.e. you cannot change the number of compartments after its creation without creating a new one (in NEURON you can set the nseg attribute on any section at any time). And changing it should not change the overall area or resistance.

The problem in your example is that you are injecting a point current via neuron.I = Iamp, which means: “inject a point current into each compartment”. Therefore with each additional compartment you also get an additional “electrode” injecting a current. Instead, you can inject it into a single place via something like:

neuron.dendrite[500*b2.um:500.1*b2.um].I = Iamp

Which would inject the current into whatever compartment happens to be in the middle of the dendrite. As a side remark: we have some changes in the pipeline that would allow you to simple write [500*b2.um] instead of this slightly weird range, but at the moment this is not supported. However, using a concrete integer index will already work:

neuron.dendrite[len(neuron.dendrite)//2].I = Iamp

Hope that makes things clear!

3

Thank you. Yes that clears things up. I have re-run the code and it works as expected.
When the morphology is mapped to the SpatialNeuron, is neuron[0] always guaranteed to be the soma or root compartment? Or is there a less ambiguous way of referring to the root section?

4

If you start your morphology declaration with morpho = b2.Soma(...), then neuron[0] will always correspond to the soma. The mapping from the morphology tree to the flat representation is deterministic, you can get the indices of individual sections by accessing the indices attribute, e.g.

print(morpho.dendrite.indices[:])

will give you the indices of the dendrite compartment. As for all morphological attributes, you can also get them directly from the SpatialNeuron (neuron.dendrite.indices[:]). There’s a subtle difference between the two accesses, though:
Accessing through the morphology will always give you the attribute for the respective section, e.g. morpho.area will give you the area of the soma. On the other hand, accessing through the SpatialNeuron will always work on the whole subtree (i.e. neuron.area will give you the area of all the compartments). We should probably make this clearer in the documentation, I now realize that this difference is probably why you expected that neuron.I = ... would set the current in the soma.

5

I now realize that this difference is probably why you expected that neuron.I = ... would set the current in the soma.

Yes, that’s why I tripped.