initial commit

README.md

@@ -1,2 +1,31 @@
-# SpikingNeuralProgramForagingInsect-PNAS
-A spiking neural program for sensorimotor control during foraging in flying insects.
+# A spiking neural program for sensorimotor control during foraging in flying insects
+
+This repository contains all accompanying code and allows to re-generate the results for the following paper:
+
+*H. Rapp, MP. Nawrot, A spiking neural program for sensorimotor control during foraging in flying insects.*
+
+If you use any parts of this code for your own work, please acknowledge us by citing the above paper.
+
+
+If you have questions or encounter any problems while using this code base, feel free to file a Github issue here and we'll be in touch !
+
+# project layout
+This project uses a mixed code base of Python and MATLAB scripts. Python and BRIAN2 is used for the spiking neural network (SNN) models and simulations thereof. The simulation results (spike trains) are dumped as numpy pickled files (NPZ) and MATLAB (MAT) files.
+
+MATLAB is used for all learning and memory experiments to train the Multispike Tempotron readout neuron on the dumped spike trains from the model simulations and for most data analysis and figures.
+
+All script files are commented and/or self-explanatory.
+
+* `./` root folder contains all Python and BASH scripts to run the SNN simulations to re-generate the data used for the paper (Note: this requires several TB of disk space and large amount of RAM !)
+* `olnet/models/` contains the BRIAN2 Mushroom body SNN model definitions
+* `olnet/plotting/` contains matplotlib scripts to plot SNN network activity
+* `matlab/` contains all the MATLAB code for fitting the readout neuron, data analysis and figures
+
+
+# Using the model
+The BRIAN2 model definition is located in the file `olnet/droso_mushroombody.py` for the model without APL neuron and `olnet/droso_mushroombody_apl.py` for the model with APL neuron. If you want to use our model for your own study, just import the model definitions from these files.
+Both models use the same interface such that you can easily swap out the implementations by modifying your `import` statement. For an example on how to use it see the `run_model` function in `mkDataSet_DrosoCustomProtocol.py`.
+
+
+# usage
+Detailed description of usage of the individual script files will be published here shortly.

figures.mplstyle

@@ -0,0 +1,9 @@
+axes.linewidth : 1
+xtick.labelsize : 8
+ytick.labelsize : 8
+axes.labelsize : 8
+lines.linewidth : 1
+lines.markersize : 2
+legend.frameon : False
+legend.fontsize : 8
+axes.prop_cycle : cycler(color=['e41a1c', '377eb8', '4daf4a', '984ea3', 'c51b7d', '4d9221', '542788', '8c510a', 'b2182b', '2166ac', '01665e'])

make_labdcond.sh

@@ -0,0 +1,30 @@
+#!/bin/bash
+PYTHON_BIN="python3.6"
+N_CPU=12
+# default sparsity condition
+$PYTHON_BIN mkDataSet_DrosoLabCondition.py --name LabCondConnectivityHighSparsityAPL_0-15-3sec \
+-N 30 --n_cpu $N_CPU -T 3 --max_pulse_duration 0.5 --min_pulse_duration 0.1 \
+--stim_noise_scale 0.004 --bg_noise_scale 0.0055 \
+--odor_ids 0 --odor_ids 15 \
+-o data/LabCondConnectivityHighSparsityAPL_0-15-3sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoLabCondition.py --name LabCondConnectivityMediumSparsityAPL_0-15-3sec \
+-N 30 --n_cpu $N_CPU -T 3 --max_pulse_duration 0.5 --min_pulse_duration 0.1 \
+--stim_noise_scale 0.004 --bg_noise_scale 0.0055 \
+--odor_ids 0 --odor_ids 15 \
+--modelParams PNperKC=8.1 \
+-o data/LabCondConnectivityMediumSparsityAPL_0-15-3sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoLabCondition.py --name LabCondWeightMediumSparsityAPL_0-15-3sec \
+-N 30 --n_cpu $N_CPU -T 3 --max_pulse_duration 0.5 --min_pulse_duration 0.1 \
+--stim_noise_scale 0.004 --bg_noise_scale 0.0055 \
+--odor_ids 0 --odor_ids 15 \
+--modelParams wPNKC=20 \
+-o data/LabCondWeightMediumSparsityAPL_0-15-3sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoLabCondition.py --name LabCondConnectivityLowSparsityAPL_0-15-3sec \
+-N 30 --n_cpu $N_CPU -T 3 --max_pulse_duration 0.5 --min_pulse_duration 0.1 \
+--stim_noise_scale 0.004 --bg_noise_scale 0.0055 \
+--odor_ids 0 --odor_ids 15 \
+--modelParams PNperKC=12 \
+-o data/LabCondConnectivityLowSparsityAPL_0-15-3sec.mat

