# Description of problem

I want to model a multicompartmental neuron, stimulate some part of it, and study spike propagation. To compare my implementation with the literature, first, I implemented Equations for a point neuron and extended it for a SpatialNeuron using the point current flag. It works fine. But dor other reasons, I want to describe my equation per surface area. Unfortunately, I fail to figure out how.

My approach was as follows:

1. make a single compartmental neuron (a spherical Soma) with an area of 1*meter**2. Inject a fixed total current (here .3*nA) using point current; and
2. make the same morphology, but this time, inject a current density of 0.3*nA/meter**2.

Since the area is 1, both should be similar as long as the equation is adjusted correctly. The results differ greatly, unfortunately.

# Minimal code to reproduce problem

I encapsulated the problem into two implementations of equation string, which you can find out below.
Equation forms are based on this paper by Schwartz. Please ignore the form of ionic channels, but pay attention to the evolution equations and their units.

import brian2 as b2
from brian2 import Soma
from brian2.units import *
import matplotlib.pyplot as plt
import numpy as np
from brian2.units.constants import faraday_constant as F
from brian2.units.constants import gas_constant as R

# Some global values, according to Schwartz 1995

El = -84*mV
Ek = -84*mV
Na_out = 154*mmole/liter

Na_in = 35*mmole/liter
gl0 = 30* nsiemens/(1*meter**2)
gKf0 = 15* nsiemens/(1*meter**2)
gKs0 = 30* nsiemens/(1*meter**2)
p_Na0 = 3.52e-9 *cm**3/second/(1*meter**2)
Cm0= 1.4*pF/meter**2
T = (25+273)*kelvin

Ri0 = 0.4*metre/siemens # defines the axial resistance
T0 = (20+273)*kelvin
Qm = 2.2**((T-T0)/(10*kelvin))
Qh = 2.9**((T-T0)/(10*kelvin))
Qn = 3**((T-T0)/(10*kelvin))
Qs = 3**((T-T0)/(10*kelvin))

# this formualtion does not use point current flag
eqs_spatial = '''
# The same equations for the whole neuron, but possibly different parameter values
# distributed transmembrane current
Im = -i_L - i_Ks -i_Kf -i_Na: amp/meter**2

i_L = gl * (v-El) : amp/meter**2

i_Ks = gKs * s * (v-Ek) : amp/meter**2
ds/dt= alpha_s*(1-s) - beta_s * s : 1
alpha_s = Qs*0.00122*23.6/exprel(-(v+12.5*mV)/(23.6*mV))/ms : Hz
beta_s = Qs*0.000739*21.8/exprel((v+80.1*mV)/(21.8*mV))/ms : Hz

i_Kf = gKf * n**4 * (v-Ek): amp/meter**2
dn/dt = alpha_n * (1-n) - beta_n * n : 1
alpha_n = Qn*0.00798*1.1/exprel(-(v+93.2*mV)/(1.1*mV))/ms : Hz
beta_n = Qn*0.0142*10.5/exprel((v+76.*mV)/(10.5*mV))/ms : Hz

i_Na = m**3 * h * p_Na * F*(F*v)/(R*T) * fraction: amp/meter**2
fraction = (Na_in - Na_out *exp(-(F*v)/(R*T)))/(1-exp(-(F*v)/(R*T))): mole/metre**3
dm/dt = alpha_m * (1-m) - beta_m * m : 1
dh/dt = alpha_h * (1-h) - beta_h * h : 1
alpha_m = Qm*1.86*10.3/exprel(-(v+18.4*mV)/(10.3*mV))/ms : Hz
beta_m = Qm*0.086*9.16/exprel((v+22.7*mV)/(9.16*mV))/ms  : Hz
alpha_h = Qh*0.0336*11.0/exprel((v+111.*mV)/(11.*mV))/ms : Hz
beta_h = Qh*2.3/(1+exp(-(v+28.8*mV)/(13.4*mV)))/ms : Hz

gl: siemens/meter**2 (constant)
gKs : siemens/meter**2 (constant)
gKf : siemens/meter**2 (constant)
p_Na: meter/second (constant)
I : amp/meter**2
'''

# but this one does by replacing the last line
eqs_spatial_point = ''.join(eqs_spatial.split('I : amp')[0])
eqs_spatial_point +='''I : amp (point current)'''

I encapsulated all neuron creation, simulation, and plotting in the following functions. The goal is the produce Fig 2Aa of the paper mentioned in the code above.

def make_neuron_and_run(eq_type='point'):
b2.start_scope()
morpho = Soma(diameter=np.sqrt(1/np.pi)*meter) # gives an area of 1 m**2

if eq_type=='point':
model = eqs_spatial_point
correct_unit = 1*nA
else:
model = eqs_spatial # per area
correct_unit = 1*nA/meter**2

neuron = b2.SpatialNeuron(morphology=morpho, model=model,
Cm=Cm0,
Ri=Ri0, # irrelevent for single compartments
method='exponential_euler')

neuron.v = El
neuron.h = 1 # this is critical since we should have no Na inactivation.
neuron.m = 0
neuron.n = .5
neuron.Cm = Cm0
neuron.gKf = gKf0
neuron.gKs = gKs0
neuron.gl = gl0
neuron.p_Na = p_Na0
neuron.I = 0*correct_unit

mon = b2.StateMonitor(neuron, ['v', 'I', 'Im'], record=True)

b2.run(2*ms)
neuron.I = 0.3*correct_unit
b2.run(2*ms)
neuron.I = 0*correct_unit
b2.run(3*ms)

return mon, neuron

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

for i in range(0,mon.v.shape[0]):
axs[0].plot(mon.t,mon.Im[i]*neuron.morphology.area[i], 'k',
alpha=(i+1.)/mon.v.shape[0],label=i, )
axs[0].set_title('Im')

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

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

# originally was desigend for a REAL multicompartmental neuron
for ax in axs:
ax.legend(ncol=6)

# What you have already tried

The formulation with point current can correctly reproduce Fig 6C:

mon, neuron = make_neuron_and_run('point')
plot_vars(mon, neuron)

but the current per area one canâ€™t.

mon, neuron = make_neuron_and_run('per_area')
plot_vars(mon, neuron)

Any input appreciated!

Hi @arashgmn ,
the (point current) flag does two things: it divides a current by the area, but it also adds it to the equation for the Im current. This is missing in your equations with the â€śper areaâ€ť formulation â€“ the current I is simply ignored, since Im does not depend on itâ€¦ When you add + I to the Im equation, it should give the same result.

1 Like

Amazing! Thanks, @mstimberg.

1 Like