[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.