[Re] Cortical network effects of subthalamic deep brain stimulation in a thalamo-cortical microcircuit model

[celinesoeiro]

Original article: AmirAli Farokhniaee and Madeleine M Lowery 2021 'Cortical network effects of subthalamic deep brain stimulation in a thalamo-cortical microcircuit model', Journal of Neural Engineering, DOI: 10.1088/1741-2552/abee50

PDF URL: https://github.com/celinesoeiro/TCM-model/blob/main/Neural_network_model_consisting_of_cortex_and_thalamus_to_simulate_effects_of_deep_brain_stimulation_in_Parkinson_s_disease.pdf
Metadata URL: https://github.com/celinesoeiro/TCM-model/blob/main/metadata.tex
Code URL: https://github.com/celinesoeiro/TCM-model/tree/main

Scientific domain: Computational Neuroscience
Programming language: Python
Suggested editor:

[rougier]

Thanks for your submission. An editor will be assigned soon.

[rougier]

@benoit-girard Can you edit this submission (else I can do it)

[benoit-girard]

Maybe you could act as an editor, so I could be one of the reviewers...

[rougier]

Ok!

[rougier]

@NikVard Can you review this submission?

[NikVard]

@rougier I'm on it!

[rougier]

Perfect ! @benoit-girard @NikVard you can start the review.

[NikVard]

I read the manuscript and the code and tried running the "DBS.py" file to generate some results. In general, the README file is well-written and along with the manuscript it guides the user on how to run the code for the first time. However, I have a few comments:

Repository

While the README file does a good job of detailing the steps to reproduce the model, I find that it is missing a few things. It would be helpful to include installation instructions in one of the following ways:

1. Provide a requirements.txt file with the necessary modules to install using pip.
2. Provide a conda .yml file, if conda was used, with the complete environment (including the Python version).
3. If using pyenv or another method, provide the Python version and a requirements.txt file to replicate your setup.

This way you guarantee that your code will be reproducible by others with minimal hassle. I had to manually go over and install the dependencies as errors appeared when running the code, which can be avoided by following of the above suggestions.

Code

I can follow the code, but when I run the DBS.py file, there is an error in line 322 in the function plot_BP_filter. From what I understand, you have made changes to the function's definition and it now takes different/fewer arguments. I believe you redefined the function to take the DBS frequency, but omitted the bandpass frequencies. Please address the issue.

There are some discrepancies between the parameters reported in the paper and the provided files. For instance, in the paper you mention that you have 540 neurons, however, by looking at the provided code and the file tcm_params.pym the default values to run the model with is using 1/10th of the reported neurons (54 total). Moreover, the default simulation time in your parameters file is only 3 seconds. As a suggestion, it would be nice to set the parameters in your code to correspond to the simulations that you report in your paper to avoid confusion when trying to replicate the results.

One more comment regarding the code for your spectal analysis. Why is the PSD function (in model_functions.py) defined this way?

def PSD(signal, fs):
    (f, S) = welch(signal, fs, nperseg=10*1024)
    
    return f, S

