The improved thalamo-cortical spiking network model of deep brain stimulation

1 Introduction

Parkinson's disease (PD) is a neurodegenerative disease associated with the loss of dopamine-producing neurons in the brain and is characterized by abnormal neural firing within the cortex and basal ganglia regions of the brain. Deep Brain Stimulation (DBS) is an established therapy for PD and is used clinically to relieve motor symptoms. The majority of research to-date on the symptoms of PD and their control by DBS has focused on the altered basal ganglia function (Wiecki and Frank, 2010; Caligiore et al., 2016). However, recent studies suggest that the key neural activity changes may be mediated and/or driven by the motor cortex (Lindenbach and Bishop, 2013). There are compelling clinical reasons to consider the importance of the motor cortex in the generation of symptoms and cortical stimulation as a potential therapeutic in PD (Caligiore et al., 2016; Arbuthnott and Garcia-Munoz, 2017). This is critical in understanding the role of antidromic activation of the cortex and alterations in cortical firing patterns in the therapeutic efficacy of DBS (Anderson et al., 2018). There has been a vast investigation in understanding DBS mechanisms of action during the last two decades. Recently, a thalamo-cortical microcircuit model of PD and DBS (TCM) was developed and introduced to the computational neuroscience community interested in network models of DBS (Farokhniaee and Lowery, 2020, 2021). Since the TCM model exhibited known neurophysiological features in PD and manifested DBS cortical effects in animal and human models despite its simplicity, it received a particular attention. Therefore, in this article, we fix a typo in the original article and suggest an improved version of the model to ensure the correct integration of the noise term. We also provide parallel computing and running on graphical processing unit to increase the computational efficiency of the model, useful for simulations with a higher number of neurons.

2 The improved model

2.1 TCM model structure

Certain alterations of synaptic weights within and between thalamus and cortex in a neural mass model of thalamocortex led to elevated beta power (~13–30 Hz oscillations) in the rats motor cortex (Reis et al., 2019), a well-known neurophysiological activity in parkinsonism. On the other hand, a rat model of PD showed the exaggerated synchronized patterns of spiking neurons in addition to the exaggerated beta power (Li et al., 2012). Inspired by these studies, the spiking neuronal network of thalamocortex in PD-like conditions was introduced (Farokhniaee and Lowery, 2020) and developed as a network model of DBS, known as TCM (Farokhniaee and Lowery, 2021). The developers of TCM introduced an important underlying mechanism of action of DBS at high frequencies that is the synaptic suppression due to short term synaptic plasticity (Farokhniaee and McIntyre, 2019) so that TCM served as a biophysically-realistic network model of DBS. As such, TCM exhibited known network effects of DBS that includes elevated beta power, exaggerated synchronized pattern of neuronal spikes, and formation of neuronal clusters such as excited and inhibited ones and optimized intensity of DBS-induced electric field to deliver the most suppression of the elevated beta power.

TCM contains 540 subthreshold noise-driven spiking neurons that obey Izhikevich neuronal dynamics (Izhikevich, 2003) inherently, connected via Tsodyks-Markram synapses (Tsodyks and Markram, 1997). The excitatory populations in the primary motor cortex were distributed into three layers: supragranular or surface (S, 100 neurons), granular or middle (M, 100 neurons), and infragranular or deep (D, 100 neurons) with a shared population of cortical inhibitory neurons (CI, 100 neurons). Thalamo-cortical relay nucleus (TCR, 100 neurons) and thalamic reticular nucleus (TRN, 40 neurons) form the excitatory and inhibitory populations of thalamus, respectively. The distribution of neurons in each substructure of the model along with complete neuron and synapse parameters are already presented in (Farokhniaee and Lowery, 2021).

The network dynamics of the model are described by the following set of equations (Equation 1). In the original paper (see Equation 2 in Farokhniaee and Lowery, 2021), i and j indices were also placed in the summation that includes postsynaptic currents (PSCs) and here we provide the correct formula where the summation must be over i' and j' only to obey the correct mathematical logic:

$\begin{array}{ll}{\dot{v}}_{ij}=0.04v_{ij}^2+5v_{ij}-u_{ij}+140+I_{ij}+{\displaystyle \sum }_{j\text{'}=1}^6{\displaystyle \sum }_{i\text{'}=1}^{N_j}{\omega }_{i\text{'}j\text{'}},{}_{ij}PS{C}_{i\text{'}j\text{'}}& (1)\end{array}$

