Complex network measures reveal optimal targets for deep brain stimulation and identify clusters of collective brain dynamics

An important question in computational neuroscience is how to improve the efficacy of deep brain stimulation by extracting information from the underlying connectivity structure. Recent studies also highlight the relation of structural and functional connectivity in disorders such as Parkinson’s disease. Exploiting the structural properties of the network, we identify nodes of strong influence, which are potential targets for Deep Brain Stimulation (DBS). Simulating the volume of the tissue activated, we confirm that the proposed targets are reported as optimal targets (sweet spots) to be beneficial for the improvement of motor symptoms. Furthermore, based on a modularity algorithm, network communities are detected as set of nodes with high-interconnectivity. This allows to localise the neural activity, directly from the underlying structural topology. For this purpose, we build a large scale computational model that consists of the following elements of the basal ganglia network: subthalamic nucleus (STN), globus pallidus (external and internal parts) (GPe-GPi), extended with the striatum, thalamus and motor cortex (MC) areas, integrating connectivity from multimodal imaging data. We analyse the network dynamics under Healthy, Parkinsonian and DBS conditions with the aim to improve DBS treatment. The dynamics of the communities define a new functional partition (or segregation) of the brain, characterising Healthy, Parkinsonian and DBS treatment conditions.

1 Introduction

Movement disorders like Parkinson’s disease and dystonia are characterised by abnormal functioning of the whole basal ganglia (BG) - thalamocortical network. In Parkinson’s disease, one of the main characteristics of the altered BG network behaviour is a synchronised abnormal β − activity (12–35 Hz) Kühn et al. [1]; Neumann et al. [2]. This enhanced BG rhythm also affects thalamic activity by sending strong inhibitory signals via GPi (internal segment of the globus pallidus). The hyperpolarisation of thalamic neurons due to increased inhibitory BG output increases the burst discharges Kim et al. [3]; Galvan et al. [4], which in turn triggers motor dysfunction through the thalamocortical pathway Kim et al. [3].

Deep brain stimulation (DBS) of the BG was shown to be an efficient treatment for movement disorders Deuschl et al. [5]; Vidailhet et al. [6], but its therapeutic mechanism is still not fully understood. The application of DBS leads to firing pattern alterations, and in particular, to disruption of hypersynchronised β-band oscillations Kühn et al. [1]; Kim et al. [3]; Crompe et al. [7]. Indeed, recordings in animal models Xu et al. [8]; Kim et al. [3]; Crompe et al. [7] and observations from computational network models Rubin and Terman [9]; Popovych and Tass [10]; So et al. [11]; Galvan and Wichmann [12] suggest that DBS in the subthalamic nucleus (STN) results in more periodic and regular firing at higher frequencies in the BG-thalamic network.

Two important questions arise concerning the structural connectivity and DBS effectiveness. The first question is how to determine those electrode positions which are the most effective for the activation of neural pathways to improve DBS outcome. The second question which naturally arises, is which differences in neural activation patterns will emerge within the brain’s structural network when simulating different conditions (i.e., Healthy, Parkinsonian and DBS). Due to the strongly heterogeneous nature of the connection topology and the stochastic and nonlinear large scale interactions of the underlying units/neurons, the emergent macroscopic behavior usually is far from trivial to predict Spiliotis and Siettos [13]; Siettos and Starke [14]; Deco et al. [15,16]; Bassett and Bullmore [17]; Bullmore and Sporns [18]; Iliopoulos and Papasotiriou [19]. Self-organisation, sustained oscillations, travelling waves, multiplicity of stationary states and spatio-temporal chaos are paradigms of the rich nonlinear behaviour at the coarse-grained systems level Spiliotis and Siettos [13]; Siettos and Starke [14]; Deco et al. [15,16]; De Santos-Sierra et al. [20]; Crowell et al. [21]; Spiliotis et al. [22], indicating thus that a precise structure–function relation remains a major open problem Deco et al. [16].

As the primary aim of this paper, we identify positions of nodes of high functional impact using network analysis. We test then the hypothesis that such high-connectivity nodes are pivotal in shaping network activity and are highly effective stimulation targets for restoring the normal network function. Notably, the electric field-based approximation of the volume of tissue activated (VTA) Butson et al. [23]; Butenko et al. [24], computed at these nodes for a monopolar DBS, overlapped with target areas previously shown to be effective against akinetic symptoms of Parkinson’s disease Dembek et al. [25].

The critical high-connectivity nodes were identified according to three different, but interrelated measures:

1) Using the measure of “clustering coefficient,” nodes are identified which form triplet interactions. This triangular interconnection allows for circular information flow and information feedback. This triplet organisation constitutes the complex level of connectivity, and is speculated to play a role in e.g. effective information distribution but also complex oscillatory network rhythm formation.

2) Using the measure of “betweenness centrality,” high-connectivity nodes are identified (also known as hubs and defined by their so-called nodal efficiency van Hartevelt et al. [26]). This nodal efficiency is related to the degree of influence of nodes in the network and can be interpreted as the amount of flow that passes through these nodes. Such hubs act as central crossroads, enhancing the ability of parallel information transfer and the functional integration in brain networksvan Hartevelt et al. [26].

3) Using the measure of “eigencentrality” (or eigenvector centrality), nodes are identified which specifically connect with other nodes of high centrality. In this way, targeting such nodes will likely influence a large population of other nodes. In many cases, centrality measures correlate strongly Li et al. [27]; nodes with extreme values of “betweenness centrality” also show high “eigencentrality” values.

The second aim of this study was to explore the relationship between anatomical structure and neural activity (i.e., functional connectivity) using modified Hodgkin-Huxley models Terman et al. [28]; Rubin and Terman [9]; So et al. [11]; Spiliotis et al. [22]. In the current study, we address this question by combining the community structures (i.e., sets of high-connectivity nodes as identified using modularity measures, Newman [29]) and a large scale biophysical model which produces virtual neural activity. We study three different conditions: Healthy, Parkinsonian state and Parkinsonian conditions with DBS. In the latter case, we extend the previous computational model Spiliotis et al. [22] in order to simulate specific spatial positions of DBS electrodes Mandali et al. [30].

Predicting DBS outcome using neural networks is not novel Rubin and Terman [9]; Popovych and Tass [10]; Spiliotis et al. [22]; Fleming et al. [31], in the cited studies, however, smaller networks were studied without taking into account connectivity structure and VTA. In the current study, we follow a different approach: initially, we integrate the high dimensional nonlinear system, which produces spatiotemporal patterns consistent with either Healthy (normal) or Parkinsonian states. Then we average the activity over the different community structures that have been previously identified using modularity network measures.

The second main result of the study shows that in all areas (the detected communities), including neocortical ones, Parkinsonian conditions alter power spectrogams, but mainly subcortical structures, with e.g., slowing of activity in the thalamus and faster activity in pathways connecting pallido-thalamic and subthalamic-pallidal nodes. Under DBS, in turn, the simulation reveals that this stimulation at high-connectivity nodes is able to restore thalamic activity, and partly also cortical one, while on the one hand, hyperdirect-pathway associated nodes remain largely unaffected, and on the other, structures in close vicinity of the electrode mainly follow the stimulation.

2 Network connectivity from data sources

To describe the structural connectivity of the network, data from different studies on human brain anatomy were utilised.

2.1 Data sources

For the current model, published data on anatomical and fiber tract positions were used from different sources: The basal ganglia nuclei and their substructures were taken from the DISTAL Atlas Ewert et al. [32]; Chakravarty et al. [33]. The Melbourne Subcortex Atlas Tian et al. [34] was used to define substructures of the thalamus, while the relevant cortical regions were selected using the Brainnetome Atlas parcellation Fan et al. [35]. Fiber tracts classified to pathways in the vicinity of the STN were taken from Petersen et al. [36], and projections of the ventral anterior nucleus to motor cortical regions, required to complete the BG-thalamo-cortical network, were extracted from the structural group connectome of 90 PPMI Parkinsons’s disease-patients GQI Marek et al. [37], post-processed in Ewert et al. [32]. All data were represented in MNI (Montreal Neurological Institute) space.

2.2 Structural connectivity of the basal ganglia-thalamo-cortical network

To investigate DBS network effects in Parkinsons’s disease, the classic circuit model of the BG-thalamo-cortical network was employed Milardi et al. [38] (Figure 1A). It consists of three inputs from motor cortical regions: the direct pathway that involves the striatum and continues as a GABAergic projection to the GPi and SNr (substantia nigra pars reticulata); the indirect pathway that also involves the striatum, but has GABAergic projections to the GPe (external segment of the globus pallidus), which in turn inhibits the STN, GPi and SNr. In addition to these pathways, there is the hyperdirect pathway through which the STN receives a direct excitatory input from the cortical areas. The glutamatergic efferents of the STN innervate the GPe, GPi and SNr. GABAergic projections of the two latter nuclei to the ventral anterior (VA) and ventral lateral (VL) regions of the thalamus represent the output of the BG circuit to the thalamo-cortical network Bosch-Bouju et al. [39].

FIGURE 1

FIGURE 1. Basal ganglia-thalamo-cortical circuits. (A): Circuit model comprising all main connections. MC/PMC—motor and premotor cortical regions, respectively. For explanations, please see the main text. (B): Simulated reduced circuit model with highlighted (bold arrows) connections possibly affected by STN-DBS. The synaptic connection between Striatum and Globus Pallidus pars externa/interna (GPe and GPi) were modelled by different constant currents. (C): Structural connectivity of the simulated network is based on the pathway atlas of human motor network constructed from multimodal data including diffusion, histological and structural MRI data, fused to a virtual 3D rendering Petersen et al. [36]. (C1): Projections from GPi to thalamus (VA nucleus) are shown in red, connections between subthalamic nucleus (STN) and GPe are shown in green, projections from STN to GPi are shown in violet. (C2): Connections between motor cortex (MC/PMC) and the thalamus projections were obtained by classifying fiber tracts from Marek et al. [37] and are shown in orange. Projections from motor cortex to STN (hyperdirect pathway) are shown in blue. Nuclei are shown in the following colours: GPe in light grey, GPi in dark grey, STN in dark orange, thalamus (VA nucleus) in yellow and motor cortex (M1) in light grey.