make_paper_figures.py

@@ -0,0 +1,30 @@
+from olnet.plotting.figures import figure1
+import numpy as np
+import matplotlib.pyplot as plt
+fileType = "png"
+
+# plot LabConditioning single-trial
+file = 'cache/LabCond_0-3-5-8-15-3sec/sim-odor-0-0-58.npz'
+mstMATFile = 'matlab/model_cache/predictions/msp_classicalLabCond-0-15.odor-0.1-sp.1/LabCond_0-3-5-8-15-3sec.mat'
+data = np.load(file)['data'][()]
+figure_1 = figure1(data, t_min=1.0, t_max=1.3, orn_range=[620,680], pn_range=[0,35], cmap='seismic',
+                   mstMatFile=mstMATFile, mstOdorIdx=0, mstTrialIdx=0, fig_size=(3.5, 6))
+figure_1.savefig("figures/system_response_labcond.{}".format(fileType), dpi=300)
+
+file = 'cache/PoisonPulse_0-3-5-8-15-10sec/sim-12-90.npz'
+mstMATFile = 'matlab/model_cache/predictions/msp_classicalLabCond-0-15.odor-0.1-sp.1/PoisonPulse_0-3-5-8-15-10sec.mat'
+data = np.load(file)['data'][()]
+figure_2 = figure1(data, t_max=8, orn_range=[640,665], pn_range=[0,25], cmap='seismic',mstMatFile=mstMATFile, mstOdorIdx=0, mstTrialIdx=12)
+figure_2.savefig("figures/system_response_poisson.{}".format(fileType), dpi=300)
+
+figure_2_alt = figure1(data, t_max=8, orn_range=-1, pn_range=[0,45], cmap='seismic',mstMatFile=mstMATFile, mstOdorIdx=0, mstTrialIdx=12)
+figure_2_alt.savefig("figures/system_response_poisson_noORN.{}".format(fileType), dpi=300, fig_size=(4.25,8))
+
+file = 'cache/GaussianCone_15-0-3-15_10sec/sim-13-27.npz'
+mstMATFile = 'matlab/model_cache/predictions/msp_classicalLabCond-0-15.odor-15.1-sp.1/Gaussian_15-0-3-15_10sec.mat'
+data = np.load(file)['data'][()]
+figure_3 = figure1(data, t_max=10, orn_range=-1, pn_range=[0,35], cmap='seismic',mstMatFile=mstMATFile, mstOdorIdx=15, mstTrialIdx=13)
+figure_3.savefig("figures/system_response_gaussian.{}".format(fileType), dpi=300)
+
+figure_3_alt = figure1(data, t_max=10, orn_range=-1, pn_range=[0,35], cmap='seismic')
+figure_3_alt.savefig("figures/system_response_gaussian_noMST.{}".format(fileType), dpi=300)

make_sequences.sh

