Data-driven and equation-free methods for neurological disorders: analysis and control of the striatum network

The striatum as part of the basal ganglia is central to both motor, and cognitive functions. Here, we propose a large-scale biophysical network for this part of the brain, using modified Hodgkin-Huxley dynamics to model neurons, and a connectivity informed by a detailed human atlas. The model shows different spatio-temporal activity patterns corresponding to lower (presumably normal) and increased cortico-striatal activation (as found in, e.g., obsessive-compulsive disorder), depending on the intensity of the cortical inputs. By applying equation-free methods, we are able to perform a macroscopic network analysis directly from microscale simulations. We identify the mean synaptic activity as the macroscopic variable of the system, which shows similarity with local field potentials. The equation-free approach results in a numerical bifurcation and stability analysis of the macroscopic dynamics of the striatal network. The different macroscopic states can be assigned to normal/healthy and pathological conditions, as known from neurological disorders. Finally, guided by the equation-free bifurcation analysis, we propose a therapeutic close loop control scheme for the striatal network.

1 Introduction and context

Complex dynamical systems of interacting units appear in nature across several disciplines. Examples of these systems are networks of coupled neurons in the brain, epidemiological networks of interacting individuals during a virus spreading, and social or economic networks of human action and perception. A common characteristic of these networks is the existence of well-defined rules for each individual entity, the so-called microscopic description, while the emergent network behaviour evolves on a different level, the macroscopic scale.

The macroscopic description, say, in the form of ordinary or partial differential equations, governs the time evolution of few macroscopic variables, which are often given by low order statistics such as densities or correlation functions. It is however very challenging, if possible at all, to derive such a macroscopic description from a microscopic model, without making assumptions about the connectivity of the system, see, e.g., (Kevrekidis and Samaey, 2009; Montbrió et al., 2015). In neuroscience, and specifically for brain networks, the microscopic description is based on the electrochemical activity of individual cells which is frequently modelled by Hodgkin-Huxley equations (Hodgkin and Huxley, 1952; Terman et al., 2002; Spiliotis et al., 2022b). These cell-neurons interact through synaptic connections, and the mathematical description results in large systems of coupled nonlinear differential equations. The heterogeneous connectivity, the nonlinear behaviour of each cell, and the stochastic environment are factors which increase the complexity of the emergent network behaviour. Existence of multiple stationary states, sustained oscillations (Deco et al., 2008; Spiliotis and Siettos, 2011; Deco et al., 2013), as well as travelling waves and spatio-temporal chaos (Laing and Chow, 2002; Bhattacharya et al., 2022; Palkar et al., 2023), are signatures of the rich nonlinear behaviour of neural networks at the macroscopic level (Deco et al., 2008; Spiliotis and Siettos, 2011; Crowell et al., 2012; Deco et al., 2013; de Santos-Sierra et al., 2014; Siettos and Starke, 2016; Spiliotis et al., 2022b).

In previous studies (Spiliotis et al., 2022a; Spiliotis et al., 2022b; Spiliotis et al., 2024) we developed a large-scale computational model of the basal ganglia network and thalamus to describe movement disorders and treatment effects of deep brain stimulation. The model of this complex network covers three areas of the basal ganglia region: the subthalamic nucleus, the globus pallidus, both pars externa and pars interna, and the thalamus and motor and pre-motor cortex. Macroscopic analysis of the network dynamics allowed us to study the differences in neural activation patterns that will emerge within the brain’s structural network when simulating different medical conditions. For example, our computational model suggests that spatio-temporal activity in the basal ganglia network shows travelling wave solutions with more varying structures in the normal state as compared to the Parkinsonian state, see (Spiliotis et al., 2022b). In addition, the macroscopic analysis yields optimal frequency ranges for deep brain stimulation as well as optimal positions for the electrodes (Spiliotis et al., 2022a).

In this work, we focus on the striatum, an essential intermediate area of the brain that connects cortical to deep brain regions. The striatum belongs to the basal ganglia area and orchestrates activities for controlling movement, decision-making, choosing actions, and those maximising reward and other psychological behaviours (Calabresi et al., 2007; Crittenden and Graybiel, 2011; Calabresi et al., 2014). The striatum integrates cortical signals to create motor activities based on experience and forthcoming selections. The significance of striatum functionality is also accentuated by its involvement in a vast number of neurological diseases ranging from Parkinson’s disease, Huntington’s disease, and dystonia to psychological disorders such as obsessive-compulsive disorder, depression, impulsivity, and attention-deficit hyperactivity disorder (Remijnse et al., 2006; Crittenden and Graybiel, 2011).

Our main aim is the development of a mathematical-computational framework to analyse the macroscopic network behaviour of the striatum area, using data from microscopic simulations of a modified Hodgkin-Huxley network of neurons. We achieve our goal by an equation-free approach (Gear et al., 2005; Kevrekidis and Samaey, 2009; Marschler et al., 2014a; Laing, 2018). We identify the mean synaptic activity as the appropriate macroscopic variable that captures the network dynamics. This is also justified from other computational and medical-clinical studies Popovych and Tass (2019); Buzsáki (2004); Parasuram et al. (2016) since neural network activity like synchronisation, is also reflected by the amplitude of the local field potential (LFP) which is modelled as an ensemble-averaged synaptic activity of neurons. The equation-free method allows to perform a numerical bifurcation and stability analysis for the macroscopic dynamics. Our analysis will reveal an interesting property which is not accessible by straightforward simulations of the network, namely, the existence of two macroscopic network states, a high activation state and an unstable low activation state. The different macroscopic states can be related to healthy and pathological conditions existing in neurological disorders. During obsessive-compulsive disorder there is an increased cortico-striatal activity (Maltby et al., 2005; Marsh et al., 2014). Our computational model also predicts this high activation solution. Additionally the model shows a second solution which provides a low activation state, leading the striatum activity to a less pathological activation. Such a state is a potential healthy target for deep brain stimulation and may result in strategies for an efficient treatment. In fact, based on our analysis we propose a closed loop macroscopic control scheme which provides better performance compared to a straightforward deep brain stimulation approach.

2 Construction of the striatum model

We extract the surface of the striatum using magnetic resonance medical data taken from a previously published atlas (Iacono et al., 2015) and transform into the MNI (Montreal Neurological Institute) coordinate system. We place neurons randomly inside this area, see Figure 1A.

Figure 1

Figure 1. Representation of the striatal network: (A) Schematic representation of the striatum area as obtained in a MNI coordinate space. Colour code represents the membrane electrical activity in mV. (B) Raster plot representation of the network activity in time (in ms) and space (index of neuron of the nucleus). Black dots represent activated neurons (i.e., time dependent action potentials passing through $-15$mV to positive values) (C) Time series of two representative medium spiny neurons (MSN) of the striatum. (D) Fourier analysis for the mean membrane activity $\bar{V}(t)=1/N{\sum }_i{V}_i(t)$ showing a $\gamma$ rythm, i.e., with main characteristic frequency above 30Hz (such rhythm appears, for instance, in the striatum during motivated behaviour and reward processing (Kalenscher et al., 2010)).