The interstructural connectivity of the network was simulated using data based on the pathway atlas of human motor network obtained from multimodal imaging, including diffusion, histological and structural MRI data, fused to a virtual 3D rendering Petersen et al. [36] or classified based on their position relative to the involved structures (thalamo-cortical projections), see Figure 1. Note that the grouping of fibers leads to an emergence of network nodes. In the current study, the simulated network (Figure 1B) does not include the substantia nigra (SN). The reason for this omission is that the dopaminergic projections of the SNc (substantia nigra pars compacta) are not myelinated, and hence less excitable (by approx. two orders of magnitude) by extracellular fields Tarnaud et al. [40], yielding a very low chance that DBS actually would affect them directly. Furthermore, to ensure homogeneity of the BG pathways, only those present in Petersen et al. [36] were employed, which did not include the striatonigral and the nigrothalamic projections. For the same reason, only the VA neurons were simulated, as they predominantly receive the pallidal output, unlike the VL nucleus that is mostly innervated by the SNr afferents Lanciego et al. [41]. Note that although the projections from the VL to the motor cortex exist, we excluded them to avoid modelling of intrinsic dynamics between subregions of the thalamus.

3 Modelling structural basal ganglia-thalamo-cortical neuronal network using complex network theory and pathway classification

Structural brain connectivity refers to the set of anatomical links (or axonal tracts) which join different brain regions. The connectivity can be described and simplified employing elements of complex network theory Bullmore and Sporns [18]; Stam and Reijneveld [42], where the neural elements at the beginning and the end of a tract serve as nodes of the network, while the anatomical tract is described as edge of the network. The topological structure of the network plays an important role in the emergent neural activity and brain functionality, however it is not well understood how the structure or topology shapes the dynamics Deco et al. [15,43]. The knowledge of the network structural properties is important since it allows to build realistic computational models and to shed light on mechanisms underlying brain functionality or dis-functionality (i.e., movement disorders), and to predict the neural dynamics on multiple scales Honey et al. [44,45]. Using the structural connectivity of Section 2, we build a directed network. The internal connectivity structure within the areas STN, GPe, GPi, thalamus and motor cortex areas had to be defined using complex network theory Stam and Reijneveld [42]; Watts and Strogatz [46].

3.1 Construction of complex network using the structural connectivity

The resulting structure of Section 2.2 is used to obtain the connectivity network in the form G = (V, E), where V is the set of nodes and E represents the set of edges. The nodes of the structural network are defined as points in three-dimensional space and correspond to the starting and ending point of a fiber tract. The resolution is set at 1mm³, meaning that if two (or more) ending (or starting) points lying within the same cube of 1mm³, they are considered as one node. The connectivity information is included in the adjacency (or connectivity) matrix A: if there is a fiber tract starting at position x = (x₁, y₁, z₁) and ending at y = (x₂, y₂, z₂) then A(x, y) = 1, otherwise A(x, y) = 0. The resulting connectivity constitutes a graph G = (V, E), where the V is the set of nodes and E is the set of edges or tracts. At the given resolution of 1mm³, the network contains 134 STN nodes, 244 and 246 GPe and GPi nodes, respectively, 833 thalamic nodes and 2070 cortical nodes.

3.2 Intrastructural small-world network connectivity

The pathway classification analysis in Section 2.2 does not contain information about the internal connectivity of each region (i.e., how the nodes are connected within a region). We model connectivity within each area using small-world structures Watts and Strogatz [46]; Bassett and Bullmore [47]; Bullmore and Sporns [18]; Stam and Reijneveld [42]; Spiliotis and Siettos [13], thus increasing overall network connectivity beyond the interstructural projections. In such small-world complex networks Mark [48]; Watts and Strogatz [46], each node interacts with its k nearest neighbours; additionally, a few randomly chosen remote connections (with a small probability p) within the area are also formed Watts and Strogatz [46]. Small-world structures are commonly used in computational neuroscience Netoff et al. [49]; Berman et al. [50]; She et al. [51]; Bassett and Bullmore [17,47]; Fang et al. [52]; De Santos-Sierra et al. [20] as a result of two main characteristics which they show: highly clustered property together with short path length Bullmore and Sporns [18]; Watts and Strogatz [46]; Newman [29], enhancing in this way the signal or rhythm propagation within the network and the synchronizability in the network.

The GPe/GPi, thalamus and MC layers were modelled as separate small-world networks. Each node increases the initial number of connections (or the degree of the node) by k = 20 degrees on average. The local internal connections lay in a distance less than 5 mm (these are the local neighbours); however, the small-world topology Watts and Strogatz [46] allows remote connections (in a distance greater than 5 mm) with a small probability p = 0.05. The choices of k and p in this model are phenomenologically extracted. These values turn out to be successful i.e. the values k and p are chosen such that the network will give high values of clustering coefficient compared to a random network (where the clustering coefficient is very low and simultaneously a low value of the characteristic path length, see also Watts and Strogatz [46]. Similar values have been used in other studies and with the chosen values, the connectivity structure resembles real-neuronal connectivity as it is shown for example, in the work of De Santos-Sierra et al. [20] and Netoff et al. [49].

For the STN, we chose a modified small world approach which results from the experimental findings of Gouty-Colomer et al. [53]; Ammari et al. [54]. The STN area is characterised by sparse connectivity, where local and remote connections coexist. Specifically, only 20% of the STN neurons develop connections (collaterals) within the other STN neurons Gouty-Colomer et al. [53]. Almost 80% of these connections are local within a distance of 200–400 μ m radius, and the other 20% are contacts which occur farther away, i.e., $>500\mu$ m. In this sense, the 20% of neurons, which form STN connections, have both local and remote connections analogous to the small-world property Spiliotis et al. [22]. Similar to the previous connectivity, in our model, only the 20% of STN neurons show an average of 25 connections each, while few of these are randomly chosen remote connections Spiliotis et al. [22]; Gouty-Colomer et al. [53]).

3.3 Network properties: Degree distributions, path distances and centralities

Network measures such as a quantification of homogeneity are used to identify structural properties of the underlying neuronal network. These measures allow to categorise structural elements according to connectivity properties (i.e., profiles), segregating them into discrete entities. The main categories are clustering and distancing measures, centralities and communities detection Bullmore and Sporns [18]; Stam and Reijneveld [42]. Another subdivision of the network measures is the local and global description. The local description refers to the individual property of the i−th node, while the ensemble over the whole set of nodes in the network defines the global description (or distribution) for the network. The statistical distribution of a network property in this paper is characterised by its mean (the first-order statistical measure).

3.3.1 Degree distribution

The degree of a node i refers to the number of edges connected to it Bullmore and Sporns [18]. In directed networks, a node has both an in-degree and out-degree, which are the numbers of in-coming and out-going edges, respectively. A high degree of connectivity (increased numbers of links) of the i−th node defines the importance of a node in the network. The degree distribution P(k) defines the probability of a randomly selected node to have specific degree k. Averaging over all the nodes of the network, we obtain the mean degree, the first characteristic of the connectivity. The degree distribution after pathway classification analysis (Section 2.1), as long as internal connectivity is not considered, follows a power law of the form.

$P\left(K\right)=c{k}^{-\gamma }(1)$

where the exponent γ was calculated as γ ≈ 2.5.

The main characteristic of a power-law degree distribution is that only a few high-connectivity nodes acting as central nodes or hubs (nodes with a high number of connections) exist, while the majority of nodes show little connectivity. The high-connectivity nodes or hubs are responsible for an effective and fast spreading of information or signals in the network. Figure 2A depicts the degree distribution of P(k) in the network, while Figure 2B shows the degree distribution of the network including internal connections. In the latter case, the network thus combines both power law properties (describing connections among nuclei) and small world characteristics (describing local connectivity within nuclei). This generates an almost symmetric distribution of P(k).

FIGURE 2

FIGURE 2. Statistical properties of the connectivity network. (A) The degree distribution (depicted in logarithmic scale) follows a power law as Eq. 1 with critical exponent γ = 2.54. (B) The augmented network, using a small world network for the internal connectivity. The resulting distribution is a combination of power laws, meaning that the central nodes still exist and act as hubs. Furthermore, each neuron is connected on average with the 20 nearest neighbours in a radius of 5 mm. (C) The distribution of path lengths in the network. The mean value was calculated as ml = 4.43. (D) The distribution of clustering coefficients. There are few nodes with a high value of clustering coefficient.

3.3.2 Paths lengths, efficacy, and clustering coefficient

In graph theory, a path is a sequence of successive steps between two nodes, assuming that it is never intersecting a single node more than once. The minimum distance (i.e., the minimum steps in case of binary networks) between two nodes defines the shortest path length. Averaging over the set of all shortest paths, we obtain the mean path length of the network:

$\bar{m}=\frac{{\sum }_{i,j}{d}_{i\to j}}{N\left(N-1\right)},(2)$

The mean path length shows the ability of the network to spread information between any two nodes. A low mean shortest path length $\bar{m}$ signifies that any two randomly chosen nodes can interchange information via just very few intermediate nodes (in our case ≈ 4 nodes).

Another similar measure which is applied in the case where there is no connecting path between two nodes (i.e., d_i→j = ∞), is the global efficiency $\bar{G}$:

$\bar{G}=\frac{{\sum }_{i,j}\frac{1}{d_{i\to j}}}{N\left(N-1\right)}.(3)$

This measure avoids calculating with infinity, since if there is no pathway between two nodes i.e., d_i→j = ∞ ⇒ 1/d_i→j = 0. The $\bar{m}$ is comparable with inverse $\bar{G}$, and according to the Cauchy inequality for the arithmetic and harmonic mean, we obtain

