We made a model of the MSN-MSN network based on a striatal connectome (Hjorth et al., 2020). We only included MSN cells of one type without distinguishing the D1 and D2 subtypes, randomly connected with inhibitory connections. Instead of the detailed cell models used in Hjorth et al. (2020) we used the single compartment cell model used in multiple previous studies (Mahon et al., 2000; Ponzi et al., 2020) including multiple ion-channels, I_Na, I_K, I_Kir, I_Af, I_As, I_Krp, and I_NaP, I_NaS which reproduces the delay to first spike property of MSNs, as can be seen in the illustrative membrane potential trace (Figure 2A).

Experimentally the distribution of the number of MSNs connected to any given MSN peaks around 100 but has a long tail (Hjorth et al., 2020). The mean appears to be around 160 (Hjorth et al., 2020). Here synapses between MSNs are established randomly with a a connection probability of 0.4. This connection probability is slightly higher than it should be, which is around 0.3 (Hjorth et al., 2020), but it keeps simulations a bit more tractable since we only need 400 cells to achieve an average of 160 connections per cell which is probably the most important quantity to respect. Thus, all simulations reported here had 400 cells. Simulations lasted for 10 min. The long simulation periods were necessary for comparison with ex-vivo Ca slice data recorded over several minutes.

The deterministic version (Fuhrmann et al., 2002) of the Tsodyks-Markram synaptic short-term plasticity model (Tsodyks and Markram, 1997) was used for the MSN-MSN synapses, as in Hjorth et al. (2020),

where the short-term depressing variable, R(t) models the fraction of resources available, and short-term facilitation is modeled by the utilization variable, S(t), and δ(t−t_sp) is the presynaptic spike at time t_sp. After a spike the synaptic current increases rapidly in proportion to the product R(t)S(t) and decays exponentially otherwise. The detailed parameters are as follows: synaptic utilization factor (U), 0.41, depression time constant (D), 222 ms, facilitation time constant (F), 1,859 ms. The synaptic exponential decay timescale was 20 ms, and the GABA reversal potential, -85 mV. These parameters were, as described in Hjorth et al. (2020), taken from Table 1 in Planert et al. (2010). Due to the necessity of very long simulations we were not able to explore multiple different synaptic plasticity parameter values. We took the mean values quoted for U and D, while for F we took the mean value plus one standard deviation. We used the slightly longer F to increase the dynamical timescale within the acceptable range, however a few simulations where F was reduced to the mean value did not seem to be much different. MSN-MSN synaptic delays were uniformly random in the range 1–3 ms, (see Supplemental material for full details).

To estimate the network parameters from slice data we performed 200 network simulations of 10 min each, where the strengths of lateral inhibition between MSNs, G_I, and driving excitation to MSNs, G_E, were varied. This latter quantity, G_E, reflects several factors including the intrinsic MSN cellular excitability, which can be controlled by various modulators such as local dopamine concentration, as well as on the strength of cortical and thalamic excitatory synaptic drive, and also on the strength of feedforward inhibition from striatal fast spiking interneurons. For any given simulation, G_E, was uniformly drawn from the range [0, 0.00002] and all MSN cells in the simulation were driven above threshold by somatic current injection (IClamp) with uniform random strength between 0.001305 and 0.001305+G_E which was fixed for the duration of a simulation. For any given simulation, the strength of lateral inhibition, G_I, was uniformly drawn from the range [0.002, 0.02] and all MSN synaptic weights were uniformly random in the range 0.3125G_I–0.9375G_I.

These G_I and G_E ranges cover the acceptable range of IPSP sizes and firing rates we expected to find in the slice data. The maximum firing rate possible occurs when inhibition, G_I, is zero and GE=2×10-5. This generates spiking at around 10 Hz, as shown in Figure 2A. MSNs only spike at higher rates than this when phasically driven in behavior and such rates are not expected to occur in slices. As shown in Figure 2B peak IPSPs at the maximum strength of inhibition, G_I = 0.02, are around 1.5 mV, which is much larger than they are expected to be, while at G_I = 0.002 peak IPSPs are around 0.1 mV which is smaller than they are expected to be. IPSPs are similar to experimentally determined ones (Planert et al., 2010; Hjorth et al., 2020) which peak between around 0.25 and 0.6 mV and decay over around 40 ms, [see Figure 8D(i, ii) in Hjorth et al., 2020]. Model simulation code is available at: https://github.com/adampdp/MSNnetwork. All network simulations are performed using the NetPyNE simulation environment (Dura-Bernal et al., 2019) on the Piz Daint Cray XC40/XC50 supercomputer of the Swiss National Supercomputing Centre (CSCS).

