The effect of deep brain stimulation on cortico-subcortical networks in Parkinson’s disease patients with freezing of gait: Exhaustive exploration of a basic model

Abstract Current treatments of Parkinson’s disease (PD) have limited efficacy in alleviating freezing of gait (FoG). In this context, concomitant deep brain stimulation (DBS) of the subthalamic nucleus (STN) and the substantia nigra pars reticulata (SNr) has been suggested as a potential therapeutic approach. However, the mechanisms underlying this approach are unknown. While the current rationale relies on network-based hypotheses of intensified disinhibition of brainstem locomotor areas to facilitate the release of gait motor programs, it is still unclear how simultaneous high-frequency DBS in two interconnected basal ganglia nuclei affects large-scale cortico-subcortical network activity. Here, we use a basic model of neural excitation, the susceptible-excited-refractory (SER) model, to compare effects of different stimulation modes of the network underlying FoG based on the mouse brain connectivity atlas. We develop a network-based computational framework to compare subcortical DBS targets through exhaustive analysis of the brain attractor dynamics in the healthy, PD, and DBS states. We show that combined STN+SNr DBS outperforms STN DBS in terms of the normalization of spike propagation flow in the FoG network. The framework aims to move toward a mechanistic understanding of the network effects of DBS and may be applicable to further perturbation-based therapies of brain disorders.

Parkinson's disease is a progressive neurodegenerative disorder with cardinal motor 2 symptoms such as axial and limb bradykinesia, rest tremor, and rigidity [1]. Many PD 3 symptoms are successfully treated pharmacologically [2] or by deep brain 4 stimulation [3,4]. However, axial symptoms, such as parkinsonian gait disorder, postural 5 instability and freezing of gait show limited response to treatment [5,6]. FoG is 6 associated with an increased risk of falls among the patients and represents a significant 7 source of morbidity. Thus, there is a need to reconcile FoG network pathophysiology 8 and efficient treatment approaches. 9 DBS is an established treatment strategy for PD patients with motor fluctuations, 10 medically refractory tremor or severe dopaminergic drug reactions [7]. Despite wide use 11 of DBS in the clinical routine, the mechanism of its action remains poorly understood. 12 DBS efficacy varies largely between patients, and particularly gait disorders and FoG 13 appear to respond better to specific stimulation patterns [5]. It was recently shown that 14 simultaneous deep brain stimulation of STN and SNr (STN+SNr DBS) outperforms 15 standard STN DBS in improving FoG symptoms [8][9][10]. However, there is a possibility 16 of worsening akinesia [8], hypomania [11], mania [12], and depression [13] with SNr 17 stimulation, partly due to projections of the SNr to the limbic system [5]. Therefore, 18 while it is an attractive alternative to induce co-stimulation of SNr with STN when 19 treating patients with DBS-resistant FoG, it is essential to assess and understand its 20 mechanistic effect on cortico-subcortical networks. 21 The effects of DBS are thought to rely not only on the specific properties of 22 individual neurons, but on the properties of large-scale brain networks [14]. The 23 STN+SNr DBS efficacy in alleviating FoG is based on the notion that gait disturbances 24 during advanced PD are associated with defective motor processing in the 25 mesencephalic locomotor region (MLR) [15]. MLR is densely interconnected with SNr, a 26 major output nucleus in the basal ganglia. Thus, a faulty output of the basal ganglia 27 system in PD might result in over-inhibition of MLR due to GABAergic projections of 28 the SNr and subsequent attenuation of locomotor activity [16]. In turn, functional 29 inhibition of STN and SNr due to DBS [17] is likely to attenuate over-inhibitory SNr 30 output and lead to restored gait and posture. To study the DBS network effects, we 31 represent a cortico-subcortical network of regions involved in FoG episode generation as 32 a directed signed graph. The network includes the basal ganglia and brainstem regions 33 as nodes, and the excitatory or inhibitory synaptic connections as signed edges. The 34 graph representation of the FoG network facilitates examining its topological properties, 35 and allows in-silico simulations of neural activity.
can observe all the possible emerging dynamical patterns in the network after 49 initializing the graph nodes with all possible initial state conditions. After some time, 50 the dynamics of the network will converge to some repeating patterns, the network 51 attractors. The change in the graph topology, due to projections silenced during DBS, 52 results in a change in the attractor space of the system. This space, in turn, can be 53 converted into a coactivation matrix, as a measure of functional connectivity. Thus, a 54 change in network topology results, via its interpretation through the excitable model, 55 in a change of the functional connectivity of the network.   Network-based computational framework to explain the effects of targeting different DBS sites. The SER model is applied to any possible initial conditions in which network nodes can be. This process leads to the emergence of different patterns of excitable dynamics, the network attractors. There are two types of attractors the system can be in: a limit cycle or a fixed point. During a limit cycle, the system repetitively goes through a set of S -susceptible, E -excited, and R -refractory states (attractor 1). By contrast, in the fixed point, all of the network nodes remain in the S state (attractor 2). The attractor space constitutes the coactivation and spike propagation matrices of the system. A change in the network topology due to a DBS target site choice affects the coactivation and spike propagation matrices. Thus, one can infer a suitable stimulation target from the changes in these matrices.
In this study, we explore the changes in dynamical landscapes emerging from the 57 FoG network via the SER dynamics during the STN and STN+SNr DBS modes. We 58 strive to detect the changes in activity propagation in the MLR leading to alleviation of 59 FoG. We hypothesize a normalization of activity propagation in the MLR due to