$\bar{m}\ge \frac{1}{\bar{G}}.(4)$

For the augmented network, the mean path length was computed to be $\bar{m}=4.435$, while the inverse global efficacy resulted in $1/\bar{G}=3.91$. Figure 2C shows the distribution of distances between any two nodes (i.e., d_i→j).

Beyond the information flow between any two nodes, information flow among three nodes in a circular path, with the first node communicating with the second, and the second with the third, but the third communicating back to the first, another quality of information is made possible, i.e. feedback information. This would enable a circuitry to act in control loops, allowing for rhythm generation. To quantify this property, we introduce the clustering coefficient, which measures the local property of a node i to form triangle motifs. The clustering coefficient of a node i is defined as ratio:

$c\left(i\right)=\frac{{\sum }_{jk}{a}_{ij}{a}_{jk}{a}_{ki}}{k_i\left(k_i-1\right)}.(5)$

The higher the number of triangles (that exist) with respect to the i − th node, the higher the clustering coefficient. Figure 2D depicts the distribution of clustering coefficients. The mean clustering coefficient is computed as $\bar{c}=0.1$. The distribution shows the existence of few nodes with high values of c.

3.3.3 Betweenness centrality

Besides the ability to generate feedback-loops, information flow within a network is governed by the degree of interconnectivity between nodes. Centrality measures are used to identify such high-interconnectivity nodes in the network. The significance of a node is related to the degree of influence which it exerts in the network. For example, the influence can be interpreted as the amount of flow which passes from this node. Important nodes act as central crossroad or hubs in the network.

“Betweenness centrality” measures the amount of influence which a node has with respect to the total information flow in the network (serving as a bridge between subgraphs, i.e., sets of nodes of the network). The “betweenness centrality” (B_c) mathematically is defined as the fraction of all shortest paths in the network that pass through a given node, specifically the “betweenness centrality” of a node i is defined as:

$B_c\left(i\right)=\sum _{j\ne i\ne k}\left(g_{jk}\left(i\right)\right)/g_{jk},(6)$

where g_jk(i) is the number of shortest paths from j to k passing from node i, and g_jk is the number of shortest paths between nodes j and k. Bridging nodes that connect disparate parts of the network often have a high ‘betweenness centrality’. Higher values of B_c(i) indicate that the node acts as a central node influencing most of the other nodes in the network. The importance of these hubs is also highlighted pathophysiologically in the sense that therapeutic intervention in, e.g., Parkinson’s disease alters both the structural and functional connectivity profile in patients (a study which, however, obviously does not have any data on the Healthy state of the network as a basis of comparisons) van den Heuvel and Sporns [55].

Figure 3 depicts the distribution of B_c of the network. Indeed, the large majority of nodes shows very low centrality. However, there are few nodes with high B_c. In Figure 3B, black filled circles depict the spatial localisation of these central nodes in the network. As can be seen in this figure, these high-centrality nodes can be found in the STN, GPe and GPi, as well as the thalamus, but also in the motor cortex. We propose these nodes as very promising targets for DBS treatment. The coordinates of these hubs in MNI space, and the brain area they belong to, are given in Table 1.

FIGURE 3

FIGURE 3. Probability distribution (A,C) and position (B,D) of nodes with high centrality measures. (A,C): Probability distribution of centrality nodes. (A) The distribution of “betweenness centrality” follows a power law distribution. Thus, $>$90% of the nodes have a “betweenness centrality” $<$20,000, and only 0.5% a centrality value of 200,000. (C) The distribution of “eigencentrality” also follows a power law distribution. Thus, again $>$ 99% of the nodes have a value $<$ 0.2, and only 0.2% an “eigencentrality” value of 0.17. (B,D) Position of high centrality nodes (i.e., the nodes with the highest B_c or E_c values in each region, usually with a B_c > 200.000/E_c > 0.17) in MNI brain space coordinates shown as black dots. Brain regions are colour-coded as given in the figure. Most of these nodes are located in STN and GPi, thalamus (Tha) and GPe, but also some in the motor cortex (MC).

TABLE 1

TABLE 1. High centrality nodes, which might have high effectiveness in the DBS treatment.

3.3.4 Eigencentrality

Beyond betweenness centrality, identifying nodes with high-connectivity, nodes with high connectivity connected to other nodes with high connectivity represent a special type of network information distribution, since nodes with such high “eigencentrality” (E_c), are hypothesised to play an important role in fast and effective signal distribution within the network.

For each node in the network, a positive number x_i is assigned. The number x_i is set to be proportional to the sum of the weights of all nodes connected to i:

$x_i={\lambda }^{-1}\sum _j{A}_{ij}{x}_j\iff Ax=\lambda x,(7)$

where λ has to be identified. The last equation shows that the element x_i is the i−th element of the eigenvector of the adjacency matrix (corresponds to the eigenvalue λ). The eigencentrality (or eigenvector centrality) is defined when one chooses λ = λ₁ the highest eigenvalue of the adjacency matrix A. Then

$E_c\left(i\right)={\lambda }_1^{-1}\sum _j{A}_{ij}{x}_j,(8)$

which gives the eigenvector centrality the nice property that it can be large either because a vertex has many neighbours or because it has important neighbours (or both) Newman [56].

3.4 Detection of communities and modularity

Networks characteristically are made up of sets of nodes (subgraphs) which are densely connected among each other within the network, and which have sparse connections to other subgraphs Newman [29]. We hypothesise that such densely connected subgraph groups (or communities) play a significant role in information processing within the network. Assigning and allocating these densely connected communities to brain structures allows to construct a modular view of the network’s dynamics Newman [29].

In this paper, the modality index identifies such densely connected communities. The modality index Newman [29] assigns a community numbers_i to each node. For example, in the case of two communities, then s_i = ±1. Here, we seek the best network partition in order to optimise the modularity function Q:

$Q=\frac{1}{4m}{s}^{T}Bs(9)$

where m = 1/2∑k_ij is the total number of edges in the network, and B_ij = A_ij−k_ik_j/2m is the resultant modularity matrix, also known as graph Laplacian matrix. In such matrices, the optimisations can be achieved using graph partitioning or spectral partitioning (eigenvalues-eigenvectors decomposition) of the matrix B Newman [29]; Leicht and Newman [57].

Figure 4 shows the communities for the augmented network as determined by the optimisation of the Q function. Using structural brain segmentation and following anatomical partitioning, the resulting communities can be assigned to distinct brain areas. Specifically, six communities emerged from the simulation as populations with 294, 473, 189, 399, 290, and 330 members, all located in MC (see also Figure 7). Three important communities were detected in BG and the thalamus, the first one with 780 members in the thalamus, the second with 293 members connecting GPi and thalamus, and the third with 283 members connecting STN and GPe (see also Figure 9). In addition, one community with 41 members was obtained in the STN itself, and further two communities with 184 and 107 members connecting MC and STN as hyperdirect patwhway (see also Figure 10). The localisation of these communities in virtual space is given in Figure 4, as projection onto the MNI coordinate space in Figure 5.

FIGURE 4

FIGURE 4. Community detection using the modularity-index algorithm Leicht and Newman [57]. The algorithm identifies 12 major areas with the number of members of densely connected communities higher than 100. Remarkably, the modularity partition generated by the simulation is in line with anatomical brain separation.

FIGURE 5

FIGURE 5. Community detection using the modularity-index algorithm Leicht and Newman [57], now projecting the same communities as in Figure 4 onto MNI space sections, showing motor cortical areas as surface projection (A), thalamus and basal ganglia as frontal cut projection (B) and STN-MC connections as hyperdirect pathway projected onto a frontal section (C). As in the previous figure, the 12 major areas are located in MC (6), thalamus (1) and basal ganglia (2), as well as three communities in locations associated with the hyperdirect pathway (one in STN in light orange, two connecting STN and MC; dark orange and green).

4 Deep brain stimulation at centrality nodes: Electric field approximation

As outlined above, the centrality nodes (defined by the higher values of the “betweenness centrality” B_c and the “eigencentrality” E_c, as well as the clustering coefficient c) can be hypothesised to be possible stimulation targets especially effective in neuromodulation. To evaluate this hypothesis, for each of the three measures (B_c, E_c and c), first the centre of mass was calculated for the four positions in or close to the STN with the highest values (see coordinates as given in Table 1). For centre of mass positions, next an approximation of the volume of tissue activated (VTA) for a conventional DBS signal (2 mA 90 μs rectangular pulse with a 130 Hz repetition rate) was calculated. The stimulation was conducted in a monopolar mode, with the active contact placed at the coordinates of the centrality nodes, and the approximation was based on the electric field magnitude Åström et al. [58], thresholded at 150 V/m. Using the simulation platform OSS-DBS Butenko et al. [24], the field was computed in a heterogeneous and anisotropic volume conductor defined with data from Zhang and Arfanakis [59] and Horn [60]; Horn et al. [61].

Since the centrality nodes were located very close to each other, the three corresponding VTAs overlapped strongly; all three being located in the dorsolateral STN (Figure 6A), which is a clinically established target for treating motor symptoms of Parkinsons’s disease Benabid et al. [62]. It must be, however, noted that the employed structural connectivity of the network (Figure 1) was inherently biased towards the dorsolateral region. Next, we estimated the structural connectivity of these VTAs using the previously described pathway atlas, but now also including fibers beyond the motor circuit. In all three cases, nearly the same fibers were “recruited” by the stimulation, namely, the motor pallido-subthalamic projections, the hyperdirect pathway descending from the primary motor cortex (upper and lower extremity), and the dorso-lateral prefrontal cortex. Importantly, for the given stimulation amplitude, a recruitment of the corticofugal pathway was not predicted, allowing to avoid capsular side-effects Tommasi et al. [63]; Xu et al. [64]. Beyond this, no activation in the pallidothalamic projections was observed.

