[Re] Mathematical analysis of depolarization block mediated by slow inactivation of fast sodium channels in midbrain dopamine neurons.

[maryupshall]

Original article: https://www.ncbi.nlm.nih.gov/pubmed/25185810

PDF URL: https://github.com/mupsh/ReScience_Qian_2014/blob/master/article/article.pdf Metadata URL: https://github.com/mupsh/ReScience_Qian_2014/blob/master/article/metadata.yaml Code URL: https://github.com/mupsh/ReScience_Qian_2014

Scientific domain: Computational Neuroscience Programming language: Python Suggested editor: Benoit Girard, Tiziano Zito

[rougier]

Thanks for your submission. I'll edit it and assign reviewers soon.

[rougier]

@rgutzen @cJarvers Could you review this submission ?

[rgutzen]

Sure. I'm happy to

[cJarvers]

I'd be happy to review this as well, but I'll be very busy throughout November. Is early December an acceptable time frame for the review? Otherwise I have to pass.

[rougier]

Thanks @rgutzen and @cJarvers. @cJarvers I think we can wait for you review in December yes. @mupsh Is that ok for you ?

[maryupshall]

Hi,

Yes, that’s ok for me.

Cheers,

Mary

On Oct 29, 2019, at 11:07 AM, Nicolas P. Rougier notifications@github.com wrote:

Thanks @rgutzen https://github.com/rgutzen and @cJarvers https://github.com/cJarvers. @cJarvers https://github.com/cJarvers I think we can wait for you review in December yes. @mupsh https://github.com/mupsh Is that ok for you ?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ReScience/submissions/issues/9?email_source=notifications&email_token=ALXX22SQIASVDPSPMWB63RLQRBGUNA5CNFSM4JFG6XY2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECQ3J2A#issuecomment-547468520, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALXX22X2CVJRBIPKWJYDB7TQRBGUNANCNFSM4JFG6XYQ.

[rgutzen]

Hi @mupsh, I think you chose a nice study to replicate and I really hope I can help improve it towards a level of rigor, I think it could have. So, here my review. You show some qualitative and some quantitative replications, but you also show some clear discrepancies that you either ignore or don't discuss sufficiently. I think there are a couple of major and minor changes that can make your replication study much stronger and more interesting and useful for the reader. In the following, I ordered my comments into general remarks, comments related to the code, comments regarding the replication of each figure, and comments related to the actual article text.

General comments

• First of all, there is only a brief introduction, which doesn't answer the interesting question of what your motivation is for replicating exactly this study? And why is this a relevant paper? What are the key insights? And why would it be beneficial for the scientific community to replicate these results?
• There is no information on the reproducibility of the original implementation. Is the original data available? Is the original code available? Only at the very end, you mention the communication with the authors. This is very relevant to follow your replication and should be moved to the beginning of the paper and elaborated on.
• Ideally, you can also ask the original authors, if they agree to publish the code they provided you alongside this replication study
• Your paper is rather short, and writing short precise papers is generally desirable. However, you mostly just state what you replicated, and I'm totally missing a description of your journey how you replicated the results. What was the approach? What are the design principles or your implementation? Which difficulties did you encounter? What did not work? These aspects, I would consider the most valuable part of your research, and the ones others can really benefit and learn from.
• Especially, when you are not able to replicate a result or figure (exactly, or at all), it is not sufficient to say that you tried, it didn't work out, and then speculate that it might be due to different parameters. There needs to be a more detailed account of what exactly you tried, why you think the replication fails, and generally, you need to present a reasonable effort in investigating and fixing the failed replication.
• For each figure I'm missing an explanation of what exactly is shown and why this is interesting. This doesn't need to be long, but it is required in order to make sense of the replication, as this is the logical basis on which anyone can then start evaluate the similarity to the original results.
• For generating the simulation results, additional simulation parameters are set. Instead of hiding them in the code, they should be reported, and if necessary, discussed (e.g time resolution, input, discarded initialization time, ...)

Code

General

• The main goal of replication studies is to validate results from other publications by reimplementation. However, another goal is also to provide understandable and easily reproducible code so that others get better insights into the computational aspects of the original study and reuse/edit/expand the analysis or part of it. The second goal requires the code to be in a certain format. If possible, it should use standards and existing software and packages which are common in the field. Furthermore, it needs to be logically structured and well documented. Writing custom code to replicate other custom code misses the opportunity of adding real value to the implementation.
• So I strongly encourage you to use standard resources for the neuron models, like for example PyNN or NeuroML
• Besides documenting the code and structuring it in a way that makes it easy to read and understand, please also explain your code in the article. This is a computational paper and addressing and discussing certain implementation decisions in the text can be very helpful. The main outline of the implementation and the design principles could even be presented in their own Section.
• For the packages you do use, make sure you also mention them in the article, state what they do and why you chose them, especially for less common packages.

  Details

• executing run.py in the provided environment return an error while producing figure 3: unable to execute 'swig': No such file or directory error: command 'swig' failed with exit status 1 After installing swig version 3.0.12 it runs through.
• For understandability, it helps to assign variables with telling names instead of using numerical values in equations without explanations.
• code lines like end_times = np.ravel(list(zip(np.arange(3, last_set + 100 + 1, 100), np.arange(100, last_set + 100 + 4, 100))))
• for some reason you named tau_hs different from the other variables as 'tah_hs'
• I strongly suggest making use of packages to handle units (like quantities or pints) at least when variables are passed to functions. This easily ensures that the right magnitudes and units are used. Once the units are check in a function, they can then be omitted for the actual computations to not slow them down. Otherwise, it can get quite bothersome to trace back which magnitude the resulting values have.
• For a better overview and less error-proneness, it is generally advisable to outsource parameter settings (into a yamls, json, or event txt file), and import them and not set them by hand in each function.
• Also the parameters list is rather fragile since it relies on an arbitrary but specific order of the parameters. A dictionary would be a better choice.