66
To study the changes in emerging dynamics induced by the system's state, we employ 67 the network shown in Fig. 2 B. The network includes the basal ganglia and brainstem 68 regions thought to be responsible for FoG episode generation. The FoG network can be 69 in the healthy, PD, STN DBS, and STN+SNr DBS configurations. The PD 70 configuration is obtained from the healthy configuration by setting the weights of the 71 edges originating from the SNc to 0. The same is done for the STN DBS and STN+SNr 72 DBS configurations based on the PD state, by setting the weights of the edges 73 originating from STN, or both STN and SNr, to 0. In the following sections, we first 74 compare the network's topological properties in healthy and PD configurations, then 75 assess the emerging dynamical patterns and corresponding coactivation matrices. 76 Finally, we analyze the activity propagation in different FoG network configurations.

77
Further details can be found in the Materials and methods section.  One aspect of the topological differences between the healthy and the PD 86 configurations is a change in the betweenness centrality. In our case, major betweenness 87 centrality changes were observed for SNc, SNr, cortex, and PPN nodes of the graph 88 ( Fig. 3 A)  Another aspect of the changes in the network topology, which subsequently defines 96 the dynamics of the system, are the changes leading to the synchronization of the nodal 97 activities. This is of interest, as PD is known to promote abnormal synchronous activity 98 of the basal ganglia nuclei [21,22]. The graph theoretical measures thought to be 99 responsible for the nodal synchronization for the topology of the FoG network in the 100 healthy and the PD state are assessed in Fig. 3.

101
For the FoG network, the mean clustering coefficient, mean betweenness centrality, 102 February 17, 2023 5/20 and average shortest path length decrease between the healthy and the PD 103 configurations, while maximal betweenness centrality and standard deviation of the 104 betweenness centrality increase (Fig. 3). These changes, despite their intricate nature 105 and interactions, usually promote synchronization.

106
The changes in the modularity were also assessed for the healthy and PD FoG 107 network configurations. We observe that the PD state leads to a larger number of 108 communities in the graph (Fig. 4). This increase in modularity could be viewed as a 109 pathological fragmentation of the FoG network in the PD configuration, which might be 110 interpreted as a disconnection syndrome similar to [23].   Table 1. The numbers of fixed points and period-3 limit cycles 118 largely vary depending on the network configuration. It should be noted that for the 119 configurations with the smaller number of period-3 limit cycles (healthy and STN DBS 120 configurations), the number of unique period-3 limit cycles is also smaller.