FIGURE 6

FIGURE 6. Recruitment of neural tissue by monopolar STN-DBS at the centrality nodes. (A): All three (overlapping) VTAs are located in the dorsolateral STN (left; green for between-centrality nodes, yellow for clustering coefficient, and blue eigencentrality) and predominantly recruit motor pallidosubthalamic projections, the hyperdirect pathway from the primary motor cortex (upper and lower extremities), as well as its branch descending from the dorso-lateral prefrontal cortex (right). Note that the latter was not used to construct the network model. For the particular stimulation protocol (2 mA), no corticofugal fibers are recruited according to the computational model. (B): The VTAs (green for between-centrality nodes, yellow for clustering coefficient, and blue eigencentrality) overlap with target areas, whose stimulation is clinically proven to improve motor symptoms, especially akinesia Dembek et al. [25]. While an overlap also with regions producing possible side-effects exists, this overlap is significantly smaller (note that the regions delineated in Dembek et al. [25] overlap as well).

Noteworthy is the spatial relation of the centrality nodes obtained by our simulation to the target spots of the STN-DBS. The VTAs significantly overlapped with STN regions shown to be effective in treating hypokinetic symptoms of Parkinsons’s disease, while a region implicated in side-effect occurrence was largely avoided Dembek et al. [25] (Figure 6B). Moreover, the VTAs of the present study contained the effective target points which were determined by projecting actual target coordinates of patients treated successfully with DBS for Parkinsons’s disease onto MNI space in another study on a large cohort Horn et al. [65]. Such a coincidence of the centrality nodes and the so-called sweet spots might explain the efficiency of the STN as a target for various neurological disorders. This relatively small nucleus is a site of convergence of various neural circuits (even though not all of them include the STN itself). Hence, its stimulation allows a wide-spread neuromodulatory intervention.

5 Modelling the dynamics of the thalamo-cortical basal ganglia circuit

After establishing the structural components and determining promising target regions in the previous sections on structural simulation, we also wished to explore the functional consequences of DBS on network activity patterns. For this, we worked on the same structural model as above, and superimposed modified Hodgkin-Huxley modelling.

In the augmented network, each node serves as a neuron (hypothesising a homogeneous neural population on the 1 mm³ cube), and the edges represent synaptic links between the neurons. Consequently, depending on the region, each node-neuron is modelled with a variation of Hodgkin-Huxley’s current-balance equations Terman et al. [28]; Rubin and Terman [9]; Hodgkin and Huxley [66]. In this section, we present the mathematical description of the neurons from each area of the basal ganglia (BG), thalamus and cortex. Next, we couple the neural activity of neurons according to structural connectivity within and between the subthalamic nucleus (STN), globus pallidus externa (GPe) and interna (GPi), thalamus (Tha) and motor cortex (MC).

5.1 Modelling and simulations of neurons in STN-GPe-GPi nuclei

The properties of neurons are expressed using the conductance-based biophysical model of Hodgkin-Huxley’s formalism as also has been used in previous work of Terman et al. [28]; Rubin and Terman [9]; Popovych and Tass [67]. The dynamics of each STN, GPe and GPi neuron are given by a current balance equation for the membrane potential: Terman et al. [28]; Bevan and Wilson [68]; Popovych and Tass [10]:

$C\frac{d{V}_i}{dt}=-I_{\text{LEAK}}-I_{\text{K}}-I_{\text{Na}}-I_{\text{Ca}}-I_{\text{T}}-I_{\text{AHP}}-I_{\text{syn}}+I_{\text{DBS}}(10)$

$\frac{d{x}_i}{dt}=\left(x_{\infty }-x_i\right)/{\tau }_{x_i}(11)$

$\frac{d{\left[{\text{Ca}}^{2+}\right]}_i}{dt}={ϵ}_1\left(-I_{\text{Ca}}-I_{\text{T}}-k_{\text{Ca}}{\left[{\text{Ca}}^{2+}\right]}_i\right),(12)$

where C is the membrane capacity, V_i is the membrane potential of the ith neuron, x_i denotes the gating variables n, h, r and x_∞ is the steady state value for the gating variables. The quantity ${[{\text{Ca}}^{2+}]}_i$ is the intracellular concentration of calcium. The exact description of the ionic currents I_LEAK, I_K, I_Na, I_Ca and the synaptic current I_syn for the STN and GP neurons are given in the supplementary material. The current I_DBS in Eq. 10 models DBS on STN neurons only and the form is described in the Supplementary Material. In the absence of DBS treatment the value is: I_DBS = 0.

5.2 Synaptic suppression of pallido-thalamic projections

Here, we model the GABAergic short term depression. The functionality of the GABAergic synapse and the resulting release of neurotransmitters is dependent on the firing history of the presynaptic neuron Zucker and Regehr [69]; Farokhniaee and McIntyre [70]. High-frequency stimulation induces suppression of GPi GABAergic synaptic transmission Farokhniaee and McIntyre [70], which in turn, leads to a thalamic activity facilitation.

The synaptic activity is defined by the activation variable s_i, which is given by Laing and Chow [71]; Ermentrout and Terman [72]; Compte et al. [73]:

$\frac{d{s}_i}{dt}=\alpha \left(1-s_i\right)H\left(V_i-{\theta }_0\right)-\beta s_i,(13)$

where H(V) is a smooth approximation of the step function, i.e., $H(V)=1/(1+e^{-(V-{\theta }_x)/{\sigma }_x})$, where α, β stands for the rate of activation and inactivation, respectively, and typically, α = O(1), β = O(ϵ) holds Laing and Chow [71]; Ermentrout and Terman [72]; Terman [74].

The inhibitory GaBAergic synaptic current for the ith neuron is given, by

$I_{i,\text{GABA}}=g_{\text{GABA}}\left(V_i-E_{\text{GABA}}\right)\sum _j{A}_{ij}{s}_j,(14)$

where A_ij has the value 1 or 0, depending on whether the neuron is connected or not. The summation is taken over all presynaptic neurons.

In case of existent synaptic suppression the GABAergic synaptic current changes to

$I_{i,\text{GABA}}=g_{\text{GABA}}\left(V_i-E_{\text{GABA}}\right)\sum _j{A}_{ij}{s}_j{P}_j,(15)$

where the factor P_j describes the probability of a neurotransmitter release (in the {ij} synapses), and follows the dynamics Benita et al. [75]:

$\begin{array}{c}\hfill \frac{d{P}_j}{dt}=\frac{P_0-P_j}{{\tau }_D}\\ \hfill P_j\left(t_{sp}\right)\to P_j\left(t_{sp}\right)A_D,\end{array}(16)$

where t_sp corresponds to the last spike-time of the presynaptic neuron and A_D is the depression factor (0 < A_D < 1), in our case, the value A_D = 0.8 was used. The value P₀ describes the steady state of P and in our case was set to 1. To simplify, when a presynaptic neuron fires at time t_sp the functionality of the synapse (the release of neurotransmitters) is reduced (suppressed by a factor A_D). In the absence of neural activity the synapse returns to a full ability of release, in a time scale 1/τ where τ = 400 ms Benita et al. [75].

5.3 Modelling neurons in the thalamus

The mathematical description of the thalamic neurons is given by the following equation.

$C\frac{d{V}_i}{dt}=-I_{\text{LEAK}}-I_{\text{K}}-I_{\text{Na}}-I_{\text{T}}-I_{\text{syn}}+I_{\text{SM}}(17)$

$\frac{d{x}_i}{dt}=\left(x_{\infty }-x_i\right)/{\tau }_x,(18)$

where C is the membrane capacity and V_i is the membrane potential of the ith neuron, while the Eq. 18 describes the first order kinetics for the gating variables h, r. The currents I_LEAK, I_K and I_Na are the ionic currents, I_T is the T-type calcium channel. The synaptic current I_syn has the form I_syn = I_GPTH + I_THTH, where the GABAergic current I_GPTH represents the inhibition of the GPi area to the thalamus, while I_THTH represents the internal excitatory or inhibitory thalamic connections. The current I_SM represents sensorimotor excitation (from motor cortex areas to thalamus). The detailed description of the ionic and synaptic currents is given in the Supplementary Material.

5.4 Modelling and simulations of neurons in the motor cortex

The motor cortex neurons MC, are described as one-compartment soma, and following the equations Pospischil et al. [76]:

$C\frac{d{V}_i}{dt}=-I_{\text{LEAK}}-I_{\text{K}}-I_{\text{Na}}-I_{\text{M}}-I_{\text{syn}}+I_{\text{app}}(19)$

$\frac{d{x}_i}{dt}=a_x\left(1-x_i\right)-b_x{x}_i(20)$

$\frac{d{p}_i}{dt}=\left(p_{\infty }-p_i\right)/{\tau }_p,(21)$

where V_i is the membrane potential, and x_i represents the gating variables for potassium and sodium current, of the ith neuron. The gating variable p_i represents the activation gate of the slow, voltage-dependent potassium current I_M. The current I_app is added to tune the oscillatory behaviour of MC neurons around 20Hz. Each MC neuron has different value of I_app which is extracted randomly from the interval [2, 3]. The synaptic activity is given from the current I_syn and the exact form is described at Supplementary Material. The whole MC area is modelled as small world network. In this network, 20% of the neurons send inhibitory signals. i.e., replicate interneurons. The cortical neurons show a regular spiking activity Pospischil et al. [76]. The exact description is given in the Supplementary Material.

5.5 Average over the detected communities macroscopic description

In this paper, we obtain a macroscopic description of the dynamics of the detected communities of Section 3.4, by averaging the mean voltage activity $\bar{V}$ of neurons over the population in the community; specifically, we define:

${\bar{V}}_x\left(t\right)=\frac{1}{N}\sum _{k=1}^N{V}_k\left(t\right).(22)$