To generate the “calcium activity” time series shown in Figure 3, first spike times are recorded from each cell. Next the spike times are turned into spike count time series using 250 ms non-overlapping time bins. Next these spike count time series are binarized according to whether there is at least one spike in each 250 ms window. This provides a method to compare network spiking activity with calcium imaging data which is known to reflect bursts of spiking, with the same 250 ms temporal resolution.

We used Simulation Based Inference to estimate network parameters from slice data. From each of the 200 ten minute simulations we extracted summary features. We first removed any cells which did not spike at least once during the 10 minute simulation. Next we converted the cell spikes rasters into binary Ca rate time series with 250 ms bins, as described above. For each of these binary time series we calculated three summary features and used these to train the SBI estimator. The first feature was the mean activity. For each cell we divided the total number of active bins by the total number of bins, then we averaged this quantity across all active cells. Next, for each cell we obtained the sequence of time intervals between consecutive active bins. If a pair of temporally adjacent bins were both active, it was counted as a zero interval. We calculated the average interval and coefficient of variation of the intervals (i.e., the interval standard deviation divided by average interval) for each cell. The second and third summary features were the averages of these quantities across all cells. We experimented with other summary features but found these performed best.

Sometimes simulations fall to fixed-point attractor states, reflecting permanent winners-take-all (WTA) states. An example is shown in Supplementary Figure 1A. Since the full permanent WTA is highly pathological and not expected to be found in experimental slices, we removed simulations which displayed this behavior before training the SBI estimator. To do this we first calculated the total fluctuation of the network activity, Supplementary Figure 1B(i). This is simply the coefficient of variation (CV) of the mean network activity time series. Unlike the three features used to train the SBI estimator, its value strongly depends on the number of active cells. If the network activity falls to a fixed point WTA state this quantity should be zero (after removing a transient period). This occurs more often when excitation, G_E, is high and inhibition, G_I, is low, Supplementary Figure 1B(i). We removed all network simulations from the SBI training with values <0.09 of this quantity. Notice that the true values for the total fluctuation calculated from the experimental slices, Supplementary Figure 1B(ii), crosses are much higher than found in any simulations, Supplementary Figure 1B(i), and much higher than the values calculated from the best fit network simulations, Supplementary Figure 1B(ii), brown circles. This is because the value of this quantity depends strongly on the number of recorded active cells which is much lower in experimental slices than in network simulations, and also varies strongly across experimental slices. We found it was not feasible to use this quantity as a summary feature for training the SBI estimator because we would need to train the SBI estimator for each experimental slice separately for its specific quantity of recorded active cells (by sampling the appropriate quantity of the cells from the network simulations).

We used Simulation Based Inference code available at: https://sbi-dev.github.io/sbi/. Simulation-based inference offers a powerful avenue for conducting Bayesian analysis without the need for direct numerical evaluation of the likelihood function. Instead, it relies on accessing simulations generated by the underlying model. The core concept of Simulation-based inference (SBI) involves creating a substantial dataset comprising pairs of model parameters and corresponding simulated data. This dataset then serves as training data for artificial neural networks (ANNs), specifically designed to approximate complex probability distributions. Here we used a Gaussian mixture model with parameters comprising means, covariance matrices, and mixture weights for each component. Training involves minimizing the negative log likelihood of the parameterized distribution evaluated at the training data. These ANNs, in turn, enable the approximation of the likelihood, facilitating the generation of posterior samples. SBI allows for the approximation of the posterior distribution of parameters given observed data, using only forward simulations. This circumvents the need to compute potentially intractable log-likelihood functions and their gradients. By incorporating a prior distribution over parameters, a stochastic simulator, and observed data, SBI yields the posterior distribution that best explains the observed data. Following training, these neural networks are then deployed to analyze experimentally observed data, providing an approximation of the posterior distribution.

