@thesamovar

Thansk for your submission, I'll assign an editor soon.

In fact, I can editing it myself unless @otizonaizit, @benoit-girard o @oliviaguest are willing to do it.
(Just add a thumb down to let me edit it)

I could edit it, especially if the final publication steps are automatized...

@benoit-girard They are not (yet). We need to rewrite both the author and editor instruction, get them approved by everybody and then things should be much more simple.

Talking of author instructions, we found it rather difficult to get the markdown->pdf compilation working. We wrote a Dockerfile and a couple of scripts that automate the process without having to install anything locally (other than docker). We'd be very happy if you wanted to make that available to other authors. Check out the readme, Dockerfile and .sh files here: https://github.com/pamelahathway/ReScience-submission/tree/Hathway-Goodman/article. This is confirmed to work on Linux and Mac, but not tried on Windows. It would be ideal to have a workflow that was accessible to everyone. I'm a Windows user myself and had to run it in a virtual machine.

@thesamovar That's great. We're also in the process of changing the submission process because you're not the only one to have had problems with md/pandoc. It will be only latex in a few days hopefully. Authors should be able to use overleaf if needed.

I'll edit you submission, sorry for the delay.

Thanks !

@damiendr @vitay You can start the review! Let's target end of June but if you need more time it's fine.

I've installed the code, just two issues to report for now:

$ python main.py
  File "main.py", line 52
    print('%s Preparing Figure 5' % time.strftime('%H:%M'))
TabError: inconsistent use of tabs and spaces in indentation