2.1 Modelling the striatum network by small word connectivity

We place in total 1995 neurons as network nodes in the striatum area. In line with medical studies (Yager et al., 2015) we assume that the vast majority of nodes (i.e., 95% of nodes) represent medium spiny neurons (MSN) while the remaining 5% of nodes are interneurons. The actual connectivity of the striatum is constructed following the idea of the small-world algorithm (Watts and Strogatz, 1998; Bassett and Bullmore, 2006; Stam and Reijneveld, 2007; Bullmore and Sporns, 2009; Spiliotis and Siettos, 2011). Small-world structures suitably model physiological networks (Netoff et al., 2004; Bassett and Bullmore, 2006; Bassett and Bullmore, 2017; de Santos-Sierra et al., 2014; Berman et al., 2016; She et al., 2016; Fang et al., 2017) since those networks are highly clustered and typically show short path lengths, enhancing in this way signal or rhythm propagation within the network and support synchronisation. Initially, each MSN is connected with $k=20$ neurons in the local neighbourhood of 5$\mathrm{m}\mathrm{m}$ diameter. Then, with a small probability $p$, for each local connection a new remote neighbour is added. For interneurons we follow the same approach, using, however, a five times higher interneuron to MSN connectivity with $k=100$ links. In this way, we obtain a network which is highly clustered with a small distance between nodes. The resulting striatal network gives a graph $G=(V,E)$, where $V$ is the set of nodes and $E$ the set of edges, i.e., the connections between neurons. The connectivity can be captured by a so-called adjacency matrix $A$, where $A_{ij}=1$ if there is a connection from node $j$ to node $i$, and $A_{ij}=0$ otherwise. Each of the nodes represents a neuron with dynamics being described by modified Hodgkin-Huxley equations. The position of the striatum in the model is based on a medical atlas, and the positions of neurons are constructed based on this information. That means each index $i$ comes with the Cartesian coordinates of the neuron and set of links. The connectivity of the neurons is constructed using the Watts and Strogatz small-world algorithm. In this type of connectivity, the nearest neurons are connected, and their activity is communicated to the nearest nodes, analogous to a graph Laplacian. Additionally, the small-world connectivity allows rare remote connections with a small probability, offering a more realistic neuronal network activity. As our model contains the actual geometric information of the position of neurons, we are finally able to model deep brain stimulation, where the position of electrodes and their spatial interaction with neighbouring neurons becomes essential (see, e.g., Eq. 13). It is the purpose of an equation-free method to reduce such a complex realistic description of the striatum to as few degrees of freedom as possible, by keeping the important dynamical signatures.

2.2 Modelling of the neuron dynamics

Our striatum network contains two types of neurons, the medium spiny neurons (MSN) representing 95% of all neurons and fast spiking neurons (FSI) which are the remaining ones. For the equations of motion of the neuron dynamics we follow (Chartove et al., 2020). It is reported therein, using models as well as experimental works, that striatal projection neurons (MSN) are capable of generating $\beta$ oscillations. In contrast, striatal fast-spiking interneurons (FSIs) are responsible for generating delta and theta rhythmicity (at 2–6 Hz). In this sense, the FS-neurons are somewhat paradigmatic for GABAeric interneurons in the striatum (that means, neurons which use neurotransmitter gamma-aminobutyric acid in synapses, mainly to inhibit other neurons), although, obviously, other types such as Somatostatin-expressing inhibitory interneurons (SOM+) exist (Melzer et al., 2017). We chose to model the striatum mainly with parvalbumin-positive fast spiking interneurons (PV+) (Melzer et al., 2017). On the one hand they are among the best characterised neurons (Tepper et al., 2010). On the other hand a recent study focusing on identifying interneurons in the striatum found that those neurons accounted for the largest group of interneurons overlapping with 5HT3-EGFP, the marker which turned out to best identify interneurons but otherwise did not very much overlap with classical markers (Muñoz-Manchado et al., 2016). In addition PV+ are found more prominently in the dorsal, whereas SOM+ are more prominently found in the ventral striatum (Zandt et al., 2024). Finally Cholinergic neurons, in turn, act only via metabotropic receptors and hence slower than the GABAergic ones.

An MSN or FS neuron at node $i$ is modelled by current balance equations for the membrane potential $V_i$

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

where $C$ is the membrane capacity. The current balance Eq. 1 contains four membrane currents, the fast sodium and potassium currents $I_{\text{Na}}$, $I_{\text{K}}$, the leak current $I_{\text{LEAK}}$. For MSN neurons an M-current $I_{\text{M}}$ occurs whereas FS neurons contain a D-current $I_{\text{D}}$ (Chartove et al., 2020). All currents follow the Hodgkin Huxley formalism Hodgkin and Huxley (1952), that means $I_X=g_X{m}_X^{n_1}{h}_X^{n_2}\cdot (V-E_X)$, where $X\in \{\text{Na,\hspace{0.17em}K,\hspace{0.17em}M,\hspace{0.17em}LEAK}\}$. The exponents $n_1$, $n_2$ represent the number of activation-inactivation channels respectively, $g_X$ is the maximum conductance of the ion and $E_X$ stands for the reverse potential for each ion. Specifically, the sodium current has three activation gates and one inactivation gate so that

$I_{\text{Na}}=g_{\text{Na}}{m}_{\text{Na}}^3{h}_{\text{Na}}\cdot \left(V_i-E_{\text{Na}}\right).(2)$

The potassium current has the form

$I_{\text{K}}=g_{\text{K}}{m}_{\text{K}}^4\cdot \left(V_i-E_{\text{K}}\right).(3)$

and the leak current reads

$I_{\text{LEAK}}=g_{\text{LEAK}}\cdot \left(V_i-E_{\text{LEAK}}\right).(4)$

Finally, the M- and D-current which enter the MSN and the FS neurons, respectively, are given by

$I_{\text{M}}=g_{\text{M}}{m}_{\text{M}}\cdot \left(V_i-E_{\text{K}}\right),I_{\text{D}}=g_{\text{D}}{m}_{\text{D}}^3{h}_{\text{D}}\cdot \left(V_i-E_{\text{D}}\right).(5)$

The gating variables $m_{\text{Na}}$, $h_{\text{Na}}$, $m_{\text{K}}$, $m_{\text{M}}$, at node $i$ each obey, following the Hodgkin-Huxley formalism (Hodgkin and Huxley, 1952; Chartove et al., 2020), an equation of the type

$\frac{d{x}_i}{dt}=a_{x_i}\left(V_i\right)-\left(a_{x_i}\left(V_i\right)+b_{x_i}\left(V_i\right)\right)x_i,(6)$