Replicating the figures

General

• Using additional ticks (and tick labels) on the axis and would make the evaluation of the figures the comparison much easier.
• Use smaller linewidth, to not cover up the smaller details.
• Explain all the elements in each figure!
• Emphasize the relevant points in the plots. Discuss which aspects of the plots are the relevant results you aim to replicate?

Figure 1a

• A legend is missing which of the curves is with and without hs.
• The code contains the following comment. # need to divide the given value by 1.5 to get the correct graph. This operation is not motivated or even mentioned in the text.
• The vertical axis misses its units
• I am missing a statement why you consider the plot replicated, like giving the agreement of characteristic points (e.g. minimum, horizontal trend in the beginning, Voltage at current drop-off, skewed parabola shape around minimum)

Figure 1b

• even when the plot is mirrored in the x-direction, there are still clear differences (absolute and relative peak sizes, the distance between peaks) in the plot which needs to be addressed!
• Without a reason why the simulation would be backward in time, and given the other differences in the plot, the resulting plot suggests that there are other discrepancies in the equations or the code which need be uncovered. In this state, I would consider the plot not replicated.
• The code contains comment # need to divide given value by 2 to get correct graph. However, this is not mentioned in the text.
• g_Na is set to 5.92 nS, although the original article reports using 9.12 nS.
• The corresponding code is quite messy and thus hard to debug. To illustrate, I try to understand how the current trace which is plotted evolves. You are plotting the clamp_current which is calculated by sodium_current(), which takes v,h,and hs as state/time-dependent arguments. These are calculated by solving the voltage_clamp ODE using ic initial condition in the time packet time. Looking into voltage_clamp(), __clamp__() is called which calls the ODE function from the parameters list (here ode_3d()). ode_3d() itself is not easy to understand and debug since it has several not intuitive equations which stretch over multiple line. But most remarkably ode_3d() does not use its time argument t, which essentially means that the whole current function does not depend on the time packet it is iterated over. The main point here is that you organize and structure the code obfuscates dependencies and makes it really hard to understand what is happening in the code.

Figure 1c

• it should be made clear that the time course of n(h), which is not a function, is fitted by a function f(h)
• "to get the same results" Figure 1C is visibly different to the original!!
• The n-h relation in the simulation varies over time. See the spread of the curve in the original plot. This is not visible in the replication and not commented.
• The characteristic point where the n-h curve crosses itself is shifted significantly from around 0.2 to around 0.5 (also due to line width not well identifiable)
• It is totally unclear how you motivate the change in the equation for tau_n.
• This clear discrepancy of the n-h curve to the original and the replotting of a fitted function f(h) unrelated to the shown n-h curve does satisfy any criteria of a replication.
• If these results are indeed not replicable the authors should state so. However, they should put a reasonable effort into an adequate replication (e.g. by motivated changes in the equations/implementation) and document this process in detail.
• For the trace of f(h) in Figure 1C, you simply re-plot a polynomial with the reported fit parameters. However, these fit results are a result of the paper and should themselves be replicated.

Figure 1D1

• The time axis is missing
• Missing comparison of essential features and why you consider the plot replicated.
• Please comment on the initial dip of the trace which is lower than the others.

Figure 1D2

• The time axis is missing
• Missing comparison of essential features and why you consider the plot replicated.
• Please comment on the increasing peak and decreasing trough values. This trend is also influenced by the change you applied to the time constant function for n, please consider this as well.

Figure 2

• missing time axis
• discrepancy in the scale of the voltage, please mention and discuss
• please discuss the shorter 'ringing' time in you simulation

Figure 3

• For a full replication, the figure is missing the part showing the instantaneous firing rate after the current step.
• For figures 3B the relation between figures 3A and the choice of hs values does not become clear. Please, motivate this in the caption.
• Figure 3C
  • This figure is missing an explanation of its elements, especially the newly introduced ones (either as legend or in the caption): solid line, dashed line, pink line, H2, red points.
  • Why did you choose to extend the hs axis with respect to the original?
  • the grey trajectory line covers the most crucial point of the plot!
  • The linewidth of the trajectory and the contraction of the region of interest due to the larger scale make the trajectory dynamics unnecessary hard to see.
• For handling the ODEs you use scipy functions and PyDSTool. Please, state this explicitly in the text and explain what exactly you use them for and why you chose them.

Figure 4

• missing h axis for Figure 4 A2
• missing arrow or other indication in which direction the current is increased. Just saying from left to right or right to left is not very understandable on its own, especially since you didn't put in the fixpoints.
• Figure B: missing explanation of the elements!
• In Figure B2 again the region around the subcritical Hopf point is covered by large red markers and overlapping, unreadable labels.
• Missing current units for Figure B
• For Figure 4C, please state which current and which voltage you plot, why this is of interest.
• Missing current axis for Figure 4 C2

Figure 5

• missing value of depolarization in all plots
• Again, there is missing an explanation of the key results which are represented in these plots, which you are trying to replicate. Without this description, a reader can a) hardly understand the figure and b) can not reasonably assess if and why the replication failed.
• Your code also doesn't allow for easy bug tracking here, due to not understandable or even misleading variable names (the solution should be a voltage but is named ic), missing clarifying comments, and convoluted structure.

Article text