On our cluster, matplotlib raises RuntimeError: Invalid DISPLAY variable. I fixed this by adding a matplotlibrc file in code/ with backend=agg. Not sure what's ReScience's policy concerning these dependencies. [EDIT: never mind, it's probably up to me to set my global matplotlibrc correctly — found that the error also happens with other python installs on the cluster]

Now python main.py runs successfully. I'll have a go at running the array job next.

A few issues with the array job submission:

• I've adapted the options in run_main_cluster.sh to run on our cluster (Torque). But there's a typo where the main script is specified as cluster/run_main_cluster.py instead of cluster/main_cluster.py.

• The README mentions changing the hardcoded paths on lines 24-25, but while I've located para.npy I'm unsure what savedir should be set to. Could these be set to relative paths within the article repository by default?

In the meantime I'll go ahead using main.py --new True.

Should we reply and make changes as comments come in or wait until reviews are finished?

I won't start my review before next week, so I think it would be better if you fix the basic runtime problems already, so I do not run into them too.

Thank you for your feedback - I fixed the tab error in code/main.py and the typo in code/cluster/run_main_cluster.sh.

In code/cluster/main_cluster.py I set the paths as the relative paths (relative to the code/cluster/ folder), but this does not work on the Imperial College cluster because the relative paths change when the code is run.
Hardcoding the paths worked best for me, but this might be different depending on the queueing system used.

The section Running on a cluster in the code/README.md file now reflects the changes made and I also changed the paths in code/cluster/analysis_main_cluster.py to relative paths.

Thank you for the modifications and explanations !

In code/cluster/main_cluster.py I set the paths as the relative paths (relative to the code/cluster/ folder), but this does not work on the Imperial College cluster because the relative paths change when the code is run.
Hardcoding the paths worked best for me, but this might be different depending on the queueing system used.

Indeed queuing systems vary and it's unfeasible to provide a script for all of them. Regarding paths, our cluster has a similar behaviour; one workaround is to cd $PBS_O_WORKDIR in the submission script — maybe that also works on your system.

I've had further trouble with the compilation of the network by Brian failing on cluster nodes. This seems to be due to concurrent access by several processes. I've tried setting directory=None in set_device() so that each process could have its own temp dir, but found out that this created 2200 relatively large temporary directories and our cluster could not cope. So I've made a few modifications to move temp files to another filesystem and run 100 jobs in parallel, times 22 iterations in series. Hopefully I should have the results tomorrow.

General impressions

First of all, I want to congratulate the authors for that replication effort — this is clearly not an easy one to tackle. It's impressive that you managed to port all the relevant details from an event-based Matlab codebase to a timestep-driven Brian implementation. This is nice work and despite some issues which I will detail below I am inclined to believe that the original article has been replicated.

I also really appreciate that you went beyond just replicating the article and had a look at the importance of some implementation details. The point about the self-stabilising properties of the RNN rule is enlightening and adds an important nuance to the original work: not just any Hebbian flavour of STDP will do!

Running the replication

I performed all the experiments on our university cluster for convenience, even the serial ones. The cluster runs Torque on Scientific Linux 7.5 (Nitrogen). Installation was straightforward and the instructions based on a conda environment let me run main.py without dependency troubles, except for a little issue with gcc which can't be too old — we have version 4.8.5 by default on our cluster and this gave a RuntimeError: Project run failed. module load gcc-6.4 solved the issue.

I sometimes had to clean the Brian output directories (STDP_standalone*) between runs, either manually or with clean=True in set_device(); otherwise there were compilation failures or illegal instruction errors. I think this is because of different CPU architectures between nodes.

Running python main.py took only a few minutes as expected.

Running python main.py --new True took a lot more than 8 hours. You must have some fast machines at Imperial! Figure 5 took 18 hours and figure 6 took another 4 hours, for a total just under a day. That was on one of the older nodes of our cluster (Xeon E5-2660 @ 2.2GHz from 2012).

The parameter exploration code was tricky to get to run. There were frequent but non-deterministic compilation or run failures on cluster nodes: it seems that parallel jobs can't share the same Brian output directory. Do you know why that wasn't an issue on your cluster? Could multiple runs potentially use the same standalone code, or does it have to be recompiled for each parameter combination?

My first solution was to set directory=None so that each job gets its own output temp dir. However, these get large (~3GB with static_arrays and results), and they are not deleted automatically. For the array job with 2200 conditions, that would require a total of 6TB and our cluster doesn't have that much space available on /tmp, as I found out the hard way ;-).

In the end here's what worked for me:

• export TMPDIR=/beegfs/username/tmp/. This moves python's tempfiles to our distributed storage and avoids leaving logs etc. in each node's /tmp.
• create a separate output directory $TMPDIR/$PBS_JOBID/$PBS_ARRAYID that I passed as an argument to set_device() in main_cluster.py and deleted at the end of each run, to prevent results from accumulating.
• reduce the parallelism by doing 100 jobs in parallel, times 22 iterations in series. Our cluster monitoring software had a fainting fit with the 2200-long array job attempt, though our admin assures me it recovers just fine after a while.

With those settings the array job took 3 to 4 hours to complete.

Here is my submission script for the record (requires a few changes to main_cluster.py to handle the iteration variable and output directory).

Code

The replication is implemented in Python with few dependencies: it uses Brian for the network simulation and Numpy for creating the input patterns. The authors used Numba to speed up some functions and Brian's standalone mode to compile the network. One advantage of Brian, for a replication, is that all model equations are clearly visible instead of being hidden in a library of standard models.

I found the code generally readable, though it is complex in some parts, especially the manipulation of spike trains in create_input.py; here a few more comments describing the algorithm implemented in each function would be welcome. The model is defined repeatedly in main_cluster.py, run_simulation.py and figure_6.py, which is not ideal, though I understand that there are slight variations each time.

Random seeding issue

The results are not deterministic: some figures are different after every run. My guess is that it's because the code sets brian2.seed(x) but not numpy.random.seed(x). On top of that it seems that the numpy seed should be set again inside a Numba function, because Numba has its own RNG.

Thus the data generated in python (ie. the input spike trains) does not use a fixed seed at the moment. This could be an issue for someone trying to validate their own replication, and also because some of your points involve picking a particular run that displays certain characteristics, and that selection cannot be reproduced reliably at the moment.

Figures

Here I will be commenting on the replication of the original figures as well as on the reproduction of the figures of the replication. ;-)

Overall, the generated figures are clear, well labeled and aesthetically pleasing.

Figure 1 is well replicated.

Figure 2 looks fine except it is missing some annotations — are these added by hand? Not a big deal though.