We used Neural Posterior Estimation (NPE) with the three summary features described above for all 151 network simulations with total fluctuation exceeding 0.09. Since results depend on the initial conditions of the training we trained the estimator thirty times on the same set of network simulations. Each time we left out a random selection of twenty simulations as a test set. We calculated the mean-squared error between the true inhibition, G_I, and excitation, G_E, parameter values for a given simulation and the peak estimated values predicted by the trained SBI estimator, (GI*-GI)2+9002(GE*-GE)2, and averaged this over the twenty left-out simulations (the factor 900 is necessary because G_I and G_E vary over different ranges). We selected the winning trained SBI with the minimum mean-squared error from the thirty obtained. Estimated posterior distributions calculated for 12 of the 20 test-simulations for this winner SBI estimator are shown in Figure 4. This trained estimator was applied to the Ca data (Figure 5), using the same summary features to obtain the posterior density over the parameters. One million samples were drawn for the posteriors shown in Figures 4, 5.

To find the best fit network simulation corresponding to each experimental slice condition we calculated the mean-squared error, (GI*-GI)2+9002(GE*-GE)2, between the peak estimated parameters for an experimental condition and the true parameters for each network simulation out of the 151 network simulations and used the simulation with the minimum of this quantity for each condition.

The current work utilizes two data sets: (1) data published by Pérez-Ortega et al. (2016), (2) unpublished data from mice (see below). Here for easy reference we briefly describe the previously published data, for full details see Pérez-Ortega et al. (2016). (CT) Control slices. Horizontal corticostriatal slices (300 μm thick) were obtained from postnatal day (20–36) Wistar rats. Microcircuit activity in the control striatum was activated using 5–8 μM NMDA as the excitatory drive. Without an excitatory drive, the striatal circuit is generally silent in the sense that significant peaks of coactive neurons are virtually absent (Jáidar et al., 2010; Plata et al., 2013a). However, in pathological states, the striatal microcircuit is very active in the absence of any excitatory drive (Jáidar et al., 2010; López-Huerta et al., 2013; Plata et al., 2013b). Pathological conditions examined here were: (DEC) Decorticated striatal microcircuits in which the cortex was detached from the striatum in horizontal slices. In these preparations NMDA was still used as excitatory drive (Pérez-Ortega et al., 2016). (PD) Parkinsonian striatal microcircuits. These preparations were obtained from the striatum of rats lesioned ipsilaterally in the substantia nigra pars compacta with 6OHDA. These preparations were tested for dopamine depletion with behavioral techniques. Rats were tested for evoked turning behavior 8 days after the surgery at postnatal day >23. Animals rotating 450 turns or more were selected for further studies. This behavior corresponds to at least 90% of dopamine depletion. NMDA was not used here as the activity is higher than in the controls in the absence of any excitatory drive (Jáidar et al., 2010; López-Huerta et al., 2013; Plata et al., 2013a,b; Pérez-Ortega et al., 2016). (DYS) L-DOPA induced dyskinetic striatal microcircuits. Animals with severe hemiparkinsonism, which increases the risk of development of LID, were selected. Chronic L-DOPA treatment was administered for 15 days. The development of Abnormal Involuntary Movements (AIMs) was monitored. Rats were considered dyskinetic when they exhibited the main AIMs after chronic application of L-DOPA. These preparations were also recorded in the absence of any excitatory drive (Pérez-Ortega et al., 2016).

Additional calcium imaging data were obtained from brain slices from mice (C57BL/6J background) from unpublished experiments using forementioned methodology and preparation according to previous works (Aparicio-Juárez et al., 2019; Calderón et al., 2022) with Calcium Orange (8.2 μM, 0.1% DMSO, 0.67% Pluronic Acid) as calcium indicator. Squared recording sites were in striatum tissue of 750 μm size during spontaneous activity. We detected the rise times of calcium imaging transients, related to burst of action potentials (Carrillo-Reid et al., 2008; García-Vilchis et al., 2019), to obtain the binary spikes raster. For every region of interest (ROIs) we obtained signals, ΔF, corresponding to the average fluorescence of the soma in every frame. Then we applied the following signal processing pipeline (Calderón et al., 2022): (a) estimation of a canonical fluorophore response as a third-degree autoregressive process, (b) noise estimation by wavelet analysis, (c) detrending of fluorescence signals by filtering, and (d) deconvolution by LASSO construction to infer spike activity above 1 standard deviation of the noise: https://github.com/vladscript/finderspiker.

