Parameter tuning

In this tutorial we will show how to tune a model to a given firing rate. We will use the RectifiedLNP model, but the same procedure can be applied to any other differentiable model.

The RectifiedLNP model is taken from the paper “Systematic errors in connectivity inferred from activity in strongly coupled recurrent circuits”. In the paper it is hand-tuned to get an interspike interval of 16 ms (62.5 Hz firing rate). We will show how to do this automatically using spikeometric.

The tuning procedure is as follows:

  1. Generate a dataset of networks with some connectivity you are interested in

  2. Select a model and initialize its parameters (you might have to try some different initializations to get stable behavior)

  3. Use the tune function to find the parameters that minimize the squared difference between the model’s firing rate and the target firing rate

If you want to fit the parameters to some other target, you can use the model as any other PyTorch module and write your own training loop. The tune() method is just a convenience function that does this for the firing rate.

Preparing the training data

We start by generating some networks that we want to fit the model on. We will use 24 networks of 100 neurons each, with the neurons organized in a ring network with inhibitory weights drawn from a mexican hat distribution.

# Create a dataset of random connectivity matrices to train on
n_neurons = 100
n_networks = 24
dataset = MexicanHatGenerator(100, a=1.0015, sigma_1=6.98, sigma_2=7.).generate(n_networks)
data_loader = DataLoader(dataset, batch_size=5, shuffle=False)

The networks look like this:

../_images/mexican_hat.png

Note that in this case the networks are not random, and all the networks have the same connectivity. It is still useful to use a dataset with several copies, since it allows us to use the DataLoader to iterate over the networks in batches. We will use 5 networks per batch to speed up the training and data generation.

The Rectified Linear-Nonlinear-Poisson model

In the RectifiedLNP model, neurons are subject to a uniform background input \(b\) in addition to the recurrent synaptic input from neighbors, weighted by \(r\). Spikes from neighboring neurons are filtered by an exponential decay with a time constant \(\tau\). For each time step, spikes are sampled from a Poisson distribution after passing the input through a thresholded scaled rectifier.

\[\begin{split}g_i(t+1) &= \sum_{j \in \mathcal{N}(i)} \mathbf{W}_{j, i} \cdot \mathbf{x}_j(t) + b \\ \mu_i(t+1) &= \lambda_0 [g_i(t+1) - \theta]_+ \\ x_i(t+1) &\sim \text{Pois}(\mu_i(t+1))\end{split}\]

We choose the same parameters as in the paper, except for the threshold \(\theta\), which is not specified in the paper and which we will tune.

model = RectifiedLNP(
    lambda_0=100.,          # Scales the output of the rectifier
    theta=1e-4,             # Threshold of the rectifier
    T=20,                   # Length of the time window for the exponential filter (ms)
    tau=10.,                # Time constant of the exponential filter (1/ms)
    dt=0.1,                 # Time step (ms)
    r=0.0025,               # Recurrent weight
    b=0.001                 # Background input
)

Let’s run the model on the test set and see what the firing rates look like:

results = torch.zeros((n_neurons*len(test)), 10000)
for i, data in enumerate(test_loader):
    results[i*data.num_nodes:(i+1)*data.num_nodes] = model.simulate(data, n_steps=10000)

# Average firing rate:  33.88600051403046 Hz
# Average ISI:  29.510713577270508 ms

We get a firing rate of about 34 Hz and an interspike interval of about 29 ms. The plots show the spikes over the first 100 ms for the first network and the ISI histogram for all the test networks.

../_images/tune_model_initial_spikes.png

This is not quite the target firing rate of 62.5 Hz. Let’s see if we can improve it by tuning the threshold parameter \(\theta\).

Tuning the model

We will use the built in tune() method to find the a threshold parameter \(\theta\) that works. The tune method takes a dataset, a target firing rate, and a set of parameter names to tune. The tunable_parameters defaults to all tunable_parameters of the model (you can can check what these are by calling model.tunable_parameters), but you can specify a subset of parameters to tune by passing a list of strings, each of which is the name of a parameter to tune. If you have a stimulus attached to the model, you can also pass "stimulus" to tune the stimulus parameters or "model" to tune the only the model parameters. It also takes optional arguments n_steps, lr and n_epochs to specify the number of steps to run the simulation for each epoch, the learning rate and the number of epochs to run the optimization for.

for data in train_loader:
    model.tune(data, firing_rate=62.5, tunable_parameters=["theta"], n_steps=10, lr=1e-6, n_epochs=100)

After the tuning procedure, the model has been updated with a new value for \(\theta\). Let’s see how the firing rate changed:

results = torch.zeros((n_neurons*len(test)), 10000)
for i, data in enumerate(test_loader):
    results[i*data.num_nodes:(i+1)*data.num_nodes] = model.simulate(data, n_steps=10000)

# New firing rate:  62.26800009608269 Hz
# New ISI:  16.059614181518555 ms

We now get firing rate 62.26 Hz and an interspike interval of 16.05 ms, which is much closer to the target.

../_images/tune_model_final_spikes.png

That concludes this tutorial. We hope you found it useful. If you want to learn how you can also tune the parameters of a stimulus, check out the Working with stimuli tutorial. The tutorial is also available as a juptyer notebook in the examples folder of the repository.