where $x_i\in \{m_{\text{Na}},h_{\text{Na}},m_{\text{K}},m_{\text{M}}\}$ stands for the respective gating variable at node $i$. The voltage dependent coefficients for the gating variables of the sodium current are given by

$a_{m_{\text{Na}}}\left(V\right)=0.32\frac{V+54}{1-e^{-\left(V+54\right)/4)}},b_{m_{\text{Na}}}\left(V\right)=0.28\frac{\left(V+27\right)}{e^{\left(V+27\right)/5}-1}$

and

$a_{h_{\text{Na}}}\left(V\right)=0.128e^{-\left(V+50\right)/18},b_{h_{\text{Na}}}\left(V\right)=\frac{4}{1+e^{-\left(V+27\right)/5}}.$

The coefficients of the activation gating for the potassium current read

$a_{m_{\text{K}}}\left(V\right)=0.032\frac{V+52}{1-e^{-\left(V+52\right)/5)}},b_{m_{\text{K}}}\left(V\right)=0.5e^{-\left(V+57\right)/40},$

and for those of the M-current we have

$a_{m_{\text{M}}}\left(V\right)=0.032\frac{V+52}{1-e^{-\left(V+52\right)/5)}},b_{m_{\text{M}}}\left(V\right)=0.5e^{-\left(V+57\right)/40}.$

The fast spiking neurons (FS) follow similar equations (Chartove et al., 2020), where instead of the M-current we use fast-activating, slowly inactivating D-current $I_{\text{D}}$, given in Eq. 5, with three activation gates and one inactivation gate, thus imposing a delay in firing upon depolarisation (Golomb et al., 2007; Chartove et al., 2020). Table 1 contains the respective parameter settings for both types of neurons.

Table 1

Table 1. Values for the conductance $g_X$ and inverse potential $E_X$ for the MSN and FS neurons.

The current $I_{\text{app}}$ in Eq. 1 is written as $I_{\text{app}}=I_0+I_{\text{DBS}}$, where $I_0$ predominantly represents a network activation current which describes the dependence of the neuronal activation due to intensity of cortical-striatal connectivity. The coupling between the neurons in Eq. 1 is described by the synaptic current $I_{\text{syn}}$. Details will be outlined in the next Section 2.3. Since our network model contains realistic spatial details about the actual neural system we would be able to model the impact of deep brain stimulation as well. Thus, the term $I_{\text{DBS}}$ representing the deep brain stimulation, enters here as an additive contribution. In our analysis we keep $I_{\text{DBS}}=0$, except the last section where we discuss the implementation of DBS in our model.

2.3 Description of the network inhibitory synaptic activity

We model the activation of a synapse using the activation variable $s_i$ for the $i$-th neuron (Compte et al., 2000; Laing and Chow, 2002; Ermentrout and Terman, 2012)

$\frac{d{s}_i}{dt}={\alpha }_X\left(1-s_i\right)H_X\left(V_i\right)-{\beta }_X{s}_i,(7)$

where $X\in \{\text{M,F}\}$ denotes whether the $i$-th neuron is a medium spiny neuron (M) or an interneuron (F), and $H_X(V)$ is a sigmoid function. The variable $s_i$ describes the activation of synapses from the pre-synaptic neuron $i$ to post-synaptic neurons. The parameters ${\alpha }_X$ and ${\beta }_X$ in Eq. 7 determine the activation and inactivation time scales, respectively, of the inhibitory synaptic connections. For MSN we choose, following (Chartove et al., 2020),

$H_{\text{M}}\left(V\right)=1+\tanh \left(V/4\right),$

with activation rates ${\alpha }_{\text{M}}=2$ and ${\beta }_{\text{M}}=1/13$. Similarly, for interneurons the expressions are given by

$H_{\text{F}}\left(V\right)=1+\tanh \left(V/10\right),$

with activation rates ${\alpha }_{\text{F}}=4$ and ${\beta }_{\text{F}}=1/13$. For a neuron $i$ of type $X\in \{\text{M,F}\}$ the total synaptic inhibition it receives from pre-synaptic neurons of type $Y\in \{\text{M,F}\}$ is given by

$I_{i,\text{GABA}}=g_{\text{X}\text{Y}}\left(V_i-E_{\text{GABA}}\right){\displaystyle \sum _{j\in Y}}{A}_{ij}{s}_j,(8)$

where $A_{ij}$ is the adjacency matrix of the graph, the summation $j\in Y$ is taken over neurons of type $Y$ and $E_{\text{GABA}}=-80\text{mV}$. The parameter $g_{\text{X}\text{Y}}$ represents the conductance between $X$ and $Y$ interactions with $X,Y\in \{\text{M},\text{F}\}$.

The synaptic current for the MSNs consists of two parts, first the sum of synaptic currents over medium spiny neurons (describing the inhibition between MSN-MSN neurons) and second, the sum over interneurons (interneurons inhibition of MSN), so that Eq. 8 yields

$I_{\mathit{syn}}=g_{\text{MM}}\left(V_i-E_{\text{GABA}}\right){\displaystyle \sum _{j\in \text{M}}}{A}_{ij}{s}_j+g_{\text{MF}}\left(V_i-E_{\text{GABA}}\right){\displaystyle \sum _{j\in \text{F}}}{A}_{ij}{s}_j.(9)$

Similarly for an interneuron the synaptic current is given by

$I_{\mathit{sys}}=g_{\text{FF}}\left(V_i-E_{\text{GABA}}\right){\displaystyle \sum _{j\in \text{F}}}{A}_{ij}{s}_j+g_{\text{FM}}\left(V_i-E_{\text{GABA}}\right){\displaystyle \sum _{j\in \text{M}}}{A}_{ij}{s}_j.(10)$

Here the first sum represents the rare case of FS-FS inhibition, while the second term governs the feedback inhibitory loop of MSN to interneurons. For the conductivity values we use $g_{\text{MM}}=g_{\text{MF}}=0.02$ and $g_{\text{FF}}=g_{\text{FM}}=0.005$.

In summary, Eqs 1, 2, 3, 4, 6, 7, 9, and 10 constitute a high dimensional heterogeneous set of coupled nonlinear differential equations defined on a graph with adjacency matrix $A$. The state of each neuron at node $i$ is described by the set of variables $(V_i,{(m_{\text{Na}})}_i,{(h_{\text{Na}})}_i,{(m_{\text{K}})}_i,{(m_{\text{M}})}_i,s_i)$. Figures 1B,C illustrates the temporal dynamics of the network.

3 Equation-free method for analysing macroscopic network behaviour

To describe the main idea in basic terms, consider a high-dimensional dynamical system, for instance the dynamics of the neural network presented in the previous section. The network model evolves in time under specified known microscopic rules, e.g., the equations of motion for each node described above. Denote by $U_t\in {\mathbb{R}}^N$ the state of the full network. Its time evolution over a time interval $T$ is given by

$U_{t+T}={\mathrm{\Phi }}_T\left(U_t\right)$