The mean voltage activity $\bar{V}$ is used for the characterisation of rhythmic activity using Fourier spectral analysis under different states (Healthy, Parkinsonian or DBS in Parkinsonian conditions). In all simulations, the Fourier power spectrum is normalised dividing by the highest absolute value. In this context, the Parkinsonian state was modelled, in brief, by increasing the activity of STN, decreasing the activity of GPe (D2 dopamine-mediated receptor activity effect in the indirect pathway), and simultaneously increasing the activation of GPi due to D1 dopamine-mediated activity in the direct pathway. The detailed description of is given in Section 3 of the Supplementary Material.

6 Collective dynamics of the structural clusters (communities)

Structural connectivity can have significant impact on the large-scale dynamics of the brain Deco et al. [43]; Papadopoulos et al. [77]; Deco et al. [15], However, the connection between anatomical-structural and functional brain connectivity is far from been trivial. Large-scale computational models and their complex nonlinear dynamics constitute an important method to explore this connection Papadopoulos et al. [77]; Schirner et al. [78]. Here, we propose a new method which correlates the structural and functional connectivity, specifically focusing on densely connected communities as identified by the modality index (Section 3.4).

6.1 Analysis of macroscopic activity in motor cortex clusters

The network analysis resulted in the identification of 6 MC areas consisting of 294, 473, 189, 399, 290 and 330 nodes, respectively. For each area, we extract the Fourier spectrum for the macroscopic variable of Eq. 22, averaging all values of each member of the detected community. The results are depicted in Figure 7. The first column depicts the six different cortical node sets (six clusters emerging from modularity analysis) in red on the virtual brain surface. The next three columns show the power spectra in three different conditions (i.e., Healthy, Parkinsonian and DBS). Under healthy conditions, the activity peaks in the low and high γ band (i.e. ≈ 80 Hz and ≈ 190–290 Hz, with one to three peaks). Under Parkinsonian conditions, these peaks are often blunted, i.e., with less obvious peaks in e.g., MC1, MC3 and MC6 (Parkinsonian conditions depicted as blue curves, and healthy conditions depicted as red curves, for comparison see also column 1 of Figure 8 which depicts the differences between Healthy and Parkinsonian cases).

FIGURE 7

FIGURE 7. Simulation of frequency spectrograms of cortical sets of nodes as average neuronal activity (power on ordinate given as relative value) in a frequency range of 0–400 Hz. The six rows depict six different cortical node sets (six clusters emerging from modularity analysis) shown in red on the virtual brain surface. The spectrograms correspond to healthy conditions (red traces), Parkinsonian conditions and DBS in a Parkinson-affected network (both as blue traces), as shown in the different columns. Traces from healthy conditions are superimposed as fine red traces. Under healthy conditions, the activity peaks in the low and high γ band (i.e., ≈80 Hz and ≈190–290 Hz, with one to three peaks). Under Parkinsonian conditions, these peaks are generally blunted, and DBS is able to reverse this at least in some instances (e.g., clusters in MC3, MC4 and MC5, where spectrograms under DBS deviate less from healthy conditions than spectrograms under Parkinsonian conditions by as much as 17–28%, taking the overall area under the curve differences as a measure; see Figure 8).

FIGURE 8

FIGURE 8. Differences of the spectrum between Park-Healthy and DBS-Healthy [subtracting the spectrum according to Eq. 23]. Each row represents one of six cortical sets of nodes as emergent from centrality measure modelling. The columns show the differences of these frequency spectra by subtracting Park vs. Healthy, DBS vs. Park and DBS vs. Healthy. The shaded area under the curve as a measure of frequency spectrum divergence is computed and depicted in columns 1, 2 and 3. The MC3, MC4 and MC5 areas in the cases of DBS vs. Healthy show a reduction in the computed area compared to Park vs. Healthy, indicating that DBS reduces the spectrum difference compared to Parkinsonian condition, although differences to the healthy condition remain.

To better estimate the differences between the conditions, we generated difference curves of the spectrograms of Figure 7 based on the following subtraction pairs: |P_X(f)|−|P_Y(f)|, where |P_X|, |P_Y|, are the power spectra of the states X, Y, respectively, and represent: healthy, Parkinsonian or DBS conditions, see Figure 8. By further calculating the area under the curve (AUC, Valor et al. [79]) for each of these difference pairs i.e.,

$E={\int }_a^b\left||P_X\left(f\right)|-|P_Y\left(f\right)|\right|df(23)$

(numbers in insets in Figure 8) one can estimate the degree of change expected to occur when changing the condition. Comparing the difference between healthy and Parkinsonian conditions [Park-Healthy, Eq. 23], one can see that the spectrograms differ by $\approx 55$ arbitrary units for most MC clusters, save in two locations (MC2 and MC6), where the values are below 50. Obviously, thus, the effect of Parkinsonian conditions is heterogeneous in the MC network, albeit within moderate boundaries; overall, the difference is in the range of 50 ± 5 (arbitrary units).

Comparing these differences now to differences DBS-Healthy, see Eq. 23, the relative effect of DBS can be gauged: In two areas, DBS actually induced more differences than the disease condition alone (MC1, MC2, with values 62 and 49, compared to 55 and 45), in three areas, DBS reduces the differences (which might be interpreted as a normalisation of activity), i.e., in MC3, MC4 and MC5 (with values 47, 40 and 44 compared to 58, 55 and 53), and in one area, there is virtually no effect of DBS regarding this measure (MC6, with a value 49, compared to 49). Again, this leads to the conclusion that the effect of DBS is heterogeneous regarding cortical activity, with alterations by +18 and +8% (positive changes meaning that the frequency spectrogram digresses even more from healthy conditions under DBS than under Parkinsonian conditions alone) occurring in some regions (MC1, MC2), and −20, −28 and −17% in others (MC3, MC4, MC5; negative values indicating that the spectrograms under DBS show less of a difference against healthy conditions than under Parkinsonian conditions without DBS), and in fact only a minimal change (+0.9%) in MC6, see Figure 8 (Supplementary Table S4). Using repetitive analyses with the data, we did not see changes of $>5\%$, while the effect sizes particularly of beneficial DBS effects in the order of 17%–28%. Therefore, we think that the differences are unlikely due to random sampling errors. Overall, just averaging these changes, this amounts to a change of cortical spectral activity by $\approx -7\%$. Where are the clusters positioned on the cortex? With MC1 and MC2, those regions where DBS seems to accentuate differences in spectral activity, we see the largest clusters forming (incidentally) a “W” or “∞” figure. With MC3, MC4 and MC5, it appears that these clusters are positioned very close to the frontal or dorsal end of the MC2 cluster, or indeed at fragmented positions of the MC1 cluster; those regions interestingly show a reduction in spectrogram difference induced by DBS. One might speculate that normalising effects of DBS are reflected in spatially fragmented changes in the cortex, but not when considering large networks. Overall, obviously, the reason for the heterogeneity in functional connectivity or for the impact of DBS remains unknown and we hope that the current analysis spurs further detailed investigations into this matter.

6.2 Analysis of macroscopic activity in basal ganglia-thalamic clusters

The second group of communities (clusters), identified by modularity measure, belongs to the basal ganglia or thalamus. Specifically, the first thalamic cluster contains 780 neurons, the second community contains 293 neurons connecting GPi and thalamus and the third cluster consists of 283 neurons connecting STN and GPe. For each cluster, we extract the Fourier spectrum for the macroscopic variable of Eq. 22. These Fourier spectra of all members of the detected communities were averaged. The results are depicted in Figure 9. The three rows depict node sets emerging from modularity analysis in the thalamus (top), in the GPi-thalamic pathway (middle) and the STN-GPe pathway (bottom), depicted as red dots on the virtual brain sections.

FIGURE 9

FIGURE 9. Simulation of frequency spectrograms of basal ganglia and thalamic sets of nodes as average neuronal activity (power on ordinate given as relative value) in a frequency range of 0–400 Hz. The three rows depict node sets emerging from modularity analysis in the thalamus (top), in the GPi-thalamic pathway (middle) and the STN-GPe pathway (bottom), depicted as red dots on the virtual brain sections. The spectrograms correspond to healthy conditions (red traces), Parkinsonian conditions and DBS in a Parkinson-affected network (both as blue traces), as shown in the different columns. Traces from healthy conditions are superimposed as fine red traces. Under healthy conditions, the activity peaks in the low γ band (i.e., ≈ 46 Hz) in the thalamus, while in the pathways projecting from GPi to the thalamus and projecting from the STN to GPe, low frequency activity in the low β band is seen (i.e. ≈ 13–14 Hz). Under Parkinsonian conditions, the situation reverses: in the thalamus, higher frequency γ activity is blunted, and low-frequency activity emerges (at ≈6 Hz, i.e., θ band), while in the cited pathways, peak frequencies are shifted to higher frequencies (≈25 and 50 Hz, i.e., in the β and low γ band). DBS clearly changes this situation: In the pathways directly connected to the STN (which is the stimulation target), the dominant frequencies are in the DBS frequency and harmonics (i.e., 130 and 260 Hz). This leads to a restoration of thalamic activity, where again the activity is peaking in the low γ band (i.e., ≈38 Hz).

Under healthy conditions (column 2 in Figure 9), in the thalamus, the activity peaks in the low γ band (i.e., ≈ 46 Hz). In the pallido-thalamic cluster (projections from GPi to thalamus; GPi-Thal), the neuronal activity is characterised by a slow β band rhythm (i.e., ≈ 13–14 Hz). The third cluster comprises neurons projecting from the STN to GPe (STN-GPe), again showing a maximum of activity in the low β band (i.e., ≈ 13–14 Hz).

