Abstract
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
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 1mm3, meaning that if two (or more) ending (or starting) points lying within the same cube of 1mm3, 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 = (x1, y1, z1) and ending at y = (x2, y2, z2) 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 1mm3, 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., 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.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
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:The mean path length shows the ability of the network to spread information between any two nodes. A low mean shortest path length 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., di→j = ∞), is the global efficiency :This measure avoids calculating with infinity, since if there is no pathway between two nodes i.e., di→j = ∞ ⇒ 1/di→j = 0. The is comparable with inverse , and according to the Cauchy inequality for the arithmetic and harmonic mean, we obtainFor the augmented network, the mean path length was computed to be , while the inverse global efficacy resulted in . Figure 2C shows the distribution of distances between any two nodes (i.e., di→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:
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 . 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” (Bc) 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:where gjk(i) is the number of shortest paths from j to k passing from node i, and gjk 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 Bc(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 Bc of the network. Indeed, the large majority of nodes shows very low centrality. However, there are few nodes with high Bc. 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
TABLE 1
| Centrality measure | Subthalamic nucleus (STN) | Globus pallidus externa (GPe) | Globus pallidus interna (Gpi) | Thalamus (Tha) | Motor cortex areas (MC) |
|---|---|---|---|---|---|
| Betweenness centrality | (15, −13, −4),(13, −14, −7), (13, −11, −6), (13, −13, −8) center: (13.5, −12.8, −6.3) | (18, −4, −7), (18, −6, −2), (19, −5, −1), (20, −6, −3) center: (18.8, −5.3, −3.3) | (18, −9, −1), (20, −9, −1), (19, −9, 0), (18, −10, −1) center: (18.8, −9.3, −0.8) | (19, −7, 13), (9, −5, 4), (9, −4, 6) (8, −4, 5) center: (11.25, −5.0, 7.0) | (16, 2, 51), (16, 1, 49), (20, 9, 46), (18, 2, 49) center: (17.5, 3.5, 48.8) |
| Eigencentrality | (15, −14, −8), (14, −15, −8), (11, −12, −7), (15, −13, −4) center: (13.8, −13.5, −6.8) | (23, −7, 1), (20, −6, 1), (22, −4, 1), (20, −2, 0) center: (21.3, −4.8, 0.8) | (20, −9, −1), (12, −1, −6), (13, −1, −6), (13, −2, −7) center: (14.5, −3.3, −5) | (19, −8, 10) (18, −8, 10), (18, −8, 9), (19, −9, 9) center: (18.5, −8.3, 9.5) | (29, −20, 34), (27, −19, 33), (28, −22, 33), (19, −1, 50) center: (25.8, −15.5, 37.5) |
| Clustering coefficient | (12, −11, −9), (17, −14, −4), (15, −10, −5), (15, −17, −7), center: (14.8, −13, −6.3) | (17, −6, 10), (21, −8, −10), (16, −3, −7), (16, −5, -12) center: (17.5, −5.5, −9.75) | (13, −4, −5), (15, −5, −6), (11, −1, −5), (14, −3, −7) center: (13.3, −3.3, −5.8) | (7, −6, 8), (7, −6, 6), (17, −13, 10) (10, −5, 3) center: (10.3, −7.5, 6.8) | (27, −7, 66), (47, −2, 37), (36, 20, 32), (44, −4, 53) center: (38.5, 1.8, 47) |
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” (Ec), 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 xi is assigned. The number xi is set to be proportional to the sum of the weights of all nodes connected to i:where λ has to be identified. The last equation shows that the element xi 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 λ = λ1 the highest eigenvalue of the adjacency matrix A. Thenwhich 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 numbersi to each node. For example, in the case of two communities, then si = ±1. Here, we seek the best network partition in order to optimise the modularity function Q:where m = 1/2∑kij is the total number of edges in the network, and Bij = Aij−kikj/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 5
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” Bc and the “eigencentrality” Ec, 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 (Bc, Ec 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
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 mm3 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]:where C is the membrane capacity, Vi is the membrane potential of the ith neuron, xi denotes the gating variables n, h, r and x∞ is the steady state value for the gating variables. The quantity is the intracellular concentration of calcium. The exact description of the ionic currents ILEAK, IK, INa, ICa and the synaptic current Isyn for the STN and GP neurons are given in the supplementary material. The current IDBS 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: IDBS = 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 si, which is given by Laing and Chow [71]; Ermentrout and Terman [72]; Compte et al. [73]:where H(V) is a smooth approximation of the step function, i.e., , 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, bywhere Aij 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 towhere the factor Pj describes the probability of a neurotransmitter release (in the {ij} synapses), and follows the dynamics Benita et al. [75]:where tsp corresponds to the last spike-time of the presynaptic neuron and AD is the depression factor (0 < AD < 1), in our case, the value AD = 0.8 was used. The value P0 describes the steady state of P and in our case was set to 1. To simplify, when a presynaptic neuron fires at time tsp the functionality of the synapse (the release of neurotransmitters) is reduced (suppressed by a factor AD). 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.where C is the membrane capacity and Vi 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 ILEAK, IK and INa are the ionic currents, IT is the T-type calcium channel. The synaptic current Isyn has the form Isyn = IGPTH + ITHTH, where the GABAergic current IGPTH represents the inhibition of the GPi area to the thalamus, while ITHTH represents the internal excitatory or inhibitory thalamic connections. The current ISM 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]:where Vi is the membrane potential, and xi represents the gating variables for potassium and sodium current, of the ith neuron. The gating variable pi represents the activation gate of the slow, voltage-dependent potassium current IM. The current Iapp is added to tune the oscillatory behaviour of MC neurons around 20Hz. Each MC neuron has different value of Iapp which is extracted randomly from the interval [2, 3]. The synaptic activity is given from the current Isyn 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 of neurons over the population in the community; specifically, we define:The mean voltage activity 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 8
To better estimate the differences between the conditions, we generated difference curves of the spectrograms of Figure 7 based on the following subtraction pairs: |PX(f)|−|PY(f)|, where |PX|, |PY|, 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.,(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 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 , 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 . 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
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
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).
Statements
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ühnAKempfFBrückeCDoyleLMartinez-TorresIPogosyanAet alHigh-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. 10.1523/jneurosci.0282-08.2008
2.
NeumannW-JStaub-BarteltFHornASchandaJSchneiderG-HBrownPet alLong term correlation of subthalamic beta band activity with motor impairment in patients with Parkinson’s disease. Clin Neurophysiol (2017) 128:2286–91. 10.1016/j.clinph.2017.08.028
3.
KimJKimYNakajimaRShinAJeongMParkAet alInhibitory basal ganglia inputs induce excitatory motor signals in the thalamus. Neuron (2017) 95:1181–96.e8. e8. 10.1016/j.neuron.2017.08.028
4.
GalvanADevergnasAWichmannT. Alterations in neuronal activity in basal ganglia-thalamocortical circuits in the parkinsonian state. Front Neuroanat (2015) 9:5. 10.3389/fnana.2015.00005
5.
DeuschlGSchade-BrittingerCKrackPVolkmannJSchäferHBötzelKet alA randomized trial of deep-brain stimulation for Parkinson’s disease. N Engl J Med Overseas Ed (2006) 355:896–908. 10.1056/nejmoa060281
6.
VidailhetMJutrasM-FGrabliDRozeE. Deep brain stimulation for dystonia. J Neurol Neurosurg Psychiatry (2013) 84:1029–42. 10.1136/jnnp-2011-301714
7.
CrompeBAristietaALebloisAElsherbinySBoraudTMalletN. The globus pallidus orchestrates abnormal network dynamics in a model of parkinsonism. Nat Commun (2020) 11:1570. 10.1038/s41467-020-15352-3
8.
XuWRussoGHashimotoTZhangJVitekJ. Subthalamic nucleus stimulation modulates thalamic neuronal activity. J Neurosci (2008) 28:11916–24. 10.1523/jneurosci.2027-08.2008
9.
RubinJTermanD. High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J Comput Neurosci (2004) 16:211–35. 10.1023/b:jcns.0000025686.47117.67
10.
PopovychOTassP. Adaptive delivery of continuous and delayed feedback deep brain stimulation - a computational study. Sci Rep (2019) 9:10585. 10.1038/s41598-019-47036-4
11.
SoRKentAGrillW. 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. 10.1007/s10827-011-0366-4
12.
GalvanAWichmannT. Pathophysiology of parkinsonism. Clin Neurophysiol (2008) 119:1459–74. 10.1016/j.clinph.2008.03.017
13.
SpiliotisKSiettosC. 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. 10.1016/j.neucom.2011.06.018
14.
SiettosCStarkeJ. Multiscale modeling of brain dynamics: From single neurons and networks to mathematical tools. WIREs Mech Dis (2016) 8:438–58. 10.1002/wsbm.1348
15.
DecoGJirsaVRobinsonPBreakspearMFristonK. The dynamic brain: From spiking neurons to neural masses and cortical fields. Plos Comput Biol (2008) 4:e1000092. 10.1371/journal.pcbi.1000092
16.
DecoGPonce-AlvarezAMantiniDRomaniGLHagmannPCorbettaM. Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J Neurosci (2013) 33:11239–52. 10.1523/jneurosci.1091-13.2013
17.
BassettDSBullmoreET. Small-world brain networks revisited. Neuroscientist (2017) 23:499–516. 10.1177/1073858416667720
18.
BullmoreESpornsO. Complex brain networks: Graph theoretical analysis of structural and functional systems. Nat Rev Neurosci (2009) 10:186–98. 10.1038/nrn2575
19.
IliopoulosAPapasotiriouI. Functional complex networks based on operational architectonics: Application on eeg-based brain–computer interface for imagined speech. Neuroscience (2022) 484:98–118. 10.1016/j.neuroscience.2021.11.045
20.
De Santos-SierraDSendiña-NadalILeyvaIAlmendralJAnavaSAyaliAet alEmergence of small-world anatomical networks in self-organizing clustered neuronal cultures. PLoS ONE (2014) 9:e85828. 10.1371/journal.pone.0085828
21.
CrowellARyapolova-WebbEOstremJGalifianakisNShimamotoSLimDet alOscillations in sensorimotor cortex in movement disorders: An electrocorticography study. Brain (2012) 135:615–30. 10.1093/brain/awr332
22.
SpiliotisKStarkeJFranzDRichterAKöhlingR. Deep brain stimulation for movement disorder treatment: Exploring frequency-dependent efficacy in a computational network model. Biol Cybern (2021) 116:93–116. 10.1007/s00422-021-00909-2
23.
ButsonCRCooperSEHendersonJMMcIntyreCC. Patient-specific analysis of the volume of tissue activated during deep brain stimulation. NeuroImage (2007) 34:661–70. 10.1016/j.neuroimage.2006.09.034
24.
ButenkoKBahlsCSchröderMKohlingRvan RienenU. OSS-DBS: Open-source simulation platform for deep brain stimulation with a comprehensive automated modeling. Plos Comput Biol (2020) 16:e1008023–18. 10.1371/journal.pcbi.1008023
25.
DembekTARoedigerJHornARekerPOehrnCDafsariHSet alProbabilistic sweet spots predict motor outcome for deep brain stimulation in Parkinson disease. Ann Neurol (2019) 86:527–38. 10.1002/ana.25567
26.
van HarteveltTJCabralJDecoGMøllerAGreenALAzizTZet alNeural 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. 10.1371/journal.pone.0086496
27.
LiCLiQVan MieghemPStanleyHEWangH. Correlation between centrality metrics and their application to the opinion model. Eur Phys J B (2015) 88:65–13. 10.1140/epjb/e2015-50671-y
28.
TermanDRubinJYewAWilsonC. Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J Neurosci (2002) 22:2963–76. 10.1523/jneurosci.22-07-02963.2002
29.
NewmanMEJ. Modularity and community structure in networks. Proc Natl Acad Sci U S A (2006) 103:8577–82. 10.1073/pnas.0601602103
30.
MandaliAChakravarthyVSRajanRSarmaSKishoreA. Electrode position and current amplitude modulate impulsivity after subthalamic stimulation in Parkinsons disease—A computational study. Front Physiol (2016) 7:585. 10.3389/fphys.2016.00585
31.
FlemingJDunnELoweryM. Simulation of closed-loop deep brain stimulation control schemes for suppression of pathological beta oscillations in Parkinson’s disease. Front Neurosci (2020) 14:166. 10.3389/fnins.2020.00166
32.
EwertSPlettigPLiNChakravartyMMCollinsDLHerringtonTMet alToward defining deep brain stimulation targets in MNI space: A subcortical atlas based on multimodal MRI, histology and structural connectivity. NeuroImage (2018) 170:271–82. 10.1016/j.neuroimage.2017.05.015
33.
ChakravartyMMBertrandGHodgeCPSadikotAFCollinsDL. The creation of a brain atlas for image guided neurosurgery using serial histological data. NeuroImage (2006) 30:359–76. 10.1016/j.neuroimage.2005.09.041
34.
TianYMarguliesDBreakspearMZaleskyA. Topographic organization of the human subcortex unveiled with functional connectivity gradients. Nat Neurosci (2020) 23:1421–32. 10.1038/s41593-020-00711-6
35.
FanLLiHZhuoJZhangYWangJChenLet alThe human brainnetome atlas: A new brain atlas based on connectional architecture. Cereb Cortex (2016) 26:3508–26. 10.1093/cercor/bhw157
36.
PetersenMVMlakarJHaberSParentMSmithYStrickPet alHolographic reconstruction of axonal pathways in the human brain. Neuron (2019) 104:1056–64.e3. e3. 10.1016/j.neuron.2019.09.030
37.
MarekKJenningsDLaschSSiderowfATannerCSimuniTet alThe Parkinson progression marker initiative (ppmi). Prog Neurobiol (2011) 95:629–35. 10.1016/j.pneurobio.2011.09.005
38.
MilardiDQuartaroneABramantiAAnastasiGBertinoSBasileGAet alThe cortico-basal ganglia-cerebellar network: Past, present and future perspectives. Front Syst Neurosci (2019) 13:61. 10.3389/fnsys.2019.00061
39.
Bosch-BoujuCHylandBParr-BrownlieL. Motor thalamus integration of cortical, cerebellar and basal ganglia information: Implications for normal and parkinsonian conditions. Front Comput Neurosci (2013) 7:163. 10.3389/fncom.2013.00163
40.
TarnaudTJosephWMartensLTangheE. 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. 10.1007/s11517-018-1799-y
41.
LanciegoJLLuquinNObesoJA. Functional neuroanatomy of the basal ganglia. Cold Spring Harbor Perspect Med (2012) 2:a009621. 10.1101/cshperspect.a009621
42.
StamCReijneveldJ. Graph theoretical analysis of complex networks in the brain. Nonlinear Biomed Phys (2007) 1:3. 10.1186/1753-4631-1-3
43.
DecoGSendenMJirsaV. 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. 10.3389/fncom.2012.00068
44.
HoneyCJSpornsOCammounLGigandetXThiranJPMeuliRet alPredicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A (2009) 106:2035–40. 10.1073/pnas.0811168106
45.
HoneyCThiviergeJ-PSpornsO. Can structure predict function in the human brain?NeuroImage (2010) 52:766–76. 10.1016/j.neuroimage.2010.01.071
46.
WattsDStrogatzS. Collective dynamics of ‘small-world’ networks. Nature (1998) 393:440–2. 10.1038/30918
47.
BassettDBullmoreE. Small-world brain networks. Neuroscientist (2006) 12:512–23. 10.1177/1073858406293182
48.
MarkN. The structure and function of complex networks. SIAM Rev (2003) 45:58. 10.1137/S003614450342480
49.
NetoffTClewleyRArnoSKeckTWhiteJ. Epilepsy in small-world networks. J Neurosci (2004) 24:8075–83. 10.1523/jneurosci.1509-04.2004
50.
BermanBSmucnyJWylieKSheltonEKronbergELeeheyMet alLevodopa modulates small-world architecture of functional brain networks in Parkinson’s disease. Mov Disord (2016) 31:1676–84. 10.1002/mds.26713
51.
SheQChenGChanR. Evaluating the small-world-ness of a sampled network: Functional connectivity of entorhinal-hippocampal circuitry. Sci Rep (2016) 6:21468. 10.1038/srep21468
52.
FangJChenHCaoZJiangYMaLMaHet alImpaired brain network architecture in newly diagnosed Parkinson’s disease based on graph theoretical analysis. Neurosci Lett (2017) 657:151–8. 10.1016/j.neulet.2017.08.002
53.
Gouty-ColomerL-AMichelFBaudeALopez-PauchetCDufourACossartRet alMouse subthalamic nucleus neurons with local axon collaterals. J Comp Neurol (2018) 526:275–84. 10.1002/cne.24334
54.
AmmariRLopezCBioulacBGarciaLHammondC. Subthalamic nucleus evokes similar long lasting glutamatergic excitations in pallidal, entopeduncular and nigral neurons in the basal ganglia slice. Neuroscience (2010) 166:808–18. 10.1016/j.neuroscience.2010.01.011
55.
van den HeuvelMPSpornsO. Rich-club organization of the human connectome. J Neurosci (2011) 31:15775–86. 10.1523/JNEUROSCI.3539-11.2011
56.
NewmanM. Networks, an introductions. New York: Oxford University Press (2010).
57.
LeichtEANewmanMEJ. Community structure in directed networks. Phys Rev Lett (2008) 100:118703. 10.1103/PhysRevLett.100.118703
58.
ÅströmMDiczfalusyEMartensHWårdellK. Relationship between neural activation and electric field distribution during deep brain stimulation. IEEE Trans Biomed Eng (2015) 62:664–72. 10.1109/TBME.2014.2363494
59.
ZhangSArfanakisK. 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. 10.1016/j.neuroimage.2018.01.046
60.
HornA. MNI T1 6thGen NLIN to MNI 2009b NLIN ANTs transform (2016). figshare. Dataset. 10.6084/m9.figshare.3502238.v1
61.
HornALiNDembekTAKappelABoulayCEwertSet alLead-dbs v2: Towards a comprehensive pipeline for deep brain stimulation imaging. NeuroImage (2019) 184:293–316. 10.1016/j.neuroimage.2018.08.068
62.
BenabidALChabardesSMitrofanisJPollakP. Deep brain stimulation of the subthalamic nucleus for the treatment of Parkinson’s disease. Lancet Neurol (2009) 8:67–81. 10.1016/S1474-4422(08)70291-6
63.
TommasiGKrackPFraixVLe BasJChabardesSBenabidAet alPyramidal tract side effects induced by deep brain stimulation of the subthalamic nucleus. J Neurol Neurosurg Psychiatry (2008) 79:813–9. 10.1136/jnnp.2007.117507
64.
XuWMiocinovicSZhangJBakerKBMcIntyreCCVitekJL. 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. 10.1016/j.expneurol.2010.08.007
65.
HornAKühnAAMerklAShihLAltermanRFoxM. Probabilistic conversion of neurosurgical dbs electrode coordinates into mni space. NeuroImage (2017) 150:395–404. 10.1016/j.neuroimage.2017.02.004
66.
HodgkinALHuxleyAF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol (1952) 117:500–44. 10.1113/jphysiol.1952.sp004764
67.
PopovychOTassP. Multisite delayed feedback for electrical brain stimulation. Front Physiol (2018) 9:46. 10.3389/fphys.2018.00046
68.
BevanMWilsonC. Mechanisms underlying spontaneous oscillation and rhythmic firing in rat subthalamic neurons. J Neurosci (1999) 19:7617–28. 10.1523/jneurosci.19-17-07617.1999
69.
ZuckerRSRegehrWG. Short-term synaptic plasticity. Annu Rev Physiol (2002) 64:355–405. 10.1146/annurev.physiol.64.092501.114547
70.
FarokhniaeeAMcIntyreCC. Theoretical principles of deep brain stimulation induced synaptic suppression. Brain Stimulation (2019) 12:1402–9. 10.1016/j.brs.2019.07.005
71.
LaingCChowC. A spiking neuron model for binocular rivalry. J Comput Neurosci (2002) 12:39–53. 10.1023/a:1014942129705
72.
ErmentroutBTermanD. Neural networks as spatio-temporal pattern-forming systems. New York: Springer (2012).
73.
CompteABrunelNGoldman-RakicPWangX-J. Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cereb Cortex (2000) 10:910–23. 10.1093/cercor/10.9.910
74.
TermanD. An introduction to dynamical systems and neuronal dynamics. Lecture Notes Maths (2005) 1860:21–68. 10.1007/978-3-540-31544-5_2
75.
BenitaJMGuillamonADecoGSanchez-VivesM. Synaptic depression and slow oscillatory activity in a biophysical network model of the cerebral cortex. Front Comput Neurosci (2012) 6:64. 10.3389/fncom.2012.00064
76.
PospischilMToledo-RodriguezMMonierCPiwkowskaZBalTFrégnacYet alMinimal hodgkin-huxley type models for different classes of cortical and thalamic neurons. Biol Cybern (2008) 99:427–41. 10.1007/s00422-008-0263-8
77.
PapadopoulosLLynnCWBattagliaDBassettDS. Relations between large-scale brain connectivity and effects of regional stimulation depend on collective dynamical state. Plos Comput Biol (2020) 16:e1008144–43. 10.1371/journal.pcbi.1008144
78.
SchirnerMMcIntoshARJirsaVDecoGRitterP. Inferring multi-scale neural mechanisms with brain network modelling. eLife (2018) 7:e28927. 10.7554/elife.28927
79.
ValorAArista RomeuEJEscobedoGCampos-EspinosaARomero-BelloIIMoreno-GonzálezJet alStudy of methionine choline deficient diet-induced steatosis in mice using endogenous fluorescence spectroscopy. Molecules (2019) 24:3150. 10.3390/molecules24173150
80.
BotMSchuurmanPROdekerkenVJJVerhagenRContarinoFMDe BieRMAet alDeep brain stimulation for Parkinson’s disease: Defining the optimal location within the subthalamic nucleus. J Neurol Neurosurg Psychiatry (2018) 89:493–8. 10.1136/jnnp-2017-316907
81.
AkramHSotiropoulosSNJbabdiSGeorgievDMahlknechtPHyamJet alSubthalamic deep brain stimulation sweet spots and hyperdirect cortical connectivity in Parkinson’s disease. NeuroImage (2017) 158:332–45. 10.1016/j.neuroimage.2017.07.012
82.
MeunierDLambiotteRBullmoreE. Modular and hierarchically modular organization of brain networks. Front Neurosci (2010) 4:200. 10.3389/fnins.2010.00200
83.
HornAReichMVorwerkJLiNWenzelGFangQet alConnectivity predicts deep brain stimulation outcome in Parkinson disease. Ann Neurol (2017) 82:67–78. 10.1002/ana.24974
84.
NeumannWSchrollHDe Almeida MarcelinoAHornAEwertSIrmenFet alFunctional segregation of basal ganglia pathways in Parkinson’s disease. Brain (2018) 141:2655–69. 10.1093/brain/awy206
Summary
Keywords
electric field volume conductor model, dynamical systems, Hodgkin Huxley neurons, complex networks, movement disorder
Citation
Spiliotis K, Butenko K, van Rienen U, Starke J and Köhling R (2022) Complex network measures reveal optimal targets for deep brain stimulation and identify clusters of collective brain dynamics. Front. Phys. 10:951724. doi: 10.3389/fphy.2022.951724
Received
24 May 2022
Accepted
29 September 2022
Published
18 October 2022
Volume
10 - 2022
Edited by
Krasimira Tsaneva-Atanasova, University of Exeter, United Kingdom
Reviewed by
Qingyun Wang, Beihang University, China
Helmut Schmidt, Institute of Computer Science, Czech Academy of Sciences Praha (Prague), Czechia
Updates
Copyright
© 2022 Spiliotis, Butenko, van Rienen, Starke and Köhling.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jens Starke jens.starke@uni-rostock.de
This article was submitted to Biophysics, a section of the journal Frontiers in Physics
Disclaimer
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.