where the so-called flow ${\mathrm{\Phi }}_T:{\mathbb{R}}^N\to {\mathbb{R}}^N$ can be obtained from the (numerical) integration of the microscopic equations of motion.

We are interested in analysing the network behaviour on a different macroscopic scale. Assume there exists a suitable low-dimensional macroscopic variable $x\in {\mathbb{R}}^d$ with $d\ll N$, which captures the emergent collective behaviour of the network dynamics. Such collective coordinates $x_t=R(U_t)$ depend on the degree of freedoms of the system, and are determined by a restriction map $R:{\mathbb{R}}^N\to {\mathbb{R}}^d$. Often one uses suitable averages for the purpose to capture the dynamics of a system at a macroscopic scale, see, e.g., (Kozma et al., 2005; Kevrekidis and Samaey, 2009; Spiliotis and Siettos, 2011). The motion of the macroscopic variable per time step is then given by

$x_{t+T}=R\left({\mathrm{\Phi }}_T\left(U_t\right)\right)=\left(R◦{\mathrm{\Phi }}_T\right)\left(U_t\right).$

The main assumption of the Equation-Free approach (Gear et al., 2003; Kevrekidis and Samaey, 2009; Marschler et al., 2014a; Sieber et al., 2018) is that a macroscopic description in the form of closed equations of motion for the collective coordinate $x_t$ exists, even if one is unable to provide an analytic derivation for such effective low-dimensional equations of motion. The success of this idea relies on the ability to construct a meaningful lifting operator $L:{\mathbb{R}}^d\to {\mathbb{R}}^N$ such that

$U_t=L\left(x_t\right)$

gives for each macrostate $x_t$ a particular microstate $U_t$ which represents all the possible microstates which are consistent with $x_t$. As with respect to the restriction map, the lifting operator needs to fulfil the obvious consistency constraint $R◦L=I$ where $I$ denotes the identity map. Whether such a lift exists depends of course on the microscopic model, the emergence of collective motion, and a suitable choice for the macroscopic variable. If all these conditions are met one obtains the macroscopic evolution law

$x_{t+T}=R\left({\mathrm{\Phi }}_T\left(L\left(x_t\right)\right)\right)=\left(R◦{\mathrm{\Phi }}_T◦L\right)\left(x_t\right)=F_T\left(x_t\right),$

and the equations of motion, i.e., $F_T$, can be derived from a numerical integration of the microscopic equations of motion over short time intervals $T$, without the need of going through cumbersome calculations (hence the notion of an equation-free approach). If successful, one can finally utilise algorithms to compute dynamical signatures, such as fixed points and their stability, as well as complete bifurcation diagrams from the macroscopic description in terms of $F_T$, see for instance (Gear et al., 2003; Kevrekidis and Samaey, 2009; Siettos, 2011). An implicit equation-free analysis (Marschler et al., 2014b; Sieber et al., 2018) could further minimise the numerical errors in the computed bifurcation diagram.

4 Equation-free analysis of the striatum model

In this section, we apply the theoretical framework of an equation-free approach (Gear et al., 2005; Kevrekidis and Samaey, 2009; Marschler et al., 2014a) to analyse the emergent network dynamics macroscopically. Initially, we describe the lifting and restriction operator as well as the timestepper construction. Then, we discuss the consequences of the resulting one-dimensional evolution equation.

4.1 Lifting and restriction operator

The mean synaptic activity of MSNs turns out to be a suitable macroscopic variable

$S_t=\frac{1}{N}{\displaystyle \sum _{i\in \text{M}}}{s}_i\left(t\right).(11)$

While such a choice looks appealing, and can be justified with hindsight, one can also support this choice by a more subtle data analysis using for instance diffusion maps, a data-driven method for dimensional reduction and manifold learning (Coifman and Lafon, 2006; Nadler et al., 2006; Laing et al., 2010; Marschler et al., 2014a; Dsilva et al., 2018). Here we skip those technical details and take Eq. 11 as our reduction map.

The crucial step to build the timestepper is the lifting operator. The construction of the lifting operator is based on two steps. First, we record a microscopic realisation of the system from a previous simulation, i.e., we store all the microscopic variables after a short period of 20 ms. Then, in the second step, we assign synaptic variables to the 1856 MSNs in the following way: Given a mean synaptic activity $S_t$, we assign synaptic variables to the 1856 MSNs by randomisation, $s_i(t)=S_t+0.05\cdot Z_i$ where $Z_i\sim N(\mathrm{0,1})$ are uncorrelated normal random variables, while we keep the other microscopic degrees of freedom unchanged. Then, we numerically integrate the coupled Hodgkin-Huxley equations for all neurons for a (macroscopically) short time $T$, to derive the new network microstate $U_{t+T}={\mathrm{\Phi }}_T(L(S_t),I_0)$. The time $T$ is chosen such that the other variables become enslaved to the mean synaptic activity $S_t$. In fact, relatively short bursts on short time scales establish such a slaving relation (Gear et al., 2005; Spiliotis and Siettos, 2011). Constructing heuristically a lifting operator which results in a microscopic initialisation closer to the unknown attractive slow manifold leads to a faster convergence to the macroscopic behaviour and smaller numerical errors of the explicit equation-free analysis. As already mentioned, an implicit equation-free analysis (Marschler et al., 2014b; Sieber et al., 2018) could even further minimise the numerical errors in the computed bifurcation diagram.

Finally, we apply the restriction operator to the new network microstate $U_{t+T}$, i.e., we compute the mean synaptic activity $S_{t+T}$ which then defines the macroscopic evolution law $S_{t+T}=F_T(S_t,I_0)$. Here we have explicitly noted a constant network activation current $I_{\mathit{app}}=I_0$ as a static parameter of our model (see Eq. 1). Figure 2 shows a graphical representation of the numerical procedure.

Figure 2

Figure 2. Equation-free construction of the timestepper: Start with the macroscopic variable $S_t$, the mean synaptic activity. Transform this value into a consistent microscopic network state $U_t\in {\mathbb{R}}^n$ through a lifting operator. We simulate the network equations for the all neurons and for a short macroscopic time $T$ to derive the new network microstate $U_{t+T}$. Finally, average the synaptic variables $s_i$ to obtain the macroscopic variable $S_{t+T}$.

Since the macroscopic variable $S_t$ changes little on the time scale $T$ we can approximate the time discrete dynamics by a time continuous first order differential equation

$S^{\prime }\left(t\right)=f\left(S\left(t\right),I_0\right)\approx \frac{F_T\left(S_t,I_0\right)-S_t}{T}.(12)$

where the right hand side $F_T(S_t,I_0)$ in the difference quotient is determined by our equation free approach.

4.2 Data-driven system identification, stability and bifurcation analysis