121
An example of an emerging dynamical landscape for the healthy network 122 configuration is shown in Fig. 5. It is noticeable that the largest proportion of all initial 123 conditions converges to a fixed point ( Fig. 5 right inset, Table 1). In a fixed point, all 124 the network nodes are constantly in the susceptible state. All the other initial 125 conditions converge to 31 different period-3 limit cycles. The limit cycle with the largest 126 basin of attraction is shown in the left inset of Fig. 5. Every node in the graph The numbers of fixed points and period-3 limit cycles are also shown in percentages with respect to the total number of initial conditions (3 12 = 531.441). The last row shows the fraction of the limit cycle with the largest basin of attraction in a limit cycle space in percentages.
repeatedly goes through a loop of S, E, and R states. The coactivation pattern of nodes 128 is different for different period-3 limit cycles. Additionally, Table 1 shows that the 129 period-3 limit cycle with the largest basin of attraction constitutes more than 30% of For the different combinations of the initial conditions of the SER dynamics, the healthy network topology converges to fixed points (right inset) and period-3 limit cycles (the largest one is shown in the left inset).
From an additional analysis, we note that for the 42% of the limit cycles in the 132 healthy network configuration, the striatum is always in a susceptible state S. By 133 contrast, this is only the case for 4% of the PD limit cycles. During the STN DBS and 134 the STN+SNr DBS, the proportion rises back to 17%.

135
To compare the emerged limit cycles for different topological configurations of the 136 network, we utilized Venn diagrams (Fig. 6). 137 We created PD and DBS topological configurations of the FoG network by changing 138 the weights of all edges coming from a certain node to 0 (see section Materials and 139 methods for further details). In a deterministic SER model, this approach is equivalent 140 to an external control exerted over the healthy, or, in the DBS case, PD FoG network 141 configuration [24]. When an external control is exerted over a network, the basin of original attractors of the network can change, or new attractors can emerge [25]. Thus, 143 a successful intervention, in our case DBS, can be viewed as the steering of the PD 144 system towards an original healthy attractor space.

145
The healthy and the PD attractor spaces are compared in Fig. 6 A. It can be seen 146 that during the exerted external control process of moving from the healthy to the PD 147 configuration, new attractors have emerged (red attractors in Fig. 6 A). In addition, the 148 basins of attraction of the original attractors have changed (right panel in Fig. 6 A).

149
Moreover, it is noticeable that healthy and PD configurations share several common 150 limit cycles. An additional analysis showed that the limit cycle with the largest basins 151 of attraction in the PD configuration does not appear in the healthy configuration.

152
From the right panel in Fig. 6 A, one can see that the fraction of the basins of 153 attraction of the unique healthy limit cycles is much smaller than the fraction of the 154 ones that appear in both the healthy and the PD configurations. Thus, for DBS, it is 155 less likely to steer the system back towards its unique healthy attractors than to the 156 attractors shared in common between the configurations, or to some new attractors.

157
In Fig. 6 B, the healthy, the PD, and the STN DBS limit cycle spaces are compared. 158 One can see that there are only two new attractors emerging from the PD attractor 159 space (in yellow). At the same time, the dynamical landscape of the STN DBS network 160 includes one of the healthy attractors, which are not the same as for the PD state (in 161 dark blue). This attractor has the striatum node always in a state S, which again 162 points out the quiescence of the striatal neurons in the healthy conditions [26]. From an 163 additional analysis, we found that in the STN DBS configuration the system has moved 164 away from the largest limit cycle in the PD configuration. additional analysis, we found that in the STN+SNr DBS configuration, the network has 173 also moved away from the largest limit cycle in the PD configuration.  The coactivation matrices in the healthy and the PD states differ (Fig. 7 A, B). The 179 pattern of activity in the healthy matrix appears to be smooth, whereas the pattern is 180 more checkerboard in the PD coactivation matrix. Fig. 7  There is a hypothesis that the normalization of excitatory signals propagation in 187 brainstem regions might be the reason for successful STN+SNr DBS during FoG [27].   In this study, we present a new network-based computational framework to uncover 201 dynamical landscapes and activity propagation in the brain networks relevant for 202 Distance to healthy state previous clinical studies which found that STN+SNr DBS outperforms STN DBS in the 209 context of FoG symptoms [8][9][10]. According to the hypothesis in [16], SNr stimulation 210 leads to temporal regularization of gait and normalization of the output to the MLR 211 and brainstem projections signaling, which is corroborated by the present results.