• Why do you cite the original paper as 'Levitan et al. (2014)' although the first author is Qian?
• Commonly in articles, you present you the work in the present tense. Here, it can make sense to refer to the original work in the past tense, but at for least your results and figures it would make sense to use the present tense. Please review, and use it consistently.
• You use 'replicated' and 'reproduced' interchangeably throughout the text. Please adhere to the definitions of ReScience use the correct terminology. Reproduction: Running the same computation on the same data; Replication: repeating a computational study in the spirit of a published protocol, but with a different realization. (https://rescience.github.io/faq/)
• You use a lot of emphasizing adjectives like 'successfully' or 'faithfully'. Instead of using these formulations, describe the replication in a more neutral and precise manner. State exactly the argumentation why and to which degree you consider a result as replicated.
• reference the location of the equation you modified in the original paper
• There is a wrong minus sign in equation (compare to the original)! It is however correct in the code.
• explain the variables in the equation
• You need to explain you motivation why you modified this time constant equation and you need to argue why exactly you did change this numeric value from 40 to 60, instead of many other changes you could have made, and why the modified version should be considered 'better' than the original.
• Figure 1 should be moved into Results section, after the first paragraph
• "All original figures, [..], were successfully replicated". Not true. At least the last figure was not replicated!
• "..., successfully showing that they had not compromised the original sodium current description." Elaborate. What constitutes 'not compromised' here?
• "seemingly identical results" Be precise how they differ, how you compared the results, and on what basis (according to which features) you consider them equal.
• 'the same but “backwards” or ”flipped”' Be more precise. Along which axis? According to which attribute you consider it being 'flipped'?
• In many sentences there are missing articles, for example:
  • ".. when depolarizing current was repeatedly applied"
  • ".. modified potassium current"
  • ".. plotted membrane current"
• "in our implementation". Either 'with' our implementation, or in our simulation.
• In the discussion section, you mention that the original code contains different parameters than the paper. Please, state also these parameters and include this in the discussion about which equation to choose, and why and how to modify it.
• Use consistent Figure names. Either 'Figure 2 B1' or 'Figure 2B1'.
• "More specifically, our replication of Figure 1B was “backwards” or ”flipped” where in our replication the sodium current amplitude in time, which is opposite in the original implementation where it decreases towards 0 in time." This sentenced is very unclear, please re-formulate.
• "Hence, we can only speculate that this may have potentially been a graphing error." There are clear qualitative differences in the plots (see above) which can not be explained by 'graphing error', please elaborate.
• "Overall, this implementation replicated the original paper."/"With the exception of Figure 6, this implementation confidently shows that the original results were correctly implemented" No, you present a partial replication and found several discrepancies. Be more precise.

[rougier]

Thanks @rgutzen for this very detailed review. @mupsh You can wait for the second review or you can start modifying your article, taking @rgutzen comments and suggestions into account.

[maryupshall]

@rgutzen Thank you for your extensive comments and we're working on changes. @cJarvers Since we're going to be making fairly extensive changes it might be easiest for us to add our changes so you don't have to go through 2 sets.

[cJarvers]

Sounds good. I'm looking forward to your revisions.

[rougier]

@cJarvers Gentle reminder

[rougier]

@mupsh Did you have time to update the manuscript?

[maryupshall]

@rougier It is still in the works. I sincerely apologize for the delay, it's been a very busy month. Considering the holidays, would it be ok if I got it back to you for the first week of January?

[rougier]

Sure. @cJarvers Would that be ok for you ?

[cJarvers]

Yes, first week of January is fine with me.

[rougier]

@mupsh Gentel reminder. Can you do the review by Monday?

[cJarvers]

Just to make sure we're all on the same page: @rougier, would you like the review by me to be done on Monday or the revisions by @mupsh? I have not done the review yet since the authors suggested I wait for their revisions first. If that is a problem, I can try to do the review by Monday.

[maryupshall]

@rougier We can get the updated manuscript by Monday, no problem.

[rougier]

@cJarvers Sorry, the message was adressed to @mupsh. @mupsh Monday would be perfect and then @cJarvers can review the revision.

[maryupshall]

Sorry for the delay in getting this out, and thanks to @rgutzen for his thorough review. Overall we've implemented most of the requested changes and believe this puts the manuscript in a better state. For clarity we've included the original comments along with our reply. @cJarvers we're ready to move on with the review now.

General Comments

• First of all, there is only a brief introduction, which doesn't answer the interesting question of what your motivation is for replicating exactly this study? And why is this a relevant paper? What are the key insights? And why would it be beneficial for the scientific community to replicate these results?

  • We’ve expanded the introduction to give context for why these neurons and depolarization block are of interest with regards to schizophrenia and the authors conclusions.

• There is no information on the reproducibility of the original implementation. Is the original data available? Is the original code available? Only at the very end, you mention the communication with the authors. This is very relevant to follow your replication and should be moved to the beginning of the paper and elaborated on.

  • We mention that they’ve provided us with the model implementation they had available.

• Ideally, you can also ask the original authors, if they agree to publish the code they provided you alongside this replication study.

  • We don’t believe this would be of interest as the code provided is only the implementation (no analysis of manipulations) of the 3d model

• Your paper is rather short, and writing short precise papers is generally desirable. However, you mostly just state what you replicated, and I'm totally missing a description of your journey how you replicated the results. What was the approach? What are the design principles or your implementation? Which difficulties did you encounter? What did not work? These aspects, I would consider the most valuable part of your research, and the ones others can really benefit and learn from.

  • We appreciate that our original manuscript was too sparse, we’ve elaborated on the replication in the modified version. We are unclear as to what you mean by design principles or approach, this could mean many things most of which we believe are outside the scope of a replication paper.

• Especially, when you are not able to replicate a result or figure (exactly, or at all), it is not sufficient to say that you tried, it didn't work out, and then speculate that it might be due to different parameters. There needs to be a more detailed account of what exactly you tried, why you think the replication fails, and generally, you need to present a reasonable effort in investigating and fixing the failed replication.

  • We provided more detail, and a justification of why the original parameters are incorrect.

• For each figure I'm missing an explanation of what exactly is shown and why this is interesting. This doesn't need to be long, but it is required in order to make sense of the replication, as this is the logical basis on which anyone can then start evaluate the similarity to the original results.

  • Added explanation of what is shown.

• For generating the simulation results, additional simulation parameters are set. Instead of hiding them in the code, they should be reported, and if necessary, discussed (e.g time resolution, input, discarded initialization time, ...)

  • We’ve elaborated on our parameters in the methods section.

Code

General

• The main goal of replication studies is to validate results from other publications by reimplementation. However, another goal is also to provide understandable and easily reproducible code so that others get better insights into the computational aspects of the original study and reuse/edit/expand the analysis or part of it. The second goal requires the code to be in a certain format. If possible, it should use standards and existing software and packages which are common in the field. Furthermore, it needs to be logically structured and well documented. Writing custom code to replicate other custom code misses the opportunity of adding real value to the implementation.

• So I strongly encourage you to use standard resources for the neuron models, like for example PyNN or NeuroML

  • While in general we agree with this statement, in this case we have to disagree. Due to the issues we’ve faced in replicating the manuscript, not having duplicate representations of the model (i.e. for time-domain solutions, and bifurcation analysis) is critical because if one implementation works but the other does not that adds another level of potential implementation error.
  • Since (to our knowledge) no neuron specification format can be parsed by any numerical bifurcation engine. Our choices are either implement the model in python and inject symbolic variables to prevent duplication which we’ve done. Or, write our own parser which would add far more problems than it would solve. For this reason, we’ve chosen to keep our implementation as is.

• Besides documenting the code and structuring it in a way that makes it easy to read and understand, please also explain your code in the article. This is a computational paper and addressing and discussing certain implementation decisions in the text can be very helpful. The main outline of the implementation and the design principles could even be presented in their own Section.

  • We agree entirely that our initial submission should not have been presented in the state that it was in. We’ve now thoroughly refactored it and documented it properly. That being said, we’re unsure what you mean by discussing design principles in the paper. Since the goal of the work is, in-essence, to validate the original work. We’re unclear as to the benefit of discussing for example OOP vs imperative implementations. Could you please provide an example of what you would like?

• For the packages you do use, make sure you also mention them in the article, state what they do and why you chose them, especially for less common packages.

  • We’ve implemented this in the methods.

Details

• executing run.py in the provided environment return an error while producing figure 3: unable to execute 'swig': No such file or directory error: command 'swig' failed with exit status 1 After installing swig version 3.0.12 it runs through.

  • We didn’t catch this omission and it’s now specified in our readme.

• For understandability, it helps to assign variables with telling names instead of using numerical values in equations without explanations.

• code lines like end_times = np.ravel(list(zip(np.arange(3, last_set + 100 + 1, 100), np.arange(100, last_set + 100 + 4, 100))))

  • Very true, this again is an oversight on our behalf, and we’ve refactored out these obscure definitions.

• for some reason you named tau_hs different from the other variables as 'tah_hs'

  • This was a typo and it’s been corrected.

• I strongly suggest making use of packages to handle units (like quantities or pints) at least when variables are passed to functions. This easily ensures that the right magnitudes and units are used. Once the units are check in a function, they can then be omitted for the actual computations to not slow them down. Otherwise, it can get quite bothersome to trace back which magnitude the resulting values have.

  • In general, yes this is true. However, in this case since the parameter scheme is consistent, except when the authors give “somewhat obviously incorrect units” (which a computational unit system would not help) such as nS/cm^2 we believe won’t really help in this case.

• For a better overview and less error-proneness, it is generally advisable to outsource parameter settings (into a yamls, json, or event txt file), and import them and not set them by hand in each function.

  • Since this model implementation is not dynamic i.e. it’s not changing. Many benefits of an external config are lost. There are two options in this case for implementing config files, either duplicating all parameters in each model’s config which we don’t want to do (see neuroML comment for duplication argument) or each specification would only be used for setting non-defaults. In which case you’d still have to check back to either the default function or the default config file which in our view is just as cumbersome as looking up the function call for non-default arguments. For this reason, we will leave our implementation as is

• Also the parameters list is rather fragile since it relies on an arbitrary but specific order of the parameters. A dictionary would be a better choice.

  • This is completely true and we’ve implemented this.

Replicating the figures

General

• Using additional ticks (and tick labels) on the axis and would make the evaluation of the figures the comparison much easier.
• Use smaller linewidth, to not cover up the smaller details.
• Explain all the elements in each figure!
• Emphasize the relevant points in the plots. Discuss which aspects of the plots are the relevant results you aim to replicate?
  • We’ve cleaned up our plots where possible and discussed them in the manuscript. We also made sure to explain all the elements of each figure in the figure captions or using a legend.

Figure 1a

• A legend is missing which of the curves is with and without hs.
• The code contains the following comment. # need to divide the given value by 1.5 to get the correct graph. This operation is not motivated or even mentioned in the text.
• The vertical axis misses its units
• I am missing a statement why you consider the plot replicated, like giving the agreement of characteristic points (e.g. minimum, horizontal trend in the beginning, Voltage at current drop-off, skewed parabola shape around minimum)
  • We’ve implemented these changes as requested.

Figure 1b

• even when the plot is mirrored in the x-direction, there are still clear differences (absolute and relative peak sizes, the distance between peaks) in the plot which needs to be addressed!
• Without a reason why the simulation would be backward in time, and given the other differences in the plot, the resulting plot suggests that there are other discrepancies in the equations or the code which need be uncovered. In this state, I would consider the plot not replicated.
• The code contains comment # need to divide given value by 2 to get correct graph. However, this is not mentioned in the text.
• The corresponding code is quite messy and thus hard to debug. To illustrate, I try to understand how the current trace which is plotted evolves. You are plotting the clamp_current which is calculated by sodium_current(), which takes v,h,and hs as state/time-dependent arguments. These are calculated by solving the voltage_clamp ODE using ic initial condition in the time packet time. Looking into voltage_clamp(), clamp() is called which calls the ODE function from the parameters list (here ode_3d()). ode_3d() itself is not easy to understand and debug since it has several not intuitive equations which stretch over multiple line. But most remarkably ode_3d() does not use its time argument t, which essentially means that the whole current function does not depend on the time packet it is iterated over. The main point here is that you organize and structure the code obfuscates dependencies and makes it really hard to understand what is happening in the code.
  • With regards to the “mirroring” upon re-reading their definition and properly re-implementing the code we’ve fixed this issue.
  • As discussed in the text there is a scaling issue which we’ve addressed as well as we could by manipulating the experimental parameters
  • Your comments about the code quality and smells are correct and we’ve refactored these. To clarify what’s happening in the call stack is that: clamp (ignore the dunder) is generic function for clamping a variable and voltage_clamp is the specific wrapper for a voltage clamp (instead of clamping another variable). The comment about not using time is subtle we did not do a good job explaining it in the code in the initial version. Since the pulses of current or clamping are discrete, they shouldn’t be implemented as a non-autonomous function, therefore it’s implemented as a sequence of solved autonomous ODEs with the initial condition of the subsequent step being the end of the current step, hence why “t” is unused. We’ve documented and refactored this in a way that we believe makes this clearer.

Figure 1c

• it should be made clear that the time course of n(h), which is not a function, is fitted by a function f(h)
• "to get the same results" Figure 1C is visibly different to the original!!
• The n-h relation in the simulation varies over time. See the spread of the curve in the original plot. This is not visible in the replication and not commented.
• The characteristic point where the n-h curve crosses itself is shifted significantly from around 0.2 to around 0.5 (also due to line width not well identifiable)
• It is totally unclear how you motivate the change in the equation for tau_n.
• This clear discrepancy of the n-h curve to the original and the replotting of a fitted function f(h) unrelated to the shown n-h curve does satisfy any criteria of a replication.
• If these results are indeed not replicable the authors should state so. However, they should put a reasonable effort into an adequate replication (e.g. by motivated changes in the equations/implementation) and document this process in detail.
• For the trace of f(h) in Figure 1C, you simply re-plot a polynomial with the reported fit parameters. However, these fit results are a result of the paper and should themselves be replicated.
  • We’ve implemented the changes as requested however to address a few specific points. And we outright say that we could not match the authors’ implementation
  • The n-h limit cycle (which is a limit cycle) appears to change over time as they’ve included a transient. And we don’t discuss this as 1) that’s standard for displaying limit cycles and 2) there are far greater issues we wanted to focus on.
  • The main issue here is that the function tau_n cannot be the correct function due to an invalid domain (see discussion in the paper) we’ve motivated our modification and furthermore show in the supplemental the result with their original equation.
  • With regards to their polynomial, we do re-fit it (terms are printed during the run for figure 1) since we’ve established that the implementation cannot be correct, the n-h space space will be quite different as well (which it is), therefore any fit polynomial will be quite different. The important part is that in using their reduced polynomial model elsewhere we replicate those results.

