ReScience
[Re] Model of thalamocortical slow-wave sleep oscillations and transitions to activated states
Original article: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6757797/
PDF URL: https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/a732939b013c7966631ef9e742c98d29ba311934/Submission/%5BRe%5DReynesAussel2024.pdf Metadata URL: https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/a732939b013c7966631ef9e742c98d29ba311934/Submission/metadata.yaml Code URL: https://github.com/Mathilde-Reynes/ReynesAussel2024/tree/a732939b013c7966631ef9e742c98d29ba311934/Model
Scientific domain: Computational Neuroscience Programming language: Python (brian2) Suggested editor: Benoît Girard Suggested reviewers: Qihong Lu; Ozan Caglayan; Pierre Enel
@benoit-girard Can you edit this submission, I think it's a perfect fit for you.
@benoit-girard @apdavison @heplesser @vitay @mstimberg @pietromarchesi @stephanmg
Can your review this submission on "Model of thalamocortical slow-wave sleep oscillations and transitions to activated states"
I'd be happy to review this, too (but certainly do not mind leaving it to the two who agreed to review earlier, either) :blush:
Thank you very much.@apdavison @heplesser @mstimberg.
@apdavison @heplesser (since you were first to answer) you can start the review (and if you can do it before Christmas that would be great) and feel free to ask questions if needed. @mstimberg Thanks again, I think I've other submission waiting for reviewers.
@apdavison @heplesser Any chance to get your reviews by next Monday (before Christmas)?
@rougier I'm afraid I don't think it will be possible before Christmas. I've pencilled this in for the first full week in January, I'm targeting January 10th.
@benoit-girard @apdavison @heplesser @vitay @mstimberg @pietromarchesi @stephanmg
Can your review this submission on "Model of thalamocortical slow-wave sleep oscillations and transitions to activated states"
terribly sorry, somehow this went under the radar for me, if help is still needed, happy to assist!
I am also aiming for the first week of 2025.
From: Stephan Grein @.> Sent: Thursday, December 19, 2024 6:20:17 PM To: ReScience/submissions @.> Cc: Hans Ekkehard Plesser @.>; Mention @.> Subject: Re: [ReScience/submissions] [Re] Model of thalamocortical slow-wave sleep oscillations and transitions to activated states (Issue #87)
@benoit-girardhttps://github.com/benoit-girard @apdavisonhttps://github.com/apdavison @heplesserhttps://github.com/heplesser @vitayhttps://github.com/vitay @mstimberghttps://github.com/mstimberg @pietromarchesihttps://github.com/pietromarchesi @stephanmghttps://github.com/stephanmg
Can your review this submission on "Model of thalamocortical slow-wave sleep oscillations and transitions to activated states"
terribly sorry, somehow this went under the radar for me, if help is still needed, happy to assist!
— Reply to this email directly, view it on GitHubhttps://github.com/ReScience/submissions/issues/87#issuecomment-2555222763, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AAICWPD4L3JYVORMQIYGJG32GL55DAVCNFSM6AAAAABNGQ7AI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKNJVGIZDENZWGM. You are receiving this because you were mentioned.Message ID: @.***>
I have started reviewing the submission, and I will post more complete comments soon; however, I would like here to make a preliminary suggestion which will both (i) help me finish the review and (ii) make the re-implemented model more reproducible.
I was able to run python Thalamo_cortical.py, which completes without error. The problem is that to reproduce the data for the figures, I need to uncomment certain lines in Thalamo_cortical.py and, for some figures, also in other files. It takes a lot of reading and checking code to make sure the correct lines have been uncommented and re-commented.
If the authors have time, it would be very helpful to put the commented-out lines into a dictionary, or into functions, so that all the changes are controlled by a single parameter (ideally provided on the command-line, e.g., python Thalamo_cortical.py Figure4). Similarly, saving the raw data should be controlled by another parameter.
Hello,
Thanks for the feedback. We are currently implementing the changes asked to facilitate the review process. I will inform you when these are available.
Best
Dear all,
The changes proposed have been implemented. The git is up to date. In the meantime, we have identified an error in Figure 10, for which the middle and bottom panels were swapped (it was just a mistake in their position as the legend didn't match the figure). This issue has been corrected, and the revised manuscript has been updated as well.
Best
@apdavison Thansk fro you review. @heplesser Gentle reminder on your review. @Mathilde-Reynes Can you post here the changes you made?
Hello, please find below the changes implemented:
- One function controls everything necessary to run the model, found within Thalamo_cortical.py;
thalamocortical_network(seed_value,analyze_speed,fig_number,raw_data,plot_figure)with:
- analyze_speed (True/False) allows the user to compute the mean and std propagation speed of cortical up-states
- fig_number (str such as "9-A1","5" etc.) allows the user to run the simulation with the parameters necessary to reproduce a specific figure from the article
- raw_data (True/False) allows the user to select if raw data (.txt files) should be saved
- plot_figure (True/False) allows the user to plot the figure corresponding to the fig_number
- Addition of .py files containing functions necessary for 1. :
- Figure_conditions.py : Designed to configure the model with specific parameters adjusted for different figures.
- Figure_plotting.py : Designed to plot a specific figure of the article.
- Figure_rawdata.py : Designed to save the raw data necessary for a specific figure.
-
py files to plot every single figure separately from the raw data were removed as everything is now integrated in Figure_plotting.py. Figure 2, 3, 4 remain as they are based off raw data from the original model.
-
Figure 10 was updated in a new file called [Re]ReynesAussel2024_erratum.pdf. As previously said, the middle and bottom panels were swapped, as they were not matching the legend below.
Read me is up to date with all the information necessary to run the model with these new changes. Best,
I commend Reynes and Aussel for their throrough and careful effort to reproduce the work of Bazhenov et al. Overall, this work already seems in quite good shape, but I think manuscript and code would benefit from some improvements.
First of all, I found it somewhat difficult to match figures in the manuscript to figures in the original paper. I would strongly suggest to mention in each figure caption the corresponding figure number in Bazhenov et al. I also wonder if it would not be easier to follow the reproduction step by step if figures were in the same order as in the original paper. I noticed that the manuscript contains excerpts of the figures in the original publication. While this definitely eases comparison, I wonder if this could raise licensing issues.
Concerning the results, I am most concerned about Fig 5 in the present manuscript. The reproduction results show approximately 20 up and down states compared to only 11 in the original paper, i.e., twice the switching frequency than in the original, not just "minimally shorter up-states" as commented on p. 14. Also, in the lower part of Fig 5, the initiation of the transition looks much sharper in the original figure than in the reproduction, with very straight opening "flanks". Furthermore, the original figure only shows broad stripes spanning the entire PY population, while the reproduction figure also includes "tongues" that include only parts of the PY population. I believe this merits further investigation. In that context, I noticed that Fig 13 in the manuscript shows a similar difference in patterns, although here related to the presence and absence of thalamocortical connections. Perhaps something similar could be underlying the differences in Fig 5?
Concerning the sustained activity initially shown in Fig 6 of the manuscript, it would be interesting to see the corresponding plot for the current implementation. I wonder if initial conditions could play a role here?
I would also like to point to a number of more detailed points.
- You frequently refer to the ModelDB version of the Bazhenov et al model, but do not provide a reference for it. Please add the ModelDB accession number and a link.
- I_{D->S} = -I_{S->D}, isn't it? It would be useful to point this out.
- I think it would help the reader if you explicitly gave I_S^{int} somewhere above Eq 3 in the same way in which you give I_D^{int}.
- Concerning Eqs 5, 6: Can the conductances g_D->S and g_S->D actually differ? Isn't this just the passive conductance connecting the dendritic and axosomatic compartments which must be symmetric for physical reasons?
- On p 6, you give an expression for the mini rate µ including a parameter F, which I think is not given anywhere, and dividing by 250 at the very end of the expression. The corresponding expression on p 8694 in the original paper has F=400 and divides by 100, not 250. In the caption to Fig 11 you then have exactly the same numbers as in the original paper.
- Does the connectivity include periodic boundary conditions or do neurons at the edge of the layer receive input only from the "inside" of the layer? Could you clarify this?
- In Sec 4.2, it would benefit the reader if the text came before Table 2, so that one had an introduction to this large (and highly commendable!) table. I am still a bit confused about what you mean by "original value" and how you then obtained the often complex expression for the current model. Could you explain this more carefully? In the very first expression for the "present model", there seems to be a stray "B" at the end of the equation.
- In Figs 2–4, I find the use of light blue vs dark blue problematic, since it provides only low contrast. Furthermore, differences are only visible if the light blue line extends beyond the dark blue—in principle difference in K could be "hiding" in Fig 2A.
- I found your
environment.ymlfile overspecified. As it includes explicit build hashes, it will at best work for the operating system on which it was generated. It is in a way useful to precisely describe the software you used for your simulations, but for users trying to use your model, you should provide a minimal environment file. I could run your code after creating a new environment with justmamba create -n ra24 python numpy scipy matplotlib brian2. - Your code has some very long lines (over 200 characters), which makes reading difficult even with today's wider screens. I would strongly suggest to break lines around 120 characters at the most (the old 80 limit seems dated). You may also want to run a linter or code formatter over the code.
- I am confused by the code you use to introduce parameter variability, e.g. https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/8de8bb5df75d3f1eb4f7fcd8d7c5ed014e4beacc/Model/Cortical_layer.py#L117-L124. If I interpret the code correctly, you here aim to ensure that at most two parameters in a row receive a random adjustment with the same sign. This seems curious to me. Why not simply add random noise to the parameters as provided by the random number generator? Was something similar done in the original model? I did not see a comment on this in the manuscript.
- On this line: https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/8de8bb5df75d3f1eb4f7fcd8d7c5ed014e4beacc/Model/Figure_plotting.py#L215 and possibly others, the string must be a raw string so that the \ will be passed correctly to LaTeX.
- I find the connectivity specifications of the form https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/8de8bb5df75d3f1eb4f7fcd8d7c5ed014e4beacc/Model/Thalamo_cortical.py#L136 difficult to comprehend. At least at some point in the code there should be an explanatory comment on how these strings work.
- I could successfully reproduce Fig 5, but when I set
fig_number = "1"or 2, the simulation ran, but in the end reported only a mean firing rate of PY neurons of 17.678 Hz. What may have gone wrong here?
So much for now. Again, I am very impressed by the thorough work you have done in digging out details and parameters from old code and various papers and built a new implementation from this!
@heplesser Thansk you very much for your very detailed review. @apdavison Do you think you can finish your review by next week? @Mathilde-Reynes Can you adress @heplesser comments and post your answers/changes here?
my apologies for the delay. I'm travelling this week, but I will try to finish my review by March 21st
The authors have done an excellent job in replicating this important study. As is usually the case in computational neuroscience (and sadly I don't think this situation has improved in the last 20 years), the published paper does not contain all the information needed to replicate the models and simulations described, so it is to the credit of Bazhenov and colleagues that they shared code on ModelDB, even if it does not correspond exactly to the precise version used in the original study.
Even with the availability of the ModelDB version, Reynes and Aussel's detective work was a non-trivial task, and the results - understanding and explaining the differences uncovered - are in general presented very clearly. I particularly appreciated Table 2, and the way the Results section was structured by the claims/assertions from the original paper.
Reproducing the results
I would like to thank the authors for the modifications they made to the code, to make it easier to re-run simulations to match the figures.
I was able to run simulations for all figures, although the RNG seed available in Thalamo_cortical.py does not seem to be the one used for the manuscript. Figure 8, in particular, looks very different, presumably because of the highly stochastic and non-linear nature of up-state onset.
Since I don't use conda, I installed the project dependencies into a plain Python virtual environment (venv). For reference, the dependencies used were as follows (with Python 3.8.0):
Brian2==2.5.1
contourpy==1.1.1
cycler==0.12.1
Cython==3.0.11
fonttools==4.55.3
importlib_resources==6.4.5
Jinja2==3.1.5
kiwisolver==1.4.7
MarkupSafe==2.1.5
matplotlib==3.7.5
mpmath==1.3.0
numpy==1.24.4
packaging==24.2
pillow==10.4.0
pyparsing==3.1.4
python-dateutil==2.9.0.post0
scipy==1.10.1
six==1.17.0
sympy==1.13.3
zipp==3.20.2
Suggested modifications
- I disagree that the full network dynamics in the replication are "very similar" to those in the original paper (p14). In particular, there are large differences in the frequency of the up-states. I suggest using "qualitatively similar" or just "similar".
- Similarly, the up-states are not "minimally shorter", and the number of oscillations is not "negligibly" larger. On the basis of Figure 5, up-states in the replication seem to be smaller than in the original by a factor of 2 or more.
- When citing using numbers, please ensure the sentence remains readable. For example "recreating [1] behavior" should be something like "recreating the behavior shown in the original paper [1]", and "using [1] C++ code" should be "using the original C++ code [1]".
- The typesetting of equations and units was not standards-compliant, which hampered the readability of the document for me. Please modify your LaTeX source; you may find this note helpful.
- Table 2 should be introduced in the text before the table itself actually appears in the document.
- I think the colormaps of membrane potentials would be better as grayscale. This would make comparison with the original figures easier.
- There were many minor issues, such as typos. Rather than enumerate them all here, I have put my annotated copy of the manuscript online here - please let me know if you are unable to access it.
@apdavison Thansk you very much for the detail review. @Mathilde-Reynes You can now proceed with the corrections of your manuscript and please post your changes/answers here (maybe one post per review addressed). If you can do it in two weeks, that would be great.
Response to @heplesser
We would like to thank the reviewer for their thorough and constructive review, and sincerely appreciate the time invested in assessing our manuscript. Their feedback raised valid points that helped us improve the clarity and robustness of both the manuscript and the code. Below we address their comments and have indicated the associated changes in our point-by-point response. We also pushed the corrected manuscript on the git; with one version with changes in red and one ‘clean’ version.
First of all, I found it somewhat difficult to match figures in the manuscript to figures in the original paper. I would strongly suggest to mention in each figure caption the corresponding figure number in Bazhenov et al. I also wonder if it would not be easier to follow the reproduction step by step if figures were in the same order as in the original paper. I noticed that the manuscript contains excerpts of the figures in the original publication. While this definitely eases comparison, I wonder if this could raise licensing issues.
• We recognize that navigating the two papers simultaneously might be indeed laborious. We followed your advice of adding, in each figure caption, the corresponding figure number from Bazhenov et al. (2002). • We acknowledge your comment regarding the proposed rearrangement of the figures. However, we believe the order logically differs as the objectives of the two manuscripts are quite distinct. Bazhenov et al. (2002) aimed to present results and their interpretations in a narrative-driven manner. In our case, we focused on reproducing their claims with a more technical approach. Therefore, we structured our manuscript so that related claims are addressed within the same subsection, rather than switching back and forth as it would have been the case should we have followed Bazhenov et al. (2002) narrative. In that sense, we believe that adhering to the original paper's order would make the manuscript less readable and harder to follow. • As regards to potential licensing issues, as per our understanding of the website of original paper’s journal’s policy (https://www.jneurosci.org/content/rights-permissions), we believe we have permission to reuse its material, under the condition that we have to specify a full journal reference, link to the version of record, as well as a copyright statement (screen capture is provided below). We added all these elements as to respect JNeurosci guidelines to reuse their material. Modifications implemented on page 1.
Concerning the results, I am most concerned about Fig 5 in the present manuscript. The reproduction results show approximately 20 up and down states compared to only 11 in the original paper, i.e., twice the switching frequency than in the original, not just "minimally shorter up-states" as commented on p. 14. Also, in the lower part of Fig 5, the initiation of the transition looks much sharper in the original figure than in the reproduction, with very straight opening "flanks". Furthermore, the original figure only shows broad stripes spanning the entire PY population, while the reproduction figure also includes "tongues" that include only parts of the PY population. I believe this merits further investigation. In that context, I noticed that Fig 13 in the manuscript shows a similar difference in patterns, although here related to the presence and absence of thalamocortical connections. Perhaps something similar could be underlying the differences in Fig 5?
• We agree that our statement “our model’s up‐states appear minimally shorter, resulting in negligibly more oscillations during a typical 30‐second simulation” underestimated the reality and we appreciate you pinpointing it. The manuscript was modified in that sense. Modifications implemented on page 15. • Regarding the ‘sharper’ up-state initiation, we agree this is a significant difference between the present and the original model that we had not explicitly mentioned in the manuscript, and that we have corrected now. Modifications implemented on page 15. • We have been further investigating these two points, as there were the two key points we struggled to fully comprehend. You are also correct in noting that Fig 13 looks more similar to the original; we had the same thought. That being said, since the original figure explicitly states the number of TC and RE cells and displays their activity—activity that does not exist when disconnected from the cortical layer—we believe this suggests that thalamocortical connections were indeed present in the original figure. Reducing the weight of the specific connections from the thalamus to the cortical layer indeed slower oscillations, but not sufficiently enough to reproduce exactly the original results. For instance, running a simulation with only 10% of the original synaptic conductance between the thalamus and the cortical layer (g_syn_ampa_tcpy = 0.00001 msiemens; g_syn_ampa_tcin = 0.00001 msiemens) lead to 0.86 Hz cortical oscillations. This is still quite far from the 0.5Hz that can be observe in some figures of the original model. Another potential approach, clearly illustrated in Figure 10 of our model, is to reduce the conductance of miniPSPs. However, in a fully connected model, while decreasing the amplitude of miniPSPs slightly lowers the frequency of up-states, it also disrupts their organization. For this reason, we do not consider this parameter a viable option for modification. We also tested an increase in AMPA PY-IN conductance, as shown in Figure 12 of our manuscript, but this adjustment did not yield results closer to the original implementation. Since the slow accumulation of the I_KCa current in the dendritic compartment is believed to facilitate the transition from an up-state to a down-state, we attempted to reduce it in order to prolong up-states and thereby slow cortical oscillations. By decreasing the conductance of I_KCa (g_kca = 0.3msiemenscm**-2 * 0.2), we observed a clear increase in the duration of cortical up-states, making them more similar to those in the original model. However, the oscillations remained too fast, with shorter down-states overall. Nonetheless, these parameters cannot explain the discrepancy between our reproduction and the original model, as we use the same parameters when running the original C++ model and have verified that the calcium dynamics are consistent. As of now, we do not have an explanation for this difference between the original model and our implementation, still we acknowledge that this significant difference is important. We have added a supplementary figure showcasing the different tests we have performed to investigate that question. Modifications implemented on page 32
Concerning the sustained activity initially shown in Fig 6 of the manuscript, it would be interesting to see the corresponding plot for the current implementation. I wonder if initial conditions could play a role here?
• When running the original model on our system, we were able to reproduce the initial sustained activity. However, we were unable to replicate this result in our own model. We suspect this discrepancy may be due to one of the following reasons: (1) differences in the compilers used (the original model is compiled with gcc from C++ code, while our model uses Cython), (2) variations in the initialization of our parameters, or (3) a potential error in the reproduction process, likely stemming from a parameter mismatch (that could also be responsible for the faster oscillations observed in our reproduction). As we believe that the absence of this initial surge in activity aligns more closely with the biological reality, we are comfortable with keeping the model as it is on that point. We have completed Figure 6 with two new graphs showing the initial activity of the various implementations of the model. Modifications implemented on page 16
I would also like to point to a number of more detailed points.
- You frequently refer to the ModelDB version of the Bazhenov et al model, but do not provide a reference for it. Please add the ModelDB accession number and a link.
• Thank you for this remark. We have added the reference for the ModelDB, including an accession number and a link as per requested.
- I_{D->S} = -I_{S->D}, isn't it? It would be useful to point this out.
• As we will develop in your point 4. below, as the conductances g_D->S and g_S->D differ, I_{D->S} does not equal -I_{S->D}. Please refer to 4. for an explanation and the relevant changes made in the manuscript.
I think it would help the reader if you explicitly gave I_S^{int} somewhere above Eq 3 in the same way in which you give I_D^{int}.
• I_S^{int} is described in Eq 4 just below Eq 3.
- Concerning Eqs 5, 6: Can the conductances g_D->S and g_S->D actually differ? Isn't this just the passive conductance connecting the dendritic and axosomatic compartments which must be symmetric for physical reasons?
• In the present model g_D->S and g_S->D actually differ. Actually, g is defined as 1/(R) where R = rs_dend for the junction between the soma and the dendrite and R=rs_soma for the junction between the dendrites and the soma. r here is the resistance between compartments (10MOhm) and s_dend and s_soma are the area of the specific compartments. Therefore, the difference in compartment sizes is the main reason why g_D->S and g_S->D differ. Indeed, for PY neurons, the dendrite area is 165 times larger than the soma, and for IN neurons, it is 50 times bigger. If the same conductance density is used for both compartments, the dendritic compartment has a much higher absolute conductance due to its larger area, and therefore a lower input resistance. Basically this asymmetry translates the fact that somas are more reactive than dendrites. Indeed here, for the same injected current, the voltage change in the soma will be significatively higher than in the dendrite. This makes the soma more sensitive to dendritic inputs. In order to clarify this, we have added in the manuscript the equations for g_D->S and g_S->D. Modifications implemented page 5
- On p 6, you give an expression for the mini rate µ including a parameter F, which I think is not given anywhere, and dividing by 250 at the very end of the expression. The corresponding expression on p 8694 in the original paper has F=400 and divides by 100, not 250. In the caption to Fig 11 you then have exactly the same numbers as in the original paper.
• You're right to point out that equation. Not only is the "250" a typo, but the equation presented here uses an exponential arrival rate, whereas the rest of the paper consistently implements a logarithmic arrival rate. We have corrected this in the manuscript, ensuring the correct equation is shown and explicitly specifying the value of F. Modifications implemented page 6
- Does the connectivity include periodic boundary conditions or do neurons at the edge of the layer receive input only from the "inside" of the layer? Could you clarify this?
• In our current implementation, we have not set a specific connectivity boundary condition. Therefore, neurons at the edges of the layer receive input from fewer neurons compared to the neurons located further away from the edges. If we are not mistaken, this aligns with the second part of your question, stating that edge neurons “receive input only from the inside of the layer.” The original paper does not discuss boundary conditions. However, in the ModelDB implementation, a function is provided to define the desired boundary condition, allowing for options such as flow boundary condition, periodic boundary condition or no boundary condition (parameter BOUND and BOUNDcx, in function b()). The ModelDB version was configured with a parameter BOUND == 0 and BOUNDcx == 0, indicating that a flow boundary condition was used. However, since the manuscript does not explicitly describe it and as per our remarks on other parameters in this paper, we cannot positively confirm whether this was the condition applied in the final figures of the manuscript. Given that we did not observe a significant impact of the boundary condition on our results, we chose not to impose a periodic boundary condition. We have added detail on boundary condition in the table of the manuscript. Modifications implemented page 7 & 10
- In Sec 4.2, it would benefit the reader if the text came before Table 2, so that one had an introduction to this large (and highly commendable!) table. I am still a bit confused about what you mean by "original value" and how you then obtained the often complex expression for the current model. Could you explain this more carefully? In the very first expression for the "present model", there seems to be a stray "B" at the end of the equation.
• Thank you for the highly relevant remark. We put the text before the table in order to facilitate comprehension as you proposed. We added an extra sentence in this text to explain more thoroughly what we develop in the table below it. We removed the stray ‘B’ which was indeed a typo. Modifications implemented page 8
- In Figs 2–4, I find the use of light blue vs dark blue problematic, since it provides only low contrast. Furthermore, differences are only visible if the light blue line extends beyond the dark blue—in principle difference in K could be "hiding" in Fig 2A.
• We replotted the figure 2-4 with different colors in order to maximize the contrast. However, we would like to point out that whatever color is used to plot, if the plots align perfectly, one would only see one of the two colors. We played with different line thickness and transparency but they only slightly facilitate the readability of the figure. In order to explicitly confirm that there exist no different in K, we added a zoomed-in plot on the right of figure, that we think combined with our sentence “Should the equations and their integration be the same, a perfect match between the two will be observed. This is the case for $I_K$ but not $I_{Na}$.”; should make our statement clearer. Modifications implemented page 12
- I found your environment.yml file overspecified. As it includes explicit build hashes, it will at best work for the operating system on which it was generated. It is in a way useful to precisely describe the software you used for your simulations, but for users trying to use your model, you should provide a minimal environment file. I could run your code after creating a new environment with just mamba create -n ra24 python numpy scipy matplotlib brian2.
• Thank you for your feedback, indeed the environment file should be minimal. We have provided a new environment file that follows this guideline. Modifications implemented on the git
- Your code has some very long lines (over 200 characters), which makes reading difficult even with today's wider screens. I would strongly suggest to break lines around 120 characters at the most (the old 80 limit seems dated). You may also want to run a linter or code formatter over the code.
• The code is now formatted using the formatter Black.
- I am confused by the code you use to introduce parameter variability, e.g. https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/8de8bb5df75d3f1eb4f7fcd8d7c5ed014e4beacc/Model/Cortical_layer.py#L117-L124. If I interpret the code correctly, you here aim to ensure that at most two parameters in a row receive a random adjustment with the same sign. This seems curious to me. Why not simply add random noise to the parameters as provided by the random number generator? Was something similar done in the original model? I did not see a comment on this in the manuscript.
• You are correct we did not explain the parameter variability. We explicitly added a paragraph on this in the reviewed manuscript. This random adjustment is strictly reproducing the one in the original model. Modifications implemented page 6
- On this line: https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/8de8bb5df75d3f1eb4f7fcd8d7c5ed014e4beacc/Model/Figure_plotting.py#L215 and possibly others, the string must be a raw string so that the \ will be passed correctly to LaTeX.
• This is now corrected.
- I find the connectivity specifications of the form https://github.com/Mathilde-Reynes/ReynesAussel2024/blob/8de8bb5df75d3f1eb4f7fcd8d7c5ed014e4beacc/Model/Thalamo_cortical.py#L136 difficult to comprehend. At least at some point in the code there should be an explanatory comment on how these strings work.
• As proposed, we have added explanatory comments for the connectivity specifications on all files with synapses (Cortical_layer.py ; Thalamus.py and Thalamo_cortical.py).
- I could successfully reproduce Fig 5, but when I set fig_number = "1" or 2, the simulation ran, but in the end reported only a mean firing rate of PY neurons of 17.678 Hz. What may have gone wrong here?
• Basically, as of this submission, whenever a user provided the function with a number that is not an actual figure that can be plotted (such as 1 or 2), a 10-seconds simulation with the base parameters (the ones of Figure 5) is ran. An additional warning message was added in the code to make the situation clearer. • The mean firing rate of pyramidal neurons is computed by dividing the total number of recorded spikes (len(R2_PYs.t)) by the number of PY neurons (N_PY) and the total simulation runtime (runtime). The mean firing rate of PY neurons is different from the slow oscillations frequency. Indeed, during an upstate, PY neurons spike multiple times (could be 15 to 25 times), thus leading to the 17Hz mean firing rate that you reported. According to us, nothing went wrong in the simulation you ran.
We would like to thank the reviewer again for the constructive comments and remarks, and are glad you can appreciate the effort we put in reproducing this model.
Response to @apdavison
We sincerely thank the reviewer for their initial review on January 10th, which prompted us to refine our code, making it more user-friendly and facilitating easier reproduction of results. We would like also to thank the reviewer for their latest review, and sincerely appreciate the time invested in assessing our manuscript. Their feedback raised valid points that helped us improve the clarity, readability and robustness of both the manuscript and the code. Below we address their comments and have indicated the associated changes in our point-by-point response. We also pushed the corrected manuscript on the git; with one version with changes in red and one ‘clean’ version.
[…] I was able to run simulations for all figures, although the RNG seed available in Thalamo_cortical.py does not seem to be the one used for the manuscript. Figure 8, in particular, looks very different, presumably because of the highly stochastic and non-linear nature of up-state onset.
• Indeed, the seed used has a significant impact on the spatiotemporal dynamics of cortical oscillations, as the miniature excitatory postsynaptic potentials (modeled as a Poisson process) influence the onset of up-states in the oscillations. Additionally, variability in some conductances may also contribute to differences, and we have added a paragraph addressing this variability in the revised manuscript as we did not mention it in the initial submission and as reviewer 1 suggested it (see modifications on page 6). However, regardless of the seed, the qualitative behavior of the model remains consistent. The RNG seed available in Thalamo_cortical.py was in fact not the one used for this specific figure of the manuscript and we apologize for this. Running the model with the conda env on the git, with the RNG seed currently available, we were not able to reproduce the figure you sent (panel A of the figure below). This same result was obtained on a Linux and a Windows machine. In order to investigate, we have installed your venv on my Windows machine and ran the model on this environment. I still obtain the same results (panel B of the figure below). Thus we could not reproduce the issue with the elements mentioned. We would like to emphasize that some compiler optimizations methods can vary floating-point precision and may lead to such differences. For example, should you allow the compiler to reorder your arithmetic operations and optimize to get faster math instructions, results can vary in complex non-linear models. Therefore, even if we ensure the seed is fixed, we cannot guarantee exact reproducibility of results—especially in the longer term or across different environments. Modifications implemented page 19
Since I don't use conda, I installed the project dependencies into a plain Python virtual environment (venv). For reference, the dependencies used were as follows (with Python 3.8.0) […].
• Thank you for providing us with your verv specifications. We have added, in the read me of the git, a note for users that would like to run the model without conda. This should help reproducibility as well. Modifications implemented on the git
Suggested modifications:
- I disagree that the full network dynamics in the replication are "very similar" to those in the original paper (p14). In particular, there are large differences in the frequency of the up-states. I suggest using "qualitatively similar" or just "similar".
- Similarly, the up-states are not "minimally shorter", and the number of oscillations is not "negligibly" larger. On the basis of Figure 5, up-states in the replication seem to be smaller than in the original by a factor of 2 or more.
• Thank you for these two very relevant remarks (1 and 2). We agree that our statements “our model’s up‐states appear minimally shorter, resulting in negligibly more oscillations during a typical 30‐second simulation” and “the reproduced model displays dynamics very similar to those described by Bazhenov et al.” totally underestimated the reality of our results and we appreciate you pinpointing it. The manuscript was modified in that sense.
Modifications implemented page 15
- When citing using numbers, please ensure the sentence remains readable. For example "recreating [1] behavior" should be something like "recreating the behavior shown in the original paper [1]", and "using [1] C++ code" should be "using the original C++ code [1]".
• Changes have been made to the manuscript to improve its fluidity and readability. Given the extensive modifications across multiple pages, we invite you to consult the ‘reviewed’ version of the manuscript, where all changes are highlighted in red. Basically, we have incorporated additional context whenever references are mentioned in the text, as per your recommendation. These changes especially apply for citations of the original paper and ModelDB.
- The typesetting of equations and units was not standards-compliant, which hampered the readability of the document for me. Please modify your LaTeX source; you may find this note helpful.
• Thank you for providing the link to this, extremely useful, note. We applied your recommendations to all our equations and units so that they are standards-compliant. We did not highlighted them in red as basically all the equations and units were changed.
- Table 2 should be introduced in the text before the table itself actually appears in the document.
• Thank you for the highly relevant remark. We put the text before the table in order to facilitate comprehension as you proposed. We added an extra sentence in this text to explain more thoroughly what we develop in the table below it. Modifications implemented page 8
- I think the colormaps of membrane potentials would be better as grayscale. This would make comparison with the original figures easier.
• The colormaps were replotted as grayscale. Modifications implemented page 22 & 15
- There were many minor issues, such as typos. Rather than enumerate them all here, I have put my annotated copy of the manuscript online here - please let me know if you are unable to access it.
• We were able to access the annotated copy of the manuscript and made all the changes you noted. We thank you for providing feedback on this aspect of our work as well, which is definitely as important as the scientific content itself.
We would like to thank the reviewer again for the very helpful comments and remarks, and are glad you can appreciate the effort we put in reproducing this model.
@Mathilde-Reynes Thanks for the detailed answer. @heplesser @apdavison Are you satisfied with the answers and corrections?
I am satisfied with the answers and corrections.
I do have one further query for the authors: now that the colormaps are in grayscale, it seems as though there is a systematic difference between the original and replicated versions (e.g., Figure 5): specifically, both up- and down-states appear substantially darker in the replications. Does this represent an actual difference in the activity levels, or is it just that the exact gray-level map is not known for the original figures?
@apdavison, Thank you for the very relevant remark. To better understand the situation, let’s refer to Figure 7 of our manuscript:
In the bottom panel of Figure 7, we see that the resting membrane potential appears to be consistent across all cell types between the original model and our implementation, even though the lines in the original paper’s plot are missing (-68 mV for PY cells, -66 mV for IN cells, -73 mV for RE, and -72 mV for TC). According to this, the voltage range during the down state is similar in both our implementation and the original model. That being said, the colors in Figure 5 during the down states should therefore match, but they do not. This suggests that either the colormap or the cbar used in both plots are not identical, as you suggested. We did attempt to find a more suitable grayscale to align with the original figures, but without success.
Your question raised an important point: why does the down state for pyramidal cells appear lighter (whiter) than that of thalamic cells in the original manuscript, even though its value is higher (-68 mV vs. -72 mV)? The down-state value of -72 mV for TC cells should logically appear lighter than the -68 mV value for PY cells, yet this is not the case. This likely indicates that two different colorbars were used across the same figure, which would also explain why we couldn’t find an appropriate color scale to replicate it. To illustrate this, we’ve re-plotted Figure 5 with two distinct colorbars for pyramidal and thalamic relay cells, using arbitrary maximum and minimum values for the colorbars. This version aligns much better with the original manuscript’s figure though it is less intuitive for readers who can no longer easily compare the two cell types. Therefore, we do not believe using two separate color scales is the best approach.
Since we think this explanation will be helpful to the reader, we have added a brief sentence in the manuscript to summarize the point discussed above: “Nonetheless, our figure still differs because the original figure likely used two different colorbars for pyramidal and thalamic cells, resulting in internal mismatches that naturally explain the differences with our figure which uses a consistent colorbar across all plots. This claims is supported by the cellular activity shown in Figure 7.” Modifications implemented page 15
@apdavison Thansk, do that means the article is now ok for you? @heplesser Are you satisfied with answers and corrections?