Here we aimed to estimate striatal MSN network model parameters from slice data. Some examples of the previously published Ca imaging slice data we used, obtained from Pérez-Ortega et al. (2016) and Serrano-Reyes et al. (2022), are reproduced in Figure 1. These slices are many minute recordings from control (CT), decorticated (DEC), Parkinsonian (PD), and dyskinetic (DYS) preparations. DEC slices have cortical-striatal transmission removed. PD slices come from rats which had previously been 6OHDA lesioned and developed Parkinsonian symptoms. DYS slices come from 6OHDA lesioned rats which were subsequently treated with L-DOPA until they developed Abnormal Involuntary Movements symptomatic of LID. CT and DEC slices are activated by the addition of NMDA, but PD and DYS slices are spontaneously active without this (see Methods and Pérez-Ortega et al., 2016 for details).

As described in Serrano-Reyes et al. (2022) the non-linear clustering algorithm UMAP has been used to visualize the dynamical structure of the Ca imaging data with time resolution 250 ms. Each time series shows all the active cells in the given slice colored according to their cluster assigned by UMAP. The lower left panels show the temporal development of the total activity across all active cells. The UMAP formalism can also be used to project the clusters into three dimensional space. The projection corresponding to each time series is shown in the far right hand panels. Each point is a cell colored according to its assigned cluster. Some distinct clusters in distinct parts of the three dimensional space can be seen. In particular the clusters appear to reactivate repetitively on exceedingly long timescales of many minutes. Sequential reactivation of clusters can be visualized using the circular projections (middle panels). Here each point represents a UMAP cluster while the colored lines between points represent temporal transitions between sequentially activated clusters (see Serrano-Reyes et al., 2022 for details of how the cluster activation levels at any given time point are determined). The lines start at the cluster of their own color and end in a differently colored cluster. The thickness represents how often the particular transition is observed.

While all slices seem to show sequentially activating cell assemblies, the different pathological conditions appear to show quite different types of patterning and cluster structures (Figure 1, Pérez-Ortega et al., 2016; Lara-González et al., 2019; Serrano-Reyes et al., 2022). Activity in the CT slice (Figure 1A), seems to alternate fairly rapidly between repeated cell assembly activations, for example the red and cyan ones. On the other hand, the DEC slice (Figure 1B), seems to show much lower activity levels with longer lived, less repetitive cluster activations. For example the grey, blue, red, dark green and purple clusters each only activate once, each for a couple of minutes. In the PD slice (Figure 1C), activity is weakly repetitive but the quantity of clusters seems somewhat reduced compared to the CT slice (Figure 1A), and the network appears to get locked into a few dominant states for long periods of many minutes, as described in several studies (Jáidar et al., 2010, 2019; Pérez-Ortega et al., 2016; Aparicio-Juárez et al., 2019; Lara-González et al., 2019). For example very slow alternation appears to exist between yellow and green clusters. Activated bursts for individual cells last longer than those in CT slice (Figure 1A), which are much shorter lived. In contrast DYS slices seem to show higher but more “fractured” activity patterns (Figure 1D). Activity is somewhat similar to CT (Figure 1A), but individual cells activate in shorter bursts with less sequential cluster repetitions. Clusters seem to switch more rapidly than CT slices but also more randomly between many different clusters, as previously described (Pérez-Ortega et al., 2016; Lara-González et al., 2019).