Figure 1D1

• The time axis is missing
• Missing comparison of essential features and why you consider the plot replicated.
• Please comment on the initial dip of the trace which is lower than the others.
  • Implemented as requested, the initial dip was a transient which has now been discarded as discussed in the MS.
  • Additionally, we had flipped the location of 1D1 and 1D2 this has been corrected.

Figure 1D2

• The time axis is missing
• Missing comparison of essential features and why you consider the plot replicated.
• Please comment on the increasing peak and decreasing trough values. This trend is also influenced by the change you applied to the time constant function for n, please consider this as well.
  • Again, this trend was a transient which has been corrected.

Figure 2

• missing time axis
• discrepancy in the scale of the voltage, please mention and discuss
• please discuss the shorter 'ringing' time in you simulation
  • Addressed & discussed.

Figure 3

• For a full replication, the figure is missing the part showing the instantaneous firing rate after the current step.
• For figures 3B the relation between figures 3A and the choice of hs values does not become clear. Please, motivate this in the caption.
• Figure 3C
  • This figure is missing an explanation of its elements, especially the newly introduced ones (either as legend or in the caption): solid line, dashed line, pink line, H2, red points.
  • Why did you choose to extend the hs axis with respect to the original?
  • the grey trajectory line covers the most crucial point of the plot!
  • The linewidth of the trajectory and the contraction of the region of interest due to the larger scale make the trajectory dynamics unnecessary hard to see.