@@ -0,0 +1,148 @@
+#!/bin/bash
+PYTHON_BIN="python3.6"
+N_CPU=12
+
+function defaultSparsity {
+# 2% sparsity - default condition
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name 'GaussianConnectivityHighSparsityAPL_15-0-3-15_10sec' \
+-N 50 --odor_ids 0 --odor_ids 3 --odor_ids 15 --n_cpu $N_CPU -T 10 --stimulus_dt 5 --bg_noise_scale 0.0055 --pulse_rate 14 \
+--gaussian 1 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+-o data/GaussianConnectivityHighSparsityAPL_15-0-3-15_10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityHighSparsityAPL_0-3-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+-o data/PoisonPulseConnectivityHighSparsityAPL_0-3-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityHighSparsityAPL_0-3-5-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+-o data/PoisonPulseConnectivityHighSparsityAPL_0-3-5-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityHighSparsityAPL_0-3-5-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+-o data/PoisonPulseConnectivityHighSparsityAPL_0-3-5-8-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityHighSparsityAPL_0-3-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+-o data/PoisonPulseConnectivityHighSparsityAPL_0-3-8-15-10sec.mat
+}
+
+function connectivityMediumSparsity {
+# 5% sparsity
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name 'GaussianConnectivityMediumSparsityAPL_15-0-3-15_10sec' \
+-N 50 --odor_ids 0 --odor_ids 3 --odor_ids 15 --n_cpu $N_CPU -T 10 --stimulus_dt 5 --bg_noise_scale 0.0055 --pulse_rate 14 \
+--gaussian 1 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=8.1 \
+-o data/GaussianConnectivityMediumSparsityAPL_15-0-3-15_10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityMediumSparsityAPL_0-3-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=8.1 \
+-o data/PoisonPulseConnectivityMediumSparsityAPL_0-3-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityMediumSparsityAPL_0-3-5-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=8.1 \
+-o data/PoisonPulseConnectivityMediumSparsityAPL_0-3-5-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityMediumSparsityAPL_0-3-5-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=8.1 \
+-o data/PoisonPulseConnectivityMediumSparsityAPL_0-3-5-8-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityMediumSparsityAPL_0-3-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=8.1 \
+-o data/PoisonPulseConnectivityMediumSparsityAPL_0-3-8-15-10sec.mat
+}
+
+
+function connectivityLowSparsity {
+# ~10% sparsity
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name 'GaussianConnectivityLowSparsityAPL_15-0-3-15_10sec' \
+-N 50 --odor_ids 0 --odor_ids 3 --odor_ids 15 --n_cpu $N_CPU -T 10 --stimulus_dt 5 --bg_noise_scale 0.0055 --pulse_rate 14 \
+--gaussian 1 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=12 \
+-o data/GaussianConnectivityLowSparsityAPL_15-0-3-15_10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityLowSparsityAPL_0-3-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=12 \
+-o data/PoisonPulseConnectivityLowSparsityAPL_0-3-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityLowSparsityAPL_0-3-5-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=12 \
+-o data/PoisonPulseConnectivityLowSparsityAPL_0-3-5-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityLowSparsityAPL_0-3-5-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=12 \
+-o data/PoisonPulseConnectivityLowSparsityAPL_0-3-5-8-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseConnectivityLowSparsityAPL_0-3-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams PNperKC=12 \
+-o data/PoisonPulseConnectivityLowSparsityAPL_0-3-8-15-10sec.mat
+}
+
+
+function weightMediumSparsity {
+# 5% sparsity
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name 'GaussianWeightMediumSparsityAPL_15-0-3-15_10sec' \
+-N 50 --odor_ids 0 --odor_ids 3 --odor_ids 15 --n_cpu $N_CPU -T 10 --stimulus_dt 5 --bg_noise_scale 0.0055 --pulse_rate 14 \
+--gaussian 1 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams wPNKC=20 \
+-o data/GaussianWeightMediumSparsityAPL_15-0-3-15_10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseWeightMediumSparsityAPL_0-3-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams wPNKC=20 \
+-o data/PoisonPulseWeightMediumSparsityAPL_0-3-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseWeightMediumSparsityAPL_0-3-5-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams wPNKC=20 \
+-o data/PoisonPulseWeightMediumSparsityAPL_0-3-5-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseWeightMediumSparsityAPL_0-3-5-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 5 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams wPNKC=20 \
+-o data/PoisonPulseWeightMediumSparsityAPL_0-3-5-8-15-10sec.mat
+
+$PYTHON_BIN mkDataSet_DrosoArtificialStim.py --name PoisonPulseWeightMediumSparsityAPL_0-3-8-15-10sec \
+-N 50 --n_cpu $N_CPU --pulse_rate 8 -T 10 --stim_noise_scale 0.004 --bg_noise_scale 0.0055 --min_pulse_duration 0.001 --max_pulse_duration 0.2 \
+--odor_ids 0 --odor_ids 3 --odor_ids 8 --odor_ids 15 \
+--gaussian 0 --gauss_mean 5 --gauss_std 1.5 --gauss_primary_odor_id 15 --gauss_rate_other 5 \
+--modelParams wPNKC=20 \
+-o data/PoisonPulseWeightMediumSparsityAPL_0-3-8-15-10sec.mat
+}

matlab/MSPTempotron.m

