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 m 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 i j = 1 if there is a connection from node j to node i , and A i j = 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.

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 β 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 d V i d t = − I LEAK − I K − I Na − I M/D − I syn + I app (1) where C is the membrane capacity. The current balance Eq. 1 contains four membrane currents, the fast sodium and potassium currents I Na , I K , the leak current I LEAK . For MSN neurons an M-current I M occurs whereas FS neurons contain a D-current I 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 ⋅ ( V − E X ) , where X ∈ { Na, K, M, 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 Na = g Na m Na 3 h Na ⋅ V i − E Na . (2) The potassium current has the form I K = g K m K 4 ⋅ V i − E K . (3) and the leak current reads I LEAK = g LEAK ⋅ V i − E LEAK . (4) Finally, the M- and D-current which enter the MSN and the FS neurons, respectively, are given by I M = g M m M ⋅ V i − E K , I D = g D m D 3 h D ⋅ V i − E D . (5)

The gating variables m Na , h Na , m K , m M , at node i each obey, following the Hodgkin-Huxley formalism (Hodgkin and Huxley, 1952; Chartove et al., 2020), an equation of the type d x i d t = a x i V i − a x i V i + b x i V i x i , (6) where x i ∈ { m Na , h Na , m K , m 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 Na V = 0.32 V + 54 1 − e − V + 54 / 4 ) , b m Na V = 0.28 V + 27 e V + 27 / 5 − 1 and a h Na V = 0.128 e − V + 50 / 18 , b h Na V = 4 1 + e − V + 27 / 5 . The coefficients of the activation gating for the potassium current read a m K V = 0.032 V + 52 1 − e − V + 52 / 5 ) , b m K V = 0.5 e − V + 57 / 40 , and for those of the M-current we have a m M V = 0.032 V + 52 1 − e − V + 52 / 5 ) , b m M V = 0.5 e − V + 57 / 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 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.

The current I app in Eq. 1 is written as I app = I 0 + I 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 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 DBS representing the deep brain stimulation, enters here as an additive contribution. In our analysis we keep I DBS = 0 , except the last section where we discuss the implementation of DBS in our model.

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) d s i d t = α X 1 − s i H X V i − β X s i , (7) where X ∈ { 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 α X and β 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 M V = 1 + tanh V / 4 , with activation rates α M = 2 and β M = 1 / 13 . Similarly, for interneurons the expressions are given by H F V = 1 + tanh V / 10 , with activation rates α F = 4 and β F = 1 / 13 . For a neuron i of type X ∈ { M,F } the total synaptic inhibition it receives from pre-synaptic neurons of type Y ∈ { M,F } is given by I i , GABA = g X Y V i − E GABA ∑ j ∈ Y A i j s j , (8) where A i j is the adjacency matrix of the graph, the summation j ∈ Y is taken over neurons of type Y and E GABA = − 80 mV . The parameter g X Y represents the conductance between X and Y interactions with X , Y ∈ { M , 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 syn = g MM V i − E GABA ∑ j ∈ M A i j s j + g MF V i − E GABA ∑ j ∈ F A i j s j . (9)

Similarly for an interneuron the synaptic current is given by I sys = g FF V i − E GABA ∑ j ∈ F A i j s j + g FM V i − E GABA ∑ j ∈ M A i j 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 MM = g MF = 0.02 and g FF = g 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 Na ) i , ( h Na ) i , ( m K ) i , ( m M ) i , s i ) . Figures 1B,C illustrates the temporal dynamics of the network.

The mean synaptic activity of MSNs turns out to be a suitable macroscopic variable S t = 1 N ∑ i ∈ M s i t . (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 ⋅ Z i where Z i ∼ N ( 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 = Φ 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 app = I 0 as a static parameter of our model (see Eq. 1). Figure 2 shows a graphical representation of the numerical procedure.

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 ′ t = f S t , I 0 ≈ F T S t , I 0 − 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.

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.

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 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.