$\begin{array}{r}\left(t-{\text{Δ}}_{j,j\text{'}}\right)+{\displaystyle \sum _{k,t_k}{u}_{jk}\delta (t-t_k)}+\xi \left(t\right)+I_{dbs}{\delta }_{j3}\\ {\dot{u}}_{ij}=a_{ij}(b_{ij}{v}_{ij}-u_{ij})\\ \text{if }{v}_{ij}\ge v_{p_{ij}}+\zeta (t) \text{then} v_{ij}\to c_{ij} \text{and} u_{ij}\to u_{ij}+d_{ij}\end{array}$

For each structure of the TCM model, j = 1, 2, 3, 4, 5, 6 corresponds to the structures S, M, D, CI, TRN and TCR, respectively. For each neuron i = 1, 2, …, N_j in structure j (where N_j is the total number of neurons in layer j) v_ij is the membrane voltage and u_ij represents the membrane recovery variable, where a_ij, b_ij, c_ij, and d_ij are each neuron parameters, with random changes that provide non-identical neurons in the network. I_ij are the bias currents with ξ(t) and ζ(t) as white gaussian noises. I_dbs is the DBS-induced intracellular transmembrane current that is added only to layer D (j = 3) by means of the Kronecker delta function. The synaptic connections deliver the PSCs to the neurons by means of a weighting matrix with ω_ij elements. $\sum _{k,t_k}{u}_{jk}\delta (t-t_k)$presents poissonian background noise, with inputs occurring at time t_k in each structure. The PSCs are the solutions I in Tsodyks-Markram dynamics given by Equation 2.

$\begin{array}{ll}\dot{u}=-\frac{u}{{\tau }_f}+U(1-u^{-})\delta (t-t_s-\text{Δ})& (2)\end{array}$

$\begin{array}{r}ẋ=-\frac{1-x}{{\tau }_d}-u^{+}{x}^{-}\delta (t-t_s-\text{Δ})\\ İ=-\frac{I}{{\tau }_s}+A{u}^{+}{x}^{-}\delta (t-t_s-\text{Δ}),\end{array}$

where x represents the fraction of available neurotransmitters after synaptic transmission and u is the fraction of available neurotransmitter resources ready to be used. t_s is the spike time, δ is the Dirac delta function, U is the increment of u produced by an incoming spike, τ_f, τ_d, .and τ_s are the decay (or recovery) time constant of variable u, x and I, respectively. A is the absolute synaptic response.

2.2 Simulation

In this opinion article we suggest a computational code that ensures embedment of all to all random connectivity by completely rewriting the code in a matrix-based computation at each time-step, that also increased the speed and efficiency of the algorithm (see for example updateTimeStep.m file in the model repository). We also identified a numerical issue in integrating the noise term in the original article computer code that is now edited and improvised by removing the multiplication of the noise term with each time step, though it was not creating a significant issue that lead to false results in the original study since this issue manifests during longer run times than the published results. As another model enhancement, when simulating the model using CPU, the parallel algorithm is utilized in updating each timestep through the whole network (multicore and multithread computing). In addition to solving by CPU power, we developed this parallelization on graphical processing unit (GPU, available only on Windows machine), useful in reducing the run time for simulations that include a higher number of neurons and populations. The simulations that run on GPU show their strength when the number of neurons is much higher than a several hundreds. Setting a desirable subpopulation to 0 will lead to deactivation of that subpopulation from the whole network.

The results of the model simulations for 12 seconds are shown in Figure 1 as raster plots. We ran the model on both Windows and Macintosh machines with success. We neglected the poissonian background noise in running the simulations. The white Gaussian noise had the mean of zero and the standard deviation of 0.5. The threshold noise was set to have the mean of zero and the standard deviation of 0.1. Figure 1A shows the effect of DBS at 20 Hz where it is seen that the synchronized activity is present after the DBS onset (second 6) whereas Figure 1B is the same network but during 130 Hz DBS in which the synchronized activity of the populations is diminished.

Figure 1

Figure 1. The raster plots of the TCM network for 12 s are presented here with the DBS onset at second 6 for different DBS frequencies (A) 20 Hz and (B) 130 Hz. The desynchronization effect of 130 Hz DBS in contrast to 20 Hz DBS is visible. The effect takes place via propagation of the signals in the network through synaptic connections.

3 Discussion

In this article we provided an opinion about the already published network model of DBS, called TCM that has shown promising network effects in both animal (rat and monkey) and human models during rest and no visual stimuli. We suggested further developments of the TCM model by fixing a typo of the equations in the original paper, improving integration of the noise term in the original model computer code and ensuring all-to-all random network connections. The typo in the formula of the original study did not affect the original simulation results and was present only in the main text and not in the computational code of the model. We also enhanced the model computational efficiency, i.e., higher simulation speed, and inclusion of parallel computing and GPU power. DBS at high frequencies (>100 Hz) particularly around 130 Hz has been used clinically for more than 30 years in alleviating PD symptoms. Its neurophysiological effects are known mostly with suppression of exaggerated beta power and desynchronization of neuronal network firing of action potentials that is elevated during Parkinsonism. Application of DBS at 130 Hz on TCM resulted in desynchronization of the spikes as shown in the raster plot, Figure 1B, in contrast to 20 Hz DBS, Figure 1A. One of the downfalls of low frequency DBS is the resonance effects and production of unwanted harmonics. These effects plus the strengths of the TCM model in explaining clinical implications of DBS at high frequencies have been well discussed already in the original article (Farokhniaee and Lowery, 2021). TCM has limitations that could be addressed in the future studies such as inclusion of naturalistic stimuli to thalamus and simulate a suprathreshold dynamics for non-rest situations. The synaptic model of the network can be substituted easily with other known models of synaptic activity. Inclusion of long-term synaptic mechanisms such as spike-timing dependent plasticity (STDP) is another aspect that can be addressed in future studies. We suggest researchers and readers to use the improvised model presented here as the platform for future developments of DBS network model. The computational code is available online (see Data Availability Statement).

Data availability statement

The improved model presented in this study can be found in the TCM_DBS repository: https://github.com/aafarokh/TCM_DBS.

Author contributions

AF: Writing – original draft, Resources, Visualization, Formal analysis, Conceptualization, Validation, Supervision, Methodology, Writing – review & editing. SA: Resources, Writing – original draft, Validation, Formal analysis, Investigation, Methodology.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

SA was employed by Didab PVT Ltd Co.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Anderson, R. W., Farokhniaee, A, Gunalan, K., Howell, B., and McIntyre, C. C. (2018). Action potential initiation, propagation, and cortical invasion in the hyperdirect pathway during subthalamic deep brain stimulation. Brain Stimul. 11, 1140–1150. doi: 10.1016/j.brs.2018.05.008

PubMed Abstract | Crossref Full Text | Google Scholar

Arbuthnott, G. W., and Garcia-Munoz, M. (2017). Are the symptoms of Parkinsonism cortical in origin? Comput. Struct. Biotechnol. J. 15, 21–25. doi: 10.1016/j.csbj.2016.10.006

PubMed Abstract | Crossref Full Text | Google Scholar

Caligiore, D., Helmich, R. C., Hallett, M., Moustafa, A. A., Timmermann, L., Toni, I., et al. (2016). Parkinson's disease as a system-level disorder. NPJ Parkinson's Dis. 2:16025. doi: 10.1038/npjparkd.2016.25

PubMed Abstract | Crossref Full Text | Google Scholar

Farokhniaee, A., and Lowery, M. M. (2020). “A thalamo-cortex microcircuit model of beta oscillations in the parkinsonian motor cortex,” in 2019 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (Berlin: IEEE), 2145–2148.

PubMed Abstract | Google Scholar

Farokhniaee, A., and McIntyre, C. C. (2019). Theoretical principles of deep brain stimulation induced synaptic suppression. Brain Stimul. 12, 1402–1409. doi: 10.1016/j.brs.2019.07.005

PubMed Abstract | Crossref Full Text | Google Scholar

Farokhniaee, A., and Lowery, M. M. (2021). Cortical network effects of subthalamic deep brain stimulation in a thalamo-cortical microcircuit model. J. Neural Eng. 18. doi: 10.1088/1741-2552/abee50

PubMed Abstract | Crossref Full Text | Google Scholar

Izhikevich, E. M. (2003). “Simple model of spiking neurons,” in IEEE Transactions on Neural Networks (New York: IEEE), 14.

Google Scholar

Li, Q., Ke, Y., Chan, D. C. W., Qian, Z.-M., Yung, K. K., Ko, H., et al. (2012). Therapeutic deep brain stimulation in parkinsonian rats directly influences motor cortex. Neuron. 76, 1030–1041. doi: 10.1016/j.neuron.2012.09.032

PubMed Abstract | Crossref Full Text | Google Scholar

Lindenbach, D., and Bishop, C. (2013). Critical involvement of the motor cortex in the pathophysiology and treatment of Parkinson's disease. Neurosci. Biobehav. Rev. 37, 2737–2750. doi: 10.1016/j.neubiorev.2013.09.008

PubMed Abstract | Crossref Full Text | Google Scholar

Reis, C., Sharott, A., Magill, P. J., van Wijk, B. C. M., Parr, T., Zeidman, P., et al. (2019). Thalamocortical dynamics underlying spontaneous transitions in beta power in parkinsonism. Neuroimage 193, 103–114. doi: 10.1016/j.neuroimage.2019.03.009

PubMed Abstract | Crossref Full Text | Google Scholar

Tsodyks, M. V., and Markram, H. (1997). The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. Proc. Natl. Acad. Sci. 94, 719–723. doi: 10.1073/pnas.94.2.719

PubMed Abstract | Crossref Full Text | Google Scholar

Wiecki, T. V., and Frank, M. J. (2010). Neurocomputational models of motor and cognitive deficits in Parkinson's disease. Prog. Brain Res. 183, 275–297. doi: 10.1016/S0079-6123(10)83014-6

PubMed Abstract | Crossref Full Text | Google Scholar