Under Parkinsonian conditions, the situation reverses: in the thalamus, higher frequency γ activity is blunted, and low-frequency activity emerges (at ≈ 6 Hz, i.e., θ band), while in the cited pallidal pathways, peak frequencies are shifted to higher frequencies (≈25 and 50 Hz, i.e., in the β and low γ band). DBS clearly changes this situation: In the pathways directly connected to the STN (which is the stimulation target), the dominant frequencies are in the DBS frequency and harmonics (i.e., 130 and 260 Hz). Importantly, this leads to a restoration of thalamic activity, where again the activity is peaking in the low γ band (i.e., ≈ 38 Hz), not quite reaching the frequency under healthy conditions, but definitely different from low-frequency Parkinsonian activity.

6.3 Analysis of macroscopic activity in clusters associated with the hyperdirect pathway

The third group of communities emerging from the simulation were clusters associated with the hyperdirect pathway, as shown in Figure 10. These include one cluster in the STN, and two clusters connecting MC and STN. The first cluster is made of 41 neurons, and the other two of 184 and 107 neurons. Under healthy conditions, in a set of STN nodes quite different from the one projecting to the GPe of the BG-thalamic clusters (compare with Figure 9), the activity is also distinctly different as it peaks at 21 Hz (compared to 14 Hz in the BG-cluster), with a much wider frequency distribution into higher frequencies. These nodes probably relate to the hyperdirect pathway, as they are close to the nodes of MC-STN connections. In the latter, under healthy conditions, peak activity emerges at ≈ 160–180 Hz, i.e. in the high γ range.

FIGURE 10

FIGURE 10. Simulation of frequency spectrograms of the STN (top row) and the hyperdirect pathway connections (two bottom rows) as average neuronal activity (power on ordinate given as relative value) in a frequency range of 0–400 Hz. The three rows thus depict node sets emerging from modularity analysis in the STN itself (top), and in nodes connecting MC and STN, i.e., two clusters of connections involving the hyperdirect pathway (middle and bottom), depicted as red dots on the virtual brain sections. The spectrograms correspond to healthy conditions (red traces), Parkinsonian conditions and DBS in a Parkinson-affected network (both as blue traces), as shown in the different columns. Traces from healthy conditions are superimposed as fine red traces. Under healthy conditions, in a set of STN nodes quite different from the one projecting to the GPe of the previous figure, the activity is also distinctly different as it peaks at 21 Hz, with a very wide frequency distribution into higher frequencies. These nodes probably relate to the hyperdirect pathway, as they are close to the nodes of the second row, i.e., MC-STN connections. In the latter, under healthy conditions, peak activity emerges at ≈ 160–180 Hz, i.e., in the high γ range. Under Parkinsonian conditions, STN nodes narrow down their spectrum to two peaks (at 44 and 86 Hz), and a loss of the wide spectral activity. Under DBS, in the STN nodes a wider distribution of frequencies appears again, albeit as an appearance of peaks at harmonic frequencies of 44 Hz (44, 88, 172 Hz). In the hyperdirect pathway related nodes, however, DBS does not change much in the activity, save appearance of a 130 Hz peak (then stimulating frequency). In the first cluster (middle), the frequency distribution otherwise remains more or less similar, and in the second cluster (bottom), the 160 Hz peak is blunted.

Under Parkinsonian conditions, STN nodes narrow down their spectrum to two peaks (at 44 and 86 Hz), and a loss of the wide spectral activity. Under DBS, in the STN nodes a wider distribution of frequencies appears again, albeit as an appearance of peaks at harmonic frequencies of 44 Hz (44, 88, 172 Hz). In the hyperdirect pathway related nodes, however, DBS does not change much in the activity, save appearance of a 130 Hz peak (then stimulating frequency). In the first cluster (see middle row in Figure 10), the frequency distribution otherwise remains more or less similar, and in the second cluster (bottom row in Figure 10), the 160 Hz peak is blunted.

7 Discussion and conclusion

In this study, we provided insights into complex network processes in order to obtain a new approach gauging the effectiveness of DBS. Additionally, we investigated the relationship between structural and functional connectivity, presenting a new in silico methodological approach to explore dynamics of brain motor area functions under healthy and Parkinsonian conditions, as well as the impact of Deep Brain Stimulation (DBS). Using a state-of-the-art network (constructed from data based on the pathway atlas of human motor network, obtained from various types of imaging, including diffusion, histological and structural MRI data, all fused to a virtual 3D rendering, Petersen et al. [36]) and integrating this network into advanced complex network measures, Bullmore and Sporns [18], we detected nodes, with high-connectivity and thus pivotal impact on the activity distribution within the network. These nodes are hypothesised to be ideal targets for DBS application. Our results based on betweenness centrality propose the following MNI coordinates (13.5, −12.8, −6.3) as optimal STN target. The following optimal points (sweet spots) are suggested in other publications: (12.5, −12.72, −5.38) in Dembek et al. [25]; (12.42, −12.58, −5.92) in Horn et al. [65]; (11.83, −11.63, −5.8) in Bot et al. [80] and (10.83, −13.31, −7.01) in Akram et al. [81], see also Table 2 in Dembek et al. [25]. Remarkably, completely different methodologies [e.g., we use graph theory, while Dembek et al. [25] use VTA with Probabilistic Stimulation Maps (PSM)] result in very similar STN sweet spots. The advantage of using a graph-theoretical approach, as in the current study, is that it is simple and takes only a few seconds to estimate the best stimulus target. Furthermore, this method could be combined in the future with already established statistical methods used e.g., Dembek et al. [25].

As a next step, we computed the volume of tissue activated (VTA) at the positions around the STN, which had emerged as pivotal nodes. Comparing these VTAs to clinically established DBS targets in Parkinsons’s disease, it is evident that the position of the nodes matches well with areas associated with alleviation of motor symptoms Benabid et al. [62]; Dembek et al. [25]; Horn et al. [65], suggesting that the high-connectivity nodes can infer potentially effective stimulation sites. Future studies should investigate whether a match between such nodes and neurosurgical targets also occurs when analysing networks constructed based on whole brain structural connectomes.

The second part of the current study addresses modelling of functional changes within a neuronal network. Specifically, here, we propose a new method to analyse network activity under different conditions (i.e., Healthy, Parkinsonian and DBS), using knowledge of detailed structural network connectivity in a large scale Basal ganglia-thalamo-cortical model. This new approach is based on the detection of network communities or modules central to activity distribution. Network analysis of structural connectivity showed that the communities or groups of highly-connected nodes can be assigned to distinct anatomical regions (Figure 4). Such a modular network organisation, in comparison to a random distribution network, clearly shows advantages like greater robustness, adaptivity, and evolvability of network function Meunier et al. [82]. Although the relation of structural/functional connectivity and how neural activity could emerge from the brain’s anatomical connections has been studied in several other experimental and computational studies Deco et al. [16,43]; Horn et al. [83], the current view of modular activity organisation is new.

Our analysis showed that in all modular areas, including neocortical ones, Parkinsonian conditions alter power spectrogams, but mainly subcortical structures, with e.g., slowing of activity in the thalamus and faster activity in pathways connecting pallido-thalamic and subthalamic-pallidal nodes. Under DBS, in turn, simulations reveal that this stimulation at high-connectivity nodes is able to restore thalamic activity, and partly also cortical one, while on the other hand, hyperdirect-pathway associated nodes remain largely unaffected. Specifically, the simulations suggest that nuclei directly involved in DBS (STN, Pallidum) mainly follow the stimulation. The thalamus, in turn, translates this into a concordant shift of its activity to the low-γ frequency band (from θ under Parkinsonian conditions), while the motor cortex, in turn, shows a discrete, and inhomogeneous response. Thus, theoretically, the thalamus also under these conditions may serve as activity gate to the cortex, while the motor cortex only adjusts in a minor way—presumably thus preserving general functionality, which would likely be lost if strong rhythmicity were to emerge in the cortex. Importantly, at least a part of these conclusions is also clinically confirmed. In a recent study Neumann et al. [84] on Parkinsons’s disease patients receiving DBS, using clinical, behavioural and fiber tracking informed computational models, the hypokinetic state depends on suppressing indirect pathway activity and not on the hyperdirect pathway. By contrast, in that study, cognitive impairment in Parkinson’s disease patients could be attributed to modulation of the hyperdirect pathway, suggesting that the hyperdirect and indirect pathways, converging in the subthalamic nucleus, are differentially involved in cognitive aspects of motor programming and kinematic gain control during motor performance.

The present study constitutes a computational approximation of the basal ganglia-thalamo-cortical network, with assumptions and limitations. Regarding the assumptions, synaptic coupling was tuned to be consistent to produce beta-band oscillatory activity within the basal ganglia, as a pathophysiological marker of Parkinson’s disease. Further, an internal connectivity in the nuclei was assumed to take the form of small world complex structures. This novel approach in basal ganglia modelling has a reasonable justification in previous publications, both modelling and experimental Netoff et al. [49]; Berman et al. [50]; She et al. [51]; Bassett and Bullmore [17,47]; Fang et al. [52]; De Santos-Sierra et al. [20]. As a limitation of the model, the exact structure of the connectivity on this microscopic level is not known, and hence also it remains to be clarified in the future how this can be analogously modelled. In the current study, the striatal input to GPe and GPi was simplified as different, but homogeneous constant currents to all neurons in GPe and GPi, but with different values between GPe and GPi, as well as healthy and Parkinsonian conditions.

As a future modelling perspective, one important topic will be the investigation of how this structural separation of the network can help to construct a low-representation model of brain activity (also called neural manifolds). The low representation will be tested under several DBS variations (i.e., functional connectivity) with respect to the parameters of DBS implantation (position, frequency, shape).

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

KS, RK, and JS contributed to conceptualisation, methodology, model analysis and investigation. KS and KB contributed to model simulations and investigation. RK, JS, and UvR reviewing and editing the manuscript, and supervision. All authors contributed to the article, writing of the original draft, and approved the submitted version.

Funding

This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—SFB 1270/2-299150580.

Conflict of interest

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

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.951724/full#supplementary-material

References

