tsodyks(2)_synapse strange behavior in paired-pulse simulations

[HDictus]

Describe the bug
The paired-pulse behavior of both tsodyks_synapse and tsodyks2_synapse is highly unpredictable. While the same parameters, simulation settings, and interval lead to the same strength (it is deterministic), the difference between intervals with a given parameterization is fully chaotic. When I do trains of multiple spikes at 10Hz the behavior of the synapse is reasonable, but paired-pulse simulations specifically seem to have a problem.

To Reproduce
recreate_bug.py :

import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
import nest


plots_dir = Path('intervalplots')
plots_dir.mkdir(exist_ok=True)

times_run = [0]
def sim(interval, nmodel, modparams, receptor_type):
    nest.ResetKernel()
    times_run[0] += 1
    nest.rng_seed = times_run[0]
    src = nest.Create('spike_train_injector')
    src.spike_times = [1., 1.+interval]
    tgt = nest.Create(nmodel, modparams)
    synapse_model = {
        'synapse_model': 'tsodyks2_synapse',
        'u': 0.5,
        'U': 0.5,
        'receptor_type': receptor_type,
        'tau_fac': 0.0,
        'tau_rec': 1600
    }
    syn = nest.Connect(src, tgt, 'one_to_one', synapse_model, return_synapsecollection=True)
    mm = nest.Create('multimeter', params={'interval': 0.25, 'record_from': ['g_1', 'V_m']})
    nest.Connect(mm, tgt, 'one_to_one')
    nest.Simulate(interval + 50)
    df = pd.DataFrame(mm.get()['events'])
    plt.clf()
    plt.plot(df['times'], df['g_1'])
    plt.savefig(plots_dir / f'{nmodel}-{interval}.png')
    return df

stds = {}
peaks = {}
for interval in (5, 10, 11, 12, 13, 14, 15, 20, 30, 40, 50, 75, 100, 150, 200, 300, 400, 600, 800):
    pks = []
    for _ in range(10):
        df = sim(interval, 'glif_cond', {'tau_syn': [0.1, 0.1]}, 1)
        pks.append(df['g_1'][df['times'] > interval].max())

    stds[interval] = np.std(pks)
    peaks[interval] = np.mean(pks)

plt.clf()
plt.plot(peaks.keys(), peaks.values())
plt.ylabel('second peak amp')
plt.xlabel('interval')
plt.savefig(plots_dir / 'peaks.png')

for interval in (5, 10, 15, 20, 30, 40, 50, 75, 100, 150, 200):
    sim(interval, 'aeif_cond_alpha_multisynapse', {'tau_syn': [0.1]}, 1)

print(stds)

This script checks that the problem occurs for more than one neuron model.

run python recreate_bug.py and check the figures in 'intervalplots'

You can play with the parameters. I found this behavior to persist regardless of parameters, so long as tau_syn << interval, and it was also there with tsodyks_synapse, and persisted after I reinstalled NEST.

Expected behavior
Given that facilitation is disabled, the amplitude of the second pulse should be most strongly reduced at small intervals, and exponentially decay back to baseline with longer intervals.

Desktop/Environment (please complete the following information):

• OS: Red Hat Enterprise Linux 9.4 (Plow)
• Shell: bash
• Python-Version: 3.11.7
• NEST version 3.8.0-post0.dev0
• Installation: built from source, no summary was shown at the end. regular install instructions did not work on the cluster, enlisted help from IT staff, installed using the commands:

module load gcc openmpi python cmake gsl boost hdf5/1.14.3-mpi py-mpi4py py-cython
python3 -m venv --system-site-packages venv
source venv/bin/activate
cmake -B build -Dwith-mpi=ON -Dwith-hdf5=ON
srun -c 18 -q debug --pty cmake --build build -j18
cmake --install build

Where in the outputs can I find the neccessary information about the install?

[clinssen]

Hi, this isn't yet a fix as such, but I do notice that the issue goes away when you set the simulation resolution to 0.01, i.e. one order of magnitude lower than the synaptic time constant. Increasing the synaptic time constant also seems to alleviate the issue. Could it be due to a too-coarse sampling of the multimeter of the recorded timeseries, since it involves a very brief peak?

[heplesser]

Hi @HDictus, here @clinssen had a good hunch: The problem is the recording interval for the multimeter. By changing the resolution to the NEST resolution with

 mm = nest.Create('multimeter', params={'interval': nest.resolution, 'record_from': ['g_1', 'V_m']})

I get the curve attached below. Does this look like you expect?

What I find disconcerting here is that multimeter accepts an interval value of 0.25 ms while the simulation resolution is 0.1 ms. It silently rounds the interval up to 0.3 ms.

[HDictus]

That does remove the wild variability, thanks. The curve you show looks like what I'd expect. However, my mistake with the multimeter must have been obscuring another issue, because on my device I get a flat line.

I've fiddled with resolution and parameters a bit and keep seeing the same result. I'm afraid I don't have much to go on to figure out what the issue is.

[heplesser]

Strange that you get a flat curve. Here is the code I ran, just for reference. It should be identical to the code you posted above except for the multimeter interval. I am running with recent master (nest-3.9.0-post0.dev0). I noticed that you had a corresponding 3.8.0-post0.dev0. Could you try with either the release version of NEST 3.9 or current master?

import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np
import nest


plots_dir = Path('intervalplots')
plots_dir.mkdir(exist_ok=True)

times_run = [0]
def sim(interval, nmodel, modparams, receptor_type):
    nest.ResetKernel()
    times_run[0] += 1
    nest.rng_seed = times_run[0]
    src = nest.Create('spike_train_injector')
    src.spike_times = [1., 1.+interval]
    tgt = nest.Create(nmodel, modparams)
    synapse_model = {
        'synapse_model': 'tsodyks2_synapse',
        'u': 0.5,
        'U': 0.5,
        'receptor_type': receptor_type,
        'tau_fac': 0.0,
        'tau_rec': 1600
    }
    syn = nest.Connect(src, tgt, 'one_to_one', synapse_model, return_synapsecollection=True)
    mm = nest.Create('multimeter', params={'interval': nest.resolution, 'record_from': ['g_1', 'V_m']})
    nest.Connect(mm, tgt, 'one_to_one')
    nest.Simulate(interval + 50)
    df = pd.DataFrame(mm.get()['events'])
    plt.clf()
    plt.plot(df['times'], df['g_1'])
    plt.savefig(plots_dir / f'{nmodel}-{interval}.png')
    return df

stds = {}
peaks = {}
for interval in (5, 10, 11, 12, 13, 14, 15, 20, 30, 40, 50, 75, 100, 150, 200, 300, 400, 600, 800):
    pks = []
    for _ in range(10):
        df = sim(interval, 'glif_cond', {'tau_syn': [0.1, 0.1]}, 1)
        pks.append(df['g_1'][df['times'] > interval].max())

    stds[interval] = np.std(pks)
    peaks[interval] = np.mean(pks)

plt.clf()
plt.plot(peaks.keys(), peaks.values(), 'o-')
plt.ylabel('second peak amp')
plt.xlabel('interval')
plt.savefig(plots_dir / 'peaks.png')

for interval in (5, 10, 15, 20, 30, 40, 50, 75, 100, 150, 200):
    sim(interval, 'aeif_cond_alpha_multisynapse', {'tau_syn': [0.1]}, 1)

print(stds)

[HDictus]

Upgrading to current master fixes the issue, thank you.