# 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:

- 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 - 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!