212
To achieve this result, we first checked topological properties of the FoG network in 213 the healthy and PD configurations. We studied graph theoretical measures and assessed 214 the corresponding FoG network synchrony in the PD configuration (Fig. 3). The PD 215 FoG network topology likely promotes synchronized activity, which additionally 216 validates the connections in the FoG network (Fig. 2 B). Synchronicity is arguably not 217 the most straightforward measure for assessing whether the graph network is closer to 218 the PD state or not. There are also small-worldness [28], global efficiency, calculation of 219 minimum spanning trees, and other approaches [29]. However, when taken alone, these 220 measures typically lead to inconsistent results across the literature [29,30] or are 221 inefficient in the case of signed [31] or directed graphs. As the directionality of the 222 graph is an essential biological parameter in our modeling framework, we needed to take 223 this property into account. In turn, the network being a signed graph is was needed for 224 the utilized SER model and provided the framework for its computational power. That 225 is why we assessed the network synchronicity, which is also an essential biological 226 property of PD [21,22]. 227 Thereafter, we studied dynamical landscapes of the FoG network in healthy, PD, 228 STN DBS, and STN+SNr DBS configurations using the discrete excitable SER model. 229 The modeling approach might seem simplistic, although it exhibits the most important 230 properties of the three-state excitatory cycle of neural elements. At the same time, the 231 amount of biological realism required for an optimal spiking model to capture the PD 232 dynamics is an open question. Our network-based SER framework is not intended to 233 render realistic neuronal firing patterns. Nevertheless, it captures essential aspects of 234 the PD and DBS dynamics. Analysis of attractors with the largest basins of attraction 235 (Table 1) reveals that for the PD FoG network configuration, the largest basin of 236 attraction constitutes 35% of the limit cycle space. This is a large increase compared to 237 the healthy configuration (15% of limit cycle space). This increase might be interpreted 238 as a sign of abnormal synchrony within the basal ganglia nuclei during PD [32]. During 239 the STN DBS, the amount of limit cycle space taken up by the attractor with the  The network-based SER framework allows exhaustive assessment and comparison of 252 the emerging dynamical landscapes (Fig. 6). This approach helps to answer one of the 253 long-standing questions in DBS research, namely the working mechanism behind it. We 254 can see that the system moves away from the most prominent attractors in the PD 255 configuration for both DBS conditions. Although moving towards the healthy 256 configuration attractors is less likely (Fig. 6 A, right panel), both DBS attractor spaces 257 incorporate one unique healthy attractor. Otherwise, the system for the DBS exhibits a 258 trend to move away to a new attractor space (yellow parts of Fig. 6 B, C). This 259 observation is in accordance with the conclusions by Wichmann and DeLong [32]. They 260 state that the system for the DBS reaches a new equilibrium, rather than induces 261 restoration of normal basal ganglia function. This new equilibrium releases the 262 downstream MLR areas from pathological basal ganglia output. 263 We continued to study the new equilibrium DBS states by exploring the coactivation 264 matrices and spike propagation flow in the FoG network. The coactivation matrices 265 revealed that combined STN+SNr DBS is closer to the healthy state in terms of 266 functional connectivity (Fig. 7). Furthermore, the checkerboard pattern of activations in 267 the PD coactivation matrix reaffirmed the fragmentation of the FoG network into more 268 communities, similar to the results in Fig. 3. The spike propagation flow normalization 269 in the MLR (Fig. 8) exclusive to the STN+SNr DBS reaffirms the hypothesis on its 270 working mechanism by [16]. Additionally, it captures the normalization of the 271 cortico-thalamic connectivity towards the healthy controls for DBS [35].

272
The healthy FoG network topological configuration governs the resulting dynamical 273 output, as network topology is the most influential parameter in the SER model. In the 274 present study we used the topology described in Fig. 2 B. The initial choice of this 275 signed graph heavily affects the observed dynamical patterns. We motivated this 276 network choice by comparing the healthy FoG network configuration with a mesoscale 277 connectome from [36] ( S1 Fig). Specifically, we compared the weights of the edges from 278 the Allen brain atlas data in [36] with the corresponding weights of the edges in the There are no edges with large weights in [36], which would correspond to a 0 weight in 284 the FoG network. Thus, the healthy FoG network configuration in Fig. 2 B depicts the 285 most significant connections in the Allen brain atlas connectome [36], which supports 286 the reliability of the chosen topology.