To reproduce these dynamics we extended a previously published model of the striatal MSN network (Ponzi and Wickens, 2010, 2013, 2022; Ponzi et al., 2020) to include short-term plasticity between the MSNs. We used a well validated MSN cell model (Ponzi et al., 2020) with a full complement of ion-channels which has been shown to accurately reproduce characteristics of MSN spiking activity, such as the long delay to first spike (Figure 2A). We investigated a 400 cell network and took the synaptic short-term plasticity and connectivity parameters from a recent detailed striatal connectome study (Hjorth et al., 2020; see Methods). We varied the two most important factors which are known to control network model cell assembly dynamics (Ponzi and Wickens, 2010, 2012, 2013, 2022; Ponzi et al., 2020). These are the strength of recurrent inhibition between MSNs, here denoted G_I, which controls the postsynaptic IPSP size when a presynaptic MSN spikes, and the MSN excitation level, here denoted G_E. Accordingly we varied these two parameters around their physiological values to investigate the dependence of striatal network dynamics on them.

Four example 10 min simulations of model generated calcium activity time series (see Methods), at different levels of G_E and G_I, are shown in Figure 3. The cell Ca activity time series are clustered and ordered using the UMAP algorithm, with exactly the same parameters as used for the slice data (Figure 1). As in the slice data various cell cluster activations can be seen. This demonstrates that this MSN network model is capable of generating exceedingly slowly varying dynamics, whereby clusters repetitively switch on timescales of many 10's of seconds. Not all 400 cells are present because cells which do not activate at all during the time period are not shown. The network simulation time series (Figure 3), tend to look denser than the slice time series (Figure 1), because many more cells are shown. The four different examples, at different levels of G_I and G_E, show quite different types of dynamical activity, assembly patterns, and clustering (see below).