Using Eq. 12 we can construct the right hand side $f$ of the macroscopic equation of motion. We perform independent parallel computations by covering the phase space with an equidistant mesh of initial values for the macroscopic variable, and the $I_0$ axis with an equidistant lattice of parameter values. We thus obtain the right hand side $f(S,I_0)$ with fairly high numerical resolution. The results for $f$ in dependence on $S$ are depicted in Figure 3, for $I_0=8$, 10, 12, 12.8, 13, and $I_0=13.2$. Despite a quite noisy neuron dynamics we obtain a rather smooth result for $f$ which shows little fluctuations. The computation of $f$ has been based on macroscopic averages over at about 2000 neurons and an ensemble average of 20 realisations, resulting in statistical errors of about $0.5\%$, in line with the data shown in Figure 3. The fixed points of the macroscopic dynamics are given by the zeros of the function $f$, while the sign of the slope at the zero determines the stability of the fixed point. The fixed point is stable for negative slope, while the fixed point is unstable for positive slope. Here stability refers to stability with respect to the macroscopic dynamics which is solely governed by the mean synaptic activity $S(t)$. While the internal microscopic dynamics is still highly complex, at the macroscopic level the motion is captured by the single scalar quantity $S(t)$. With our equation free approach we have been able to determine the right hand side of the macroscopic equation of motion (22), see Figure 3. Thus, the zeros of this right hand side and the sign of the slope allow us to determine the location and the macroscopic linear stability of the macroscopic stationary state.

Figure 3

Figure 3. Equation-free system identification: For different values of parameter $I_0$, we construct numerically the right hand side $f$ in dependence on the mean synaptic activity $S$. Red dots indicate the zeros of $f$, i.e., the fixed point solutions of the macroscopic dynamics. Clearly, the right hand side shows one fixed point at $I_0=8$, two fixed points in the range $I_0\in [\mathrm{10,13.2}]$ and finally, no fixed point for $I_0>13.2$.

We observe that the shape of right hand side $f$ is smooth and the graph shifts down, as the value of parameter $I_0$ increases. As a consequence the number of macroscopic fixed points changes. For $I_0=8$ we obtain one stable fixed point. As the value of $I_0$ increases, e.g., for $I_0=10$, we see two fixed points, an unstable one at small values $S^{*}=0.09$ and a second stable one at $S^{*}=0.73$ with high network activation. For increasing values of the parameter $I_0$ these two fixed points still persist with stability properties unchanged, but at a critical value close to $I_0=13.2$ the two fixed points collide and disappear in a saddle-node bifurcation. We can condense this information in a bifurcation diagram, see Figure 4. There are two branches of steady state solutions. The high neural activation solutions are stable (solid red line in Figure 4) while the low activation branch is unstable (dashed blue line in Figure 4). These two branches bifurcate in a saddle node bifurcation at $I_{\text{CRIT}}=13.19$. In general, an increased intensity of the current $I_0$ changes the rhythmicity and the density of activation. In the pathological case which corresponds to high activation, neurons exhibit spiking activity with variable periods (i.e., non-constant period between two spikes), and some neurons appear to show brief intervals of synchronised activity, preceded and followed by non-synchronous firing. Such synchrony could either be due to transient common activation via network inputs (e.g., inhibition of fast-spiking neurons), or it could actually occur by chance with this tonic firing at a relatively high frequency. The equation-free method remarkably reveals also an unstable low neural activation branch. Such an unstable state is not accessible by direct numerical simulations of the network model, it is a genuine outcome of the equation free approach. In terms of the microscopic dynamics such a state corresponds to an invariant saddle in the full phase space containing all microscopic degrees of freedom. For a potential neurophysiological interpretation of this state we recall that during the pathological case of obsessive-compulsive disorder, there is a hyperactivity of the striatum network. Thus, the stable high-activation branch of the solution can be seen as a pathological condition. The unstable low activation state that cannot be reached in direct simulations is nevertheless accessible by control techniques, such as closed loop deep brain stimulations. When successful, stabilising this unstable low activation state will produce a therapeutic effect on the striatum network hyperactivity.

Figure 4

Figure 4. Macroscopic system analysis of the striatum network: Bifurcation diagram as obtained from the equation-free analysis of the striatum network. The network activation current $I_{\mathit{app}}=I_0$ is used as the bifurcation parameter, and the mean synaptic activity of MSNs acts as macroscopic variable. Solid line (red) are stable fixed points, dashed line (blue) are unstable fixed points. The two branches disappear in a saddle-node bifurcation at $I_{\text{CRIT}}=13.19$. The insets show temporal simulations of mean synaptic activity S, for $I_0=\mathrm{8,10,12}$. Simulations converge to the upper stable branch of the bifurcation diagram.

5 Discussion

The recently invented field of network physiology aims at inferring dynamical interactions in complex biological or medical systems from observed data. With its inherently interdisciplinary intention this field aims to understand, based on data analysis, modelling approaches, or clinical practice, how diverse biological or physiological sub-systems interact from the cellular microscopic to the phenomenological macroscopic level, to explain diverse physiological phenomena, such as healthy or unhealthy states (see, e.g., (Schöll, 2022) for a recent editorial). Looking at the emerging field of network physiology from an equation free perspective has the potential to add an additional facet to this area of research. An equation free approach aims at uncovering the complex dynamical behaviour at a macroscopic level without the need to reconstruct the complex underlying microscopic dynamical network, thus addressing a main goal of network physiology from the outset. We have showcased a computational framework to analyse biophysical neuronal network models, and we applied the method to the striatum area. Based on a realistic mathematical model for the microscopic dynamics of the striatum we have been able to detect relevant macrostates and their dynamical features using an equation-free approach. One major contribution of this research work is that the method bridges the different levels of spatio-temporal scales, the microscopic ones where the physics of neurons is known and the macroscopic ones where the analysis is performed. The activity of neurons and the individual synaptic activity is given using the Hodgkin-Huxley formalism, which constitutes the microscopic description of the model. The network connects these neurons and produces a macroscopic or emergent behaviour with different spatio-temporal properties. Importantly, our equation-free approach allows us to study this emergent behaviour in detail, i.e., to perform stability and bifurcation analysis. The synaptic activity shows steady behaviour, which corresponds to the high network activity, the upper branch of solution in Figure 4, while the corresponding spectrum of the mean membrane activity shows a characteristic peak at the gamma band (see as well Figure 1D). Several other studies also analyse the macroscopic network activity or the emergent network behaviour (Fesce, 2024; Kromer and Tass, 2024; Venkadesh et al., 2024). Additionally, in (Kromer and Tass, 2024), a detailed study of mean synaptic activity, including synaptic plasticity, is performed. The proposed equation-free approach can be applied to these works containing multiple spatio-temporal scales. For example, by studying synaptic plasticity, one can extract critical values of synaptic strength, which contribute to phase transition in the macroscopic network dynamics (Marschler et al., 2014a).

