fwm_behaviour_1.ipynb

fwm_utils.py

import collections
import numpy as np


Population = collections.namedtuple('Population', 'n_neurons t_rc t_ref gains biases')
SimStep = collections.namedtuple('SimStep', 'voltage refs spikes')


def generate_population(n_neurons, t_rc, t_ref, max_rate):
    """Generate a population for simulation using the given parameters."""
    # Generate the max_rates for the population, and the intercepts (in range
    # +/- 1)
    max_rates = np.random.uniform(0, max_rate, size=n_neurons)
    intercepts = np.linspace(-1., 1., n_neurons, endpoint=False)
    z = 1. / (1 - np.exp((t_ref - 1./max_rates)/t_rc))

    gains = (1. - z)/(intercepts - 1.)
    biases = 1. - gains*intercepts

    # Generate and return the population
    return Population(n_neurons, t_rc, t_ref, gains, biases)


def simulate_population_step(population, in_currents, pstep, dt):
    """Simulate a population of neurons given the input current and the
    previous simulation step.
    """
    # Add the biases to the input currents
    in_currents += population.biases

    # Calculate dV and the subsequent voltage
    dv = (in_currents - pstep.voltage)*dt / population.t_rc
    dv[pstep.refs < 0] = 0
    voltage = pstep.voltage + dv
    voltage[voltage < 0] = 0.

    # Decrease refractory periods
    refs = pstep.refs
    refs[refs > 0.] -= 1

    # Deal with spikes
    spikes = voltage > 1.
    voltage[spikes] = 0.
    pets = np.random.rand(population.n_neurons) < (voltage - 1.)/dv
    refs[spikes] = int((population.t_ref) / dt)
    refs[np.all(np.vstack([spikes, pets]), axis=0)] -= 1

    # Return the simulation step
    return SimStep(voltage, refs, spikes)


def simulate_feedforward(population, inputs, weight_matrix, dt=0.001,
                         filter_tc=0.005):
    """Simulate a feedforward connection from the spiking inputs to the
    population of neurons.

    `inputs` is expected to be a m-by-t matrix of input values.
    """
    assert weight_matrix.shape == (inputs.shape[0], population.n_neurons)
    spike_times = [list([-1.]) for n in range(population.n_neurons)]

    # Create a series of filter holder values
    synaptic_filters = np.zeros(population.n_neurons)
    f_m1 = np.exp(-dt/filter_tc)
    f_m2 = 1. - f_m1

    # For each simulation step (i.e., each column of the input) apply the input
    # to the population.
    pstep = SimStep(np.zeros(population.n_neurons),
                    np.zeros(population.n_neurons), [])
    for t, inp in enumerate(inputs.T):
        # Apply any spikes to the synapses to get the current contributions for
        # this simulation step.
        ins = inp.dot(weight_matrix)
        synaptic_filters *= f_m2
        synaptic_filters += f_m1 * ins

        # Simulate the neurons for this step
        pstep = simulate_population_step(
            population, np.copy(synaptic_filters), pstep, dt)

        # Save the spike times
        for n in np.where(pstep.spikes)[0]:
            spike_times[n].append(t * dt)

    # Return the spike times
    return spike_times


def get_interspike_times(spike_times):
    """Get a list of interspike times for each neuron."""
    interspike_times = list()

    for spikes in spike_times:
        interspike_times.append(
            [spikes[n] - spikes[n-1] for n in range(1, len(spikes))]
        )

    return interspike_times