• For handling the ODEs you use scipy functions and PyDSTool. Please, state this explicitly in the text and explain what exactly you use them for and why you chose them.
  • We’ve added the figure components as discussed. We’ve also discussed their use. However, we don’t discuss why we’ve used them as really there is no other viable option. We are not aware of any other established ode solvers in python asides from scipy odeint. And there are a few XPPAUT / AUTO python interfaces however they’re all at least several years without maintenance and are often undocumented, hence PyDSTool is the only viable option. It seems slightly unnecessary to state we used “x” because there’s nothing else.
  • For the relation between figures 3A and the choice of hs values in 3B, these come from the “distinct” regimens in A2. We don’t address this as the goal is to show the change in stability over the depolarization block hence the parameters are somewhat arbitrary. Furthermore, this a scientific discussion which we feel is outside the domain of the replication.
  • We’ve addressed the other display points.
  • Added explanation of the elements in the figure caption.

Figure 4

• missing h axis for Figure 4 A2
• missing arrow or other indication in which direction the current is increased. Just saying from left to right or right to left is not very understandable on its own, especially since you didn't put in the fixpoints.
• Figure B: missing explanation of the elements!
• In Figure B2 again the region around the subcritical Hopf point is covered by large red markers and overlapping, unreadable labels.
• Missing current units for Figure B
• For Figure 4C, please state which current and which voltage you plot, why this is of interest.
• Missing current axis for Figure 4 C2
  • Added explanation of the elements in the Figure caption.
  • Fixed and added explanation of the elements.
  • With regards to the bifurcation we are somewhat limited in terms of control of elements and features (such as plot order). That said we’ve done what we can to maximize readability.

