Injecting current spike in LIF neurons

Support
1

Hi everyone!

I’m fairly new to Brian so I hope I’m posting this correctly. Please let me know if that’s not the case!

I’m simulating a basic network N of n neurons that are not connected to each other. The neurons are LIF and are set to fire using a bias current i_offset . I’m trying to simulate events that cause a sudden spike in the bias current to a random neuron in the network. To do so I’ve tried making a network operation that takes the i_offset of a certain neuron and adds a spike to it using a function of t:

 @network_operation(dt=500 * ms)
  def event():
        neuron = np.random.randint(low=0, high=n, dtype=int)
        i_offset_neuron = N.i_offset[neuron]
        base_i_offset = i_offset_neuron
        new_i_offset = 'base_i_offset + (10 / (1 + 1000 * t / second)) * nA'
        N.i_offset[neuron] = new_i_offset

The code does run but I don’t see any change in spiking behavior of the neuron (in either spike or statemonitor). Can anyone point out why this isn’t working as intended?
The goal would be to have a spike added to the existing i_offset. Also, I want to be able to quickly change the width and height of the spike (therefore I’m not using dirac delta functions).

Best,

Joris

2

Hi @JFRWKievits and welcome to the Brian forum :wave: . Could you please also post the definition of your NeuronGroup to see how i_offset is applied? Thanks!

3

Hi @mstimberg !

Of course!
Here are the equations and Network setup. The gaussian function just provides a list of len(N) with normally distributed values around some ideal parameter value.

   eqs = """
    R : ohm
    i_offset : amp
    tau : second
    Vr : volt
    Vt : volt
    dv/dt = (R*i_offset-v)/tau : volt (unless refractory)
    """

    # The neurons with threshold voltage, reset voltage and refractoriness
    N = NeuronGroup(n, eqs, threshold='v > Vt', reset='v = Vr', refractory=t_ref, method='euler', name='neurons')
    N.Vr = gaussian(Vr_ideal, Vr_ideal * f, len(N))
    N.Vt = gaussian(Vt_ideal, Vt_ideal * f, len(N))
    N.tau = gaussian(tau_ideal, tau_ideal * f, len(N))
    N.R = N.tau / capacitance
    N.i_offset = 20 * nA

I hope this helps!

4

Hi again. So I hope I understand correctly what you want to do: every 500ms, you want to randomly select a neuron and have it’s i_offset current jump and then decay as an inverse of time, correct?

In your current network_operation, you assign to N.i_offset. An assignment to a variable sets its current value. In your case this is a string referring to time t, but this string will be evaluated once when you do the assignment. According to your model equations, i_offset is a parameter, not a dynamical variable, so its value will not change after that assignment. Also, note that t refers to the time since the start of the simulation, so at t=500*ms, you’ll increase your i_offsetonly by 0.02nA.

There are several ways you can implement a dynamically decaying current spike. The easiest is to have it decay exponentially by defining it as a differential equation that decays back to the base value with time constant tau_i (changing this time constant allows you to change the width of the spike).

tau_i = 2*ms
base_i_offset = 20*nA
eqs = """
  R : ohm
  di_offset/dtau_i = (base_i_offset - i_offset)/tau_i : amp
  tau : second
  Vr : volt
  Vt : volt
  dv/dt = (R*i_offset-v)/tau : volt (unless refractory)
  """

# The neurons with threshold voltage, reset voltage and refractoriness
N = NeuronGroup(n, eqs, threshold='v > Vt', reset='v = Vr', refractory=t_ref, 
                method='euler', name='neurons')

You can then make the i_offset value jump up every 500ms. You could do this via a network_operation, but it might be more straightforward to pre-generate the indices of the neurons that should receive this spike, i.e. something like this:

n_spikes = 10
spike_trigger = SpikeGeneratorGroup(n, np.random.randint(n, size=n_spikes), 500*ms*(np.arange(n_spikes) + 1)
trigger = Synapses(spike_trigger, N, on_pre='i_offset += 10*nA')

If it is important to have an \frac{1}{1 + \alpha t} shape of the current, then you can put it directly into the definition of the current:

base_i_offset = 20*nA
eqs = """
  R : ohm
  i_offset = base_i_offset + (10 / (1 + 1000 * (t - last_injection) / second)) * nA : amp
  last_injection : second  # time of the last current injection for this neuron
  tau : second
  Vr : volt
  Vt : volt
  dv/dt = (R*i_offset-v)/tau : volt (unless refractory)
  """
N = ...
N.last_injection = -1e9*second

Now, instead of directly increasing the current each time, you’ll just update the last_injection time:

spike_trigger = SpikeGeneratorGroup(...)
trigger = Synapses(spike_trigger, N, on_pre='last_injection = t')

I hope this all makes sense and I did not misunderstand what you are trying to do!

5

@mstimberg this is fantastic, thank you so much!

I have implemented it as a network_operation and it works beautifully.

Best,

Joris