@@ -0,0 +1,91 @@
+% MSPTEMPOTRON(ts, t_i, w, V_thresh, V_rest, tau_m, tau_s) - multi-spike tempotron neuron model
+%   ts: time vector
+%   t_i: input pattern as cell array of spike times for each synapse
+%   w: synaptic efficiencies / weights
+%   V_thresh: spike threshold potential
+%   V_rest: resting potential
+%   tau_m: membrane time constant
+%   tau_s: synapse time constant
+
+function [v_t,t_out,t_out_idx,v_unreset,V_thresh, V_rest, V_0, tau_m, tau_s] = MSPTempotron(exp_fn, ts, t_i, w, V_thresh, V_rest, tau_m, tau_s)
+
+   % MATLAB vararg parsing boilerplate
+   if nargin < 6
+      tau_m = 0.020;
+   end
+
+   if nargin < 7
+      tau_s = 0.005; 
+   end
+
+
+   eta = tau_m/tau_s;
+   V_0 = eta^(eta/(eta-1)) / (eta - 1);     % normalizing constant for syn. currents
+
+   v_t = zeros(1, length(ts)) + V_rest;     % init V(t) with resting membrane potential 
+   t_out= [];                               % output spike times
+   t_out_idx = [];                          % indices of ts vector where output spikes occour
+   t_sp_idx = 1;
+
+   v_t(t_sp_idx:end) = (0 .* v_t(t_sp_idx:end));    % membrane potential V(t)
+
+   % simulate neuron
+   for i=1:length(w)
+       v_sub = msp_tempotron_kernel(exp_fn, ts, t_i{i}, tau_m, tau_s, V_0);
+       v_t = v_t + (w(i).*v_sub);    
+   end
+
+   v_unreset = v_t; % save unresetted membrane potential
+
+   % determine output spike times & perform soft-reset of V(t)
+   while (~isempty(t_sp_idx))
+       % reached V_threshold, soft-reset & emit spike time
+        above = v_t > V_thresh; % 1 and 0 for event / non event
+        crossings = diff(above);
+        idx = find(crossings>0)+0;
+        t_sp_idx = idx(idx ~= t_sp_idx);
+
+        if (~isempty(t_sp_idx))
+            t_sp_idx = t_sp_idx(1);
+            t_out = unique([t_out ts(t_sp_idx)]);
+            t_out_idx = unique([t_out_idx t_sp_idx]);
+            v_reset = (V_thresh .* exp(-(ts(t_sp_idx:end)- ts(t_sp_idx))/tau_m));% + (v_t(t_sp_idx) - V_thresh);
+            v_t(t_sp_idx:end) = v_t(t_sp_idx:end) - v_reset;
+        end
+   end
+end
+
+
+% kernel to compute synaptic input current
+function [v] = msp_tempotron_kernel(exp_fn, t, t_i, tau_m, tau_s, V_0)
+    exp_kernel = @(x, tau) exp(-(x)/tau);
+    step_fn = @(x) (sign(x) + 1) / 2;
+    v = zeros(1, length(t));
+    for i=1:length(t_i)
+       tmp = step_fn(t-t_i(i)) .* (t-t_i(i));
+       v = v + (exp_kernel(tmp, tau_m) - exp_kernel(tmp, tau_s));
+    end
+
+    v = v .* V_0;
+end
+
+function [v] = msp_tempotron_kernel_opt(exp_fn, t, t_i, tau_m, tau_s, V_0)
+    exp_kernel = @(x, tau) exp(-(x)/tau);
+    step_fn = @(x) (sign(x) + 1) / 2;
+    v = zeros(1, length(t));
+    if ~isempty(t_i)
+        tmp = zeros(length(t_i), length(t));
+        for i=1:length(t_i)
+           %tmp = heaviside(t-t_i(i)) .* (t-t_i(i));
+           tmp(i, :) = step_fn(t-t_i(i)) .* (t-t_i(i));
+        end
+
+        %tmp_idx = find(tmp > 0);
+        %tmp_2 = tmp;
+        %tmp_2(tmp_idx) = exp_kernel(tmp_2(tmp_idx), tau_s);
+        %tmp(tmp_idx) = exp_kernel(tmp(tmp_idx), tau_m);
+        v = v + sum(exp_kernel(tmp, tau_m) - exp_kernel(tmp, tau_s));
+        %v = v + sum(tmp - tmp_2);
+    end
+    v = v .* V_0;
+end