1. Kühn A, Kempf F, Brücke C, Doyle L, Martinez-Torres I, Pogosyan A, et al. High-frequency stimulation of the subthalamic nucleus suppresses oscillatory β activity in patients with Parkinson’s disease in parallel with improvement in motor performance. J Neurosci (2008) 28:6165–73. doi:10.1523/jneurosci.0282-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Neumann W-J, Staub-Bartelt F, Horn A, Schanda J, Schneider G-H, Brown P, et al. Long term correlation of subthalamic beta band activity with motor impairment in patients with Parkinson’s disease. Clin Neurophysiol (2017) 128:2286–91. doi:10.1016/j.clinph.2017.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Kim J, Kim Y, Nakajima R, Shin A, Jeong M, Park A, et al. Inhibitory basal ganglia inputs induce excitatory motor signals in the thalamus. Neuron (2017) 95:1181–96.e8. e8. doi:10.1016/j.neuron.2017.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Galvan A, Devergnas A, Wichmann T. Alterations in neuronal activity in basal ganglia-thalamocortical circuits in the parkinsonian state. Front Neuroanat (2015) 9:5. doi:10.3389/fnana.2015.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Deuschl G, Schade-Brittinger C, Krack P, Volkmann J, Schäfer H, Bötzel K, et al. A randomized trial of deep-brain stimulation for Parkinson’s disease. N Engl J Med Overseas Ed (2006) 355:896–908. doi:10.1056/nejmoa060281

CrossRef Full Text | Google Scholar

6. Vidailhet M, Jutras M-F, Grabli D, Roze E. Deep brain stimulation for dystonia. J Neurol Neurosurg Psychiatry (2013) 84:1029–42. doi:10.1136/jnnp-2011-301714

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Crompe B, Aristieta A, Leblois A, Elsherbiny S, Boraud T, Mallet N. The globus pallidus orchestrates abnormal network dynamics in a model of parkinsonism. Nat Commun (2020) 11:1570. doi:10.1038/s41467-020-15352-3

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Xu W, Russo G, Hashimoto T, Zhang J, Vitek J. Subthalamic nucleus stimulation modulates thalamic neuronal activity. J Neurosci (2008) 28:11916–24. doi:10.1523/jneurosci.2027-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Rubin J, Terman D. High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J Comput Neurosci (2004) 16:211–35. doi:10.1023/b:jcns.0000025686.47117.67

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Popovych O, Tass P. Adaptive delivery of continuous and delayed feedback deep brain stimulation - a computational study. Sci Rep (2019) 9:10585. doi:10.1038/s41598-019-47036-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. So R, Kent A, Grill W. Relative contributions of local cell and passing fiber activation and silencing to changes in thalamic fidelity during deep brain stimulation and lesioning: A computational modeling study. J Comput Neurosci (2012) 32:499–519. doi:10.1007/s10827-011-0366-4

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Galvan A, Wichmann T. Pathophysiology of parkinsonism. Clin Neurophysiol (2008) 119:1459–74. doi:10.1016/j.clinph.2008.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Spiliotis K, Siettos C. A timestepper-based approach for the coarse-grained analysis of microscopic neuronal simulators on networks: Bifurcation and rare-events micro- to macro-computations. Neurocomputing (2011) 74:3576–89. doi:10.1016/j.neucom.2011.06.018

CrossRef Full Text | Google Scholar

14. Siettos C, Starke J. Multiscale modeling of brain dynamics: From single neurons and networks to mathematical tools. WIREs Mech Dis (2016) 8:438–58. doi:10.1002/wsbm.1348

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Deco G, Jirsa V, Robinson P, Breakspear M, Friston K. The dynamic brain: From spiking neurons to neural masses and cortical fields. Plos Comput Biol (2008) 4:e1000092. doi:10.1371/journal.pcbi.1000092

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Deco G, Ponce-Alvarez A, Mantini D, Romani GL, Hagmann P, Corbetta M. Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J Neurosci (2013) 33:11239–52. doi:10.1523/jneurosci.1091-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bassett DS, Bullmore ET. Small-world brain networks revisited. Neuroscientist (2017) 23:499–516. doi:10.1177/1073858416667720

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Bullmore E, Sporns O. Complex brain networks: Graph theoretical analysis of structural and functional systems. Nat Rev Neurosci (2009) 10:186–98. doi:10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Iliopoulos A, Papasotiriou I. Functional complex networks based on operational architectonics: Application on eeg-based brain–computer interface for imagined speech. Neuroscience (2022) 484:98–118. doi:10.1016/j.neuroscience.2021.11.045

PubMed Abstract | CrossRef Full Text | Google Scholar

20. De Santos-Sierra D, Sendiña-Nadal I, Leyva I, Almendral J, Anava S, Ayali A, et al. Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures. PLoS ONE (2014) 9:e85828. doi:10.1371/journal.pone.0085828

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Crowell A, Ryapolova-Webb E, Ostrem J, Galifianakis N, Shimamoto S, Lim D, et al. Oscillations in sensorimotor cortex in movement disorders: An electrocorticography study. Brain (2012) 135:615–30. doi:10.1093/brain/awr332

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Spiliotis K, Starke J, Franz D, Richter A, Köhling R. Deep brain stimulation for movement disorder treatment: Exploring frequency-dependent efficacy in a computational network model. Biol Cybern (2021) 116:93–116. doi:10.1007/s00422-021-00909-2

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Butson CR, Cooper SE, Henderson JM, McIntyre CC. Patient-specific analysis of the volume of tissue activated during deep brain stimulation. NeuroImage (2007) 34:661–70. doi:10.1016/j.neuroimage.2006.09.034

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Butenko K, Bahls C, Schröder M, Kohling R, van Rienen U. OSS-DBS: Open-source simulation platform for deep brain stimulation with a comprehensive automated modeling. Plos Comput Biol (2020) 16:e1008023–18. doi:10.1371/journal.pcbi.1008023

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Dembek TA, Roediger J, Horn A, Reker P, Oehrn C, Dafsari HS, et al. Probabilistic sweet spots predict motor outcome for deep brain stimulation in Parkinson disease. Ann Neurol (2019) 86:527–38. doi:10.1002/ana.25567

PubMed Abstract | CrossRef Full Text | Google Scholar

26. van Hartevelt TJ, Cabral J, Deco G, Møller A, Green AL, Aziz TZ, et al. Neural plasticity in human brain connectivity: The effects of long term deep brain stimulation of the subthalamic nucleus in Parkinson’s disease. PLOS ONE (2014) 9:e86496–13. doi:10.1371/journal.pone.0086496

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Li C, Li Q, Van Mieghem P, Stanley HE, Wang H. Correlation between centrality metrics and their application to the opinion model. Eur Phys J B (2015) 88:65–13. doi:10.1140/epjb/e2015-50671-y

CrossRef Full Text | Google Scholar

28. Terman D, Rubin J, Yew A, Wilson C. Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J Neurosci (2002) 22:2963–76. doi:10.1523/jneurosci.22-07-02963.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Newman MEJ. Modularity and community structure in networks. Proc Natl Acad Sci U S A (2006) 103:8577–82. doi:10.1073/pnas.0601602103

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Mandali A, Chakravarthy VS, Rajan R, Sarma S, Kishore A. Electrode position and current amplitude modulate impulsivity after subthalamic stimulation in Parkinsons disease—A computational study. Front Physiol (2016) 7:585. doi:10.3389/fphys.2016.00585

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Fleming J, Dunn E, Lowery M. Simulation of closed-loop deep brain stimulation control schemes for suppression of pathological beta oscillations in Parkinson’s disease. Front Neurosci (2020) 14:166. doi:10.3389/fnins.2020.00166

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ewert S, Plettig P, Li N, Chakravarty MM, Collins DL, Herrington TM, et al. Toward defining deep brain stimulation targets in MNI space: A subcortical atlas based on multimodal MRI, histology and structural connectivity. NeuroImage (2018) 170:271–82. doi:10.1016/j.neuroimage.2017.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Chakravarty MM, Bertrand G, Hodge CP, Sadikot AF, Collins DL. The creation of a brain atlas for image guided neurosurgery using serial histological data. NeuroImage (2006) 30:359–76. doi:10.1016/j.neuroimage.2005.09.041

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Tian Y, Margulies D, Breakspear M, Zalesky A. Topographic organization of the human subcortex unveiled with functional connectivity gradients. Nat Neurosci (2020) 23:1421–32. doi:10.1038/s41593-020-00711-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Fan L, Li H, Zhuo J, Zhang Y, Wang J, Chen L, et al. The human brainnetome atlas: A new brain atlas based on connectional architecture. Cereb Cortex (2016) 26:3508–26. doi:10.1093/cercor/bhw157

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Petersen MV, Mlakar J, Haber S, Parent M, Smith Y, Strick P, et al. Holographic reconstruction of axonal pathways in the human brain. Neuron (2019) 104:1056–64.e3. e3. doi:10.1016/j.neuron.2019.09.030

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Marek K, Jennings D, Lasch S, Siderowf A, Tanner C, Simuni T, et al. The Parkinson progression marker initiative (ppmi). Prog Neurobiol (2011) 95:629–35. doi:10.1016/j.pneurobio.2011.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Milardi D, Quartarone A, Bramanti A, Anastasi G, Bertino S, Basile GA, et al. The cortico-basal ganglia-cerebellar network: Past, present and future perspectives. Front Syst Neurosci (2019) 13:61. doi:10.3389/fnsys.2019.00061

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Bosch-Bouju C, Hyland B, Parr-Brownlie L. Motor thalamus integration of cortical, cerebellar and basal ganglia information: Implications for normal and parkinsonian conditions. Front Comput Neurosci (2013) 7:163. doi:10.3389/fncom.2013.00163

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Tarnaud T, Joseph W, Martens L, Tanghe E. Dependence of excitability indices on membrane channel dynamics, myelin impedance, electrode location and stimulus waveforms in myelinated and unmyelinated fibre models. Med Biol Eng Comput (2018) 56:1595–613. doi:10.1007/s11517-018-1799-y

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lanciego JL, Luquin N, Obeso JA. Functional neuroanatomy of the basal ganglia. Cold Spring Harbor Perspect Med (2012) 2:a009621. doi:10.1101/cshperspect.a009621

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Stam C, Reijneveld J. Graph theoretical analysis of complex networks in the brain. Nonlinear Biomed Phys (2007) 1:3. doi:10.1186/1753-4631-1-3

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Deco G, Senden M, Jirsa V. How anatomy shapes dynamics: A semi-analytical study of the brain at rest by a simple spin model. Front Comput Neurosci (2012) 6:68–7. doi:10.3389/fncom.2012.00068

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Honey CJ, Sporns O, Cammoun L, Gigandet X, Thiran JP, Meuli R, et al. Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A (2009) 106:2035–40. doi:10.1073/pnas.0811168106

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Honey C, Thivierge J-P, Sporns O. Can structure predict function in the human brain? NeuroImage (2010) 52:766–76. doi:10.1016/j.neuroimage.2010.01.071

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Watts D, Strogatz S. Collective dynamics of ‘small-world’ networks. Nature (1998) 393:440–2. doi:10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Bassett D, Bullmore E. Small-world brain networks. Neuroscientist (2006) 12:512–23. doi:10.1177/1073858406293182

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Mark N. The structure and function of complex networks. SIAM Rev (2003) 45:58. doi:10.1137/S003614450342480