You have hardcoded the parameter nperseg to 10*1024. If you were aiming to obtain an estimate using Welch's method, this will not achieve what you wanted. With the above code, the parameter nfft will be set to the value 10*1024 and you will have a single segment for your estimate. A better way to do this is to actually use Welch's method as originally intended to obtain the averaged periodograms. If you do so, then 10*1024 is a very large window for obtaining an averaged PSD estimate. Assuming you are using the default time step for your simulations at 10 ms (taken from the parameters file dt = 100/ms # time step of 10 ms), a window with a length of 1 second (to achieve a frequency resolution if 1 Hz in your PSDs) would only be 100 sample points. To smoothen out your PSD estimates, you should set the parameter noverlap to get overlapping windows.

I would also like to suggest one more modification. The model currently plots figures as it runs, which are opened and block execution but never saved. Instead, opt for saving the figures in a dedicated "results" directory (i.e. using the code fig.savefig("out.png")). This way, the figures after each simulation are stored on-disk and are easier to inspect.

Paper

• First paragraph of "Methods": "Neurons of the S, M, S, and TC..." -> "Neurons of the S, M, D, and TC..."
• Table 1: you mention that the model contains 540 neurons in total, however, from the table, we see there are 600 neurons. From the code, I believe your TR should have 40 inhibitory neurons.
• In the last paragraph of page 3, you mention that the noise term was "scaled so that the mean firing rates of the neurons are compatible with experimental recordings;". Which recordings are you referring to? You can add the experimental reference from the original work here (I believe it's this one: "Li Q, Ke Y, Chan D C W, Qian Z M, Yung K K L, Ko H, Arbuthnott G and Yung W H 2012 Therapeutic deep brain stimulation in parkinsonian rats directly influences motor cortex Neuron 76 1030–41"
• In the first paragraph of page 4, you mention that "The scale of this threshold noise was set to one-third of the scale of ξ(t)". What does it mean "the scale of this threshold noise was set to one-third of ξ(t)"? Was it arbitrarily selected? In the original model, the mean of the threshold noise was set to half of the additive white noise.
• The paragraph just below Table 2 is quite unclear. I think it would benefit greatly from rephrasing, especially the part starting with "for the excitatory neurons..."
• In Section 2.3 "Thalamo-cortical microcircuit model", you mention that "The values of synaptic strengths were adjusted to produce a simulated LFP in cortical layer D that resembles the LFP recorded under normal conditions". Could you specify how the original values were adjusted? How "close" was the simulated LFP to the recorded LFP and how was the similarity estimated?
• In Section 2.4 "Stimulation Protocol", you describe the "DBS off" vs the "DBS on" conditions in a simulation. In the original work, stimulation was applied for a total time T. Here, you stimulate for a time of T/3, which is 5 seconds for a simulation time of T=15 s. In the provided code of the original paper, they used a simulation time of 10 (+1) seconds to reach a steady state. Was it consistent?
• It would have also been interesting to replicate the original authors' reported increase in beta-band power in response to low-frequency DBS (i.e. <40 Hz) and examine the results.

Results

Figure 2: It would be more convienient if you could have a direct comparison with the matrices of the original work. There are also some minor discrepancies between the two, which you could address later in the "Discussion" section.

Figure 3: This is an important figure and it really deserves more "real-estate". Don't hesitate to allocate more space to it!

Figure 6: It is really difficult to see the differences between all conditions in this case. I am missing the key points of this figure. If you wanted to highlight the same frequencies as Figure 4, maybe you could indicate it on the figure directly, or zoom in the frequency range of interest.

Minor language notes

• Above equation 1: "... its dynamics is given..." -> "... its dynamics are given..."
• Change "one third" to "one-third"
• In section 3 "Implementation details": "separeted file" -> separate, also not hyphenated properly. In LaTeX, you can manually avoid it in such cases by explicitly stating the hyphenation (i.e. se-pa-ra-te).
• Section 3.3 "split and conquer" -> "divide and conquer"
• After section 3.2, your paragraphs are no longer properly spaced. I believe you have a misaligned bracket in your Tex file somewhere in the beginning of 3.2
• There is an errand "Figure" in Section 4 "Results".

Conclusion

Overall, the authors have translated the original model from MATLAB to Python. The model is nicely documented, however, some changes need to be made to make it run without issues. In it's current state, the file DBS.py produces an error at line 322 and does not complete its execution.

Using their reproduced model, the authors have also replicated some of the main findings of the work by Farokhniaee et al. (2021). Notably, they have shown that high-frequency DBS at 130 Hz can reduce beta synchronization in parkinsonian states, as measured from the simulated cortical LFP.

I am attaching the PDF with comments and corrections for convenience. Minor comments: Some sentences and typos to be addressed (marked in red).

There is an opportunity to model the effects of low-frequency DBS (<40 Hz) and verify the claims of the original authors and it is something that can be added easily to the current work. I believe that by doing so, the authours would further improve the scientific utility of the current work.

Neural_network_model_consisting_of_cortex_and_thalamus_to_simulate_effects_of_deep_brain_stimulation_in_Parkinson_s_disease.pdf

[celinesoeiro]

Thank you for taking the time to review our work, @NikVard. I will carefully review your comments and address them.

[celinesoeiro]

@NikVard, we thank the reviewer for the time and effort devoted to our manuscript. We found his comments and suggestions helpful in improving the paper. We revised the manuscript, answering all points raised by the reviewer. Below, we reproduce the reviewer’s comments in boldface and our replies are in regular font. The PDF file was also updated and can be found in the in this link.

In Section 2.4 "Stimulation Protocol", you describe the "DBS off" vs the "DBS on" conditions in a simulation. In the original work, stimulation was applied for a total time T. Here, you stimulate for a time of T/3, which is 5 seconds for a simulation time of T=15 s. In the provided code of the original paper, they used a simulation time of 10 (+1) seconds to reach a steady state. Was it consistent?

• In the original code (https://github.com/ModelDBRepository/266941/blob/main/ThalamoCortical_Microcircuit_PD_DBS_ModelDB.m), in line 190, one can see that when DBS is applied the authors divide the total time T in 3 equal parts. Regarding the +1 second to reach a steady state, we applied it at first but the results with or without it were very similar. Therefore we believe it is safe to say the results are consistent.

It would have also been interesting to replicate the original authors' reported increase in beta-band power in response to low-frequency DBS (i.e. <40 Hz) and examine the results.

• Done, with a frequency of 25 hz.

The model currently plots figures as it runs, which are opened and block execution but never saved. Instead, opt for saving the figures in a dedicated "results" directory (i.e. using the code fig.savefig("out.png”)).

• Done

Updates in the PSD function to smoothen out your PSD estimates, you should set the parameter noverlap to get overlapping windows. Remove the nperseg.

• Done

Include installation instructions in the README Provide a conda .yml file, if conda was used, with the complete environment (including the Python version). Provide a requirements.txt file with the necessary modules to install using pip.

• Done

Update the parameters in tcm_params.py to match the ones in the paper

• Done

I can follow the code, but when I run the DBS.py file, there is an error in line 322 in the function plot_BP_filter

• Thank you for pointing that out. We believe now you can run the code without this error.

In Section 2.3 "Thalamo-cortical microcircuit model", you mention that "The values of synaptic strengths were adjusted to produce a simulated LFP in cortical layer D that resembles the LFP recorded under normal conditions". Could you specify how the original values were adjusted? How "close" was the simulated LFP to the recorded LFP and how was the similarity estimated?

• In the original paper, the authors used Bayesian model comparison (BCM) to adjust the strengths for both normal and Parkinsonian conditions. For the Parkinsonian condition, they further adjusted the inhibitory strength from CI to D to increase the level of synchronous activity in layer D neurons. We added this explanation to the manuscript together with the cited references in the original article. In our reimplementation, we used the synaptic strengths given in the original code.

In the first paragraph of page 4, you mention that "The scale of this threshold noise was set to one-third of the scale of ξ(t)". What does it mean "the scale of this threshold noise was set to one-third of ξ(t)"? Was it arbitrarily selected? In the original model, the mean of the threshold noise was set to half of the additive white noise.

• This is due to a discrepancy between the original code and the original article. While in the original article the authors write “half of the additive white noise”, in the code they use one-third of the additive white noise (see lines 147 and 148 of the original code in ModelDB). We used the information given in the code because this provided a better fit to the plots in the original article. We added a comment about this discrepancy to the manuscript.

In the last paragraph of page 3, you mention that the noise term was "scaled so that the mean firing rates of the neurons are compatible with experimental recordings;". Which recordings are you referring to? You can add the experimental reference from the original work here (I believe it's this one: "Li Q, Ke Y, Chan D C W, Qian Z M, Yung K K L, Ko H, Arbuthnott G and Yung W H 2012 Therapeutic deep brain stimulation in parkinsonian rats directly influences motor cortex Neuron 76 1030–41”.

• Done

The paragraph just below Table 2 is quite unclear. I think it would benefit greatly from rephrasing, especially the part starting with "for the excitatory neurons..."

• We agree it was quite confuse. We have rephrased the paragraph, and hope it is clearer now.

Table 1: you mention that the model contains 540 neurons in total, however, from the table, we see there are 600 neurons. From the code, I believe your TR should have 40 inhibitory neurons.

• Thanks for pointing out this mistake to us. We have corrected the number of TR neurons in Table 1.

First paragraph of "Methods": "Neurons of the S, M, S, and TC..." -> "Neurons of the S, M, D, and TC..."

• Done

[rougier]

@celinesoeiro Thanks a lot for these corrections. @benoit-girard Can you upload your review in a few days?

[NikVard]

@celinesoeiro Thank you for addressing the earlier comments and for the thorough description of the steps you took in response. After your edits, I can now run the code without issues and I will make some more comments that could be addressed.

Code

1. First, thank you for addressing my comment and providing a requirements file and instructions to install the necessary python modules. I would like to point out that the conda environment you provided is OSX specific, and conda failed to find the relevant packages on other systems (linux and windows). You can follow the instructions from this link to make a non-platform-specific .yaml file that everyone can run if you wish.

2. As this is a computational work, could you consider specifying explicitly in the text the simulation parameters that you used to produce the figures? I believe there may be some errors in the code that need to be addressed, which I will detail in the following points. It is also good practice to include the parameters in the text to encourage the reproducibility of your work by others who read the paper alone. Could you kindly specify the simulation parameters in a table?

3. I noticed something strange in the declaration of your parameters. In the tcm_params.py file, the default time step is initialized to 100/ms; that is 100ms, not 10ms as stated in the comments. Can you clarify if you ran all simulations with a time step of 10ms or 100ms?

4. Unless I am mistaken, there are also some follow-up issues with the initialization of other variables. Following point 3 from above, in the tcm_params.py file you define dt as 100/ms, which, as stated above, evaluates to 0.1 seconds. Then, you proceed to define your sampling frequency as samp_freq = int(ms/dt) # sampling frequency in Hz, which gives a value of 10000. This corresponds to a 10kHz sampling frequency (and a dt of 0.1 ms). This is inconsistent with what you probably intended to do, as with a dt of 0.1 seconds (100ms equivalently) you would only have a sampling frequency of 10Hz. This is evident also by evaluating the expression you wrote in terms of units, as ms / (100/ms) gives units of time squared, not Hz. The proper definition for your sampling frequency should be int(1/dt).

5. Following up from point 4, when defining your simulation steps you write sim_steps = int(T/dt) # number of simulation steps. This evaluates to 150000, as you divide T (which you have initialized earlier to 15000, and given in ms) with dt (which is 0.1, and I assume is given in seconds). In reality, with a dt of 0.1s (100ms) you only have 10 steps per second, which would evaluate to sim_steps = 150 (10 steps per second and 15 seconds of simulation time). Even if we assume a 10kHz sampling rate and ignore the above, if we follow the initialization of the sim_steps to 150000, then a single step does not correspond to a ms but 0.1ms, making all of the delays right below that (all the td_... variables) inconsistent. Overall, your comments and the declared variables create some confusion and could be improved.

6. Moreover, some of the figures that are produced by running the DBS.py file are strange. For example, in the figure LFP - 80Hz.png the x-axis I believe gives simulation steps instead of time. Furthermore, the figure with the LFP looks a bit strange (see the following figure). . This figure was generated by running the DBS.py with the default parameters, and a DBS frequency of 150HZ (default value). Could you modify the code to be consistent with the rest of the plots?

7. On another note, the simulations take an incredibly long time of around 2 hours to complete on a Linux computer with a 10th gen 16-core Intel i7 processor with 64GB of RAM. Can you confirm that you are indeed running each simulation with the intended parameters? Also, if possible, could you specify the time it takes to run a complete simulation on your end, along with the technical specs of the system you used?

8. I noticed that in your results you are generating plots with the PSDs for DBS applied at 130Hz and 180Hz, however, the x-axis is clipped at 100Hz, thus I cannot verify the results. If you consider the results useful, please modify the x-axis so that we can see the results.

9. Finally, in the model_functions.py file, the function PSD that computes the periodogram using Welch's method (before plotting it using the plot_PSD_DBS function) and is being called during the DBS.py execution (lines 323, 328) also has hard-coded values for the nperseg parameter (line 222, nperseg=10*1024). Taken directly from the file:

def PSD(signal, fs):
    (f, S) = welch(signal, fs, nperseg=10*1024)

    return f, S

This should also be fixed as for the plotting functions.

Minor language notes

I noticed that in your PDF, the spacing between paragraphs is not always consistent. Specifically, almost all paragraphs do not have vertical spacing, but the first paragraph in page 2 is different, as is the last paragraph of page 9. Is this done on purpose or is it an errand LaTeX error?

I am also listing some typos that I came across for your convenience. I would advise to proofread your text using a syntax / grammar checker (Grammarly, for example), to catch any others that I might have missed.

1. In section 3 Implementation details - "...separete file" -> "separate file".

2. In section 3.3 The TCM model - "...to devide and conquer..." -> "divide"

Conclusion

I would like to thank the authors, as they took into account all earlier comments and made significant modifications to their work, implementing a range of suggestions and explaining their approach in a detailed and convincing manner. Importantly, they successfully replicated the claims of the original work regarding the effects of low-frequency DBS (<40 Hz) by using a stimulation protocol at 25Hz. Overall, the work represents a thorough replication of an interesting model to investigate DBS for PD.

[rougier]

@NikVard Thansk fro your second review
@benoit-girard Can you still do the review?