Abstract
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., (; 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 (; 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 (; Spiliotis and Siettos, 2011; ), as well as travelling waves and spatio-temporal chaos (; ; Palkar et al., 2023), are signatures of the rich nonlinear behaviour of neural networks at the macroscopic level (; Spiliotis and Siettos, 2011; ; ; ; 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 (; ; ). 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; ).
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 (; ; 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); ; 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 () and transform into the MNI (Montreal Neurological Institute) coordinate system. We place neurons randomly inside this area, see Figure 1A.
FIGURE 1
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; ; Stam and Reijneveld, 2007; ; Spiliotis and Siettos, 2011). Small-world structures suitably model physiological networks (Netoff et al., 2004; ; ; ; ; She et al., 2016; ) 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 neurons in the local neighbourhood of 5 diameter. Then, with a small probability , 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 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 , where is the set of nodes and the set of edges, i.e., the connections between neurons. The connectivity can be captured by a so-called adjacency matrix , where if there is a connection from node to node , and 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 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 (). It is reported therein, using models as well as experimental works, that striatal projection neurons (MSN) are capable of generating 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 is modelled by current balance equations for the membrane potential where is the membrane capacity. The current balance Eq. 1 contains four membrane currents, the fast sodium and potassium currents , , the leak current . For MSN neurons an M-current occurs whereas FS neurons contain a D-current (). All currents follow the Hodgkin Huxley formalism , that means , where . The exponents , represent the number of activation-inactivation channels respectively, is the maximum conductance of the ion and stands for the reverse potential for each ion. Specifically, the sodium current has three activation gates and one inactivation gate so thatThe potassium current has the formand the leak current readsFinally, the M- and D-current which enter the MSN and the FS neurons, respectively, are given by
The gating variables , , , , at node each obey, following the Hodgkin-Huxley formalism (; ), an equation of the typewhere stands for the respective gating variable at node . The voltage dependent coefficients for the gating variables of the sodium current are given byandThe coefficients of the activation gating for the potassium current readand for those of the M-current we haveThe fast spiking neurons (FS) follow similar equations (), where instead of the M-current we use fast-activating, slowly inactivating D-current , given in Eq. 5, with three activation gates and one inactivation gate, thus imposing a delay in firing upon depolarisation (; ). Table 1 contains the respective parameter settings for both types of neurons.
TABLE 1
| Parameters and | MSN | FS |
|---|---|---|
| 80 | ||
| 100 | ||
| 1.3 | - | |
| - | ||
| −67 mV | −70 mV | |
| −100 mV | −90 mV | |
| 50 mV | 50 mV | |
| −100 mV | - | |
| - | −90 mV |
Values for the conductance and inverse potential for the MSN and FS neurons.
The current in Eq. 1 is written as , where 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 . 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 representing the deep brain stimulation, enters here as an additive contribution. In our analysis we keep , 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 for the -th neuron (; ; )where denotes whether the -th neuron is a medium spiny neuron (M) or an interneuron (F), and is a sigmoid function. The variable describes the activation of synapses from the pre-synaptic neuron to post-synaptic neurons. The parameters and in Eq. 7 determine the activation and inactivation time scales, respectively, of the inhibitory synaptic connections. For MSN we choose, following (),with activation rates and . Similarly, for interneurons the expressions are given bywith activation rates and . For a neuron of type the total synaptic inhibition it receives from pre-synaptic neurons of type is given bywhere is the adjacency matrix of the graph, the summation is taken over neurons of type and . The parameter represents the conductance between and interactions with .
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
Similarly for an interneuron the synaptic current is given byHere 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 and .
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 . The state of each neuron at node is described by the set of variables . 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 the state of the full network. Its time evolution over a time interval is given bywhere the so-called flow 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 with , which captures the emergent collective behaviour of the network dynamics. Such collective coordinates depend on the degree of freedoms of the system, and are determined by a restriction map . Often one uses suitable averages for the purpose to capture the dynamics of a system at a macroscopic scale, see, e.g., (; ; Spiliotis and Siettos, 2011). The motion of the macroscopic variable per time step is then given byThe main assumption of the Equation-Free approach (; ; 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 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 such thatgives for each macrostate a particular microstate which represents all the possible microstates which are consistent with . As with respect to the restriction map, the lifting operator needs to fulfil the obvious consistency constraint where 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 lawand the equations of motion, i.e., , can be derived from a numerical integration of the microscopic equations of motion over short time intervals , 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 , see for instance (; ; 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 (; ; 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 variableWhile 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 (; Nadler et al., 2006; Laing et al., 2010; Marschler et al., 2014a; ). 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 , we assign synaptic variables to the 1856 MSNs by randomisation, where 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 , to derive the new network microstate . The time is chosen such that the other variables become enslaved to the mean synaptic activity . In fact, relatively short bursts on short time scales establish such a slaving relation (; 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 , i.e., we compute the mean synaptic activity which then defines the macroscopic evolution law . Here we have explicitly noted a constant network activation current as a static parameter of our model (see Eq. 1). Figure 2 shows a graphical representation of the numerical procedure.
FIGURE 2
Since the macroscopic variable changes little on the time scale we can approximate the time discrete dynamics by a time continuous first order differential equationwhere the right hand side 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 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 axis with an equidistant lattice of parameter values. We thus obtain the right hand side with fairly high numerical resolution. The results for in dependence on are depicted in Figure 3, for , 10, 12, 12.8, 13, and . Despite a quite noisy neuron dynamics we obtain a rather smooth result for which shows little fluctuations. The computation of has been based on macroscopic averages over at about 2000 neurons and an ensemble average of 20 realisations, resulting in statistical errors of about , in line with the data shown in Figure 3. The fixed points of the macroscopic dynamics are given by the zeros of the function , 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 . While the internal microscopic dynamics is still highly complex, at the macroscopic level the motion is captured by the single scalar quantity . 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
We observe that the shape of right hand side is smooth and the graph shifts down, as the value of parameter increases. As a consequence the number of macroscopic fixed points changes. For we obtain one stable fixed point. As the value of increases, e.g., for , we see two fixed points, an unstable one at small values and a second stable one at with high network activation. For increasing values of the parameter these two fixed points still persist with stability properties unchanged, but at a critical value close to 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 . In general, an increased intensity of the current 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
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 (; ; Venkadesh et al., 2024). Additionally, in (), 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 () extracting coordinates for the striatal neurons, on modified Hodgkin-Huxley equations for medium spiny neurons (MSN) and fast-spiking neurons (FSN) (; ), and on complex network structures for neuronal connectivity (Netoff et al., 2004; ; ; ; ; She et al., 2016; ). 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 .
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; ; 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 , comes with its corresponding spatial position, given by a vector . The impact of the current applied in DBS and acting on the -th neuron depends on the distance between the position of the electrode in real space, denoted by , and the position of the -th neuron. It is modelled byThis quantity enters the equation for the -th neuron in an additive way, , 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 denotes the angular frequency, the period, and the phase shift the duration of the pulse, with abbreviating the Heaviside function (, if , else ). The function 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 and , respectively. Figure 5 shows the application of DBS at , with the DBS amplitude and frequency and . We observe that DBS induces strong synchronisation in the neural activity of striatum.
FIGURE 5
Thanks to the equation-free framework we are now able to design a macroscopic proportional feedback controller for DBS. For instance, for , the equation-free analysis showed the existence of one unstable fixed point at mean synaptic activity . We use the amplitude of DBS, that means , as control variable which is adjusted due to linear proportional feedbackwhere 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 and for frequency 200Hz. We observe that after switching on the feedback control (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
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.
Statements
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
1
BassettD.BullmoreE. (2006). Small-world brain networks. Neuroscientist12, 512–523. 10.1177/1073858406293182
2
BassettD. S.BullmoreE. T. (2017). Small-world brain networks revisited. Neurosci.23, 499–516. PMID: 27655008. 10.1177/1073858416667720
3
BermanB.SmucnyJ.WylieK.SheltonE.KronbergE.LeeheyM.et al (2016). Levodopa modulates small-world architecture of functional brain networks in Parkinson’s disease. Mov. Disord.31, 1676–1684. 10.1002/mds.26713
4
BhattacharyaS.BrincatS. L.LundqvistM.MillerE. K. (2022). Traveling waves in the prefrontal cortex during working memory. PLOS Comput. Biol.18, e1009827. 10.1371/journal.pcbi.1009827
5
BlomstedtP.SjöbergR. L.HanssonM.BodlundO.HarizM. I. (2013). Deep brain stimulation in the treatment of obsessive-compulsive disorder. World Neurosurg.80, e245–e253. 10.1016/j.wneu.2012.10.006
6
BullmoreE.SpornsO. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci.10, 186–198. 10.1038/nrn2575
7
BuzsákiG. (2004). Large-scale recording of neuronal ensembles. Nat. Neurosci.7, 446–451. 10.1038/nn1233
8
CalabresiP.PicconiB.TozziA.Di FilippoM. (2007). Dopamine-mediated regulation of corticostriatal synaptic plasticity. Trends Neurosci.30, 211–219. Fifty years of dopamine research. 10.1016/j.tins.2007.03.001
9
CalabresiP.PicconiB.TozziA.GhiglieriV.Di FilippoM. (2014). Direct and indirect pathways of basal ganglia: a critical reappraisal. Nat. Neurosci.17, 1022–1030. 10.1038/nn.3743
10
ChartoveJ. A. K.McCarthyM. M.Pittman-PollettaB. R.KopellN. 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. 10.1371/journal.pcbi.1007300
11
CoifmanR. R.LafonS. (2006). Diffusion maps. Appl. Comput. Harmon. Analysis21, 5–30. Special Issue: Diffusion Maps and Wavelets. 10.1016/j.acha.2006.04.006
12
CompteA.BrunelN.Goldman-RakicP.WangX.-J. (2000). Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cereb. Cortex10, 910–923. 10.1093/cercor/10.9.910
13
CrittendenJ.GraybielA. (2011). Basal ganglia disorders associated with imbalances in the striatal striosome and matrix compartments. Front. Neuroanat.5, 59. 10.3389/fnana.2011.00059
14
CrowellA.Ryapolova-WebbE.OstremJ.GalifianakisN.ShimamotoS.LimD.et al (2012). Oscillations in sensorimotor cortex in movement disorders: an electrocorticography study. Brain135, 615–630. 10.1093/brain/awr332
15
DecoG.JirsaV.RobinsonP.BreakspearM.FristonK. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol.4, e1000092. 10.1371/journal.pcbi.1000092
16
DecoG.Ponce-AlvarezA.MantiniD.RomaniG. L.HagmannP.CorbettaM. (2013). Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J. Neurosci.33, 11239–11252. 10.1523/JNEUROSCI.1091-13.2013
17
de Santos-SierraD.Sendiña-NadalI.LeyvaI.AlmendralJ. A.AnavaS.AyaliA.et al (2014). Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures. Plos One9, 1–8. 10.1371/journal.pone.0085828
18
DsilvaC. J.TalmonR.CoifmanR. R.KevrekidisI. G. (2018). Parsimonious representation of nonlinear dynamical systems through manifold learning: a chemotaxis case study. Appl. Comput. Harmon. Analysis44, 759–773. 10.1016/j.acha.2015.06.008
19
ErmentroutB.TermanD. (2012). Neural networks as spatio-temporal pattern-forming systems. New York: Springer.
20
FangJ.ChenH.CaoZ.JiangY.MaL.MaH.et al (2017). Impaired brain network architecture in newly diagnosed Parkinson’s disease based on graph theoretical analysis. Neurosci. Lett.657, 151–158. 10.1016/j.neulet.2017.08.002
21
FesceR. (2024). The emergence of identity, agency and consciousness from the temporal dynamics of neural elaboration. Front. Netw. Physiology4, 1292388. 10.3389/fnetp.2024.1292388
22
GearC.KaperT.KevrekidisI.ZagarisA. (2005). Projecting to a slow manifold: singularly perturbed systems and legacy codes. SIAM J. Appl. Dyn. Syst.4, 711–732. 10.1137/040608295
23
GearC. W.HymanJ. M.KevrekididP. G.KevrekidisI. G.RunborgO.TheodoropoulosC. (2003). Equation-free, coarse-grained multiscale computation: enabling mocroscopic simulators to perform system-level analysis. Commun. Math. Sci.1, 715–762. 10.4310/cms.2003.v1.n4.a5
24
GolombD.DonnerK.ShachamL.ShlosbergD.AmitaiY.HanselD. (2007). Mechanisms of firing patterns in fast-spiking cortical interneurons. PLOS Comput. Biol.3, e156. 10.1371/journal.pcbi.0030156
25
HodgkinA. L.HuxleyA. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiology117, 500–544. 10.1113/jphysiol.1952.sp004764
26
IaconoM. I.NeufeldE.AkinnagbeE.BowerK.WolfJ.Vogiatzis OikonomidisI.et al (2015). Mida: a multimodal imaging-based detailed anatomical model of the human head and neck. PLOS ONE10, e0124126. 10.1371/journal.pone.0124126
27
KalenscherT.LansinkC. S.LankelmaJ. V.PennartzC. M. A. (2010). Reward-associated gamma oscillations in ventral striatum are regionally differentiated and modulate local firing activity. J. Neurophysiol.103, 1658–1672. 10.1152/jn.00432.2009
28
KevrekidisI.SamaeyG. (2009). Equation-free multiscale computation: algorithms and applications. Annu. Rev. Phys. Chem.60, 321–344. 10.1146/annurev.physchem.59.032607.093610
29
KozmaR.PuljicM.BalisterP.BollobásB.FreemanW. J. (2005). Phase transitions in the neuropercolation model of neural populations with mixed local and non-local interactions. Neurocomputing92, 367–379. 10.1007/s00422-005-0565-z
30
KromerJ. A.TassP. A. (2024). Coordinated reset stimulation of plastic neural networks with spatially dependent synaptic connections. Front. Netw. Physiology4, 1351815. 10.3389/fnetp.2024.1351815
31
LaingC.ChowC. (2002). A spiking neuron model for binocular rivalry. J. Comput. Neurosci.12, 39–53. 10.1023/a:1014942129705
32
LaingC. R. (2018). On the application of “equation-free modelling” to neural systems. J. Comput. Neurosci.20, 5–23. 10.1007/s10827-006-3843-z
33
LaingC. R.FrewenT.KevrekidisI. G. (2010). Reduced models for binocular rivalry. J. Comput. Neurosci.28, 459–476. 10.1007/s10827-010-0227-6
34
MaltbyN.TolinD. F.WorhunskyP.O’KeefeT. M.KiehlK. A. (2005). Dysfunctional action monitoring hyperactivates frontal-striatal circuits in obsessive-compulsive disorder: an event-related fmri study. NeuroImage24, 495–503. 10.1016/j.neuroimage.2004.08.041
35
MarschlerC.Faust-EllsässerC.StarkeJ.Van HemmenJ. (2014a). Bifurcation of learning and structure formation in neuronal maps. EPL108, 48005. 10.1209/0295-5075/108/48005
36
MarschlerC.SieberJ.BerkemerR.KawamotoA.StarkeJ. (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. 10.1137/130913961
37
MarshR.HorgaG.ParasharN.WangZ.PetersonB. S.SimpsonH. B. (2014). Altered activation in fronto-striatal circuits during sequential processing of conflict in unmedicated adults with obsessive-compulsive disorder. Biol. Psychiatry75, 615–622. Obsessive-Compulsive Disorder and the Connectome. 10.1016/j.biopsych.2013.02.004
38
MelzerS.GilM.KoserD. E.MichaelM.HuangK. W.MonyerH. (2017). Distinct corticostriatal gabaergic neurons modulate striatal output neurons and motor activity. Cell Rep.19, 1045–1055. 10.1016/j.celrep.2017.04.024
39
MontbrióE.PazóD.RoxinA. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X5, 021028. 10.1103/PhysRevX.5.021028
40
MontgomeryE. B.GaleJ. 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. 10.1016/j.neubiorev.2007.06.003
41
Muñoz-ManchadoA.FoldiC.SzydlowskiS.SjulsonL.FarriesM.WilsonC.et al (2016). Novel striatal gabaergic interneuron populations labeled in the 5ht3a(egfp) mouse. Cereb. Cortex26, 96–105. 10.1093/cercor/bhu179
42
NadlerB.LafonS.CoifmanR. R.KevrekidisI. G. (2006). Diffusion maps, spectral clustering and reaction coordinates of dynamical systems. Appl. Comput. Harmon. Analysis21, 113–127. Special Issue: Diffusion Maps and Wavelets. 10.1016/j.acha.2005.07.004
43
NetoffT.ClewleyR.ArnoS.KeckT.WhiteJ. (2004). Epilepsy in small-world networks. J. Neurosci.24, 8075–8083. 10.1523/JNEUROSCI.1509-04.2004
44
PalkarG.WuJ.-y.ErmentroutB. (2023). The inhibitory control of traveling waves in cortical networks. PLOS Comput. Biol.19, e1010697. 10.1371/journal.pcbi.1010697
45
ParasuramH.NairB.D’AngeloE.HinesM.NaldiG.DiwakarS. (2016). Computational modeling of single neuron extracellular electric potentials and network local field potentials using lfpsim. Front. Comput. Neurosci.10, 65. 10.3389/fncom.2016.00065
46
PopovychO.TassP. (2019). Adaptive delivery of continuous and delayed feedback deep brain stimulation - a computational study. Sci. Rep.9, 10585. 10.1038/s41598-019-47036-4
47
RemijnseP. L.NielenM. M. A.van BalkomA. J. L. M.CathD. C.van OppenP.UylingsH. B. M.et al (2006). Reduced orbitofrontal-striatal activity on a reversal learning task in obsessive-compulsive disorder. Archives General Psychiatry63, 1225–1236. 10.1001/archpsyc.63.11.1225
48
Rodriguez-RomagueraJ.MonteF. H. M. D.QuirkG. J. (2012). Deep brain stimulation of the ventral striatum enhances extinction of conditioned fear. Proc. Natl. Acad. Sci.109, 8764–8769. 10.1073/pnas.1200782109
49
RubinJ.TermanD. (2004). High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J. Comput. Neurosci.16, 211–235. 10.1023/B:JCNS.0000025686.47117.67
50
SchöllE. (2022). Network physiology, insights in dynamical systems: 2021. Front. Netw. Physiol.2, 1–3. 10.3389/fnetp.2022.961339
51
SheQ.ChenG.ChanR. (2016). Evaluating the small-world-ness of a sampled network: functional connectivity of entorhinal-hippocampal circuitry. Sci. Rep.6, 21468. 10.1038/srep21468
52
SieberJ.MarschlerC.StarkeJ. (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. 10.1137/17M1126084
53
SiettosC.StarkeJ. (2016). Multiscale modeling of brain dynamics: from single neurons and networks to mathematical tools. Wiley Interdiscip. Rev. Syst. Biol. Med.8, 438–458. 10.1002/wsbm.1348
54
SiettosC. I. (2011). Equation-free multiscale computational analysis of individual-based epidemic dynamics on networks. Appl. Math. Comput.218, 324–336. 10.1016/j.amc.2011.05.067
55
SpiliotisK.ButenkoK.StarkeJ.van RienenU.KöhlingR. (2024). Towards an optimised deep brain stimulation using a large-scale computational network and realistic volume conductor model. J. Neural Eng.20, 066045. 10.1088/1741-2552/ad0e7c
56
SpiliotisK.ButenkoK.van RienenU.StarkeJ.KöhlingR. (2022a). Complex network measures reveal optimal targets for deep brain stimulation and identify clusters of collective brain dynamics. Front. Phys.10. 10.3389/fphy.2022.951724
57
SpiliotisK.SiettosC. (2011). A timestepper-based approach for the coarse-grained analysis of microscopic neuronal simulators on networks: bifurcation and rare-events micro-to macro-computations. Neurocomputing74, 3576–3589. 10.1016/j.neucom.2011.06.018
58
SpiliotisK.StarkeJ.FranzD.RichterA.KöhlingR. (2022b). Deep brain stimulation for movement disorder treatment: exploring frequency-dependent efficacy in a computational network model. Biol. Cybern.116, 93–116. 10.1007/s00422-021-00909-2
59
StamC.ReijneveldJ. (2007). Graph theoretical analysis of complex networks in the brain. Nonlinear Biomed. Phys.1, 3. 10.1186/1753-4631-1-3
60
SzalisznyóK.SilversteinD. N. (2021). Computational predictions for ocd pathophysiology and treatment: a review. Front. Psychiatry12, 687062. 10.3389/fpsyt.2021.687062
61
TepperJ. M.TecuapetlaF.KoósT.Ibáñez-SandovalO. (2010). Heterogeneity and diversity of striatal gabaergic interneurons. Front. Neuroanat.4 (150), 1–18. 10.3389/fnana.2010.00150
62
TermanD.RubinJ.YewA.WilsonC. (2002). Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J. Neurosci.22, 2963–2976. 10.1523/JNEUROSCI.22-07-02963.2002
63
VenkadeshS.ShaikhA.ShakeriH.BarretoE.Van HornJ. D. (2024). Biophysical modulation and robustness of itinerant complexity in neuronal networks. Front. Netw. Physiology4, 1302499. 10.3389/fnetp.2024.1302499
64
WattsD.StrogatzS. (1998). Collective dynamics of 'small-world' networks. Nature393, 440–442. 10.1038/30918
65
WidgeA. S.ZorowitzS.BasuI.PaulkA. C.CashS. S.EskandarE. N.et al (2019). Deep brain stimulation of the internal capsule enhances human cognitive control and prefrontal cortex function. Nat. Commun.10, 1536. 10.1038/s41467-019-09557-4
66
WuH.HarizM.Visser-VandewalleV.ZrinzoL.CoenenV. A.ShethS. A.et al (2021). Deep brain stimulation for refractory obsessive-compulsive disorder (ocd): emerging or established therapy?Mol. psychiatry26, 60–65. 10.1038/s41380-020-00933-x
67
YagerL.GarciaA.WunschA.FergusonS. (2015). The ins and outs of the striatum: role in drug addiction. Neuroscience301, 529–541. 10.1016/j.neuroscience.2015.06.033
68
ZandtM. V.FlanaganD.PittengerC. (2024). Sex differences in the distribution and density of regulatory interneurons in the striatum. BioRxiv [Preprint]. 10.1101/2024.02.29.582798
Summary
Keywords
network physiology, equation free method, complex network dynamics, obsessive compulsive disorders, control of neurological disorders
Citation
Spiliotis K, Köhling R, Just W and Starke J (2024) Data-driven and equation-free methods for neurological disorders: analysis and control of the striatum network. Front. Netw. Physiol. 4:1399347. doi: 10.3389/fnetp.2024.1399347
Received
11 March 2024
Accepted
16 July 2024
Published
07 August 2024
Volume
4 - 2024
Edited by
Eckehard Schöll, Technical University of Berlin, Germany
Reviewed by
Oleksandr Popovych, Helmholtz Association of German Research Centres (HZ), Germany
Ali Foroutannia, University of Canberra, Australia
Rossella Rizzo, University of Palermo, Italy
Updates
Copyright
© 2024 Spiliotis, Köhling, Just and Starke.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Konstantinos Spiliotis, konstantinos.spiliotis@uni-rostock.de; Wolfram Just, wolfram.just@uni-rostock.de
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.