Next, we used Simulation Based Inference (SBI) with sequential neural posterior estimation (SNPE; see https://sbi-dev.github.io/sbi/) to fit the MSN network parameters to Ca slice data. A mixture of Gaussians density estimator was trained to map summary features calculated from MSN network generated time series data to the network parameters (see Methods). We performed 200 network simulations of 10 min each, where the strengths of lateral inhibition, G_I, and excitation, G_E, were varied. Some simulations were rejected due to highly pathological behavior (see Methods). We trained the SBI density estimator on the remaining MSN network simulations after leaving out a further 20 simulations as a test set. Figure 4 shows the SBI parameter estimation for twelve simulations randomly chosen from this test set. The blue lines show the marginal posterior distributions for lateral inhibition, GI*, and cortical excitation, GE*, (here and in the following ^* denotes estimated value, as opposed to actual simulation parameter value) the color plots show the joint posterior density, and the red lines and points indicate the true parameter values, G_I and G_E, for the given network simulation. In most cases the joint posterior density is quite sharply peaked around the true parameter values, and close to zero elsewhere, despite the wide range of different parameter settings shown, including simulations close to the borders of the prior parameter ranges. This suggests the SBI parameter estimation procedure is highly effective.

Finally, we calculated the same summary features from the Ca slice data and used the trained SBI density estimator to map these features to corresponding MSN network parameters. Figure 5 illustrates the results of applying this procedure to the Ca slice data. Here the blue lines indicate the marginal posterior density estimates and the red lines and points indicate the maxima of the posterior estimates. The posterior densities for the four control slices (Figure 5A), are all very sharply peaked. Remarkably, they all indicate very similar tightly clustered levels of estimated inhibition GI*, of around 0.008. These maximal posterior estimate values are shown in Figure 6A, CT. Interestingly the estimated inhibition level GI*, is very close to its physiologically correct value. Indeed, as shown in Figure 2B, peak IPSP sizes around 0.45 mV, the physiologically observed value (Planert et al., 2010; Hjorth et al., 2020), occur between about G_I = 0.006 and G_I = 0.01.

The observation that best fit network models had IPSPs of about the physiologically correct size is highly non-trivial because the SBI density estimator was trained on a wide range of inhibition levels from G_I = 0.002 to G_I = 0.02, which generate a large spread in IPSP sizes (Figure 2B). To investigate this further, we used SBI to estimate the MSN network model parameters from unpublished Ca slices from a different cohort of multiple mice, under the same NMDA activated control conditions. The estimated posterior densities from two of these six slices are shown in Figure 5B, “MICE.” Again the densities are sharply peaked. Figure 6A, MICE shows maximal posterior estimated inhibition GI*, levels in this extra group of six Ca slices. Estimated GI* levels are again fairly tightly clustered between about GI*=0.007 and GI*=0.011 corresponding to peak IPSP sizes of around 0.4 and 0.7 mV which are again close to, but slightly larger than, physiologically observed values (Planert et al., 2010; Hjorth et al., 2020).

Excitation levels for the control groups (Figure 6B, CT, MICE), are more widely spread. The CT group includes three quite closely clustered around GE*=1.25×10-5, while one has lower excitation. The MICE group includes five clustered with very low excitation, around GE*=1.5×10-6, but one with higher excitation. Unlike inhibition, it is not possible to compare these estimated values with empirical physiological cortical-striatal (or thalamic-striatal) EPSP sizes because we are here investigating a slice preparation. All slices are activated by the bath addition of NMDA, since the striatum is purely inhibitory and shows almost no activity without it. Therefore we did not include cortical-striatal synapses or spiking cortical input to the striatum network model. MSN cells are excited purely by somatic current injection (see Methods) modeling the bath addition of NMDA. Moreover, MSNs can change their intrinsic excitability levels depending on various modulatory and pathological factors so that they can spike with different rates for any given level of cortical (or thalamic) excitation. Our parameter G_E also reflects these effects. However, since the mean activity is one of the features used to estimate network parameters (see below), providing that the estimated level of inhibition is roughly correct, as we in fact found it to be, the excitation level should also be reasonable since between them these two factors determine the mean spiking rate. Differences in NMDA application could also cause large differences in excitation levels across the two groups, CT and MICE, and within groups. It is to be expected that excitation will be more variable across experiments than inhibition since the inhibition level is a structural characteristic of the striatum while the excitation level can strongly depend on the experimental conditions. During behavior in live animals, levels of excitation arising from cortical and thalamic driving would also be expected to vary strongly, while the strength of lateral inhibition between MSNs is thought to be more determined by fixed synaptic properties and should be less variable.

To see if our findings in control slices are not simply a chance result of the SBI estimation procedure, we also estimated network parameters from the pathological slice preparations. Posterior distributions generated by the trained SBI density estimator for two DEC slices, two PD slices, and two DYS slices, obtained from previously published work (Pérez-Ortega et al., 2016; Serrano-Reyes et al., 2022), are also shown in Figure 5. Densities are again quite sharply peaked, except for the decorticated preparations (Figure 5C). Interestingly posterior densities tend to be quite similar within each group, but can differ quite strongly between groups. Remarkably, like the control slices (Figure 6A, “CT,” “MICE”), we find that maximal posterior inhibition values GI* (Figure 6A, “DEC,” “PD,” and “DYS”), are quite strongly clustered within each group, but differ between most groups. Peak excitation values, GE* (Figure 6B), are more widely spread but clustering within groups, with different levels of excitation between groups, is still evident.

DEC preparations (Figure 5C), show very high inhibition, GI*, and very low excitation, GE*, as confirmed by peak values for these quantities [Figures 6A, B, “DEC”; the inhibition values, GI*, for the two samples (Figure 6A, “DEC”), are indistinguishable in the figure]. Although colocalized across samples, their estimated densities are much broader than the other preparations (Figure 5C), and they are also close to the borders of the parameter priors we used to train the SBI, suggesting lower confidence in the estimated peak values.

Estimated densities from PD preparations (Figure 5D), are as sharply peaked as control ones (Figure 5A), with similar peak excitation values, GE* (Figure 6B, “PD”), but peak inhibition values, GI* (Figure 6A, “PD”), are lower. In fact estimated peak inhibition values for both PD slice preparations (Figure 6A, “PD”), are lower than all ten estimated peak inhibition values for control slices, CT and MICE (Figure 6A, “CT,” “MICE”), as well as all four DEC and DYS slices (Figure 6A, “DEC,” “DYS”), suggesting that reduced inhibition in 6OHDA slices is a significant finding.

Estimated densities from DYS preparations (Figure 5E), are also quite sharply peaked, although in one sample estimated excitation is on the border of the prior parameter range used to the train the SBI density estimator. Peak inhibition values (Figure 6A, “DYS”), are clustered, but in contrast to the PD slices, are similar to those from control animals (Figure 6A, “CT,” “MICE”). On the other hand, both DYS peak estimated excitation values (Figure 6B, “DYS”) are higher than all of the CT, MICE, PD, and DEC slices (Figure 6B, “CT,” “MICE,” “PD,” and “DEC”), suggesting their excess excitation is also a significant finding.

Given that estimated parameter values were quite strongly clustered within groups, but varied between groups (Figures 6A, B), we wondered what network dynamics these parameter values predict for each of the groups. For each group we calculated the average of the maximal estimated inhibition, GI* and maximal estimated excitation, GE*, levels across all the slices in the group, and searched our 200 network simulations for the simulation with most closely matching parameter values for each group using the mean-square distance (see Methods). The “winner” simulations for CT, DEC, PD, and DYS preparations are the ones shown in Figure 3.

The cluster activation patterns shown by these network simulations admit many similarities with the corresponding experimental slice data (Figure 1). Compared to the best fit CT simulation (Figure 3A, left panel), the best fit DEC simulation (Figure 3B, left panel), shows much sparser activity. Cell clusters are still seen but they activate much more rarely with large silent periods, and without strongly apparent sequential patterns. Some cell clusters, such as the dark blue colored one around cell index number 150 (Figure 3B, left panel), activate strongly only once during the 6 min period shown. Cluster transition graphs (Figure 3B, middle panel), appear to show fewer strong “loops” than the best fit CT simulation (Figure 3A, middle panel). For example the best fit CT simulation (Figure 3A, middle panel), shows a fairly strong recurrent “circuit” composed by the dark blue cluster and purple cluster. These observations are to a certain degree evident in the experimental data (Figures 1A, B). Several fairly strong recurrent loops can be seen in the CT experimental slice (Figure 1A, middle panel), but not in the DEC preparation (Figure 1A, middle panel), while the DEC activity is much sparser than CT activity with large silent gaps between cluster activations, as described above.

Compared to the best fit CT simulation (Figure 3A, left panel), the best fit PD simulation (Figure 3C, left panel), shows larger, longer-lived persistent, and “locked-in” clusters. When such clusters are active the rest of the network is mostly silent. The transition graph is dominated by a strong recurrent loop (Figure 3C, middle panel, red, and green). These findings are recapitulated in the experimental slices. The PD experimental transition graph (Figure 1C, middle panel), is also dominated by a strong recurrent loop between yellow and green clusters, and the PD clusters (Figure 1C, left panel), are much longer lived than the CT experimental slice (Figure 1A, left panel). Finally the best fit DYS simulation displays higher activity levels with many shorter fractured bursts (Figure 3D, left panel), and many weak cluster transitions (Figure 3D, middle panel), compared to the best fit CT simulation (Figure 3A), and again these features are also found in the experimental DYS slice (Figure 1D).

The resemblance of sequential cluster structure of the best fit network simulations to the corresponding experimental slices is a highly non-trivial finding because we did not use any clustering related features to train the SBI density estimator. In fact, because the number of recorded cells varies strongly across experimental slice preparations, we were only able to use features derived from single cell activity, rather than cross-cell correlation based features (see Methods). The features we used were: the mean activity, the mean interval between subsequent activations, and the coefficient of variation (CV) of the intervals between subsequent activations. These quantites were averaged across all active cells to provide the three features for SBI parameter estimation. If the individual cell activation time series in a given network simulation, or in a given experimental slice, are temporally shifted with respect to each other, these feature values remain unchanged. This demonstrates that the sequential cluster structure in the population can be revealed by single cell activity characteristics in this network model.

The dependence of these three features on network model inhibition, G_I and excitation, G_E, parameters, is shown in Figures 6C–E(i). The mean activity [Figure 6C(i)], decreases with increasing inhibition G_I, and at higher inhibition also with decreasing excitation G_E, albeit not as strongly. The mean interval between activations [Figure 6D(i)], increases with increasing inhibition and more weakly with increasing excitation. The behavior of the interval CV [Figure 6E(i)], is non-monotonic (Ponzi and Wickens, 2013, 2022; Ponzi et al., 2020). It increases as inhibition is decreased from high values, before a transition around G_I = 0.006 occurs and it decreases suddenly to very low values, or zero. In fact the transition occurs when the network changes from a strongly fluctuating state with high interval CV, to a WTA fixed point state with low or zero interval CV. In the WTA state the network activity is dominated by a set of regularly firing cells which permanently suppress all the others into silence. An example of a such a “pathological” simulation in the WTA regime is shown in Supplementary Figure 1A.

We found the three features [Figures 6C–E(i)], represent the inhibition, G_I and excitation, G_E, space in different ways, and are sufficient to train the SBI density estimator to accurately map features to parameters (Figure 4). The colored circles in Figures 6C–E(i) show the estimated inhibition GI* and estimated excitation GE*, for the different experimental slices, averaged across the slices in each condition. Most interestingly we find that the two control conditions, CT, and MICE, and the dyskinetic condition, DYS [Figures 6C–E(i)], blue, orange, and purple circles) all reside in the regime with high interval CV just above the transition to the WTA state [Figure 6E(i)], blue, orange, and purple circles). The strong fluctuations found in this regime are generated by the switching coherent cell assembly dynamics seen in the best fit network simulations for CT and DYS conditions (Figures 3A, D). In previous work (Ponzi and Wickens, 2012, 2013, 2022; Ponzi et al., 2020) we have shown that this regime generates complex dynamics which is optimal for neuronal computation (see Discussion). On the other hand the decorticated condition, DEC [Figures 6C–E(i), green circle] resides on the border of the G_I, G_E, parameter space, at very low excitation and very high inhibition and, as described above, its true parameter values may lie outside the prior parameter distribution. Its very low excitation and very high inhibition are responsible for the short lived bursts and very long silent periods found in its corresponding best fit network simulation for the DEC condition (Figure 3). Intriguingly, the reduction in estimated inhibition GI*, we found in the PD experimental slices [Figures 6C–E(i), red circle] situates the PD experimental slices much closer to the border of the pathological WTA regime [Figure 6E(i), red circle]. It is this proximity to the WTA regime which produces the metastable cell assemblies which silence the rest of the network for extended periods, as the system transiently visits the vicinity of WTA states.