Our realistic microscopic model was based on an FDA-approved state-of-the-art human atlas (Iacono et al., 2015) extracting coordinates for the striatal neurons, on modified Hodgkin-Huxley equations for medium spiny neurons (MSN) and fast-spiking neurons (FSN) (Hodgkin and Huxley, 1952; Chartove et al., 2020), and on complex network structures for neuronal connectivity (Netoff et al., 2004; Bassett and Bullmore, 2006; Bassett and Bullmore, 2017; de Santos-Sierra et al., 2014; Berman et al., 2016; She et al., 2016; Fang et al., 2017). Depending on the parameters, the network model produces patterns which can be associated with healthy or pathological conditions, reflected by low or high synaptic activity. In clinical studies of obsessive-compulsive disorder (Maltby et al., 2005; Marsh et al., 2014) a dysfunctional hyperactivity of the frontal-striatal circuits is observed similarly to the high activation state we obtain in our model for increasing the intensity of the cortico-striatal current $I_0$.

Within an equation-free approach we were able to investigate the crucial macroscopic behaviour for the mean synaptic activity. Such an analysis not just reproduces the dynamically stable high activity branch, but also shows an unstable low activity state which is inaccessible by direct simulations of the model. Such unstable dynamical states could be promising targets for treating pathological conditions.

Deep brain stimulation (DBS) of the striatum has evolved as a promising therapy for patients with severe and resistant forms of obsessive compulsive disorders (OCD) and mental impairments (Rodriguez-Romaguera et al., 2012; Blomstedt et al., 2013; Widge et al., 2019; Wu et al., 2021). While there exist different computational approaches modelling DBS for OCD, see for instance (Szalisznyó and Silverstein, 2021), we can utilise our realistic large scale dynamical system to obtain insights about pathological neural activity during OCD. Since our model has been based on the realistic spatial structures of the striatum each neuron, labelled by an index $i$, comes with its corresponding spatial position, given by a vector ${\vec{x}}_i$. The impact of the current applied in DBS and acting on the $i$-th neuron depends on the distance between the position of the electrode in real space, denoted by ${\vec{x}}_E$, and the position of the $i$-th neuron. It is modelled by

