# Matching a spike-response-like model in Brian

**URL:** <https://brian.discourse.group/t/matching-a-spike-response-like-model-in-brian/398>\
**Category:** Science, Projects, and Showcases\
**Tags:** equations\
**Created:** [2 June 2021 16:03 UTC](https://brian.discourse.group/t/matching-a-spike-response-like-model-in-brian/398 "2021-06-02T16:03:51Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![oliverdibb](https://avatars.discourse-cdn.com/v4/letter/o/9e8a1a/32.png) [@oliverdibb](https://brian.discourse.group/u/oliverdibb)\
**Post date:** [2 June 2021 16:03 UTC](https://brian.discourse.group/t/matching-a-spike-response-like-model-in-brian/398/1 "2021-06-02T16:03:51Z")

</div>

Hi!

First - I’ve placed this in the Science section because it’s mainly about model choice, but it also revolves around implementation in Brian, so apologies if it’s misplaced.

I’m trying to reproduce a [paper](https://academic.oup.com/cercor/article/30/3/952/5547080) which uses a spike-response-like neuron model. In the paper, the membrane potential of an excitatory neuron i is given as

\begin{align} u\_i(t) = \sum\_{j} w^{inpE}\_{ij} y\_j(t) + \sum\_{j} w^{EE}\_{ij}z\_j(t) - \sum\_{j}w^{IE}\_{ij}h(t) + E\_{gen} + E\_i, \end{align}

where the PSPs are given as

\begin{align} y(t) = \sum\_{t\_f} \epsilon(t - t\_f) \end{align}

for presynaptic spike times \{t\_f\}, and response kernel  
\begin{align} \epsilon(t) = K\left( e^{-t/\tau\_f} - e^{-t/\tau\_r}\right), \end{align}  
where the constant K is chosen for a peak value of 1. The quantities E\_{gen} and E\_i represent generic excitability of the excitatory population and the random individual excitability of neuron i (drawn lognormally) respectively.  
Summed postsynaptic potentials z and h are defined identically, but for different presynaptic spike times (from excitatory and inhibitory afferents respectively).

The original model has a stochastic firing criterion (without an adaptive threshold), with an instantaneous firing rate given by the transfer function  
\begin{align} r(t) = r\_0 \exp\left(u(t)/4\right). \end{align}

The original paper was implemented in NEST, and I’m fairly certain they used the model `pp_psc_delta` given [here](https://nest-simulator.readthedocs.io/en/latest/models/pp_psc_delta.html). Model parameters are set to the mean values in [this](https://link.springer.com/article/10.1007/s10827-006-7074-5) paper. Note that with the given transfer function, this model would correspond to the `pp_psc_delta` model with `c_1 = 0`.

In replicating the model in Brian, I’m using the equations given in the documentation:

```auto
dv/dt = (A*s - v + exc_gen + exc_i) / tau_f : 1
ds/dt = -s / tau_r : 1

```

with presynaptic spike update

```auto
s_post += w.

```

Here, `A` is the scaling factor (a function of the time constants \tau\_f, \tau\_r, given in the documentation) ensuring a peak PSP value of 1, and `exc_gen` and `exc_i` are the quantities E\_{gen} and E\_i.

Since thresholds can only be given as explicit spike conditions, I used that (given in the NEST model documentation) the instantaneous spike probability of the NEST model is given as

P(\text{spike}) = 1 - \exp(-r h) for instantaneous firing rate r and simulation time step h.  
I matched this with the instantaneous probability of firing for a model with threshold \theta = V\_{\text{th}} + \xi, for soft threshold V\_{\text{th}} and noise \xi from a distribution with pdf \rho(\xi):  
P(\text{spike}) = \int\_{\mathbb{R}} H(V - V\_{th} - \xi)\rho(\xi) d\xi

Below, see the v - P(\text{spike}) curve from the original model (blue) vs v - P(\text{spike}) curves from a matched normal distribution (orange, dashed) and a matched skewnormal distribution (green, dashdot), along with the estimated PDF of the CDF given by the v-P(\text{spike}) curve from the original model.

 ![threshold_noise](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/1563c3a8a549d3d74e5d249d0d6b6cd6e780882f.png)

The threshold is implemented one of two ways: either I use a Gaussian noisy threshold with  
`threshold = 'v > theta + sigma*randn()'` (for \theta, \sigma of the matched normal distribution), or a skewnorm threshold implemented as a network operation, where threshold values are generated as variates from the skewnormal distribution matched to the estimated PDF shown above.

The model is not behaving as it should, and in going through the network model piece by piece to check for errors, I’ve noticed a few areas where I’m not certain about choices made in matching the model to the paper.

First, I’m not sure how the spike and reset mechanism of spike-response models changes the details of implementation in Brian. In my current implementation, the membrane potential is reset to  
`v = exc_gen + exc_I` after a spike. However in some sources (e.g. in the Jolivet et al. paper linked above), spike-response models are implied to not explicitly reset the membrane potential after a spike.

Second, I’ve assumed that the v-P(\text{spike}) curves are CDFs and matched the threshold noise distributions to the distribution whose CDF matches the given v-P(\text{spike}) curve. This assumption might be wrong, but I don’t feel confident on why it might or might not be.

Finally, on the implementation level - are there any ways to implement a skewnorm distributed threshold (or in general, any non-normally distributed noisy threshold) without preventing the use of standalone mode? The current implementation through network operations does not allow this.

I hope this post hasn’t been too rambling - any comments or input on the steps I’ve taken to match the model would be greatly appreciated.

Thanks!

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [16 June 2021 14:59 UTC](https://brian.discourse.group/t/matching-a-spike-response-like-model-in-brian/398/2 "2021-06-16T14:59:16Z")

</div>

Hi. Sorry for the very late reply, you might have solved some of the issues already… I have to say I am having a bit of a hard time making a link between the description in the paper and the one in the nest documentation (and of both with more common integrate-and-fire models)…  
Regarding the threshold: if r = r\_0 \exp(u/4) and P(\mathrm{spike}) = 1 - \exp(-r h) then I think the most straightforward implementation in Brian would be to use

```auto
threshold='rand() < 1 - exp(-r_0*exp(u/4)*dt)'

```

As the threshold definition, instead of introducing a noisy threshold.  
Regarding the membrane potential: if I understand the model correctly, the membrane potential is an instantaneous function of the inputs (+ some bias term) and not a dynamical variable defined by a differential equation. So the equations would be something along the lines of

```auto
'''
u = psp + exc_gen + exc_i : 1
dpsp/dt = (A*s - psp)/tau_f : 1
ds/dt = -s/tau_r : 1
'''

```

From their paper, I don’t think there’s anything like a reset (or a reset/AHP kernel used in other spike-response models), there’s just a refractory period. Since this period is different for each neuron, it should be a parameter in the equations (e.g. `ref_period : second (constant)`), and then used as `refractoriness='ref_period'`. This assumes that each neuron has a fixed refractory period, their text is a bit ambiguous about that.

Let me know whether this gets you started!