Figures 6C–E(ii), colored crosses show the true values of the three features calculated directly from the experimental slices. The brown circles [Figures 6C–E(ii), brown circles], show the feature values calculated directly from the best-fit network simulations in each condition, shown in Figure 3. True features values for the two control conditions, CT and MICE [Figures 6C–E(ii), blue, orange] are very close to the feature values indicated by their corresponding estimated excitation, GE*, and inhibition GI*, levels [Figures 6C–E(i), blue, orange circles], and to the values directly calculated from the best-fit simulations [Figures 6C–E(ii), brown circles]. This confirms that our model is able to accurately reproduce control slice data.

True feature values calculated from the experimental slices in pathological conditions [Figures 6C–E(ii), green, red, and purple], however, can sometimes stray a little from the values indicated by their corresponding estimated excitation GE*, and inhibition GI*, levels [Figures 6C–E(i), green, red, and purple circles] and from their best fit simulation values [Figures 6C–E(ii), brown circles]. True experimental mean activity levels [Figure 6C(ii), colored crosses] are always in good correspondence with their corresponding best fit simulation values Figure 6C(ii), brown circles]. But in the DEC condition, the true mean interval [Figure 6D(ii), green], is much shorter than the corresponding best-fit simulation mean interval [Figure 6D(ii), brown circle]. The true interval CV for DEC [Figure 6E(ii), green] is also lower than the corresponding best-fit interval CV [Figure 6E(ii), brown circle]. For the DYS condition, the true mean interval[Figure 6D(ii), purple], is a bit lower than its corresponding best-fit estimated value [Figure 6D(ii), brown circle]. On the other hand in the PD condition the true interval CV [Figure 6E(ii), red], is actually higher than its corresponding best-fit estimated value [Figure 6E(ii), brown circle]. Indeed the interval CV in the experimental PD condition acquires a very high value, around five, which is both much higher than the other experimental conditions, and higher than found in any of our network simulations [Figure 6E(i)]. And indeed the PD interval CV estimated by SBI [Figure 6E(ii), brown circle], is close to the maximum value, around four, which can be found amongst the available network simulations (see Discussion).