287
Overall, we propose a practical computational framework that simultaneously 288 captures the essential aspects of PD and DBS, explores the new STN+SNr DBS 289 approach and confirms its efficiency in the context of FoG.

291
Selecting the best therapeutic strategy for a specific symptom in PD can be 292 burdensome. In this study, we propose a straightforward network-based computational 293 framework to determine the most suitable DBS target for treating FoG. This approach 294 is based on a discrete SER model and requires only preliminary knowledge of the 295 cortico-subcortical network topology. We also showed that STN+SNr DBS is more 296 beneficial for patients with FoG due to the normalization of the spike propagation flow 297 in the MLR, which confirms our initial hypothesis. In the future, we aim to explore 298 additional DBS targets and network regions, specifically the limbic system, to check the 299 influence of STN+SNr DBS on mood and anxiety.

301
The FoG network 302 The freezing of gait network (Fig. 2 B) was used to compare possible network effects 303 induced by STN DBS and STN+SNr DBS. These effects were then used to identify the 304 possible working mechanisms behind DBS and to infer a suitable stimulation mode for 305 FoG. The network consists of the basal ganglia, the motor cortex, the thalamus, the 306 pedunculopontine nucleus, the pontine reticular formation, the cuneiform nucleus, and 307 the locus coeruleus. The choice of the network regions and connections was based on an 308 analysis of the anatomical literature [37][38][39]. A schematic for locomotion control [37] 309 was a primary resource used for the construction of the network. The brainstem 310 circuitry was also partially taken from the Mouse Brain Connectivity Atlas provided by 311 Allen Institute [36,40] in order to expand the network. being the working mechanism behind the DBS [17]. The hypothesis is based on the 326 similarity between the effects observed after surgical lesions of the brain regions and 327 their high-frequency DBS. According to this idea, DBS is thought to function via 328 inhibition of the neurons in the vicinity of the stimulating electrode. In our model, this 329 phenomenon was taken to be equivalent to the disconnection of the DBS target region. 330 Thus, in the case of STN DBS, the weights of all the edges originating from the STN 331 node were set to 0. During the combined STN+SNr DBS, the weights of the edges 332 originating from both the STN node and the SNr node were set to 0.

333
Graph theoretical measures 334 Changes in network topology can affect the resulting network dynamics [42][43][44]. The 335 dynamics can become more or less synchronized among various network nodes upon the 336 changes in graph theoretical measures, which assess different aspects of network 337 topology. Here we focus on the measures often thought to affect node 338 synchronization [45][46][47]: clustering coefficient, mean betweenness centrality, maximal 339 betweenness centrality, and average shortest path length.

340
Clustering coefficient corresponds to the density of connections between a node's 341 neighbors [42]. When the clustering coefficient is high, the node is densely connected 342 with its neighbors and forms a local cluster. The mean clustering coefficient represents 343 the average of clustering coefficients across all the graph nodes. A decrease in the mean 344 clustering coefficient usually promotes synchronization [47].

345
Betweenness centrality reflects the node's influence on the activity flow in the 346 graph [48]. It is based on the node's ability to participate in the network's shortest path. 347 This participation results in the node controlling a network's activity [42]. The maximal 348 betweenness centrality of the graph corresponds to the largest betweenness centrality 349 among the nodes, whereas the mean betweenness centrality corresponds to the mean of 350 the betweenness centrality across nodes [45]. For the betweenness centrality measures, 351 the relationship with synchronization is complex. Synchronization tends to increase 352 with the maximum betweenness centrality [47] while in contrast it decreases with the 353 mean betweenness centrality [45]. However, there are more complex cases where one 354 should observe the heterogeneity of the betweenness centrality of the network [46,47]. 355 We measure heterogeneity of the betweenness centrality by assessing its standard 356 deviation across nodes. More homogeneous distributions of centralities tend to promote 357 more synchronization [46]. On the other hand, the authors in [47] report the opposite, 358 indicating that this measure cannot solely assess synchronization in complex networks. 359 Average shortest path length relates to the efficiency of the information processing in 360 a network and corresponds to the average shortest path length between every pair of 361 nodes in the graph [47]. Usually, a decrease in the average shortest path length 362 promotes synchronization [45,47].