Figure 5

• missing value of depolarization in all plots
• Again, there is missing an explanation of the key results which are represented in these plots, which you are trying to replicate. Without this description, a reader can a) hardly understand the figure and b) can not reasonably assess if and why the replication failed.
• Your code also doesn't allow for easy bug tracking here, due to not understandable or even misleading variable names (the solution should be a voltage but is named ic), missing clarifying comments, and convoluted structure.
  • Again, we agree the code was confusing and should not have been shipped in that state. It’s been refactored. While no-one should have been expected to understand the code as is. The “ic” (lines 46 to 49) was again used as the initial condition for the next step (due to the discrete conductance jump) and the state (i.e. the current solution) was concatenated with the growing full solution. The new implementation makes this much clearer.

Article text

• Why do you cite the original paper as 'Levitan et al. (2014)' although the first author is Qian?

  • This was a copying error; we’ve fixed the references. @rougier we can fix the repo name, is it O.K. to edit this submission to reflect this?

• Commonly in articles, you present you the work in the present tense. Here, it can make sense to refer to the original work in the past tense, but at for least your results and figures it would make sense to use the present tense. Please review, and use it consistently.

• You use 'replicated' and 'reproduced' interchangeably throughout the text. Please adhere to the definitions of ReScience use the correct terminology. Reproduction: Running the same computation on the same data; Replication: repeating a computational study in the spirit of a published protocol, but with a different realization. (https://rescience.github.io/faq/)

• You use a lot of emphasizing adjectives like 'successfully' or 'faithfully'. Instead of using these formulations, describe the replication in a more neutral and precise manner. State exactly the argumentation why and to which degree you consider a result as replicated.

• reference the location of the equation you modified in the original paper

• There is a wrong minus sign in equation (compare to the original)! It is however correct in the code.

• explain the variables in the equation

• You need to explain you motivation why you modified this time constant equation and you need to argue why exactly you did change this numeric value from 40 to 60, instead of many other changes you could have made, and why the modified version should be considered 'better' than the original.

• Figure 1 should be moved into Results section, after the first paragraph

• "All original figures, [..], were successfully replicated". Not true. At least the last figure was not replicated!

• "..., successfully showing that they had not compromised the original sodium current description." Elaborate. What constitutes 'not compromised' here?

• "seemingly identical results" Be precise how they differ, how you compared the results, and on what basis (according to which features) you consider them equal.

• 'the same but “backwards” or ”flipped”' Be more precise. Along which axis? According to which attribute you consider it being 'flipped'?

• In many sentences there are missing articles, for example:

  • ".. when depolarizing current was repeatedly applied"
  • ".. modified potassium current"
  • ".. plotted membrane current"

• "in our implementation". Either 'with' our implementation, or in our simulation.

• In the discussion section, you mention that the original code contains different parameters than the paper. Please, state also these parameters and include this in the discussion about which equation to choose, and why and how to modify it.

• Use consistent Figure names. Either 'Figure 2 B1' or 'Figure 2B1'.

• "More specifically, our replication of Figure 1B was “backwards” or ”flipped” where in our replication the sodium current amplitude in time, which is opposite in the original implementation where it decreases towards 0 in time." This sentenced is very unclear, please re-formulate.

• "Hence, we can only speculate that this may have potentially been a graphing error." There are clear qualitative differences in the plots (see above) which can not be explained by 'graphing error', please elaborate.

  • We’ve addressed the above concerns.

• "Overall, this implementation replicated the original paper."/"With the exception of Figure 6, this implementation confidently shows that the original results were correctly implemented" No, you present a partial replication and found several discrepancies. Be more precise.

  • Upon re-reading the definitions of replication, it also notes that we should be able to change something which “shouldn’t matter” and see if the conclusions are the same. In this spirit while we believe we’ve identified some parameter issues, and unit presentation issues the qualitative behaviour is the same once corrected and importantly the conclusion remains the same. As such we will still refer to this as replicated unless we’ve misunderstood. @rougier do you have any thoughts on this matter?

[rougier]

Impressive ! @cJarvers can we set a deadline for your review ? Like two weeks, would that be ok ?

[cJarvers]

@rougier Yes, two weeks should be fine.

[cJarvers]

Okay, @mupsh , here's my review:

This review is relative to commit https://github.com/mupsh/ReScience_Levitan_2014/commit/d5c641ffd14d5740e032e2864083215ac49dec73

Reproducibility

I ran the code using Python 3.6.10, swig 4.0.1, PyDSTool 0.90.3, numpy 1.18.1, scipy 1.3.0 and matplotlib 3.1.1.

The code produced a warning message ("Test function going crazy: class Hopf_Bor(BordermMethod, BiAltMethod)"), but otherwise ran without problems and produced the figures as shown in the paper.

Replication

The core results of the original paper are replicated. This requires some major parameter changes, so I'm not sure whether I share your optimistic conclusion (that "the 2D and 3D models were implemented correctly in the original work"), but I find the explanations of the discrepancies and your overall reasoning clear enough to accept this as a replication - thanks to your work I'm sufficiently confident that the conclusions of the original article were correct, even if there were errors in the implementation or the reporting.

Clarity of Article

In general, I find the article is now a good deal clearer than before the revision. The explanations of what each figure shows and what adjustments were necessary to replicate each are easy to follow.

Unfortunately, the article is not self-contained. In order to understand the explanations of the replication, it is necessary to know the original article, which is behind a paywall. I could not find any guidelines about how self-contained the article should be (do you know any, @rougier?) and I guess most readers interested in a replication would know the original. Still, I think it would be good have a section that explains the model itself. This would enable to reader to better understand what is going on, without going back to the original article.

I also concur with @rgutzen that adding a section that explains the code would be helpful. This should not include overly general things like OOP vs. imperative programming (and I don't think that was meant with "design principles") and should be more aimed towards giving the reader a starting point to exploring/understanding the code. This would essentially be a more elaborate version of the first paragraph of your methods section (explaining why and how you used which package), including some pointers like the fact that you factorized your code into the differential equations for all three models (diff_eq.py), the steady-state equations and time constants (gating.py), the membrane currents (current.py) and nullclines (nullclines.py).

Another thing that could improve the clarity of the article would be to list all the changes made to the original model and between different figure in one central place, e.g., a table. Unless I missed something, these changes are:

• the changed formula for tau_n (used everywhere)
• the value of g_Na in Figure 1A
• the value of g_Na and the pulse width in Figure 1B
• the newly fitted polynomial coefficients in Figure 1C; the fitted values should be reported in the text
• the reversal potentials for NMDA and AMPA receptors in Figure 6
• the conductances g_AMPA and g_NMDA in Figure 6A and 6B (different values)

Especially for the sodium conductance, it would be good to know what value you used for different simulations, the one from 1A or 1B.

Otherwise, I noticed a few smaller issues:

• Use of tenses is sometimes inconsistent. For example, in the second paragraph one sentence starts with "Previous work had", whereas the next starts with "Previous work has".
• The formulation "This new model included, or exclused a slow inactivation term" might confuse readers. It's not clear that this relates to the 2D and 3D versions mentioned before.
• The legend in Figure 1A covers a small part of the curve.
• The "ringing" in Figure 2 is hardly visible. It might help to make the figure bigger and possibly decrease the line width.
• In the text, the term "ringing" is not explained, except by referring to Figure 2. A brief definition might help.
• In the explanation for Figure 4, you mention that you use a different naming convention (SNP vs. LPC) but do not explain why.
• The last sentence of section 2.4 is incomplete.

Clarity of Code

The code is well structured and commented. The separation into simulation code, runners for the individual figures and plot setup works well and helped me find my way around the code.

In two places, this wasn't quite true: The diff_eq.py module contains a lot of stuff, some of which is not directly related to the differential equations and could for example go to an experiments.py module (e.g. the clamp experiments or IV curves). The other place is the functions for figures 3c and 4b, due to the use of pyDSTool. I don't really know how to separate that out better, though.

I also agree with @rgutzen that using packages for handling / checking units would be useful. This is especially true since some of the differences from the original paper involve scaling, so it would be good to make sure this is not due to some implicit unit conversion error somewhere.

One minor point I noticed: the explanation of parameter parameter_name in diff_eq.pulse is only "Parameter to"; there's probably something missing there.

Summary

Overall, I'm quite happy with the paper. I think some improvements are possible: (1) explaining the model in the text to make the article more self-contained, (2) listing changes and parameter values in one place / table and (3) adding explicit unit handling to the code.

[rougier]

thanks @cJarvers. @mupsh can you address the comments ?

@mupsh Don't forget to change the title!

@cJarvers I think you're right and we may need to add a section for aiming at self-contained article when the original article is not available (behind a paywall). Could you propose some changes authors instruction based on your experience ?

Note I found this link which seems to be free: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4254877/ (or is it just me ?)

[maryupshall]

@rougier We're working on implementing the requested changes.

[cJarvers]

@rougier True, the ncbi-link is free. Question is whether green OA like this solves the issue, since readers still need to find the free version. Maybe adding a link to the replication article would work. Anyways, I'll open a separate issue to suggest some changes to the author instructions.

[maryupshall]

Sorry again for the delay in getting this out, and thanks to @cJarvers for his thorough review. As before, we've implemented most of the requested changes and believe this puts the manuscript in an even better state. For clarity we've included the original comments along with our reply.

Reproducibility

I ran the code using Python 3.6.10, swig 4.0.1, PyDSTool 0.90.3, numpy 1.18.1, scipy 1.3.0 and matplotlib 3.1.1.

The code produced a warning message ("Test function going crazy: class Hopf_Bor(BordermMethod, BiAltMethod)"), but otherwise ran without problems and produced the figures as shown in the paper.

• That message is normal. pyDSTool does not have a silent mode or verbosity setting. We have made a note of this in the README.

Replication

The core results of the original paper are replicated. This requires some major parameter changes, so I'm not sure whether I share your optimistic conclusion (that "the 2D and 3D models were implemented correctly in the original work"), but I find the explanations of the discrepancies and your overall reasoning clear enough to accept this as a replication - thanks to your work I'm sufficiently confident that the conclusions of the original article were correct, even if there were errors in the implementation or the reporting.

Clarity of Article

In general, I find the article is now a good deal clearer than before the revision. The explanations of what each figure shows and what adjustments were necessary to replicate each are easy to follow.

Unfortunately, the article is not self-contained. In order to understand the explanations of the replication, it is necessary to know the original article, which is behind a paywall. I could not find any guidelines about how self-contained the article should be (do you know any, @rougier?) and I guess most readers interested in a replication would know the original. Still, I think it would be good have a section that explains the model itself. This would enable to reader to better understand what is going on, without going back to the original article.

• We added a brief paragraph explaining the model itself give that the original article is open access.

I also concur with @rgutzen that adding a section that explains the code would be helpful. This should not include overly general things like OOP vs. imperative programming (and I don't think that was meant with "design principles") and should be more aimed towards giving the reader a starting point to exploring/understanding the code. This would essentially be a more elaborate version of the first paragraph of your methods section (explaining why and how you used which package), including some pointers like the fact that you factorized your code into the differential equations for all three models (diff_eq.py), the steady-state equations and time constants (gating.py), the membrane currents (current.py) and nullclines (nullclines.py).

• We have implemented an explanation of our codebase from a big-picture standpoint.

Another thing that could improve the clarity of the article would be to list all the changes made to the original model and between different figure in one central place, e.g., a table. Unless I missed something, these changes are:

the changed formula for tau_n (used everywhere) the value of g_Na in Figure 1A the value of g_Na and the pulse width in Figure 1B the newly fitted polynomial coefficients in Figure 1C; the fitted values should be reported in the text the reversal potentials for NMDA and AMPA receptors in Figure 6 the conductances g_AMPA and g_NMDA in Figure 6A and 6B (different values)

• We've added a table with a summary of the changes. Because we doe not use our own fit values in figure 1C we have elected not to report them. The fit values are reported in the output of the program for interested parties but we believe it would confuse things to report it in the manuscript as we have the author's reported parameters, the parameters we got form the authors in their XPP file, and this would be a third set to report. Considering the risk of confusion coupled with the fact that we generate them, but do not use them, we did not include our fits.

Especially for the sodium conductance, it would be good to know what value you used for different simulations, the one from 1A or 1B.

Otherwise, I noticed a few smaller issues:

Use of tenses is sometimes inconsistent. For example, in the second paragraph one sentence starts with "Previous work had", whereas the next starts with "Previous work has". The formulation "This new model included, or exclused a slow inactivation term" might confuse readers. It's not clear that this relates to the 2D and 3D versions mentioned before. The legend in Figure 1A covers a small part of the curve. The "ringing" in Figure 2 is hardly visible. It might help to make the figure bigger and possibly decrease the line width. In the text, the term "ringing" is not explained, except by referring to Figure 2. A brief definition might help. In the explanation for Figure 4, you mention that you use a different naming convention (SNP vs. LPC) but do not explain why. The last sentence of section 2.4 is incomplete.

• We've addressed the above concerns in the manuscript and the article. With regards to bifurcation names: the same mathematical description often has many "English" names. In this case the two bifurcations are the same; pyDSTool outputs LPC as text instead of SNP which AUTO uses. pyDSTool additionally does not have the flexibility to rename the labels documented. Considering the equivalence, we did not feel it necessary to explain the "grammatical difference".

Clarity of Code

The code is well structured and commented. The separation into simulation code, runners for the individual figures and plot setup works well and helped me find my way around the code.

In two places, this wasn't quite true: The diff_eq.py module contains a lot of stuff, some of which is not directly related to the differential equations and could for example go to an experiments.py module (e.g. the clamp experiments or IV curves). The other place is the functions for figures 3c and 4b, due to the use of pyDSTool. I don't really know how to separate that out better, though.

I also agree with @rgutzen that using packages for handling / checking units would be useful. This is especially true since some of the differences from the original paper involve scaling, so it would be good to make sure this is not due to some implicit unit conversion error somewhere.

• We've refactored our more experimental functions into experiments.py. Additionally, we have implemented units where prior to computation units are forced into mV, uA/cm^2, or mS/cm^2 (depending on the natural unit used) and the unit is discarded. We experimented with doing unit checking more robustly i.e. running one iteration of the ODE to ensure unit matching. Unfortunately based on the amount of jumping between ODE and "analysis" type code we do, that became unwieldy and less readable. We believe our approach finds a balance.

• With regards to the bifurcation ODE (3c/4b) we agree entirely. Unfortunately, due to the level of fine-tuning required for each plot there was not an easy refactor for this. Form a non-refactor standpoint the code presented maps fairly nicely onto the natural steps someone familiar with XPPAUTO might use to reproduce them. For this reason, we've left the code in a more "scientific" state.

One minor point I noticed: the explanation of parameter parameter_name in diff_eq.pulse is only "Parameter to"; there's probably something missing there.

• Fixed.

[cJarvers]

Thanks for your reply, @mupsh. I'm quite happy with how you addressed my comments.

I ran the code again (with slightly different versions of the dependencies: python 3.8.2, scipy 1.4.1, matplotlib 3.2.1, sympy 1.5.1, Pint 0.11, PyDSTool 0.91) and everything behaves as expected.

There are only two minor things I noticed in the article:

• In the very first line, the word "dopaminergic" stretches over the right border. Not sure how hard that is to fix.
• Reference 11 seems to be incomplete.

@rougier After these minor things are fixed, I'd recommend acceptance.

[rougier]

@cJarvers @rgutzen Thanks for your very nice and extensive reviews. I therefore accept the manuscript that will be published soon hopefully. Congratulations @mupsh 🎉 !

I will make a PR to your repository to complete your metadata. In the meantime, can you register your code repository at https://www.softwareheritage.org/save-and-reference-research-software/. You should obtain a swid that you can post here.

[maryupshall]

@rougier Great news, thank you very much! We've made the minor fixes, I've registered the code as requested and obtained the following SWH-ID: https://archive.softwareheritage.org/swh:1:rev:e1040fe3fbaefef4166964772efa08ada444db6e;origin=https://github.com/mupsh/ReScience_Qian_2014.git/

[rougier]

I created a PR (https://github.com/mupsh/ReScience_Qian_2014/pull/10) for small changes to your repo before publication. Can you check it?

[rougier]

Ok, I just saw PR has been merged. I'll try to publish the paper by the end of the week. Don't hesitate to spam me if it's not done :)

[maryupshall]

@rougier I might be mistaken, but spamming you just in case :)

[rougier]

It has been published and online ! Sorry for the delay. See https://zenodo.org/record/3842360