Utilising activity patterns of a complex biophysical network model to optimise intra-striatal deep brain stimulation

Data sourcesThe structural Magnetic Resonance Imaging (MRI) data are taken from the Human Connectome Project (HCP)51. These MRI data were laid over another atlas, the MNIPD25 atlas, which is most commonly used for surgery planning52,53,54 with the help of Lead-DBS, a toolbox for atlas connectivity and segregation analyses55. The structural composition striatum, in turn, is taken from the MIDA atlas48 and translated into the MNI (Montreal Neurological Institute) space by using the segmented MRI data as a reference.Modelling intrastriatal connectivity using complex networksIn contrast to other models of the striatum, we leverage a graph network modelling the connections of neurons under the hypothesis the graph properties are linked to biological activities. Thus, four main graph indices characterise the structural connectivity of the network, i.e. degree distribution (number of connections of each neuron), connection efficacy (minimal distance between nodes), clustering coefficients (numbers of locally interconnected triplets) and the degree of betweenness-centrality, i.e. the number of neurons serving as high-density hubs. The knowledge of the topological structure of the network plays an important role in emergent neural activity and functionality42,56. Furthermore, knowledge of structural details allows us to investigate the mechanisms involved in neural functionality or dysfunctionality. Network structural properties, in turn, can be identified using network statistical measures as listed above (i.e. degree distribution average path length, clustering coefficients, centrality38,49,57). In our case, the positions of nodes were defined by the data source analysis. The following section will describe how these nodes are connected to each other.Construction of complex network using the structural connectivityInitially, data source analysis defined a network consisting of 1,995 nodes. In this network, we assume that the majority of the nodes (95%) represent MSN neurons, and the remaining 5% are interneurons, consistent with58. The connectivity of the striatum is constructed following the idea of the small-world algorithm40,49,50,57,59: initially, each MSN neuron is connected with \(k=20\) neurons in the vicinity of 5 \(\textrm{mm}\) (local connections). Next, with a small probability p, and for each local connection, a new remote neighbour is added (e.g. \(p=0.05\)). For interneurons, we follow the same structure, using, however, a five times denser interneuron to MSN connectivity (i.e. \(k=100\)). In this way, the resulting network is highly clustered (like a lattice structure), additionally with a small distance between nodes (like random networks).Unlike traditionally used lattice connectivity models or all-to-all models, small-world structures better represent the physiological networks as a result of two main characteristics44,50,60,61,62,63,64: they are highly clustered, and typically show short path lengths38,49,59,65, enhancing in this way signal or rhythm propagation within the network, and the synchronisability of the network. We constructed the striatum network phenomenologically (no exact knowledge of the individual neuron connectivity in the striatum is known). Nevertheless, the choices of local connections per node k and the probability of random remote connections p in the model are successfully chosen so that the network shows higher values of clustering coefficients than a random network (where the clustering coefficients have very small values). In this sense, our choice produces more realistic connectivity features (compared, e.g., to fully connected or randomly connected models), which are closer to biological (brain) systems.The resulting striatal network is represented as graph \(G=(V,\, E)\), where V is the set of nodes and E represents the set of edges, i.e. connections. The nodes of the structural network are defined as points in three-dimensional space, as a single neuron with spatial coordinates. The connectivity can be represented with the adjacency (or connectivity) matrix A: if two neurons at positions \({\textbf {x}}=(x_1,y_1,z_1)\) and \({\textbf {y}}=(x_2,y_2,z_2)\) are linked then \(A({\textbf {x}},{\textbf {y}})=1\), otherwise \(A({\textbf {x}},{\textbf {y}})=0\). In the next section, we provide tools that allow us to extract the connectivity properties of the striatal network.Network measures characterising striatal connectivityNetwork measures59 reveal important properties of the structure, e.g. the expected number of connections for each node of the network, the centrality of nodes (importance of nodes), the average shortest path (expected number of steps between any two nodes), and, importantly, the detection of communities (or modules) and how these organise the network by creating dynamical patterns42,56. A network measure can express a specific property of a node (e.g. the clustering property of one node) and appears as an individual property of the ith node, while the resulting distribution over the whole network defines the global-macroscopic description (e.g. the statistical distribution of clustering).Degree distributionThe degree of a node i refers to the number of links connected to it59. In directed networks, a node has both an in-degree and out-degree, numbering in-coming and out-going edges, respectively. A high degree of connectivity (increased numbers of links) of the ith node defines the importance of a node in the network. The degree distribution P(k) defines the probability of a randomly selected node having a specific degree k.Path lengths, efficacy and clustering coefficientThe minimum number of steps between two nodes in the network (in the case of binary networks such as the one described here) defines the shortest path length \(d_{i\rightarrow {j}}\) between node i and j. Averaging over the set of all shortest paths, we obtain the mean path length of the network59:$$\begin{aligned} \overline{m}=\frac{\sum _{i,j} d_{i\rightarrow {j}}}{N(N-1)}. \end{aligned}$$
(1)
The mean path length shows the ability of the network to spread information (signal activity) between any two nodes. A low mean shortest path length \(\overline{m}\) highlights that any two randomly chosen nodes can interchange information very fast.If no pathway exists between node i and j then \(d_{i\rightarrow j}=\infty\), and this pair of nodes is excluded from Eq. (1). A similar measure reflecting pathway length, which, however, avoids divisions by 0, is the global efficacy \(\overline{l}\)59:$$\begin{aligned} \overline{l}=\frac{\sum _{i,j} \frac{1}{d_{i\rightarrow {j}}}}{N(N-1)}. \end{aligned}$$
(2)
now \(d_{i\rightarrow j}=\infty \implies 1/d_{i\rightarrow j}=0\). The inverse global efficacy \(1/\overline{l}\) is thus reflecting the mean shortest path.Another measure that characterises the local connectivity is the clustering coefficient. It measures the proportion of triangle loops that exist in a node, expressing a feedback mechanism that enhances the rhythm generation. Specifically, the clustering coefficient of a node i is defined as ratio:$$\begin{aligned} c(i) =\frac{\sum _{jk} a_{ij}a_{jk}a_{ki}}{k_i(k_i-1)}. \end{aligned}$$
(3)
The higher the number of triangles (that exist) with respect to the ith node, the higher the clustering coefficient.Betweenness centralityAnother   significant measure which quantifies the importance of a node is ‘betweenness centrality’. The term ‘centrality’ is related to the degree of influence which this particular node exerts within the network. ‘Betweenness centrality’ thus measures the amount of influence, which a node has with respect to the total information flow in the network. Important nodes that connect different subgraphs in the network (i.e., act as a bridge) show high betweenness centrality. The betweenness centrality BC(i) of the ith node is mathematically defined as the fraction of all shortest paths in the network that pass through the node, that is,$$\begin{aligned} BC(i)=\sum _{j\ne i \ne k}(g_{jk} (i))/g_{jk}\,, \end{aligned}$$
(4)
where \(g_{jk} (i)\) is the number of shortest paths from j to k passing over i, and \(g_{jk}\) is the number of shortest paths between nodes j and k. Bridging nodes that connect different subsets of the network often have a high betweenness centrality. Higher values of BC(i) indicate that the node acts as a central hub. The importance of these hubs is also highlighted pathophysiologically as such hubs are ideally suited as targets of therapeutic intervention, i.e. for DBS38.Detection of communities and modularityNetworks characteristically contain subsets (subgraphs) that have dense internal connectivity (connectivity among nodes in the subset) and sparse connections to other subgraphs65. The partition of the network into densely connected subgraphs (or communities) plays a significant role in information processing within the network, and it is also related to different biological functions of the area (e.g. striatum)66. Assigning and allocating these densely connected communities to brain structures allows the construction of a modular view of the dynamics of the network67,68.The modality index identifies such densely connected communities. The modality index65 assigns a community number \(s_i\) to each node. For example, if there are two communities, then \(s_i=\pm 1\). Here, we seek the best network partition in order to maximise the modularity function Q:$$\begin{aligned} Q=\frac{1}{4m}{s}^TBs \end{aligned}$$
(5)
where \(m=1/2 \sum k_{ij}\) is the total number of edges in the network, and \(B_{ij}=A_{ij}-k_i k_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 B65,69.Modelling and simulation of MSN networksModelling MSN neuronsThe dynamics of each MSN neuron are modelled by current balance equations for the membrane potential:46$$\begin{aligned} C\frac{dV_i}{dt}&=-I_{\text {LEAK}}-I_{\text {K}}-I_{\text {Na}}-I_{\text {M}}-I_{\text {syn}}+I_{\text {app}}\, , \end{aligned}$$
(6)
$$\begin{aligned} \frac{dx_i}{dt}&=(x_{\infty }-x_i)/\tau _{x_i}\, , \end{aligned}$$
(7)
where C is the membrane capacity, and \(V_i\) is the membrane potential of the ith neuron. The current balance Eq. (6) contains four membrane currents46, the fast sodium and potassium currents \(I_{\text {Na}}\) and \(I_{\text {K}}\), the leak current \(I_{\text {LEAK}}\), and an M-current \(I_{\text {M}}\). All ionic currents follow the Hodgkin–Huxley formalism47: \(I=g_X m^{n_1}_X h^{n_2}_X \cdot (V-E_X)\), where the exponents \(n_1, n_2\) represent the number of activation-inactivation channels, respectively, \(g_X\) is the maximum conductance of the X ion, and \(E_X\) stands for the reversal potential for each X ion (\(X \in \{\text {Na, K}\}\)). Specifically, the sodium current has three activation gates and one inactivation gate, that is, \(I_{\text {Na}}=g_{\text {Na}} m^3_{\text {Na}} h_{\text {Na}} \cdot (V-E_{\text {Na}})\). The potassium current has the form \(I_{\text {K}}=g_{\text {K}} m^4_{\text {K}} \cdot (V-E_{\text {K}})\) (i.e., \(n_1 = 4,n_2=0\)). Finally, the M-current and leak currents follow \(I_{\text {M}}=g_{\text {M}} m_{\text {M}} \cdot (V-E_{\text {K}})\) and \(I_{\text {LEAK}}=g_{\text {LEAK}} \cdot (V-E_{\text {LEAK}})\), respectively (see also Table 1).The variable \(x_i\) denotes the gating variables \(m_{\text {X}}\) and \(h_{\text {X}}\). Following the Hodgkin–Huxley formalism46,47, the function \(x_{\infty }\) is given by \(x_{\infty }=\frac{\alpha _l}{\alpha _l+\beta _l}\) and the time by \(\tau _{X}=\frac{1}{\alpha _l+\beta _l}\), \(l \in m_{\text {X}}, h_{\text {X}}\). For the sodium current and for the gating variable \(m_{\text {X}}\), we obtain$$\begin{aligned} a_m(V)=0.32\frac{V+54}{1-e^{-(V+54)/4)}}, \hspace{0.5cm} b_m(V)=0.28\frac{(V+27)}{e^{(V+27)/5}-1}; \end{aligned}$$
(8)
similarly, for the gating variable \(h_{\text {X}}\):$$\begin{aligned} a_h(V)=0.128e^{-(V+50)/18}, \hspace{0.5cm} b_h(V)=\frac{4}{1+e^{-(V+27)/5}}. \end{aligned}$$
(9)
For the potassium current, with only one activation gating: \(m_{\text {k}}\)$$\begin{aligned} a_m(V)=0.032\frac{V+52}{1-e^{-(V+52)/5)}}, \hspace{0.5cm} b_m(V)=0.5 {e^{-(V+57)/40}}. \end{aligned}$$
(10)
For the M-current, we obtain$$\begin{aligned} a_m(V)=0.032\frac{V+52}{1-e^{-(V+52)/5)}}, \hspace{0.5cm} b_m(V)=0.5 {e^{-(V+57)/40}}. \end{aligned}$$
(11)
The current \(I_{\text {app}}\) is written as \(I_{\text {app}}=I_{0}+I_{\text {DBS}}\) in Eq. (7), where \(I_0\) represents a network activation current, describing the dependence of the neuronal activation due to dopamine receptor activation, or due to intensity of cortical-striatal connectivity. The current \(I_{\text {DBS}}\) specifically models deep brain stimulation, and it is applied on any element within the reach of the stimulation electrodes (the stimulation is modelled as declining exponentially, see Eq.(12)). The mathematical description is given by the form:$$\begin{aligned} I_{\text {DBS}}=A_{\text {DBS}}e^{-\frac{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2}{\sigma ^2}}H(\sin (2\pi t/T_{\text {DBS}})\cdot (1-H(\sin (2\pi (t+\delta _{\text {DBS}})/T_{\text {DBS}})), \end{aligned}$$
(12)
while in the absence of DBS treatment \(I_{\text {DBS}}=0\).The synaptic currents \(I_{\text {syn}}\) for MSN neurons can be written as a sum: \(I_{\text {syn}}=I_{\text {MSN} \rightarrow \text {MSN}}+I_{\text {FS} \rightarrow \text {MSN}}\), where \(I_{\text {MSN} \rightarrow \text {MSN}}\) is the inhibitory synaptic current between MSN neurons, carried by fast-spiking neurons, whose exact description is given in the next section.Modelling Fast Spiking (FS) neuronsThe dynamics of each FS neuron are modelled by current balance equations for the membrane potential:46:$$\begin{aligned} C\frac{dV_i}{dt}&=-I_{\text {LEAK}}-I_{\text {K}}-I_{\text {Na}}-I_{\text {D}}-I_{\text {syn}}+I_{\text {app}}. \end{aligned}$$
(13)
$$\begin{aligned} \frac{dx_i}{dt}&=(x_{\infty }-x_i)/\tau _{x_i} , \end{aligned}$$
(14)
where C is the membrane capacity, and \(V_i\) is the membrane potential of the ith FS neuron. The current balance Eq. (14) contains four membrane currents46:\(I_{\text {Na}}=g_{\text {Na}} m^3_{\text {Na},\infty } h_{\text {Na}} (V-E_{\text {Na}})\), \(I_{\text {K}}=g_{\text {K}} m^2_{\text {K}} (V-E_{\text {K}})\), \(I_{\text {LEAK}}=g_{\text {LEAK}} (V-E_{\text {LEAK}})\), while the fast-activating, slowly inactivating dendritic potassium D-current has the form46,70: \(I_{\text {D}}=g_{\text {D}} m^3_{\text {D}} h_{\text {D}} (V-E_{\text {K}})\), with three activation gates and one inactivation gate (i.e. \(n_1=3, n_2=1\), see also Table 1), thus imposing a delay in firing upon depolarisation. For the sodium current \(I_{\text {Na}}\), and for the gating variable m, we obtain:$$\begin{aligned} m_{\infty }=\frac{1}{1+e^{-(V+24)/11.5}}, \end{aligned}$$
(15)
and for the sodium inactivation$$\begin{aligned} h_{\infty }=\frac{1}{1+e^{-(V+58.3)/6.7}}. \end{aligned}$$
(16)
For the potassium \(I_{\text {K}}\)current, and for the activation m variable, we use the equation:$$\begin{aligned} m_{\infty }=\frac{1}{1+e^{-(V+12.4)/6.8}}. \end{aligned}$$
(17)
Finally, for the D-current with activation and inactivation variables, we use$$\begin{aligned} m_{\infty }=\frac{1}{1+e^{-(V+50)/20}}, \end{aligned}$$
(18)
and for the inactivation$$\begin{aligned} h_{\infty }=\frac{1}{1+e^{-(V+70)/6}}. \end{aligned}$$
(19)
Table 1 The currents for medium spiny neurons (MSN) and fast-spiking neurons (FS).Description of the network inhibitory synaptic activityThe coupling between the neurons in eqns. (6) and (13) is described by the synaptic current \(I_{\text {syn}}\). Initially, we model the activation of a synapse using the activation variable \(s_i\) (for the ith neuron), which is given by71,72,73:$$\begin{aligned} \frac{ds_i}{dt} = \alpha (1-s_i)H(V_i)-\beta s_i, \end{aligned}$$
(20)
where the function H(V) is a smooth approximation of the step (Heaviside) function \(H_\text {step}\) (i.e. \(H_\text {step}(x)=1, x>0\) and \(H_\text {step}(x)=0, x<0\).) The variable \(s_i\) describes the activation of synapses from the pre-synaptic neuron i to the post-synaptic neuron j. The form of function H is given by:$$\begin{aligned} H(V)=1+\tanh (V/10), \end{aligned}$$
(21)
The parameters \(\alpha , \beta\) in Eq. (20) are related to the activation and inactivation time scales, respectively, of the inhibitory (GABA-ergic) synaptic connections. In cases of MSN – MSN and MSN – FSI interactions, the activation rates in equation (20) are \(\alpha =4, \beta =1/13 \approx 0.08\). Similarly, for FS – MSN and FS – FS interactions, the activation rates in equation (20) are \(\alpha =4, \beta =1/11\).For each ith neuron in the network, the total synaptic inhibition which it receives from the pre-synaptic neurons is:$$\begin{aligned} I_{i,\text {GABA}}=g_{\text {X} \text {Y}}(V_i-E_{\text {GABA}})\sum _{j}{A_{ij}s_{j}}, \end{aligned}$$
(22)
with \(E_{\text {GABA}}=-80mV\). The matrix element \(A_{ij}\) has the value 1 or 0, depending on whether neurons i and j are connected or not. In this way, it resembles the modified Watts and Strogatz (WS) small-world topology49,57,59,60,63,74,75,76. The summation is taken over all presynaptic neurons. The parameter \(g_{\text {X}\text {Y}}\) represents the conductance between X and Y interactions \(X, Y \in \{\text {MSN}, \text {FS}\}\).Modelling the connectivity within the striatumThe synaptic current for the MSN neurons is given by:$$\begin{aligned} I_{\text {syn}}=I_{\text {MSN} \rightarrow \text {MSN}}+I_{\text {FS} \rightarrow \text {MSN}}. \end{aligned}$$
(23)
The current \(I_{\text {MSN} \rightarrow \text {MSN}}\) indicates the inhibition between MSN-MSN neurons, while the second term \(I_{\text {FS} \rightarrow \text {MSN}}\) represents the interneuronal inhibition. Taken together, the mathematical form of the synaptic current for MSN neurons is:$$\begin{aligned} I_{i,\text {MM}}= g_{\text {MM}}(V_i-E_{\text {GABA}})\sum _{j}{A_{ij}s_{j}}+g_{\text {FM}}(V_i-E_{\text {GABA}})\sum _{j}{A_{ij}s_{j}}, \end{aligned}$$
(24)
where, again, the element \(A_{ij}\) has the value 1 or 0, depending on whether neurons i and j are connected or not, while the sum is taken over all presynaptic neurons.Similarly, for the FS neurons, the synaptic current is analysed as a sum:$$\begin{aligned} I_{\text {syn}}=I_{\text {FS} \rightarrow \text {FS}}+I_{\text {MSN} \rightarrow \text {FS}}. \end{aligned}$$
(25)
The current \(I_{\text {FS} \rightarrow \text {FS}}\) represents the rare case of FS-FS inhibition, while the second term \(I_{\text {MSN} \rightarrow \text {FS}}\) imitates the feedback inhibitory loop of MSN to interneurons. Then, the mathematical form of the synaptic current for each FS neuron is:$$\begin{aligned} I_{i,\text {FS}}= g_{\text {FF}}(V_i-E_{\text {GABA}})\sum _{j}{A_{ij}s_{j}}+g_{\text {MF}}(V_i-E_{\text {GABA}})\sum _{j}{A_{ij}s_{j}}. \end{aligned}$$
(26)
Restoring normal striatal activity by optimising DBS positionAs already emphasised, striatal neuronal activity is not only involved in major tasks such as movement control but also in decision-making, reward behaviour and other cognitive/emotional tasks, with behavioural control being driven by the ventral parts of the striatum1,2. Under pathological conditions (e.g. obsessive-compulsive disorder), abnormal striatal activity has been reported; decreased dopamine levels in conjunction with aberrant cortico-striatal interactions are thought to lead to reduced striatal network activity7,27,46.In order to quantify striatal activity, we define the mean network activity as a macroscopic variable77:$$\begin{aligned} \langle a\rangle _T(t)=\frac{1}{T} \int _{t}^{t+T}{n_A([t, t+T])} dt, \end{aligned}$$
(27)
where N is the number of neurons in the population, \(n_A([t, t+T])\) is the number of spikes (summed over all neurons in the population) that occur between t and \(t+T\), and where T is a small macroscopic time interval (10 ms). The values obtained were then divided by the number of bins (100); thus, \(\langle a\rangle _T\) actually stands for the number of spikes per 0.1 ms bin within the entire population of neurons (2000 neurons). In order to define the mean activity as spikes per neuron in Hz, we ultimately multiply the values by 5 in two steps. In the initial step, the value is multiplied by 10,000 (0.1 ms * 10,000) to obtain a value per second. Subsequently, the value is divided by 2,000 (the number of neurons) in order to obtain the value per neuron (in Hz). Thus, we define the average network activation rate in Hz as \(\langle a\rangle _T\). Low values of the macroscopic network rate indicate low striatal output, characteristic of a disturbed and abnormally low dopamine, or low intensity of corticostriatal activity.Another macroscopic variable that we explore is the mean membrane voltage \(\overline{V}\) of neurons in the network; specifically, we define:$$\begin{aligned} \overline{V}_x(t)= \frac{1}{N}\sum _{k=1}^{N}{V_k(t)}. \end{aligned}$$
(28)
The mean voltage activity \(\overline{V}\) (indirectly related to the local field potential (LFP) in the case of supra-threshold values resulting in spiking) is utilised for the characterisation of synchronised rhythm (through Fourier spectrum) under different states (healthy, abnormal or abnormal plus DBS).Optimising DBS parameters using macroscopic quantities of the striatal networkDifferences in striatal targeting areas, as well as different intensity and frequency values of the DBS signal, result in differences in distant network activation. In the model, we vary position, stimulation intensity and frequency, resulting in the parameter vector \(r=(x_0,y_0,z_0, A_{\text {DBS}},f_{\text {DBS}}) \in \mathbb {R}^{5}\) to estimate optimal DBS outcome; see Eq. (12). The effectiveness of DBS is evaluated using three objective functions. The objective functions defined below indicate the impact of DBS, i.e. the ability of DBS to restore neuronal activity to a healthy state.Optimising with respect to mean network activityThe first objective function is based on the mean network activity \(\langle a \rangle _T\) Eq. (27). We define the objective function as the difference between healthy and DBS mean activities, i.e.$$\begin{aligned} \Phi _1(r):=|| \langle a\rangle _{T, X_\text {healthy}}- \ \langle a\rangle _{T,X_\text {DBS}} ||_{L^1} =\int _{a}^{b}{\big |\langle a\rangle _{T, X_\text {healthy}}(t)- \ \langle a\rangle _{T,X_\text {DBS}}(t)\big |}dt, \end{aligned}$$
(29)
where a, b are the times of activation and inactivation of DBS (usually \(a=0\) and \(b=1\) sec), \(r=(x_0,y_0,z_0,A_{\text {DBS}},f_{\text {DBS}}) \in \mathbb {R}^{5}\) is the DBS parameter vector, i.e. the position \((x_0,y_0,z_0)\) of the DBS electrode, the amplitude \(A_{\text {DBS}}\) and the frequency f of the pulse. The values of the model parameters, \(r \in \mathbb {R}^{5}\) were estimated numerically by minimising the residual, i.e.:$$\begin{aligned} \mathop {\mathrm {arg\,min}}\limits _{r \in \mathbb {R}^{5}} \Phi _1(r):= \left\{ r: \min \Phi _1(r),\ r \in \mathbb {R}^{5} \right\} . \end{aligned}$$
(30)
Minima of the difference function \(\Phi _1\), (similarly  for \(\Phi _2\) and \(\Phi _3\)) , i.e. values of the objective function close to 0 imply \(\langle a \rangle_{T, X_\text {healthy}} \approx \langle a \rangle_{T, X_\text {DBS}}\), which it is interpreted as effective DBS action restoration of normal striatal activity.Optimising with respect to network rhythmicityThe second objective function \(\Phi _2\) is based on the frequency spectrum produced by the model. For this, we define the Fourier transform of the mean activity \(\overline{V}\) i.e.$$\begin{aligned} X_r(f)= \int _{-\infty }^{\infty }{\overline{V}(t)e^{-j f t}dt}, \end{aligned}$$
(31)
where j is the imaginary unit. The power spectrum is defined as \(|P(f)|=|X(f)|^2\). To estimate the effectiveness of DBS, we calculate differences curves 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. X, Y in this case represent the conditions under DBS X, and healthy state Y.The objective function \(\Phi _2\) is then defined as the area under the curve (AUC,78) for each of these difference pairs, i.e.$$\begin{aligned} \Phi _2(r):= \int _{a}^{b} \Big | \ |P_X(f)|-|P_Y(f)|\ \Big | df, \end{aligned}$$
(32)
where a, b are the frequency range (i.e. [a,b]=[0 300]Hz) and \(r=(x_0,y_0,z_0,A_{\text {DBS}},f_{\text {DBS}}) \in \mathbb {R}^{5}\) is the DBS parameter vector. Similar to the previous case, i.e. \(\Phi _1\), the values of the model parameter, \(r \in \mathbb {R}^{5}\) were estimated numerically by minimising the residual, that is,$$\begin{aligned} \mathop {\mathrm {arg\,min}}\limits _{r \in \mathbb {R}^{3}} \Phi _2(r):= \left\{ r: \min \Phi _2(r),\ r \in \mathbb {R}^{5} \right\} . \end{aligned}$$
(33)
We can also combine the aforementioned objective functions to obtain a third optimisation scheme that takes into account both macroscopic characteristics (i.e. rate and rhythmicity) of the network activity. We will present this scheme in the next subsection.Optimising using a combination of network rhythmicity and firing rateWe define a combination of objectives functions eqs. (32) and (29):$$\begin{aligned} \Phi _3(r):= \Phi _1(r)+\gamma \cdot \Phi _2(r), \end{aligned}$$
(34)
where \(\gamma\) is a scaling factor that balances (reduces or increases) the importance of the phase spectrum in the optimisation process. In this instance, \(\gamma\) is adapted by iterative optimisation steps with decreasing sizes of \(\gamma\). At the chosen value, the objective function \(\Phi _3\) combines both \(\Phi _1, \Phi _2\) optimally, i.e. considers both the rate and the rhythmicity effect.The minimisation problem was solved in all cases using a nonlinear least-squares solver (the MATLAB function lsqnonln). The step size tolerance was set to \(\texttt {tol}(X)=0.001\) and the function tolerance was set to \(\texttt {tolF}=0.1\). The maximum number of iterations was set to 100.Figure 2Connectivity properties of the striatal network. Blue histograms stand for the MSN neurons, while red for FS neurons. (A, B) Distribution of the number of connections leaving a node (out-degree distribution). Clearly, MSN neurons show a relatively symmetrical spread of data with relatively sparse connections (mean 25 connections), unlike fast-spiking GABAergic interneurons (minority; mean 110 connections). (C, D) Distribution of efficacy (the number of steps as maximal distance between any two neurons, here given as reciprocal value to avoid divisions by 0). (E, F) The distribution of clustering coefficients, i.e. the number of locally interconnected neuronal triplets. MSN neurons show a relatively high number of triangle loops, see Eq. (3), while FS neurons show a lower number of triangles. (G, H) The distribution of betweenness centrality (BC), measuring the number of paths passing from a given node. MSN neurons have a lower value of betweenness centrality, while FS neurons show almost 10 times higher values of BC.Figure 3Modularity: Community detection in the striatal network. Communities were extracted using an iterative process of branching groups of neurons fulfilling two criteria: (a) dense connectivity among members and (b) sparse connectivity to the other communities. This iterative branching stops when an optimum is reached in any branch. In this way, the striatal network is partitioned into separate subgraphs, using a commonly used modularity-index algorithm69. Using a boundary condition that at least 180 neurons should be within one community (9 % of the population), in the present model, the algorithm identifies six communities, which remain stable with repetitive (20 times) realisations. As the figure shows, three communities are located in the caudate nucleus (top, red, yellow and violet hues) and three in the putamen (bottom, blue and green hues). The ten nodes with the highest betweenness-centrality (hubs, black circles) are equally distributed between the caudate nucleus (n = 5) and the putamen (n = 5).Sensitivity analysisThe last part of our theoretical analysis (Sensitivity Analysis) is further investigating the model parameters. Specifically, the aim is to examine the significance or the sensitivity of the emergent network behaviour with respect to vital parameters such as the intensity of cortico-striatal connectivity and the connectivity conductance between MSN neurons. Using Sobol’ indices as order parameters, we are partitioning the variance of the output into fractions according to the parameter’s input contribution. The Sobol’ indices lie in the interval [0, 1]. A Sobol’ index close to 1 reflects a more significant influence of the parameter on the model’s response.For the global sensitivity analysis, i.e., to quantify the effects of the input random variables in the variance of the response of the model, we use Sobol’ indices79, based on functional decomposition applied to the variance. The implementation and use of the method are straightforward. In the system of interest, we consider \(\textbf{X}\) as a random input vector following a certain probability density function \(f_\textbf{X}\), and \(\textbf{Y}\) as the response of this system.The total variance of a model’s response is denoted by \(\textrm{var}(\textbf{Y})\), and the conditional variances, which consider the contribution of one or more parameters, are denoted by \(\textrm{var}[\mathbb {E}(\textbf{Y}|\textbf{X}_k)]\), where \(\textbf{X}_k\) is the input vector with k the parameters, with \(k\in \mathbb {N}\). From Sobol’ decomposition and its orthogonality, we obtain the total variance of \(\textbf{Y}\) as the sum of the conditional variances, i.e.,$$\begin{aligned} \textrm{var}(\textbf{Y}) = \sum \limits _{k} \textrm{var}[\mathbb {E}(\textbf{Y}|\textbf{X}_k)], \end{aligned}$$
(35)
where \(\mathbb {E}\) is the expectation. In this sense, Sobol’ indices are defined by$$\begin{aligned} S_k = \frac{\textrm{var}[\mathbb {E}(\textbf{Y}|\textbf{X}_k)]}{\textrm{var}(\textbf{Y})}, \end{aligned}$$
(36)
such that their total sum equals one. Therefore, a Sobol’ index can have a value within the interval [0,1]. The closer the sensitivity index of a parameter approaches the value 1, the greater its influence on the response of the model.The aim is to examine the significance of a particular variable by measuring the proportion of variance in the Quantity of Interest (QoI) for which it is responsible. For this, we compute the first-order Sobol’ index, which quantifies the share of variance in the output due to the examined parameter. In addition, the higher-order index quantification considers the interaction of all studied variables.To compute Sobol’ indices, we choose the Polynomial Chaos Expansion (PCE), which has proven to be very efficient in dealing with uncertainties. The PCE method is an efficient alternative to the Monte-Carlo methods, being much faster in obtaining similar results, as long as the number of uncertain variables is less than 20, see, e.g.,80. PCE can significantly decrease the number of simulations while providing an accurate approximation of the model’s response.The PCE method uses an approximation of a system’s random model response \(\textbf{Y}\), assuming that \(\textbf{Y}\) has finite variance. The representation of \(\textbf{Y}\) can be given by$$\begin{aligned} \textbf{Y} = \sum \limits _{n=0}^{N_p-1}c_n\Psi _n(\textbf{X}), \end{aligned}$$
(37)
with \(N_p\) as the number of expansion factors, \(c_n\) as the coefficients, and \(\Psi _n\) as a multidimensional generalised PC basis defined in a Hilbert space \(L^2(\textbf{X},f_\textbf{X})\). The coefficients are computed following$$\begin{aligned} c_n = \frac{\langle \textbf{Y}, \Psi _n(\textbf{X})\rangle }{\langle \Psi _n,\Psi _n\rangle }, \quad \forall \, n=0,1,…, N_p-1. \end{aligned}$$
(38)
The equation above is obtained by taking the inner product of Eq. (37) and \(\Psi _n\) and using the orthogonality of the basis. The coefficients are calculated with a pseudo-spectral projection method.In particular, since we are using the uniform distribution \(f_\textbf{X}\) of the parameters, \(\Psi _n\) is based on Legendre polynomials, which are orthogonal with respect to the uniform distribution. These polynomials are obtained with a three-term recurrence formula,$$\begin{aligned} P_0(x)= & {} 1, \end{aligned}$$
(39)
$$\begin{aligned} P_1(x)= & {} x,\end{aligned}$$
(40)
$$\begin{aligned} (n+1)P_{n+1}(x)= & {} (2n+1)xP_n(x)-nP_{n-1}(x), \end{aligned}$$
(41)
where \(P_n(x)\), \(n\in \mathbb {N}\), denotes the Legendre polynomials.Once the simulations are complete, both an approximation of the model output and the calculated statistics of the QoIs based on the PCE are available.Many statistical measures can then be obtained directly from the PCE representation, such as the mean and variance of the model response. The Sobol’ indices can also be computed based on the polynomial chaos decomposition of the model81 and are called PC-based Sobol’ indices.