$I_{\text{DBS}}=A_{\text{DBS}}\exp \left(-\Vert {\vec{x}}_i-{\vec{x}}_E{\Vert }^2/{\sigma }^2\right)H\left(\sin \left({\omega }_{\text{DBS}}t\right)\right)\cdot \left(1-H\left(\sin \left({\omega }_{\text{DBS}}t+{\delta }_{\text{DBS}}\right)\right)\right.(13)$

This quantity enters the equation for the $i$-th neuron in an additive way, $I_{\text{app}}=I_0+I_{\mathit{DBS}}$, see Eq. 1. The impact of DBS decreases super-exponentially with respect to the distance from the electrode. Far from the electrode, the intensity current is almost zero, having no impact on the neuron action, and different electrode positions result in different network activations. Here ${\omega }_{\text{DBS}}=2\pi /T_{\text{DBS}}$ denotes the angular frequency, $T_{\text{DBS}}$ the period, and the phase shift ${\delta }_{\text{DBS}}$ the duration of the pulse, with $H$ abbreviating the Heaviside function ($H(x)=1$, if $x>0$, else $H(x)=0$). The function $H$ takes into account that DBS is not applied as a plain harmonic, but in terms of periodic pulses. The crucial amplitude and frequency in DBS is determined by the two parameters $A_{\text{DBS}}$ and ${\omega }_{\text{DBS}}$, respectively. Figure 5 shows the application of DBS at ${\vec{x}}_E=(-\mathrm{9,9,5})$, with the DBS amplitude and frequency $A_{\text{DBS}}=200$ and ${\omega }_{\text{DBS}}=2\pi /130$. We observe that DBS induces strong synchronisation in the neural activity of striatum.

Figure 5

Figure 5. Deep brain stimulation (DBS) on the striatum model: Simulation of the network model with the current Eq. 13 added to the network equations.

Thanks to the equation-free framework we are now able to design a macroscopic proportional feedback controller for DBS. For instance, for $I_0=10$, the equation-free analysis showed the existence of one unstable fixed point at mean synaptic activity $S^{*}=0.08$. We use the amplitude of DBS, that means $A_{\text{DBS}}$, as control variable which is adjusted due to linear proportional feedback

$\frac{d{A}_{\text{DBS}}}{dt}=-K_p\left(S-S^{*}\right)(14)$

where $K_p$ denotes the gain of the control. By choosing the gain appropriately we aim at driving the system towards the low activation state. Figure 6 shows the application of DBS at the point ${\vec{x}}_E=(-\mathrm{9,9,5})$ and for frequency 200Hz. We observe that after switching on the feedback control ($t>150$ms) the macroscopic activity gets closer to the healthy low activation state, see Figure 6B, and that synchronisation is destroyed in favour of a desynchronised state, see Figure 6A. In general, explaining the mechanism of DBS and how it acts in the evolved brain network is still a mystery. For example, in Parkinson’s disease, it is unclear whether DBS suppresses or enhances the neural activity of the targeted areas (Rubin and Terman, 2004; Montgomery and Gale, 2008). In Figure 6, we present two stages of DBS: the first 150 ms without a control scheme and the second part after 150 ms. While DBS without control induces synchronised activity of neurons such a synchronised state is suppressed when control is turned on. In that respect the closed-loop DBS results in realistic patterns closer to healthy conditions.

Figure 6

Figure 6. Closed loop control scheme for DBS on the striatum network: Application of DBS with constant amplitude along the line of Eq. 13 for $t<150\text{ms}$ (red). Close loop control scheme for DBS, using Eq. 14, adjusting the DBS amplitude by linear proportional feedback for $t>150\text{ms}$. (A) Raster plot for $n=500$ randomly chosen neurons. Black dots represent activated neurons (i.e., time dependent action potentials passing through −15 mV towards positive values. (B) The mean synaptic activity for DBS without control (red), and DBS with linear proportional feedback (blue).

There are still considerable unknowns for a successful application of DBS such as the anatomical targets of stimulation, optimal stimulation parameters like amplitude and frequency of stimulation, as well as long-term effects of stimulation. In obsessive compulsive disorders hyperactive frontal-striatal activity has been reported (Maltby et al., 2005; Marsh et al., 2014). We conjecture that this hyperactivity is qualitatively similar to the stable upper branch solution as depicted in the bifurcation diagram Figure 4. Since our network model allows for properly modelling the network activation current a corresponding equation-free analysis of the model may then provide some answers to the open questions raised above. Our successful simple showcase provides evidence that such an ambitious program may succeed.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

Author contributions

KS: Formal Analysis, Investigation, Methodology, Software, Writing–original draft, Writing–review and editing. RK: Conceptualization, Investigation, Methodology, Writing–review and editing. WJ: Methodology, Writing–review and editing. JS: Conceptualization, Funding acquisition, Methodology, Supervision, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - SFB 1270/2 - 299150580 - Collaborative Research Centre ELAINE.

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.

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, D. S., and Bullmore, E. T. (2017). Small-world brain networks revisited. Neurosci. 23, 499–516. PMID: 27655008. doi:10.1177/1073858416667720

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., Brincat, S. L., Lundqvist, M., and Miller, E. K. (2022). Traveling waves in the prefrontal cortex during working memory. PLOS Comput. Biol. 18, e1009827. doi:10.1371/journal.pcbi.1009827

PubMed Abstract | CrossRef Full Text | Google Scholar

Blomstedt, P., Sjöberg, R. L., Hansson, M., Bodlund, O., and Hariz, M. I. (2013). Deep brain stimulation in the treatment of obsessive-compulsive disorder. World Neurosurg. 80, e245–e253. doi:10.1016/j.wneu.2012.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G. (2004). Large-scale recording of neuronal ensembles. Nat. Neurosci. 7, 446–451. doi:10.1038/nn1233

PubMed Abstract | CrossRef Full Text | Google Scholar

Calabresi, P., Picconi, B., Tozzi, A., and Di Filippo, M. (2007). Dopamine-mediated regulation of corticostriatal synaptic plasticity. Trends Neurosci. 30, 211–219. Fifty years of dopamine research. doi:10.1016/j.tins.2007.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Calabresi, P., Picconi, B., Tozzi, A., Ghiglieri, V., and Di Filippo, M. (2014). Direct and indirect pathways of basal ganglia: a critical reappraisal. Nat. Neurosci. 17, 1022–1030. doi:10.1038/nn.3743

PubMed Abstract | CrossRef Full Text | Google Scholar

Chartove, J. A. K., McCarthy, M. M., Pittman-Polletta, B. R., and Kopell, N. J. (2020). A biophysical model of striatal microcircuits suggests gamma and beta oscillations interleaved at delta/theta frequencies mediate periodicity in motor control. PLOS Comput. Biol. 16, e1007300–e1007330. doi:10.1371/journal.pcbi.1007300

PubMed Abstract | CrossRef Full Text | Google Scholar

Coifman, R. R., and Lafon, S. (2006). Diffusion maps. Appl. Comput. Harmon. Analysis 21, 5–30. Special Issue: Diffusion Maps and Wavelets. doi:10.1016/j.acha.2006.04.006

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Crittenden, J., and Graybiel, A. (2011). Basal ganglia disorders associated with imbalances in the striatal striosome and matrix compartments. Front. Neuroanat. 5, 59. doi:10.3389/fnana.2011.00059

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

de Santos-Sierra, D., Sendiña-Nadal, I., Leyva, I., Almendral, J. A., Anava, S., Ayali, A., et al. (2014). Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures. Plos One 9, 1–8. doi:10.1371/journal.pone.0085828

CrossRef Full Text | Google Scholar

Dsilva, C. J., Talmon, R., Coifman, R. R., and Kevrekidis, I. G. (2018). Parsimonious representation of nonlinear dynamical systems through manifold learning: a chemotaxis case study. Appl. Comput. Harmon. Analysis 44, 759–773. doi:10.1016/j.acha.2015.06.008

CrossRef Full Text | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Fesce, R. (2024). The emergence of identity, agency and consciousness from the temporal dynamics of neural elaboration. Front. Netw. Physiology 4, 1292388. doi:10.3389/fnetp.2024.1292388

PubMed Abstract | CrossRef Full Text | Google Scholar

Gear, C., Kaper, T., Kevrekidis, I., and Zagaris, A. (2005). Projecting to a slow manifold: singularly perturbed systems and legacy codes. SIAM J. Appl. Dyn. Syst. 4, 711–732. doi:10.1137/040608295

CrossRef Full Text | Google Scholar

Gear, C. W., Hyman, J. M., Kevrekidid, P. G., Kevrekidis, I. G., Runborg, O., and Theodoropoulos, C. (2003). Equation-free, coarse-grained multiscale computation: enabling mocroscopic simulators to perform system-level analysis. Commun. Math. Sci. 1, 715–762. doi:10.4310/cms.2003.v1.n4.a5

CrossRef Full Text | Google Scholar

Golomb, D., Donner, K., Shacham, L., Shlosberg, D., Amitai, Y., and Hansel, D. (2007). Mechanisms of firing patterns in fast-spiking cortical interneurons. PLOS Comput. Biol. 3, e156. doi:10.1371/journal.pcbi.0030156

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiology 117, 500–544. doi:10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Iacono, M. I., Neufeld, E., Akinnagbe, E., Bower, K., Wolf, J., Vogiatzis Oikonomidis, I., et al. (2015). Mida: a multimodal imaging-based detailed anatomical model of the human head and neck. PLOS ONE 10, e0124126. doi:10.1371/journal.pone.0124126

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalenscher, T., Lansink, C. S., Lankelma, J. V., and Pennartz, C. M. A. (2010). Reward-associated gamma oscillations in ventral striatum are regionally differentiated and modulate local firing activity. J. Neurophysiol. 103, 1658–1672. doi:10.1152/jn.00432.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Kevrekidis, I., and Samaey, G. (2009). Equation-free multiscale computation: algorithms and applications. Annu. Rev. Phys. Chem. 60, 321–344. doi:10.1146/annurev.physchem.59.032607.093610

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozma, R., Puljic, M., Balister, P., Bollobás, B., and Freeman, W. J. (2005). Phase transitions in the neuropercolation model of neural populations with mixed local and non-local interactions. Neurocomputing 92, 367–379. doi:10.1007/s00422-005-0565-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Kromer, J. A., and Tass, P. A. (2024). Coordinated reset stimulation of plastic neural networks with spatially dependent synaptic connections. Front. Netw. Physiology 4, 1351815. doi:10.3389/fnetp.2024.1351815

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2018). On the application of “equation-free modelling” to neural systems. J. Comput. Neurosci. 20, 5–23. doi:10.1007/s10827-006-3843-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R., Frewen, T., and Kevrekidis, I. G. (2010). Reduced models for binocular rivalry. J. Comput. Neurosci. 28, 459–476. doi:10.1007/s10827-010-0227-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltby, N., Tolin, D. F., Worhunsky, P., O’Keefe, T. M., and Kiehl, K. A. (2005). Dysfunctional action monitoring hyperactivates frontal-striatal circuits in obsessive-compulsive disorder: an event-related fmri study. NeuroImage 24, 495–503. doi:10.1016/j.neuroimage.2004.08.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Marschler, C., Faust-Ellsässer, C., Starke, J., and Van Hemmen, J. (2014a). Bifurcation of learning and structure formation in neuronal maps. EPL 108, 48005. doi:10.1209/0295-5075/108/48005

CrossRef Full Text | Google Scholar

Marschler, C., Sieber, J., Berkemer, R., Kawamoto, A., and Starke, J. (2014b). Implicit methods for equation-free analysis: convergence results and analysis of emergent waves in microscopic traffic models. SIAM J. Appl. Dyn. Syst. 13, 1202–1238. doi:10.1137/130913961

CrossRef Full Text | Google Scholar

Marsh, R., Horga, G., Parashar, N., Wang, Z., Peterson, B. S., and Simpson, H. B. (2014). Altered activation in fronto-striatal circuits during sequential processing of conflict in unmedicated adults with obsessive-compulsive disorder. Biol. Psychiatry 75, 615–622. Obsessive-Compulsive Disorder and the Connectome. doi:10.1016/j.biopsych.2013.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Melzer, S., Gil, M., Koser, D. E., Michael, M., Huang, K. W., and Monyer, H. (2017). Distinct corticostriatal gabaergic neurons modulate striatal output neurons and motor activity. Cell Rep. 19, 1045–1055. doi:10.1016/j.celrep.2017.04.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5, 021028. doi:10.1103/PhysRevX.5.021028

CrossRef Full Text | Google Scholar

Montgomery, E. B., and Gale, J. T. (2008). Mechanisms of action of deep brain stimulation (dbs). Neurosci. Biobehav. Rev. 32, 388–407. Special section: Neurobiology of Deep Brain Stimulation: Innovations in Treatment and Basal Ganglia Function. doi:10.1016/j.neubiorev.2007.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Manchado, A., Foldi, C., Szydlowski, S., Sjulson, L., Farries, M., Wilson, C., et al. (2016). Novel striatal gabaergic interneuron populations labeled in the 5ht3a(egfp) mouse. Cereb. Cortex 26, 96–105. doi:10.1093/cercor/bhu179

PubMed Abstract | CrossRef Full Text | Google Scholar

Nadler, B., Lafon, S., Coifman, R. R., and Kevrekidis, I. G. (2006). Diffusion maps, spectral clustering and reaction coordinates of dynamical systems. Appl. Comput. Harmon. Analysis 21, 113–127. Special Issue: Diffusion Maps and Wavelets. doi:10.1016/j.acha.2005.07.004

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Palkar, G., Wu, J.-y., and Ermentrout, B. (2023). The inhibitory control of traveling waves in cortical networks. PLOS Comput. Biol. 19, e1010697. doi:10.1371/journal.pcbi.1010697

PubMed Abstract | CrossRef Full Text | Google Scholar

Parasuram, H., Nair, B., D’Angelo, E., Hines, M., Naldi, G., and Diwakar, S. (2016). Computational modeling of single neuron extracellular electric potentials and network local field potentials using lfpsim. Front. Comput. Neurosci. 10, 65. doi:10.3389/fncom.2016.00065

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Remijnse, P. L., Nielen, M. M. A., van Balkom, A. J. L. M., Cath, D. C., van Oppen, P., Uylings, H. B. M., et al. (2006). Reduced orbitofrontal-striatal activity on a reversal learning task in obsessive-compulsive disorder. Archives General Psychiatry 63, 1225–1236. doi:10.1001/archpsyc.63.11.1225

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Romaguera, J., Monte, F. H. M. D., and Quirk, G. J. (2012). Deep brain stimulation of the ventral striatum enhances extinction of conditioned fear. Proc. Natl. Acad. Sci. 109, 8764–8769. doi:10.1073/pnas.1200782109

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Schöll, E. (2022). Network physiology, insights in dynamical systems: 2021. Front. Netw. Physiol. 2, 1–3. doi:10.3389/fnetp.2022.961339

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Sieber, J., Marschler, C., and Starke, J. (2018). Convergence of equation-free methods in the case of finite time scale separation with application to deterministic and stochastic systems. SIAM J. Appl. Dyn. Syst. 17, 2574–2614. doi:10.1137/17M1126084

CrossRef Full Text | Google Scholar

Siettos, C., and Starke, J. (2016). Multiscale modeling of brain dynamics: from single neurons and networks to mathematical tools. Wiley Interdiscip. Rev. Syst. Biol. Med. 8, 438–458. doi:10.1002/wsbm.1348

PubMed Abstract | CrossRef Full Text | Google Scholar

Siettos, C. I. (2011). Equation-free multiscale computational analysis of individual-based epidemic dynamics on networks. Appl. Math. Comput. 218, 324–336. doi:10.1016/j.amc.2011.05.067

CrossRef Full Text | Google Scholar

Spiliotis, K., Butenko, K., Starke, J., van Rienen, U., and Köhling, R. (2024). Towards an optimised deep brain stimulation using a large-scale computational network and realistic volume conductor model. J. Neural Eng. 20, 066045. doi:10.1088/1741-2552/ad0e7c

CrossRef Full Text | Google Scholar

Spiliotis, K., Butenko, K., van Rienen, U., Starke, J., and Köhling, R. (2022a). Complex network measures reveal optimal targets for deep brain stimulation and identify clusters of collective brain dynamics. Front. Phys. 10. doi:10.3389/fphy.2022.951724

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Szalisznyó, K., and Silverstein, D. N. (2021). Computational predictions for ocd pathophysiology and treatment: a review. Front. Psychiatry 12, 687062. doi:10.3389/fpsyt.2021.687062

PubMed Abstract | CrossRef Full Text | Google Scholar

Tepper, J. M., Tecuapetla, F., Koós, T., and Ibáñez-Sandoval, O. (2010). Heterogeneity and diversity of striatal gabaergic interneurons. Front. Neuroanat. 4 (150), 1–18. doi:10.3389/fnana.2010.00150

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Venkadesh, S., Shaikh, A., Shakeri, H., Barreto, E., and Van Horn, J. D. (2024). Biophysical modulation and robustness of itinerant complexity in neuronal networks. Front. Netw. Physiology 4, 1302499. doi:10.3389/fnetp.2024.1302499

CrossRef Full Text | Google Scholar

Watts, D., and Strogatz, S. (1998). Collective dynamics of 'small-world' networks. Nature 393, 440–442. doi:10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

Widge, A. S., Zorowitz, S., Basu, I., Paulk, A. C., Cash, S. S., Eskandar, E. N., et al. (2019). Deep brain stimulation of the internal capsule enhances human cognitive control and prefrontal cortex function. Nat. Commun. 10, 1536. doi:10.1038/s41467-019-09557-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H., Hariz, M., Visser-Vandewalle, V., Zrinzo, L., Coenen, V. A., Sheth, S. A., et al. (2021). Deep brain stimulation for refractory obsessive-compulsive disorder (ocd): emerging or established therapy? Mol. psychiatry 26, 60–65. doi:10.1038/s41380-020-00933-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yager, L., Garcia, A., Wunsch, A., and Ferguson, S. (2015). The ins and outs of the striatum: role in drug addiction. Neuroscience 301, 529–541. doi:10.1016/j.neuroscience.2015.06.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Zandt, M. V., Flanagan, D., and Pittenger, C. (2024). Sex differences in the distribution and density of regulatory interneurons in the striatum. BioRxiv [Preprint]. doi:10.1101/2024.02.29.582798

CrossRef Full Text | Google Scholar