363
Modularity assesses the extent to which a network can be segregated into connected 364 neighborhood [42]. Here it was evaluated using the Leiden community detection method 365 for directed and weighted graphs [49].

366
The SER model 367 We studied all possible dynamical patterns which emerge in the healthy, PD, and DBS 368 FoG networks. Thus, for computational effort, a minimal excitable dynamical system, a 369 discrete excitable SER model [18,19] was used. In the SER model, the S stands for 370 susceptible, E excited, and R refractory states. In the deterministic version of the 371 model, the node goes through the aforementioned states. If the node is in a susceptible 372 state S and if the sum of the weights of the edges originating from excited neighbouring 373 nodes and converging on the target note is larger than 0, then in the following discrete 374 time step, this node will become excited (E). From the E state, the node will always go 375 to a refractory state R. In turn, from the R state, the node will always go to the state 376 S in the following discrete time step. 377 We applied the SER model over all possible initial conditions of the FoG network for 378 the healthy, PD, and DBS network configurations to explore and compare emerging 379 dynamical patterns. At the beginning, every node of the FoG network could be in one 380 of the states: S, E, or R. After considering all possible combinations of these initial 381 conditions, every possible dynamical pattern which can emerge in the FoG network for a 382 certain topological configuration was explored. After that, the dynamical landscapes 383 were compared between the healthy, PD, STN DBS, and STN+SNr DBS FoG network 384 configurations.

385
As the SER model is a three-state model, there are 3 n possible patterns of activity 386 in a network with n nodes. These patterns of activity represent the states of the 387 network. Thus, the network dynamics can be represented as a time series of the states. 388 As the amount of the states is finite for a system with a finite number of nodes, the 389 dynamics of the network eventually converges to one of two types of behavior, the 390 attractors. The first type of attractor is a fixed point. The system reaches a fixed point 391 after some transient time when all nodes in the network remain in the state S. The 392 other attractor type is a limit cycle. In this case, all the nodes in the network will go 393 through a series of the S, E, and R states with a period of 3k, where k ∈ N . Attractors 394 could be compared via their basin of attraction, which is a set of initial states that 395 converge onto a particular attractor [24]. Thus, to compare the dynamical landscapes in 396 the various FoG network states, we compared the emerging attractors and their basins 397 of attraction for different FoG network configurations and the SER dynamics after the 398 time T = 100 with the transient time t = 40.

399
Additionally, we studied coactivation and activity propagation patterns of the 400 network in different configurations. Coactivation matrices were obtained from the time 401 series of the network states by summing up simultaneously occurring spiking events (E 402 states).

403
To calculate the activity propagation flow in the FoG network under different 404 topological configurations, we calculate 1-norm distances between coactivation matrices 405 corresponding to the FoG network projections. We assessed the activity propagation 406 flow differently depending on whether the projection in the FoG network was excitatory 407 or inhibitory.

408
To study excitatory projections, we used the shifted version of the coactivation 409 matrices ( Fig. 9 A). It is obtained by shifting the activity of the projection target region 410 one time step further from the activity source region. This way, if the excitatory state 411 E co-occurs in the projection target and the source region, it means that the excitation 412 of the target follows the excitation of the source. In our model, this was the case for the 413 excitatory projections.

414
To assess activity propagation flow along the inhibitory projections of the FoG 415 network, we used the usual coactivation matrices (Fig. 9 B). If the excitatory state E 416 occurred simultaneously in the projection target and the projection source region, it 417 meant that the excitation of the source was followed by the target in the refractory 418 state R on the next time step. This would be true only after the transient time when 419 the system has reached its stable dynamical pattern. In our model, the excitation of the 420 source followed by the target region entering the refractory state wwas equivalent to 421 inhibition. Thus, we utilized the coactivation matrices to assess the spike propagation 422 along the inhibitory projections.

423
After obtaining the shifted and synchronous coactivation matrices, we calculated excitatory projections which exist in the FoG network (Fig. 8).