NMODL alpha-synapse (missing at_time)

[espenhgn]

Hi,
I'm trying to define the common alpha synapse model in the Arbor NMODL dialect and encountered some issue. I'm basing this on the NEURON version found here: https://github.com/neuronsimulator/nrn/blob/99a6a370ef143b3c681eda1deea68bc7e2b8f21d/src/nrnoc/syn.mod

After adding t as to the PARAMETER block and removing i (nA) from the ASSIGNED block,
my problem is that Arbor don't define the at_time function used in the BREAKPOINT block, hence I fail to determine the time relative to synapse activation (t-onset).

What is the correct way to deal with such discontinuities that occur at time of activations?

Edit: The functionality of at_time is not super well described, but it's defined as: https://github.com/neuronsimulator/nrn/blob/99a6a370ef143b3c681eda1deea68bc7e2b8f21d/src/nrncvode/cvodeobj.cpp#L746.
I assume at_time(onset) will evaluate to True at synapse activation times in this case.
Edit 2: The NEURON book from 2005 notes that at_time() has become "deprecated" (Ch. 9.4.2.1). Funny.

My version of the file is:

COMMENT
  an synaptic current with alpha function conductance defined by
        i = g * (v - e)      i(nanoamps), g(microsiemens);
        where
         g = 0 for t < onset and
         g = gmax * (t - onset)/tau * exp(-(t - onset - tau)/tau)
          for t > onset
  this has the property that the maximum value is gmax and occurs at
    t = delay + tau.

 Rewritten for Arbor from:
 https://github.com/neuronsimulator/nrn/blob/99a6a370ef143b3c681eda1deea68bc7e2b8f21d/src/nrnoc/syn.mod

ENDCOMMENT

NEURON {
	POINT_PROCESS alphasyn
	RANGE onset, tau, gmax, e, i
	NONSPECIFIC_CURRENT i
}
UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(uS) = (microsiemens)
}

PARAMETER {
	onset=0 (ms)
	tau=.1 (ms)	<1e-3,1e6>
	gmax=0 	(uS)	<0,1e9>
	e=0	(mV)
        t
}

ASSIGNED { v (mV) g (uS) }

BREAKPOINT {
  if (gmax) {
    at_time(onset)
  }
  g = gmax * alpha( (t - onset)/tau )
  i = g*(v - e)
}

FUNCTION alpha(x) {
	if (x < 0 || x > 10) {
		alpha = 0
	}else{
		alpha = x * exp(1 - x)
	}
}

I have this added to the mechanisms/default folder as alphasyn.mod and added alphasyn to the corresponding CMakeLists.txt file, reinstalling Arbor simply using pip (obviously now failing though).

The code I set up to check my implementation is:

class Recipe (arbor.recipe):
    def __init__(self, cell, times=[[10.]], weights=[1E-3], synapse_type='expsyn'):
        super().__init__()
     
        self.the_cell = cell
        self.vprobe_id = (0, 0)
        self.cprobe_id = (0, 1)
        self.the_props = arbor.neuron_cable_properties()
        self.the_cat = arbor.default_catalogue()
        self.the_props.register(self.the_cat)

        assert(len(times) == len(weights)), 'len(times) != len(weights)'
        self.times = times
        self.weights = weights
        self.synapse_type = synapse_type
        
    def num_cells(self):
        return 1

    def num_sources(self, gid):
        return 0

    def num_targets(self, gid):
        return len(self.times)

    def connections_on(self, gid):
        return []

    def event_generators(self, gid):
        events = []
        for i, (w, t) in enumerate(zip(self.weights, self.times)):
            events += [arbor.event_generator(f'{i}', w, 
                                             arbor.explicit_schedule(t))]        
        return events
    
    def cell_kind(self, gid):
        return arbor.cell_kind.cable

    def cell_description(self, gid):
        return self.the_cell

    def global_properties(self, kind):
        return self.the_props
    
    def probes(self, gid):
        return [
            arbor.cable_probe_membrane_voltage_cell(),
            arbor.cable_probe_point_state(0, self.synapse_type, 'g')
        ]
    
# define segment tree
tree = arbor.segment_tree()
tree.append(arbor.mnpos, 
            arbor.mpoint(0, 0, 0, 10), 
            arbor.mpoint(0, 0, 10, 10), tag=1)

# define morphology (needed for arbor.place_pwlin)
morphology = arbor.morphology(tree)

# define relative locations (segment midpoints)
nseg = 1
locations = (np.arange(nseg) + 0.5) / nseg

# Label dictionary
labels = arbor.label_dict()
labels['all'] = '(all)'

# decor
decor = arbor.decor()

# set initial voltage, temperature, axial resistivity, membrane capacitance
decor.set_property(
    Vm=-65,  # Initial membrane voltage [mV]
    tempK=300,  # Temperature [Kelvin]
    rL=1000,  # Axial resistivity [Ω cm]
    cm=0.01,  # Membrane capacitance [F/m**2]
)