In the caption, the dashes around "in this example" look a bit short (use --- in pandoc's markdown to make a —)

For figures 3A and 4A, I understand that you were able to re-run the original matlab code with your own inputs and pass the results to your plotting routines, hence the exact same spike patterns in 4A and 4B. Correct me if I'm wrong.

Figure 3B looks fine, but shows wide variations between runs (seeding issue). Rarely, it entirely fails to find the pattern (expected).

Figure 4 looks fine appart from the seeding issues. Minor issue: some of the white circles are almost entirely obscured by dark circles, thus the apparent density may not reflect the data — maybe some transparency would give better results.

I found the caption to be quite obscure at first: the plot looks like a spike raster, but it says it's a plot of the weights. Only after reading the original paper did I understand that this was a plot of input spikes, coloured by the corresponding weight. The caption could be amended to make this clearer.

Figure 5 can be created in two ways: fig_5_new() and main_cluster.py.

The restricted version created by fig_5_new() is missing the dt=1e-6 condition. Otherwise I could find no obvious discrepancies given the limited number of runs. There are a few cosmetic problems.

The full version created by the cluster has some issues. The legend is missing the 'orig paper' label. The plots also show only one of the two conditions: dt=1e-4 or 1e-6. But while the legend says 1e-4, the plain black lines seem to correspond better to the dotted lines in the replication paper, which labels them 1e-6. Looking into the logs from the cluster run, I find dt = 100. us.

Figure 6 is missing its markers, and the variance seems to be half as large as reported for the 1e-4 condition. But the figure still validates the point about the effect of the timestep, with low variance and the expected mean at dt=1e-6.

I'm unable to reproduce figure 7 reliably. For figures 7C and 7D, this seems to be due to the seeding issue. Most of the runs are a lot worse than the one shown in the article — which sort of proves your point, actually!

7B sometimes looks like the one you show, and sometimes RNN converges to a value > 0.5 instead, as shown below.

7A looks very different: the curves are always straights lines. There must be another problem here.

Figure 8 is well reproduced.

Model description

The authors chose not to repeat the model equations that didn't differ from Masquelier et al (2008), for instance those for the STDP rule. I am personally fine with this.

Case of Δt=0

The replication says:

If the learning rule is modified so that the synaptic weight is increased as well as depressed when input neuron and output neuron spikes happen during the same time step (effectively causing a smaller weight increase), the pattern is found earlier

The variant seems to be implemented in main_cluster.py line 306:

stdp_on_post = ('''LTDtrace = -aminus
               wi = clip(wi + LTPtrace + 0.5*(aminus-aplus)*(lastspike_pre==lastspike_post), wmin, wmax)
               LTPtrace = 0''')

Therefore it will be used to generate figure 5 when run on the cluster.

But the text isn't very clear about the scope of this variant, as it first mentions "the rest of the paper" and only then figure 5 in passing. Maybe you could be more explicit about using the variant to generate figure 5 specifically? Or should the variant have been commented out for generating figure 5 too?

Regarding the variant code, I'm not sure why you take the mean of LTP and LTD in this case, rather than the just sum?

Other results from the original paper

While I am reasonably confident that the replication does reproduce the original work, adding a replication of fig. 1 (original) would confirm the equivalence of the input spike trains.

Another figure wich could be nice to replicate is fig. 4 (original), showing the output firing patterns over time. Though in this case the equivalence is implied to some degree by your fig. 3 (fig. 5 original) which shows the development of selective responses.

Other remarks

The following remarks are things that popped into my mind while reading the work — not comments to address for the review!

I'm curious about the fact that fig. 4 in the original paper shows the membrane potential relaxing to baseline over the exact duration of the pattern. I guess this is because the later spikes moved onto the LTD side of the output spike earlier in the run, and were therefore potentiated less?

I also wonder about the importance of the after-spike hyperpolarisation for learning. Is it critical for stabilising the learning, does it have to match the duration of the pattern or the time constants of the EPSPs, and would the model still work with a shorter or longer period?

Thanks @damiendr for the excellent and detailed review. We'll reply in detail later, but for now just wanted to make a quick comment about Brian standalone mode.

Your problems with the generated code directories of Brian are unavoidable with the current release version of Brian unfortunately, but I've created an issue to sort this out for the next release (brian-team/brian2#973). In the meantime we could sort it out by explicitly creating a temporary directory (https://docs.python.org/2/library/tempfile.html#tempfile.mkdtemp) and then deleting the directory after we've got the data back (more or less what you did). In the release version of Brian the code needs to be recompiled for each parameter change (something I'm also working on fixing in Brian at the moment). By default, we use -march=native compilation flag, so it will also need recompilation on each processor architecture.

Thank you for taking the time and effort to get the parameter exploration to run and for providing your submission script, @damiendr!

As a response to the small issues you reported, I made a few changes to the figures. All other comments will be addressed after the second review.

Running the replication

Dan already commented on the Brian standalone mode and how to manage the directories until the next release version of Brian.

Our cluster runs gcc version 5.4.0 by default, it's good to know that that older versions might create an issue. I added a line into the .sh file to load gcc 5.4.0.

I only tried main.py --new True on my work computer and laptop and not on the cluster, so I wasn't aware that it took so long on older nodes, I will make a note of that in the README.

Figures

Figure 5:

I fixed the cosmetic problems in figure_5_new() and the labeling in main_cluster.py.

The reason the 1e-6 condition is missing in both fig_5_new() and the cluster version is that a time sep of 0.001 ms (1e-6) instead of 0.1 ms (1e-4) takes about 100 times longer to complete. Since the results from the two time steps is not that large, we thought it might be sufficient to run the larger time steps (0.1 ms).

It is of course possible to run them at 0.001 ms as well, the time step can be changed in line 28 of run_simulation.py and line 33 of main_cluster.py.

Figure 6:

I fixed the missing markers.

Figure 7A:

The plotting resolution was wrong, this has now been corrected. The plotting resolution for these two figures is slightly different, but the general increases/decreases should be clearly visible.

Figure 7B:

I cannot reproduce the increase of weights for RNN that you show. Did this occur often?

Case of Δt=0

You are right, the normal STDP rule was commented out and the modified rule was commented in - I changed this now. This might also be the reason why the performance was a bit better at pattern frequency 0.5 and looks more like the dotted line in the article figure (dt=1e-6 s). While this was a mistake, it shows that this modification only has an impact in very specific circumstances.

Response to Reviewer 1

Below we comment on the review by @damiendr and the changes made to the article and the code in response.

Running the replication

I performed all the experiments on our university cluster for convenience, even the serial ones. The cluster runs Torque on Scientific Linux 7.5 (Nitrogen). Installation was straightforward and the instructions based on a conda environment let me run main.py without dependency troubles, except for a little issue with gcc which can't be too old — we have version 4.8.5 by default on our cluster and this gave a RuntimeError: Project run failed. module load gcc-6.4 solved the issue.

The README now states that gcc version 5.4.0 or higher is required.

Running python main.py --new True took a lot more than 8 hours. You must have some fast machines at Imperial! Figure 5 took 18 hours and figure 6 took another 4 hours, for a total just under a day. That was on one of the older nodes of our cluster (Xeon E5-2660 @ 2.2GHz from 2012).

A note was made in the README file that runs might take longer.

I sometimes had to clean the Brian output directories (STDP_standalone*) between runs, either manually or with clean=True in set_device(); otherwise there were compilation failures or illegal instruction errors. I think this is because of different CPU architectures between nodes.

[...]

Running large array jobs on a cluster can be difficult to get to run. Issues are difficult to foresee and to debug, so for the moment the code remains as it is, but thank you very much for your solutions and notes - they will definitely be helpful for someone else with a similar cluster!

Code

I found the code generally readable, though it is complex in some parts, especially the manipulation of spike trains in create_input.py; here a few more comments describing the algorithm implemented in each function would be welcome. The model is defined repeatedly in main_cluster.py, run_simulation.py and figure_6.py, which is not ideal, though I understand that there are slight variations each time.

Comments were added to the functions in create_input.py. Most of them were adapted from the Matlab code of the original article.

Random seeding issue

The results are not deterministic: some figures are different after every run. My guess is that it's because the code sets brian2.seed(x) but not numpy.random.seed(x). On top of that it seems that the numpy seed should be set again inside a Numba function, because Numba has its own RNG.

Thus the data generated in python (ie. the input spike trains) does not use a fixed seed at the moment. This could be an issue for someone trying to validate their own replication, and also because some of your points involve picking a particular run that displays certain characteristics, and that selection cannot be reproduced reliably at the moment.

Thank you for pointing this out, the random seed issue is now corrected - each run with a certain seed is now deterministic. Figure 7CD now shows a similar behaviour to the one in the article with seed 28 as an example.

Figures

Figure 2 looks fine except it is missing some annotations — are these added by hand? Not a big deal though.

Yes, the annotations are done by hand.

In the caption, the dashes around "in this example" look a bit short (use --- in pandoc's markdown to make a —)

This has been corrected.

For figures 3A and 4A, I understand that you were able to re-run the original matlab code with your own inputs and pass the results to your plotting routines, hence the exact same spike patterns in 4A and 4B. Correct me if I'm wrong.

We used the input spikes created by the Matlab code and ran it both in the original code and in ours, therefore the spike patterns are the same.

Figure 3B looks fine, but shows wide variations between runs (seeding issue). Rarely, it entirely fails to find the pattern (expected).

Indeed, with every seed, Figure 3B will look slightly different.

Figure 4 looks fine appart from the seeding issues. Minor issue: some of the white circles are almost entirely obscured by dark circles, thus the apparent density may not reflect the data — maybe some transparency would give better results.

We acknowledge that it is not immediately apparent that this figure represents a raster plot due to the high density of spikes and the overlaps. Since we were replicating the figure in the original paper (Figure 6a), we did not want to change the format of the plot. A colorbar was added to make the colouring of the dots clearer.

I found the caption to be quite obscure at first: the plot looks like a spike raster, but it says it's a plot of the weights. Only after reading the original paper did I understand that this was a plot of input spikes, coloured by the corresponding weight. The caption could be amended to make this clearer.

The caption was changed to make it clearer that it shows a raster plot of spikes.

Figure 5 can be created in two ways: fig_5_new() and main_cluster.py.

The restricted version created by fig_5_new() is missing the dt=1e-6 condition. Otherwise I could find no obvious discrepancies given the limited number of runs. There are a few cosmetic problems.

The full version created by the cluster has some issues. The legend is missing the 'orig paper' label. The plots also show only one of the two conditions: dt=1e-4 or 1e-6. But while the legend says 1e-4, the plain black lines seem to correspond better to the dotted lines in the replication paper, which labels them 1e-6. Looking into the logs from the cluster run, I find dt = 100. us.

The cosmetic problems in figure_5_new() and the labeling in main_cluster.py were corrected.

The reason the 1e-6 condition is missing in both fig_5_new() and the cluster version is that a time sep of 0.001 ms (1e-6) instead of 0.1 ms (1e-4) takes about 100 times longer to complete. Since the results from the two time steps is not that large, we thought it might be sufficient to run the larger time steps (0.1 ms).

Figure 6 is missing its markers, and the variance seems to be half as large as reported for the 1e-4 condition. But the figure still validates the point about the effect of the timestep, with low variance and the expected mean at dt=1e-6.

Markers were added. The variance might be lower due to the random seeds not being fixed (runs were not the same), but the point is still the same.

I'm unable to reproduce figure 7 reliably. For figures 7C and 7D, this seems to be due to the seeding issue. Most of the runs are a lot worse than the one shown in the article — which sort of proves your point, actually!

This has now been addressed and Figure 7CD now shows a behaviour similar to the one shown in the article (seed=28 in this example, which is the default seed for Figure 7CD).

7B sometimes looks like the one you show, and sometimes RNN converges to a value > 0.5 instead, as shown below.

If this behaviour only occurred once, it might have been one of the rare failed runs, for example now seed=0 will fail and show this behaviour.

7A looks very different: the curves are always straights lines. There must be another problem here.

The plotting resolution was wrong, this has now been corrected. The plotting resolution for these two figures is slightly different in the article figure, but the general increases/decreases should be clearly visible.

Model description

The authors chose not to repeat the model equations that didn't differ from Masquelier et al (2008), for instance those for the STDP rule. I am personally fine with this.

We will leave it like it is unless Reviewer 2 has a different opinion.

Case of Δt=0

The replication says:

If the learning rule is modified so that the synaptic weight is increased as well as depressed when input neuron and output neuron spikes happen during the same time step (effectively causing a smaller weight increase), the pattern is found earlier

The variant seems to be implemented in main_cluster.py line 306:

stdp_on_post = ('''LTDtrace = -aminus
               wi = clip(wi + LTPtrace + 0.5*(aminus-aplus)*(lastspike_pre==lastspike_post), wmin, wmax)
               LTPtrace = 0''')

Therefore it will be used to generate figure 5 when run on the cluster.

But the text isn't very clear about the scope of this variant, as it first mentions "the rest of the paper" and only then figure 5 in passing. Maybe you could be more explicit about using the variant to generate figure 5 specifically? Or should the variant have been commented out for generating figure 5 too?

This was a mistake, it should have run the standard STDP equations, this was changed in the main_cluster.py script.

Using the modified STDP rule might be the reason why the performance was a bit better at pattern frequency 0.5 and therefore looks more like the dotted line in the article figure (which corresponds to dt=1e-6 s). Apart from that condition the dotted and the solid lines are very similar.

The modification only comes into play at an output neuron spike, when it falls into the same time window as some input neuron spikes. In the pattern frequency = 0.5 condition there are more instances of the pattern and therefore (if successful) more output neuron spikes. It is therefore interesting that modifying the STDP rule for the case of Δt=0 increases the performance. While for Figure 5 the standard STDP rule should have been run (and was run for the article figures), it shows that this modification only has an impact in very specific circumstances.

Regarding the variant code, I'm not sure why you take the mean of LTP and LTD in this case, rather than the just sum?

The article text was changed to explain the modified learning rule better.

In essence, the modification accounts for the fact that input neuron and output neuron spikes fall into the same time window in a time-based simulation, which would not happen in an event-based simulation. If there is one input neuron spike and one output neuron spike in the same time window, there is a 50% chance of the input neuron spike being before the output neuron spike and vice versa. The modification reflects this and adds the mean of LTP and LTD in those cases.

Other results from the original paper

While I am reasonably confident that the replication does reproduce the original work, adding a replication of fig. 1 (original) would confirm the equivalence of the input spike trains.

Another figure wich could be nice to replicate is fig. 4 (original), showing the output firing patterns over time. Though in this case the equivalence is implied to some degree by your fig. 3 (fig. 5 original) which shows the development of selective responses.

main.py and run_simulation.py now also create Figure 9 (original Figure 4) and Figure 10 (original Figure 1) as supplemental figures (they have not been added to the article). Note that Figure 9 only looks like the original figure if the pattern finding was successful and that the two spikes per pattern in original Figure 4b only rarely occur.

Other remarks

The following remarks are things that popped into my mind while reading the work — not comments to address for the review!

I'm curious about the fact that fig. 4 in the original paper shows the membrane potential relaxing to baseline over the exact duration of the pattern. I guess this is because the later spikes moved onto the LTD side of the output spike earlier in the run, and were therefore potentiated less?

Yes this is correct, the spikes just after the output spike are mostly depressed, therefore the voltage is slow to reach levels close to threshold again during the pattern.

I also wonder about the importance of the after-spike hyperpolarisation for learning. Is it critical for stabilising the learning, does it have to match the duration of the pattern or the time constants of the EPSPs, and would the model still work with a shorter or longer period?

Yes, there seems to be a minimum time between spikes necessary for the model to work.

Increasing the negative afterpotential to create longer periods between spikes works just as well despite lower output neuron firing rates that result from the large afterpotential. Decreasing the size of the afterpotential leads to a faster increase in voltage after a spike and therefore higher firing rates, which at some point leads to a failure to find the pattern (very high firing rates lead to an overall depression of all weights and a failure to become specific the pattern). We did not attempt to quantify this (how small the afterpotential can be) but it would be interesting to see whether there it is a relationship with parameters such as the time constants of the EPSPs.

Sorry for the delayed review. The first review was already quite extensive and I agree with all the points (which are already solved anyway), so I will just make a couple of additional remarks.

The submission is of extremely high quality, be it the correctness/exhaustiveness of the reimplementation, the description of the reimplementation, or the added value of the additional experiments performed on the original model. What I particularly appreciate is the understanding of the conditions and mechanisms allowing the Masquelier model to work or not. I tried last year to reproduce that same model using another neurosimulator (ANNarchy) and failed miserably (it detected patterns, but there were a lot of false positive, i.e. spikes between two patterns. There was not enough depression). Now I understand why: I took dt=0.5 ms, wrongly understood what they meant with nearest spike approximation (I implemented what you call NN, not RNN, so it ended up with too many potentiations) and did not care about simultaneous pre and post spikes. It had no chance to work... Thanks to your work I will try again.

Running the replication

I ran the preloaded simulations on an arch system without anaconda (gcc=8.1.1 python=3.6.5 numpy=1.14.5 matplotlib=2.2.2 numba=0.38.0 brian2=2.1.3.1) and it ran without problems. It is good to know that your code is not bound to the exact versions mentioned in environment.yml and still usable with newer versions. Perhaps you could use >= instead of = in that file and the README for users from the future? I am not sure how long anaconda ships specific package versions.

I ran the complete simulation on a Mint desktop using anaconda and the suggested package versions (gcc=5.4.0). I do not have a cluster at disposal so I could not try that (or would it work on a single machine?). As the first reviewer managed to do it, I guess it is OK. The simulation took around eight hours as expected. It was not a high-end machine either (Xeon X5675, 12 GB RAM, 2011) so I do not know what happened so it lasts 24 hours for @damiendr. The simulation was single-threaded all along.

Code

The code is clear and well structured. It could have deserved a bit more comments, but for who is already familiar with brian, it should not be a problem to understand what it is doing.

Article

The paper is very clear and well written. The results from Masquelier et al. 2008 were successfully reproduced (within a very reasonable margin due to random initializations and a fundamentally different implementation - SRM vs. LIF).

I am fine with not repeating the equations from Masquelier which have not changed, but I find Eqs 1,2,3 confusing. $x$ is the EPSP for a single synapse, but it is not clear how multiple presynaptic spikes influence the postsynaptic potential. A sum over afferent spikes is missing. Also, the synaptic weights $w$ should intervene somewhere: either $x$ is incremented from $w$ instead of one, or it is multiplied by it in the equation for $u$.

One should also specify the numerical method used to solve the ODEs: Euler, exact, exponential? It plays surely a role regarding the right values for dt.

It would be nice to add some comments on how to go from a SRM model to a discretized ODE-based model. Perhaps recap the idea of SRM, and that the kernels they use are just alpha functions which can be very well modeled by two coupled ODEs. In SRM, they cut the kernels at +/- 7*tau. Is something similar necessary with ODEs, or are the variables already so close to zero that it does not matter?

Effect of learning rule at $\Delta t = 0$: "it differs from the implementation used in the majority of modelling papers". In my experience (I may be wrong), papers mostly ignore that problem and rely on the choices made by the underlying simulator. Brian treats them as pre-post (but because of RNN, no? with the "classical" NN, brian would treat them as both pre-post and post-pre). ANNarchy as both a pre-post and a post-pre, unless a flag unless_post is set to have only post-pre (if the spikes occur at the same time, there is no causality between pre and post firings - pre did not make post fire -, so the weight should not be increased). I do not know what NEST does. You do not have to discuss this in the article, this sentence just made me think.

Typos: (version 5de5d9a)

• LIF neuron: We converted the kernels of the posystnaptic potential
• Implementation details, version of STDP learning rule: The choice of THE learning rule (or? I am not a native English speaker...)

The new figures 9 and 10 are nice, do you plan to ship a supplementary file too, or do they only appear when running the code?

@damiendr @vitay Thanks for your very detailed review. I think @damiendr comments have already been adressed. @pamelahathway @thesamovar Can you address the remaining comments and questions?

After this, @damiendr and @vitay, could you tell me if you accept the submission for publication?

Response to Reviewer 2

Thank you for your insightful comments, @vitay.

Here are our answers to the second review and the code and article with the incorporated changes.

Sorry for the delayed review. The first review was already quite extensive and I agree with all the points (which are already solved anyway), so I will just make a couple of additional remarks.

The submission is of extremely high quality, be it the correctness/exhaustiveness of the reimplementation, the description of the reimplementation, or the added value of the additional experiments performed on the original model. What I particularly appreciate is the understanding of the conditions and mechanisms allowing the Masquelier model to work or not. I tried last year to reproduce that same model using another neurosimulator (ANNarchy) and failed miserably (it detected patterns, but there were a lot of false positive, i.e. spikes between two patterns. There was not enough depression). Now I understand why: I took dt=0.5 ms, wrongly understood what they meant with nearest spike approximation (I implemented what you call NN, not RNN, so it ended up with too many potentiations) and did not care about simultaneous pre and post spikes. It had no chance to work... Thanks to your work I will try again.

These were also errors we made when we started, the original paper is not very clear on the learning rule and does not use discrete time steps and therefore does not comment on this.

Running the replication

I ran the preloaded simulations on an arch system without anaconda (gcc=8.1.1 python=3.6.5 numpy=1.14.5 matplotlib=2.2.2 numba=0.38.0 brian2=2.1.3.1) and it ran without problems. It is good to know that your code is not bound to the exact versions mentioned in environment.yml and still usable with newer versions. Perhaps you could use >= instead of = in that file and the README for users from the future? I am not sure how long anaconda ships specific package versions.

Since we are not completely sure that future versions will work, we will leave the = , but a few lines were added to the README explaining that later versions probably also work.

I ran the complete simulation on a Mint desktop using anaconda and the suggested package versions (gcc=5.4.0). I do not have a cluster at disposal so I could not try that (or would it work on a single machine?). As the first reviewer managed to do it, I guess it is OK. The simulation took around eight hours as expected. It was not a high-end machine either (Xeon X5675, 12 GB RAM, 2011) so I do not know what happened so it lasts 24 hours for @damiendr. The simulation was single-threaded all along.

It would work on a single machine, but would take days if not weeks, so it definitely is best run in parallel on a cluster.

Code

The code is clear and well structured. It could have deserved a bit more comments, but for who is already familiar with brian, it should not be a problem to understand what it is doing.

More comments were added to the Brian code.

Article

The paper is very clear and well written. The results from Masquelier et al. 2008 were successfully reproduced (within a very reasonable margin due to random initializations and a fundamentally different implementation - SRM vs. LIF).

I am fine with not repeating the equations from Masquelier which have not changed, but I find Eqs 1,2,3 confusing. $x$ is the EPSP for a single synapse, but it is not clear how multiple presynaptic spikes influence the postsynaptic potential. A sum over afferent spikes is missing. Also, the synaptic weights $w$ should intervene somewhere: either $x$ is incremented from $w$ instead of one, or it is multiplied by it in the equation for $u$.

It is true that we missed adding the weights to $x$, which is now corrected: In the event of a presynaptic spike at neuron $i$, $x$ is increased by $w_{i}$ ($x\leftarrow x+w_{i}$).

In Brian, x is a variable of the postsynaptic neuron (the equations eqs in the code are used for the postsynaptic neuron N2), which represents the sum of all incoming EPSPs. Each incoming spike from synapse $i$ increases the postsynaptic variable $x$ by its weight $w_{i}$ via the synapse equation stdp_on_pre = '''x_post += deltax*wi''' (deltax=1), which is called at every presynaptic spike. Since $x$ is incremented by $w_{i}$ and not set to a certain value, multiple $w_{i}$ can be added to the existing trace when multiple spikes happen in succession from the same synapse.

One should also specify the numerical method used to solve the ODEs: Euler, exact, exponential? It plays surely a role regarding the right values for dt.

Since ODEs in this case there is an exact solution to the equations, we used the exact (in the code it is called 'linear') method for integrating the differential equations and added this to the text.

We ran 10 repetitions or each 'linear', 'euler', 'exponential_euler', 'rk2' and 'rk4' integration methods. While the exact output neuron spike timings differed slightly (due to slight variations in potentiations and depressions), the overall performance was the same for all methods except for the exponential Euler method, which fails about 50% of the time. Since this method is not an appropriate choice as an integration method for the ODE's used in this model, this is not surprising.

It would be nice to add some comments on how to go from a SRM model to a discretized ODE-based model. Perhaps recap the idea of SRM, and that the kernels they use are just alpha functions which can be very well modeled by two coupled ODEs.

This was added to the text:

"Brian on the other hand uses differential equations to model the system parameters and evaluates those equations for each time step. The SRM kernels can be modelled using alpha functions [@Brette2007], so we converted the postsynaptic potential and the spike afterpotential into the following differential equations:"

There is a detailed explanation of how to go from SRM to ODE's in the Brian documentation: http://brian2.readthedocs.io/en/stable/user/converting_from_integrated_form.html

In SRM, they cut the kernels at +/- 7*tau. Is something similar necessary with ODEs, or are the variables already so close to zero that it does not matter?

Using ODE's eliminates the necessity to use such cut-offs as in the original paper. However, the variables are already 0 at +/-7*tau.

Effect of learning rule at $\Delta t = 0$: "it differs from the implementation used in the majority of modelling papers". In my experience (I may be wrong), papers mostly ignore that problem and rely on the choices made by the underlying simulator. Brian treats them as pre-post (but because of RNN, no? with the "classical" NN, brian would treat them as both pre-post and post-pre). ANNarchy as both a pre-post and a post-pre, unless a flag unless_post is set to have only post-pre (if the spikes occur at the same time, there is no causality between pre and post firings - pre did not make post fire -, so the weight should not be increased). I do not know what NEST does. You do not have to discuss this in the article, this sentence just made me think.

We removed that half sentence since most modelling papers probably do ignore this effect.

Typos: (version 5de5d9a)

• LIF neuron: We converted the kernels of the posystnaptic potential

This was corrected, thank you.

• Implementation details, version of STDP learning rule: The choice of THE learning rule (or? I am not a native English speaker...)

We think this is correct, so we will leave out the 'the'.

The new figures 9 and 10 are nice, do you plan to ship a supplementary file too, or do they only appear when running the code?

The two figures do not add much more information about the equivalence of the two models, so they will only appear when running the code.

Thanks for your detailed reply. I am satisfied with the changes and recommend the acceptance of the submission.

Thank you for your replies. I have been able to re-run the updated code on our cluster and I can confirm that the results for figure 5 now match. The fact that the modified rule caused the original discrepancy is interesting — I did not expect the Δt=0 case to make so much difference in practice. I learned something here about these corner cases of neural simulations, which is good!

The modifications for the other figures and the text also address my remarks, and my recommendation is therefore to accept the article.

Many, many thanks @vitay and @damiendr for your reviews and comments!

Congratulation @thesamovar, the article is accepted for publication. It will take a few days until it appears online. Don't hesitate to ping me if it is not the case.

Article has been published and will soon appeared at https://rescience.github.io/read/

Repository has been archived at https://github.com/ReScience-Archives/Hathway-Goodman-2018
See also https://github.com/ReScience/ReScience/blob/master/Volume%204%20-%20Issue%201.md