`Equations`: how to reporduce `point current` in a multicompratmental neurons?

Support
,
1

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.

fig2Aa

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)

point

but the current per area one can’t.

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

area

Any input appreciated!

2

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
3

Amazing! Thanks, @mstimberg.

1 Like