# set passive mechanism all over
pas = arbor.mechanism('pas/e=-65')  # passive mech w. leak reversal potential (mV)
pas.set('g', 0.0001)  # leak conductivity (S/cm2)
decor.paint('"all"', pas)

# create synapse
times = [[10.]]
weights = [1E-3]
syn_loc_set = '(location 0 0.5)'
synapse_type = 'alphasyn'
synapse_params = {'gmax': 1E-3, 'tau': 2}
decor.place(syn_loc_set, arbor.mechanism(synapse_type, synapse_params), '0')

# number of CVs per branch
policy = arbor.cv_policy_fixed_per_branch(nseg, '"all"')
decor.discretization(policy)

# create cell and set properties
cell = arbor.cable_cell(tree, labels, decor)

# create single cell model
model = arbor.single_cell_model(cell)

# instantiate recipe with cell
recipe = Recipe(cell, times, weights, synapse_type)

# instantiate simulation
context = arbor.context()
domains = arbor.partition_load_balance(recipe, context)
sim = arbor.simulation(recipe, domains, context)

# set up sampling on probes
schedule = arbor.regular_schedule(0.1)
v_handle = sim.sample(recipe.vprobe_id, schedule, arbor.sampling_policy.exact)
c_handle = sim.sample(recipe.cprobe_id, schedule, arbor.sampling_policy.exact)

# run simulation for 500 ms of simulated activity and collect results.
sim.run(tfinal=40)

# extract time, V_m and I_m for each compartment
V_m_samples, V_m_meta = sim.samples(v_handle)[0]
I_c_samples, I_c_meta = sim.samples(c_handle)[0]

# prep recorded data for plotting
time = I_c_samples[:, 0]
I_c = I_c_samples[:, 1:].ravel()  # synapse conductance
V_m = V_m_samples[:, 1:].ravel() # postsynaptic potential

fig, axes = plt.subplots(2, 1, figsize=(10, 10), sharex=True)
axes[0].plot(time, I_c)
axes[1].plot(time, V_m)

[espenhgn]

Hmm, this may be a non issue and a waste of time. The exp2syn should well mimic the alpha synapse if tau2 is fractionally larger than tau1

[halfflat]

The limitation in Arbor mechanisms (as they currently are defined) is that we don't explicitly expose time as a parameter for use in BREAKPOINT. One can work around this though by adding a state variable u say which has derivative 1, e.g.

STATE {
    u
}
INITIAL {
    u = -onset
}
BREAKPOINT {
     SOLVE state METHOD cnexp
     g = gmax * alpha(u / tau)
     i = g*(v-e)
}
DERIVATIVE state {
    u' = 1
}

We don't support at_time; NEURON uses it to indicate to a variable time step solver that there is a discontinuity at a particular point in time but we only have a solver that respects the specified max dt in the simulation, and outside of event handling, we don't account for discontinuities in the synapse model.

Exact timing can be recovered though by removing onset from the model, and instead using an explicit event (from an event generator in the recipe) to trigger the synapse.

While the description in terms of the exponential function will work, a more natural definition of the alpha model though would be to use the differential form. The provided exp2syn synapse is close, as you noted, but doesn't apply in the case where the two time constants are equal. But t*exp(-t) is a solution to the second order differential equation y'' +2y' + y = 0, so we should be able to model it with a first order ODE with two state variables. I haven't verified the following, but I think this should be the equivalent to the exp2syn with equal time constants:

PARAMETER {
    tau = 0.1
    e = 0
}
STATE {
    u w
}
INITIAL {
    u = 0
    w = 0
}
BREAKPOINT {
    SOLVE state METHOD sparse
    i = u*(v-e)
}
DERIVATIVE state {
    u' = w/tau
    w' = -(2*w+u)/tau
}
NET_RECEIVE(weight) {
    w = w+tau*weight
}

This is the event-based version, where the connection weight plays the role of gmax, and the effect of successive events is accumulative.

[espenhgn]

Thanks! Your suggestion produces the correct temporal profile, but I believe the correct factor in the NET_RECEIVE part is exp(1) not tau that's all (at least to a couple digits precision):

The full .mod file is then

COMMENT
Bla Bla
ENDCOMMENT

NEURON {
	POINT_PROCESS alphasyn
	RANGE tau, e, i
	NONSPECIFIC_CURRENT i
}

PARAMETER {
    tau = 0.1
    e = 0
}

STATE {
    g w
}

INITIAL {
    g = 0
    w = 0
}

BREAKPOINT {
    SOLVE state METHOD sparse
    i = g * (v - e)
}

DERIVATIVE state {
    g' = w / tau
    w' = -(2 * w + g) / tau
}

NET_RECEIVE(weight) {
    w = w + exp(1) * weight
}

[espenhgn]

Think this can be marked as "answered"