CrossRef Full Text | Google Scholar

49. Netoff T, Clewley R, Arno S, Keck T, White J. Epilepsy in small-world networks. J Neurosci (2004) 24:8075–83. doi:10.1523/jneurosci.1509-04.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Berman B, Smucny J, Wylie K, Shelton E, Kronberg E, Leehey M, et al. Levodopa modulates small-world architecture of functional brain networks in Parkinson’s disease. Mov Disord (2016) 31:1676–84. doi:10.1002/mds.26713

PubMed Abstract | CrossRef Full Text | Google Scholar

51. She Q, Chen G, Chan R. Evaluating the small-world-ness of a sampled network: Functional connectivity of entorhinal-hippocampal circuitry. Sci Rep (2016) 6:21468. doi:10.1038/srep21468

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Fang J, Chen H, Cao Z, Jiang Y, Ma L, Ma H, et al. Impaired brain network architecture in newly diagnosed Parkinson’s disease based on graph theoretical analysis. Neurosci Lett (2017) 657:151–8. doi:10.1016/j.neulet.2017.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Gouty-Colomer L-A, Michel F, Baude A, Lopez-Pauchet C, Dufour A, Cossart R, et al. Mouse subthalamic nucleus neurons with local axon collaterals. J Comp Neurol (2018) 526:275–84. doi:10.1002/cne.24334

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Ammari R, Lopez C, Bioulac B, Garcia L, Hammond C. Subthalamic nucleus evokes similar long lasting glutamatergic excitations in pallidal, entopeduncular and nigral neurons in the basal ganglia slice. Neuroscience (2010) 166:808–18. doi:10.1016/j.neuroscience.2010.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

55. van den Heuvel MP, Sporns O. Rich-club organization of the human connectome. J Neurosci (2011) 31:15775–86. doi:10.1523/JNEUROSCI.3539-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Newman M. Networks, an introductions. New York: Oxford University Press (2010).

Google Scholar

57. Leicht EA, Newman MEJ. Community structure in directed networks. Phys Rev Lett (2008) 100:118703. doi:10.1103/PhysRevLett.100.118703

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Åström M, Diczfalusy E, Martens H, Wårdell K. Relationship between neural activation and electric field distribution during deep brain stimulation. IEEE Trans Biomed Eng (2015) 62:664–72. doi:10.1109/TBME.2014.2363494

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Zhang S, Arfanakis K. Evaluation of standardized and study-specific diffusion tensor imaging templates of the adult human brain: Template characteristics, spatial normalization accuracy, and detection of small inter-group FA differences. NeuroImage (2018) 172:40–50. doi:10.1016/j.neuroimage.2018.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Horn A. MNI T1 6thGen NLIN to MNI 2009b NLIN ANTs transform (2016). figshare. Dataset. doi:10.6084/m9.figshare.3502238.v1

CrossRef Full Text | Google Scholar

61. Horn A, Li N, Dembek TA, Kappel A, Boulay C, Ewert S, et al. Lead-dbs v2: Towards a comprehensive pipeline for deep brain stimulation imaging. NeuroImage (2019) 184:293–316. doi:10.1016/j.neuroimage.2018.08.068

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Benabid AL, Chabardes S, Mitrofanis J, Pollak P. Deep brain stimulation of the subthalamic nucleus for the treatment of Parkinson’s disease. Lancet Neurol (2009) 8:67–81. doi:10.1016/S1474-4422(08)70291-6

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Tommasi G, Krack P, Fraix V, Le Bas J, Chabardes S, Benabid A, et al. Pyramidal tract side effects induced by deep brain stimulation of the subthalamic nucleus. J Neurol Neurosurg Psychiatry (2008) 79:813–9. doi:10.1136/jnnp.2007.117507

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Xu W, Miocinovic S, Zhang J, Baker KB, McIntyre CC, Vitek JL. Dissociation of motor symptoms during deep brain stimulation of the subthalamic nucleus in the region of the internal capsule. Exp Neurol (2011) 228:294–7. doi:10.1016/j.expneurol.2010.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Horn A, Kühn AA, Merkl A, Shih L, Alterman R, Fox M. Probabilistic conversion of neurosurgical dbs electrode coordinates into mni space. NeuroImage (2017) 150:395–404. doi:10.1016/j.neuroimage.2017.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol (1952) 117:500–44. doi:10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Popovych O, Tass P. Multisite delayed feedback for electrical brain stimulation. Front Physiol (2018) 9:46. doi:10.3389/fphys.2018.00046

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Bevan M, Wilson C. Mechanisms underlying spontaneous oscillation and rhythmic firing in rat subthalamic neurons. J Neurosci (1999) 19:7617–28. doi:10.1523/jneurosci.19-17-07617.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Zucker RS, Regehr WG. Short-term synaptic plasticity. Annu Rev Physiol (2002) 64:355–405. doi:10.1146/annurev.physiol.64.092501.114547

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Farokhniaee A, McIntyre CC. Theoretical principles of deep brain stimulation induced synaptic suppression. Brain Stimulation (2019) 12:1402–9. doi:10.1016/j.brs.2019.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Laing C, Chow C. A spiking neuron model for binocular rivalry. J Comput Neurosci (2002) 12:39–53. doi:10.1023/a:1014942129705

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Ermentrout B, Terman D. Neural networks as spatio-temporal pattern-forming systems. New York: Springer (2012).

Google Scholar

73. Compte A, Brunel N, Goldman-Rakic P, Wang X-J. Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cereb Cortex (2000) 10:910–23. doi:10.1093/cercor/10.9.910

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Terman D. An introduction to dynamical systems and neuronal dynamics. Lecture Notes Maths (2005) 1860:21–68. doi:10.1007/978-3-540-31544-5_2

CrossRef Full Text | Google Scholar

75. Benita JM, Guillamon A, Deco G, Sanchez-Vives M. Synaptic depression and slow oscillatory activity in a biophysical network model of the cerebral cortex. Front Comput Neurosci (2012) 6:64. doi:10.3389/fncom.2012.00064

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Pospischil M, Toledo-Rodriguez M, Monier C, Piwkowska Z, Bal T, Frégnac Y, et al. Minimal hodgkin-huxley type models for different classes of cortical and thalamic neurons. Biol Cybern (2008) 99:427–41. doi:10.1007/s00422-008-0263-8

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Papadopoulos L, Lynn CW, Battaglia D, Bassett DS. Relations between large-scale brain connectivity and effects of regional stimulation depend on collective dynamical state. Plos Comput Biol (2020) 16:e1008144–43. doi:10.1371/journal.pcbi.1008144

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Schirner M, McIntosh AR, Jirsa V, Deco G, Ritter P. Inferring multi-scale neural mechanisms with brain network modelling. eLife (2018) 7:e28927. doi:10.7554/elife.28927

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Valor A, Arista Romeu EJ, Escobedo G, Campos-Espinosa A, Romero-Bello II, Moreno-González J, et al. Study of methionine choline deficient diet-induced steatosis in mice using endogenous fluorescence spectroscopy. Molecules (2019) 24:3150. doi:10.3390/molecules24173150

CrossRef Full Text | Google Scholar

80. Bot M, Schuurman PR, Odekerken VJJ, Verhagen R, Contarino FM, De Bie RMA, et al. Deep brain stimulation for Parkinson’s disease: Defining the optimal location within the subthalamic nucleus. J Neurol Neurosurg Psychiatry (2018) 89:493–8. doi:10.1136/jnnp-2017-316907

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Akram H, Sotiropoulos SN, Jbabdi S, Georgiev D, Mahlknecht P, Hyam J, et al. Subthalamic deep brain stimulation sweet spots and hyperdirect cortical connectivity in Parkinson’s disease. NeuroImage (2017) 158:332–45. doi:10.1016/j.neuroimage.2017.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Meunier D, Lambiotte R, Bullmore E. Modular and hierarchically modular organization of brain networks. Front Neurosci (2010) 4:200. doi:10.3389/fnins.2010.00200

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Horn A, Reich M, Vorwerk J, Li N, Wenzel G, Fang Q, et al. Connectivity predicts deep brain stimulation outcome in Parkinson disease. Ann Neurol (2017) 82:67–78. doi:10.1002/ana.24974

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Neumann W, Schroll H, De Almeida Marcelino A, Horn A, Ewert S, Irmen F, et al. Functional segregation of basal ganglia pathways in Parkinson’s disease. Brain (2018) 141:2655–69. doi:10.1093/brain/awy206

PubMed Abstract | CrossRef Full Text | Google Scholar